[B-Plane] B-평면 타켓팅 - 2

행성 간 임무를 수행하는 동안 우주비행체는 섭동력의 영향을 받아 원하는 임무 궤도에서 벗어날 수 있으므로 실제 B-벡터도 설계한 값과 차이가 발생한다. 이러한 오차를 조정하기 위해서 슈팅방법(shooting method) 기반의 B-평면 타겟팅 방법을 사용하여 우주비행체의 속도벡터를 조정하는데 이를 궤적수정기동(TCM, trajectory correction manuever)이라고 한다.

 

 

목표로 하는 B-벡터를 ${\overrightarrow{B}}^{\star }$, 현재의 B-벡터를 $\overrightarrow{B}$ 라고 하고 테일러 시리즈를 이용하여 전개한 후 1차 절삭하면 B-벡터의 오차 $\mathrm{\Delta }\overrightarrow{B}={\overrightarrow{B}}^{\star }-\overrightarrow{B}$ 와 속도벡터 섭동 $\mathrm{\Delta }\overrightarrow{v}$ 의 관계를 다음과 같이 근사적으로 표현할 수 있다.

 

$$\begin{array}{}\text{(1)}& \mathrm{\Delta }\overrightarrow{B}=\frac{\partial \overrightarrow{B}}{\partial \overrightarrow{v}}\mathrm{\Delta }\overrightarrow{v}\end{array}$$

 

위 식을 각각의 좌표계의 축 성분으로 표현하면 다음과 같다.

 

$$\begin{array}{r}\text{(2)}& & \left[\begin{array}{c}\mathrm{\Delta }{B}_T\\ \mathrm{\Delta }{B}_R\end{array}\right]=\left[\begin{array}{ccc}\frac{\partial B_T}{\partial v_x}& \frac{\partial B_T}{\partial v_y}& \frac{\partial B_T}{\partial v_z}\\ \frac{\partial B_R}{\partial v_x}& \frac{\partial B_R}{\partial v_y}& \frac{\partial B_R}{\partial v_z}\end{array}\right]\left[\begin{array}{c}\mathrm{\Delta }{v}_x\\ \mathrm{\Delta }{v}_y\\ \mathrm{\Delta }{v}_z\end{array}\right]\\ \text{or }& \mathrm{\Delta }{\mathbf{B}}^s=J\mathrm{\Delta }{\mathbf{v}}^a\end{array}$$

 

여기서 $\mathrm{\Delta }{B}_T=B_T^{\star }-B_T,\mathrm{\Delta }{B}_R=B_R^{\star }-B_R$ 은 B-벡터를 B-평면 좌표계인 TRS좌표계 $\{s\}$ 로 표현한 값이고 $\mathrm{\Delta }{v}_x=v_x^{\star }-v_x,\mathrm{\Delta }{v}_y=v_y^{\star }-v_y,\mathrm{\Delta }{v}_z=v_z^{\star }-v_z$ 는 속도벡터를 행성중심관성좌표계 $\{a\}$ 로 표현한 값임에 주의해야 한다. $J$ 는 $2\times 3$ 인 자코비안 행렬이다.

위 식에 의하면 속도벡터 오차 $\mathrm{\Delta }{\mathbf{v}}^a$ 는 다음과 같이 추정할 수 있다.

 

$$\begin{array}{}\text{(3)}& \mathrm{\Delta }{\mathbf{v}}^a=J^T(J{J}^T{)}^{-1}\mathrm{\Delta }{\mathbf{B}}^s\end{array}$$

 

식 (2)에 의하면 미지수가 수식의 갯수보다 많으므로 수식을 더 추가하는 방법을 사용하는 것이 일반적이다. 주로 최근접 도달 시간(TCA, time to closest approach)이나 이심율(eccentricity)에 관한 식을 추가적으로 사용한다.

만약 TCA 식을 추가한다면,

 

$$\begin{array}{}\text{(4)}& \left[\begin{array}{c}\mathrm{\Delta }{B}_T\\ \mathrm{\Delta }{B}_R\\ \mathrm{\Delta }{t}_{TCA}\end{array}\right]=\left[\begin{array}{ccc}\frac{\partial B_T}{\partial v_x}& \frac{\partial B_T}{\partial v_y}& \frac{\partial B_T}{\partial v_z}\\ \frac{\partial B_R}{\partial v_x}& \frac{\partial B_R}{\partial v_y}& \frac{\partial B_R}{\partial v_z}\\ \frac{\partial t_{TCA}}{\partial v_x}& \frac{\partial t_{TCA}}{\partial v_y}& \frac{\partial t_{TCA}}{\partial v_z}\end{array}\right]\left[\begin{array}{c}\mathrm{\Delta }{v}_x\\ \mathrm{\Delta }{v}_y\\ \mathrm{\Delta }{v}_z\end{array}\right]\end{array}$$

 

이고, 이심율 식을 추가한다면 다음과 같이 된다.

 

$$\begin{array}{}\text{(5)}& \left[\begin{array}{c}\mathrm{\Delta }{B}_T\\ \mathrm{\Delta }{B}_R\\ \mathrm{\Delta }e\end{array}\right]=\left[\begin{array}{ccc}\frac{\partial B_T}{\partial v_x}& \frac{\partial B_T}{\partial v_y}& \frac{\partial B_T}{\partial v_z}\\ \frac{\partial B_R}{\partial v_x}& \frac{\partial B_R}{\partial v_y}& \frac{\partial B_R}{\partial v_z}\\ \frac{\partial e}{\partial v_x}& \frac{\partial e}{\partial v_y}& \frac{\partial e}{\partial v_z}\end{array}\right]\left[\begin{array}{c}\mathrm{\Delta }{v}_x\\ \mathrm{\Delta }{v}_y\\ \mathrm{\Delta }{v}_z\end{array}\right]\end{array}$$

 

식 (4)와 (5)에 의하면 자코비안 $J$ 는 $3\times 3$ 행렬이 된다.

 

$$\begin{array}{}\text{(6)}& \mathrm{\Delta }{\mathbf{B}}^s=J\mathrm{\Delta }{\mathbf{v}}^a\end{array}$$

 

따라서 속도벡터 오차 $\mathrm{\Delta }{\mathbf{v}}^a$ 는 다음과 같이 추정할 수 있다.

 

$$\begin{array}{}\text{(7)}& \mathrm{\Delta }{\mathbf{v}}^a=J^{-1}\mathrm{\Delta }{\mathbf{B}}^s\end{array}$$

 

식 (2)와 (7)은 B-벡터의 오차 $\mathrm{\Delta }{\mathbf{B}}^s$ 와 $\mathrm{\Delta }{\mathbf{v}}^a$ 가 선형 관계를 갖는다고 가정했으므로 테일러 시리즈의 고차 항에 의한 오차가 발생한다. 따라서 이러한 절삭 오차를 줄이기 위한 반복 계산 방법이 필요하다 (https://pasus.tistory.com/328). 추정된 속도벡터 오차 $\mathrm{\Delta }{\mathbf{v}}^a$ 를 현재 속도벡터 ${\mathbf{v}}^a$ 에 더해주면 속도벡터를 새로운 값으로 업데이트할 수 있다.

 

$$\begin{array}{}\text{(8)}& {\mathbf{v}}^a& \leftarrow {\mathbf{v}}^a+\mathrm{\Delta }{\mathbf{v}}^a\end{array}$$

 

다음 메트랩 코드는 식 (5)의 자코비안을 수치적으로 계산하는 함수다. 이 함수는 앞선 게시글 (https://pasus.tistory.com/365)에 있는 매트랩 함수 rv2bplane.m 이 필요하다.

 

bplane_jacobian.m

 

function J = bplane_jacobian(r_vec, v_vec_ini, mu, dv)
%
% J = bplane_jacobian(r_vec, v_vec_ini, mu, dv)
% need rv2bplane.m
%
% input
%   r_vec, v_vec_ini : position and starting velocity vector (km, km/s)
%   mu: gravitational parameter of target planet (km^3/s^2)
%   dv: perturbation amount
%
% output
%   J: Jacobian
%
% (c) st.watermelon

J = zeros(3, 3);

% finite differences with +/- vx
v_vec= v_vec_ini + [+dv 0 0]';
[bplane, hyper] = rv2bplane(r_vec, v_vec, mu);
BT_plus = bplane.BT;
BR_plus = bplane.BR;
e_plus = hyper.e;
v_vec= v_vec_ini + [-dv 0 0]';
[bplane, hyper] = rv2bplane(r_vec, v_vec, mu);
BT_minus = bplane.BT;
BR_minus = bplane.BR;
e_minus = hyper.e;
J(:, 1) = 1/(2*dv) * [(BT_plus-BT_minus)
           (BR_plus-BR_minus) 
          (e_plus-e_minus) ];

% finite differences with +/- vy
v_vec= v_vec_ini + [0 +dv 0]';
[bplane, hyper] = rv2bplane(r_vec, v_vec, mu);
BT_plus = bplane.BT;
BR_plus = bplane.BR;
e_plus = hyper.e;
v_vec= v_vec_ini + [0 -dv 0]';
[bplane, hyper] = rv2bplane(r_vec, v_vec, mu);
BT_minus = bplane.BT;
BR_minus = bplane.BR;
e_minus = hyper.e;
J(:, 2) = 1/(2*dv) * [(BT_plus-BT_minus)
           (BR_plus-BR_minus) 
           (e_plus-e_minus) ];

% finite differences with +/- vz
v_vec= v_vec_ini + [0 0 +dv]';
[bplane, hyper] = rv2bplane(r_vec, v_vec, mu);
BT_plus = bplane.BT;
BR_plus = bplane.BR;
e_plus = hyper.e;
v_vec= v_vec_ini + [0 0 -dv]';
[bplane, hyper] = rv2bplane(r_vec, v_vec, mu);
BT_minus = bplane.BT;
BR_minus = bplane.BR;
e_minus = hyper.e;
J(:, 3) = 1/(2*dv) *[(BT_plus-BT_minus)
           (BR_plus-BR_minus)
           (e_plus-e_minus) ];

end

 

 

 

이제 B-평면과 슈팅방법을 이용하여 궤적수정기동(TCM)을 수행해 보자. 시나리오는 다음과 같다.

경사각이 30도인 달 궤도에 진입하려고 하며 목표로 하는 근지점은 고도 1,000km 로서 달의 반경을 고려하면 2,737 km 이다. 기준 위치벡터와 속도벡터를 계산하기 위하여 다음과 같은 궤도요소를 갖는 달의 진입 쌍곡선궤도를 가정한다.

 

$$\begin{array}{}\text{(9)}& r_p=2737,e=1.1,i={30}^0,\mathrm{\Omega }={45}^0,\omega ={60}^0,\theta =-{146.0905}^0\end{array}$$

 

달의 중력파라미터가 $\mu =4903k{m}^3/s^2$ 이므로 변환함수 (https://pasus.tistory.com/288)를 이용하면 기준 위치벡터와 속도벡터를 다음과 같이 계산할 수 있다.

 

$$\begin{array}{}\text{(10)}& {\mathbf{r}}_{ref}^a=\left[\begin{array}{c}3,503.47\\ -37,139.74\\ -32,922.45\end{array}\right]km,{\mathbf{v}}_{ref}^a=\left[\begin{array}{c}-0.32022\\ 0.37905\\ 0.28548\end{array}\right]km/s\end{array}$$

 

여기서 위첨자 $a$ 는 달중심관성좌표계를 의미한다.

 

% data
mu = 4903; % Moon's gravitation constant
rp = 1000+1737; % periapsis
e = 1.1; % eccentricity
i = 30*pi/180; % inclination
W = 45*pi/180; % RAAN
w = 60*pi/180; % AOP
the = -2.54976; % true anomaly (rad)
a = rp/(1-e); % SMA

[r_vec_ref,v_vec_ref] = coe2rv(a,e,i,W,w,the,mu);

 

이 경우 쌍곡선 초과속도(hyperbolic excess velocity) 벡터 ${\overrightarrow{v}}_{\mathrm{\infty }}$ 와 최근접 도달 시간(TCA) 및 $B_T,B_R$ 은 다음과 같이 계산된다.

 

$$\begin{array}{r}\text{(11)}& & {\mathbf{v}}_{\mathrm{\infty }}^a=\left[\begin{array}{c}-0.22998\\ 0.28610\\ 0.21069\end{array}\right]km/s,t_{TCA}=92,554.66sec\\ & B_T=12,524.17km,B_R=-677.97km\end{array}$$

 

[bplane, hyper] = rv2bplane(r_vec_ref, v_vec_ref, mu);
vinf_vec_target = bplane.S_hat*hyper.vinf; % hyperbolic excess velocity vector
tca_ref = hyper.tca; % TCA

% target
[bplane, hyper] = vri2bplane(vinf_vec_target, rp, i, mu);
BT_target = bplane.BT;
BR_target = bplane.BR;
rp_target = hyper.rp;
e_target = hyper.e;

 

여기서 ${\mathbf{v}}_{\mathrm{\infty }}^a$ 와 $r_p$, 그리고 경사각 $i$ 를 설계값으로 하면 식 (11)의 $B_T,B_R$ 가 주 목표점이 되며 $e$ 또는 $t_{TCA}$ 도 목표점으로 추가될 수 있다. 여기서는 목표점을 $B_T,B_R,e$ 로 정한다.

 

target = [BT_target BR_target e_target]';

 

섭동력의 영향을 받아 원하는 임무 궤도에서 벗어난 우주비행체의 위치와 속도벡터를 모사하기 위해 기준 위치벡터에 $\pm 100km$ 의 무작위 섭동량을 가하고, 기준 속도벡터에는 $\pm 0.02km/s$ 의 무작위 섭동량을 가한다.

 

r_vec = r_vec_ref + (rand(3,1)*200 - 100); % +-100 km
v_vec = v_vec_ref + (rand(3,1)*0.04 - 0.02); % +-0.02 km/sec

 

그러면 $B_T,B_R,e$ 에는 다음 크기의 오차가 생긴다 (이 오차는 매트랩 코드를 실행할 때 마다 다른 값으로 나온다).

 

$$\begin{array}{}\text{(12)}& \mathrm{\Delta }{B}_T=2071.28km,\mathrm{\Delta }{B}_R=573.34km,\mathrm{\Delta }e=0.015\end{array}$$

 

슈팅방법을 이용하면 4번의 반복계산 끝에 다음과 같은 결과를 얻을 수 있다.

 

corrected errors 
 BT: 0.000000, 
 BR: -0.000000, 
 e: 0.000000 

initial v vector = 
 -0.334542 
 0.375924 
 0.302105 

corrected v vector = 
 -0.318094 
 0.375469 
 0.292377 

initial and corrected rp =  2124.936346, 2737.000000 
initial and corrected e =  1.084944, 1.100000
initial and corrected i =  31.764900, 30.929764 
reference tca = 92554.660858 
initial and corrected tca =  90046.819988, 92611.192631

 

TCM에 의해서 $B_T,B_R,e$ 오차는 모두 $0$ 이 되는 것을 알 수 있다. 또한 근지점 오차도 $0$ 이 된다. 이 때 TCM에 사용된 속도벡터의 증분은 다음과 같다.

 

$$\begin{array}{}\text{(13)}& \mathrm{\Delta }{\mathbf{v}}^a=\left[\begin{array}{c}0.0206\\ -0.0179\\ -0.0099\end{array}\right]km/s\end{array}$$

 

참고로 동일한 조건에서 식 (3)과 같이 $B_T,B_R$ 만 목표점으로 삼으면 결과는 다음과 같다.

 

target BT, BR 
 12524.172767 
 -677.971612 

initial errors 
 BT: -164.140425, 
 BR: -3650.502267, 
 e: -0.007716 

corrected errors 
 BT: -0.000000, 
 BR: -0.000000, 
 e: -0.001504 

initial v vector = 
 -0.324528 
 0.385272 
 0.272323 

corrected v vector = 
 -0.318547 
 0.376297 
 0.293459 

initial and corrected rp =  2946.047133, 2756.517497 
initial and corrected e =  1.107716, 1.101504
initial and corrected i =  30.823518, 30.975053 
reference tca = 92554.660858 
initial and corrected tca =  92792.920043, 92386.956344

 

이 경우 $B_T,B_R$ 오차는 $0$ 이 되었으나 이심율 오차는 남아있는 것을 알 수 있다. 또한 근지점 오차도 $0$ 이 되지 않았다.

다음은 B-평면 타겟팅 알고리즘이다. 여기서 $\overrightarrow{r}$ 과 $\overrightarrow{v}$ 는 유주비행체의 현재 위치와 속도벡터, ${\overrightarrow{v}}_{\mathrm{\infty }}^{\star },i^{\star },r_p^{\star }$ 는 각각 쌍곡선 초과속도 벡터, 경사각, 근지점 거리의 설계값이다.

0. 입력: ${\overrightarrow{v}}_{\mathrm{\infty }}^{\star },i^{\star },r_p^{\star }$ 및 $\overrightarrow{r},\overrightarrow{v}$

1. ${\overrightarrow{v}}_{\mathrm{\infty }}^{\star },i^{\star },r_p^{\star }$ 부터 $B_T^{\star },B_R^{\star },e^{\star }$ 를 계산한다.

2. $\overrightarrow{r},\overrightarrow{v}$ 로부터 현재의 $B_T,B_R,e$ 를 계산한다

3. $\mathrm{\Delta }{B}_T\mathrm{\Delta }{B}_R\mathrm{\Delta }e$ 를 계산하고 $\overrightarrow{v}$ 를 업데이트한다.

 

 

 

다음은 전체 매트랩 코드다.

 

% B-Plane Targeting example
% need rv2bplane.m, vri2bplane.m, coe2rv
% (c) st.watermelon

% data
mu = 4903; % Moon's gravitation constant
rp = 1000+1737; % periapsis
e = 1.1; % eccentricity
i = 30*pi/180; % inclination
W = 45*pi/180; % RAAN
w = 60*pi/180; % AOP
the = -2.54976; % true anomaly (rad)
a = rp/(1-e); % SMA

[r_vec_ref,v_vec_ref] = coe2rv(a,e,i,W,w,the,mu);
[bplane, hyper] = rv2bplane(r_vec_ref, v_vec_ref, mu);
vinf_vec_target = bplane.S_hat*hyper.vinf; % hyperbolic excess velocity vector
tca_ref = hyper.tca; % TCA
i_ref = hyper.i;

% target
[bplane, hyper] = vri2bplane(vinf_vec_target, rp, i, mu);
BT_target = bplane.BT;
BR_target = bplane.BR;
rp_target = hyper.rp;
e_target = hyper.e;

target = [BT_target BR_target e_target]';

% perturbations in position and velocity vectors
r_vec = r_vec_ref + (rand(3,1)*200 - 100); % +-100 km
v_vec = v_vec_ref + (rand(3,1)*0.04 - 0.02); % +-0.02 km/sec

v_vec_ini = v_vec;
% iteration
dv = 0.0001; % for numerical Jacobian computation 
for iter=1:10
    iter
    
    % 1. compute BT, BR and e from current position and velocity vectors
    [bplane, hyper] = rv2bplane(r_vec, v_vec, mu);
    BT = bplane.BT;
    BR = bplane.BR;
    e = hyper.e;
    current = [BT BR e]';

    % 2. compute Jacobian
    J = bplane_jacobian(r_vec, v_vec, mu, dv);

    % 3. TCM
    correc = inv(J)*(target-current);
    v_vec = v_vec + correc;
    if norm(correc)<1e-8
        break
    end

end

% results
% target
fprintf("target BT, BR and e \n %f \n %f \n %f \n\n", target);

% initial errors
[bplane, hyper] = rv2bplane(r_vec, v_vec_ini, mu);
rp_ini = hyper.rp;
e_ini = hyper.e;
i_ini = hyper.i;
tca_ini = hyper.tca;

fprintf("initial errors \n BT: %f, \n BR: %f, \n e: %f \n\n", ...
    BT_target-bplane.BT, BR_target-bplane.BR, e_target-hyper.e);

% corrected errors
[bplane, hyper] = rv2bplane(r_vec, v_vec, mu);
fprintf("corrected errors \n BT: %f, \n BR: %f, \n e: %f \n\n", ...
    BT_target-bplane.BT, BR_target-bplane.BR, e_target-hyper.e);

% initial and corrected v vectors
fprintf("initial v vector = \n %f \n %f \n %f \n\n", v_vec_ini);
fprintf("corrected v vector = \n %f \n %f \n %f \n\n", v_vec);

% some comparision
fprintf("initial and corrected rp =  %f, %f \n\n", rp_ini, hyper.rp);
fprintf("initial and corrected e =  %f, %f\n\n", e_ini, hyper.e);
fprintf("initial and corrected i =  %f, %f \n\n", i_ini*180/pi, hyper.i*180/pi)
fprintf('reference tca = %f \n', tca_ref)
fprintf("initial and corrected tca =  %f, %f \n\n", tca_ini, hyper.tca);