部分フィードバック線形化を用いた無人航空機のモデル予測型着 陸制御
○周楊 大橋あすか 鷹羽浄嗣 （立命館大学）
Model Predictive Landing Control of an Unmanned Aerial Vehicle via Partial Feedback Linearization
∗ Y. Zhou, A. Ohashi and K. Takaba (Ritsumeikan Univ.)
Abstract–This paper introduces a model predictive control approach of an unmanned aerial vehicle (UAV) with the aid of a feedback linearization. As is well known, the feedback linearization is one of the effective techniques to cope with the nonlinearity of dynamical systems. Since the UAV is an underactuated nonlinear system, it is impossible to exactly linearize the dynamics of the UAV. Therefore, we take an approach to linearize only the translational motion, and then apply the linear optimal control to it. However, the UAV is easily affected by wind disturbances in an actual environment. A model predictive control is proposed to cope with the disturbances. We apply this approach to a landing control of the UAV to a moving ground vehicle. The effectiveness of the proposed method is verified by numerical simulations.
Key Words:Unmanned aerial vehicle, Model predictive control, Feedback linearization, Landing control
1 はじめに
近年，クワッドコプターやマルチコプターなどの小 型無人航空機(Unmanned Aerial Vehicle(UAV))が，災 害探索，宅配，セキュリティ監視及び空撮などの応用 で活躍されている．UAVの自律制御あるいは自動制御 は， これらの応用に不可欠な技術の一つである．
UAVの自動制御技術に関しては，モデル予測制御 (Model Predictive Control(MPC))[1]，最小時間制御[2]
などいくつかの最適化ベースの手法が提案されている．
Villarreal, Rossiter, and Shin [1]は，適応型ラゲールベー スのモデル予測制御をUAVの姿勢安定化に適用した．
Ritzet al.[2]は，最小時間制御をUAVの軌道生成問題 に適用した．
将来のタスクの複雑さを考慮すると，UAVと地上車 両(Ground Vehicle(GV))の協調制御が重要である．協 調制御システム設計の一つの例として，移動GVに対 するUAVの自律着陸制御問題が挙げられる．この問題 には，二つの課題がある．一つは，非線形劣駆動シス テムに対処することである．よく知られているように，
フィードバック線形化は動的システムの非線形性に対処 する有効な手法の一つである．UAVは劣駆動システム であるため，厳密なフィードバック線形化をUAVのダ イナミクスに適用することが不可能である．もう一つ は，UAVをGVに安全に着陸させることである．Zhou
and Takaba[3]は，最適な着陸軌道を得るためには，部
分的なフィードバック線形化によるUAVの最適制御を 提案した．ただし，UAVは風の影響を受けやすいため，
これらの外乱に対処できる制御手法が必要である．よ く知られているように，MPCは外乱に対処する有効な
手法の一つである．
MPCは制御対象のモデルを用いて将来の動作を予測 する高度な最適化ベースの制御手法である．MPC法に 基づくUAVのいくつかの研究がある．例えば，高い計 算コストを必要とする非線形MPC[4]と近似線形化を
用いるMPC[5]などがある．
本稿では，MPCベースの部分フィードバック線形化 手法をUAVに適用した．この手法は，UAVの着陸軌 道をオンラインで修正することで，風の外乱による悪 影響を低減させることができる．また，部分フィード バック線形化を用いているため，非線形MPCと近似線 形化によるMPCと比べ，それぞれ計算コストと線形化 誤差が減少できる．さらに，この制御手法を三次元運 動に拡張した．最後に，提案した手法の有効性を数値 シミュレーションによる検証する．
2 問題定式化
Fig. 1は着陸制御タスクの概念図を示している．所
望の進入降下角を維持し，UAVは移動GVの上に着陸 するはずである．
Fig. 1: Landing Control of the UAV to the Moving GV
2B2-5
2.1 モデリング 2.1.1 座標システム
Fig. 2はUAVの座標系を示している．地上座標系は
地面に固定されたマーカーの重心を原点とする座標系 である．北向きをX軸，西向きをY 軸，地球の中心か ら外向きをZ 軸とする．機体座標系はUAVの重心を 原点とする座標系である．機体前方をx軸，左向きを y軸，上向きをz軸とする．回転角ロール，ピッチ及 びヨーはそれぞれϕ，θ，ψで表される．機体座標系を 地上座標系に変換する回転行列R(η)は次式になる．
R(η) =R(ψ)R(θ)R(ϕ)
ここで，R(ψ)，R(θ)及びR(ϕ)はそれそれx，y及びz 軸を中心に回転する回転行列である．
Fig. 2: Coordinate Systems
Fig. 3: Geometry of the Landing Control Problem 2.1.2 GVモデル
Fig. 3 では，GV が X 軸に対して角度 κの方向
に一定の速度 VT で移動すると仮定する．XB(t) = (XT(t), YT(t))をGVの位置，hをGVの高さ，XB0
をGVの初期位置とする．GVのモデルは次のように 表される．
X˙B(t) =VT，XB0 =XB(0) 2.1.3 UAVモデル
UAVの運動方程式[6]は次式で与えられる．
mξ¨+mgEZ =R(η)EZu Jη¨ =τ−C(η,η) ˙˙ η
ここで，uはUAVの機体座標系のz方向の総推力であ り，ほかの記号の定義はTable 1に示す．
Table 1:物理量とパラメータ
m UAVの質量
ξ= (X, Y, Z)T 地上座標系における位置 g 重力加速度g= 9.81[m/s2] EZ= (0,0,1)T Z方向の単位ベクトル
η= (ϕ, θ, ψ)T ロール角，ピッチ角，ヨー角 J=diag(Ix, Iy, Iz) 慣性モーメント行列
τ= (τϕ, τθ, τψ)T 揚力によるトルク
C 航空力学トルク
運動方程式は次のように書き換えることができる．
mX¨ =u(sinϕsinψ+ cosϕcosψsinθ) mY¨ =u(cosϕsinθsinψ−cosψsinϕ) mZ¨ =ucosθcosϕ−mg
Ixϕ¨=τϕ
Iyθ¨=τθ
Izψ¨=τψ
ここで，τϕ，τθ及びτψはトルク入力を表す．上記の方 程式をu1=u/m，u2 =τϕ/Ix，u3 =τθ/Iy，u4 = τψ/Izで正規化すると，次式を得る．
X¨ =u1(sinϕsinψ+ cosϕcosψsinθ) (1) Y¨ =u1(cosϕsinθsinψ−cosψsinϕ) (2) Z¨ =u1cosθcosϕ−g (3)
ϕ¨=u2 θ¨=u3 ψ¨=u4
2.2 問題設定
モデリングに基づき，最適着陸制御問題が次のよう に定式化する．その幾何図形はFig. 3に示す．
2.2.1 GVに対するUAVの進入降下角
UAVには前カメラが内蔵されていると想定する．安 全や監視のため，GVをUAVの前カメラの視野内に保 つことが望まれる．UAVの進入降下角を所望の値に保 つことにより，この目的を達成したいと考えている．定 数γを所望の進入降下角とし，V = [−VX −VY δVT] を所望の進入降下角に対応する軌道の方向ベクトルと し，F = [X−XT Y −YT Z−h]をUAVとGVの相 対位置ベクトルとする．実際の進入降下角がγと完全 に等しいとき，V とF は平行になる．従って，次の関
係式が得られる．
|V ×F|=|V||F|sinρ= 0
ここで，δ= tanγ，ρはV とF のなす角である.従っ て，|V ×F|2を十分に小さくすることで，進入降下角 をγに近づけることができる．
2.2.2 終端制約
着陸制御の最後に，UAVの運動がGVの運動と一致 する必要があるため，終端制約は次のように表される．
X(L) =XT(L)，X(L) =˙ VX，X¨(L) = 0，X(3)(L) = 0 Y(L) =YT(L)，Y˙(L) =VY，X(L) = 0，¨ X(3)(L) = 0 Z(L) =h，Z(L) = 0，˙ Z¨(L) = 0，Z(3)(L) = 0
ψ=κ，ψ˙= 0 (4)
ここで，Lは終端時刻である. 2.2.3 風外乱
実際の状況では，UAVは風の影響を受けやすいため，
その悪影響を低減させることが重要である．
本稿における制御目標は，風外乱が存在する場合で も進入降下角を所望の値に近づけながら，UAVを移動 GVの上に着陸させることである．この目的を達成す るために，モデル予測制御(MPC)手法を適用する．
3 コントローラー設計
MPCの主な考え方は，特定の有限時間最適制御問題 を繰り返し解き，各時刻で得られた最適な入力列の最 初の値を適用する．最適制御を計算する際に，MPCは 制御対象のモデルを用いてシステムの未来の動作を予 測する．MPC法により，UAVの着陸軌道をオンライ ンで修正し，風外乱による悪影響を低減することがで きる．
有限ホライゾンの最適制御問題の直接的な解法には，
計算的に要求が厳しい反復非線形最適化が含まれるこ とと，線形近似と反復の早期終了による数値誤差が含 まれることに注意されたい．本稿では，MPCと部分 フィードバック線形化を組み合わせ，近似誤差なしの 最適制御を得るための簡単なアルゴリズムを示す．
3.1 部分フィードバック線形化
フィードバック線形化は動的システムの非線形性に 対処する有効な手法である．なぜなら，この手法は非 線形ダイナミクスを近似なしの線形ダイナミクスに変 換するからである[7]．ただし，UAVは劣駆動システ ム，すなわち，システムの自由度に対して入力の数が
少ないシステムであるため，全ての状態変数を線形化 することができない．
本稿では，[8]に沿って線形化されたモデルを採用す る．ここで，並進運動が非線形フィードバックにより，
線形化される．また，u1の出現を式(1)， 式(2)及び式 (3)の高次導関数まで遅らせ，ほかの入力が現れること を期待する必要がある．ここで，u1をu¯1による駆動 される二重積分器の出力に等しく設定する．
¨ u1= ¯u1
表記の一貫性のために，ほかの入力も次のように設定 する．
u2= ¯u2
u3= ¯u3
u4= ¯u4
ここで，複雑化を避けるために，入力の数を出力の数 と同じに設定する．また，UAVの絶対位置とヨー角を 制御したいので，出力関数は次のように選ばれる．
y=h(¯x) = [X Y Z ψ]T
得られた拡張システムは下の形式の方程式で述べる．
˙¯
x=f(¯x) +
∑4 i=1
gi(¯x)¯ui (5)
ここで，
¯
x= [X X Y˙ Y Z˙ Z ϕ˙ ϕ θ˙ θ ψ˙ ψ u˙ 1 u˙1]T また，
f(¯x) = [ ˙X W1u1 Y W˙ 2u1 Z W˙ 3u1−g ϕ˙ 0 ˙θ 0 ˙ψ 0 ˙u1 0]T g1(¯x) = [0 0 0 0 0 0 0 0 0 0 0 0 0 1]T
g2(¯x) = [0 0 0 0 0 0 0 1 0 0 0 0 0 0]T g3(¯x) = [0 0 0 0 0 0 0 0 0 1 0 0 0 0]T g4(¯x) = [0 0 0 0 0 0 0 0 0 0 0 1 0 0]T
W1= sinϕsinψ+ cosϕcosψsinθ W2= cosϕsinθsinψ−cosψsinϕ W3= cosθcosϕ
非干渉化された線形システムに変換するためには，次 の状態フィードバック制御則が求められる．
u¯ =α(¯x) +β(¯x)v (6) α(¯x)は4×1のベクトル，β(¯x)は4×4の行列である．
補助制御入力vは後で定義する．
非線形システム(5)に対して，入力が現れるまでの 微分階数ri(i= 1,2,3,4)は次のように表される．
r1=r2=r3= 4，r4= 2 以上より，次の式が得られる．
[X(4) Y(4) Z(4) ψ]¨T=b(¯x) + ∆(¯x) ¯u ここで，
∆(¯x) =
Lg1Lrf1−1h1(¯x) ... Lg4Lrf1−1h1(¯x)
... ... ...
Lg1Lrf4−1h4(¯x) ... Lg4Lrf4−1h4(¯x)
b(¯x) =
Lrf1h1(¯x) ...
Lrf4h4(¯x)
行列∆(¯x)はu1 ̸= 0,−π/2 < ϕ < π/2,−π/2 < θ <
π/2を満たす任意の動作点で非特異である．この場合，
α(¯x) =−∆−1(¯x)b(¯x) β(¯x) = ∆−1(¯x)
を持つ状態フィードバック制御則は干渉問題を解決す ることができる．以上より，次の式が得られる．
[X(4) Y(4) Z(4) ψ]¨T=v
さらに，拡張されたシステム(5)の次元は次の条件 r1+r2+r3+r4= 14
を満たすため，システムは動的フィードバックにより，
適切な座標変換で線形可制御なシステムに変換するこ とができる．座標変換は次のように与えられる．
x1=h1(¯x) =X， x8=L3fh2(¯x) =Y(3) x2=Lfh1(¯x) = ˙X， x9=h3(¯x) =Z x3=L2fh1(¯x) = ¨X， x10=Lfh3(¯x) = ˙Z x4=L3fh1(¯x) =X(3)，x11=L2fh3(¯x) = ¨Z x5=h2(¯x) =Y， x12=L3fh3(¯x) =Z(3) x6=Lfh2(¯x) = ˙Y， x13=h4(¯x) =ψ x7=L2fh2(¯x) = ¨Y， x14=Lfh4(¯x) = ˙ψ 次に，拡張された状態ベクトルxと補助制御入力vは 次のように定義する．
x= [X X˙ X X¨ (3) Y Y˙ Y Y¨ (3) Z Z˙ Z Z¨ (3) ψ ψ]˙T v= [v1 v2 v3 v4]T
よって，次の線形時不変状態方程式が得られる．
˙
x=Ax+Bv (7)
ここで，
A=
A1 O4×4 O4×4 O4×2 O4×4 A1 O4×4 O4×2 O4×4 O4×4 A1 O4×2
O2×4 O2×4 O2×4 A2
，B=
B1
B2
B3
B4
A1=
0 1 0 0
0 0 1 0
0 0 0 1
0 0 0 0
，A2= [0 1
0 0 ]
，B1=
0 0 0 0
0 0 0 0
0 0 0 0
1 0 0 0
B2=
0 0 0 0
0 0 0 0
0 0 0 0
0 1 0 0
，B3=
0 0 0 0
0 0 0 0
0 0 0 0
0 0 1 0
，B4=
[0 0 0 0
0 0 0 1
]
(A, B)は可制御であることに注意されたい．
3.2 縮小ホライゾンMPC
本稿では，いわゆる縮小ホライゾンMPC(shrinking- horizonMPC(SH-MPC))法(see e.g. [9])を採用してい る. SH-MPCはMPCの一種であり，終端時刻Lが固 定され，各更新時刻で予測ホライゾンの長さを減らし，
特定の有限時間最適制御問題を解く手法である．
有限時間最適制御を定式化するためには，評価関数 を導入する．更新間隔とステップはそれぞれ∆tとkで 表される．制御間隔はN 個に分け，L=N∆tで表さ れる．次に，更新時刻tkは以下のようになる．
tk=k∆t (k= 0,1,2, . . . , N −1)
このモデル予測着陸制御に対して，評価関数を次のよ うに定める．
Jk =1 2
∫ L tk
(q(|V ×F|2+v(t)TΩv(t))dt
ここで，q >0と正定値行列Ω∈R4×4は重みである．
右辺第1項は，UAVの進入降下角をγに近づけること を目的とする．第2項は，制御入力の急激な変化を避 けることを目的とする．
以上より，移動GVへのUAVのモデル予測着陸制御 は，式(7)と終端制約(4)の下で，各更新時刻tkで評 価関数Jkを最小化する最適制御問題として定式化され る．表記を簡単にするために，次のように表される．
xk[τ] :=x(tk+τ) また，ほかの変数も同じ表記で表される．
評価関数は次のように書き換えることができる．
Jk =1 2
∫ L tk
(x(t)TQx(t)
+ 2s(t)x(t) +R(t) +v(t)TΩv(t))dt
=1 2
∫ L−tk 0
(xk[τ]TQxk[τ]
+ 2sk[τ]xk[τ] +Rk[τ] +vk[τ]TΩvk[τ])dτ
ここで，
Q=
q(δ2V2
T+V2
Y) O1×3 −qVX VY O1×3 qδVX VT O1×5 O3×1 O3×3 O3×1 O3×3 O3×1 O3×5
−qVX VY O1×3 q(δ2V2 T+V2
X) O1×3 qδVY VT O1×5 O3×1 O3×3 O3×1 O3×3 O3×1 O3×5 qδVX VT O1×3 qδVY VT O1×3 qV2
T O1×5 O5×1 O5×3 O5×1 O5×3 O5×1 O5×5
sk[τ] =
q(−δ2V2
TXT−V2
YXT+VX VY YT−δVX VT h) 0
0 0 q(−δ2V2
TYT−V2
XYT+VX VY XT−δVY VT h) 0
0 0 q(−V2
Th)−δVY VT YT−δVX VT XT) 0
0 0 0 0
T
Rk[τ] =V2 T h2+δ2V2
T(X2 T+Y2
T) + 2δVT h(VY YT+VX XT) +V2
X Y2 T+V2
Y X2
T−2VX VY XT YT
3.3 MPC設計
評価関数Jkと終端制約(4)はそれぞれ以下のように 表される．
Jk = 1 2
∫ L−tk 0
(xk[τ]TQxk[τ]
+ 2sk[τ]xk[τ] +Rk[τ] +vk[τ]TΩvk[τ])dτ， (8)
xk[L−tk] =xL:=
[XT k[L−tk] VX 0 0 YT k[L−tk] VY 0 0 h 0 0 0 κ 0]T (9)
ここで，
Q=
q(δ2V2
T+V2
Y) O1×3 −qVX VY O1×3 qδVX VT O1×5 O3×1 O3×3 O3×1 O3×3 O3×1 O3×5
−qVX VY O1×3 q(δ2V2 T+V2
X) O1×3 qδVY VT O1×5 O3×1 O3×3 O3×1 O3×3 O3×1 O3×5 qδVX VT O1×3 qδVY VT O1×3 qV2
T O1×5 O5×1 O5×3 O5×1 O5×3 O5×1 O5×5
sk[τ] =
q(−δ2V2
TXT−V2
YXT+VX VY YT−δVX VT h) 0
0 0 q(−δ2V2
TYT−V2
XYT+VX VY XT−δVY VT h) 0
0 0 q(−V2
Th)−δVY VT YT−δVX VT XT) 0
0 0 0 0
T
Rk[τ] =V2 T h2+δ2V2
T(X2 T+Y2
T) + 2δVT h(VY YT+VX XT) +V2
X Y2 T+V2
Y X2
T−2VX VY XT YT
従って，モデル予測着陸制御は，線形システム(7)に対 し，終端制約(9)の下で，各更新時刻tkでJkを最小化 する最適制御問題とみなすことができる．
この最適制御問題を解くために，ハミルトニアンH を定義する．
H(x,λ,v, tk+τ) := 1
2(xk[τ]TQxk[τ] + 2sk[τ]xk[τ]
+Rk[τ] +vk[τ]TΩvk[τ]) +λTk[τ](Axk[τ] +Bvk[τ]) ここで，λk[τ]∈R14は随伴ベクトルである．
次の定理は最小原理としてよく知られている(see e.g.
[10])．
Theorem 評価関数JKを最小化する最適制御入力 vk[τ](0≤τ≤L−tk)が存在するとし，対応する軌道 をxk[τ]とする．このとき，ベクトル値関数λk[τ]が存 在して，以下が成り立つ．
˙
xk[τ] =Axk[τ] +Bvk[τ]
xk[0] =x(tk) (10)
xk[L−tk] =xL
λ˙k[τ] =− (∂Ho
∂x )T
(x,λ, tk+τ) Ho(x,λ, tk+τ) = min
v H(x,λ,v, tk+τ) (11) ここで，制御入力vは制約をかかていない，H も滑ら かであるため，式(11)の最小値は次式で表される．
∂H
∂v (x,λ,v, tk+τ) = 0 (12) また，最小原理に基づいて最適制御問題を解くことは，
境界条件xk[0] =x(tk)とxk[L−tk] =xLの下で，以下 の連立方程式を解くことと同じである．
[
˙ xk[τ]
λ˙k[τ]
]
= [
A −BΩ−1BT
−Q −AT ] [
xk[τ]
λk[τ]
] +
[ 0
−sk[τ]T ]
(13) 式(12)により，最適補助制御vk[τ]は次のように与え られる．
vk[τ] =−Ω−1BTλk[τ] (14) ここで，λk[0]は未知であり，xk[0]とxk[L−tk]は既 知であることに注意されたい．最適軌道を取得するに は，λk[0]を見つける必要がある．
式(13)から，次のようになる．
[ xk[τ]
λk[τ] ]
= [
E11[τ] E12[τ]
E21[τ] E22[τ]
] [ xk[0]
λk[0]
] +
[ l1k[τ]
l2k[τ]
]
(15) ここで，
E[τ] = [
E11[τ] E12[τ]
E21[τ] E22[τ]
]
= exp(
[
A −BΩ−1BT
−Q −AT ]
τ)
lk[τ] = [
l1k[τ]
l2k[τ]
]
=
∫ τ 0
E[τ−ζ]
[ 0
−sk[ζ]T ]
dζ
次に，τ =L−tkで，以下の式が得られる．
xL =E11[L−tk]xk[0] +E12[L−tk]λk[0] +l1k[L−tk] E12[L−tk]は正則ならば，λk[0]は次のように求めら れる．
λk[0] =E12[L−tk]−1(xL−E11[L−tk]xk[0]−l1k[L−tk]) (16)
上記の議論から，着陸制御に対する縮小ホライゾンMPC のアルゴリズムは以下のとおりである．
Step 0: E12[L−tk]がk= 0,1, . . . , N−1に対して正則になる ように，L,∆t及びNを選ぶ．
Step 1: 式(16)により，初期値λk[0]を計算する．
Step 2: k= 0を設定する．
Step 3: ステップ1のλk[0]と式(10)のxk[0]は既知であるため，式 (15)により，最適軌道(xk[τ],λk[τ])(0≤τ≤L−tk)を取得 することができる．次に，最適補助入力vk[τ](0≤τ≤L−tk) は式(14)により与えられる．
Step 4: 得られたvk[τ]をv(tk+τ)に置き換える．次に，式(6) を解くことにより，最適制御入力が得られる．
Step 5: 最適制御入力の最初の間隔の値のみを適用する．すなわち，
ui(tk+τ) =uik[τ](i= 1,2,3,4)，τ∈[0,∆t)
Step 6: k=N−1の場合，終了する(着陸が完成).それ以外の場
合はkをk:=k+ 1に更新し，ステップ2に戻る．
提案されたMPCシステムの有限時間最適制御のブロッ ク線図はFig. 4に示す．
Dynamics of UAV ഥ
𝒖 = 𝛼 ഥ𝒙 + 𝛽 ഥ𝒙 𝒗
∫ ∫
ത 𝑢1 ത 𝑢2 ത 𝑢3 ത 𝑢4
𝑢2
𝑢1 𝑢3 𝑢4 ሶ𝑢1
𝝀𝑘[𝜏] ሶ𝐱𝑘[𝜏] = 𝐴𝐱𝑘[𝜏] − 𝐵Ω−1𝐵T𝝀𝑘[𝜏]
ሶ𝝀𝑘[𝜏] = −𝑄𝐱𝑘[𝜏] − 𝐴T𝝀𝑘[𝜏] − 𝒔𝑘[𝜏]T 𝐱𝑘[0] = 𝐱(𝑡𝑘), 𝝀𝑘[0] = 𝜆(𝑡𝑘)
𝒗𝑘[𝜏]
Linearized Dynamics
Calculation of the optimal control
−Ω−1𝐵T
𝐱 𝑡𝑘
𝒙
Fig. 4: Finite-horizon Optimal Control in the MPC-based Landing Control System
4 シミュレーション
本節では，MATLAB/Simulinkを用い，提案された 制御手法の有効性を数値シミュレーションによる検証 する．
本稿では，大きさが一定の突風を風外乱として扱う．
ここでは，t= 4から6[s]まで，一定の突風FW がUAV のX方向に発生すると仮定する．
シミュレーションでは，GVはX軸に対してπ/4の 方向で移動すると仮定する．また，t≥6[s]のとき，GV
の速度が0.3から0.5[m/s]に変化すると仮定する．シ
ミュレーションでのパラメータ設定はTable 2に示す．
Table 2:シミュレーション設定
ピッチ角の初期値 θ(0) = 0.04[rad]
更新間隔 ∆t= 0.5[s]
入力の初期値 [u1(0) ˙u1(0)]T= [g/cos(0.04) 0]T
UAVの質量 m= 0.5[kg]
慣性モーメント Ix, Iy, Iz= 0.035[kg·m2] GVの初期速度 VT= 0.3[m/s]
終端時刻 L= 15[s]
GVの初期位置 [XT(0), YT(0)] = [3 2]
GVの高さ h= 0[m]
座標変換されたシステムの初期状態 x(0) = [0 0.2g∗tan(0.04) 0 0 0 0 0 4 0 0 0 0 0]T 座標変換されたシステムの終端状態 x(L) = [XT(L)VX0 0YT(L)VY 0 0h0 0 0π/4 0]T
所望の進入降下角 γ=π/4[rad]
重み Ω=diag([0.1 0.1 0.1 0.1]), q= 1
(a)FW = 0[N]
(b)FW = 0.3[N]
(c) FW = 0.5[N]
Fig. 5: Simulation Results (Relative Trajectories of UAV to GV)
0 5 10 15 -3
-2 -1 0
(a)X(t)−XT(t)
0 5 10 15
-2 -1 0
(b)Y(t)−YT(t)
0 5 10 15
-1 0 1 2 3 4
(c)Z(t)−h
Fig. 6: Simulation Results (Relative Positions(FW = 0.5[N]))
Fig. 5はUAVとGVの相対位置の軌道を示している．
すなわち，(X(t)−XT(t),Y(t)−YT(t),Z(t)−h)．原 点に戻るとは，UAVがGVに着陸したことを意味する．
UAVが青い線上にある場合は，UAVが所望の進入降下 角を保てることを意味する．Fig. 6は風外乱が0.5[N]
の場合の相対位置の応答を示している．Figs. 5-6から，
0.5[N](VW ind ≃ 3.16 [m/s])までの風外乱が存在する 場合でも，UAVは移動GVにうまく着陸することがわ かる．ただし，進入降下角は外乱によりわずかに変動 する．
0 5 10 15
4.8 4.9 5 5.1
(a) Thrustu
0 5 10 15
-4 -2 0 2 4 6
10-3
(b) Torqueτϕ
0 5 10 15
-0.01 0 0.01
(c) Torqueτθ
Fig. 7: Simulation Results (Control Inputs)
Fig. 7は制御入力の応答を示している．ホバリング
状態(u= 0.5g ≃4.9 [N],τϕ =τθ = 0 [N·m])と比べ ると，制御入力の変動は比較的小さいため，許容範囲 以内である．
回転角の応答はFig. 8に示している．図より，回転 角は評価関数Jkに含まれていないが，FW = 0.5[N]の 場合でも回転角の変動は許容範囲以内であることがわ かる．
0 5 10 15 -0.04
-0.02 0 0.02 0.04 0.06
(a) Roll Angleϕ
0 5 10 15
-0.1 -0.05 0 0.05
(b) Pitch Angleθ
Fig. 8: Simulation Results (Rotation Angles)
上記の数値シミュレーション結果からまとめると，提 案されたMPC法が，風外乱が存在する場合でも，所望 の進入降下角で移動GVへのUAVの着陸制御に有効で あることを確認した．
注意:風外乱のもう一つの典型的なタイプは持続的な 外乱である．そのような外乱は，UAVが終端時刻に目 標着陸地点から逸脱する原因となる可能性がある．持 続的な外乱に対処するためには，UAVを所定の小さい 領域に着陸させるための領域着陸制御システムを設計 する必要がある．この問題は本研究の範囲を超えてい るため，将来の作業として扱う．
5 おわりに
本稿では，UAVの着陸制御問題に対し，部分フィー ドバック線形化を用いた縮小ホライゾンモデル予測制 御法を提案した．提案された手法の利点としては，風 外乱が存在する場合でも，UAVが移動GVの上に着陸 することができる．また，部分フィードバック線形化 により，近似誤差なしで計算コストを低減することが できる．さらに，提案された手法を三次元運動の場合 に拡張した．シミュレーション結果から，提案された 手法の有効性を検証した．
今後の課題として，実験を通じて提案されたMPC法 の実用性を検証する．また，線形モデルの高次導関数 を直接測定することができないため，MPC法を状態オ ブザーバと組み合わせる必要がある．
参考文献
[1] O. J. G. Villarreal, J. A. Rossiter, and H. Shin:
Laguerre-based adaptive MPC for attitude stabiliza- tion of quad-rotor.UKACC 12th Int. Conf. on Control (CONTROL). 2018.
[2] R. Ritz, M. Hehn, S. Lupashin, and R. D’Andrea:
Quadrocopter performance benchmarking using op- timal control. Proc. of 2011 Int. Conf. Intelligent Robots and Systems. 2011.
[3] Y. Zhou, K. Takaba: Optimal landing control of an un- manned aerial vehicle via partial feedback lineariza- tion. Proc. of 2019 Int. Conf. on Advanced Mecha- tronic Systems (ICAMechS), 2019, pp. 218–223.
[4] A. Zanelli, G. Horn, G. Frison, and M. Diehl: Nonlin- ear model predictive control of a human-sized quadro- tor.Proc. of the European Control Conference (ECC).
2018.
[5] K. Benkhoud and S. Bouallegue: Model predictive control design for a convertible quad tilt-wing UAV.
Proc. of 2016 4th Int. Conf. on Control Engineering Information Technology (CEIT-2016) Tunisia. 2016.
[6] L. R. G. Carrillo, A. E. D. Lopez, R. Lozano and C.
Pegard: Quad Rotorcraft Control Vision-Based Hov- ering and Navigation. Springer-Verlag. 2013.
[7] A. Isidori: Nonlinear Control Systems. Springer- Verlag. 1989.
[8] V. Mistler, A. Benallegue, N. K. M’Sirdi: Exact Lin- earization and Non-interacting Control of a 4 Rotors Helicopter via Dynamic Feedback, Proc. 10th IEEE Int. Workshop on Robot-Human Interactive Commu- nication, 2001.
[9] S. S. Farahani, R. Majumdar, V. S. Prabhu, and S.
Soudjani: Shrinking horizon model predictive control with signal temporal logic constraints under stochastic disturbances.IEEE Trans. on Automatic Control, vol.
64, no. 8, 2019, pp. 3324–3331.
[10] D. Liberzon: Calculus of Variations and Optimal Control Theory: A Concise Introduction. Princeton Univ. Press. 2011.
