[CR3BP] 헤일로 궤도 (Halo Orbit) 계산

헤일로 궤도(halo orbit)는 라그랑지 포인트 L1, L2, L3 포인트를 중심으로 형성되는 3차원 주기궤도(periodic orbit)이다.

앞서 살펴본 주기궤도의 조건 (https://pasus.tistory.com/277)에 따라 헤일로 궤도는 (x-z) 평면에 대해 대칭이고, (x-z) 평면을 직각으로 통과한다. 따라서 시간 \(t_0\) 의 초기조건과 주기 \(T\) 의 반인 시간 \(T/2\) 에서의 상태벡터는 다음과 같아야 한다.

 

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

 

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

 

 

헤일로 궤도는 선형화 운동 방정식에서는 존재하는 않는 비선형적인 궤도다. 선형 방정식의 해는 다음과 같았다 (https://pasus.tistory.com/273).

 

\[ \begin{align} x(t) &= -A_x \cos ( \omega_p t+\phi) \tag{2} \\ \\ y(t) &= \kappa A_x \sin (\omega_p t+\phi) \\ \\ z(t) &= A_z \sin (\omega_v t+ \psi) \end{align} \]

 

여기서 진폭을 충분히 키워서 선형화 범위를 넘어서게 한다면 비선형 효과로 인하여 \(\omega_p=\omega_v\) 가 되도록 만들 수 있다. 이 때는 \(A_x\)와 \(A_z\) 는 독립적이지 않고 다음 식으로 주어지는 대수관계를 갖는다.

 

\[ l_1 A_x^2+l_2 A_z^2+\Delta =0 \tag{3} \]

 

따라서 전통적으로 헤일로 궤도를 계산할 때는 \(A_z\) 만 지정한다. 이에 대한 자세한 사항은 Richardson의

논문
"Halo Orbit Formulation for the ISEE-3 Mission, David L. Richardson, J. of Guidance and Control, 1980"

보고서
"Periodic and Quasi Periodic Halo Orbits in the Earth-Sun/Earth-Moon Systems, David L. Richardson, FDSS Contract #NNG10CP02C, Advanced Mission Design, 2010"

을 참고하면 된다.

리야프노프 궤도와 달리 헤일로 궤도는 초기조건에 대한 추측값(initial guess)을 얻기 위해서 선형화 방정식 대신에 Richardson이 논문에서 제안한 3차 근사식을 적용한다. 헤일로 궤도의 계산 절차는 다음과 같다.

1. Richardson의 3차 근사식을 이용하여 초기조건의 추측값을 계산한다.

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

 

\[ \begin{align} \dot{\mathbf{x}}(t) &= \begin{bmatrix} \dot{x} \\ \dot{y} \\ \dot{z} \\ \dot{v}_x \\ \dot{v}_y \\ \dot{v}_z \end{bmatrix} = \mathbf{f}( \mathbf{x}(t)) \tag{4} \\ \\ &= \begin{bmatrix} v_x \\ v_y \\ v_z \\ 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 } \\ - \frac{(1-\mu) z}{r_1^3 }- \frac{\mu z}{ r_2^3 } \end{bmatrix} \end{align} \]

 

여기서

 

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

 

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

 

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

X_dot = zeros(6,1);

% unpack
x = X(1);
y = X(2);
z = X(3);
vx = X(4);
vy = X(5);
vz = X(6);

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

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

 

3. 시간 \(t_1\) 에서 \(v_x=0, \ v_z=0\) 인지 확인한다. 아니라면 x 및 z축 위치 및 y축 초기속도를 조정하기 위한 미분보정(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)\) 에서 선형화한 미분방정식의 시스템 행렬로서

 

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

 

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

 

\[ \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 } \\ \\ \bar{U}_{zz} (t) &= -\frac{3(1-\mu) z^2}{r_1^5 }-\frac{3 \mu z^2}{r_2^5 }+ \frac{(1-\mu)}{r_1^3 }+ \frac{\mu}{ r_2^3 } \\ \\ \bar{U}_{xz} (t) &= - \frac{3(1-\mu)(x+\mu) z}{r_1^5 }- \frac{3\mu (x+\mu-1)z}{r_2^5 } \\ \\ \bar{U}_{yz} (t) &= -\frac{3(1-\mu)yz}{r_1^5 }-\frac{3 \mu yz}{r_2^5 } \end{align} \]

 

식 (5)의 미분방정식을 시간 \(t_1\) 까지 수치적분하여 \(\Phi (t_1, t_0 )\) 을 계산한다. 상태천이행렬 \(\Phi (t, t_0 )\) 는 \( 6 \times 6\) 행렬이므로 \(36\) 개의 스칼라 미분방정식이 나오고 식 (4) 또한 \(6\) 개의 스칼라 미분방정식이므로 \(\Phi (t_1, t_0 )\) 을 계산하기 위해서는 총 \(42\) 개의 미분방정식을 수치적분해야 한다.

 

 

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

 

function Phi_dot = state_trans_ode_3d(t, Phi, mu)
%
% CR3BP state transition matrix differential equation 3D version
% input:
%       t
%       Phi: 36x1 for state transition matrix
%             6x1 for cr3bp state
%       mu
% output
%       Phi_dot: 36x1 for state transition matrix dot
%                 6x1 for cr3bp state
%
% coded by st.watermelon

Phi_dot = zeros(42,1);

% unpack
x = Phi(37);
y = Phi(38);
z = Phi(39);
vx = Phi(40);
vy = Phi(41);
vz = Phi(42);

% parameter computation
r1 = sqrt((x+mu)^2 + y^2 + z^2);
r2 = sqrt((x+mu-1)^2 + y^2 + z^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;
Uzz = -3*(1-mu)*z^2/r1^5 - 3*mu*z^2/r2^5 + (1-mu)/r1^3 + mu/r2^3;
Uxz = -3*(1-mu)*(x+mu)*z/r1^5 - 3*mu*(x+mu-1)*z/r2^5;
Uyz = -3*(1-mu)*y*z/r1^5 - 3*mu*y*z/r2^5;

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

PHI = reshape(Phi(1:36),6,6);
PHI_dot = A * PHI;

Phi_dot(1:36) = reshape(PHI_dot,36,1);
Phi_dot(37) = vx;
Phi_dot(38) = vy;
Phi_dot(39) = vz;
Phi_dot(40) = 2*vy + x - (1-mu)*(x+mu)/r1^3 - mu*(x+mu-1)/r2^3;
Phi_dot(41) = -2*vx + y - (1-mu)*y/r1^3 - mu*y/r2^3;
Phi_dot(42) = -(1-mu)*z/r1^3 - mu*z/r2^3;

 

4. 초기 z축 위치 \(z_0\) 를 고정시켰을 때 (물론 \(x_0\) 를 고정시킬 수도 있다), \(x_0\) 와 \(v_{y0}\) 는 다음 연립방정식을 풀어서 계산할 수 있다.

 

\[ \begin{align} 0 &= \Phi_{21} \delta x_0+ \Phi_{25} \delta v_{y0}+v_{y1} \delta t_1 \tag{8} \\ \\ \delta v_{x1} &= \Phi_{41} \delta x_0+ \Phi_{45} \delta v_{y0}+ \dot{v}_{x1} \delta t_1 \\ \\ \delta v_{z1} &= \Phi_{61} \delta x_0+ \Phi_{65} \delta v_{y0}+ \dot{v}_{z1} \delta t_1 \end{align} \]

 

여기서 \(\Phi_{ij}\) 는 상태천이행렬 \(\Phi (t_1, t_0 )\) 의 성분이다. \(\delta v_{x1}=-v_{x1}, \ \delta v_{z1}=-v_{z1}\) 이므로 식 (8)의 해는 다음과 같다.

 

\[ \begin{bmatrix} \delta x_0 \\ \delta v_{y0} \\ \delta t_1 \end{bmatrix} = \begin{bmatrix} \Phi_{21} & \Phi_{25} & v_{y1} \\ \Phi_{41} & \Phi_{45} & \dot{v}_{x1} \\ \Phi_{61} & \Phi_{65} & \dot{v}_{z1} \end{bmatrix}^{ -1} \begin{bmatrix} 0 \\ -v_{x1} \\ -v_{z1} \end{bmatrix} \tag{9} \]

 

식 (9)에서 \((\delta x_0, \ \delta v_{y0})\) 을 계산한 후 초기값을 다음식으로 업데이트한다.

 

\[ \begin{align} x_0 \ & \gets \ x_0+ \delta x_0 \tag{10} \\ \\ v_{y0} \ & \gets \ v_{y0}+ \delta v_{y0} \end{align} \]

 

다시 2번항으로 돌아가서 시간 \(t_1\) 에서 \(v_x=0, \ v_z=0\) 이 될 때까지 반복하면 헤일로 궤도가 계산된다.

여기서 \( y=0, \ v_x=0, \ v_z=0 \) 의 기준은 \( |y| \le 10^{-11} \sim 10^{-15} \), \( |v_x, \ v_z | \le 10^{-8} \sim ~10^{-15}\) 일 정도로 매우 정밀하게 계산해야 한다. 2번항과 4번항에 대한 매트랩 코드는 다음과 같다.

 

function [X0, t1, flag] = dp_cr3bp_3d(X0_guess, mu)
%
% Differential correction for CR3BP 3D version
% input:
%       X0_guess: initial guess from linear dynamics
%               = [x0 0 z0 0 vy0 0]'
%       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, z, vx, vy, vz and state transition matrix at y=0
    tspan = [0 10];
    Phi0(1:36) = reshape(eye(6),36,1);
    Phi0(37:42) = X0;

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

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

    PHI_t1 = reshape(Phi_t1(1:36), 6,6);
    x1 = Phi_t1(37);
    y1 = Phi_t1(38);
    z1 = Phi_t1(39);
    vx1 = Phi_t1(40);
    vy1 = Phi_t1(41);
    vz1 = Phi_t1(42);

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

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

    % update initial condition x0 and vy0
    PP = [PHI_t1(2,1) PHI_t1(2,5) vy1;
          PHI_t1(4,1) PHI_t1(4,5) vxdot1;
          PHI_t1(6,1) PHI_t1(6,5) vzdot1 ];

    deltv = PP\[0; -vx1; -vz1];

    % learning rate
    eta = 1-0.9^iter;
    X0(1,1) = X0(1,1) + eta*deltv(1);
    X0(5,1) = X0(5,1) + eta*deltv(2);

end

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

end % end of function


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

 

 

 

다음 그림은 \(A_z=0.05\) 로 지정했을 때 Richardson의 3차 근사식 기반 초기조건의 추측값과 미분보정을 통한 초기값을 이용하여 궤도를 시뮬레이션 한 것이다. 미세한 초기조건의 차이가 주기궤도의 구현 여부를 결정함을 알 수 있다.

 

 

 

\(A_z \gt 0\) 이면 'Northern' 궤도, \(A_z \lt 0\) 이면 'Southern' 궤도라고 한다.

 

 

다음 그림은 \(A_z=0.001\) 부터 \(A_z=0.5\) 까지 변화해 가면서 계산한 지구-달 시스템의 L1 포인트 헤일로 궤도들의 집합이다. \(A_z\) 가 작을 때는 리야프노프 궤도와 비슷하지만 \(A_z\) 가 클 수록 (x-y) 평면 밖으로 큰 편위를 갖는다는 것을 알 수 있다.

 

 

 

다음 그림은 지구-달 시스템의 L1 과 L2 포인트 헤일로 궤도들의 집합을 함께 그린 것이다.