計測自動制御学会論文集

特集 「移動体の誘導制御」

Vol. 10, No. 18, 153/159 (2011)

繰り返し学習制御を用いた航空機線形モデルにおける空力微係数の同定
藤森篤*・大原伸介*

Identification of Aerodynamic Derivatives in Aircraft Linear Model Using Iterative Learning
ControlAtsushi FUJIMORI *, Shinsuke OH-HARA*
This paper presents a system identification technique for continuous-time state-space system using the iterative learning
control. The transfer function parameters are regarded as functions with respect to the state-space parameters which will beidentified. The relationship between the state-space parameters and the response error is explicitly derived. The proposed
technique is applied to identification of aerodynamic derivatives of an aircraft linear model.
Key Words: system identification, iterative learning control, aerodynamic derivatives

１．はじめに線 形 時 不 変 モ デ ル の シ ス テ ム 同 定 法 は ， Ljung1) や時間系のシステム同定法では，入力と出力信号の微分値を必要とする 9)．特に出力信号の微分値は観測ノイズに敏感であSöderström ら 2)の教科書に纏められているように，最小二乗り，ノイズ強度が大きいときは正しくパラメータの推定がで(LS) 法をベースとしたさまざまな方法が体系化されてきていきないことがある．一方，繰り返し学習制御を用いた同定法

る．LS 法によるシステム同定法の基本的な考え方は，実シスでは，基準信号の微分値を用いるので，同定パラメータに対テムに内在するすべての動特性を励起させるように入力信号する観測ノイズの直接の影響を避けている．また，この同定を印加し，そのときの出力信号を入手する．この同じ入力信法はシステムの動特性の不十分な励起に対してロバストであ号を，同定したい線形モデルにも入力した際に得られるであることや，不安定システムの同定にも適用可能であるというろう予測出力との差 (予測誤差) の二乗値を最小とするように利点を有している．杉江らは，文献 5)～8) で伝達関数モデル線形モデルのパラメータを決定しようとする方法である．その同定手順を示した．しかしながら，状態空間(SS) モデルのこで扱うモデルは離散時間パルス伝達関数を基本としている． パラメータの同定には対応していない．
上記に示した一連の信号の印加とサンプリングによって同定

本研究は，SS パラメータの同定にも適用可能な繰り返し学計算が行われるため，これらが妥当になされているならば，習制御を提案することが目的である．基本的な考え方は，

すなわち，すべての動特性が十分励起されてノイズの少ない「TF パラメータは SS パラメータの関数である」とみなして，出力信号を入手すれば，LS 法により適切なシステム同定が可能である．
最近，繰り返し学習制御 3), 4)を用いた同定法が杉江らによって提案されている 5)～8)．この方法では，信号の印加とサンプリングを繰り返して実行し，同定したいパラメータを更新させて推定値を得ようとする方法であり，連続時間の伝達関誤差ベクトルを SS パラメータの線形表現で表す．これを SSパラメータの更新則に用いることにより，繰り返し学習制御しながら SS パラメータを同定する．提案法の応用として航空機の線形モデルにおける空力微係数の同定に本手法を適用する．なお，文献 6) などに示されている同定パラメータの収束性を改善する方策についても言及する．

数 (TF) モデルの同定法として近年注目を集めている．この方法では，信号の印加とサンプリングを繰り返すので，不慮の２．同定問題の設定

信号の印加とサンプリングにはあまり影響を受けずに同定が行うことができる．ところで，これまで提案されてきた連続本研究で扱う同定したいシステムは，次式で与えられる 1

入力 1 出力線形時不変 (LTI) システムである．
* 山梨大学大学院 医学工学総合研究部 機械システム工学
 x (t )  Ap ( ) x (t )  B p ( )u (t )専攻 甲府市武田 4-3-11
(E-mails: {afujimori, sohhara}@yamanashi.ac.jp)
 y (t )  C p ( ) x(t )  D p ( )u (t )  v(t )* Faculty of Mechanical Systems Engineering, University of
Yamanashi, Takeda 4-3-11, Kofu
(Received December 7, 2010)
TRIA018/11/1018 © 2010 SICE(1)

uf(k)(t)D(k)(p)rs(t)N(k)(p)

v (t)
+ (k)
ub (t)+ u (t)(k)r(k)(t) +K(p)-D*(p)(k)

y*(t)

Y(p)

v(k)(t)(k)+ + y (t) + e (t)
(k)

P(p)

uf*(t)

rs(t)D(k)(p)uf(k)(t) + u(k)(t)ub(k)(t)

+P(p)

(k)
- e(k)(t)+ + y (t) +-K(p)Fig. 1

An iterative learning control system for
identification

N(k)(p)x(t) は n 次元状態変数ベクトル，u(t)は制御入力，y (t)は観測

N*(p)r(k)(t) ++r*(t)

S(p)

y*(t)

-出力，v(t)は y(t)に含まれる観測雑音である．η はシステム行列に含まれる同定したいパラメータにより構成される q 次元Fig. 2 An equivalent iterative learning control system

ベクトルであり，以後，SS パラメータベクトルと呼ぶ．式(1)は，可制御かつ可観測であるとする．このとき， (1)式の３．応答誤差とパラメータ更新

伝達関数はP( p) N ( p)b p    b0
 n m
D ( p ) p  an 1 p n 1    a0m

３．１(2)と表わされる．ただし，n≧m．また，q≦m+n+1 であるとすK(p)は制御系全体の安定性を確保するために付加されたフィる．p は微分オペレータを表わし，ードバック制御器であり，その構造は未知でも構わない．

d l u (t )p u (t ) 
dt ll

ub(k)(t)はフィードバック入力，uf(k)(t)は規範信号 rs(t)を推定中(3)の伝達関数多項式 D(k)(p)を通して作られるフィードフォワーを意味する．bi , aj (i=0,1,..., m, j=0,1,..., n-1) は伝達関数

ド入力，r(k)(t)は N(k)(p)を通して作られた y(k)(t)に対する繰りの分子および分母多項式の係数であり，これらを集めて以下返し k 回目での目標値である．同定においては，y(k)(t)とのベクトルを定義する．

r(k)(t)の差から応答誤差 e(k)(t)を構成し，これを小さくするよ a  [a0  an 1 ]T   n , b  [b0  bm ]T   m 1
 a1  [Ta1]  Tn 1,   [ ] TbT Ta

同定のための繰り返し学習制御系

システム同定のための繰り返し学習制御系を Fig. 1 に示す．

m  n 1

うにパラメータを更新しようとする．もし，制御対象が安定でフィードバック制御器を特に必要としない場合は，フィー(4)ドバック入力 ub(k)(t)は不要となる．

以後，θ を TF パラメータベクトルと呼ぶ．繰り返し学習制

Fig. 1 の各信号は，以下のように表現できる．

御による同定では，システムの出力があらかじめ設定した目標値応答となるように制御を実行しながらパラメータの同定を実施する．この目標値を作り出す規範信号を rs(t)として与える．この rs(t)は滑らかで最高 n 回微分可能であるとし，以下のベクトルを定義する．
l (t )   rs (t )Tdrs (t )d rs (t ) l 1   (l  n)
dt
dt l l(5)

たとき，サンプリング時間 Ts 間隔で N+1 組の信号(t  0, Ts , , NTs )

(6)

(8)

r ( k ) (t )  N ( k ) ( p )rs (t )  mT (t )b( k )(9)u f (t )  D* ( p) rs (t )  nT (t ) a*1(10)r * (t )  N * ( p )rs (t )  mT (t ) b*(11)e( k ) (t )  y ( k ) (t )  r ( k ) (t )定することである．

 Y ( p )u (fk ) (t )  S ( p)r ( k ) (t )  S ( p )v ( k ) (t )本研究では，制御の実行時のデータサンプリングと η の更新を交互に繰り返し行いながら η を同定しようとする方法で(12)と表される．ただし，あるので，繰り返し k 回目での信号，推定パラメータやそれを用いた伝達関数多項式に( )u (fk ) (t )  D ( k ) ( p)rs (t )  nT (t ) a(1k )と表す．繰り返し k 回目での応答誤差はによって (1) 式での SS パラメータベクトル η を精度よく同(k)(7)

TF パラメータが真値であるとき，本研究の目的は，繰り返し学習制御を有限時間内で実行し{n (t ), y (t )}ub( k ) (t )  K ( p)(r ( k ) (t )  y ( k ) (t ))*，また真値を( ) なる上付き添え字を付けて表すものとする．

154

S ( p )  (1  P ( p ) K ( p )) 1(13)Y ( p )  (1  P ( p ) K ( p )) 1 P ( p )(14)

v(k)(t)=0 のとき，e(k)(t)=0 は (12) 式で添え字(k)を * に変更し0
 g S (0)g (T )g S (0)GS   S s g S ( NTs ) g S (( N  1)Ts )0 
 0  g S (0) (21)0  0 
 gY (0) (22)繰り返し学習制御系を Fig. 2 に示す．Fig. 1 に比べて，新た0 gY (0)g
T()gY (0)
GY   Y s gY ( NTs ) gY (( N  1)Ts )に破線の信号が付加されている．しかしながら，これらは式

このとき， (16a) 式はサンプリングデータによるベクトルでたものに対応するから，
y* (t )  Y ( p )u *f (t )
 S ( p)N * ( p) *D ( p )rs (t )  S ( p )r * (t )D* ( p )(15)が成立する．このことを踏まえて，Fig. 1 の制御系と等価な(15)により応答誤差 e(k)(t)において相殺されてしまうので，e( k )  GY (u (fk )  u*f )  GS (r ( k )  r* )  GS v ( k )Fig. 2 は Fig. 1 と等価なシステムとなっている．これらを含めた応答誤差 e(k)(t)は次式で表される．

 Y ( p )(u (t )  u (t ))*f S ( p )(r ( k ) (t )  r * (t ))  S ( p)v ( k ) (t )
 Y ( p )nT (t )( a(1k )   a*1 ) S ( p )mT (t )( b( k )  b* )  S ( p )v ( k ) (t )(16a)u*f   n a*1 ,r ( k )   mb( k ) ,r*   mb*(24)e( k )  GY  n ( a(1k )   a*1 )  GS  m ( b( k )   b* )  GS v ( k )(16b)

 GY  n 1 ( a( k )   a* )  GS  m ( b( k )  b* )  GS v ( k )存関係を陽に表現したものである．すなわち，e(k)(t)を小さ  ( ( k )   * )  GS v ( k )くすることは TF パラメータを真値により近づけることを意味する．これを具体的に実現するために，同定したいパラメータ数より十分多いサンプリングデータにより応答誤差を表現し，ついでパラメータ空間への射影を導く．
サンプリングデータによるベクトル表現

(25)

ただし，

  [GS  m

GY  n 1 ]   ( N 1)( m  n 1)(26)と表現できる．

サンプリング時間 Ts で採取した N+1 個の信号を扱いたい

３．３(k)

SS パラメータと応答誤差ベクトル(25) 式は，応答誤差と TF パラメータの誤差の依存関係でため，たとえば出力 y (t)に対して
y ( k )  [ y ( k ) (0)u (fk )   n a(1k ) ,を (23) 式に代入すると，(16b) 式は，応答誤差 e(k)(t)と TF パラメータ θa1(k) , θb(k)の依

３．２(23)と近似できる．また，e( k ) (t )  y ( k ) (t )  y* (t )  (r ( k ) (t )  y* (t ))
(k )
f


y ( k ) (Ts )  y ( k ) ( NTs )]T   N 1

あるが，これを SS パラメータとの依存関係へ拡張し，応答(17)

誤差を低減することは SS パラメータの推定値が真値へ収束することをもたらすことを示す．そこで，TF パラメータベを定義する．ほかの信号に対してもサンプリングデータによるベクトル表現を同様なボールド体として表現する． l (t )クトル θ は SS パラメータベクトル η の関数であるとみなしに対して次式の行列を定義する．

たとき，θ を  周りで展開した 1 次近似は

 lT (0)   rs (0) l     T ( NT )  s  l
 r ( NT )s sd l rs (0) dt l (18)
ldrs ( NTs )d rs ( NTs ) 
dtdt ldrs (0)dt

 ( )   ( ) 

ただし，N≫l とすれば，ほとんどの場合
rank l  l  1(27)と表される．ここで，    * ,    ( k ) と置いて θ*を η(k)周りで展開した 1 次近似表現は *   ( k )   ( k ) ( *   ( k ) )(28)

ただし，(19) (k ) が満たれる．

 ( ( k ) ) T(29)これを (25) 式に用いると，S(p)，Y(p)をサンプリング時間 Ts の 0 次ホールドで離散化したとき，擬似的インパルス入力列(i  0)1
uimp (iTs )  
0 (i  1,, N ) ( )(   ) Te( k )   ( k ) ( ( k )   * )  GS v ( k )(30)となる． (30) 式の両辺に(ΛΨ(k))T を掛けて

(20)

ε ( k )  M ( k ) ( ( k )   * )  ν ( k )を S(p)，Y(p)に加えたときの出力値をそれぞれ gS(t)，gY(t)とただし，し，以下の行列を定義する．

155(31)ε ( k )  ( ( k ) )T e( k ) ,ν ( k )  (  ( k ) )T GS v ( k ) ,M ( k )  ( ( k ) )T ( ( k ) )Norm of response error

(32)

0.35

0.3

ここで( k  0)(33)

0.25||e(k)||rank( ( k ) )  qが成り立てば M(k)は正則ゆえ，ν(k)=0 のとき，0.2

ε(k )  0

 (k )   *

(34)0.15がいえる．このことから， (31) 式により ε(k) (e(k))の低減を図れば SS パラメータの同定が可能である．
３．４

0.10

SS パラメータの更新則

102030

本節では， (34) 式を実現する SS パラメータの更新則を示す．杉江らはいくつかの更新則を提案しているが，本論文で4050kFig. 3 Norm of response error (Case 1)

は，文献 6)での更新則に修正を施しパラメータの収束速度を早めた更新則を提案する．繰り返し学習制御による η(k)の更新則を次式で与える．

Estimated SS parameters ( k 1)   ( k )  H ( k )ε ( k )10(35)

0

ただし，H ( k )   ( k ) ( M ( k ) ) 1(36)−10α(k)は繰り返し回数 k に依存した可変ゲインであり，(0      1)

−30(37)ηi (k ) :   −20−40とする非減少ゲインである．一例として，α(k)を線形的に増−50加させるように  k (k )   kr( k  ke )(38)−60

η1−70

η2η3−800( k  ke )1020304050kFig. 4 Estimated SS parameters (Case 1)と与える．ただし，ke(>0)は  ( ke )   を満たす繰り返し回数である． (31), (35), (36) 式より，

 ( k 1)  (1   ( k ) ) ( k )   ( k ) ( *  ( M ( k ) ) 1 ν ( k ) )(39)

本節では，前節までに示した繰り返し学習制御による SS(k)ν =0 ならば， (37) 式より
k   : (k )   *パラメータ同定法を航空機の縦の状態空間モデルに適用して，(40)空力微係数の同定を試みる．釣り合い位置からの微小擾乱をが保証される．文献 6)では，

(k )仮定して導かれた航空機の縦の状態空間表現は，

1k 1(41)(k) x (t )  Ap ( ) x(t )  B p ( )u (t ) y (t )  C p x(t )  v(t )(k)と与えている．これを用いた η の更新は， (39) 式より ηに関する「固有値」が繰り返し回数 k が大きくなるにしたがって 1 に近づくので，η(k)の更新速度が低下する．また，カルマン・フィルタを用いた場合 6)でも同様の傾向がある．このことは，繰り返し回数の増加をもたらし，場合によっては実質的に真値へ到達できないことがある． (37) 式は逆に繰り返し回数 k が増えるにしたがって固有値を小さくさせ収束速度の低下を防いでいる．

４．航空機線形モデルの空力微係数の同定

156(42)と書ける 10)．ただし，x  [u  
1q ]T , u   ey  

1Ap  E M ,B p  E L , C p  [0 1 1 0] ,

10E0
000U0100  M  0 0 0 
0,, L  Z e 0
1  M  e   [ X

Feedforward inputs
8
uf*(t)

u (t), u (t) [deg]M  ]T(43)と設定した．これらの空力微係数の真値は，k=3k=15k=50

6Z *  [6.515 73.45 4.371]T(44)4(k)
f*fである．目標値 rs(t)は，20GR ( p ) 26( p  2)6(45)なる規範モデルに単位ステップを印加して作られるステップ応答を用いた．rs(t)の微分値は GR(p)の SS 表現から求めるこ−2とが出来るが，ここではサンプリングデータ rs(t) (t =0, Ts, …,−402468NTs)を用いた数値微分により求め， (18) 式の Φl (l=n, m) を構10t [s]成した．また， (29) 式の Ψ(k)を求める際にも数値微分を利用した．観測出力信号には，雑音信号比(NSR)が 20%となる正Fig. 5 Feedforward input responses (Case 1)

規白色雑音が混入するとした．ただし，NSR の定義は

NSR=Outputs

1|| v(t ) |||| y (t ) ||(46)である．この飛行条件では縦運動は安定系であるので，制御器 K(p)は不要である．出力信号のサンプリングはサンプリ
ング時間 Ts=0.01 [sec]，データ数 N=1000 とした．同定パラy*(t), y(k)(t) [deg]0

メータ更新則の α(k)は (38) 式で与え，その上下限値を α=0.7，−1

  0.9 , また ke=50 とした．SS パラメータの初期値は Xα 以外は負となることから，−2 (0)  1 1 1T−3y*(t)k=3
k=15
k=50−4−502

468(47)と与えた．

50 回の繰り返し学習制御を実施した際の応答誤差のノル
ム，SS パラメーの推移を Fig. 3, Fig. 4 に示す．応答誤差ノ10ルムは単調に減少し，SS パラメータは漸近的に収束していt [s]き，k≧25 でほぼ真値へ収束しているのがわかる．k=50 でのSS パラメータはFig. 6 Output responses (Case 1)
 Xu
 0M  Zu
 M u g cos 00 g sin 00Xu0Z
M

 (50)  [6.489 72.18 4.2834]T0 1 U  Zq 
M q (48)で，どの空力微係数も真値に対して 2%以内の誤差となった．
繰り返し回数 k=3, 15, 50 でのフィードフォワード入力，出力信号の応答を Fig. 5, Fig. 6 に示す．フィードフォワード入状態変数 x の成分は，飛行方向速度 u，ピッチ角 θ，迎角 α，力，出力信号においても繰り返し回数が増加するにしたがっ

ピッチ角速度 q，制御入力はエレベータ舵角 δe ，出力は経て，同定パラメータが真値であるときの応答にそれぞれ近づ路角であり，すべて釣り合い位置からの変動量である．U はいていき，k=50 に於いてはほとんど一致した応答を示して飛行速度，Θ0 は定常状態でのピッチ角である．行列 E, L, Mいる．

中の Xu, Zα, …は安定軸で記述された空力微係数で，本研究

４．２

ではこれらの指定した空力微係数によって同定するべき SS

SS パラメータベクトル η をパラメータベクトル η を構成する．対象とする航空機は，文献 10) にある可変特性研究機としてその構造パラメータの数値を用い，飛行条件は高度 4000 [m]， 飛行速度 U=100 [m/s]Z eM  e ]T(49)と設定した．これらの空力微係数の真値は，とした．
４．１

  [ X

同定ケース２ *  [6.515 4.952 3.116]T

同定ケース１

SS パラメータベクトル η を

157(50)Norm of response error

Norm of response error

10.80.90.7
0.8
0.70.5(k)||e ||||e(k)||

0.6

0.4

0.6
0.5
0.4

0.3
0.3
0.2

0.2

0.1
0102030400.1
050

1020304050

4050kkFig. 9 Norm of response error (Case 3)

Fig. 7 Norm of response error (Case 2)

Estimated SS parameters

Estimated SS parameters

10806−10−20

2−30

ηi

ηi

4−400−2−4

η1
η2−50

η

η1−60

η

η2−70

η53

η3−60102030

40−80
050

4

η61020kFig. 10 Estimated SS parameters (Case 3)

Fig. 8 Estimated SS parameters (Case 2)と設定した．これらの空力微係数の真値は，である．同定条件は α=0.9 とした以外は同定ケース 1 と同様(k)とした(よって α =0.9 と固定した)．応答誤差のノルム，SS *  [6.515 0.2202 73.45

パラメーの推移を Fig. 7, Fig. 8 に示す．応答誤差ノルムは単1.534 4.371 0.2976]T調に減少し，SS パラメータは漸近的に収束していき，k≧25(53)である．同定条件は α=0.8 とした以外は同定ケース 1 と同様でほぼ収束しているのがわかる．k=50 での SS パラメータは (50)  [6.479 5.906 3.162]T30kとした．誤差のノルム，SS パラメーの推移を Fig. 9, Fig.10(51)に示す．50 回の繰り返し学習制御を実施した際の SS パラメで，Zδe の推定値には約 20%の誤差を含んだものとなった．

ータはこの例にあるように空力微係数の中には同定し易い微係数と (50)  [7.851 0.2295 71.15

そうでないものがあり，後者の微係数を精度よく同定する対1.213 4.261 0.2342]T(54)

策が必要である．
４．３となり，各パラメータの推定誤差は 20.5, 4.2, 3.1, 20.9, 3.5,同定ケース３21.3%である．このケースでも，同定しやすい微係数とそう

SS パラメータベクトル η を  [ XZuZ

ZqM

M  ]Tでないものがある． (42) 式の伝達関数の分母子多項式は 4(52)次，分子多項式は 3 次となるので，最高 8 つのパラメータが同定可能である 1)．それ以上のパラメータを一度に同定したい場合は，多変数化する必要がある．

158

５．おわりに

Journal of Dynamics Systems, Measurement, and Control,
129-4, 825/836 (2007)

本研究では，杉江らによって開発されてきた繰り返し学習

8） T. Sugie: Identification of linear continuous-time systems

制御を用いた連続時間伝達関数 (TF) モデルのシステム同定

based iterative learning control, In Recent Advances in

法を状態空間 (SS) モデル表現のパラメータの同定への拡張

Learning and Control, Springer-Verlag, London, 205/218

について示した．そこでは，「TF パラメータは SS パラメー(2008)タの関数である」とみなして，応答誤差を SS パラメータの線形表現で表した．また，同定パラメータの更新則に改良を

9） H. Unbehauen and G. P. Rao: Continuous-time approaches

加えて推定速度の向上を図った．提案法の応用として航空機

to system identification - A survey, Automatica, 26-1, 23/35の縦運動の線形モデルにおける空力微係数の同定に本手法を(1990)
10） 加藤, 大屋, 柄沢: 航空機力学入門, 東京大学出版会

適用し，その有効性と現方法における解決するべき問題点を(1982)指摘した．

参考文献
1） L. Ljung: System Identification - Theory for the User, 2nd

[著藤森

篤者紹介]

(正会員)

1988 年名古屋大学大学院博士後期課程満了(航空工学専攻)．

edn, Prentice Hall, Upper Saddle River, NJ. (1999)

89 年静岡大学工学部エネルギー機械工学科講師，90 年同2） T. Söderström and P. Stoica: System Identification, Prentice

機械工学科助教授．2006 年山梨大学大学院医学工学総合Hall, Upper Saddle River, NJ. (1989)

研究部教授．飛行制御系設計，移動ロボットに関する研

3） Z. Bien and J. X. Xu: Iterative Learning Control: Analysis,

究に従事．工学博士．日本航空宇宙学会, 日本機械学会，

Design, Integration and Applications, Kluwer Academic

AIAA 各会員．

Press, Boston (1998)
4） Y. Chen and C. Wen: Iterative Learning Control:(Convergence, Robustness and Applications), no. 248 of大原伸介(正会員)

Lecture Notes in Control and Information Sciences,

2006 年大阪大学大学院博士後期課程修了(電子制御機械工

Springer-Verlag, London (1999)学専攻)．同年国際レスキューシステム研究機構研究員，

5） 酒井, 杉江: 射連続時間システム同定のための耐雑音性

2007 年電気通信大学産学官連携研究員 , 同特任助教．

を有する反復学習制御, 計測自動制御学会論文集, 42-5,

2009 年山梨大学大学院医学工学総合研究部助教．拘束シ543/550 (2006)

ステムの制御，移動ロボットに関する研究に従事．博士(工学)．日本機械学会，IEEE 各会員．

6） 酒井, 杉江: 射影型反復学習による連続時間閉ループ同定法, 計測自動制御学会論文集, 43-10, 877/882 (2007)

7） T.-H. Kim, X. Zheng and T. Sugie: Noise tolerant iterative
learning control and identification for continuous-timesystems with unknown bounded input disturbance, ASME,

159