최적유도법칙 (Optimal Guidance Law): 최종 속도 미설정

이전 포스트(https://pasus.tistory.com/293)와 유사한 문제를 풀어본다. 차이점은 최종 시간에서 속도벡터에 관한 제약조건이 없는 경우이다. 편의상 운동 방정식을 다시 쓴다.

 

\[ \begin{align} & \dot{\mathbf{r}}= \mathbf{v} \tag{1} \\ \\ & \dot{\mathbf{v}}= \mathbf{g}(\mathbf{r})+\mathbf{a} \end{align} \]

 

여기서 \(\mathbf{r}\) 과 \(\mathbf{v}\) 는 각각 관성좌표계에 대한 위치벡터와 속도벡터를 나타낸다. \(\mathbf{a}\) 는 제어 가속도, \(\mathbf{g}(\mathbf{r})\) 은 비행체 또는 미사일에 작용하는 중력 가속도로서 위치의 함수이다.

비용함수와 제약조건은 다음과 같다.

 

\[ \begin{align} & J=\frac{1}{2} \int_{t_0}^{t_f} \mathbf{a}^T \mathbf{a} \ dt \tag{2} \\ \\ & \ \ \mathbf{r}(t_0 )=\mathbf{r}_0, \ \ \ \ \mathbf{r}(t_f )=\mathbf{r}_f \tag{3} \\ \\ & \ \ \mathbf{v}(t_0 )= \mathbf{v}_0 \end{align} \]

 

랑데부가 아닌 요격 임무의 경우에는 최종 속도에 관한 제약조건을 두지 않는 경우가 대부분으로서 식 (3)과 같이 최종 속도는 자유(free)로 둔다.

 

 

식 (1)과 (2), (3)으로 주어지는 최적제어 문제는 예제(https://pasus.tistory.com/257, https://pasus.tistory.com/260)에서 풀었던 최종상태가 부분적으로 설정된 LQR 문제와 거의 동일하지만, 여기서 다시 자세히 풀어보도록 한다.

 

 

이전 포스트와 마찬가지로 식 (1)에서 중력 가속도를 상수벡터로 가정하겠다. 그러면 식 (1), (3)은 다음과 같이 상태변수 방정식으로 표현할 수 있다.

 

\[ \dot{\mathbf{x}} =A \mathbf{x}+B( \mathbf{a}+\mathbf{g}) \tag{4} \]

여기서

\[ \begin{align} & A= \begin{bmatrix} 0 & I \\ 0 & 0 \end{bmatrix} , \ \ \ \ B=\begin{bmatrix} 0 \\ I \end{bmatrix} \\ \\ & \mathbf{x}= \begin{bmatrix} \mathbf{r} \\ \mathbf{v} \end{bmatrix} , \ \ \ \ \mathbf{x}(t_0 )= \begin{bmatrix} \mathbf{r}_0 \\ \mathbf{v}_0 \end{bmatrix} \end{align} \]

 

이다. 최종 상태변수의 제약조건은 다음과 같다.

 

\[ \begin{align} & C \mathbf{x}(t_f )= \mathbf{r}_f \tag{5} \\ \\ & C = \begin{bmatrix} I & 0 \end{bmatrix} \end{align} \]

 

이 문제를 풀기 위해 먼저 해밀토니안(Hamiltonian) 함수를 정의한다.

 

\[ \mathcal{H}= \frac{1}{2} (\mathbf{a}^T \mathbf{a})+ \lambda^T \left(A \mathbf{x}+B( \mathbf{a}+\mathbf{g}) \right) \tag{6} \]

 

그러면 다음과 같이 상태변수와 코스테이트(costate) 미분 방정식을 얻을 수 있다.

 

\[ \begin{align} \dot{\mathbf{x}} &= \frac{\partial \mathcal{H}}{ \partial \lambda} =A \mathbf{x}+B(\mathbf{a}+\mathbf{g}) \tag{7} \\ \\ -\dot{\lambda} &= \frac{\partial \mathcal{H}}{ \partial \mathbf{x}}=A^T \lambda \tag{8} \end{align} \]

 

코스테이트 미분 방정식은 상태변수와 독립이므로 다음과 같이 쉽게 해를 구할 수 있다.

 

\[ \lambda (t)=e^{A^T (t_f-t) } \lambda (t_f) \tag{9} \]

 

상태변수와 코스테이트 미분 방정식을 풀기 위해서는 경계조건이 필요하다. 여기서는 초기 시간과 초기 상태변수, 그리고 최종 시간이 정해졌으므로 \(dt_0=0, \ dx(t_0 )=0, \ dt_f=0\) 이 된다. 따라서 최적제어의 필요조건을 정리한 표(https://pasus.tistory.com/258)의 경계조건에 의하면,

 

\[ C^T \nu -\lambda (t_f )=0 \tag{10} \]

 

를 얻을 수 있다. 여기서 \(\nu\) 는 라그랑지 곱수이다. 식 (10)을 (9)에 대입하면,

 

\[ \lambda (t)=e^{A^T (t_f-t) } C^T \nu \tag{11} \]

 

이 된다. 정정조건(stationary condition)은 다음과 같다.

 

\[ 0= \frac{\partial \mathcal{H}}{\partial \mathbf{a}} =B^T \lambda + \mathbf{a} \tag{12} \]

 

따라서 최적제어는 식 (12)와 (11)로부터 다음과 같이 계산할 수 있다.

 

\[ \begin{align} \mathbf{a}(t) &= -B^T \lambda (t) \tag{13} \\ \\ &= -B^T e^{A^T (t_f-t) } C^T \nu \end{align} \]

 

식 (13)을 식 (7)에 대입하면 상태변수 미분 방정식은 다음과 같이 된다.

 

\[ \dot{\mathbf{x}}=A\mathbf{x}+B \left( \mathbf{g}-B^T e^{A^T (t_f-t) } C^T \nu \right) \tag{14} \]

 

위 식을 풀면 \(\mathbf{x}(t)\) 를 다음과 같이 계산할 수 있다.

 

\[ \begin{align} \mathbf{x}(t) &= e^{A(t-t_0 )} \mathbf{x}(t_0 )+ \int_{t_0}^t e^{A(t-\tau)} \ d \tau \ B \mathbf{g} \tag{15} \\ \\ & \ \ \ \ \ \ \ -\int_{t_0}^t e^{A(t-\tau)} BB^T e^{A^T (t_f-\tau) } \ d \tau \ C^T \nu \end{align} \]

 

위 식에 의하면 \(t=t_f\) 에서 \(\mathbf{r}(t_f )\) 를 다음과 같이 구할 수 있다.

 

\[ \begin{align} \mathbf{r}(t_f) &= C e^{A(t_f-t_0 )} \mathbf{x}(t_0 )+ C \int_{t_0}^{t_f} e^{A(t_f-\tau)} \ d \tau \ B \mathbf{g} \tag{16} \\ \\ & \ \ \ \ \ \ \ - C \int_{t_0}^{t_f} e^{A(t_f-\tau)} BB^T e^{A^T (t_f-\tau) } \ d \tau \ C^T \nu \end{align} \]

 

식 (16)는 시간 \(t_0\) 에서 상태변수 \(\mathbf{x}(t_0 )\) 가 주어졌을 때 \(\mathbf{r}(t_f )\) 를 계산하는 식이다. 이 식은 임의의 시간 \(t \le t_f\) 에서 상태변수 \(\mathbf{x}(t)\) 가 주어졌을 때에도 성립해야 한다. 따라서 식 (16)은 다음과 같이 수정할 수 있다.

 

\[ \begin{align} \mathbf{r}(t_f) &= C e^{A(t_f-t )} \mathbf{x}(t )+ C \int_{t}^{t_f} e^{A(t_f-\tau)} \ d \tau \ B \mathbf{g} \tag{17} \\ \\ & \ \ \ \ \ \ \ - C \int_{t}^{t_f} e^{A(t_f-\tau)} BB^T e^{A^T (t_f-\tau) } \ d \tau \ C^T \nu \end{align} \]

 

이제 이 값이 \(\mathbf{r}(t_f )= \mathbf{r}_f\) 가 되도록 \(\nu\) 를 계산해야 한다. 여기서 \(P(t)\) 를 다음과 같이 정의하면,

 

\[ P(t)= C \int_t^{t_f} e^{A(t_f-\tau)} BB^T e^{A^T (t_f-\tau) } \ d \tau \ C^T \tag{18} \]

 

\(\nu\) 를 다음과 같이 계산할 수 있다.

 

\[ \nu=-P^{-1} (t) \left( \mathbf{r}_f-C \int_t^{t_f} e^{A(t_f-\tau)} \ d\tau \ B \mathbf{g}-C e^{A(t_f-t)} \mathbf{x}(t) \right) \tag{19} \]

 

여기서 만약 \(P^{-1} (t)\) 이 존재하지 않는다면 이 최적제어 문제는 해가 존재하지 않는다.

 

 

이제 식 (19)를 식 (13)에 대입하면 최종적으로 최적제어를 계산할 수 있다.

 

\[ \mathbf{a}(t)=B^T e^{A^T (t_f-t) } C^T P^{-1} (t) \left( \mathbf{r}_f-C \int_t^{t_f} e^{A(t_f-\tau)} \ d\tau \ B \mathbf{g}-C e^{A(t_f-t)} \mathbf{x}(t) \right) \tag{20} \]

 

이제 식 (20)을 구성하는 행렬을 차례로 계산해 보자. 먼저 행렬지수함수는 다음과 같이 계산된다.

 

\[ e^{A(t_f-t)} =I+A(t_f-t)= \begin{bmatrix} I & (t_f-t)I \\ 0 & I \end{bmatrix} \tag{21} \]

 

행렬 \(P(t)\) 는 정의와 식 (18)에 의해서 다음과 같이 계산된다.

 

\[ P(t) = \frac{t_{go}^3}{3} I \tag{22} \]

 

여기서 \(t_{go}=t_f-t\) 는 잔여시간(time-to-go)을 나타낸다. 식 (21)과 (22)를 이용하면 식 (20)의 구성 부분을 다음과 같이 계산할 수 있다.

 

\[ \begin{align} & B^T e^{A^T (t_f-t) } C^T P^{-1} (t)= \frac{3}{t_{go}^3 } I \tag{23} \\ \\ & C \int_t^{t_f} e^{A(t_f-\tau)} \ d\tau \ B \mathbf{g}= \frac{t_{go}^2}{2} \mathbf{g} \tag{24} \\ \\ & r_f-C \int_t^{t_f} e^{A(t_f-\tau)} \ d\tau \ B \mathbf{g}-Ce^{A(t_f-t)} \mathbf{x}(t) \tag{25} \\ \\ & \ \ \ \ \ \ \ = \mathbf{r}_f- \frac{t_{go}^2}{2} \mathbf{g}-\mathbf{r}(t)-t_{go} \mathbf{v} (t) \end{align} \]

 

식 (23)~(25)를 식 (20)에 대입하면 최적제어는 다음과 같다.

 

\[ \begin{align} \mathbf{a}(t) &= \frac{3}{t_{go}^2} \left( \mathbf{r}_f- \left( \mathbf{r}(t)+t_{go} \mathbf{v}(t) \right) \right)-\frac{3}{2} \mathbf{g} \tag{26} \\ \\ &= \frac{3}{t_{go}^2 } \left( \mathbf{r}_f- \left( \mathbf{r}(t)+t_{go} \mathbf{v}(t)+ \frac{t_{go}^2}{2} \mathbf{g} \right) \right) \\ \\ &= \frac{3}{t_{go}^2} ZEM \end{align} \]

 

여기서 ZEM(Zero-Effort-Miss)는 다음과 같이 정의된다.

 

\[ ZEM = \mathbf{r}_f- \left( \mathbf{r}(t)+t_{go} \mathbf{v}(t)+ \frac{t_{go}^2}{2} \mathbf{g} \right) \tag{27} \]

 

정의에 의하면 ZEM는 현재 시간 \(t\) 에서 임무 종료 시간 \(t_f\) 까지 제어 가속도가 \(0\) 일 때, 즉 zero-effort 일 때 임무 종료 시에 예측되는 거리 오차를 뜻한다.

식 (26)을 최종 속도가 설정되지 않은 경우의 최적 피드백 유도법칙(optimal feedback guidance law)이라고 한다.

이제, 식 (26)을 식 (15)에 대입하면 다음과 같이 된다.

 

\[ \begin{align} \mathbf{x}(t_f ) &=e^{A(t_f-t)} \mathbf{x}(t)+ \int_t^{t_f} e^{A(t_f-\tau)} \ d\tau \ B \mathbf{g} \tag{28} \\ \\ & \ \ \ \ \ \ \ -\int_t^{t_f} e^{A(t_f-\tau)} BB^T e^{A^T (t_f-\tau) } \ d\tau \ C^T \nu \\ \\ &= \begin{bmatrix} I & t_{go} I \\ 0 & I \end{bmatrix} \mathbf{x}(t)+ \begin{bmatrix} \frac{t_{go}^2}{2} \mathbf{g} \\ t_{go} \mathbf{g} \end{bmatrix} \\ \\ & \ \ \ \ \ \ \ + \begin{bmatrix} \frac{t_{go}^3}{3} I & \frac{t_{go}^2}{2} I \\ \frac{t_{go}^2}{2} I & t_{go} I \end{bmatrix} \begin{bmatrix} I \\ 0 \end{bmatrix} \frac{3}{t_{go}^3 } I \left( \mathbf{r}_f-\frac{t_{go}^2}{2} \mathbf{g}-\mathbf{r}(t)-t_{go} \mathbf{v}(t) \right) \end{align} \]

 

위 식으로부터 최종 시간에서 위치벡터와 속도벡터를 구하면 다음과 같다.

 

\[ \begin{align} \mathbf{r}(t_f ) &= \mathbf{r}(t)+t_{go} \mathbf{v}(t)+ \frac{t_{go}^2}{2} \mathbf{g}+ \left( \mathbf{r}_f-\frac{t_{go}^2}{2} \mathbf{g}-\mathbf{r}(t)-t_{go} \mathbf{v}(t) \right) \tag{29} \\ \\ &= \mathbf{r}_f \\ \\ \mathbf{v}(t_f ) &= \mathbf{v}(t)+t_{go} \mathbf{g}+ \frac{3}{2t_{go} } \left( \mathbf{r}_f- \frac{t_{go}^2}{2} \mathbf{g}-\mathbf{r}(t)-t_{go} \mathbf{v}(t) \right) \tag{30} \\ \\ &= \frac{3}{2t_{go}} (\mathbf{r}_f-\mathbf{r}(t))-\frac{1}{2} \mathbf{v}(t)+ \frac{1}{4} t_{go} \mathbf{g} \end{align} \]

 

계산에서 보듯이 식 (26)으로 주어지는 최적제어를 사용하면 \(\mathbf{r}(t_f )=\mathbf{r}_f\) 를 만족한다.

식 (26) 대신에 식 (31)의 제어를 사용하면 어떨까. 직관적으로는 식 (31)도 타당할 것 같다.

 

\[ \mathbf{a}(t)= \frac{3}{t_{go}^2 } \left( \mathbf{r}_f- \left( \mathbf{r}(t)+t_{go} \mathbf{v}(t) \right) \right) - \mathbf{g} \tag{31} \]

 

이 경우는 현재 시간 \(t\) 에서 예상되는 최종 위치벡터는 다음과 같이 계산된다.

 

\[ \begin{align} \mathbf{r}(t_f ) &= \mathbf{r}(t)+t_{go} \mathbf{v}(t)+ \frac{t_{go}^2}{2} \mathbf{g}+ \left( \mathbf{r}_f-\frac{t_{go}^2}{3} \mathbf{g}-\mathbf{r}(t)-t_{go} \mathbf{v}(t) \right) \tag{32} \\ \\ &= \mathbf{r}_f + \frac{t_{go}^2}{6} \mathbf{g} \end{align} \]

 

최종 임무 시간에 접근함에 따라 \(t_{go} \to 0\) 이 되므로 이 경우도 점근적으로 \(\mathbf{r}(t_f )=\mathbf{r}_f\) 이 된다.