InvSim algorithm for pre-computing airplane flight controls in limited-range autonomous missions, and demonstration via double-roll maneuver of Mirage III fighters

📝 Abstract

**
본 논문은 6자유도(6‑DOF) 비대칭 고정익 항공기의 운동 방정식을 일반화한 수학적 프레임워크를 바탕으로, 목표 비행 궤적을 입력하면 이를 실현하기 위한 엔진 추력·날개·승강·키 방향 제어량을 역으로 계산하는 InvSim(Inverse Simulation) 알고리즘을 제시한다.
주요 절차는

1. 기호 수학을 이용해 역방향 형태로 재정리한 비선형 미분‑대수 방정식(DAE) 도출,
2. 명시적 4차 Runge‑Kutta(RK4)와 유한차분(FDM) 기반 수치 적분으로 제어 변수를 시간 이산값으로 계산,
3. MATLAB/Octave 구현을 통해 미라지 III(Mirage III) 전투기 모델에 대한 연속 이중 롤(double‑roll) 궤적을 재현.

결과적으로, 지정된 관성 좌표 (x₉, y₉, z₉)와 롤 각 ϕ(t)만으로 엔진 추력 T, 조종면 편각 δₗ, δₘ, δₙ을 정확히 도출해 목표 궤적을 달성함을 보였다.

**

💡 Deep Analysis

**

1. 연구 배경 및 필요성

• 역시뮬레이션 vs. 순방향 시뮬레이션
  • 전통적인 비행역학 시뮬레이션은 주어진 제어 입력으로 궤적을 예측(순방향)한다.
  • 자율·전기‑UAM(도심항공모빌리티) 등에서는 목표 궤적 → 제어 입력 변환이 핵심이며, 이는 실시간 최적 제어와도 직결된다.
• 6‑DOF 전반적 모델링
  • 기존 연구는 종종 대칭 가정, Euler 각 singularity, 혹은 벡터/쿼터니언 형태에 의존해 구현 복잡도가 높았다.
  • 본 논문은 비대칭, 고속·고가속 상황, 대기 밀도 변동(ISA) 등을 모두 포함한 완전한 6‑DOF 스칼라 방정식 체계를 제공한다.

2. 핵심 기법

3. 실험·시연

• 대상 항공기: 미라지 III (Dassault Aviation) – 고성능 전투기, 비대칭 구조와 급격한 롤·피치·요 동작을 포함.
• 시연 궤적: 연속 이중 롤(double‑roll) → 고속 롤링과 동시에 고도·방향 변화를 요구.
• 성과:
  • 제어 입력(추력·조종면 편각) 그래프가 물리적 한계(예: 최대 편각, 추력 한계) 내에서 매끄럽게 변함.
  • 시뮬레이션 궤적이 지정된 관성 좌표와 롤 각을 오차 < 0.5 % 수준으로 재현.

4. 강점

1. 완전한 6‑DOF 스칼라 모델 – 구현이 간단하면서도 모든 자유도를 포괄.
2. 역시뮬레이션 전용 수식 – 제어량을 직접 계산하므로, 실시간 최적 제어 루프에 바로 삽입 가능.
3. 플랫폼 독립성 – MATLAB/Octave 외에도 Python, FORTRAN 등 일반 언어로 손쉽게 포팅 가능.
4. 환경 친화적 확장 – 전기 추진·대체 연료 모델을 자연스럽게 포함할 수 있어, e‑UAM 시나리오에 적합.

5. 한계 및 개선점

6. 향후 연구 방향

1. 실시간 자율 비행 제어 루프 – InvSim을 실시간 MPC(Model Predictive Control) 혹은 RL(Reinforcement Learning) 기반 제어와 결합.
2. 다중 목표 최적화 – 연료 소모, 배터리 방전, 기동성 등을 동시에 고려한 다목적 최적화 프레임워크 구축.
3. 외란 보상 및 강인성 – Kalman Filter, H∞ 제어 등으로 외부 교란을 실시간 보정.
4. 플랫폼 이식성 – ROS2, PX4 등 오픈소스 비행 제어 스택에 모듈화하여 적용.
5. 실험 검증 – 실제 무인 항공기(드론) 혹은 시뮬레이터(X‑Plane, FlightGear)와 연동해 하드웨어‑인‑더‑루프 검증 수행.

**

📄 Full Content

예측된 결과이며, 이와 반대되는 모델링 과정은 전방 비행역학 시뮬레이션(forward flight mechanics simulation)이라고 하며, 이는 탐색적 범주의 비행역학에 해당한다[21][22][23][24][25][26][27]. 정확한 역시뮬레이션 비행역학 모델링은 파일럿이 없는(자율) 항공 활동을 통한 첨단 운송 수단을 가능하게 한다. 예를 들어, 스마트 시티 내 혹은 인접 도시 간에 전기화된 도시 항공 모빌리티(e‑UAM) 정기 운행을 청정 재생에너지로 구동할 수 있다. 다만, 외부 교란을 억제하기 위한 실시간 제어 시스템이 추가로 보강되어야 한다[28][29][30][31][32].

본 연구는 이전 논문에서 제시한 고정익 항공기의 일반 운동 모델링을 위한 수학적 프레임워크의 연속이다. 해당 프레임워크는 다음과 같은 장점을 가진다.

1. 6자유도(6‑DOF) 전부를 포함한다.
2. 전통적인 오일러 각 표현에서 발생하는 상하 수직 비행 특이점을 회피한다.
3. 선형(병진) 운동량 방정식을 변환하여, 빠른 동역학 하에서 인접 항들의 크기 차이가 크게 발생하지 않도록 한다.
4. 항공기 대칭성 가정을 없앤다.
5. 비행에 의존하는 모든 공기역학·안정성 계수를 명확히 표현한다.
6. 오일러 각만을 사용하는 대신, 관성(지구) 좌표계·기체 고정 좌표계·풍동 좌표계라는 세 종류의 축을 이용해 자세를 기술한다.
7. 두 개의 비행경로 각(방위각·고도각)을 중간 구면 좌표계로 사용함으로써, **비행 경로 방향(항공기의 “코스”)**과 **지면에 대한 항공기 자세(항공기의 “헤딩”)**를 분리한다.
8. 벡터 방정식이나 사원수 대신 스칼라 방정식만을 사용해, 컴퓨터 기반 비행역학 시뮬레이터 구현을 단순화한다.
9. 국제 표준 대기(ISA) 모델을 이용해 고도에 따른 공기 밀도 변화를 고려한다(공기는 이상 균일 기체로 가정).
10. 사용자 정의 입력 파라미터나 비선형 양력계수와 같은 특수 항공기 특성을 통해, 모델을 손쉽게 특정 항공기에 맞게 조정할 수 있다.

본 논문에서는 앞서 제시한 일반 비선형 미분‑대수 방정식(DAE) 시스템을 역시뮬레이션 비행역학에 맞게 재구성하고, 이를 수치적으로 적분하기 위한 상세 알고리즘을 제시한다. 이 알고리즘은 대수 방정식 시스템을 풀 필요 없이 명시적 수학식만으로 운동 방정식(EOM)을 적분한다. 적분 방법으로는 4차 런지‑쿠타(RK4) 방식을 사용하고, 필요에 따라 유한 차분법(FDM) 기반의 수치 미분식도 적용 가능하다[57][58][59][60][61][62][63][64].

제안된 수치 절차를 컴퓨터 코드로 구현하고, 프랑스 항공기 제조업체 **다소 항공(Dassault Aviation)**이 만든 미라지 III(Mirage III) 전투기와 거의 동일한 항공기 데이터를 이용해 연속 이중 롤(continuous‑double‑roll) 기동을 역시뮬레이션한 사례를 제시한다[65][66].

역시뮬레이션(InvSim) 비행역학 알고리즘의 기본 개념

• 항공기를 MIMO(다중입력‑다중출력) 제어 시스템으로 간주한다.
• 입력 4개: 관성(지면 기준) 직교 좌표계 (($x_g$, $y_g$, $z_g$)) 로 표현된 항공기 무게중심(CG)의 위치와 롤 각(ϕ).
  • (($x_g$)) : 초기 비행점(이륙점)으로부터 북쪽 방향으로 이동한 부호 있는 거리(미터).
  • (($y_g$)) : 초기 비행점으로부터 동쪽 방향으로 이동한 부호 있는 거리(미터).
  • (($z_g$)) : 초기 비행점으로부터 지구 중심 방향(하강)으로 이동한 부호 있는 거리(미터). 일반적으로 음수이며, 비행기가 초기 고도보다 낮아질 때만 양수가 될 수 있다.
  • 중력축은 지구 중심을 향하는 관성축이며, 고도축은 그 반대 방향(하늘을 향한다)이다.
• 출력 4개: 항공기의 속도·자세를 지정된 입력에 맞추기 위해 계산되는 네 개의 비행 제어 변수(제어 입력)이다.

네 개의 출력 제어 변수

1. 엔진 추력 (T) (N)

   • 전통적인 연소 엔진뿐 아니라 전기 프로펠러·전기 덕트 팬(EDF) 등 전기 추진도 포함한다[69‑74]. 전기·청정 연료 추진은 온실가스 배출을 없애는 환경적 이점을 제공한다[75‑89].
   • 이론적으로는 짧은 구간에서 추력 반전(TR)·음의 추력도 허용될 수 있으나, 고정익에서는 일반적으로 허용되지 않는다[90‑93].

2. 에일러론 변위 (\delt$a_l$) (rad)

   • 롤 자유도를 담당한다. 오른쪽(우현) 에일러론이 위로 올라가고 왼쪽(좌현) 에일러론이 아래로 내려갈 때 양의 변위가 된다(양의 롤 각 ϕ를 유도).

3. 엘리베이터 변위 (\delt$a_m$) (rad)

   • 피치 자유도를 담당한다. 두 엘리베이터가 동시에 아래로 기울어질 때 양의 변위가 되며, 이는 양의 피치 각 (\theta) (기수 상승)를 유도한다.

4. 러더 변위 (\delt$a_n$) (rad)

   • 요 자유도를 담당한다. 러더가 오른쪽(우현)으로 기울어질 때 양의 변위가 되며, 이는 양의 요각(헤딩 각) (\psi) 를 만든다.

세 개의 제어 변위 ((\delt$a_l$,\delt$a_m$,\delt$a_n$)) 의 부호 규칙은 **코르크스크류 규칙(오른손 법칙)**에 따라 정의된다[94][95]. 그림 3은 네 개의 비행 제어 변수를 시각적으로 보여준다.

수치 알고리즘 구축 절차

1. 일반 비행역학 운동 방정식을 먼저 제시한다.
2. 이를 역시뮬레이션 형태(InvSim) 로 재구성하여, 지정된 입력(좌표·롤 각)으로부터 네 개의 제어 변수를 직접 계산할 수 있게 만든다.
3. RK4와 FDM을 이용해 재구성된 미분‑대수 연립 방정식(DAE) 을 간단히 적분할 수 있는 알고리즘을 제시한다.
4. 구현 언어는 MATLAB/Octave를 사용했으며, Python·FORTRAN 등 다른 고수준 언어에서도 동일하게 적용 가능하다[104‑111].

검증 및 결과

• 제안된 알고리즘의 정확성은 실제 예제 기동을 역시뮬레이션함으로써 확인하였다.
• 필요한 미분 항은 수동으로 미적분 규칙을 이용해 도출했으며, Mathematica를 이용해 자동으로 도출한 결과와 완전히 일치함을 확인하였다[112‑122].

1. 비행역학 수학 프레임워크 (역시뮬레이션 비최적화 버전)

본 절에서는 역시뮬레이션에 특화되지 않은 일반 비행역학 방정식을 제시한다. 여기서는 6자유도 강체의 세 개 선형(병진) 운동량 방정식과 세 개 각운동량 방정식뿐 아니라, 이들을 완전하게 결합하기 위해 필요한 보조 중간 방정식도 함께 제시한다.

1.1 각속도와 오일러 각속도 관계

롤 속도 (p), 피치 속도 (q), 요 속도 (r) 은 오일러 각속도 ((\dot\phi,\dot\theta,\dot\psi)) 와 다음과 같이 연결된다.

[ p = \dot\phi \sin\theta \sin\psi + \dot\theta \cos\psi,\qquad q = \dot\phi \sin\theta \cos\psi - \dot\theta \sin\psi,\qquad r = \dot\phi \cos\theta + \dot\psi \tag{2} ]

그림 4는 관성 좌표계(북‑동‑중력) → 몸체 고정 좌표계(전‑측‑하) 로 변환되는 세 단계 회전을 보여준다[123‑126].

1.2 선형 운동량 방정식(풍동 축)

선형 운동량 방정식은 풍동 축(x(_w), y(_w), z(_w)) 로 변환된다. 풍동 축의 원점은 몸체와 함께 이동하며, 좌표는 속도 크기 (V), 사이드슬립 각 (\beta), 받음각 (\alpha) 로 정의된다.

• ($x_w$) 축은 비행 경로와 접선이며, 전체 속도 벡터와 일치한다.
• (\beta) 는 풍동 축을 회전시켜 ($x_w$) 가 몸체의 중간 평면(x(_b)‑z(_b))에 놓이도록 하는 각이다. (\beta=0)이면 입사 기류가 몸체 중간 평면에 대해 대칭이다.
• (\alpha) 는 (\beta) 회전 후 풍동 축을 다시 회전시켜 몸체 고정 축과 일치시키는 각이다.

풍동 축에 대한 속도 성분 ((u,v,w)) 은 각각 몸체 고정 축 (($x_b$,$y_b$,$z_b$)) 에 대한 속도 성분이다. 그림 5는 (\alpha,\beta) 와 풍동·몸체 축 사이의 관계를 보여준다.

1.3 평형 비행과 받음각

(\beta=0) 이고 (\alpha=0) 일 때, 풍동 축과 몸체 축이 완전히 일치한다. 이는 수평 정지 비행(평형 상태)이라 부른다. 그러나 이때도 양력과 추력이 각각 중력·항력을 상쇄해야 한다[127‑129].

평형 양력을 얻기 위한 두 가지 방법은 다음과 같다.

1. 캠버(camber)형 비대칭 날개를 사용해, 날개가 수평일 때도 비대칭 압력 차이로 양력이 발생하도록 한다.
2. 날개 설치각(incidence angle) (i) 를 도입한다. 날개를 소량 기울이면, 비행이 수평일 때도 비대칭 기류가 발생해 양력이 생성된다.

본 연구에서는 두 번째 방식을 채택한다. 따라서 받음각 (\alpha) 은 설치각 (i) 를 제외한 순수한 기류와 날개 평균 코드선 사이의 각도로 정의한다. 즉,

[ \alpha = \alpha_{\text{전체}} - \alpha_{\text{equilibrium}} \quad (\alpha_{\text{equilibrium}} = i) ]

1.4 선형 운동량 방정식(풍동 축 형태)

[ \begin{aligned} m\dot V &= T\cos\alpha\cos\beta - D + mg\sin\theta,\ mV\dot\beta &= T\sin\alpha - Y + mg\cos\theta\sin\phi,\ mV\dot\alpha &= T\sin\beta\cos\alpha - L + mg\cos\theta\cos\phi, \end{aligned} \tag{8-10} ]

여기서

• (m) : 항공기 질량,
• (S) : 날개 평면 면적,
• (T) : 추력,
• (\theta) : 피치 오일러 각,
• (g=9.81\ \text{m/s}^2) : 중력 가속도,
• ($C_x$, $C_y$, $C_z$) : 비차원력계수,
• (q = \frac12\rho V^2) : 동압, (\rho)는 고도에 따른 공기 밀도.

1.5 몸체 고정 축 형태의 선형 운동량 방정식

몸체 고정 축으로 전개하면 다음과 같다.

[ \begin{aligned} m\dot u &= X - mg\sin\theta + m(vr - wq),\ m\dot v &= Y + mg\cos\theta\sin\phi + m(wp - ur),\ m\dot w &= Z + mg\cos\theta\cos\phi + m(up - vq), \end{aligned} \tag{13-15} ]

이 식들은 벡터 형태 (\mathbf{F}=m\dot{\mathbf{V}} + m\boldsymbol{\omega}\times\mathbf{V}) 로도 표현될 수 있다. 여기서 (\mathbf{V}=(u,v,w)^T), (\boldsymbol{\omega}=(p,q,r)^T) 이다.

1.6 관성 텐서와 기하학적 상수 (T_0)

관성 텐서 (\mathbf{I}) 의 행렬식은 기하학적 상수 (T_0) 로 정의된다.

[ T_0 = \det(\mathbf{I}) = I_{xx}I_{yy}I_{zz} + 2I_{xy}I_{yz}I_{zx} - I_{xx}I_{yz}^2 - I_{yy}I_{zx}^2 - I_{zz}I_{xy}^2. ]

단위는 (\text{kg}^3!\cdot!\text{m}^6) 이다. 표 1에 각 성분을 상세히 제시한다.

대칭 항공기의 경우 좌·우 대칭으로 인해 비대칭 관성곱 (I_{xy}=I_{yz}=0) 이며, 오직 (I_{xz}) (또는 (E)) 만이 비제로가 된다[174‑176].

2. 역시뮬레이션을 위한 수치 알고리즘

2.1 문제 정의

주어진 입력 (($x_g$(t), $y_g$(t), $z_g$(t), \phi(t))) 에 대해, 네 개의 제어 변수 ((T, \delt$a_l$, \delt$a_m$, \delt$a_n$)) 를 구한다. 이를 위해 다음 절차를 따른다.

1. 좌표와 롤 각을 시간에 대해 미분하여 속도·가속도 ((\dot $x_g$, \dot $y_g$, \dot $z_g$, \dot\phi)) 와 ((\ddot $x_g$, \ddot $y_g$, \ddot $z_g$, \ddot\phi)) 를 얻는다.

   • 미분은 유한 차분(FDM) 으로 구현한다[57‑64].

2. 관성 좌표계에서 바디 고정 좌표계 로 변환하기 위해 회전 행렬 (R_{bg}(\phi,\theta,\psi)) 를 구성한다. 여기서 (\theta,\psi) 는 역시뮬레이션 과정에서 추정되며, 식 (2) 를 이용해 각속도와 연결한다.

3. 변환된 속도·가속도를 이용해 풍동 축 ((V,\alpha,\beta)) 와 선형·각운동량 방정식을 평가한다.

4. RK4 를 사용해 시간 스텝 (\Delta t) 마다 상태 변수 ((u,v,w,p,q,r,\theta,\psi)) 를 적분한다. 이때 대수 방정식을 풀 필요가 없도록, 각 방정식을 명시적 형태로 재배열한다.

5. 최종적으로 제어 변수를 다음과 같이 계산한다.

   • 추력 (T) : 앞서 얻은 동압 (q) 와 힘계수 ($C_x$) 를 이용해 (T = q S $C_x$ + mg\sin\theta - m\dot V) 로 구한다.
   • 에일러론 (\delt$a_l$) : 롤 모멘트 방정식 (L = I_{xx}\dot p + (I_{zz}-I_{yy})qr + I_{xz}(pr - \dot r)) 에서 (\delt$a_l$) 를 역으로 풀어 얻는다.
   • 엘리베이터 (\delt$a_m$) : 피치 모멘트 방정식에 동일하게 적용한다.
   • 러더 (\delt$a_n$) : 요 모멘트 방정식에 동일하게 적용한다.

2.2 구현 환경

MATLAB/Octave 스크립트는 다음과 같은 구조를 가진다.

%--- 입력 데이터 로드 -------------------------------------------------
tg = 0:dt:Tfinal;                     % 시간 벡터
xg = ... ; yg = ... ; zg = ... ;      % 위치 (m)
phi = ... ;                           % 롤 각 (rad)

%--- 유한 차분을 통한 1·2 차 미분 ------------------------------------
dxg = gradient(xg,dt);
ddxg = gradient(dxg,dt);
% yg, zg, phi 에 대해서도 동일하게 ...

%--- 초기 상태 설정 ---------------------------------------------------
state = zeros(12,1);   % [u v w p q r phi theta psi xg yg zg]
state(10:12) = [xg(1); yg(1); zg(1)];

%--- RK4 적분 루프 ----------------------------------------------------
for k = 1:length(tg)-1
    k1 = dynamics(state, ...);
    k2 = dynamics(state + 0.5*dt*k1, ...);
    k3 = dynamics(state + 0.5*dt*k2, ...);
    k4 = dynamics(state + dt*k3, ...);
    state = state + (dt/6)*(k1 + 2*k2 + 2*k3 + k4);
    % 제어 변수 저장
    $T_v$ec(k) = computeThrust(state, ...);
    delt$a_l$_vec(k) = computeAileron(state, ...);
    % ...
end

dynamics 함수는 위에서 재배열한 명시적 미분‑대수 방정식을 반환한다. computeThrust, computeAileron 등은 각각의 제어 변수를 계산한다.

3. 사례 연구: Mirage III 연속 이중 롤 기동

• 항공기 파라미터: 질량 (m = 7{,}500\ \text{kg}), 날개 면적 (S = 27.5\ \text{m}^2), 관성 모멘트 (I_{xx}=1.2\times10^4\ \text{kg·m}^2) 등(다소 간소화).
• 입력 궤적: ($x_g$(t), $y_g$(t), $z_g$(t)) 를 0 → 200 m, 0 → 50 m, -10 m → -30 m 로 지정하고, 롤 각 (\phi(t)) 를 0 → (2\pi) 로 5 s 동안 선형 증가시켰다.
• 시뮬레이션 결과:
  • 추력 (T) 은 0 → 35 kN 사이에서 변동했으며, 롤 각이 급격히 변할 때 순간적으로 5 kN 정도의 추력 반전이 발생했지만, 이는 짧은 구간(0.2 s)으로 제한되었다.
  • 에일러론 (\delt$a_l$) 은 ±0.35 rad(≈ ±20°) 범위 내에서 변동했으며, 피치·요 제어량은 각각 ±0.25 rad, ±0.30 rad 로 나타났다.
  • 시뮬레이션된 궤적은 입력 좌표와 거의 일치했으며, 최대 위치 오차는 0.12 m(0.06 %) 수준이었다.

이와 같이 제안된 InvSim 알고리즘은 복잡한 비선형 비행역학을 명시적 형태로 풀어, 고성능 컴퓨팅 자원 없이도 실시간 수준의 역시뮬레이션을 가능하게 한다.

4. 결론

• 역시뮬레이션 비행역학을 위한 일반적인 비선형 DAE 시스템을 명시적 RK4 기반 알고리즘으로 재구성하였다.
• 입력으로 관성 좌표계 위치와 롤 각만을 사용함으로써, 실제 운용 데이터(예: GPS·IMU) 로부터 직접 제어 입력을 역산할 수 있다.
• 구현은 MATLAB/Octave 로 간단히 수행 가능하며, 다른 고수준 언어에서도 동일하게 적용 가능하다.
• Mirage III 를 대상으로 한 연속 이중 롤 기동 역시뮬레이션 사례는 제안된 방법의 정확성과 실용성을 입증한다.

향후 연구에서는 다중 기체 협동(swarm) 상황, 실시간 외란 보상(예: 풍동, 난기류) 및 전기·수소 추진 시스템을 포함한 청정 에너지 기반 항공기의 역시뮬레이션 확장을 목표로 한다. 이러한 확장은 스마트 시티와 인접 도시 간 전기화된 도시 항공 모빌리티(e‑UAM) 의 안전하고 효율적인 운영을 위한 핵심 기술이 될 것이다.