섭동력을 받는 램버트 문제의 보정 해

램버트 문제의 해(https://pasus.tistory.com/316)는 두 위치 \(\mathbf{r}_1\) 과 \(\mathbf{r}_2\) 사이를 비행하는 데 걸리는 시간 \(\Delta t\) 가 주어졌을 때, 두 위치를 연결하는 이체문제 (two-body problem) 궤적(trajectory)을 계산한다. 하지만 램버트 문제에서 고려하지 않았던 섭동력(perturbation)으로 인하여 궤적이 목표로 한 위치 \(\mathbf{r}_2\) 에 도달하지 못할 때는 어떻게 해야 할까.

 

 

일반적인 섭동력 (J2 섭동력, 태양복사압력, 항력, 달 또는 태양 등의 제3의 중력 등)의 경우 이러한 오차 거리(miss distance)가 크지 않기 때문에, 출발 위치 \(\mathbf{r}_1\) 에서 램버트 솔버가 산출한 속도벡터 \(\mathbf{v}_1\) 을 보정하여 최종 위치 \(\mathbf{r}_2\) 에서 오차 거리를 줄일 수 있다. 이와 같은 방법을 슈팅방법(shooting method)이라고 한다(https://pasus.tistory.com/276).

슈팅방법은 기본적으로 경계값 문제를 초기값 문제로 바꾸어 해를 구하는 방법으로서, 초기값으로부터 산출된 최종값과 목표값 간의 차이를 이용하여 원래 초기값을 수정한다.

최종 위치 오차 \(\Delta \mathbf{r}_2\) 는 출발 위치에서의 속도벡터 섭동 \(\Delta \mathbf{v}_1\) 을 이용하여 다음과 같이 근사적으로 표현할 수 있다.

 

\[ \begin{align} \Delta \mathbf{r}_2= \frac{d \mathbf{r}_2}{d \mathbf{v}_1 } \Delta \mathbf{v}_1 \tag{1} \end{align} \]

 

여기서 미분 \(\frac{d \mathbf{r}_2}{d \mathbf{v}_1 }\) 는 자코비안(Jacobian)으로서 \(3 \times 3\) 행렬인데, 출발 위치에서의 초기 속도변화에 대한 최종 위치의 변화를 나타내므로 '민감도' 행렬이라고도 한다.

식 (1)에 의하면 최종 위치벡터 오차에 대한 초기 속도벡터 오차 \(\Delta \mathbf{v}_1\) 은 다음과 같이 추정할 수 있다.

 

\[ \begin{align} \Delta \mathbf{v}_1= \left( \frac{d \mathbf{r}_2}{d \mathbf{v}_1 } \right)^{-1} \Delta \mathbf{r}_2 \tag{2} \end{align} \]

 

추정된 속도벡터 오차를 속도벡터 \(\mathbf{v}_1\) 에서 빼주면 \(\mathbf{v}_1\) 을 새로운 값으로 업데이트할 수 있다.

 

\[ \begin{align} \mathbf{v}_1 & \gets \mathbf{v}_1- \Delta \mathbf{v}_1 \tag{3} \\ \\ &= \mathbf{v}_1- \left( \frac{d \mathbf{r}_2}{d \mathbf{v}_1 } \right)^{-1} \Delta \mathbf{r}_2 \end{align} \]

 

식 (2)의 자코비안 또는 민감도 행렬 \( \frac{d \mathbf{r}_2}{d \mathbf{v}_1}\) 의 구성 성분을 쓰면 다음과 같다.

 

\[ \begin{align} \frac{d \mathbf{r}_2}{d \mathbf{v}_1} = \mathcal{J} = \begin{bmatrix} \frac{ \partial x_2}{\partial u_1} & \frac{ \partial x_2}{\partial v_1} & \frac{ \partial x_2}{\partial w_1} \\ \frac{ \partial y_2}{\partial u_1} & \frac{ \partial y_2}{\partial v_1} & \frac{ \partial y_2}{\partial w_1} \\ \frac{ \partial z_2}{\partial u_1} & \frac{ \partial z_2}{\partial v_1} & \frac{ \partial z_2}{\partial w_1} \end{bmatrix} \tag{5} \end{align} \]

 

여기서 \(\mathbf{r}_2= \begin{bmatrix} x_2 & y_2 & z_2 \end{bmatrix}^T\), \(\mathbf{v}_1= \begin{bmatrix} u_1 & v_1 & w_1 \end{bmatrix}^T\) 로서 각 벡터를 ECI 좌표계로 표현한 것이다.

이제 J2 섭동력 하에서도 궤적이 시간 \(\Delta t\) 후에 목표 위치 \(\mathbf{r}_2\) 로 정확히 도달할 수 있도록, 램버트가 산출한 초기 속도벡터 \(\mathbf{v}_1\) 을 슈팅방법을 이용하여 보정해 보도록 하겠다.

J2 섭동력을 포함한 램버트 문제는 다음과 같다.

 

\[ \begin{align} & \frac{d^2 \mathbf{r} }{dt^2}+ \frac{\mu}{r^3} \mathbf{r}= \mathbf{a}_p \tag{5} \\ \\ & \mathbf{a}_p= \frac{3}{2} J_2 \frac{\mu}{r^5} R_e^2 \begin{bmatrix} x \left( 5 \left( \frac{z}{r} \right)^2-1 \right) \\ y \left( 5 \left( \frac{z}{r} \right)^2-1 \right) \\ z \left( 5 \left( \frac{z}{r} \right)^2-3 \right) \end{bmatrix} \\ \\ & \mathbf{r}(t_1 )= \mathbf{r}_1, \ \ \ \mathbf{r}(t_2 )= \mathbf{r}_2, \ \ \ \Delta t=t_2-t_1 \end{align} \]

 

여기서 \(\mathbf{r}= \begin{bmatrix} x & y & z \end{bmatrix}^T\), \(r=|\mathbf{r}|\), \(J_2=0.0010826269\) 이며 \(\mathbf{a}_p\) 는 J2 섭동 가속도이다.

식 (4)에 의하면 자코비안 \(\mathcal{J}\) 는 유한 차분을 사용하여 한 번에 한 열씩 수치적으로 계산할 수 있다. 방법은 다음과 같다.

\(\mathcal{J}\) 의 첫번째 열을 계산하기 위하여 \(\mathbf{v}_1\) 의 첫번째 성분에만 \(+\Delta\) 만큼의 조그마한 (예를 들면 \(0.01 km/s\)) 섭동을 준다.

 

\[ \begin{align} \mathbf{v}_1^+= \mathbf{v}_1+ \begin{bmatrix} +\Delta \\ 0 \\ 0 \end{bmatrix} \tag{6} \end{align} \]

 

그리고 \(\mathbf{r}_1\) 과 \(\mathbf{v}_1^+\) 을 초기 조건으로 삼아 식 (5)를 수치적분하여 시간 \(\Delta t\) 후의 위치벡터 \(\mathbf{r}_2^+\) 를 계산한다. 이 단계를 슈팅 단계라고 한다.

그런 후 이번에는 \(\mathbf{v}_1\) 의 첫번째 성분에 반대 방향으로 \(-\Delta\) 만큼의 섭동을 준다.

 

\[ \begin{align} \mathbf{v}_1^+= \mathbf{v}_1+ \begin{bmatrix} -\Delta \\ 0 \\ 0 \end{bmatrix} \tag{7} \end{align} \]

 

마찬가지로 \(\mathbf{r}_1\) 과 \(\mathbf{v}_1^-\) 을 초기 조건으로 삼아 식 (4)를 수치적분하여 시간 \(\Delta t\) 후의 위치벡터 \(\mathbf{r}_2^-\) 를 계산한다. 그러면 자코비안 \(\mathcal{j}\) 의 첫번째 열을 다음과 같이 계산할 수 있다.

 

\[ \begin{align} \mathcal{J}_1= \frac{ \mathbf{r}_2^+- \mathbf{r}_2^-}{2 \Delta } \tag{8} \end{align} \]

 

동일한 방법으로 \(\mathbf{v}_1\) 의 두번째 성분과 세번째 성분에만 각각 동일한 섭동을 주어서 자코비안 \(\mathcal{J}\) 의 두번째 열과 세번째 열을 계산하면 자코비안 \(\mathcal{J}\) 의 전체 행렬을 구할 수 있다.

 

 

자코비안을 구한 후에는 식 (3)으로 초기 속도벡터를 보정하면 된다. 그리고 다시 새로운 속도벡터를 이용하여 자코비안을 다시 계산하고 속도벡터를 보정하는 프로세스를 오차 거리가 일정 수준 이하로 작아질 때 까지 계속 반복한다.

 

 

수치 예를 들어보자. 램버트 문제의 경계조건은 다음과 같다.

 

\[ \begin{align} & \mathbf{r}(t_1 )= \mathbf{r}_1= \begin{bmatrix} 5,000 & 10,000 & 2,100 \end{bmatrix}^T \ \ km \tag{9} \\ \\ & \mathbf{r}(t_2)= \mathbf{r}_2= \begin{bmatrix} -14,600 & 2,500 & 7,000 \end{bmatrix}^T \ \ km \\ \\ & \Delta t=t_2-t_1=3,600 \ sec \ \ (1 \ hr) \end{align} \]

 

위 문제를 범용변수 기반 램버트 알고리즘(https://pasus.tistory.com/316)으로 푼 초기 속도는 다음과 같다.

 

\[ \begin{align} \mathbf{v}(t_1 )= \mathbf{v}_1= \begin{bmatrix} -5.992494639666394 \\ 1.925363415280893 \\ 3.245636528490488 \end{bmatrix} \ \ km/s \tag{10} \end{align} \]

 

식 (5)에서 \(\mathbf{a}_p=0\) (no J2) 으로 하고 \(\mathbf{r}_1\) 과 \(\mathbf{v}_1\) 을 초기 조건으로 하여 최종 위치벡터를 계산하고 오차 거리를 구하면 \(5.954239 \times 10^{-9} \ m\) 로서 램버트 솔버가 매우 정확한 해를 구한 것을 알 수 있다. 하지만 J2를 고려했을 때, 동일한 초기 조건으로 최종 위치벡터를 계산하고 오차 거리를 구하면 \(538.5933 \ m\) 로서 오차 거리가 상당히 커진 것을 볼 수 있다.

이 오차 거리를 감소시키기 위하여 초기 속도벡터를 슈팅 방법으로 보정하면 결과는 다음과 같다. 우선 3회의 반복횟수 만에 초기 속도벡터가 다음 값으로 수정되었다.

 

\[ \begin{align} \mathbf{v}_{1_{corrected}}= \begin{bmatrix} -5.992359526670324 \\ 1.924215530723804 \\3.246086959795341 \end{bmatrix} \ \ km/s \tag{11} \end{align} \]

 

보정된 속도와 원래 속도의 크기 차이는 \(0.1240477 \ m/s\) 에 불과할 만큼 매우 작다. 하지만 이러한 속도 차이가 상당한 크기의 오차 거리를 유발한 것이다.

이제 보정된 속도벡터를 초기값으로 하여 최종 위치벡터를 계산하고 오차 거리를 구하면 \(2.746990 \times 10^{-9} \ m\) 로서 슈팅 방법으로 수행한 보정이 매우 정확하다는 것을 알 수 있다.

다음은 자코비안 행렬 (또는 민감도 행렬)을 계산하는 알고리즘을 매트랩 코드로 구현한 것이다.

 

jacobian_shooting.m

 

function J = jacobian_shooting(odefun, y0ini, tof, dv)

% J = jacobian_shooting(odefun, y0ini, tof, dv)
%        given odefun, y0in, tof, dv
%        compute 3x3 Jacobian for shooting method
%       
% Input  odefun: ODE function
%        y0ini: initial position and velocity vector in km, km/s
%             = [x y z vx vy vz]'
%        tof: time of flight in sec
%        dv: small velocity amoint used for finite differences in km/s
% Output J: Jacobian
%
% all in inertial frame (ECI)
%
% coded by st.watermelon

J = zeros(3, 3);

tspan = [0:tof/100:tof]';
OPTIONS = odeset('RelTol',3e-14,'AbsTol',1e-15);

% finite differences with +/- vx
y0 = y0ini + [0 0 0 +dv 0 0]';
[t,y] = ode113(odefun, tspan, y0, OPTIONS);
rplus = y(end,1:3)';
y0 = y0ini + [0 0 0 -dv 0 0]';
[t,y] = ode113(odefun, tspan, y0, OPTIONS);
rminus = y(end,1:3)';
J(:, 1) = (rplus-rminus)/(2*dv);

% finite differences with +/- vy
y0 = y0ini + [0 0 0 0 +dv 0]';
[t,y] = ode113(odefun, tspan, y0, OPTIONS);
rplus = y(end,1:3)';
y0 = y0ini + [0 0 0 0 -dv 0]';
[t,y] = ode113(odefun, tspan, y0, OPTIONS);
rminus = y(end,1:3)';
J(:, 2) = (rplus-rminus)/(2*dv);

% finite differences with +/- vz
y0 = y0ini + [0 0 0 0 0 +dv]';
[t,y] = ode113(odefun, tspan, y0, OPTIONS);
rplus = y(end,1:3)';
y0 = y0ini + [0 0 0 0 0 -dv]';
[t,y] = ode113(odefun, tspan, y0, OPTIONS);
rminus = y(end,1:3)';
J(:, 3) = (rplus-rminus)/(2*dv);

 

다음은 식 (5)의 운동 방정식을 구현한 코드다.

 

orbitbody_J2.m

 

function ydot=orbitbody_J2(t,y)

% two-body orbital equation with J2
% coded by st.watermelon

mu=398600;
Rearth=6378;
J2=1.0826269e-3;

r=y(1:3);
v=y(4:6);

normr=norm(r);
rhat=r/normr;

acc=-(mu/(normr*normr))*rhat;

% J2
tmp = 5*r(3)^2/normr^2;
acc=acc + 3/2*J2*mu/normr^5*Rearth^2* ...
          [r(1)*(tmp-1); r(2)*tmp-1; r(3)*(tmp-1)];

%acc=acc - (mu*J2*Rearth*Rearth/2)* ...
%    ((6*z/(normr^5))*khat+((3/(normr^4))-(15*z*z/(normr^6)))*rhat);
ydot=zeros(6,1);
ydot(1:3)=v;
ydot(4:6)=acc;

 

다음은 수치 예제를 실행한 파일이다.

 

orbit_simul.m

 

% 
% J2-perturbed Lambert correction with shooting method
% coded by st.watermelon

clear all


% data
r1=[5000 10000  2100]';
r2=[-14600  2500  7000]';
tof=60*60; % 1 hour

OPTIONS = odeset('RelTol',3e-14,'AbsTol',1e-15);

% 1. Lambert solution
[v1, v2] = lambert(r1, r2, tof);

% 2. check solution accuracy  
y0=[r1;v1];
[t,y]=ode113(@orbitbody, [0:60:tof]',y0,OPTIONS); 

r2_computed = y(end,1:3)';
err_is = r2_computed-r2;

fprintf('Lambert solver accuracy \n');
fprintf('final error distance = %e m \n\n', norm(err_is)*100);


% 3. check again with J2 perturbation
y0=[r1;v1];
[t,yJ2]=ode113(@orbitbody_J2, [0:60:tof]',y0,OPTIONS);

r2_J2computed = yJ2(end,1:3)';
errJ2_is = r2_J2computed-r2;

fprintf('J2-perturbed Lambert solver accuracy \n');
fprintf('final error distance = %e m \n\n', norm(errJ2_is)*100);

% 4. correction with shooting method
v1_J2 = v1;
for iter=1:10
    y0=[r1;v1_J2];
    [t,y]=ode113(@orbitbody_J2, [0:60:tof]',y0,OPTIONS);
    J = jacobian_shooting(@orbitbody_J2, [r1; v1_J2], tof, 0.01);
    correc = inv(J)*(y(end,1:3)'-r2);
    v1_J2 = v1_J2 - correc;
    if norm(correc)<1e-8
        break
    end
end
v1_J2
fprintf('shooting iteration = %d \n', iter);
fprintf('initial velocity differences = %e m/s \n\n', norm(v1-v1_J2)*100);


% 5. check corrected accuracy
y0=[r1;v1_J2];
[t,yJ2_c]=ode113(@orbitbody_J2, [0:60:tof]',y0,OPTIONS);

r2_J2c = yJ2_c(end,1:3)';
errJ2c_is = r2_J2c-r2;

fprintf('corrected solution accuracy \n');
fprintf('final error distance = %e m \n', norm(errJ2c_is)*100);