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

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

 

 

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

 

\[ \begin{align} \Delta \vec{B}= \frac{ \partial \vec{B} }{\partial \vec{v} } \Delta \vec{v} \tag{1} \end{align} \]

 

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

 

\[ \begin{align} & \begin{bmatrix} \Delta B_T \\ \Delta B_R \end{bmatrix} = \begin{bmatrix} \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{bmatrix} \begin{bmatrix} \Delta v_x \\ \Delta v_y \\ \Delta v_z \end{bmatrix} \tag{2} \\ \\ \mbox{or } \ \ \ & \Delta \mathbf{B}^s =J \Delta \mathbf{v}^a \end{align} \]

 

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

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

 

\[ \begin{align} \Delta \mathbf{v}^a=J^T (JJ^T )^{-1} \Delta \mathbf{B}^s \tag{3} \end{align} \]

 

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

만약 TCA 식을 추가한다면,

 

\[ \begin{align} \begin{bmatrix} \Delta B_T \\ \Delta B_R \\ \Delta t_{TCA} \end{bmatrix} = \begin{bmatrix} \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{bmatrix} \begin{bmatrix} \Delta v_x \\ \Delta v_y \\ \Delta v_z \end{bmatrix} \tag{4} \end{align} \]

 

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

 

\[ \begin{align} \begin{bmatrix} \Delta B_T \\ \Delta B_R \\ \Delta e \end{bmatrix} = \begin{bmatrix} \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{bmatrix} \begin{bmatrix} \Delta v_x \\ \Delta v_y \\ \Delta v_z \end{bmatrix} \tag{5} \end{align} \]

 

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

 

\[ \begin{align} \Delta \mathbf{B}^s =J \Delta \mathbf{v}^a \tag{6} \end{align} \]

 

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

 

\[ \begin{align} \Delta \mathbf{v}^a=J^{-1} \Delta \mathbf{B}^s \tag{7} \end{align} \]

 

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

 

\[ \begin{align} \mathbf{v}^a & \gets \mathbf{v}^a+ \Delta \mathbf{v}^a \tag{8} \end{align} \]

 

다음 메트랩 코드는 식 (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{align} r_p=2737,\ \ \ e=1.1,\ \ \ i=30^0, \ \ \ \Omega=45^0, \ \ \ \omega=60^0, \ \ \ \theta = -146.0905^0 \tag{9} \end{align} \]

 

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

 

\[ \begin{align} \mathbf{r}_{ref}^a= \begin{bmatrix} 3,503.47 \\ -37,139.74 \\ -32,922.45 \end{bmatrix} \ km, \ \ \ \ \ \mathbf{v}_{ref}^a=\begin{bmatrix} -0.32022 \\ 0.37905 \\ 0.28548 \end{bmatrix} \ km / s \tag{10} \end{align} \]

 

여기서 위첨자 \(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) 벡터 \(\vec{v}_\infty\) 와 최근접 도달 시간(TCA) 및 \(B_T, \ B_R\) 은 다음과 같이 계산된다.

 

\[ \begin{align} & \mathbf{v}_\infty^a= \begin{bmatrix} -0.22998 \\ 0.28610 \\ 0.21069 \end{bmatrix} \ km/s, \ \ \ \ \ t_{TCA}=92,554.66 \ sec \tag{11} \\ \\ & B_T=12,524.17 \ km, \ \ \ \ \ B_R=-677.97 \ km \end{align} \]

 

[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}_\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 100 \ km\) 의 무작위 섭동량을 가하고, 기준 속도벡터에는 \( \pm 0.02 \ km/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{align} \Delta B_T=2071.28 \ km, \ \ \ \ \ \Delta B_R=573.34 \ km, \ \ \ \ \ \Delta e=0.015 \tag{12} \end{align} \]

 

슈팅방법을 이용하면 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{align} \Delta \mathbf{v}^a= \begin{bmatrix} 0.0206 \\ -0.0179 \\ -0.0099 \end{bmatrix} \ km/s \tag{13} \end{align} \]

 

참고로 동일한 조건에서 식 (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-평면 타겟팅 알고리즘이다. 여기서 \(\vec{r}\) 과 \(\vec{v}\) 는 유주비행체의 현재 위치와 속도벡터, \(\vec{v}_\infty^\star, \ i^\star, \ r_p^\star\) 는 각각 쌍곡선 초과속도 벡터, 경사각, 근지점 거리의 설계값이다.

0. 입력: \(\vec{v}_\infty^\star, \ i^\star, \ r_p^\star\) 및 \(\vec{r}, \ \vec{v}\)

1. \(\vec{v}_\infty^\star, \ i^\star, \ r_p^\star\) 부터 \(B_T^\star, \ B_R^\star, \ e^\star\) 를 계산한다.

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

3. \( \Delta B_T \ \Delta B_R \ \Delta e \) 를 계산하고 \(\vec{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);