[CR3BP] 리야프노프 궤도 (Lyapunov Orbit) 계산

리야프노프 궤도(Lyapunov orbit)는 (x-y) 평면에서 라그랑지 포인트 L1, L2, L3를 중심으로 공전하는 주기궤도(periodic orbit)이다. 앞서 살펴본 주기궤도의 조건 (https://pasus.tistory.com/277)에 따라 리야프노프 궤도는 x축에 대해 대칭이고, x축을 직각으로 통과한다.

따라서 시간 \(t_0\) 의 초기조건과 주기 \(T\) 의 반인 시간 \(T/2\) 에서의 상태벡터는 다음과 같아야 한다.

 

\[ \mathbf{x}(t_0 )= \begin{bmatrix} x(t_0) \\ y(t_0) \\ \dot{x}(t_0) \\ \dot{y}(t_0) \end{bmatrix} = \begin{bmatrix} x_0 \\ 0 \\ 0 \\ v_{y0} \end{bmatrix}, \ \ \ \ \mathbf{x}(T/2)= \begin{bmatrix}x(T/2) \\ 0 \\ 0 \\ v_y (T/2) \end{bmatrix} \tag{1} \]

 

리야프노프 궤도를 구하기 위해서는 반주기(\( T/2\)) 시간에서 x축을 직각으로 통과하도록 초기조건 \((x_0, v_{y0})\) 을 잘 선정하는 것이 핵심이다.

 

 

리야프노프 궤도의 계산 절차는 다음과 같다.

 

 

1. 선형화 운동모델의 주기궤도로부터 초기조건의 추측값(initial guess)을 계산한다 (https://pasus.tistory.com/273).

 

\[ \mathbf{x}(t_0 )= \begin{bmatrix} x(t_0) \\ y(t_0) \\ \dot{x}(t_0) \\ \dot{y}(t_0) \end{bmatrix} = \begin{bmatrix} x_e-A_x \\ 0 \\ 0 \\ \omega_p \kappa A_x \end{bmatrix} \tag{2} \]

 

여기서 \(x_e\) 는 라그랑지 포인트 L1, L2, L3의 위치이고 \(A_x \gt 0\) 은 궤도의 x 축 반경이며

 

\[ \begin{align} \kappa &= \frac{\omega_p^2+(1+2c_2 )}{2 \omega_p } \tag{3} \\ \\ c_2 &= \frac{(1-\mu)}{|x_e+ \mu|^3 } + \frac{\mu}{ |x_e+\mu-1|^3 } \\ \\ \omega_p^2 &= \frac{2-c_2+ \sqrt{9c_2^2-8c_2 }}{2} \end{align} \]

 

이다. 이에 해당하는 매트랩 코드는 다음과 같다.

 

% cr3bp linear dynamics for Lypapunov orbit
% input:
%       Lpoint: Lagrange point 1 or 2 (not 3)
%       Ax: x-axis amplitude
%       param
% output:
%       x0: initial x-position
%       vy0: initial y_dot
%       wp: frequency
%       kappa: Ay = kappa*Ax
%
% coded by st.watermelon

% unpack
mu = param.mu;
if Lpoint == 1
    xe = param.L1; % x-position of Lagrange point
else
    xe = param.L2;
end


% parameter computation
c2 = (1-mu)/(sqrt((xe+mu)^2))^3 + mu/(sqrt((xe+mu-1)^2))^3;
wp2 = 1/2 * (2 - c2 + sqrt(9*c2*c2-8*c2));
wp = sqrt(wp2);
kappa = (wp2+1+2*c2)/(2*wp);

% initial guess
x0 = xe - Ax;
vy0 = wp*kappa*Ax;

 

2. 초기조건의 추측값 식 (2)를 이용하여 CR3BP의 무차원 비선형 운동방정식을 \(y=0\) 이 될 때까지 수치적분한다. \(y=0\) 일 때의 시간을 \(t_1\) 으로 표기한다.

 

\[ \begin{align} \dot{\mathbf{x}}(t) &= \begin{bmatrix} \dot{x} \\ \dot{y} \\ \dot{v}_x \\ \dot{v}_y \end{bmatrix} = \mathbf{f}(\mathbf{x}(t)) \tag{4} \\ \\ &= \begin{bmatrix} v_x \\ v_y \\ 2v_y+x-\frac{(1-\mu)(x+\mu)}{r_1^3 }- \frac{\mu (x+\mu-1)}{ r_2^3} \\ -2v_x+y- \frac{(1-\mu)y}{r_1^3 }- \frac{\mu y}{ r_2^3 } \end{bmatrix} \end{align} \]

 

여기서

 

\[ \begin{align} & r_1= \sqrt{ (x+\mu)^2+y^2 } \\ \\ & r_2= \sqrt{ (x+\mu-1)^2+y^2 } \end{align} \]

 

이다. CR3BP 비선형 운동방정식의 매트랩 코드는 다음과 같다.

 

function X_dot = cr3bp_ode_2d(t, X, mu)
% CR3BP nonlinear differential equation 2D version
% input:
%       t
%       X: [x, y, vx, vy]
%       mu
% output:
%       X_dot
%
% coded by st.watermelon

X_dot = zeros(4,1);

% unpack
x = X(1);
y = X(2);
vx = X(3);
vy = X(4);

% paramter computation
r1 = sqrt((x+mu)^2 + y^2);
r2 = sqrt((x+mu-1)^2 + y^2);

% xvec_dot computation
X_dot(1) = vx;
X_dot(2) = vy;
X_dot(3) = 2*vy + x - (1-mu)*(x+mu)/r1^3 - mu*(x+mu-1)/r2^3;
X_dot(4) = -2*vx + y - (1-mu)*y/r1^3 - mu*y/r2^3;

 

3. 시간 \(t_1\) 에서 \(v_x=0\) 인지 확인한다. 아니라면 y축 초기속도 \(v_{y0}\) 를 조정하기 위한 미분보정(differential correction)을 실시해야 한다. 미분보정을 위해서 다음과 같이 상태천이행렬 미분방정식을 만든다 (https://pasus.tistory.com/274).

 

\[ \frac{d}{dt} \Phi (t, t_0 ) = A(t) \Phi (t, t_0 ), \ \ \ \ \Phi(t_0, t_0 )=I \tag{5} \]

 

여기서 \(A(t)\) 는 식 (4)를 수치적분한 궤적 \(\mathbf{x}(t)\) 에서 선형화한 미분방정식의 시스템 행렬로서

 

\[ A(t)= \left. \frac{ \partial \mathbf{f} }{\partial \mathbf{x} } \right|_{\mathbf{x}(t)} = \begin{bmatrix} 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \\ -\bar{U}_{xx} (t) & -\bar{U}_{xy} (t) & 0 & 2 \\ -\bar{U}_{xy} (t) & -\bar{U}_{yy} (t) & -2 & 0 \end{bmatrix}_{\mathbf{x}(t)} \tag{6} \]

 

이며, 각 성분은 다음과 같이 계산된다.

 

\[ \begin{align} & \bar{U}_{xx} (t)= -1- \frac{3(1-\mu) (x+\mu)^2}{r_1^5 }- \frac{3 \mu (x+\mu-1)^2}{r_2^5 } \tag{7} \\ & \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ + \frac{(1-\mu)}{r_1^3 } + \frac{\mu}{r_2^3 } \\ \\ & \bar{U}_{yy} (t)= -1- \frac{3(1-\mu) y^2}{r_1^5 }- \frac{3 \mu y^2}{r_2^5 }+ \frac{(1-\mu)}{r_1^3 }+ \frac{\mu}{ r_2^3 } \\ \\ & \bar{U}_{xy} (t)= - \frac{3(1-\mu)(x+\mu)y}{r_1^5 }- \frac{3 \mu (x+\mu-1) y}{r_2^5 } \end{align} \]

 

식 (5)의 미분방정식을 시간 \(t_1\) 까지 수치적분하여 \( \Phi(t_1, t_0 )\) 을 계산한다. 상태천이행렬 \(\Phi (t, t_0 )\) 는 \(4 \times 4\) 행렬이므로 \(16\) 개의 스칼라 미분방정식이 나오고 식 (4) 또한 \(4\) 개의 스칼라 미분방정식이므로 \(\Phi (t_1, t_0 )\) 을 계산하기 위해서는 총 \(20\) 개의 미분방정식을 수치적분해야 한다. 이에 대한 해당하는 매트랩 코드는 다음과 같다.

 

function Phi_dot = state_trans_ode_2d(t, Phi, mu)
%
% CR3BP state transition matrix differential equation 2D version
% input:
%       t
%       Phi: 16x1 for state transition matrix
%             4x1 for cr3bp state
%       mu
% output
%       Phi_dot: 16x1 for state transition matrix dot
%                 4x1 for cr3bp state
%
% coded by st.watermelon

Phi_dot = zeros(20,1);

% unpack
x = Phi(17);
y = Phi(18);
vx = Phi(19);
vy = Phi(20);

% parameter computation
r1 = sqrt((x+mu)^2 + y^2);
r2 = sqrt((x+mu-1)^2 + y^2);

Uxx = -1 - 3*(1-mu)*(x+mu)^2/r1^5 - 3*mu*(x+mu-1)^2/r2^5 + (1-mu)/r1^3 + mu/r2^3;
Uyy = -1 - 3*(1-mu)*y^2/r1^5 - 3*mu*y^2/r2^5 + (1-mu)/r1^3 + mu/r2^3;
Uxy = -3*(1-mu)*(x+mu)*y/r1^5 - 3*mu*(x+mu-1)*y/r2^5;

% Phi_dot computation
A = [0     0     1  0;
     0     0     0  1;
    -Uxx  -Uxy   0  2;
    -Uxy  -Uyy  -2  0];

PHI = reshape(Phi(1:16),4,4);
PHI_dot = A * PHI;

Phi_dot(1:16) = reshape(PHI_dot,16,1);
Phi_dot(17) = vx;
Phi_dot(18) = vy;
Phi_dot(19) = 2*vy + x - (1-mu)*(x+mu)/r1^3 - mu*(x+mu-1)/r2^3;
Phi_dot(20) = -2*vx + y - (1-mu)*y/r1^3 - mu*y/r2^3;

 

4. 초기속도 \(v_y0\) 는 다음 연립방정식을 풀어서 계산할 수 있다.

 

\[ \begin{align} 0 &= \Phi_{24} \delta v_{y0}+ v_{y1} \delta t_1 \tag{8} \\ \\ \delta v_{x1} &= \Phi_{34} \delta v_{y0}+ \dot{v}_{x1} \delta t_1 \end{align} \]

 

여기서 \(\Phi _{ij}\) 는 상태천이행렬 \(\Phi (t_1, t_0 )\) 의 성분이다. 식 (8)에서 \(v_{y0}\) 를 계산하면 다음과 같다.

 

\[ \begin{align} \delta v_{y0} &= \left( \Phi_{34}- \frac{\dot{v}_{x1}}{v_{y1}} \Phi_{24} \right)^{-1} \delta v_{x1} \tag{9} \\ \\ &= - \left( \Phi_{34}- \frac{\dot{v}_{x1} }{v_{y1} } \Phi_{24} \right)^{-1} v_{x1} \end{align} \]

 

식 (9)에서 \(v_{y0}\) 를 계산한 후 초기값을 \(v_{y0} \ \gets \ v_{y0}+ \delta v_{y0}\) 로 업데이트한다. 다시 2번항으로 돌아가서 시간 \(t_1\) 에서 \(v_x=0\) 이 될 때까지 반복하면 리야프노프 궤도가 계산된다.

여기서 \(y=0\) 과 \(v_x=0\) 의 기준은 \( |y| \le 10^{-11} \sim 10^{-15}, \ \ |v_x | \le 10^{-8} \sim 10^{-15}\) 일 정도로 매우 정밀하게 계산해야 한다.

 

 

2번항과 4번항에 대한 매트랩 코드는 다음과 같다.

 

function [X0, t1, flag] = dp_cr3bp_2d(X0_guess, mu)
%
% Differential correction for CR3BP 2D version
% input:
%       X0_guess: initial guess from linear dynamics
%               = [x0 0 0 vy0]'
%       mu
% output
%       X0: corrected initial condition
%       t1: corrected half period (T/2)
%
% coded by st.watermelon

% setup
max_iter = 100; % maximum DC iteration
absTol = 1e-15;
relTol = 3e-14;
X0 = X0_guess;

for iter = 1: max_iter
    
    fprintf("iter %2d \n", iter)

    % compute t1 and x, vx, vy, and state transition matrix at y=0
    tspan = [0 10];
    Phi0(1:16) = reshape(eye(4),16,1);
    Phi0(17:20) = X0;

    opt2 = odeset('RelTol',relTol,'AbsTol',absTol,'Events',@events);
    [t, Phi, t1, Phi_t1, ie] = ode113(@(t,Phi) state_trans_ode_2d(t,Phi,mu), tspan, Phi0, opt2);

    t1 = t1(end);
    Phi_t1 = Phi_t1(end,:);

    PHI_t1 = reshape(Phi_t1(1:16), 4,4);
    x1 = Phi_t1(17);
    y1 = Phi_t1(18);
    vx1 = Phi_t1(19);
    vy1 = Phi_t1(20);

    r1 = sqrt((x1+mu)^2 + y1^2);
    r2 = sqrt((x1+mu-1)^2 + y1^2);
    vxdot1 = 2*vy1 + x1 - (1-mu)*(x1+mu)/r1^3 - mu*(x1+mu-1)/r2^3;

    % update initial condition
    dvy0 = -vx1/(PHI_t1(3,4)-PHI_t1(2,4)*vxdot1/vy1);
    
    % learning rate
    eta = 1-0.9^iter;
    X0(4,1) = X0(4,1) + eta*dvy0;

    if abs(vx1) <= 1e-11
        break;
    end

end

if iter == max_iter
    flag = 1;
else
    flag = 0;
end

end % end of function


function [value,isterminal,direction] = events(t,Phi)
%
value = Phi(18); % event y=0
isterminal = 1; % stop integration
direction = -1; % negative direction

 

다음은 선형 방정식으로 계산한 궤도와 비선형 방정식으로 계산한 궤도이다. 궤도의 크기가 커질수록 차이가 발생함을 알 수 있다. 리야프노프 궤도의 크기는 궤도에너지에 따라 달라지며 이에 따라 궤도 주기도 달라지므로 임무 요구조건에 따라 궤도 크기를 설계해야 한다.

 

 

다음은 선형화 방정식으로 리야프노프 궤도의 초기값을 획득한 후, x 좌표값을 \(0.001\) 씩 더해 가면서 계산한 지구-달 시스템의 L1과 L2 포인트 리야프노프 궤도들의 집합이다.

 

 

궤도의 크기가 작을 때는 선형화 방정식에서 계산된 궤도의 초기값을 미분보정의 초기 추측값으로 사용해도 잘 수렴하지만 궤도가 커질 수록 초기조건의 추측값에 대해서 미분보정이 매우 민감해 지기 때문에, x 좌표값 \(x_0\) 을 서서히 증가시키며 y 축 속도의 초기 추측값으로 전 단계에서 수렴된 궤도의 초기값 \(v_{y0}\) 을 사용했다.