2.1 운동방정식

이 절에서는 수직축 2자유도 운동 모델로 가정한 지능화 탄약의 상태를 지배하는 운동방정식을 정의한다. 지능화 탄약은 운용개념 상 포발사 후 일정 조건에 이르면 조종날개를 전개하여 유도비행하는 상태천이모델이기 때문에, 운동방정식을 다음과 같이 두 단계로 표현할 수 있다.

여기서 D₁[N], D₂[N], L[N]은 각각 항력과 양력, m[kg]은 탄의 비행 중 질량, g[m/s²]는 중력가속도, v[m/s]는 비행속도, γ[rad]는 비행경로각, r[m], h[m]는 각각 비행 사거리와 고도이다. 조종날개 전개 전에는 운동방정식 (1)에 따라 항력 D₁과 중력에 영향을 받아 상태를 변화시키며 비행한 후, 조종날개를 전개 및 구동하면 운동방정식 (2)에 따라 항력 D₂와 양력 L이 운동방정식의 입력으로서 작용하여 유도명령을 적용하여 비행하는 시스템이다.

2.2 공력 모델링

(1), (2)와 같이 표현된 운동방정식의 입력에 해당하는 D₁, D₂, L의 공력 모델은 다음과 같다. 단, 탄의 형상이 횡단면을 기준으로 대칭인 점을 고려하여 다음과 같이 모델링하였다.

이때, S는 공력환산계수, ρ는 공기밀도로서 고도 h에 의존하는 함수이며, ρ₀와 K_ρ는 대기권의 고도 구간에 의존하는 상수이다. 이후 사거리 최대화 최적유도를 수행함에 있어서, 지능화 탄약의 유도비행 시에는 트림점(trim) 부근에서 비행하는 것을 가정하기 때문에 비행 상태에 따라 비행기의 수직축 유도명령을 제한할 필요가 있다. 이때, 다음 장에서 후술하는 최적화문제에서 유도명령 제한 구속조건을 구성할 때 상태변수를 포함하지 않는 단순한 형태를 취하기 위해서는 수직축 유도명령, 즉 운동방정식의 입력(u_L)을 양력 계수(C_L)로 정하는 것이 유용하다. 이러한 모델링을 통해 C_L=u_L이 운동방정식 (2)에 작용하는 유일한 입력변수라는 것을 알 수 있으며, 이 양력계수 명령이 이 연구에서 최적화할 변수에 해당한다.

C_D10, C_D20는 공기의 물리적 성질을 고려했을 때 마하수에 의존하는 함수이며, 이 연구에 사용된 C_D10, C_D20의 공력 데이터는 그림. 1의 ○와 △가 나타내는 것과 같다. 고속 비행체의 항력 모델링에 관한 Bridges의 연구에 따르면, C_D20와 같이 마하 1을 경계선으로 step 함수와 유사한 형태로 변화하는 항력을 다음과 같이 모델링할 수 있다⁽¹⁵⁾.

단, M은 마하수, a₂₀, b₂₀, c₂₀, d₂₀는 C_D20의 항력 모델링 계수이다. 또한, C_D10는 마하 1을 기점으로 항력계수의 값이 상승한 후 마하수가 증가함에 따라 점차적으로 감소하는 형태로, Bridges의 연구결과를 참고하여 다음과 같이 모델링하였다.

이때, a₁₀, b₁₀, c₁₀, d₁₀, fvd₁₀, ed₁₀, nd₁₀은 C_D10항력 모델링 계수이다. C_D10, C_D20와 마찬가지로 K_D 또한 공력의 성질을 엄밀하게 고려한다면 마하수에 관한 함수로 표현되어야 하지만, C_D20+K_DC_L²에 관한 공력 데이터를 분석한 결과, 트림점이 존재하는 범위 내에서는 상수로 지정하여도 이 연구의 목적을 저해하지 않는 수준으로 공력 데이터를 모사할 수 있다고 판단하여, 마하수에 무관한 상수로 정의하였다.

운동방정식에서 마하수 M은 속도 v와 고도 h의 함수로서, 음속 a(h)의 모델링을 이용하여 다음과 같이 나타낼 수 있다.

단, $T_a=\left\{\begin{array}{ll}{T}_0+\lambda h,& \text{ when }h\leqq 11000\\ T_0+11000\lambda ,& \text{ when }h>11000\end{array}\right.$

마지막으로 운동방정식의 상태변수 x를 x=[v γ r h]^T로 정의하고 양력계수 C_L을 운동방정식에 대한 입력 u_L로 표현하면, 공력 모델링 (3)을 이용하여 탄의 운동방정식 식(1), 식(2)를 각각 다음과 같이 다시 기술할 수 있다.

단, $K_a=\frac{1}{2}{\rho }_{0}S$

3.1 최적화 문제의 정의

2장에 기술한 것과 같이 지능화 탄약의 운동방정식은 2가지로 표현되며, 이는 사거리 최적화를 수행하기 위해 다점 경계치문제(Multi-Point Boundary Value Problem)를 풀어야 한다는 것을 의미한다. 그러나 식 (7)과 같이 조종날개 전개 전의 운동방정식 f₁은 입력 u_L에 의존하지 않기 때문에, 온전히 초기 상태변수와 시간의 함수(단, 시간에 대해 non-explicit)이다. 따라서 조종날개 전개 직전과 직후의 상태변수에 대하여 최적화 문제의 경계조건을 조종날개 전개시각 t_i에 대하여 정의하면, 조종날개 전개 후의 운동방정식 f₂의 관점으로 최적화 문제를 다음과 같은 Mayer 형태로 기술할 수 있다.

여기서, $t _ { 0 } = 0$ 은 포발사 시각이며, $t _ { i }$는 조종날개의 전개시각, $t _ { f }$는 탄착시각으로 $t _ { i }$, $t _ { f }$는 $U _ { L }$과 같은 최적화 변수에 해당한다. 이 문제에서는, 시각 $t _ { i }$에서의 경계조건으로서 $f _ { i }$에 관한 시간의 함수로 부가함으로써 다점 경계치 문제를 2점 경계치문제로 간주할 수 있다는 점에 주목하자. 그런데 입력변수 $u _ { L }$에 공력학적인 한계를 고려하여 부등식 구속조건 $\left| u _ { L } \right| \leqq u _ { L \max }$을 부가해야 한다. 이러한 경우, 폰트랴긴의 최대값 원칙의 접근방식을 사용할 수도 있으나, 부등식 구속조건을 $u _ { L }$에 대한 새로운 입력 $u _ { d }$를 정의하여 다음과 같은 등식 구속조건으로 변환하는 방법을 사용할 경우 $u _ { L }$값에 따른 경우의 수를 고려하지 않아도 되며, 최적화 프로그래밍의 범용성을 높일 수 있다⁽¹⁶⁾.

이때, $u = \left[ u _ { L } , u _ { d } \right] ^ { \top }$로 정의하였고, $u _ { Lmax } $는 조종날개 전개 시 탄약의 속도를 고려하였을 때 마하수 $M (v,h)$에 무관한 상수로 결정하였다. 따라서 운동방정식의 입력에 대한 구속조건 (10)은 상태변수와 무관계한 오로지 입력 $u$만의 함수이다.

결론적으로 목적함수 $J$에 대하여 입력에 대한 구속조건 $C(u)$, 경계조건 $\chi \left( x _ { i } \right)$, $\psi \left( x _ { f } \right)$, 그리고 운동방정식 $\dot { x } = f _ { 2 }$를 병합하여 수정된 목적함수 $\overline { J }$는 다음과 같이 표현할 수 있다.

단, $H = \lambda ^ { \top } f _ { 2 } + \mu ^ { \top } C$

여기서 $\mathcal { V }$, $\eta$, $\mu$, $\lambda$는 라그랑쥬 미정계수이며, $\mathcal { V }$는 스칼라 상수, $\eta$는 (4x1)차원 벡터 상수, $\mu$, $\lambda$는 시간에 대한 (4x1)차원 벡터함수이다.

3.2 사거리 최대화 문제의 수치적 해법

(8)과 같이 표현된 지능화 탄약의 운동방정식은 비선형 시스템이기 때문에 (9)로 표현된 최적화 문제의 해를 해석적으로 구하기가 어려우므로, 이 연구에서는 수치적 해법을 이용하여 최적화를 수행한다. 수치 알고리듬으로서, 최적해까지의 수렴속도는 빠르지 않으나 수학적으로 적용이 용이한 steepest descent법을 적용하였다⁽¹⁷⁾.

(11)과 같이 정의된 수정된 목적함수에 대하여, 각 변수들의 미소변화량에 대한 $\overline { J }$의 미소변화량은 다음과 같이 표현된다.

단, 위 식에서 $\mathcal { V }$, $\eta$, $\mu$, $\lambda$와 같은 라그랑쥬 미정계수의 미소변화량의 계수는 과 동일하기 때문에 생략하였다. $\overline { J }$가 최소점이 되는 필요조건은 (12)의 미소변화량의 계수가 0이 되는 것이므로, 최적화를 위한 해석적 필요조건은 경계조건 및 구속조건 (9), (10)을 포함하여 다음과 같이 나타낼 수 있다.

그런데 수치적으로 문제를 풀 경우, 최적화하여야 할 변수는 $u$, $t _ { i }$, $t _ { f }$이므로 (14)의 좌변은 $\overline { J }$를 최대화하는, $u$, $t _ { i }$, $t _ { f }$에 대한 steepest ascent에 해당한다. 따라서 (12)를 만족시키기 위하여 ((14)를 $\delta u$, $\delta t _ { i }$, $\delta t _ { f }$를 다음과 같이 정의하여 대체한다.

단, $\alpha _ { u }$, $\alpha _ { t _ { i } }$, $\alpha _ { t _ { f } }$는 충분히 작은 상수이다. (15)와 같이 $u$, $t _ { i }$, $t _ { f }$의 값이 변동하였을 때 $\overline { J }$의 값이 가장 급격하게 감소하는 방향의 기울기를 $\delta u$, $\delta t _ { i }$, $\delta t _ { f }$에 대입함으로써 극소점(해)을 찾을 수 있다. 다음으로 라그랑쥬 미정계수를 구하기 위하여 $u$, $t _ { i }$, $t _ { f }$에 대한 구속조건 및 경계조건의 미소변화량에 대한 식을 표현하면 다음과 같다.

(15)를 (16)에 대입한 후 (13)과 연계하면 라그랑쥬 미정계수는 다음과 같이 나타낼 수 있다.

단, $\lambda _ { \nu _ { 1 } } = \Phi _ { \lambda } [ 0,0,0,1 ] ^ { \top }$, $\lambda _ { \nu _ { 2 } } = \Phi _ { \lambda } [ 0,0,-1,0 ] ^ { \top }$, $\Phi _ { \lambda } = \Phi _ { \lambda } \left( t ^ { - } , t _ { f } \right)$는 역시간 $t ^ { - }$에 대한 $\lambda$의 천이행렬(transition matrix),

$\mu _ { \nu _ { 1 } } = - \frac { \left( \frac { \partial C } { \partial u } \right) \left( \frac { \partial f } { \partial u } \right) ^ { \top } \lambda _ { \nu _ { 1 } } \nu } { \left( \frac { \partial C } { \partial u } \right) \left( \frac { \partial C } { \partial u } \right) ^ { \top } }$, $\mu _ { \nu _ { 2 } } = \frac { - \frac { d C } { \alpha _ { u } } - \left( \frac { \partial C } { \partial u } \right) \left( \frac { \partial f } { \partial u } \right) ^ { \top } \lambda _ { \nu _ { 2 } } } { \left( \frac { \partial C } { \partial u } \right) \left( \frac { \partial C } { \partial u } \right) ^ { \top } }$

여기까지 정의한 수식을 이용하여 최적화 해를 구하는 알고리듬은 다음과 같으며, 이를 그림. 2와 같은 순서도로 도식화하였다.

(a) 적당한 입력 $u$, $t _ { i }$, $t _ { f }$를 정한다.

(b) 운동방정식 로 시간방향으로 상태량을 구한 뒤, (13)을 이용하여 역시간 방향으로 $\lambda$의 천이행렬 $\Phi _ { \lambda }$를 구한다.

(c) 라그랑쥬 미정계수 $\nu$, $\mu$, $\lambda$, $\eta$를 차례대로 계산한다.

(d) (15)에 의하여 $\alpha _ { u }$, $\alpha _ { t _ { i } }$, $\alpha _ { t _ { f } }$를 변화시키며 $u \leftarrow u + \delta u$, $t _ { i } \leftarrow t _ { i } + \delta t _ { i }$, $t _ { f } \leftarrow t _ { f } + \delta t _ { f }$와 같이 갱신하여 (9)의 $J$를 최소로 하는 $u ^ { * }$, $t _ { i } ^ { * }$, $t _ { f } ^ { * }$를 찾는다(직선탐색).

(e) 충분히 작은 크기의 $\epsilon _ { J }$에 대하여 $| d \overline { J } | < \epsilon _ { J }$, $\psi \left( x _ { f } \right) = 0$의 조건을 동시에 만족하면 알고리듬을 종료하고, 만족하지 못하면 $u \leftarrow u ^ { * }$, $t _ { i } \leftarrow t _ { i } ^ { * }$, $t _ { f } \leftarrow t _ { f } ^ { * }$와 같이 갱신하고 (b)로 돌아가서 알고리듬을 반복한다.

4.1 조종날개 전개시각 최적화의 효과

그림. 3~그림. 6에 나타난 4개의 비행궤적의 초기조건은 $v _ { 0 }$= 750[m/s], $\gamma _ { 0 }$=45[deg], $r _ { 0 }$=0[m], $h _ { 0 }$=0[m]으로 동일하며, 다음 네 가지 경우의 비행조건을 가정하였다.

1. 발사 시부터 탄착 시까지 조종날개를 전개하지 않음

2. 조종날개를 발사 후 10[sec]에 전개하여 최적 유도비행 (조종날개 전개시점 최적화 미수행)

3. 조종날개를 발사 후 19.66[sec]에 전개, 최적 유도비행 (조종날개 전개시점 최적화)

4. 조종날개를 발사 후 35[sec]에 전개, 최적 유도비행 (조종날개 전개시점 최적화 미수행)

조종날개 전개시각이 19.66[sec]인 경우는 유도명령 $u _ { L }$와 함께 조종날개 전개시각 를 최적화 변수로 하여 사거리 최대화를 수행한 결과이며, 다른 두 경우는 유도명령 $u _ { L }$만을 최적화한 결과이다. 조종날개를 전개한 세 경우는 조종날개를 전개하지 않은 경우에 비교하여 모두 약 200%의 사거리를 비행하였으며, 그 중에서도 조종날개 전개시각 최적화를 수행한 경우가 가장 긴 사거리를 비행하였다. 전개시각 10[sec]의 경우, 그림. 4, 그림. 5(a)에서 확인할 수 있듯이 조종날개를 일찍 전개하여 확보한 속도를 이용하여 큰 양력을 받아 다른 두 경우에 비하여 높은 고도까지 비행한 반면 속도가 빠르게 줄었으며, 전개시각 35[sec]의 경우 늦은 조종날개 전개가 고도 확보 실패로 이어져 조종날개를 전개한 경우 중 가장 짧은 사거리를 비행하였다. 그림. 6에서 조종날개 전개 직후의 유도명령 $u _ { L }$은 세 경우에 각각 다른 양상을 보이지만, 비행궤적, 양력과 항력 프로파일은 정점 이후 세 경우가 거의 유사한 궤적을 보인다.

표 1을 보면, 발사 고각을 60[deg], 75[deg]로 설정한 경우의 최적화를 수행한 경우에도 비슷한 경향의 결과를 얻었다는 것을 알 수가 있다. 이 경우에도 조종날개 전개시점 최적화를 수행한 경우가 그렇지 않은 경우보다 긴 사거리를 비행하였으므로, 보다 긴 사거리를 비행하기 위해서는 무조건 고도가 높은 탄도 정점 부근에서 조종날개를 전개하기보다, 적정한 시각에 조종날개를 전개해야 한다는 결론을 내릴 수 있다.

또한, 표 1의 최적화 결과는 조종날개 전개시각 최적화를 한 경우가 가장 효율적인 비행을 했다는 것을 알려준다. 표 1이 제시하는 최적화 결과의 초기조건은 모든 경우에 대하여 포구속도 750 [m/s]로 같기 때문에 다음과 같이 정의되는 역학적 에너지 $E _ { m }$이 동일하다.

그러나 조종날개 전개시각을 최적화한 경우와 그렇지 않은 경우는 탄착 시($h=0$[m])의 비행속도가 다르며, 이것은 최종적인 $E _ { m }$이 서로 다르다는 것과 동일한 의미이다. 다른 식으로 표현하면, 발사 시의 역학적 에너지 $E _ { m0 }$와 탄착 시의 역학적 에너지 $E _ { mf }$를 정의했을 때, 비행 전체에서 소비된 역학적 에너지 $\Delta E _ { m } = E _ { m 0 } - E _ { m f }$가 세 경우에 대해 모두 다르다. 한편, 이동체의 이동 효율을 나타내는 무차원 지표인 Cost of Transport(이하 CoT로 표기)는 위에서 정의된 변수를 이용하여 다음과 같이 정의된다⁽¹⁸⁾.

따라서, 최종 CoT가 낮은 비행(이동)을 했을 때 그 물체는 효율적인 비행을 한 것을 의미한다. 이를 고려했을 때 표 1의 결과가 보여주듯이, 조종날개 전개시각 최적화를 수행했을 때가 CoT가 가장 낮아 효율적으로 비행했다고 할 수 있다. 이 경향은 다른 고각에서 발사한 경우에도 비슷하며, 발사 고각이 75[deg]인 경우가 대개 소비된 역학적 에너지 $\Delta E _ { m }$는 적은데도 불구하고 비행한 사거리가 짧기 때문에 CoT는 높은 것을 알 수 있다. 최적화 문제에서 역학적 에너지는 최적화 변수로 고려하지 않았음에도 불구하고 이러한 결과를 얻은 사실로부터 미루어 보았을 때, 가능한 한 역학적 에너지를 소비하지 않으면서 효율적으로 비행하는 것이 지능화 탄약의 사거리 최대화로 이어지는 전략이라는 것을 알 수 있다.

4.2 발사 고각에 따른 최적화 결과

그림. 7~그림. 14에 나타난 비행궤적의 초기조건은 3개의 상태변수 $v _ { 0 }$=750[m/s], $r _ { 0 }$=0[m], $h _ { 0 }$=0[m] 에 대해서는 동일하나, 각 발사 고각이 각각 $\gamma _ { 0 }$= 45[deg], 60[deg], 75[deg]로 다른 경우에 대하여 $u _ { L }$ 최적화 및 조종날개 전개시각 최적화를 수행한 결과이다. 직관적인 관점에서 최대사거리 획득을 위한 조종날개 전개시각은 탄도 비행으로부터 이어지는 고도 확보를 위하여 정점 부근이라고 생각할 수 있으나, 이 결과는 그에 반하여, 심지어는 발사고각 45[deg]로 사격한 저고도에서도 조종날개 전개시점이 상당히 정점 이전에 위치한 것을 알 수 있다. 그림. 9를 보면, 발사고각이 낮을 때는 조종날개를 일찍 전개함으로써 고도의 손실을 감수하며 일찍 양력을 얻어 사거리를 증가시키는 것이 총 사거리 상으로 더 큰 이득을 볼 수 있기 때문이라고 사료된다. 반면, 발사고각이 클 때는 조종날개 전개로 인한 항력 증가가 큰 사거리 대비 고도 손실을 야기하기 때문에, 조종날개를 전개하지 않은 형상으로 정점 부근에 이른 후 조종날개를 전개하는 것이 최적의 비행 시나리오라는 결론을 얻을 수 있다. 이것은 다른 포구속도에서 탄약을 발사했을　경우에도 동일하게 찾아볼 수 있는 경향으로, 표 2의 결과를 보면 모든 포구속도의 경우에 발사고각이 증가할수록 최적 조종날개 전개시각이 정점 이전에서 정점 부근으로 이동하는 것을 알 수 있다. 단, $t _ { APEX }$는 각 경우의 사거리 최대화 비행에서 정점고도를 획득했을 때의 시간이다.

또한, 4.1절의 경우와 마찬가지로, 그림. 9, 그림. 10을 보면, 발사고각과 조종날개 전개시각이 다르기 때문에 초기 앙력 및 항력 프로파일과 전체 궤적에서 $u _ { L }$의 양상은 다르지만 비행궤적이 안정된 이후의 양력 및 항력 프로파일, 양항비 프로파일이 세 경우가 모두 비슷한 수준인 것을 알 수 있다.

한편, 발사부터 탄착까지 사거리에 대한 역학적 에너지 순간 소비량 $\frac { d E _ { m } } { d r }$는 다음과 같이 계산될 수 있다.

식(20)은 $E _ { m }$ 순간소비량이 사거리에 대한 순간 이동효율(순간 CoT)이라는 것을 의미하고 있다. 그림. 11은 이와 같이 정의된 $E _ { m }$ 순간소비량과 비행 시 양항비를 비교하여 보여주는데, 사거리에 따른 양항비의 추이가 $E _ { m }$ 순간소비량과 밀접한 관계를 갖고 있음을 알 수 있다. 조종날개 전개 직후에는 양항비가 높아짐에 따라 $E _ { m }$ 순간소비량도 늘어나지만, 이내 양항비가 안정된 비행구간에 들어서면서는 서로 반비례의 관계가 되어, 입력 $u _ { L }$이 $u _ { Lmax }$에 포화되어 양항비가 감소하는 순간에는 $E _ { m }$ 순간소비량이 잠시 높아지다가, 이내 안정된 양항비를 되찾으면 다시 $E _ { m }$ 순간소비량도 안정하는 것을 알 수 있다. 따라서 비행 구간에 따라 사거리 최적화를 위하여 역학적 에너지를 소비하는 양상이 다르며, 활공비행 시에는 양항비를 높게 유지하도록 역학적 에너지를 소비하는 것이 사거리 최대화를 위한 전략이라고 판단된다.

흥미로운 점은, 그림. 12~그림. 13에서 알 수 있듯이 세 경우가 모두 발사 고각이 다름에도 불구하고, 역학적 에너지 $E _ { m }$가 탄착 지점에 이를수록 비슷한 경로로 수렴하여, 탄착 시에 도달한 최종 역학적 에너지 $E _ { mf }$의 값이 동일하다는 것이다. 이 경향은 모든 포구속도에 대하여 동일하게 관찰되는 경향인데, 표 3에 따르면 발사 고각이 다름에도 불구하고 발사 시 역학적 에너지, 탄착 시 역학적 에너지, 따라서 소비된 역학적 에너지가 모두 동일한 수준이다. 단, 그림. 14와 같이 CoT 관점에서 보았을 때는 발사 고각 60[deg]인 경우가 제일 효율적인 비행을 하였는데, 이는 고각 60[deg]에서 최대사거리 값이 가장 크기 때문이다. 또한, 표 3에 따르면 같은 고각에서 다른 포구속도로 발사했을 경우에 탄착 시의 역학적 에너지는 같은 수준인데, 포구속도가 빠를수록 비행에서 소비되는 역학적 에너지의 총량이 사거리에 비하여 커지기 때문에, 비행 효율이 떨어진다는 것을 그림 14에서 알 수 있다.