双線形特性を持つ 3-DOF Helicopter のロバスト制御
M2013SC011 長屋秋馬
指導教員：高見勲

1

はじめにヘリコプタは救助, 軍事, 医療, 報道, 輸送と多くの場面で用いられている. 最近では小型ヘリコプタの発展により, 海外大手通販サイトでは小型ヘリコプタを利用した無人宅配サービスの導入を予定しており, 今もなお活躍の場面を広げている. ヘリコプタの自動操縦に関し, 現在はパイロットの負担を軽減する目的でホバリングや旋回などの部分的な運用はされている. 本研究ではヘリコプタの離陸や着陸といった動作は考慮せず, ホバリング状態の機体が前進する動作のみの自動操縦を考える. ヘリコプタの特徴の 1 つとして人員や貨物の有無による積載量の変化があり, これに対しロバスト性を持ったコントローラの設計を考える.
本研究で用いる実験装置は前後移動, 上下移動, ピッチングのみ可動であり, 自由度が 3 である. 実際のヘリコプ
タの自由度 6 からは複雑さは減少しているもののヘリコプタ特有の非線形性は有している. ヘリコプタはその推進力に非線形性があり, 推進力はヘリコプタの仰角と揚力の積で表される. これは双線型システムと呼ばれ, これまでいくつもの研究が行われており, 動的状態フィードバックによる厳密な線形化法 [1] や, 2 段階線形化による制御系設計 [2], パラメータ同定を伴う非干渉制御 [3] などの手法が考えられてきた. 筆者も本学卒業論文において状態フィードバックによる双線形システムの安定化手法 [4]を基にヘリコプタシステムへの適用と評価関数の導入を行ってきた.
本研究の特徴としては 2 つが挙げられる. 1 つは, ヘリコプタが持つ双線形性を無視することなくモデル化及び制御系設計を行い, 理論的な安定性を保証する. もう 1 つは貨物質量の変動に対するロバスト性を付与し, 貨物を搭載することによる性能の低下を軽減する. 本論文ではモデリング, 制御系設計, シミュレーションについて記述する. モデリングでは Lagrange の運動方程式を用いて貨物の有無に対応した数式モデルを算出する. 制御系設計ではポリトープ表現を用いて貨物の有無に対しロバスト性を持ったコントローラを設計するシミュレーションは貨物のある場合とない場合の 2 つの場合を行い, 従来法と提案法を比較して質量の変化に対するロバスト性の有無を検証する.

2モデリング

実験機概略図を図 1 に示す. 概略図内 A, B, O はそれぞれヘリコプタ本体中心, カウンターウェイト, タワーである. ヘリコプタ本体にはロータが 2 つ付いており, 適切な電圧を与えることで揚力を得ることができる. 電気モータでの制御を可能とするために点 O を挟んだ点 B にカウンターウェイトを乗せて相対的なヘリコプタ本体の質量を軽減している. 本研究ではヘリコプタ本体に貨物質量 Mc

Travel axis

λCounter
WeightBθhPitchaxis

Fb
BackRotor
Mb gθwOA
FfLoadMc gεMw gρElevation axis

Mf g
FrontRotor

図 1 実験機概略図が搭載されたモデルを考える. 本モデリングにおいて貨物の揺れや貨物を吊るす紐の伸縮, 貨物によるヘリコプタ本体の重心の変動は考えないとする. 図のようにトラベル角を λ, エレベーション角を ϵ, ピッチング角を ρ とする. それぞれ機体の前後移動, 上下移動, 姿勢を示している. モ
デリングで用いるパラメータを定義する. Kf をロータの揚力定数 [N/V], Mw , Mf , Mb , Mc をそれぞれカウンターウェイト, 前ロータ, 後ろロータ, 貨物の質量 [kg], La , Lw ,Lh をそれぞれ OA 間, OB 間, 点 A と各ロータ間の長さ[m], Jλ (ϵ), Jϵ , Jρ をそれぞれ各回転方向の慣性モーメン
ト [kg · m2 ], θh を点 A から点 O への仰角 [rad], θw を点O から点 B への仰角 [rad], g を重力加速度 [m/s2 ] とする.
ここでヘリコプタの総質量 Mh を Mh = Mf + Mb + Mcとしてモデリングを簡略化する. 前後ロータの入力電圧をVf , Vb として入力を u(t) = u⋆ (t) + uc = [Vf Vb ]T とする.
ここで uc は, 貨物を持たない時 (Mc = 0) に定常状態を維持するために必要な電圧であり, 時間 0[t] より初期入力として与え続ける. 状態変数を xp (t) = [λ ϵ ρ λ̇ ϵ̇ ρ̇]T とし, Lagrange の運動方程式より以下の運動方程式が得られた. シミュレーションはこの運動方程式を用いて行う.

Jλ (ϵ)λ̈ = −Kf La (Vf + Vb ) sin ρ+2Mh L2a sin(ϵ − θh ) cos(ϵ − θh )λ̇ϵ̇+2Mw L2w sin(ϵ − θw ) cos(ϵ − θw )λ̇ϵ̇

(1)Jϵ ϵ̈ = Kf La (Vf + Vb ) cos ρ−Mh L2a sin(ϵ − θh ) cos(ϵ − θh )λ̇2−Mw L2w sin(ϵ − θw ) cos(ϵ − θw )λ̇2−Mh gLa cos(ϵ − θh )+Mw gLw cos(ϵ − θw )Jλ ρ̈ = Kf Lh (Vf − Vb )

(2)
(3)

式 (1) 右辺の第 2 項, 第 3 項はコリオリ力の影響である.
式 (2) 右辺の第 2 項, 第 3 項は図 (1) 内 λ 軸を中心に円周外方向にかかる遠心力である. 式 (2) 右辺の第 4 項, 第 5

項はヘリコプタ総質量, カウンターウェイトにかかる重力による影響である. ここで uc を次式とおくことで重力による影響を打ち消す初期電圧を与える. なお, ヘリコプタの総質量は変動要因 Mc によって一意に求まらないため,ここでは Mc = 0 として貨物を吊るさないモデルに対する初期電圧 uc を求める.
[]
1 V0uc =2 V0ggLw Mw
V0 =cos(θw )(Mf + Mb ) cos(θh ) −Kf
Kf La(4)(5)Jλ (ϵ) = Mh L2a cos2 (ϵ − θh ) + Mw L2w cos2 (ϵ − θw ) (6)ϵ は常に零近傍であると仮定し, また θh , θw も零近傍であるため, Jλ (ϵ) ∼= J˜λ = Mh L2a + Mw L2w と近似する. 以上の近似とフィードフォワード入力の導入を行い, 式 (1),(2) を以下のように簡略化する.
Jλ λ̈ = −Kf La (Vf + Vb + V0 ) sin ρ(7)Jϵ ϵ̈ = Kf La (Vf + Vb )(8)

状態変数を xp (t) = [λ ϵ ρ λ̇ ϵ̇ ρ̇]T とし, 式 (3), (7), (8)より以下の状態方程式を得られた.

pp ppp0Ap =[Ø3·30

Ø4·2
Kf La /Jϵ
Kf Lh /Jρ

−Kf LaV0 /JλØ2·6Kf La /Jϵ−Kf Lh /Jρ[]Kf La sin ρ Ø3·211
Np (ρ) = −J˜λ
Ø2·2
Bp =C p = [ I23

制御系設計

3.1拡大系

Ø2·4 ]

I3000吸引領域]
(11)

(12)
(13)

P = conv{x(1) , x(2) , . . . , x(p) }(19)xT(i) P x(i) ≤ c,(20)i = 1, . . . , p

E = {ξ ∈ Rn | ξ T P ξ ≤ c}

P ′ = {ξ ∈ Rn | aTk ξ ≤ 1, k = 1, , , q}
−1

− aTk P −1 ak ≥ 0,k = 1, . . . , qx˙e (t) = Ae xe (t) + Be u⋆ (t) + Ne (ρ)u⋆ (t) + Br r(t)(14)y(t) = Ce xe (t)
[]
[]
Ap0
BpAe =, Be =(15)−Cp 00[][]Np (ρ)0
Ne (ρ) =, Br =(16)0I(17)(22)(23)

よって式 (20), 式 (23) が成り立ち, P ′ における漸近安定が保証されるとき, 初期値が P に内包される系の状態変数は 0 < t において常に P ′ に内包されるといえる. 本研
究で設計されるコントローラはポリトープ表現を用いてx ∈ P ′ における漸近安定を保証する.
本研究では式 (24) で示される ρ の範囲において安定性を保証する.
ρ̂min ≤ ρ ≤ ρ̂maxρ̂min = −π/8ρ̂max = π/8{(21)E の外周は V (x(t)) = c の等位線であるため, 初期値がP に内包され, V̇ (x(t)) < 0(t > 0) が常に成立するとき
V (x(t)) < c (t > 0) で x(t) → 0 に収束することが言える. また, 式 (22) で定義される超平面で囲われた範囲 P ′について, 式 (23) が成立するとき E は P ′ に内包される.

c(10)(18)

リアプノフ関数を時不変な正定対称行列 P を用いてV (xe (t)) = xTe (t)P xe (t) と表すとき, 常に V (x(t)) >0, V̇ (x(t)) < 0 が言えるならば系が漸近安定である. 初期値の取りうる範囲を式 (19) で表される凸多面体 P に内包されると仮定し, 式 (20) が成立するとき, 式 (21) で定義される P は E に内包される [5]. ここで x(i) は P の頂点であり, conv(·) は x(i) を頂点とした凸包の操作である.

]

本研究では出力を目標値に追従させるために λ と ϵ に積分器を付加した. 観測出力 y(t) と目標値 r(t) の偏差を∫te(t) = r(t) − y(t) とし, 偏差の積分を w(t) = 0 e(t)dt とする. 拡大系の状態変数を xe (t) = [xTp (t) wT (t)]T とす
ると拡大系は次式のようになる.

Ce = [Cp 0]

3.2(9)y(t) = Cp x(t)[u⋆ (t) = Ke xe (t)= Kxp (t) + Gw(t)∼0
式 (1) より, 系はゆっくりと動くと仮定し, λ̇ϵ̇ ∼= 0, λ̇2 =のように近似を行う. λ 方向の慣性モーメント Jλ (ϵ) は次の式で表され, ϵ の角度によって大きさが変化する.

{ ẋ (t) = A x (t) + B u⋆ (t) + N (ρ)u⋆ (t)式 (18) のような状態フィードバック形式の積分型コントローラにより, 系が安定であることが言えるとき定常偏差なく制御量 y(t) を目標値 r(t) に追従可能であることが知られている.

(24)

式 (16) で表される Ne ρ に関し, Ne (ρ̂m in) から Ne (ρ̂m ax)は単調増加である. その為, Ne (ρ̂min ), Ne (ρ̂max ) の両端点で安定が保証されるとき, それらに内分する Ne (ρ) に対しても安定が保証される.

3.3

貨物の変化に対するロバスト H2 制御器

貨物の質量 Mc の変化に対しロバスト性を持つコントローラを設計する. 本研究では 0 ≤ M c ≤ 5[kg] の範囲でのシステムの安定をポリトープ表現を用いて保証する.
行列 Ae , Be , Ne (ρ) の中には Mc が入っているため, Mcの変動によってこれら行列は影響を受ける. M c = 0 の時の {Ae , Be , Ne (ρ)} は {Ae(1) , Be(1) , Ne(1) (ρ)}, M c = 5

の時は {Ae(2) , Be(2) , Ne(2) (ρ)} と改めて定義する. 以下のLMI 条件式のもと二分探索を用いて γ 2 の最小化を行う.
min : γ 2
subject to : X ≻ 0

評価関数を用いてシミュレーション結果を比較する. 評価関数を式 (28) で定義する.
∫ ∞J=xTe C2T C2 xe + uT D22 D2 udt(28)
0c−1 ≻ 0[He(Ae(n) X + (Be(n) + Ne(n) (ρ̂))Y )
C2 X
D2 Y]
XC2T Y T D2T≺ 0, (n = 1, 2)−I00−I

評価関数 J で評価した結果を表 1 に示す. 表 (1) より, 貨物が搭載されていない場合においてノンロバストコントローラの性能が最も良いことがわかる. これはロバストコントローラが保守的になっている為である. 以上より,ロバストコントローラが設計通り, 貨物の重量に対するロバスト安定性を持つことがわかる.

(ρ̂ = ρ̂min , ρ̂max )[]ZBrT≻0BrX
γ 2 − trace(Z) > 0

Nonrobust
RobustX − c−1 x(i) xT(i) ≻ 0, (i = 1, ..., 8)T
c−1 − aTk Xak ≥ 0, (k = 1, 2)

0kg
0.00240.00325kg0.0041

表 1 評価関数 J4

シミュレーション

重み行列 C2 , D2 を式 (25), 式 (26) として前章の LMIを解く. システムの安定を保証するゲイン Ke を Ke =Y X −1 として算出すると式 (27) となった.
C2 = (diag [ 111111111000 ] · 10−5 ) 2 (25)12

D2 = (diag [ 1 1 ] · 10−4 )[−5.0 20.2
59.4
−16.0
Ke =5.0
20.2 −59.4
16.0]
14.7
−0.8 7.0
−14.7
0.8
7.0

(26)
29.1
29.1(27)

ヘリコプタに貨物が搭載されていない場合 (Mc = 0) と5kg の貨物が搭載されている場合 (Mc = 5) の 2 つのシミュレーションを行い, 貨物の有無による安定性の違いを検証する. 前章で設計した貨物の質量 Mc に対するロバストコントローラと貨物の質量 Mc に対しロバスト性を持たないノンロバストコントローラを比較する. なお, ノンロバストコントローラは Mc = 0 の場合を想定して設計した. ヘリコプタが浮いて ϵ = 0 で静止している状態を定常状態とし, 定常状態から目標値 λ = −2π をステップ入力で与える. 実際の実験装置ではロータの逆トルクが発生するが本シミュレーションではこの影響を考えない.
貨物が搭載されていない場合 (Mc = 0) でのシミュレー
ション結果を図 2, 3, 4, 5, 6 に示す. 図 3 より, ノンロバストコントローラはロバストコントローラに比べ大きく振れている. しかしながら, 図 2, 4, 5, 6 より, ノンロバストコントローラが早く収束していることが見て取れる.
これに対し, 貨物が搭載されている場合 (Mc = 5) でのノンロバストコントローラを用いたシミュレーションでは数値が発散し, 目標値に追従することはなかった. 貨物が搭載されている場合 (Mc = 5) でのロバストコントロー
ラを用いたシミュレーション結果を図 7, 8, 9, 10, 11 に示す. 図 7, 8, 9 より, システムは安定していることがわ
かり, Mc = 0 時のシミュレーションに見られた振動も無くなっている.

5終わりに

本研究では貨物の質量に関し変動要因を持ち非線形なヘリコプタシステムに対し, 状態フィードバックでロバスト安定性を保証するコントローラの設計を行った. シミュレーションにより提案法と従来法を比較することで従来法が受ける貨物の影響を提案法で軽減できることを示した. 本研究に用いられた双線形システムは線形項に対し,双線形項が持つ影響がごくわずかであり, 双線形項を考慮することなく制御を行うことは可能だった. その為, 今後の課題としては他の双線形システムを用いた理論検証などを行い, 双線形項が持つ影響の強いシステムに対しても本研究が適応可能であることを示すことが挙げられる.

参考文献[1] 井村, 家木, 佐伯, 和田: 動的状態フィードバックによる厳密な線形化法を用いたツインローターヘリコプ
ターモデルの基礎実験日本機械学会論文集 (C 編), 66
巻, 648 号, 2630/2637, 2000.
[2] 佐伯, 和田, 井村, 坂上: 2 段階線形化法に基づくツインローターヘリコプターモデルの飛行制御系設計日
本機械学会論文集 (C 編), 67 巻, 656 号, 1038/1045,2001.
[3] 石飛, 木下, 西: 3 自由度モデルヘリコプタに対するパラメータ同定を伴う非干渉制御日本機械学会論文集(C 編), 70 巻, 699 号, 3186/3191, 2004.
[4] F. A. Mato, C. Cosentino, A. Merola: Stabilizationof Bilinear Systems via Linear State Feedback Control Circuits and Systems II: Express Briefs, IEEE
Transactions on, vol.56, issue.1 pp. 76-80, 2009.
[5] S. Boyd, L. E. Ghaoui, E. Feron, V. Balakrishnan:
Linear Matrix Inequalities in System and Control
Theory.
Society for Industrial and Applied Mathematics,
1994.00
Nonrobust controller
Robust controller
Reference

Robust controller
Reference-1-2

λ[rad]

λ[rad]-2-4-3-4-5-6-6-8-7020406080

10002040time[s]6080100time[s]図 2 Mc = 0 時 λ

図 7 Mc = 5 時 λ

0.06

0.04
Nonrobust controller
Robust controller
Reference

0.04

Robust controller
Reference
0.03

ε[rad]

ε[rad]

0.02
0.02

00.01

0-0.02
-0.01

-0.04

-0.02
020406080

100020time[s]4060

図 3 Mc = 0 時 ϵ0.04

Robust controller

Nonrobust controller
Robust controller
0.03

ρ[rad]

ρ[rad]

0.20.100.020.01-0.10-0.2-0.010204060

80

100020time[s]4060214.59f

10V [V]

215

8214

213.5
Nonrobust controller
Robust controller

Robust controller
213204060

80

1000time[s]50100time[s]図 5 Mc = 0 時 Vf

図 10 Mc = 5 時 Vf

10214.5

98b215

V [V]

11214

213.5
Nonrobust controller
Robust controller

70100

図 9 Mc = 5 時 ρ117080time[s]図 4 Mc = 0 時 ρV f[V]100

図 8 Mc = 5 時 ϵ0.3

V b [V]

80time[s]Robust controller
21320406080time[s]図 6 Mc = 0 時 Vb100

050time[s]図 11 Mc = 10 時 Vb100