計測自動制御学会東北支部 40 周年記念学術講演会
2004 年 12 月 22 日，23 日 日本大学工学部

講演番号：B03

既約分解表現に基づいた物理パラメータの新しい同定法と故障診断への応用○黒沢忠輝（八戸高専）

大日方五郎

川合忠雄（名古屋大学）佐藤勝俊（八戸高専）

A New Identification Method of A Physical Parameter Based on Left Co-prime
Factorization and Its Application to Fault Diagnosis*T. Kurosawa (Hachinohe National College of Technology), G. Obinata, T. Kawai(Nagoya University) and K. Sato (Hachinohe National College of Technology)
Abstract－ Many methods have been proposed for identifying linear systems. On the other hand, a few

methods are known for estimating unknown physical parameters, which are not necessarily enough toexpress the subject dynamics. We can use such a method for estimating the parameter value of a brittleelement, which means a kind of fault diagnosis. We have proposed the new identifying method forpresuming only one physical parameter in a system. The new method is based on Youla parameterization,and the parameterization is applied to derive an input-output relation with respect to the target parameters.
In this paper, It examines that this technique applies to fault diagnosis.
Key Words:1Identification, Left Co-prime Factorization, Linear Fractional Transformation, Fault Diagnosis

はじめに関数であり，N0(s)，D0(s)はGp0(s)の既約因子である．
すなわち，次式で表される対象とするシステムの動特性を正確に表す数式モN (s)デルを得ることは，故障診断，制御系設計，状態量
G p 0 (s) = 0(2)D0 ( s )の推定などモデルベースの多くの技術の前提として重要である．これはモデリングやシステム同定と呼2.2 同定問題とその解法ばれる問題として，非常に多くの研究の中で取り上未知パラメータを推定する方法の仮定と手順を以げられてきた[1]．これらで提案された同定法には，そ下に示す．
れを使用する様々な状況に対応できる方法が揃って仮定１ 対象システムは１入力１出力系とする．
いるが，その中で質量や電気抵抗値などの物理パラ仮定２ Gp0は安定とする．すなわち未知パラメー
メータの推定にターゲットを絞った方法は多くない． タの基準値p0をGp0が安定になるように選べる．
機械システムの故障診断や経年変化の把握では，特仮定３ パラメータpの同定にはGp0の左既約分解定の物理パラメータの変動値を知りたいし，車の積に対応した表現(Fig.1)を用いるが，U(s)，V(s)はBezout
載量のように使用状況で変化する物理パラメータを等式[5]を満足することは求めない．またU(s)，V(s)は把握したいということがある．これに対し筆者らは安定プロパであるとするが，V-1(s)が安定プロパであ
線形連続時間系の従来用いられてきた同定の枠組のることは求めない．
ひとつ[2][3]を示し，その構造を利用して指定したパラ（Step1）線形分数変換(以後LFTと呼ぶ)により，未知メータだけを同定する新しい方法を提案している[4]． パラメータベクトルpに対する偏差δpi＝pi－p0iは，ほ本研究はこれに引き続き，この手法が故障診断や異とんど全てのパラメータについてFig.2のように整理常検知に応用することを検討する．また応用する際できることは知られている[6]．未知パラメータの基準の問題点を明らかにし，その対処法を与え，例を用値p0を設定し，真値からの偏差δpi＝pi－p0iを対象シスいて故障診断の有効性を示す．
テムGpから分離する．

2

同定問題の定式化

2.1

既約分解に基づくシステム表現

N0

モデルに含まれる未知パラメータをベクトルpで表し，対象システムをGp(s)と表すこととする．物理パラメータの基準値が事前に知られており，これをp0 とする．Gp(s)とGp0(s)はともに漸近安定であり，Gp(s)は次のように表現できると仮定する．
G p (s) =N 0 ( s) + V ( s) R( s)
D0 ( s ) − U ( s ) R ( s )Gp (s)

uu1 +y1+ zy
D0 -1

RxV(1)

ここに，N0(s)，D0(s)，U(s)，V(s)は安定プロパな伝達

u2 ++ y2U

Fig. 1 Description of Gp based on co-prime factorization

33

GpyzuyMq1m1ux
δpk1xzq2m2

δp(a) Uncertain system(b) LFTG p (δ p ) = FA ( M , Σ)

Fig.3 Quarter Car Suspension Model

3-1 モデリングおよび同定計算法
対象システムを Fig.3 に示す．本システムは線形モデルとすると運動方程式は次式で表される．

(3)

ここで，M，Σは次のような行列である．
⎡ y⎤⎡u ⎤ ⎡ M 11 M 12 ⎤ ⎡u ⎤ ⎫⎥ ⎢ ⎥⎪⎢ x ⎥ = M ⎢ z ⎥ = ⎢M
⎣ ⎦⎣ ⎦ ⎣ 21 M 22 ⎦ ⎣ z ⎦ ⎪⎪0 ⎤ ⎬⎡δ p1⎪⎢⎥z = Σx = ⎢%⎥x ⎪⎢⎣ 0δ pn ⎥⎦ ⎪⎭d 2 q1dq dq= − k1 (q1 − q2 ) − µ ( 1 − 2 )dt 2dtdt2d qdq dq
m2 22 = − k2 ( q2 − u ) + k1 (q1 − q2 ) + µ ( 1 − 2 )
dt
dt
dt

m1(4)

D0 = M ⎫⎪U = M M 22 , R = −δ pi ⎬⎪V = M 12 − M 11M 21−1M 22⎭

x = [ q1 q1 q2Tq2 ] ⎫⎪⎪x = ( A0 + Lδ kH ) x + Bu ⎬⎪y = Cx⎪⎭−121N 0 = M 11M ,−121(5)この書き換えができるためには，M21が可逆である必要がある．すなわちD0が可逆であることを意味する．
（Step3） 書き換えられた Fig.1 において N0(s)，D0(s)，U(s)，V(s)が安定プロパであれば，u1，y1，u2，y2 が入出力 u，y から計算できる．
x( s ) = V ( s )u + U ( s ) y(6)z ( s ) = D0 ( s ) y − N 0 ( s )u(7)である．したがって対象システムの入出力 u，y を観測し，式(6)，(7)によって x，z を算出し，x から z への伝達特性として通常の方法により R=diag(δp)を同定することができる．
例えば未知パラメータ R の同定は次のような方法が考えられる．推定誤差の評価をT

J = ∫ ( zi (t ) − δ pi x(t )) 2 dt0のように定義すると，最小化の必要条件
TTdJ= 2δ pi ∫ xi2 (t )dt − 2∫ zi (t ) xi (t ) dt = 000d δ pi(8)(9)0i

T03i2i

ここで⎡ 0⎢ k⎢ − 10⎢ m1A0 = ⎢⎢ 0⎢ k10⎢⎣⎢ m21−µ0k10
m10 ⎤ ⎫⎪µ ⎥⎥ ⎪m1 ⎥ ⎪⎥, ⎪001 ⎥ ⎪k + k2µµ⎥ ⎪− 10− ⎥ ⎪⎪m2
m2m2 ⎦⎥ ⎬⎪k2 T
B = [0 0 0] , C = [1 0 0 0], ⎪m2⎪⎪
11 T⎪
L = [0 −0] ,⎪m1m2⎪H = [1 0 −1 0]⎪⎭m1(13)である．式(12)は(14)のように δ k をゲインとするフィードバックシステムに書き換えることができる．すなわち⎡ y ⎤ ⎡ M 11 M 12 ⎤ ⎡u ⎤⎥⎢ ⎥⎢ x ⎥ = ⎢ M⎣ ⎦ ⎣ 21 M 22 ⎦ ⎣ z ⎦⎡ C ( sI − A0 ) −1 B C ( sI − A0 ) −1 L ⎤ ⎡u ⎤=⎢−1−1 ⎥ ⎢ ⎥⎣ H ( sI − A0 ) B H ( sI − A0 ) L ⎦ ⎣ z ⎦Ti(12)x = A0 x + Bu + Lz
y = Cx, x = Hx, z = δ kx

より次式が得られる．

∫ z (t ) x (t )dtδp =∫ x (t )dt(11)［問題設定］サスペンションのばねが経年または疲労により劣化するものとする．このとき，基準値 k10と偏差δk を用いて k1＝k10+δk とすると，状態方程式は次式のように表すことができる．

（Step2） 次に，Fig.2(b)をFig.1に書き換える．本研
究では，簡単のため行列Σの対角項にδpiが1度しか現れないと仮定する．式(4)中の伝達関数を用い，Fig.1中の伝達関数は次式で与えられる．
−121k2uFig. 2 Pulling out uncertain parameterµ(10)計算例

例として自動車の 1/4 サスペンションモデル[7]を用い，着目した物理パラメータの同定計算を示すとともに故障診断への応用について検討を行う．

(15)の LFT 表現に帰着する．式(15)を用いて式(5)を計算すると D0-1 はプロパではない．そこで式(2)を満たすために，安定プロパな伝達関数 W を用いて次式のよ34

うな変換を行い，式(16)の左辺および D0-1 が安定プロパとなるよう次のような変換を行う．
D 0 U V ⎤⎦ = W [ N 0⎡⎣ N 0

D0 U V ](16)

本問題の場合，式(5)から得られる D0 は分母が 4 次，分子が 2 次となるため，ここで W を1
(Ts + 1) 2

W=(17)とした．ここで式(17)は２次形式のローパスフィルタであり，1/T が遮断周波である．これは D0 を安定プロパにする際に導入する W にフィルタリング効果を持たせる事ができる意味を持つ．この式(17)の効果についてシミュレーションによる検証を行った．対象
システムの真値は m1=350[kg] ， m2=35[kg] ， k1=2 ×104[N/m]，k2=1.9×105[N/m]，µ=320[Ns/m]とし，サスペンションのばねの基準値 k10＝1.9×104[N/m]からの偏差δk を推定する．対象システムへの入力 u は白色ノイズとし，観測する出力 y に S/N＝0.01 のノイズが含まれるものとして，100 組のデータについて推定を行った． Fig.4 に，推定結果に及ぼす設計パラメー
タ T の影響を示す．実システムの２次の固有周波数はω1＝77.5 [rad/sec]であり，この付近以下に遮断周波数 1/T を合わせることで，推定結果の 100 組の平均値 Ave.は不偏性を保っており，また 1/T が固有角周波数近辺のときに最も良い推定結果が得られる．しかしこれより大きな値を採用すると，推定値の不偏性が失われることがわかる．

4

故障診断への応用

4.2 故障診断システムの構築(Step1) m1， k1，それぞれの同定器を設計する．
(Step2)まずはじめに m1 を同定することにより自動車の使用状況を把握する．
(Step3)Step2 で得られた m1 をもとに k1 を同定する．
この監視システムを構築しシミュレーションを行った．実プラントの車重 m1 は，ある区間で 50[kg]増加した後に元に戻るものとし，ばね k1 は劣化時間の短縮のため 100[sec]間に 0.1％の劣化が生じるものとした．入力は白色ノイズとし，入力および出力は測定できるものとする．Fig.5 に計算結果を示す．車重m1，ばね k1 の同定結果は実プラントの変化に対し良好に追従しており，質量 m1 が変化してもばね k1 の推定には影響がないことがわかる．
4.3 複数のパラメータ変動に対する考察
複数パラメータが変動するとした問題設定のように，想定したパラメータ以外の真値が全てわかっているとするには実際上の無理があると思われるので，想定した以外のパラメータは小さな不確かさを有していると考え，提案した同定方法を故障診断に適用する方法を述べる．
1) 故障が想定される箇所に対応したパラメータをFig.2 で示した方法によって抜きだし，LFT 表現を求める．
2) 不確かさがあると考えられるパラメータについては次のようなパラメータ感度によりあらかじめ故障診断に大きな影響を与える事項をチェックする．3.1 節を例にすると，質量 m1 の変化に対するばね k1 の変化の影響を次の感度関数で表す．
∆m1m∂m kS = lim 1 = 1 1∆k → 0 ∆k∂k1 m11k1

システムのモデル化が完了した後，同定器の基準値を真値にすれば推定偏差は零になる．この場合において，同定器に偏差が現れた場合にはパラメータの変化すなわち故障や劣化などの異常とみることができる．そこで故障の箇所はあらかじめ想定できると仮定し，前章の 1/4 サスペンションモデルを用い，故障診断システムの構築を試みる．
4.1 問題設定
サスペンションのばね k1 は使用年数とともにゆっくりと劣化する．また車重 m1 は乗車人数や積載物など使用状況により随時変化する．自動車の使用状況を把握しながらばねの「へたり」を監視したい．

=

(18)

∂m1 ∂G k1
1 ∂G k1=∂G ∂k1 m1 ∂G ∂k1 m1∂m1500

Plant Parameter m1
m1 Identifier

m1 [kg]450
400
3503020100204020010

80100Plant Parameter k1
k1 Identifier

020000-1019990-20-300.160Time [sec]k1 [N/m]Error [%]

300

Ave.
Max.
Min.

1

101001/T [rad/sec]

19980
01000

204060Time [sec]Fig.5 Identification Result

Fig.4 Identification to k1

3580100

影響が大きくなる周波数領域を確認するため，参考文献
Fig.6 に式(18)の周波数応答を示す．また，状態(1) L. Ljung, System Identification Theory for the User
変数の単位の取り方で不確かさの影響を抑制でSecond Edition, Prentice-Hall, (1999).
きる場合には，単位の調整を行う．
3） 式(5)によって Fig.1 に対応した表現を求め，x と(2) P. M. J. Van den Hof, R. J. P. Schrama, Identificationz が入出力 u，y から計算できるようにする．こand Control -Closed Loop Issues, Vol.31, No.12,の際 N0(s)，D0(s)，U(s)，V(s)を決定する際の自pp.1751-1770, (1995).
由度である M を 2)で確認した不確かさの影響を(3) S. Dasgupta, B. D. O. Anderson, A Parameterization
考慮して決定する（複数のパラメータに対してfor the closed-loop Identification of Nonlinear Time
M の要素を独立に選択できることに注意する）．
Varying Systems, Automatica, Vol.32, No.149,
4） 式(10)に含まれる 2 つの積分を実時間でかつデpp.1349-1360, (1996).
ィジタルで行うことになるが，これについては
やはり 2)で確認した不確かさの影響を考慮する． (4) 大日方，黒沢，川合，「事前情報を活用する物理パラメータの同定方法（第一報：新しい同定方法
5 おわりにの提案）」，日本機械学会論文集 C 編，70 巻，691 号，pp714－719，(2004)．
得られた結果をまとめると以下のようになる．
(1) 着目したパラメータだけを同定する方法を提案(5) K. Zhou, J. C. Doyle, K. Glover, Robust and Optimal
した．
Control, Prentice-Hall, (1996).
(2) 提案した同定法の枠組みを成立させるためには，定義するいくつかの伝達関数を安定プロパなも(6) S. Boyd, C. Barratt, Linear Controller Design のにする必要があるが，その際に導入する任意
Limits of performance -, Prentice Hall, (1991).
性のある伝達関数がフィルタとして働き推定結(7) N.H. McClamroch, State Models of Dynamic
果に及ぼす影響を例題によって明らかにした．
Systems, Springer-Verlag New York, (1980).
(3) いくつかのパラメータに変化が起きるような状
況で本同定法を使用した場合，パラメータ推定値は他のパラメータ変化の影響を受けるが，その際の問題点を明らかにし，例を用いて故障診断の有効性を示した．

0Gain [dB]-50-100-150-200-250-300 1
10210Frequency [rad/sec]

103300Phase [deg]]

250200150
10050 110210Frequency [rad/sec]

103Fig.6 Sensitivity function36