計測自動制御学会東北支部

第 209 回研究集会(2003.6.12)

資料番号209-3

事前情報を活用する物理パラメータの同定
A Method of Identification for Physical Parameter using Advance Information○黒沢忠輝*，大日方五郎**，川合忠雄**，佐藤勝俊*○Tadateru Kurosawa*，Goro Obinata**，Tadao Kawai**，Katsutoshi Sato**八戸工業高等専門学校，**名古屋大学大学院工学研究科*Dept. Mechanical Eng., Hachinohe National College of Technology,**Graduate School of Eng., Nagoya University
キーワード：同定(Identification)，線形分数変換(Linear Fractional Transformation)，左既約分解表現(Left Co-primeFactorization)
連絡先：〒039-1192

青森県八戸市大字田面木字上野平 16-1

八戸工業高等専門学校

機械工学科

黒沢忠輝，Tel:(0178)27-7272，Fax:(0178)27-7275，E-mail:kuro-m@hachinohe-ct.ac.jp-------------------------------------------------

本論文では，線形連続時間系の従来用いられてきた

１．はじめに同定の枠組のひとつ (4)(5) を示し，その構造を利用して対象とするシステムの動特性を正確に表す数式モデ

指定したパラメータだけを同定する新しい同定方法をルを得ることは，故障診断，制御系設計，状態量の推提案する．また，この手法が故障診断などに有効であ定・予測などのモデルベースの多くの技術の前提としることを示す．て重要である．これは，モデリングやシステム同定と呼ばれる問題として，非常に多くの研究の中で取り上げられてきた (1)〜 (3)．これらで提案された同定法には，２．同定問題の定式化

それを使用する様々な状況に対応できる方法が揃って本章では，まず 2.1 節で基準対象システムが知られいるが，その中で質量や電気抵抗値などの物理パラメているという条件下で，その基準対象からの変動分をータの推定にターゲットを絞った方法は多くない．機同定する方法を記述し，そのあとの 2.2 節で本論文の械システムの故障診断や経年変化の把握では，特定の目的である指定したパラメータだけを同定する方法を物理パラメータの変動値を知りたいし，車の積載量の説明する．

ように使用状況で変化する物理パラメータを把握したいということがある．さらに従来の方法には，次に述べるような事態に対応できる方法がないように思われ

2.1

既約分解に基づくパラメータ表現

る．例えば，力学系の運動方程式中の物理パラメータ同定の対象となるシステムは同定しようとするパラを同定する問題を考える時，同定をする前に質量の値メータを係数として含む 1 入出力の伝達関数で表現では正確に知られているが粘性係数は漠然としか存在範きると仮定する．このモデルに含まれるパラメータを囲しかわからないというように，推定しようとするパベクトル p で表し，対象システムを Gp(s)と表すこととラメータごとに事前情報が異なる場合が多くある．とする．物理パラメータの基準値が事前に知られており，ころが従来，このような事前情報を陽に活用し推定のこれを p0 とする．Gp(s) と Gp0(s)はともに漸近安定で精度を上げるための一般的な方法は知られていない．

あるとし，は次のように表現できると仮定する．

−1−u
N0

Gp (s)
u1 +y1+ zy

GpyD0 -1zyux+ y2

U

N 0 ( s) + V ( s) R( s)
D0 ( s ) − U ( s ) R ( s )Fig.2 Pulling out Uncertain Parameter

(1)

2.2

ここに， N 0 ( s ) ， D0 ( s ) ， U ( s ) ， V ( s ) は安定プロパな伝達関数であり， N 0 ( s ) ， D0 ( s ) は Gp0(s)の既約因子でN0 (s)D0 ( s )(2)

同定問題とその解法

未知パラメータを推定する方法の仮定と手順を以下に示す．

ある．すなわち，

G p 0 ( s) =(b) LFT(a) Uncertain System

Fig.1 Left-coprime factorization-based description of G pG p (s) =z

δpx
u2 +u

δp xRVM

仮定１

対象システムは１入力 1 出力系とする．

仮定２

Gp0 は安定とする．すなわち未知パラメー

タの基準値 p0 を Gp0 が安定になるように選べる．

である．もし U ( s ) ， V ( s ) が次の Bezout 方程式を満た

仮定３

パラメータ p の同定には Gp0 の左既約分解に対応した表現(Fig.1)を用いるが，U ( s ) ，V ( s ) は式(1)すなら，N 0 ( s )U ( s ) + D0 ( s )V ( s ) = 1(3)を満足することは求めない．また U ( s ) ，V ( s ) は安定プロパであるとするが， V −1 ( s ) が安定プロパであること

式(1)は， R ( s ) が安定プロパであるとき，コントローラは求めない．

U ( s ) / V ( s ) によって安定化される全てのシステムを表（ Step1 ）(6)

未知パラメータベクトルの基準値 p0 は与すことが知られており ， N 0 ( s ) ， D0 ( s ) ，U ( s ) ，V ( s )えられ，仮定 2 は満足されるとする(仮定 2 は後の例でを既知として R ( s ) を同定することは，システム Gp0(s)わかるように，それほど厳しい条件ではない)．未知パがコントローラ U ( s ) / V ( s ) によって制御されている閉ラメータベクトル p に対しての偏差 δ pi = pi − p0i は，ループの状態での同定法の表現として用いられてきた．

ほとんど全てのパラメータについて Fig.2 のようにしかしこの論文では， U ( s ) ， V ( s ) が式 (3) を満たすこ

LFT を用いて整理できることは知られている(7)．対象とは仮定せず，また閉ループ系での同定問題をあつかシステム Gp から未知パラメータの基準値 p0 からの偏うものでもない．ただし対象システムの左既約分解表

差 δ pi = pi − p0i を分離する．

現は利用する．ここで，既約分解と線形分数変換につG p (δ p ) = FA ( M , Σ )いての知られている結果をまとめておく．
結果１

R ( s ) が安定であるとき，式 (1) の Gp(s) はここで， M ， Σ は次のような行列である．

U ( s ) / V ( s ) によって安定化されるすべての伝達関数の yu   M 11 M 12  u    
 x  = M  z  = M    21 M 22   z  
0  δ p1
x z = Σ x = %  0δ pn  

クラスを表す．
結果２下線形分数変換（ Lower Linear FractionalTransformation，以後，下 LFT または FA と記述する）は次式で表される．

FA ( M , R) = M 11 + M 12 R( I − M 22 R ) −1 M 21

(4)

（Step2）(8)

次に，Fig.2(b)を Fig.1 に書き換える．本論

文では，簡単のため行列 Σ の対角項に δ pi が 1 度しか現れないと仮定する(一般には δ pi が繰り返し現れ，Σ

ここで，M
M =  11
 M 21(7)M 12 
M 22 である．Fig.1 のブロック線図は，G p ( R ) = FA ( M , R )(5)のサイズは未知パラメータ数より大きくなる ) ．式 (8)中の伝達関数を用いて，Fig.1 中の伝達関数は、(6)で表すことができる．

−2−N 0 = M 11M 21−1 ,U = M 21−1M 22 
R = −δ pi

D0 = M 21−1 ,V = M 12 − M 11M 21−1M 22 ,(9)のように与えられる．この書き換えができるためには，Gp0

M 21 が可逆である必要がある．これは D0 が可逆であるk0ことを意味することに注意しよう．

ys（Step3） 書き換えられた Fig.1 において N 0 ( s ) ，D0 ( s ) ，.
y−1sU ( s ) ， V ( s ) が安定プロパであれば， u1 ， y1 ， u2 ， y2−1..
yd0が入出力 u ， y から計算できる．このとき，

x( s ) = V ( s )u + U ( s ) y(10)z ( s ) = D0 ( s ) y − N 0 ( s )u(11)Gp0 =測し，式(10)，(11)によって x ， z を算出し， x から z+u++z

Block Diagram

Fig.3である．したがって対象システムの入出力 u ， y を観への伝達特性として通常の方法により R = diag(δ pi ) をδdx-m 0−11m0 s 2 + d 0 s + k0(18)であり，安定な系となっている．

同定することができる．例えば未知パラメータ R の同次に D0 ( s ) が可逆かつ安定プロパとなり，かつ式 (2)定は次のような方法が考えられる．推定誤差の評価をを満足するような既約因子 N 0 ( s ) ，D0 ( s ) ，U ( s ) ，V ( s )TJ = ∫ ( zi (t ) − δ pi x(t )) 2 dt0(12)を求める．これは唯一には定まらないが， D0−1 が安定プロパとなるように例えばのように定義すると，最小化の必要条件
TTdJ= 2δ pi ∫ xi2 (t )dt − 2 ∫ zi (t ) xi (t ) dt = 000dδ pi1
m0 s 2 + d 0 s + k0 ,()Ds=
0(Ts + 1) 2(Ts + 1) 2sV ( s ) = 0,U (s) =(Ts + 1) 2
N 0 ( s) =(13)より次式が得られる．と選ぶことができる．ここで T は正の実数で任意であT∫ z (t ) x (t )dtδp =∫ x (t )dti0i

T0i(19)(14)

2iるが，ノイズなどの関わりで推定精度に影響を及ぼすので，適切に選択する必要がある(その選択方法は後述する)．以上の設計により，観測した入力 u と出力 y を用いて式(10)，(11)から偏差 δ d への入力 x および出力z が得られ，式(14)により δ d が一意に求められる．

３．簡単な例

続いて，観測ノイズの影響を考慮した場合について検討する．入力 u または出力 y に観測ノイズが含まれ3.1る場合， x および z にノイズが含まれ，δ p は不変推定着目した 1 つのパラメータの同定量とならない．そこで式(19)における設計パラメータ

次式で表される簡単な１自由度系
G p (s) =y1=u ms 2 + ds + k

T の効果をシミュレーションにより明らかにする．対(15)象システムは式(15)で表され，既知パラメータを m0 = 1 ，k0 = 1 とし，また推定するパラメータの真値 d = 0.2 との各パラメータのうち質量 m ，剛性 k の値はわかってする．本システムの固有角周波数は ω 0 = 1 であることいるが，減衰 d が未知としよう．

は自明である．推定値を d̂ = d0 + δ d とし，真値に対し既知パラメータ m = m0 ， k = k0 とし， d を基準値10 ％の誤差をもつ基準値 d 0 = 0.18 からの偏差 δ d を

d 0 > 0 と偏差 δ d を用いて次のように表す．
d = d0 + δ d(16)uN0

u1 ++式(16)を用いて式(15)をブロック線図で表すと Fig.3 のD0 -1zδdようになる．偏差 δ d の入力と出力を x および z とすると，入出力関係より，次式のように表せる．

y = FA (G p 0 , δ d )u z = δ dxyy1xV(17)ここで基準システム Gp0 は

u2
Fig.4

−3−++Uy2

Measurement Noiseny推定する．対象システムへの入力 u は白色ノイズとし，Fig.5 中，Ave.は推定値 100 事例の誤差の平均値を表し，出力 y には Fig.4 に示すようにノイズ n が含まれる．

Max.，Min.は 100 事例中の誤差の最大および最小値でここで n は白色ノイズとする．Fig.5 に，S/N 比をある．図より対象システムの固有角周波数 ω 0 の前後にS/N =2E{n }E{ y 2 }(20)

良好な範囲があり，±1%未満の推定結果を得ることができる．また， 1/T ＝3[rad/sec]あたりまで，ほぼ推定とし，S / N ＝0.01 となる入力 u と出力 y の 100 組のデ値の不偏性が成立していることがわかる．このようにータを作り，これらのデータについて，推定結果に及既約分解について，事前情報に基づく基準システムにぼす設計パラメータ T の影響を示す．シミュレーショ対して適切なフィルタリング要素を既約因子として与ンでの微分方程式の解法では，固有角周波数に対してえることにより，推定の向上が見込まれる．

十分に連続系としての性質が現れ，かつ計算精度が保また，式(14)を用いて δd を計算する際の時間間隔 T は他パラメータの変化による着目したパラメータ推定への影響たれるように配慮して，時間刻みを１ [msec]とした．

3.210[sec]として，１ [msec]刻みのデータを加算して積分着目した以外のパラメータの変化による影響を調べ値を近似した．横軸は式 (19) における 1/T ，すなわちるため，以下のような検証を行った．対象システムは２次のローパスフィルタ 1/(Ts + 1) 2 の遮断周波数であ式(15)とし，Fig.6 のように m ， d ， k それぞれの推定り，縦軸は真値に対する推定値の誤差を表す．また，器を設計し，同時にそれぞれの推定を行う．基準シス

u15Error [%]

105y

Plant

Ave.
Max.
Min.

0-5

Identifier
for m

Identifier
for dIdentifierfor kδmδdδk-10-15
0.01

0.1

1101001/T [rad/sec]
Fig.5

Identification to d

Fig.60.2

δm

δm0.2

0.00.0-0.2
0.2

0.0-0.2
0.20.0δkδk0.0-0.20.2δdδd-0.20.2-0.2

Real Time Identification

0102030400.0-0.2050

10Time [sec](a)203040Time [sec]Parameter Varying to m(b)Fig.7

Identification Result−4−

Parameter Varying to k50

テム Gp0 を Gp0= Gp，すなわち実システム初期値と同一で表される．ここで，とした場合の逐次推定結果を Fig.7 に示す．同定の計I 4× 4  Ο 4× 4 Ο 4×1  , B =  −1  , 
A=
−1−1
−M f K f −M f C f  M f Ff  
C = [ 0 0 0 1 Ο1× 4 ] ,
M f = diag[m1 m2 m3 m4 ],−c200 c1 + c2 −c
−c30 
c2 + c32, 
Cf = 
 0−c3c3 + c4 −c4 −00cc
44 −k200  
 k1 + k2
 −k− k30  
k 2 + k32, 
Kf = 
 0− k3
k3 + k 4 − k 4  −k40k4  
 0
Ff = [0 0 0 1]T

算は，1 組の(x, z)の離散値が追加されるごとに 2[sec]前の離散値を 1 組捨てるようにして δm ， δd ， δk を計算し（2[sec]の時間窓），各パラメータの時間変化もつかまえるようにした． Fig7 では，0.5[sec]毎の各推定値を時系列に表示している(直線補間表示)． Fig.7(a)は t＝20[sec]においてパラメータ m をステップ状に変化させた場合，Fig.7(b)は t＝20[sec] においてパラメー
タ k をステップ状に変化させた場合である．どちらの場合においても変化したパラメータに対する推定器は変化に追従した推定を行っているが，それ以外の推定器では推定値が変動することがわかる．この結果は，着目していないパラメータ変化が起きたときに推定結(22)

果が影響を受けることを示している．着目していないである．

パラメータの変動では不規則な変化が見られるので，［問題設定］未知パラメータを k3 とする． k3 は式(22)あらかじめ変化パターンが知られている場合には，そ中の Kf にあるように複数の行および列に存在する．の変化パターンと比較することによって想定している

k3 = k30 + δ k とすると，式(21)はうちどのパラメータに変化がおきたかを知ることがで Af 0 + δ Af
Pf ( s ) =  C fきると思われる．

４．多自由度系の構造パラメータ同定(23)となる．よって式(22)はI 4× 4  Ο 4× 4
Af 0 = ,−1−1
−M f K f 0 −M f C f Ο 4× 4  Ο 4× 4
δ Af = 
,−1M
δ
K
Ο
−ff
4× 4 

前章までに未知パラメータが１つの項に存在する場合について，１自由度系を例として同定法を述べた．
しかし実際には多自由度系のようにパラメータが状態方程式や伝達関数の複数項に存在する場合が多い．そ k1 + k2
 −k2
Kf0 =  0 0こでこのような対象に本手法を用いた場合について例を示す．多自由度系の例として，建築構造物モデル(8)を用い，Fig.8 に示す．文献では地震等に対する構造物− k2k2 + k30− k300000 δ k3
δKf = 0 −δ k3
00の振動制御問題を扱っているが，本研究ではアクチュエータを同定のための入力装置とし，構造パラメータのみを文献のものと同一とし，Table.1 に示す．文献より対象システム pf(s)は0− k30
k30 + k4
− k40−δ k3
δ k300
0 00 
0 
0  ,−k4  k4  (24)となる．次に式(24)中のδAf を，次式のように表す．

 A B
Pf ( s ) = C 0 x4

Bf 
0 (21)uδ Af = Lδ k3 H(25)

ここで，
m4
c4x3

k4x2

c3

k3x1

c2

k2

c1

k1

Table.1

m3

m2

m1

Fig.8 Dynamical model of 4-degree-freedom Structure−5−

Primary Structure

m1

1.62[kg]

m2

1.48[kg]

m3

1.48[kg]m4

2.20[kg]ki (i=1,…,4)

2600

[N/m]

ci (i=1,…,4)

0.08[Ns/m]

て積分値を近似した．Fig.9 に，推定結果に及ぼす設計T 
110  L = Ο1× 4 0 −m2 m3 H = [ 0 1 −1 0 Ο1× 4 ](26)

パラメータ T の影響を示す．実システムの 1 次の固有周波数は ω1 = 13.2 [rad/sec]であり，この付近以下に遮断周波数 1/T を合わせることで，推定結果の 100 組の平である．すなわち

均値 Ave.は不偏性を保っており，また 1/T が固有角周x = A0 x + Bu + Lz 
y = Cxx = Hx, z = δ kx 

波数近辺のときに最も良い推定結果が得られる．しか(27)

しこれより大きな値を採用すると，推定値の不偏性が失われることがわかる．のように δ k3 をゲインとするフィードバックシステムに書き換えることができ，
 y   M 11
 x  =  M   21

５．おわりにM 12  u M 22   z 

本研究によって得られた結果をまとめる．

(28)

 C ( sI − A0 ) B C ( sI − A0 ) L  u =−1−1   
 H ( sI − A0 ) B H ( sI − A0 ) L   z 
−1−1

(1)

着目したパラメータだけを同定するための方法を提案した．の LFT 表現に帰着する．式(28)を用いて式(9)を計算す

(2)

−10

提案した同定法の枠組みを成立させるためには，ると D はプロパではない．そこで式(2)を満たすため定義するいくつかの伝達関数を安定プロパなものに，安定プロパな伝達関数 W を用い， N 0 ， D0 ， U ，にする必要があるが，その際に導入する任意性のV が安定プロパで，かつ D 0−1 が安定プロパとなるよう次ある伝達関数がフィルタとして働き推定結果に及式のような変換を行う．

ぼす影響を例題によって明らかにした．

 N 0D 0 U V  = W [ N 0

D0 U V ]

(29)

(3)

いくつかのパラメータに変化が起きるような状
況で本同定法を使用した場合，パラメータ推定値

本問題の場合，式(9)から得られる D0 は分母が 5 次，は他のパラメータ変化の影響を受けるが，着目し分子が 8 次となるため， W をているパラメータの変動パターンが知られている

W=1(Ts + 1)3

場合にはそのパラメータに変動が起きたかを識別(30)として， D0−1 および D0 が安定かつプロパになるようにした．シミュレーションは対象システムへの入力 u は白色ノイズとし，出力 y に S / N ＝0.01 のノイズが含ますることができることがわかった．

文

献(1) L. Ljung, System Identification Theory for the User

れるものとして，100 組のデータについて推定を行っ

Second Edition, Prentice-Hall, (1999).

た．微分方程式の解法では，4 次までの固有角周波数(2) 片山,システム同定入門，朝倉書店，(1994)．

に対して十分に連続系としての性質が現れ，かつ計算(3) 足立 修一：MATLAB による制御のためのシステ

精度が保たれるように配慮し，時間刻みを 0.1[msec]

ム同定，東京電機大学出版局，(1996)．とした．また式(14)を用いて δk3 を計算する際の時間間(4) P. M. J. Van den Hof, R. J. P. Schrama, Identification

隔 T は 1[sec]として，0.1[msec]刻みのデータを加算しand Control -Closed Loop Issues, Vol.31, No.12,
pp.1751-1770, (1995).

1510Error [%](5) S. Dasgupta, B. D. O. Anderson, A Parameterization

Ave.
Max.
Min.

5for the closed-loop Identification of Nonlinear Time
Varying

Systems,Automatica,

Vol.32,

No.149,

pp.1349-1360, (1996).

0(6) K. Zhou, J. C. Doyle, K. Glover, Robust and Optimal-5

Control, Prentice-Hall, (1996).

-10(7) S. Boyd, C. Barratt, Linear Controller Design - Limits-150.01(8) 西村 , 野波 , 中田 , アクティブ動吸振器による多of performance -, Prentice Hall, (1991).
0.1

Fig.9

1

101001/T [rad/sec]

自由度構造物の極配置指定領域指定 H ∞ 制御 , 機

Identification to k3

論 C, Vol.61, No.584, pp.1351-1358, (1995).

−6−