[CR3BP] 주기궤도의 매니폴드 계산

라그랑지 포인트(Lagrange point) 에 대한 안정성(stability) 판별과 부분공간(subspace)의 계산 (https://pasus.tistory.com/272)과 유사하게 주기궤도(periodic orbit) 상에 있는 임의의 포인트에 대해서도 안정성 판별과 부분공간을 계산할 수 있다.

주기궤도(periodic orbit) 상에 고정된 포인트에서 계산된 모노드로미 행렬 (monodromy matrix)은 궤도 상에 있는 포인트마다 서로 다른 값을 가지므로 고유벡터(eigenvector)는 달라진다. 반면에 고유값(eigenvalue)은 궤도를 따라 일정하게 유지되는데, 이 때문에 고유값을 '주기궤도의 고유값' 이라고 하며 주기궤도의 한 속성으로 본다.

 

 

이에 대해 자세히 알아보기 위하여 먼저 주기(period)가 \(T\) 인 주기궤도 방정식을 선형화한 미분방정식 (https://pasus.tistory.com/285)을 분석해 보자.

 

\[ \delta \dot{\mathbf{x}} (t)=A(t) \delta \mathbf{x}(t), \ \ \ \ \ A(t)=A(t+T) \tag{1} \]

 

여기서 시스템 행렬 \(A(t)\) 의 기본행렬(fundamental matrix)과 상태천이행렬(state transition matrix)을 각각 \(\Psi(t)\) 와 \(\Phi (t, \tau)\) 라고 하자 (https://pasus.tistory.com/274). 그러면 다음 식이 성립한다.

 

\[ \Phi (t_0+t, t_0 )= \Psi (t_0+t) \Psi^{-1} (t_0 ) \tag{2} \]

 

모노드로미 행렬을 \(M\) 이라고 하면, \(M=e^{\bar{A} T}\) 인 상수 행렬 \(\bar{A}\) 와 다음 식을 만족하는 주기 \(T\) 인 주기 행렬 \(P(t)=P(t+T)\) 가 존재한다.

 

\[ P(t)= \Psi (t) e^{-\bar{A}t} \tag{3} \]

 

식 (3)을 (2)에 대입하면 다음과 같이 된다.

 

\[ \begin{align} \Phi (t_0+t, t_0 ) &= \Psi (t_0+t) \Psi^{-1} (t_0 ) \tag{4} \\ \\ &= P(t_0+t) e^{\bar{A}(t_0+t) } e^{-\bar{A}t_0 } P^{-1} (t_0 ) \\ \\ &= P(t_0+t) e^{\bar{A}t} P^{-1} (t_0) \end{align} \]

 

행렬 \(\bar{A}\) 의 고유값 분해(eigen decomposition)를 다음과 같다고 하면,

 

\[ \bar{A} =Q \Lambda Q^{-1} \tag{5} \]

 

여기서 \(\Lambda\) 는 \(\bar{A}\) 의 고유값 (Floquet exponent)으로 구성된 (블록)대각행렬이고 \(Q\) 는 고유벡터를 열로 하는 행렬이다. 그러면,

 

\[ e^{\bar{A}T} =Q e^{\Lambda T} Q^{-1} \tag{6} \]

 

이 된다. 식 (6)을 (4)에 대입하고 \(t=T\) 에서 계산하면 다음과 같이 된다.

 

\[ \Phi (t_0+T,t_0 )=P(t_0+T) Q e^{\Lambda T} Q^{-1} P^{-1} (t_0) \tag{7} \]

 

여기서 \(P(t)=P(t+T)\) 이고 \(\Phi (t_0+T, t_0 )=M\) 이므로 위 식은 다음과 같이 된다.

 

\[ \begin{align} M &= P(t_0 )Q e^{\Lambda T} Q^{-1} P^{-1} (t_0 ) \tag{8} \\ \\ &= \bar{P}(t_0 ) e^{\Lambda T} \bar{P}^{-1} (t_0) \end{align} \]

 

식 (8)에 의하면,

 

\[ M \bar{P}(t_0 )= \bar{P} (t_0 ) e^{\Lambda T} \tag{9} \]

 

가 되므로 \(e^{\Lambda T}\) 는 모노드로미 행렬 \(M\) 의 고유값 (Floquet multiplier)으로 구성된 (블록)대각행렬이 되고 \(\bar{P}(t_0 )\) 는 그에 대응하는 고유벡터를 열로 하는 행렬이 된다. 따라서 고유값은 상수이므로 주기궤도의 모든 포인트에서 일정하며, 고유벡터는 \(t_0\) 의 함수이므로 궤도 상에 있는 포인트마다 서로 다른 값을 갖게 된다.

참고로 \(M\) 의 고유값을 \(\mu\), \(\bar{A}\) 의 고유값을 \(\gamma\) 라고 하면 다음 식이 성립한다.

 

\[ \mu = e^{\gamma T} \tag{10} \]

 

모노드로미 행렬에서 \(|\mu | \lt 1\) 인 고유값은 안정한 고유값으로서 이에 해당하는 고유벡터는 안정한 부분공간을 형성하며, 안정한 부분공간에서 시작하는 궤적은 \(t \to \infty\) (순방향 전파)함에 따라 주기궤도에 점근적으로 접근한다. 반면 \(|\mu | \gt 1\) 인 고유값은 불안정한 고유값으로서 이에 해당하는 고유벡터는 불안정한 부분공간을 형성하며, 불안정한 부분공간에서 시작하는 궤적은 \(t \to -\infty \) (역방향 전파)함에 따라 주기궤도에 점근적으로 접근한다. 글로벌 매니폴드(manifold)는 안정 부분공간과 불안정 부분공간을 이용하여 계산할 수 있다 (https://pasus.tistory.com/234).

 

 

CR3BP에서 모노드로미 행렬의 고유값은 역쌍으로 나오므로 주기궤도가 불안정한 매니폴드를 가지고 있으면 안정한 매니폴드도 존재한다. 또한 모노드로미 행렬의 6개의 고유값 중 2개는 항상 1이 된다 (https://pasus.tistory.com/285).

우주비행체가 안정적인 매니폴드에 있으면 결국 주기궤도에 도달하게 되고 불안정한 매니폴드에 있으면 궤도에서 점점 멀어지게 된다. 이 속성은 라그랑지 포인트에 대한 우주임무 설계에 매우 중요하다. 특히 안정한 매니폴드는 탄도적으로 궤도에 접근하는 자연스러운 경로를 제공하기 때문에 지구 파킹궤도에서 라그랑지 포인트의 주기궤도로 이동하기 위한 천이궤도(transfer orbit) 설계에 활용된다.

고유벡터는 궤도 상에 있는 포인트마다 방향이 바뀌므로 안정한 부분공간과 불안정한 부분공간의 방향도 궤도를 따라가며 달라진다. 또한 고유벡터에 음수를 곱해도 고유벡터가 되므로 매니폴드는 서로 반대인 두 방향으로 궤도에 접근(안정) 하거나 이탈(불안정)하게 된다.

 

 

궤도를 따라 임의의 한 포인트에 대해서 모노드로미 행렬이 얻어지면 궤도를 따라 다른 임의의 포인트와 관련된 고유벡터를 계산해야 하는데, 이에는 두가지 방법이 있다. 첫째는 새로운 포인트에서 모노드로미 행렬을 다시 구한 후 고유벡터를 계산하는 방법이고, 둘째는 상태천이행렬을 이용하는 방법이다.

임의의 포인트 \(\bar{\mathbf{x}} (t_0 )\) 에서 계산된 안정한 고유벡터를 \(Y^s (t_0)\) 로, 불안정한 고유벡터를 \(Y^u (t_0)\) 라 하면, 포인트 \(\bar{\mathbf{x}} (t_i )\) 에서의 안정 및 불안정한 고유벡터는 상태천이행렬을 이용하여 각각 다음과 같이 계산할 수 있다.

 

\[ \begin{align} Y^s (t_i) &= \Phi (t_i, t_0 ) Y^s (t_0) \tag{11} \\ \\ Y^u (t_i) &= \Phi (t_i, t_0 ) Y^u (t_0) \end{align} \]

 

주기궤도 \(\bar{\mathbf{x}} (t)\) 가 주어졌을 때, 매니폴드를 계산하는 절차는 아래와 같다.

먼저 주기궤도 상의 임의의 포인트 \(\bar{\mathbf{x}} (t_i )\) 에서 안정한 고유벡터와 불안정한 고유벡터 방향으로 섭동을 준다.

 

\[ \begin{align} X^s (t_i ) &= \bar{\mathbf{x}} (t_i )+ d \ Y^s (t_i ) \tag{12} \\ \\ X^u (t_i ) &= \bar{\mathbf{x}} (t_i )+ d \ Y^u (t_i ) \end{align} \]

 

여기서 \(d\) 는 주기궤도로부터의 섭동량이다. \(d\) 의 크기는 선형화 가정의 범위내에 있을 만큼 작야야 하지만 너무 작으면 매니폴드 계산시간이 길어질 수 있다. 어떤 연구에 의하면 \(Y^s (t_i )\) 와 \(Y^u (t_i )\) 대신에 고유벡터의 위치 성분만으로 정규화한 안정 및 불안정 고유벡터 \(V^s (t_i )\) 와 \(V^u (t_i )\) 를 사용하기도 한다.

 

\[ \begin{align} X^s (t_i ) &= \bar{\mathbf{x}} (t_i )+ d \ V^s (t_i ) \tag{13} \\ \\ X^u (t_i ) &= \bar{\mathbf{x}} (t_i )+ d \ V^u (t_i ) \end{align} \]

여기서

\[ \begin{align} & V^s (t_i )= \frac{ Y^s (t_i ) }{ \sqrt{x_s^2+y_s^2+z_s^2} } \\ \\ & V^u (t_i )= \frac{ Y^u (t_i ) }{ \sqrt{x_u^2+y_u^2+z_u^2} } \\ \\ & Y^s (t_i )= [x_s \ \ y_s \ \ z_s \ \ v_{xs} \ \ v_{ys} \ \ v_{zs} ]^T \\ \\ & Y^u (t_i )= [x_u \ \ y_u \ \ z_u \ \ v_{xu} \ \ v_{yu} \ \ v_{zu} ]^T \end{align} \]

 

이다. \(d\) 와 \(-d\) 를 모두 사용하여 매니폴드의 초기조건을 계산한 후, 시간적으로 순방향으로 수치적분을 수행하여 주기궤도의 불안정한 매니폴드인 \(W^{u+}\) 와 \(W^{u-}\) 의 궤적을 생성한다. 안정한 매니폴드인 \(W^{s+}\) 와 \(W^{s-}\) 의 궤적은 시간적으로 역방향으로 수치적분을 수행하여 찾을 수 있다.

 

 

다음 그림은 L1 리야프노프 궤도 상에 있는 50개 포인트에서 계산한 매니폴드이다. 파란색은 안정한 매니폴드, 빨간색은 불안정한 매니폴드이다.

 

 

다음 그림은 L1 헤일로 궤도 상에 있는 50개 포인트에서 계산한 매니폴드이다. 파란색은 안정한 매니폴드, 빨간색은 불안정한 매니폴드이다.

 

 

다음은 매니폴드의 초기조건을 계산하기 위한 매트랩 코드이다.

 

function [X0Wsp,X0Wsn,X0Wup,X0Wun] = halo_manifold_init(x0,T,n_points,mu)
%
% x0W = halo_manifold_init(x0, T, n_points, stable, dir)
% compute initial point of Lyapunov orbit manifold
% input:
%       x0: initial condition of the selected Lyapunov orbit
%       T: period of the orbit
%       n_points: number of points on the orbit for computing manifold
%       mu
% output
%       X0W =initial point on manifold of the fixed point on the orbit
%           =[ x, y, z, vx, vy, vz] (n_points x 6)
%       X0Wsp: manifold of stable, positive branch 
%       X0Wsn: manifold of stable, negative branch
%       X0Wup: manifold of unstable, positive branch
%       X0Wun: manifold of unstable, negative branch
%
% coded by st.watermelon

N = 6;
X0W = zeros(n_points,N);

% compute state transition matrix(STM) at n_points
absTol = 1e-15;
relTol = 3e-14;

tspan = linspace(0,T,n_points+1);
Phi0(1:N*N) = reshape(eye(N),N*N,1);
Phi0(N*N+1:N*N+N) = x0;

opt1 = odeset('RelTol',relTol,'AbsTol',absTol);
[t, Phi] = ode113(@(t,Phi) state_trans_ode_3d(t,Phi,mu), tspan, Phi0, opt1);

Phi_t = Phi(:, 1:N*N); % Phi(t,t_0)
x_t = Phi(:, N*N+1:end); % x(t)

% compute monodromy matrix
M = reshape(Phi_t(end,1:N*N), N,N);

% compute stable/unstable eigenvalues and eigenvectors of monodromy
[se,ue,Ws,Wu] = e_space(M);

eps = 1e-6;
% compute initial points on manifold at t0
% stable
d = eps/norm(Ws(1:3,1));
x0w = x0 + d*Ws(:,1);
X0Wsp(1,:) = x0w';
x0w = x0 - d*Ws(:,1);
X0Wsn(1,:) = x0w';

% unstable
d = eps/norm(Wu(1:3,1));
x0w = x0 + d*Wu(:,1);
X0Wup(1,:) = x0w';
x0w = x0 - d*Wu(:,1);
X0Wun(1,:) = x0w';

% compute initial points using eigenvectors at n_points with STM
for kk=2:n_points
    PHI_t = reshape(Phi_t(kk,1:N*N), N,N);

    % stable 
    Wst = PHI_t*Ws(:,1);
    d = eps/norm(Wst(1:3,1));
    x0w = x_t(kk,:)' + d*Wst;
    X0Wsp(kk,:) = x0w';
    x0w = x_t(kk,:)' - d*Wst;
    X0Wsn(kk,:) = x0w';
    
    % unstable
    Wut = PHI_t*Wu(:,1);
    d = eps/norm(Wut(1:3,1));
    x0w = x_t(kk,:)' + d*Wut;
    X0Wup(kk,:) = x0w';
    x0w = x_t(kk,:)' - d*Wut;
    X0Wun(kk,:) = x0w';

end

end


function [se,ue,Ws,Wu] = e_space(M)

[mm,nn] = size(M);
[V, D] = eig(M);

eps = 1e-7;
s = 0; u = 0; c = 0; % number of stable, unstable, center eigenvalues
for k = 1:mm
    if (abs(D(k,k))-1)<-eps % stable |lambda|<1
        s = s+1;
        se(s) = D(k,k);
        Ws(:,s) = V(:,k);
    elseif (abs(D(k,k))-1)>eps % unstable |lambda|>1
        u = u+1;
        ue(u) = D(k,k);
        Wu(:,u) = V(:,k);
    else % marginal |lambda|=1
        c = c+1;
        ce(c) = D(k,k);
        Wc(:,c) = V(:,k);
    end
end

end