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

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

 

$$\begin{array}{r}\text{(1)}& & \dot{\mathbf{r}}=\mathbf{v}\\ & \dot{\mathbf{v}}=\mathbf{g}(\mathbf{r})+\mathbf{a}\end{array}$$

 

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

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

 

$$\begin{array}{r}\text{(2)}& & J=\frac{1}{2}{\int }_{t_0}^{t_f}{\mathbf{a}}^T\mathbf{a}dt\\ \text{(3)}& & \mathbf{r}(t_0)={\mathbf{r}}_0,\mathbf{r}(t_f)={\mathbf{r}}_f\\ & \mathbf{v}(t_0)={\mathbf{v}}_0\end{array}$$

 

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

 

 

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

 

 

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

 

$$\begin{array}{}\text{(4)}& \dot{\mathbf{x}}=A\mathbf{x}+B(\mathbf{a}+\mathbf{g})\end{array}$$

여기서

$$\begin{array}{rl}& A=\left[\begin{array}{cc}0& I\\ 0& 0\end{array}\right],B=\left[\begin{array}{c}0\\ I\end{array}\right]\\ \\ & \mathbf{x}=\left[\begin{array}{c}\mathbf{r}\\ \mathbf{v}\end{array}\right],\mathbf{x}(t_0)=\left[\begin{array}{c}{\mathbf{r}}_0\\ {\mathbf{v}}_0\end{array}\right]\end{array}$$

 

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

 

$$\begin{array}{r}\text{(5)}& & C\mathbf{x}(t_f)={\mathbf{r}}_f\\ & C=\left[\begin{array}{cc}I& 0\end{array}\right]\end{array}$$

 

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

 

$$\begin{array}{}\text{(6)}& \mathcal{H}=\frac{1}{2}({\mathbf{a}}^T\mathbf{a})+{\lambda }^T(A\mathbf{x}+B(\mathbf{a}+\mathbf{g}))\end{array}$$

 

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

 

$$\begin{array}{r}\text{(7)}& \dot{\mathbf{x}}& =\frac{\partial \mathcal{H}}{\partial \lambda }=A\mathbf{x}+B(\mathbf{a}+\mathbf{g})\\ \text{(8)}& -\dot{\lambda }& =\frac{\partial \mathcal{H}}{\partial \mathbf{x}}=A^T\lambda \end{array}$$

 

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

 

$$\begin{array}{}\text{(9)}& \lambda (t)=e^{A^T(t_f-t)}\lambda (t_f)\end{array}$$

 

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

 

$$\begin{array}{}\text{(10)}& C^T\nu -\lambda (t_f)=0\end{array}$$

 

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

 

$$\begin{array}{}\text{(11)}& \lambda (t)=e^{A^T(t_f-t)}{C}^T\nu \end{array}$$

 

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

 

$$\begin{array}{}\text{(12)}& 0=\frac{\partial \mathcal{H}}{\partial \mathbf{a}}=B^T\lambda +\mathbf{a}\end{array}$$

 

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

 

$$\begin{array}{r}\text{(13)}& \mathbf{a}(t)& =-B^T\lambda (t)\\ & =-B^T{e}^{A^T(t_f-t)}{C}^T\nu \end{array}$$

 

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

 

$$\begin{array}{}\text{(14)}& \dot{\mathbf{x}}=A\mathbf{x}+B(\mathbf{g}-B^T{e}^{A^T(t_f-t)}{C}^T\nu )\end{array}$$

 

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

 

$$\begin{array}{r}\text{(15)}& \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}\\ & -{\int }_{t_0}^t{e}^{A(t-\tau )}B{B}^T{e}^{A^T(t_f-\tau )}d\tau C^T\nu \end{array}$$

 

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

 

$$\begin{array}{r}\text{(16)}& \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}\\ & -C{\int }_{t_0}^{t_f}{e}^{A(t_f-\tau )}B{B}^T{e}^{A^T(t_f-\tau )}d\tau C^T\nu \end{array}$$

 

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

 

$$\begin{array}{r}\text{(17)}& \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}\\ & -C{\int }_t^{t_f}{e}^{A(t_f-\tau )}B{B}^T{e}^{A^T(t_f-\tau )}d\tau C^T\nu \end{array}$$

 

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

 

$$\begin{array}{}\text{(18)}& P(t)=C{\int }_t^{t_f}{e}^{A(t_f-\tau )}B{B}^T{e}^{A^T(t_f-\tau )}d\tau C^T\end{array}$$

 

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

 

$$\begin{array}{}\text{(19)}& \nu =-P^{-1}(t)({\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))\end{array}$$

 

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

 

 

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

 

$$\begin{array}{}\text{(20)}& \mathbf{a}(t)=B^T{e}^{A^T(t_f-t)}{C}^T{P}^{-1}(t)({\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))\end{array}$$

 

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

 

$$\begin{array}{}\text{(21)}& e^{A(t_f-t)}=I+A(t_f-t)=\left[\begin{array}{cc}I& (t_f-t)I\\ 0& I\end{array}\right]\end{array}$$

 

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

 

$$\begin{array}{}\text{(22)}& P(t)=\frac{t_{go}^3}{3}I\end{array}$$

 

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

 

$$\begin{array}{r}\text{(23)}& & B^T{e}^{A^T(t_f-t)}{C}^T{P}^{-1}(t)=\frac{3}{t_{go}^3}I\\ \text{(24)}& & C{\int }_t^{t_f}{e}^{A(t_f-\tau )}d\tau B\mathbf{g}=\frac{t_{go}^2}{2}\mathbf{g}\\ \text{(25)}& & 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)\\ & ={\mathbf{r}}_f-\frac{t_{go}^2}{2}\mathbf{g}-\mathbf{r}(t)-t_{go}\mathbf{v}(t)\end{array}$$

 

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

 

$$\begin{array}{r}\text{(26)}& \mathbf{a}(t)& =\frac{3}{t_{go}^2}({\mathbf{r}}_f-(\mathbf{r}(t)+t_{go}\mathbf{v}(t)))-\frac{3}{2}\mathbf{g}\\ & =\frac{3}{t_{go}^2}({\mathbf{r}}_f-(\mathbf{r}(t)+t_{go}\mathbf{v}(t)+\frac{t_{go}^2}{2}\mathbf{g}))\\ \\ & =\frac{3}{t_{go}^2}ZEM\end{array}$$

 

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

 

$$\begin{array}{}\text{(27)}& ZEM={\mathbf{r}}_f-(\mathbf{r}(t)+t_{go}\mathbf{v}(t)+\frac{t_{go}^2}{2}\mathbf{g})\end{array}$$

 

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

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

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

 

$$\begin{array}{r}\text{(28)}& \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}\\ & -{\int }_t^{t_f}{e}^{A(t_f-\tau )}B{B}^T{e}^{A^T(t_f-\tau )}d\tau C^T\nu \\ \\ & =\left[\begin{array}{cc}I& t_{go}I\\ 0& I\end{array}\right]\mathbf{x}(t)+\left[\begin{array}{c}\frac{t_{go}^2}{2}\mathbf{g}\\ t_{go}\mathbf{g}\end{array}\right]\\ \\ & +\left[\begin{array}{cc}\frac{t_{go}^3}{3}I& \frac{t_{go}^2}{2}I\\ \frac{t_{go}^2}{2}I& t_{go}I\end{array}\right]\left[\begin{array}{c}I\\ 0\end{array}\right]\frac{3}{t_{go}^3}I({\mathbf{r}}_f-\frac{t_{go}^2}{2}\mathbf{g}-\mathbf{r}(t)-t_{go}\mathbf{v}(t))\end{array}$$

 

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

 

$$\begin{array}{r}\text{(29)}& \mathbf{r}(t_f)& =\mathbf{r}(t)+t_{go}\mathbf{v}(t)+\frac{t_{go}^2}{2}\mathbf{g}+({\mathbf{r}}_f-\frac{t_{go}^2}{2}\mathbf{g}-\mathbf{r}(t)-t_{go}\mathbf{v}(t))\\ & ={\mathbf{r}}_f\\ \\ \text{(30)}& \mathbf{v}(t_f)& =\mathbf{v}(t)+t_{go}\mathbf{g}+\frac{3}{2t_{go}}({\mathbf{r}}_f-\frac{t_{go}^2}{2}\mathbf{g}-\mathbf{r}(t)-t_{go}\mathbf{v}(t))\\ & =\frac{3}{2t_{go}}({\mathbf{r}}_f-\mathbf{r}(t))-\frac{1}{2}\mathbf{v}(t)+\frac{1}{4}{t}_{go}\mathbf{g}\end{array}$$

 

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

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

 

$$\begin{array}{}\text{(31)}& \mathbf{a}(t)=\frac{3}{t_{go}^2}({\mathbf{r}}_f-(\mathbf{r}(t)+t_{go}\mathbf{v}(t)))-\mathbf{g}\end{array}$$

 

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

 

$$\begin{array}{r}\text{(32)}& \mathbf{r}(t_f)& =\mathbf{r}(t)+t_{go}\mathbf{v}(t)+\frac{t_{go}^2}{2}\mathbf{g}+({\mathbf{r}}_f-\frac{t_{go}^2}{3}\mathbf{g}-\mathbf{r}(t)-t_{go}\mathbf{v}(t))\\ & ={\mathbf{r}}_f+\frac{t_{go}^2}{6}\mathbf{g}\end{array}$$

 

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