2.1 문제 정의

본 장에서는 1단 고체 추진기관이 탑재된 중단거리 유도탄의 교전 시나리오를 기반으로 충돌각 제어 궤적 최적화 문제를 정의한다. 그림. 1은 유도탄 발사부터 표적 요격까지의 비행궤적을 나타낸다.

여기서 $x _{0} , y _{0} ,V _{0} , \gamma _{0}$는 각각 발사 시점에서의 거리, 고도, 속도, 경로각을 나타내며, $x _{f} , y _{f} , V _{f} , \gamma _{f}$는 표적 요격 시점의 거리, 고도, 속도, 경로각을 나타낸다. 또한, $y _{max}$는 유도탄 운용 상에서 허용된 최대고도를 나타낸다. 이와 같은 유도탄의 비행궤적은 식(1)에 기술된 유도탄의 종방향 운동방정식에 의하여 결정된다.

여기서, $x, y, V, \gamma, \alpha, \alpha _{c}$는 각각 유도탄의 거리, 고도, 속도, 경로각, 받음각, 받음각 명령을 나타낸다. 또한, $m, g, \tau _{c}$는 유도탄의 질량과 중력가속도, 받음각 명령에 대한 시정수를 나타내고, $T, L, D$는 유도탄에 작용하는 추력, 양력, 항력을 나타낸다. 본 연구에서는 추진 구간에서 일정 크기의 추력 $T _{0}$가 발생하고, 연소구간 동안 유도탄의 무게가 일정하게 감소한다고 가정하였다. 그리고, 유도탄에서 발생하는 양력과 항력은 다음과 같이 받음각의 1차 및 2차 함수로 근사화하였다.

여기서, $\rho, S _{ref}$는 공기 밀도와 유도탄의 기준면적을 나타내고, $C _{L _{0}} , C _{L _{\alpha 1}} , C _{D _{0}} , C _{D _{\alpha 1}} , C _{D _{\alpha 2}}$는 양력 및 항력과 관련된 공력계수들을 나타낸다. 만일 비행 중 발생하는 유도탄의 받음각이 작다고 가정하고, 유도탄의 상태변수 및 제어변수를 각각 $p=[x, y, V, \gamma, \alpha] ^{T}$와 $u= \alpha _{c}$로 선정하면, 추진-활공 유도탄의 비선형 운동방정식은 다음과 같이 기술된다. 단, 상태변수의 개수를 최소화하고 최적화 문제 설정을 용이하게 하도록, 유도탄의 질량은 상태변수로 선정하지 않았다.

2.2 최적화 문제 설정

본 연구에서 고려하는 충돌각 제어 최적화 문제는 고도 및 받음각 구속조건을 만족하며 유도탄의 종말 경로각을 만족하는 최소에너지 비행궤적을 산출하는 문제이다. P0는 2.1절에서 기술한 식들을 기반으로 한 비선형 최적제어 문제를 나타낸다.

P0 : $\min \frac { 1 } { 2 } \int _ { t _ { 0 } } ^ { t _ { f } } \alpha _ { c } ^ { 2 } d t$

• dynamics

$\dot { p } = f _ { p } ( p ) + B _ { p } ( p ) u$

• state constraint

$0 \leq y \leq y _ { \max }$

• control constraint

$- \alpha _ { \max } \leq u \leq \alpha _ { \max }$

• initial conditions

$p \left( t _ { 0 } \right) = \left[ x \left( t _ { 0 } \right) , y \left( t _ { 0 } \right) , V \left( t _ { 0 } \right) , \gamma \left( t _ { 0 } \right) , \alpha \left( t _ { 0 } \right) \right] ^ { T }$

• final conditions

$x \left( t _ { f } \right) = x _ { f } , y \left( t _ { f } \right) = y _ { f } , \quad \gamma \left( t _ { f } \right) = \gamma _ { f }$

위와 같은 최적제어 문제를 수치적으로 해결하기 위해서는 단조 증가하거나 단조 감소하고, 경계값이 정해진 변수를 독립변수로 설정해야 한다. 본 연구에서는 고려하는 문제는 최종 시간이 고정되지 않았기 때문에, 시간 대신 거리를 독립변수로 선정하고 다음 식을 이용하여 운동방정식을 변환하였다.

이에 따라 유도탄의 상태변수는 $q = [ y , V , \gamma , \alpha ] ^ { T }$로 재선정되고, 식(3)의 운동방정식은 식(5)의 운동방정식으로 변환된다. 단, 유도탄의 추력과 질량은 컨벡스 프로그래밍을 위한 이산화를 수행하며 각각 식(9) 및 식(10)과 같이 거리에 대한 형태로 근사화하여 변환된다.

독립변수 변경에 따라 P0은 다음과 같이 P1로 변경된다.

P1 : $\min \frac { 1 } { 2 } \int _ { x _ { 0 } } ^ { x _ { f } } \left( \frac { \alpha _ { c } ^ { 2 } } { V \cos \gamma } \right) d x$

• dynamics

$q ^ { \prime } = f _ { q } ( q ) + B _ { q } ( q ) u$

• state constraint

$0 \leq y \leq y _ { \max }$

• control constraint

$- \alpha _ { \max } \leq u \leq \alpha _ { \max }$

• initial conditions

$q \left( x _ { 0 } \right) = \left[ y \left( x _ { 0 } \right) , V \left( x _ { 0 } \right) , \gamma \left( x _ { 0 } \right) , \alpha \left( x _ { 0 } \right) \right] ^ { T }$

• final conditions

$y \left( x _ { f } \right) = y _ { f } , \quad \gamma \left( x _ { f } \right) = \gamma _ { f }$

3.1 운동방정식 선형화 및 신뢰구간 설정

최적화 문제에서 $k$번째로 산출된 결과를 $q ^ {k}$라고 가정하자. 운동방정식 (5)는 부분 선형화를 기반으로 다음과 같이 선형화할 수 있다.

여기서 $a _ {ij}$, $b _ {ij}$는 $A$행렬과 $B$행렬의 $i$번째 행과 $j$번째 열의 값을 나타낸다. 부분 선형화 기법은 선형화를 수행하면서 제어 입력에 대한 정보가 요구되지 않기 때문에, 진동 특성을 억제하고 수렴 속도를 높이는 것으로 알려져 있다. 그리고 선형화는 하나의 평형점에 대한 근사화 결과를 산출하므로 다음과 같은 신뢰 구간을 추가적으로 설정하여 타당성을 확보하는 것이 필요하다.

3.2 운동방정식 이산화

본 연구에서 수행하는 컨벡스 최적화 문제는 독립변수인 $x$를 초기값($x _ {0}$)과 최종값($x _ {f}$) 사이를 동일한 $N$개의 구간으로 나누어 최적화 결과를 도출한다. 운동방정식은 미분방정식으로 적절한 적분 기법을 사용하여 각 이산화 지점에서의 산술적인 구속조건 형태로 변환해야 한다. 사다리꼴 기법을 기반으로 이산화를 수행하면 다음과 같다⁽¹⁶⁾.

유도탄의 질량과 추력은 연소 종료 시점인 $t _ {T}$를 기준으로 값이 결정되기 때문에, 다음과 같이 $t _ {T}$ 시점을 포함하는 $i _ {T}$번째 구간$\left[ x _ { i _ { T } - 1 } , x _ { i _ { T } } \right]$)을 산출하는 것이 필요하다.

$i _ { T }$번째 구간을 기준으로 유도탄의 질량과 추력은 다음과 같이 기술된다.

3.3 컨벡스 최적화 문제 정의

지금까지 기술된 결과를 바탕으로 비선형 최적제어 문제 P1은 다음과 같은 컨벡스 최적화 문제 P2로 변환된다.

P2 : $\min \frac { 1 } { 2 } \sum _ { i = 1 } ^ { N } k _ { i } \eta _ { i }$ where $\kappa _ { i } = \frac { 2 \Delta x } { \left( V _ { i - 1 } \cos \gamma _ { i - 1 } + V _ { i } \cos \gamma _ { i } \right) }$

• dynamics

$H _ { q } = F$

• state constraint

$0 \leq y _ { i } \leq y _ { \max }$ $\left\| q - q ^ { k } \right\| \leq r _ { q }$

• control constraint

$- \alpha _ { \max } \leq u _ { i } \leq \alpha _ { \max }$ $u _ { i } ^ { 2 } \leq \eta _ { i }$

• initial conditions

$q \left( x _ { 0 } \right) = \left[ y \left( x _ { 0 } \right) , V \left( x _ { 0 } \right) , \gamma \left( x _ { 0 } \right) , \alpha \left( x _ { 0 } \right) \right] ^ { T }$

• final conditions

$y \left( x _ { N } \right) = y _ { f } , \quad \gamma \left( x _ { N } \right) = \gamma _ { f }$

여기서 $\eta $는 2차 원뿔 프로그래밍(Second-Order Cone Programming) 형태로 기술하기 위하여 추가된 여유변수(slack variable)을 나타낸다. $\kappa _{i}$는 각 구간에서의 입력에 대한 가중치로 이전 단계에서 산출된 최적해로부터 계산된 값을 사용한다.

4.1 순차 컨벡스 프로그래밍 결과

다음 그림. 2은 허용 고도값 1000m에 대한 순차 컨벡스 프로그래밍의 결과를 보여준다. 총 10번의 연속적인 최적화를 수행한 후, 충돌각 및 고도 제한조건을 만족하는 수렴해를 도출하였다. 특히, 초기 비행궤적이 최종 비행궤적과 매우 상이함에도 불구하고 빠르게 수렴하는 것을 확인할 수 있다. 이와 같이 순차 컨벡스 프로그래밍은 초기 값에 매우 강건한 특성을 지니고 있기 때문에, 비선형 프로그래밍보다 빠르게 수렴 해를 도출하는 장점을 지니고 있다.

그림. 3은 속도, 경로각, 받음각에 대한 초기 및 최종 상태 결과를 나타낸다. 속도의 경우에는 추진-활공 구간을 고려하지 않고 초기값과 최종값을 선형적으로 연결하는 초기 상태 예측이 주어졌음에도 최종 결과에는 추진-활공 구간이 반영되어 산출되었음을 확인할 수 있다. 또한, 경로각은 초기 발사 시점과 최종 요격 시점에서의 값을 유지한 채 최적화가 수행되었음을 확인할 수 있다. 마지막으로 받음각은 전 구간에서 3도 이내로 발생하며 초기 구간에서 다른 구간보다 크게 발생함을 확인할 수 있다.

다음 그림. 4은 허용 고도 변화에 따른 유도탄의 비행궤적을 나타낸다. 허용 고도가 1200m일 경우에는 제한 값보다 작은 약 1100m 정도까지 고도가 상승하지만, 허용 고도가 1000m, 800m, 600m로 감소한 경우에는 허용 고도까지 상승한 후 표적을 요격하는 결과를 확인할 수 있다.

다음 표 2는 허용 고도 변화에 따른 목적함수 값의 변화를 나타낸다. 고도 제한 값까지 상승하지 않는 1400m, 1200m에서는 동일한 값을 나타내지만, 고도 제한으로 인하여 비행궤적이 변하는 경우에는 그에 따라 최소에너지가 증가하고 있음을 확인할 수 있다.

4.2 결과 검증

본 절에서는 순차 컨벡스 프로그래밍으로 산출된 결과의 타당성을 확인하기 위하여 동일한 문제에 대하여 비선형 프로그래밍(GPOPS-II) 및 참고문헌 [1]⁽¹⁾의 최적제어 기반의 유도법칙으로 산출된 결과를 비교하였다⁽¹⁾. 단, 최적제어 기반의 유도법칙은 조종지연을 고려하지 않고 유도한 결과를 사용하였다. 또한, 고도에 대한 구속조건을 배제한 결과를 비교하기 위하여 컨벡스 및 비선형 프로그래밍의 허용 고도 제한을 1200m로 설정하여 고도 제한에 걸리지 않도록 하였다. 다음 그림. 5는 세 가지 방법으로 산출된 결과를 나타낸다.

세 가지 다른 방법으로 산출된 비행궤적은 전반적으로 유사한 결과를 도출하지만, 최대고도 측면에서 약간의 차이가 발생함을 확인할 수 있다. 특히, 최적제어 기반의 유도법칙은 수치적으로 산출된 다른 기법들에 비하여 속도가 낮은 초기 구간에서의 명령이 크게 발생하고 그에 따라 고도가 상승함을 볼 수 있다. 표 3을 살펴보면 수치 기반의 최적화 기법들은 유도탄의 속도 변화에 대한 부분이 고려되어 일정 속도를 가정하여 산출된 최적제어 기반의 유도기법에 비하여 에너지가 적게 소모됨을 확인할 수 있다.

다음 표 4는 컨벡스 프로그래밍과 비선형 프로그래밍의 계산시간을 비교한 결과를 나타내는데, 컨벡스 프로그래밍이 비선형 프로그래밍에 비하여 빠르고 안정적으로 해에 수렴하는 것을 확인할 수 있다. 컨벡스 프로그래밍은 MOSEK을 사용하여 수행하였고, 코드 최적화를 수행하지 않았음에도 10회의 최적화에 대하여 약 6.068초의 평균 계산시간이 소요되었다. 반면, GPOPS-II의 경우 tolerance 값과 최적화 방식에 따라 계산시간에 차이가 크게 발생하였는데, tolerance 0.001에 대하여 SNOPT와 IPOPT을 사용한 결과 각각 17.34초와 983.83초가 소요되었다. (여기서, mesh는 hp로 적용하였다.) 결과에서 나타난 것과 같이 GPOPS-II는 문제에 따른 설정 값을 적절하게 선정하는 것이 계산시간에 매우 중요하게 작용하는 반면, 컨벡스 프로그래밍의 경우에는 초기 상태변수 예측 이외의 다른 설정 값에 크게 영향을 받지 않는 장점이 있다.

다음 그림. 6은 허용 고도 제한에 따른 순차 컨벡스 프로그래밍과 비선형 프로그래밍을 비교한 결과를 나타낸다. 허용 고도는 1000m로 설정한 후, 두 최적화 결과를 산출하였다. 허용 고도 구속조건을 만족하는 유사한 비행궤적을 나타냄을 확인할 수 있다.