部分空間法を用いた 2 自由度ヘリコプタのシステム同定
M2011MM059 佐々浩貴
指導教員：高見勲

1

はじめに

2 自由度ヘリコプタは, 不安定系であり, 非線形性の強いダイナミクスを持つ. また, 多入力多出力であることから, 線形化誤差, システム自身の特性変動によるモデルパラメータの変動など, モデル化誤差が生じ, シミュレーションと実験装置との出力結果に誤差が生じるため, モデリングが困難である. そこで, ヘリコプタのモデリングに対して有効であるのがシステム同定である. システム同定とは対象とする動的システムの入出力データの測定値からモデリングを行う方法であり, 対象と「同一である」ことを証明できるような「数式モデル」を作成することである[3].2 自由度ヘリコプタは 2 入力 2 出力システムであるため, 多入力多出力にも対応した部分空間法 (N4SID 法) によって数式モデルを決定する [2]. 実験装置に部分空間法を適用し, 実験装置の振る舞いに近いモデルを構築することで, どのような制御入力を与えて良いかがわかり, 最適な制御系の設計が行えると考える.
本研究では部分空間同定法の精度を検証するため, 図 1のヘリコプタの数式モデルを用いてシステム同定を行う.
推定したモデルがプラントモデルのボード線図と一致することを実証する. また, 検証で用いた入力信号と周期を実験装置にも適用し, 推定したモデルの同定検証を行う.

2

図 2 experiment and simulationて取れるのは, 実験結果とシミュレーション結果が異なっているということである. このことを疑問に思い, プラントモデルが間違っているのではないかと考えた. そこでシステム同定が必要であると考えた. 本研究で扱う 2 自由度
ヘリコプタは 2 入力 2 出力の MIMO システムであるため, 多入力多出力にも対応した部分空間法を適用することにした. 式 (1)〜(4) は 2 自由度ヘリコプタの数式モデルである.

同定対象:2 自由度ヘリコプタ

本研究では図 1 に示すヘリコプタをモデリングした数式を用いる.2 自由度ヘリコプタはピッチ方向とヨー方向に運動する機器である. この実験装置はメインモーター,テールモータへの入力電圧をそれぞれ Vm,p ,Vm,y , ピッチ角 θ, ヨー角 ψ を観測出力とする 2 入力 2 出力の MIMOシステムである.

0 0
A=
 00000010
B− jeq,p +mpheli l2cm0Bcm



0000

Kyp2jeq,y +mheli lcm

Kpy
2jeq,p +mheli lcm
Kyy
2jeq,y +mheli lcm


B=
Kpp
 jeq,p +mheli lcm
2C=

D=4

動機

図 2 は同定対象とシミュレーションを比較した図である. 縦軸がピッチ角, 横軸が時間となっている. 図 2 から見(2)]1001000000[3 (1)− jeq,y +myheli l2[

図 1 実験機の図

01000(3)](4)

部分空間同定法

1980 年代後半に提案された部分空間法は, インパルス応答のような特殊な応答データではなく, 通常の入出力データから状態空間モデルの同定を行うことが出来る. そのため, 不安定系や閉ループ系などへも適用することが可能である. また, 多入力多出力である MIMO システムに対しても同定可能である. 本研究では部分空間法 (あるいはN4SID 法) を用いた [2].

4.1

拡大可観測行列の推定4.2

N ×N 行列

多入力多出力システムの状態空間表現を式 (5),(6) のように表現する.

x(t + 1) = Ax(t) + Bu(t) + w(t)(5)y(t) = Cx(t) + Du(t) + v(t)(6)

TT −1Π⊥UU T = I − U (U U )
Yr (t) = TT −1U Π⊥
U =0
U T = U − U U (U U )





y(t)
u(t)


y(t + 1) 
 u(t + 1) 


U(t)
=r..
..


.
.


y(t + r − 1)u(t + r − 1)

Yr (t) = Or x(t) + Sr Ur (t) + V (t)

(7)

⊥⊥
Y Π⊥
U T = Or XΠU T + V ΠU T

最後の項である V を消すことを考える. そのために,s×N行列 Φ を用いる.ΦT を式 (18) に右側から掛けて, 消去すると式 (19) が得られる.

G=C
 CA 
Or = ..
.

r−1CA

5(8)

1TY Π⊥UT ΦN
Sr = 
DCB..
.

0D..
.

CAr−2 BCAr−3...
...
..
.
..
.

00..
.
CB

0
0 .. . 
D(9)本研究で用いる部分空間アルゴリズムをまとめる.

1. 入出力から式 (20) を構成する.

+ Cw(t + k − 2) + v(t + k − 1)Ĝ = W1 GW2 = U SV T ≈ U1 S1 V1T(10)u(k), y(k); k = 1, 2, ..., N + r − 1 が利用できると仮定し,
式 (11)〜(14) を用いる.
Y = [Yr (1) Yr (2) · · · Yr (N )](11)

X = [x(t) x(2) · · · x(N )](12)

U = [Ur (1) Ur (2) · · · Ur (N )](13)

V = [V (1) V (2) · · · V (N )]

(14)

(20)

(21)

最後の近似式は特異値のうちで最も重要な n 個を用いて, そのほかを 0 にすることで得られる.

3. n×n 行列 R を選び,(pr×n) 行列Ôr = W1−1 U1 R を定義する. 式 (22),(23) について解き,Â,Ĉ を得る.
Ĉ = Ôr (1 : p, 1 : n)

(22)Ôr (p + 1 : pr, 1 : n) = Ôr (1 : p(r − 1), 1 : n)Â(23)
式 (22) の右辺は行列の部分行列を指定している.

式 (7),(11)〜(14) より式 (15) が得られる.

Y = Or X + S r U + V

1TY Π⊥UT ΦN

T −1を選び, 式
2. 重み行列 W1 = I と W2 = (ΦΠ⊥UT Φ )(21) のように特異値分解する.U, V はそれぞれ正規直交行列であり,S は対角成分以外は零で, 対角成分にG の特異値を持った行列である.U1 , V1T はそれぞれrp×n, n×α 行列である. S1 は n×n 行列である.

である. また,V (t) はV (t) = CAk−2 w(t) + CAk−3 w(t + 1) + · · ·(19)部分空間法のアルゴリズム

G=(18)雑音項の除去と表すことが出来る. このとき Or ,Sr は

(17)

Π⊥
U T を右側から式 (15) に掛けることで,U とその項を消す.

4.3

とおくと,(16)を構成する. ここで式 (16) の行列は行列 U と直交な射影である. すなわち

ここで,U 項の除去

4. 線形回帰問題から B̂,D̂ を推定する.
arg min

次に, 拡大可観測行列 Or を推定することを考える. そのために,U とその項を除くとともに, 雑音項 V を取り除くことも考える.

1 ∑∥y(t) − Ĉ(qI − Â)−1 Bu(t) − Du(t)N t=1N(15)
B,D,x0−Ĉ(qI − Â)−1 x0 δ(t)∥2(24)式 (24) から B̂, D̂ を推定する [1][2].

6

安定化

7.1

このプラントモデルは積分特性を持っていることを考慮する必要がある. そこで, プラントモデルを 4 次系から2 次系に表現し直し, 外部入力信号 d と出力 θ̇, ψ̇ を用いて同定を行うことを考えた.2 次系で表現した行列は式 (25)〜(27) である.

[Ap =B00− jeq,y +myheli l2cm

同定対象のステップ応答における立ち上がり時間の間に 5〜8 サンプル点が入る間隔をサンプリング周期とする.
また, 同定対象のステップ応答が定常値の 95 ％ に達成する時間を T95 としたとき, 式 (28) を満たす T を用いる [3].

1
1T95 ≤ T ≤ T95
154]− jeq,p +mpheli l2

周期の選定(28)(25)

Bcm[Bp =Kpp2jeq,p +mheli lcm
Kyp2jeq,y +mheli lcm[]
10Cp =

Kpy2jeq,p +mheli lcm
Kyy
2jeq,y +mheli lcm

01(26)][Dp =]

0000(27)

また, 制御器を設けることで安定化させ, 入力信号を与えることにする [4]. 図 3 のようなブロック線図を用いた.
図中の Ad, Bd, Cd 行列は式 (25)〜(27) を離散化したものである.
d++w(nT)u(nT)

Bd+1/z+x(nT)

Cd++y(nT)Ad
Kd

図 3 Block diagram of control system7

図 6 step response
図 6 より同定対象のステップ応答における立ち上がり時間に約 1〜2 秒かかっている. 以上のことから T = 0.1と定めた.

8

同定モデルの検証

本研究では部分空間同定法の検証のため,式 (25)〜(27)で求めた 2 次系のプラントモデルから N4SID 法にて同定を行った. 図 7〜10 は得られた行列式と同定対象のプラントモデルとのボード線図を比較したものである. 実線が同定対象のプラントモデルであり, 点線が同定から導き出した数式モデルを表わしている [5].

入力の選定

同定入力として Random Noise(RN) および疑似白色二値信号である M 系列信号をメインモーター, テイルモー
ターにそれぞれ 0v〜1v,-1.5v〜1.5v の電圧を印加した [5].
加えて,1000 個の入出力信号を収集した. また, 前述したプラントより, モデル次数を 2 次とした.
図 7 Bode diagram of u1-θ

図 4 u1 の M 系列信号

図 8 Bode diagram of u1-ψ
図 5 u2 の M 系列信号

図 12 ﬁgure virtual plant and plant

図 9 Bode diagram of u2-θ

図 13 ﬁgure virtual plant and plant
図 10 Bode diagram of u2-ψ以上の図 7〜10 から観測ノイズを考量した場合においてもボード線図がほぼ合致していることから部分空間法によるシステム同定の精度が高いことが実証できた.

9

実験

第 7 章で用いた外部入力と周期を用いて実験装置に適用した. また, 図 11 は外部入力信号 d と出力 θ̇, ψ̇ を用いて同定を 行った推定モデルの出力結果と実験装置の出力θ̇, ψ̇ を表している. 実線が実験装置の振る舞いを示しており, 破線が N4SID 法から得られたモデルの振る舞いを示している.

図 12,13 は入力信号 d と出力 θ̇, ψ̇ を用いて 2 次系で表現し, 積分することで実験機との出力結果を近づけることが出来た. 図 11 では θ̇ と ψ̇ の一致率はそれぞれ 30%,27%である. その時点での同定精度を上げることでさらに出力結果を近づけることができるのではないかと考える.

10成果

本研究の成果を以下にまとめる.

• 2 自由度ヘリコプタへの部分空間法の適用• 入力 (M 系列信号) の選定• 周期の選定• 安定化を行ったシステム同定• ボード線図の一致

参考文献[1] L.Ljung：『System identiﬁcation-Theory for theUser』．Prentice Hall,Englewood Cliﬀs, pp340-351,
1987.
[2] 足立修一,『MATLAB による制御のための上級システム同定』：東京電機大学出版局,1999[3] 足立修一,『MATLAB による制御のためのシステム同定』：東京電機大学出版局,1999

図 11 angular velocity of θ and ψ次に推定したモデルの出力を積分し,θ, ψ を導出した. 図12,13 はそれぞれ推定したモデルと実験装置の θ, ψ を表わしている. 実線が実験装置の出力を示しており, 破線がN4SID 法から得られたモデルの出力を示している.

[4] 松葉, 牛田, 奥：『3 自由度運動を行う小型無人ヘリコ
プタの閉ループ部分空間同定によるモデリングとモデ
ルベース飛行制御』．日本機械学会論文集, 78 巻, 第792 号, 2012.
『周波数領域におけるシステム同[5] 足立修一, 室井秀夫：定の性能評価』．第 51 回自動制御連合講演会論文集,2008.