• 대한전기학회
Mobile QR Code QR CODE : The Transactions of the Korean Institute of Electrical Engineers
  • COPE
  • kcse
  • 한국과학기술단체총연합회
  • 한국학술지인용색인
  • Scopus
  • crossref
  • orcid

  1. (Dept. of Electrical and Electronic Engineering, Incheon National University, Korea.)



Unmanned aerial vehicle, Multirotor aerial vehicle, Fault diagnosis and isolation algorithm, Extended Kalman filter, Interacting multiple model

1. 서 론

UAV(Unmanned aerial vehicle)는 지상에서 운용하는 UGV(Unmanned ground vehicle)보다 더 높은 자유도를 가져 다양한 영역에 적용할 수 있지만 더 많은 외란에 영향을 받고, 탑재할 수 있는 장비 무게에 제한이 있다. 하지만 모터, 센서 그리고 배터리 등의 부품 소형화 개발이 진행됨에 따라 UAV 활용 기술 연구개발이 촉진되고 있다. 이 중 multi rotor를 가진 MAV(Multirotor Aerial Vehicle)는 다른 기체보다 높은 안정도 그리고 수직 이착륙이 가능한 이점이 있어 많은 연구가 진행되고 있다. MAV의 대표적인 활용 예로 quadrotor가 있으며 농업(1), 운송(2), 탐색(3) 그리고 구조(4)에 활용하기 위한 연구개발이 수행되고 있다.

하지만 quadrotor는 비행 중 모터가 파손되면 기체 힘의 불균형으로 추락할 수 있으며, 이로 인해 인명 재산 피해가 발생할 수 있다. 이러한 문제를 방지하고 큰 하중 무게를 견디기 위해 더 많은 모터를 가진 MAV 연구를 지향하고 있다. 하지만 많은 모터를 가졌다고 해도 고장 상황에 이를 활용한 알고리즘이 개발되지 않는다면 기체의 고장에 대해 큰 효과를 기대하기 어렵다. 이와 같은 이유로 기체 고장 시 정상 부품을 활용하거나 비상 행동을 취할 수 있도록 고장을 인식하고 이를 분류할 수 있는 알고리즘 개발이 요구된다.

고장을 분류하는 방법에는 신호 기반, 데이터 기반 그리고 모델기반 고장 분류방법으로 나누어질 수 있다. 신호 기반 고장 분류방법에는 WT(Wavelet Transform)와 FFT(Fast Fourier Transform)로 주파수 영역에서 신호를 분석하는 방법이 있으며(5), 데이터 기반 방법은 고장 분류 모델을 학습시켜 고장과 정상을 분류한다. 데이터 기반 고장 분류방법으로 SVM(Support Vector Machine) 또는 LSTM (Long Short-Term memory)을 적용한 방법이 널리 사용되고, (6)은 관성 센서 데이터를 SVM에 적용하여 고장을 분류했으며, (7)은 PCA(Principal Component Analysis)와 LSTM을 사용하여 향상된 연산으로 센서 고장을 분류하였다. 모델기반 고장 분류방법은 target 모델기반으로 필터를 모델링하여 고장을 분류하는 방법이 선호된다. (8)은 드론의 추력 고장을 시스템에 bias로 모델링 하여 이중 칼만필터로 고장을 예측하였으며 (9-10)은 드론의 센서 고장과 차량 EPS의 고장을 사전에 설계된 다중 필터 모델 고장을 분류하였다.

시스템의 신호를 분석하여 고장을 분류하는 FFT는 동작특성이 다양한 시스템에 대해서는 유효할 수 있으나, octorotor와 같이 동작 특성이 단조로운 시스템에는 제한된다. 데이터 기반 고장 분류는 딥러닝 알고리즘으로 높은 신뢰도의 고장 분류를 수행하지만, 고장 분류 모델을 설계하기 위한 정밀한 데이터 수집과정이 요구된다. 이와 같은 특성들을 고려하여 MAV기체의 고장을 분류하기 위해 모델기반 방법을 채택하였으며, 단일 모델로는 다양한 상태(정상, 고장)의 시스템 동작을 나타낼 수 없으므로 다중 모델을 결합하는 상호 다중 모델 IMM(Interacting Multiple Model)을 활용하였다.

IMM은 단일 필터 모델로 예측할 수 없는 target 모델의 상태(state)를 예측하기 위해 복수의 필터 모델을 Markov chain rule에 따라 결합하는 방법으로(11) 빠르고 정확한 예측을 보여준다. IMM은 역동적으로 움직이는 target 모델 추적에 활용할 수 있으며(12), 앞서 언급했듯이 (9-10)과 같이 target 모델에 고장으로 변화된 시스템 특성을 감지하여 고장을 분류할 수도 있다. 이를 기반으로 MAV 기체의 프로펠러 고장이 나타나 있는 시뮬레이션에 IMM을 적용하고 프로펠러 고장을 분류하였다.

고장 분류는 MATLAB/Simulink 시뮬레이션 상에서 진행되었으며, 고장을 기체 동역학(13)의 변화로 간주하여 나타내었다. 또한, 드론의 주요 부품인 모터와 배터리를 모델링하여 더욱 현실에 가까운 기체의 동작을 관측하였는데, 모터는 brushless 모터를 사용하는 기체를 고려하여 모델링을 했고 배터리는 SoC(State of Charge)연구에 흔히 사용되는 ECM(Equivalent Circuit Model)으로 모델링을 수행하였다. 본 논문에서는 동역학 기반으로 MAV 기체의 프로펠러 고장을 나타내고, 복수의 필터 모델을 Markov chain rule에 따라 결합하는 IMM으로 고장을 감지하고 분류하였다.

2. Octorotor system dynamic

2.1 Multirotor Aerial Vehicle frame

그림. 1. 드론 구조 및 좌표계

Fig. 1. MAV structure and frame

../../Resources/kiee/KIEE.2022.71.5.744/fig1.png

Octorotor의 동체 동역학에 확인하기 위해서는 우선 그림 1의 관성 좌표계(inertial frame)와 기체 좌표계(body frame)에 대해 이해해야 한다. 관성 좌표계는 정지한 좌표계에서 기체의 움직임을 관측할 수 있고, 기체 좌표계는 기체 중심에서 기체의 움직임을 관측한다. 기체의 상태는 시뮬레이션 운용자의 기호에 따라 관성 좌표계 또는 기체 좌표계로 자유롭게 설정할 수 있으며 본 논문에서 GPS와 IMU 센서 활용을 근거로 하여 기체의 관측 좌표계(frame)를 선정하였다.

GPS와 IMU는 MAV 기체에서 활용하는 가장 대표적인 센서로 각각의 센서의 관측값은 관성 좌표계와 기체 좌표계로 분류할 수 있다. GPS와 IMU는 기체의 위치와 자세를 관성 좌표계로 측정하며 가속도, 각 가속도 그리고 각속도 측정값은 IMU에서 기체 좌표계 기준으로 관측한다. 이를 고려하여 기체의 움직임 X, Y, Z와 자세 Φ, Θ, Ψ[rad]를 관성 좌표계로 기체의 가속도 $\dot{V}_{x}^{B}$, $\dot{V}_{y}^{B}$, $\dot{V}_{z}^{B}$와 각속도p, q, r[rad/s]를 기체 좌표계로 나타냈다.

시뮬레이션 운용자가 기체상태 관측 좌표계를 선정하였으면, 기체-관성 좌표계의 변환 행렬을 이용하여 기체의 상태를 선정된 좌표계로 변환해 준다. 식 (1), (2), (3)은 각 축을 기준으로 관성 좌표계를 기체 좌표계로 나타내는 변환행렬이고, 기체의 자세 Φ, Θ, Ψ에 따라 삼각 함수 cosine(c)과 sine(s)으로 나타낸다.

(1)
$\mathrm{R}_{\mathrm{X}}(\phi)=\left[\begin{array}{ccc}1 & 0 & 0 \\ 0 & \mathrm{c} \phi & -\mathrm{s} \phi \\ 0 & \mathrm{s} \phi & \mathrm{c} \phi\end{array}\right]$

(2)
$\mathrm{R}_{\mathrm{y}}(\theta)=\left[\begin{array}{ccc}c \theta & 0 & -\mathrm{s} \theta \\ 0 & 0 & 0 \\ \mathrm{~s} \theta & 0 & c \theta\end{array}\right]$

(3)
$\mathbf{R}_{z}(\psi)=\left[\begin{array}{ccc}\mathrm{c} \psi & -\mathrm{s} \psi & 0 \\ \mathrm{~s} \psi & \mathrm{c} \psi & 0 \\ 0 & 0 & 1\end{array}\right]$

(4)
$\mathrm{R}_{\mathrm{BE}}=\mathrm{R}_{\mathrm{z}}(\psi) \mathbf{R}_{\mathrm{y}}(\theta) \mathrm{R}_{\mathrm{x}}(\phi)$ $\mathrm{R}_{\mathrm{BE}}=\left[\begin{array}{ccc}c \theta c \psi & \mathrm{s} \phi_{\mathrm{s}} \theta{c} \psi-\mathrm{c} \phi_{\mathrm{s}} \psi & \mathrm{c} \phi{\mathrm{s}} \theta \mathrm{c} \psi+\mathrm{s} \phi{\mathrm{S}} \psi \\ c \theta \mathrm{s} \psi & \mathrm{s} \phi{\mathrm{s}} \theta_{\mathrm{s}} \phi+c \phi{\mathrm{c}} \psi & \mathrm{c} \phi{\mathrm{s}} \theta{\mathrm{s}} \psi-\mathrm{s} \phi{\mathrm{c}} \psi \\ -\mathrm{s} \theta & \mathrm{s} \phi \cdot \theta & \mathrm{c} \phi \boldsymbol{c} \theta\end{array}\right]$

(5)
$\mathrm{T}_{\mathrm{BE}}=\left[\begin{array}{ccc}1 & \mathrm{s} \phi \mathrm{t} \theta & \mathrm{c} \phi \mathrm{t} \theta \\ 0 & \mathrm{c} \phi & -\mathrm{s} \phi \\ 0 & \frac{\mathrm{s} \phi}{\mathrm{c} \theta} & \frac{\mathrm{c} \phi}{\mathrm{c} \theta}\end{array}\right]$

각 축의 변환 행렬 식 (1), (2), (3)들을 결합하여 관성 좌표에서 기체의 좌표로 변환하는 식 (4)를 유도할 수 있으며, 아래 첨자$_{BE}$ 기체 좌표에서 관성 좌표로 변환을 의미한다. 회전 운동에서 관성 좌표계와 기체 좌표계의 관계는 식 (5)와 같이 나타낼 수 있다.

2.2 Octorotor dynamic

기체의 동역학을 분석하기 위해 기체 프레임을 강체(rigid body)로 가정한 후, 회전 운동 하는 강체의 동작을 해석할 수 있는 Newton-Euler 방정식을 적용한다(13).

(6)

$\mathrm{V}^{\mathrm{B}}=\left[V_{\mathrm{x}}, V_{\mathrm{y}}, V_{z}\right]^{\mathrm{T}}, \omega^{\mathrm{B}}=\lceil p, q, r\rceil^{T}$

$\mathrm{~F}_{\mathrm{ext}}=m\left(\dot{\mathrm{V}}^{\mathrm{B}}+\omega^{B} \times \mathrm{V}^{\mathrm{B}}\right)$

기체의 병진 운동을 나타내는 식 (6)의 첨자 는 기체 좌표계에서 관측한 값을 의미하며, 식 (6)에 대한 유도는 Appendix에 나타냈다.

그림. 2. 드론 상위, 하위 제어기 구조

Fig. 2. MAV upper/lower controller structure

../../Resources/kiee/KIEE.2022.71.5.744/fig2.png

$V^{B}$[m/s], $\omega^{B}$는 기체 좌표계에서 관측한 각 축의 속도와 회전 속도를 의미하며 p, q, r은 roll, pitch, yaw의 기체의 자세 변화율 m은 기체의 질량을 나타낸다. 식 (6)을 통해 외력 $F_{ext}$가 기체의 병진 운동 동역학에 주는 영향을 알 수 있고 기체의 회전 운동 또한 Newton-Euler 방정식으로 식 (7)과 같이 풀 수 있다(13).

(7)
$\begin{aligned} \tau_{\text {ext }} &=\mathrm{J} \dot{\omega}^{\mathrm{B}}+\omega^{\mathrm{B}} \times \mathrm{J} \omega^{\mathrm{B}} \\ \mathrm{J} &=\left[\begin{array}{ccc}J_{\mathrm{xx}} & 0 & 0 \\ 0 & J_{\mathrm{yy}} & 0 \\ 0 & 0 & J_{z z}\end{array}\right] \end{aligned}$

식 (6)과 유사하게 기체의 회전 운동 동역학 식 (7) 또한 기체에 작용하는 외력 $\tau_{ext}$[N·m]로 동작하며, J[kg·$m^{2}$]는 x, y, z 축으로 서로 상관관계가 없는 기체의 관성 모멘트이다. 식 (7)에 대한 유도 또한 Appendix에 나타냈으며, 식 (6)(7)로 정리한 기체 동역학은 식 (8)과 같이 나타낼 수 있다.

(8)

${\left[\begin{array}{cc}\mathrm{M} & 0_{3 \times 3} \\ 0_{3 \times 3} & \mathrm{~J}\end{array}\right]\left[\begin{array}{c}\dot{V}^{\mathrm{B}} \\ \dot{\omega}^{\mathrm{B}}\end{array}\right]+\left[\begin{array}{c}\omega^{\mathrm{B}} \times \mathrm{M} V^{\mathrm{B}} \\ \omega^{\mathrm{B}} \times \mathrm{M} \omega^{\mathrm{B}}\end{array}\right]=\left[\begin{array}{c}\mathrm{F}_{\text {ext }} \\ \tau_{\text {ext }}\end{array}\right] }$

$\mathrm{M}=\left[\begin{array}{ccc}m & 0 & 0 \\ 0 & m & 0 \\ 0 & 0 & m\end{array}\right]$

식 (8)의 은 x, y, z축의 힘의 작용을 나타내기 위해 무게 m을 각 축에 대해 augment 한 것이다. 기체에 작용하는 외력은 모터의 추력과 회전력, 중력 그리고 gyroscopic effect로 기체 동작에 영향을 준다. Gyroscopic effect는 회전하는 강체에 작용한 힘을 해석한 것이며, 이를 포함하여 기체의 상태 방정식을 나타내면 식 (9)와 같다.

(9)

$\dot{\mathrm{V}}^{\mathrm{B}}=\left[\begin{array}{c}0 \\ 0 \\ \frac{T}{m}\end{array}\right]+g\left[\begin{array}{c}\mathrm{s} \theta \\ -\mathrm{s} \phi \boldsymbol{c} \theta \\ -\mathrm{c} \psi \mathrm{c} \theta\end{array}\right]-\mathrm{S}\left(\omega^{\mathrm{B}}\right) \mathrm{V}^{\mathrm{B}}$

$\dot{\omega}^{\mathrm{B}}=\mathrm{J}^{-1}\left(\tau-\mathrm{S}\left(\omega^{\mathrm{B}}\right) \mathrm{J} \omega^{\mathrm{B}}-\mathrm{S}\left(\omega^{\mathrm{B}}\right)\left[\begin{array}{c}0 \\ 0 \\ J_{r} \Omega\end{array}\right]\right)$

$J_{\mathrm{r}} \Omega=\sum_{\mathrm{D}=1}^{8}-1^{\mathrm{D}} J_{\mathrm{r}, \mathrm{D}} \Omega_{\mathrm{D}}, \quad \tau=\left[\tau_{x} \tau_{y} \tau_{z}\right]^{\mathrm{T}}$

$\mathbf{u}=\left[\begin{array}{llll}T & \tau_{\mathrm{x}} & \tau_{\mathrm{y}} & \tau_{\mathrm{z}}\end{array}\right]^{\mathrm{T}}$

S($\omega^{B}$) skew matrix는 기체 동역학을 해석하기 위한 유용한 도구이며, T[N]는 기체의 추력 그리고 g는 중력[m/s$^{2}$]이다. $J_{r,D}$[kg·m$^{2}$]과 Ω$_{D}$[rad/s]는 p 번째 모터의 관성 모멘트와 모터의 회전 속도로 식 (9)에서 gyroscopic effect를 계산하는 데 활용한다. 식 (9)의 기체의 가속도와 회전 각속도 외의 이동 거리와 자세는 식 (9)의 동역학에 식 (4), (5)의 회전 변환을 적용한 뒤 적분을 수행하면 유도할 수 있다. 이를 통해 추력과 회전력이 전달되었을 때의 기체의 움직임을 시뮬레이션으로 구현할 수 있게 되었고, 이제 모터의 회전 속도와 동역학의 관계를 정리하기 위해 추력 T와 각 축의 회전력 $\tau$을 기체 동역학의 입력 u로 나타내었다. 기체 동역학 입력 u는 모터의 회전 속도, 위치 그리고 각 비례 계수에 의해 결정되며 모터와 식 (10)의 관계를 가진다.

(10)

$\Omega_{\text {square }}=\left[\begin{array}{llllll}\Omega_{1}^{2} & \Omega_{2}^{2} & \Omega_{3}^{2} & \Omega_{4}^{2} & \Omega_{5}^{2} & \Omega_{6}^{2} & \Omega_{7}^{2} & \Omega_{\mathrm{S}}^{2}\end{array}\right]$

$A=\left[\begin{array}{cccccc}b & b & b & b & b & b \\ b l & \frac{\sqrt{2}}{2} b l & 0 & -\frac{\sqrt{2}}{2} b l-b l-\frac{\sqrt{2}}{2} b l & 0 & \frac{\sqrt{2}}{2} b l \\ 0 & -\frac{\sqrt{2}}{2} b l-b l-\frac{\sqrt{2}}{2} b l & 0 & \frac{\sqrt{2}}{2} b l & b l & \frac{\sqrt{2}}{2} b l \\ -d & d & -d & d & -d & d\end{array}\right]$

$\mathrm{u}=\mathrm{A} \Omega_{\text {square }}$

b[N·s$^{2}$/rad$^{2}$]는 추력 상수, l[m]은 기체의 무게중심에서 모터까지 길이, d[N·m·s$^{2}$]는 drag constant 그리고 는 모터의 속도를 기체의 추력과 회전력으로 변환하는 변환 행렬이다. 이 중 b는 각 모터의 추력과 연동되는 매개체로 모터의 고장으로 정상적인 추력을 발생할 수 없을 때 조정할 수 있는 매개 변수(parameter)이다.

2.3 MAV controller

그림 2의 모터의 상위 제어기는 모터의 지령(reference) 회전 속도, 하위제어기는 모터의 전압을 조절하여 모터를 적정속도로 회전시킨다. 하위제어기에서 전달하는 모터의 회전 속도는 식 (10)의 과정으로 MAV 기체를 동작시키며, 기체의 동작은 상위-하위제어기로 피드백되어 기체가 목표 지점에 도달할 수 있도록 다시 모터의 지령속도와 전압을 제어한다.

위치 제어기(X-Y position controller)는 식 (11)과 같이 기체를 지령 위치 $x_{ref}$, $y_{ref}$로 이동시키기 위해 기체에 필요한 자세 값 $Φ_{ref}$, $Θ_{ref}$를 계산하며, 본 논문에서는 안정적인 비행을 위해서 지령 각도를 ±3°로 제한을 두었다.

(11)
$\begin{aligned} {\left[\begin{array}{l}\theta_{\text {ref }} \\ \phi_{\mathrm{ref}}\end{array}\right] } &=\left[\begin{array}{cc}c \psi & -\mathrm{s} \psi \\ \mathrm{s} \psi & \mathrm{c} \psi\end{array}\right]\left(\left[\begin{array}{cc}K_{\mathrm{Dx}} & 0 \\ 0 & K_{\mathrm{Dy}}\end{array}\right]\left[\begin{array}{c}V_{\mathrm{x}} \\ V_{\mathrm{y}}\end{array}\right]\right.\\ &\left.+\left[\begin{array}{cc}K_{\mathrm{Px}} & 0 \\ 0 & K_{\mathrm{Py}}\end{array}\right]\left[\begin{array}{ll}x_{\mathrm{ref}} & -x \\ y_{\mathrm{ref}} & -y\end{array}\right]\right) \end{aligned}$

식 (11)과 운용자가 입력해준 자세와 고도 지령 값은 각각 자세와 고도 제어기(Attitude & Altitude control)에 전달하여 기체의 추력과 회전력 지령 값을 계산하기 위해 사용되며 그 결과는 식 (12)와 같다.

(12)

$\mathrm{u}_{\mathrm{ref}}=\left[\begin{array}{l}T_{\mathrm{ref}} \\ \tau_{\mathrm{ref}}\end{array}\right]=\left[\begin{array}{cc}-\mathrm{K}_{\mathrm{D} z} & 0_{1 \times 3} \\ 0_{1 \times 1} & -\mathrm{K}_{\mathrm{D} \varphi}\end{array}\right]\left[\begin{array}{c}\dot{z} \\ \omega^{\mathrm{B}}\end{array}\right]$

$+\left[\begin{array}{ll}\mathrm{K}_{\mathrm{P} z} & 0_{1 \times 3} \\ 0_{1 \times 1} & \mathrm{~K}_{\mathrm{P} \varphi}\end{array}\right]\left[\begin{array}{c}z-z_{\mathrm{ref}} \\ \varphi_{\mathrm{ref}}-\varphi\end{array}\right]+\left[\begin{array}{c}g \\ 0_{3 \times 1}\end{array}\right]$

$\varphi=\left[\begin{array}{ll}\phi \theta & \psi\end{array}\right]$

$\mathrm{K}_{\mathrm{D} \varphi}=\left[J_{\mathrm{xx}} K_{\mathrm{D} \phi} J_{\mathrm{yy}} K_{\mathrm{D} \theta} \quad J_{z z} K_{\mathrm{D} \psi}\right]$

$\mathrm{K}_{\mathrm{P} \varphi}=\left[J_{\mathrm{xx}} K_{\mathrm{P} \phi} J_{\mathrm{yy}} K_{\mathrm{P} \theta} J_{z z} K_{\mathrm{P} \psi}\right]$

추력과 회전력 지령을 유도하기 위해 PD제어기를 적용하였으며, $\mathrm{K}_{\mathrm{D}, \varphi}$, $\mathrm{K}_{\mathrm{P}, \varphi}$는 회전력 지령 값을 계산하기 위한 제어 이득 $\mathrm{K}_{\mathrm{D}, \mathrm{Z}}$, $\mathrm{K}_{\mathrm{P}, \mathrm{Z}}$는 추력 지령을 계산하기 위한 제어 이득이다. 식 (12) 과정 이후 하위제어기가 추력과 회전력을 전달받아 모터를 제어하기 위해서는 $u_{ref}$를 모터의 지령 회전 속도로 변환해야 한다. Quadrotor와 같이 모터의 개수가 4개인 경우에는 식 (10)의 변환 행렬이 4x4 크기로 이에 대한 역변환(inverse)을 바로 이용할 수 있으나 8개의 모터를 가진 octorotor의 경우 정방행렬이 되지 않기 때문에 이를 바로 이용할 수 없다. 이를 해결하기 위해 식 (13)의 제한을 모터 지령에 추가하여 변환 행렬의 크기를 줄여야 한다(13).

(13)

$\Omega_{1, \text { ref }}=\Omega_{3, \text { ref }}, \Omega_{2, \text { ref }}=\Omega_{4, \text { ref }}$

$\Omega_{5, \text { ref }}=\Omega_{7, \text { ref }}, \Omega_{6, \text { ref }}=\Omega_{8, \text { ref }}$

$\mathbf{A}_{\mathrm{f}}=\left[\begin{array}{cccc}2 b & 2 b & 2 b & 2 b \\ b l & 0 & -b l & 0 \\ -b l & -\sqrt{2} b l & b l & \sqrt{2} b l \\ -2 d & 2 d & -2 d & 2 d\end{array}\right]$

$\Omega_{\text {sauare } 1,2,5,6, \text { ref }}=\mathrm{A}_{\mathrm{f}}^{-1} \mathrm{u}_{\mathrm{ref}}$

$A_{f}$는 식 (10)의 대해 식 (13)의 제약을 고려한 모터 속도를 추력으로 변환하는 행렬로, 기체의 지령 입력 $u_{ref}$를 8개의 모터의 지령 회전 속도 제곱 값으로 변환한다. 식 (13)의 속도 지령 정보가 모터 동역학과 하위제어기에 전달되면 안정된 비행을 수행하도록 모터가 동작한다. 시뮬레이션 상에서 모터의 동작을 나타내기 위해 PMSM과 배터리를 모델링 하였으며(14), PMSM은 sensorless(15)와 SVPWM신호(16)로 제어하였다.

2.4 Propeller fault

프로펠러는 모터의 회전을 기체의 추력으로 변환하는 주요 매개체로 프로펠러의 길이는 추력 T와 식 (14)의 관계를 가진다(17).

(14)
$\begin{aligned} T &=C_{\mathrm{T}} \rho D^{4} \Omega^{2} \\ &=b \Omega^{2} \text { where } b=C_{\mathrm{T}} \rho D^{4} \end{aligned}$

$C_{T}$은 추력의 비례 계수, 공기 밀도 그리고 D는 프로펠러 길이로 이 계수들을 추력 상수 와 밀접한 관계를 갖는데, (18)에서 시뮬레이션으로 유인 항공기의 엔진 터빈 블레이드 파손을 발생시키고 추력의 감소를 확인하였다. 이를 기반으로 본 논문은 프로펠러 파손을 추력 상수 b의 감소로 나타내어 시뮬레이션 상에서 고장을 모델링하고 상호 다중 모델을 적용하여 발생한 고장을 분류하였다.

3. Interacting Multiple Model(IMM)

Target 모델의 동역학이 변화하면 단일 필터로는 변환된 모델의 상태를 예측하기 제한된다(9). 하지만 복수의 필터 모델을 사용하는 상호 다중 모델 IMM을 적용하면 target 모델 출력과 유사한 예측을 수행하는 필터 모델에 가중치를 부여하여 단일 필터보다 정확한 예측을 수행할 수 있다(9-10).

그림 3은 IMM의 구조로, IMM은 Mixing, Filtering, Mode probability update 그리고 Merging 단계로 target 모델의 상태를 예측할 수 있다(10). 또한, 각 필터 모델의 가중치는 필터 모델과 target 모델의 유사도를 나타내는 mode probability로 target의 모드(정상 또는 고장)를 구분할 수 있다.

그림. 3. 상호 다중 모델 개요

Fig. 3. Introduction of Interacting Multiple Model

../../Resources/kiee/KIEE.2022.71.5.744/fig3.png

그림 3block에 나타나 있는 과정은 식 (15)부터 식 (19)로 나타내었다. Mixing 단계에서는 식 (15)와 같이 사전에 정의된 전이 확률 $π_{i|j}$와 각 필터 모드의 확률(mode probability) $u_{i,k-1}$을 결합하여 mixing probability $u_{i|j,k}$를 계산한다.

(15)

$\mu_{\mathrm{i} \mid \mathrm{j}, \mathrm{k}}=\frac{\pi_{\mathrm{i} \mid \mathrm{j}} \mu_{\mathrm{i}, \mathrm{k}-1}}{\bar{c}_{\mathrm{j}}}$

$\bar{c}_{\mathrm{j}}=\sum_{\mathrm{i}=0}^{\mathrm{N}=8} \pi_{\mathrm{i} \mid \mathrm{j}} \mu_{\mathrm{i}, \mathrm{k}-1}$

i와 j는 k-1과 k step의 필터 모델 지표(index)로 지표 0은 정상 모드 필터를 지표 1부터 8까지는 8개 프로펠러의 고장 모드 필터를 나타낸다. 식 (15)의 mixing probability 는 각 필터의 상호 관계를 정보이며, 식 (16)에서 이전 필터의 state $\hat{\mathbf{X}}_{\mathrm{k}-1}^{\mathrm{i}}$와 error covariance $\hat{\mathbf{P}}_{\mathrm{k}-1}^{\mathrm{i}}$를 활용해 k step 필터의 입력 state $\mathbf{X}_{\mathrm{k}-1}^{\mathrm{0j}}$와 error covariance $\mathbf{P}_{\mathrm{k}-1}^{\mathrm{0j}}$를 계산한다.

(16)

$\mathbf{x}_{\mathrm{k}-1}^{0 \mathrm{j}}=\sum_{\mathrm{i}=0}^{\mathrm{N}=8} \hat{\mathrm{x}}_{\mathrm{k}-1}^{\mathrm{i}} \mu_{\mathrm{i} \mid \mathrm{j}, \mathrm{k}}$

$\mathrm{P}_{\mathrm{k}-1}^{0 \mathrm{j}}=\sum_{\mathrm{i}=0}^{\mathrm{N}=8} \mu_{\mathrm{i} \mid \mathrm{j}, \mathrm{k}}\left[\hat{\mathrm{P}}_{\mathrm{k}-1}^{\mathrm{i}}+\right.$

$\left.\left(\hat{\mathrm{x}}_{\mathrm{k}-1}^{\mathrm{i}}-\mathbf{x}_{k-1}^{0 \mathrm{j}}\right)\left(\hat{\mathrm{x}}_{\mathrm{k}-1}^{\mathrm{i}}-\mathrm{x}_{\mathrm{k}-1}^{0 \mathrm{j}}\right)^{\mathrm{T}}\right]$

Filtering 단계에서는 식 (9), (10) 동역학이 고려된 정상 모드와 고장 모드의 필터들이 각 필터 모델을 기반으로 prediction과 correction을 식 (17)과 같이 수행하며, 기체의 모델이 비선형성을 가지기 때문에 각 필터는 EKF(Extended Kalman Filter)를 기반으로 동작한다.

(17)

Prediction

$x_{k}^{j}=A_{MAV, j} {x}_{k-1}^{0j}+B_{MAV,j} u_{k-1}$

$P_{k}^{j}=F_{MAV,j} P_{k-1}^{0j} F_{MAV,j}^{T}+Q_{MAV}$

Correction

$S_{j,k} = H_{MAV,j} P_{k}^{j} H_{MAV,j}^{T}+R_{MAV}$

$K_{j,k} = P_{k}^{j} H_{MAV,j}^{T} S_{j,k}^{-1} $

$\hat{x}_{k}^{j} = x_{k}^{j}+K_{j, k} (y_{k}-H_{MAV,j} x_{k}^{j})$

$\hat{P}_{k}^{j} = (I_{12 \times 12}-K_{j,k} H_{MAV, j}) P_{k}^{j}$

$Q_{MAV}$, $R_{MAV}$는 기체에 작용하는 시스템 노이즈 공분산, $A_{MAV,j}$, $B_{MAV,j}$는 식 (9)를 통해 유도된 기체 시스템 방정식, $F_{MAV,j}$, $H_{MAV,j}$는 기체 시스템과 출력의 Jacobian matrix로 각각 모델 지표 j에 따라 정상, 고장 모델 특성을 갖는다. 계산된 Kalman gain 로 $K_{j,k}$각 필터의 특성(정상 또는 고장)에 따라 상태 $\hat{x}_{k}^{j}$와 error covariance $\hat{P}_{k}^{j}$가 예측되며, 이 과정으로 계산되는 상태 예측값($H_{MAV,j}$ $x_{k}^{j}$)과 센서 error covariance($S_{j,k}$)는 식 (18)에서 각 필터와 target의 유사도를 나타내는 $A_{j,k}$를 계산하는 데 사용한다. 식 (18)에서는 $A_{j,k}$를 정규화하여 mode probability $u_{j,k}$를 update하고 이후 update $u_{j,k}$된 와 기체상태 임계값 ε을 비교하여 시스템의 고장 여부를 나타낸다.

(18)

$\Lambda_{\mathrm{j}, \mathrm{k}}=\frac{\exp \left(-\frac{1}{2} \mathrm{E}_{\mathrm{j}, \mathrm{k}}^{\mathrm{T}} \mathrm{S}_{\mathrm{j}, \mathrm{k}}^{-1} \mathrm{E}_{\mathrm{j}, \mathrm{k}}\right)}{\sqrt{\left|2 \pi \mathrm{S}_{\mathrm{j}, \mathrm{k}}\right|}}$

$c=\sum_{\mathrm{j}=0}^{\mathrm{N}=\mathrm{s}} \Lambda_{\mathrm{j}, \mathrm{k}} \bar{c}_{\mathrm{j}}, \quad \mu_{\mathrm{j}, \mathrm{k}}=\frac{\Lambda_{\mathrm{j}, \mathrm{k}} \bar{c}_{\mathrm{j}}}{c}$

$\mathrm{E}_{\mathrm{j}, \mathrm{k}}=\mathrm{y}_{\mathrm{k}}-\mathrm{H}_{\mathrm{MAV},{ }_{\mathrm{j}}} \mathrm{x}_{\mathrm{k}}^{\mathrm{j}}$

MAV system mode $= \begin{cases}\text { Propeller } j \text { fault } & \mu_{\mathrm{j}=1,2, \ldots,, \mathrm{s}, \mathrm{k}}>\epsilon \\ \text { Normal } & \text { other wise }\end{cases}$

마지막으로 식 (19) Merging에서는 식 (17)식 (18)의 상태 $\hat{x}_{k}^{j}$와 각 필터 mode probability $\hat{u}_{k}^{j}$를 결합하여 상태 예측 기댓값 $\hat{x}_{k}$를 연산한다. 다수의 필터 모델의 예측을 식 (19)로 결합한 상태 기댓값은 시스템 변동 시 단일 필터 예측값보다 더 높은 정확도를 보일 수 있다(10).

(19)
$\hat{\mathbf{X}}_{\mathrm{k}}=\sum_{\mathrm{j}=0}^{\mathrm{N}=8} \hat{\mathbf{x}}_{\mathrm{k}}^{\mathrm{j}} \mu_{\mathrm{j}, \mathrm{k}}$

본 논문에서는 IMM으로 프로펠러의 고장을 분류하기 위해 기체의 동역학과 추력 상수 b를 기반으로 1개의 정상 모드와 8개의 고장 모드 필터를 설계하고, 고장을 검출하였다.

4. Propeller fault simulation setting

4.1 MAV simulation setting

고장 분류 시뮬레이션은 MATLAB/Simulink에서 정상적으로 호버링비행 하는 MAV 기체에 고장을 발생시켰고, 호버링을 하기위한 기체, 제어기 이득은 표 1, 2의 값을 사용하였다. 표 1의 기체의 파라미터는 (13)을 참조하였으며, 제어 이득 표 2는 시뮬레이션 수행하여 조정하였다.

표 1. 정상 모드 기체 파라미터

Table 1. Normal mode MAV parameter

$J_{xx}$ : 0.0429

$J_{yy}$ : 0.0429

$J_{zz}$ : 0.0748

m : 2.371

$J_{r}$ : 1.945$e^{-4}$

$l$ : 1.0

$b_{normal}$ : 8.54858$e^{-6}$

$d$ : 1.205$e^{-6}$

표 2. 기체 상위/하위 제어기 이득

Table 2. MAV upper/lower controller gain

$K_{Dx}$ : -30.57

$K_{Dy}$ : -30.57

$K_{Dz}$ : -24

$K_{Px}$ : 36.25

$K_{Py}$ : -36.25

$K_{Pz}$ : 36

$K_{DΦ}$ : 20

$K_{DΘ}$ : 20

$K_{DΨ}$ : 20

$K_{PΦ}$ : 100

$K_{PΘ}$ : 100

$K_{PΨ}$ : 100

표1, 2 파라미터를 기반 기체의 정상 비행을 시뮬레이션으로 확인하였고, 시뮬레이션에 고장을 발생시켜 IMM의 고장 분류를 검증하였다.

4.2 Propeller fault and IMM setting

프로펠러는 기체에 추력을 발생시키는 주요 부품으로 프로펠러 파손은 기체의 추력 손실로 이어진다. 프로펠러 파손이 발생했을 때의 데이터를 수집하고 고장 분류 성능을 확인하기 위해 정상 호버링 중인 기체에 추력을 감소시키고 IMM을 사용하여 발생된 고장을 분류하였다. 프로펠러의 고장을 모사하기 위해 표 1의 추력 상수를 8개 프로펠러의 지표 p를 통해 $b_{p}$로 확장하였다. 본 논문에서는 1번째 프로펠러의 추력 상수를 변화시켜 고장 분류를 수행하였고 1번 외의 프로펠러 추력 b는 식 (20)과 같이 정상으로 설정하였다.

(20)

$b_{\mathrm{D}} \in\{2,3, \cdots, 8\}=b_{\text {normal }}$

$b_{1}=F_{\mathrm{k}} b_{\text {normal }}$

$F_{\mathrm{k}} \in\left\{\begin{array}{llllllll}1 & 0.9 & 0.7 & 0.5 & 0.3 & 0.1 & 0\end{array}\right\}$

프로펠러의 정상부터 탈조까지 모사하기 위해 고장의 지표 $F_{K}$를 1(정상)부터 0(완전 파손)까지 변화시켜 시뮬레이션을 수행하였고 각각은 추력의 0~100%감소를 나타낸다. 이후 고장이 발생한 프로펠러를 구분하기 위해 그림 3 IMM의 필터들의 추력 상수 b를 식 (21)로 설정하면 그림 4의 과정으로 기체의 고장을 분류할 수 있다.

(21)

Filter 0: $\quad b_{p \in\{1,2, \cdots, 8\}}=b_{\text {normal }}$

Filter 1: $b_{1}=0.5 b_{\text {normal }}, \quad b_{\mathrm{p}} \in\{2,4, \cdots, 8\}=b_{\text {normal }}$

Filter 2: $b_{2}=0.5 b_{\text {normal }}, \quad b_{p} \in\{1,3, \cdots, 8\}=b_{\text {normal }}$

$\vdots$

Filter 8: $b_{8}=0.5 b_{\text {normal }}, \quad b_{\mathrm{D} \in\{1,2, \cdots, 7\}}=b_{\text {normal }}$

그림. 4. 기체 고장 분류 과정

Fig. 4. MAV Fault detection and isolation Process

../../Resources/kiee/KIEE.2022.71.5.744/fig4.png

식 (21)은 정상과 고장이 발생한 프로펠러를 구분하기 위해 설계된 9개 필터의 추력 상수이다. 고장이 발생한 프로펠러를 분류하기 위해 각 프로펠러 지표 p와 연동되는 정상 추력 상수 를 낮춰 각 모드의 필터를 모델링 하였고, p=0을 제외한 식 (21)의 Filter p는 p번째 프로펠러 추력 상수만 감소 된 모델을 의미한다.

식 (21)의 필터 고장을 50%로 설정하여 고장 필터 모델을 설계한 이유는 효율적으로 프로펠러의 고장을 분류하기 위함이다. 만일 탈조 상태 즉, 추력 감소 100%로 고장 필터 추력을 설정하면 추력이 30%감소한 고장이 발생할 때 탈조 상태보다 정상 상태에 가까운 결과를 나타내어 정상으로 인식된다. 이와 같은 이유로 정상과 탈조의 중간인 0.5$b_{normal}$로 고장 필터의 추력 상수를 설정하여 고장을 검출하였다. 마지막으로 상호 다중 모델을 동작시키기 위해 전이 확률$π_{i|j}$를 설정해야 하며 본 논문에서는 그림 5과 같은 전의 확률을 활용했다.

그림. 5. 상호 다중 모델 전이 확률

Fig. 5. Interacting Multiple model transition probability

../../Resources/kiee/KIEE.2022.71.5.744/fig5.png

전이 확률은 각 필터 모델의 상호 관계를 나타내며, 고장은 정상에서 발생하고 고장이 발생한 뒤 다른 고장으로 전환될 확률은 없다는 가정을 그림 5의 전이 확률을 통해 나타냈다.

5. Simulation result

5.1 Simulation verification

기체에 고장을 발생시키고 고장을 분류하기 전에 시뮬레이션이 올바르게 설계되었는지 검증해야 한다. 본 논문에서는 기체의 정상, 고장 모드에서 원형 경로와 호버링 비행을 그림 6, 7, 8과 같이 실행해 보고 각 경우의 기체 데이터를 확인함으로써 시뮬레이션을 검증을 수행하였다.

그림. 6. 정상 기체의 원형 경로 비행

Fig. 6. Normal mode MAV circular trajectory flying

../../Resources/kiee/KIEE.2022.71.5.744/fig6.png

그림. 7. 정상과 고장 기체 경로

Fig. 7. Normal mode fault MAV trajectory

../../Resources/kiee/KIEE.2022.71.5.744/fig7.png

그림 6은 정상 모드에서 기체의 동작을 확인하기 위해 10m 직경의 원형 경로 비행 시뮬레이션을 실행시킨 결과이다. 초기 시작점에서 원형 경로에 도달하는 과정에서 경로를 벗어나지만, 경로에 도달한 후에는 이탈하지 않고 비행하는 것을 확인할 수 있다. 이를 통해 정상 기체가 원형 경로를 비행할 수 있도록 정상 모드 시뮬레이션이 올바르게 설계되었음을 확인할 수 있다.

그림 7은 정상 기체 호버링 비행 중에 8개 프로펠러에 고장을 발생시킨 후 경로 데이터이며 각 프로펠러 고장 시 기체의 동작을 분석함으로써 시뮬레이션을 검증할 수 있다. 그림 7의 라벨 Fault의 p는 8개 프로펠러 중 고장이 발생한 프로펠러 지표이며 감소된 추력의 크기는 각 프로펠러 모두 같고 50 감소로 설정하였다.

8개의 프로펠러 고장 중 프로펠러 1번의 고장 경로를 확인해 보면 +Y로 이동하는 것을 확인할 수 있는데 이는 +Y를 지향하는 1번 프로펠러가 추력이 감소하여 기체가 기울어졌기 때문이다. 이와 같은 현상은 다른 프로펠러에서도 동일하게 발생하며 이 중 1번과 마주 보는 위치에 설치된 5번 프로펠러 고장 경로는 –Y로 이동하는 것을 확인할 수 있다. 이러한 대칭적인 비행경로는 1번 5번 외에 다른 마주 보는 프로펠러 고장에서도 확인 가능하며 시뮬레이션 설계를 검증해준다.

마지막으로 고장의 크기가 변화할 때의 시뮬레이션 데이터를 비교해 보기 위해 1개의 프로펠러의 고장 수준을 변화시켜 시뮬레이션을 실행시켰다. 시뮬레이션을 검증하기 위한 데이터 중 고장이 발생한 프로펠러의 회전 속도를 비교해 보았는데 이는 프로펠러 추력이 감소할수록 모터가 고장을 보상하기 위해 더 빠른 회전을 수행하는 것을 관측하기 위함이다.

그림. 8. 고장 변화에 따른 1번 프로펠러 속도 데이터

Fig. 8. 1st propeller speed data depends on level of fault

../../Resources/kiee/KIEE.2022.71.5.744/fig8.png

그림 8은 1번 프로펠러 추력이 정상부터 70%까지 감소한 데이터이다. 70%보다 추력이 감소한 고장은 모터의 속도가 800[rad/s]로 제한되어있으므로 70%고장 발생 후 모터 속도와 거의 동일한 결과를 보여 생략했다.

고장이 발생한 시간인 2초 이후 데이터를 확인해 보면 정상부터 고장의 심각할수록 모터의 제어기에서 더 빠른 회전을 수행하여 감소한 추력을 보상하려는 동작을 관측할 수 있다. 이를 통해 제어기가 고장을 보상하기 위해 동작하는 것을 확인할 수 있지만, 모터가 추력을 보상하기 위해 더 빠른 회전을 수행해도 그림 7과 같이 기체의 편향된 동작이 관측된다. 이는 PD 제어기로는 프로펠러가 파손된 기체를 제어할 수 없으므로 고장 감지, 분류 후 적절한 대응을 수행할 수 있는 알고리즘의 필요성을 나타내며, 본 논문은 이를 개발하기 위한 첫 단계로 고장 분류알고리즘을 연구하였다.

5.2 Fault detection algorithm verification

그림 6, 7, 8의 검증된 기체 기반으로 1번 프로펠러에 대해 고장을 발생시켜 보았으며 그림 9-a부터 9-f까지 고장의 크기를 변화시켜 가며 IMM 고장 분류를 확인한 결과이다.

그림. 9-a. 고장 진단 확률, 1번 프로펠러 추력 10 감소

Fig. 9-a. Fault diagnose probability 1st propeller thrust 10 decrease

../../Resources/kiee/KIEE.2022.71.5.744/fig9_a.png

그림. 9-b. 고장 진단 확률, 1번 프로펠러 추력 30 감소

Fig. 9-b. Fault diagnose probability 1st propeller thrust 30 decrease

../../Resources/kiee/KIEE.2022.71.5.744/fig9_b.png

그림. 9-c. 고장 진단 확률, 1번 프로펠러 추력 50 감소

Fig. 9-c. Fault diagnose probability 1st propeller thrust 50 decrease

../../Resources/kiee/KIEE.2022.71.5.744/fig9_c.png

그림. 9-d. 고장 진단 확률, 1번 프로펠러 추력 70 감소

Fig. 9-d. Fault diagnose probability 1st propeller thrust 70 decrease

../../Resources/kiee/KIEE.2022.71.5.744/fig9_d.png

그림. 9-e. 고장 진단 확률, 1번 프로펠러 추력 90 감소

Fig. 9-e. Fault diagnose probability 1st propeller thrust 90 decrease

../../Resources/kiee/KIEE.2022.71.5.744/fig9_e.png

그림. 9-f. 고장 진단 확률, 1번 프로펠러 추력 100 감소

Fig. 9-f. Fault diagnose probability 1st propeller thrust 100 decrease

../../Resources/kiee/KIEE.2022.71.5.744/fig9_f.png

안정된 호버링 비행을 수행하고 있는 기체 프로펠러에 고장이 발생하였을 때 IMM의 고장 분류를 그림 9를 통해 확인했으며 각 데이터의 라벨링은 그림 7과 같다. 초기 IMM의 mode probability는 9개의 필터를 고려하여 필터마다 1/9로 설정하였고 고장을 인식하기 위한 식 (18)의 기체상태 임계값 ε을 0.65로 설정하였다.

시뮬레이션의 고장이 발생하기 2초 전에는 정상 비행이 관측되기 때문에 정상 필터의 mode probability가 상승한다. 하지만 프로펠러 1번에 고장이 인가된 2초 이후에는 고장의 크기에 따라 mode probability 값에 변화가 발생하며, 고장 인가 후 고장 인식 시간을 표 3을 통해 나타내었다.

표 3. 그림 9에 대한 고장 감지 시간

Table 3. For figure 9, 1st propeller fault detection time

9-a

9-b

9-c

9-d

9-e

9-f

time [s]

-

1.74

0.38

0.21

0.16

0.14

기체에 고장이 가장 작게 발생한 9-a는 고장을 인지하지 못하고, 추력이 30% 감소한 9-b는 1.74초 9-c에서 9-f까지의 고장은 0.5초 이내의 빠른 고장 인지 성능을 보인다. 이러한 차이는 IMM이 식 (18) 정상, 고장 필터 모델 예측과 고장 기체 동작의 유사도를 기반으로 mode probability를 update하기 때문이다. 그림 9-a 10% 추력 감소를 고장으로 인식하지 못한 이유는 고장이 발생해도 기체의 움직임이 정상 모드와 유사도가 더 높아 정상모드로 분류된 것이며, 9-b부터 9-f까지는 정상 모드보다 고장 모드에 유사도가 더 높아 고장으로 인식된 것이다. 하지만 9-b는 고장 인식 과정에서 상대적으로 지연된 결과를 나타내는데 이는 추력 30%감소 고장이 다른 심각한 고장보다 더 높은 정상 필터 유사도를 가지기 때문이다. 다른 고장보다 높은 정상 모드 유사도는 식 (18)의 mode probability 정규화 과정에서 고장 모드 가중치 update를 지연시켜 그림 9-b와 같은 지연된 고장 인식 결과를 나타낸다. 이에 대해 현재 식 (18)의 임계값과 식 (21)의 정상, 고장 필터 모델로 구분할 수 있는 고장 범위를 그림 1011로 확인해 보았다.

그림. 10. 고장 진단 확률, 1번 프로펠러 추력 24% 감소

Fig. 10. Fault diagnose probability 1st propeller thrust 24% decrease

../../Resources/kiee/KIEE.2022.71.5.744/fig10.png

그림. 11. 고장 진단 확률, 1번 프로펠러 추력 23% 감소

Fig. 11. Fault diagnose probability 1st propeller thrust 23% decrease

../../Resources/kiee/KIEE.2022.71.5.744/fig11.png

그림 1011은 추력 감소가 24%, 23% 발생했을 때의 고장 진단확률이며 고장 발생 시점과 같은 시뮬레이션 설정은 그림 9와 같지만 보다 지연된 반응을 보여 30초 동안 확률 변화를 관측했다. 그림 10의 고장 감지에 소요된 시간은 9.57초이며 이는 그림 9-b처럼 정상 모드와의 유사도가 높아 발생한 결과이다. 그림 11은 고장 모드의 확률의 반전이 발생하는 것이 확인되지만 그 크기가 식 (18)에서 지정한 고장 임계값보다 낮아 고장으로 인식하지 못했다. 이를 통해 본 논문에서 사용한 고장 임계값과 고장 필터 모델은 추력이 24% 감소한 고장에 대해 분류할 수 있음을 확인했다. 비록 추력이 10%, 23% 감소한 고장들은 본 논문에서 구분하지 못했지만, 식 (18)의 고장 임계값과 식 (21) 고장 필터 모델을 조절하면 IMM으로 고장을 분류할 수 있다.

6. Conclusion

MAV 기체에 모터와 배터리를 세부적으로 모델링 함으로써 더 현실적인 기체의 상태를 관측하는 시뮬레이션을 설계하였다. 설계된 시뮬레이션을 검증하기 위해 정상모드일 때 원형 경로 그리고 고장 모드의 호버링 비행을 수행해보고 원형 경로 추적 능력과 고장 발생 후 기체의 편향된 상태 그리고 대칭성이 관측됨으로써 유효한 기체 동역학 시뮬레이션이 설계되었음을 확인하였다. 이러한 시뮬레이션은 기체에 적용할 알고리즘의 성능을 실 테스트 전에 확인할 수 있어 더 안전한 알고리즘 개발을 할 수 있다. 본 논문에서는 시뮬레이션에 프로펠러 고장을 인가하고 IMM 고장 분류알고리즘을 유효성을 확인해 보았다. 고장 분류알고리즘은 고장을 인지하고 발생한 고장을 분류하여 고장에 강인한 제어알고리즘 개발의 첫 단계가 되며 IMM은 다양한 수준의 고장이 발생하였을 때 대부분 빠르고 정확한 고장 인지 및 분류 성능을 보였다. 하지만 발생한 고장이 정상 필터 모델과 많은 유사도를 보일수록 고장 인식에 지연이 발생하고 고장 필터보다 정상 필터 모델에 더 높은 유사도를 나타내면 고장 인식 기능이 제한된다. 이러한 문제는 현실의 고장 데이터를 취득 후 적절한 고장 필터 모델과 고장 임계값을 설정함으로써 개선할 수 있는 문제이기 때문에 IMM고장 분류는 높은 신뢰성을 가지는 결과를 나타낼 수 있다.

향후 복수의 프로펠러 고장 분류를 수행할 수 있는 고장 필터 모델을 연구할 것이며, 또한 시뮬레이션으로 수집된 고장 데이터로 LSTM, CNN(Convolution Neural Network) 그리고 SVM과 같은 데이터 기반 고장 분류알고리즘을 분석하고 IMM과 비교해 봄으로써 더 신뢰성 있는 고장 분류알고리즘 개발을 수행할 것이다.

7. Appendix

7.1 Skew matrix

Skew matrix는 기체의 동역학을 해석할 때 사용하는 유용한 tool로 식 (22)과 같은 특성이 있으며, 첨자 E, B는 관성 좌표계와 기체 좌표계의 상태를 나타낸다.

(22)

$\mathrm{S}^{\mathrm{T}}\left(\omega^{\mathrm{B}}\right)+S\left(\omega^{\mathrm{B}}\right)=0$

$\omega^{\mathrm{B}}=[p, q, r]^{\mathrm{T}}, \mathrm{S}\left(\omega^{\mathrm{B}}\right)=\left[\begin{array}{ccc}0 & -r & q \\ r & 0 & -p \\ -q & p & 0\end{array}\right]$

$\dot{\mathrm{R}}_{\mathrm{BE}}=\mathrm{S}\left(\omega^{\mathrm{E}}\right) \mathrm{R}_{\mathrm{BE}}$

$\mathrm{S}\left(\omega^{\mathrm{B}}\right) \omega^{\mathrm{B}}=\omega^{\mathrm{B}} \times \omega^{\mathrm{B}}$

$\dot{\mathrm{R}}_{\mathrm{BE}}=\mathrm{R}_{\mathrm{BE}} \mathrm{S}\left(\omega^{\mathrm{B}}\right) \mathrm{R}_{\mathrm{BE}}^{\mathrm{T}}=\mathrm{S}\left(\boldsymbol{R}_{\mathrm{BE}} \omega^{\mathrm{B}}\right)$

7.2 MAV rigid body linear motion model derive

기체의 선형 동역학은 식 (22)의 skew matrix의 특성을 이용하여 식 (23)와 같이 순차적으로 풀어 유도할 수 있다.

(23)

$\mathrm{P}^{E}=[X Y Z]^{\mathrm{T}}, \omega^{\mathrm{B}}=[p, q, r]^{\mathrm{T}}$

$\ddot{\mathrm{P}}^{\mathrm{E}}=\frac{d}{d t} \dot{\boldsymbol{P}}^{\mathrm{E}}=\frac{d}{d t}\left(\boldsymbol{R}_{\mathrm{BE}} \dot{\boldsymbol{P}}^{\mathrm{B}}\right)$

$\ddot{\mathrm{P}}^{\mathrm{E}}=\dot{\mathrm{R}}_{\mathrm{BE}} \dot{\mathrm{P}}^{\mathrm{B}}+\mathrm{R}_{\mathrm{BE}} \ddot{\mathrm{P}}^{\mathrm{B}}$

$\ddot{\mathrm{P}}^{\mathrm{E}}=\mathrm{S}\left(\omega^{\mathrm{E}}\right) \mathrm{R}_{\mathrm{BE}} \dot{\mathrm{P}}^{\mathrm{B}}+\mathrm{R}_{\mathrm{BE}} \ddot{\mathrm{P}}^{\mathrm{B}}$

$\ddot{\mathrm{P}}^{\mathrm{E}}=\mathrm{S}\left(\boldsymbol{R}_{\mathrm{BE}} \omega^{\mathrm{B}}\right) \mathrm{R}_{\mathrm{BE}} \dot{\mathrm{P}}^{\mathrm{B}}+\mathrm{R}_{\mathrm{BE}} \ddot{\mathrm{P}}^{\mathrm{B}}$

$\ddot{\mathrm{P}}^{\mathrm{E}}=\mathrm{R}_{B E}^{\mathrm{E}} S\left(\omega^{\mathrm{B}}\right) \mathrm{R}_{B E}^{\mathrm{T}} \mathrm{R}_{\mathrm{BE}} \dot{\mathrm{P}}^{\mathrm{B}}+\mathrm{R}_{\mathrm{BE}} \ddot{\mathrm{P}}^{\mathrm{B}}$

$\ddot{\mathrm{P}}^{\mathrm{E}}=\mathrm{R}_{\mathrm{BE}}\left(\mathrm{S}\left(\omega^{\mathrm{B}} \dot{\mathrm{P}}^{\mathrm{B}}\right)+\ddot{\mathrm{P}}\right)$

$\ddot{\mathrm{P}}^{\mathrm{E}}=\mathrm{R}_{\mathrm{BE}}\left(\dot{\mathrm{V}}^{\mathrm{B}}+\omega^{\mathrm{B}} \times \mathrm{V}^{\mathrm{B}}\right)$

7.3 MAV rigid body linear motion model derive

회전 동역학 또한 기체의 병진 운동과 마찬가지로 식 (24)의 순차적인 과정으로 회전력과 기체의 회전 상태와의 관계를 해석할 수 있다.

(24)
$\begin{aligned} \ddot{\boldsymbol{\theta}}^{\mathrm{E}} &=\frac{d}{d t} \dot{\boldsymbol{\theta}}^{\mathrm{E}}=\frac{d}{d t}\left(\mathrm{~T}_{\mathrm{BE}} \omega^{\mathrm{B}}\right) \\ \ddot{\boldsymbol{\theta}}^{\mathrm{E}} &=\left(\dot{\mathrm{T}}_{\mathrm{BE}} \omega^{\mathrm{B}}\right)+\left(\boldsymbol{T}_{\mathrm{BE}} \dot{\omega}^{\mathrm{B}}\right) \\ \ddot{\boldsymbol{\theta}}^{\mathrm{E}} &=\mathrm{T}_{\mathrm{BE}}\left(\dot{\omega}^{\mathrm{B}}+\omega^{\mathrm{B}} \times \boldsymbol{\omega}^{\mathrm{B}}\right) \end{aligned}$

Acknowledgements

This work was supported in part by the National Research

Foundation of Korea (NRF) grant funded by the Korea government (MSIT) under Grant NRF-2020R1F1A1071547.

References

1 
U. R. Mogili, B. B. V. L. Deep, 2018, Review on Application of Drone Systems in Precision Agriculture, Procedia Comput. Sci, Vol. 133, pp. 502-508DOI
2 
E. Frachtenberg, 2019, Practical Drone Delivery, Comput, Vol. 52, No. 15, pp. 53-57DOI
3 
Y. Granesh, R. Raju, R. Hegde, 2015, Surveilance Drone for Landmin Detection, 2015 Interanational Conference on Advanced Computing and Communications(ADCOM), pp. 33-38DOI
4 
J. Cui, 2015, Drones for cooperative search and rescue in pos-disaster situation, 2015 IEEE 7th International Conference on Cybernetics and Intelligent Systems (CIS) and IEEE in Robotics Automation and Mechatronics (RAM), Vol. , No. , pp. 167-174DOI
5 
I. Atoui, H. Meradi, R. Boulkroune, R. Saidi, A. Grid, 2013, Fault detection and diagnosis in rotating machinery by vibration monitoring using FFT and Wavelet techniques, 2013 8th International Workshop on Systems Singnal Processing and their Application(WoSSPA), pp. 401-406DOI
6 
I. Atoui, H. Meradi, R. Boulkroune, R. Saidi, A. Grid, 2013, Fault detection and diagnosis in rotating machinery by vibration monitoring using FFT and Wavelet techniques, 2013 8th International Workshop on Systems Singnal Processing and their Application(WoSSPA), pp. 401-406DOI
7 
B. Wang, M. Jiang, D. Liu, 2020, Real-Time Fault Detection for UAV Based on Model Acceleration Engine, IEEE Transactions on Instrumentation and Measurement, Vol. 69, No. 12, pp. 9505-9516DOI
8 
M. H Amoozgar, A. Chmseddine, Y. Zhang, 2013, Experimental Test of a Two-stage Kalman Filter for Actuactor Fault Detection and Diagnosis of Unmmaned Quadrotor Helicopter, J Intell Robot Syst, Vol. 70, pp. 107-117DOI
9 
L. Cork, R. Walker, 2007, Sensor Fault Detection for UAVs using a Nonlinear Dynamic Model and the IMM-UKF Algorithm, Information, decision and control, pp. 230-235DOI
10 
C. Jeong, C. Park, C. M. Kang, 2021, Fault Diagnosis Algorithm Using Parallel Connected Interacting Mutiple Model for Electric Power Steering Systems, Transactions of KSAE, Vol. 29, No. 10, pp. 943-950Google Search
11 
H. A. P. Blom, Y. Bar-Shalom, 1988, The Interacting Mutiple Model Algorithm for Systems with Markovian Switching Coefficients, IEEE Transaction on Automatic Control, Vol. 33, No. 8DOI
12 
S. H. Jeong, W. S Lee, Y. S. Kang, 2013, Neighboring Vehicle Maneuver Detection using IMM Algorithm for ADAS, Journal of Institute of Control Robotics and systems, Vol. 19, No. 8, pp. 718-724DOI
13 
N. Osmić, M. Kurić, I. Petrović, 2016, Detailed octorotor modeling and PD control, 2016 IEEE International conference on Systems Man and Cybernetics(SMC), pp. 2182-2189DOI
14 
G. L. Plett, 2015, Battery Management Systems, Volume I: Battery Modeling. NorwoodGoogle Search
15 
H. Nak, M. O. Gűlbahce. M. Gokasan, A. F. Ergence, 2015, Performance investigation of extended Kalman filter based observer for PMSM using in washing machine applications, 2015 9th International Conference on Electrical and Electronics Engineering(EKLECO), pp. 618-623DOI
16 
K. V. Kumar, P. A. Micahel, J. P. John, S. S. Kumar, 2010, Simulation and comparison SPWM and SVPWM control for three phase inverter, ARRN Journal of Engineering and Applied Sciences, Vol. 5, No. 7, pp. 61-74Google Search
17 
E. Kuantama, D. Craciun, R. Trace, 2016, Quadcopter Propeller Design and Performance Analysis, In New Advances in Mechanisms, Mechanical Transmissions and Robotics. Mechanisms and Machine Science; Springer: Cham, Vol. 46, pp. 269-277DOI
18 
H. Liu, M. H. Che Man, K. H. Low, 2021, UAV airborne collision to manned aircraft engine: Damage of fan blades and resultant thrust loss, Aerospace Science and Technology, Vol. 113DOI

저자소개

박채훈(Chaehun Park)
../../Resources/kiee/KIEE.2022.71.5.744/au1.png

He received the B.S degree in Electrical Engineering from Incheon National University, Incheon, South Korea in 2021.

Currently, he is pursuing the M.S degree in the same institution.

His research interests include FDI, model predictive control and reinforcement learning.

정철민(Cheolmin Jeong)
../../Resources/kiee/KIEE.2022.71.5.744/au2.png

He received the B.S degree in Electrical Engineering from Incheon National University, Incheon, South Korea in 2020

Currently, he is pursuing the M.S degree in the same institution.

His research interest include optimal control, linear control, tracking, system FDI, SLAM, autonomous vehicle, autonomous robot, reinforcement learning

강창묵(Chang Mook Kang)
../../Resources/kiee/KIEE.2022.71.5.744/au3.png

He received the B.S and Ph.D degrees in Electrical Engineering from Hanyang University, Seoul, South Korea, in 2012 and 2018 respectively.

He was a Senior Engineer with Agency for Defense Development, Daejeon, Korea, from 2018 to 2019.

He is currently an Assistant Professor with the Department of Electric Engineering, Incheon National University, Incheon, South Korea.

His research interest include linear system, optimal control. autonomous vehicle and artificial intelligent.