計測自動制御学会東北支部第229回研究集会(2006.6.9)
資料番号229-4

既約分解表現を利用した構造物の故障診断
A Fault Diagnosis of Structure Using Co-prime Factorization○米田祐輔*，黒沢忠輝*，佐藤勝俊*，大日方五郎**○Yuusuke Maita*, Tadateru Kurosawa*, Katsutoshi Satou*, Gorou Obinata****

八戸高専,**名古屋大学

Hachinohe National College of Technology, **Nagoya University

キーワード：システム同定(System Identification)，故障診断(Fault Diagnosis)，既約分解表現(Co-prime Factorization)，線形分数変換(Linear Fractional Transformation)

連絡先：〒039-1192 青森県八戸市田面木字上野平16–1 八戸工業高等専門学校 機械工学科
黒沢忠輝，Tel.:(0178)27-7272，Fax.:(0178)27-7275，E-mail:kuro-m@hachinohe-ct.ac.jp

1. はじめにム中の指定したパラメータだけを同定する新しい開ループ同定方法を提案した2)．この手法

橋梁，高速道路などの老朽化および輸送機は従来のヘルスモニタリング技術のような膨器や大規模産業設備などの安全かつ効果的な大なセンサ計測は必要なく，構造物全体の動延命が社会的な重要課題となってきている背特性を表す一定時間長の入出力信号１組が得景から，構造物や設備機器の状態変化や健全られればよい．本研究では塔状建築構造物を性を客観的かつ経済的に監視できる方法の確対象として，加振入力に対する構造物の応答立が急務となっている．近年，構造物の健全を得る加振実験を行い，提案する同定法の推性を常時にモニターしようとする構造物のヘ定精度向上および高速化のための入力の選定

ルスモニタリング1)と呼ばれる概念が提唱され多くの研究が行われている．
しかしながら，従来あるひずみ計などの電気的なセンサや非指針について検討を行った。

破壊検査技術によるモニタリングなどは，大2. 同定問題の定式化

型構造物を常時監視するには膨大なセンサ数

2.1 左既約分解表現への帰着が必要であることや，その耐久性，信頼性お同定の対象となるシステムは同定しようとよび運営・維持管理費用などの問題点が山積するパラメータを係数として含む１入出力のしている．これに対し筆者らは線形連続時間伝達関数で表現できると仮定する．このモデ系における閉ループ構造の表現方法の一つでルに含まれるパラメータをベクトル p で表し，ある既約分解表現の構造を利用して，システ対象システムをGp(s)と表すこととする．物理-1-

2.2 同定問題とその解法パラメータの基準値が事前に知られており，これをp0とする．Gp(s)とGp0(s)はともに漸近安未知パラメータを推定する方法の仮定と手定であるとし，Gp(s)は次のように表現できる順を以下に示す．と仮定する．

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

N ( s) + V ( s) R( s)G p (s) = 0
D0 ( s ) − U ( s ) R ( s )2) Gp0(s)は安定とする．すなわち未知パラメ(1)ータの基準値p0 をGp0(s)が安定になるよ

ここでN0(s)，D0(s)，U(s)，V(s)，R(s)は安定プうに選べる．

ロパな伝達関数であり，N0(s)，D0(s)はGp0(s)3) パラメータpの同定にはGp0(s)の左既約分の既約因子である．すなわち，解に対応した表現(図１)を用いるが，U(s)，N (s)G p 0 ( s) = 0
D0 ( s )V(s)は式(3)を満足することは求めない．

(2)

またN0(s)，D0(s)は安定プロパであるとすである．もしU(s)，V(s)が次のBezout方程式

N 0 ( s )U ( s ) + D0 ( s )V ( s ) = 1

るが，V-1(s)が安定プロパであることは求めない．

(3)（Step1）未知パラメータベクトルの基準値p0を満たすなら，式(1)は，R(s)が安定プロパでは与えられ，仮定2は満足されるとする．p0あるとき，コントローラU(s)／V(s)によって安に対しての偏差δpは，ほとんど全てのパラメ定化される全てのシステムを表すことが知らータについて図2のように線形分数変換

れており3)，N0(s)，D0(s)，U(s)，V(s)を既知と(Linear Fractional Transformation，以後LFTと呼してR(s)を同定することは，システムGp0(s)がぶ)を用いて整理できることは知られている4)．

コントローラU(s)／V(s)によって制御されて対象システムGp(s)から未知パラメータの基いる閉ループの状態での同定法の表現として準値p0からの偏差δpを分離する．

用いられてきた．しかしこの論文ではU(s)，G p (δ p) = FA ( M , Σ)V(s)が式(3)を満たすことは仮定せず，また閉

ここで M ，Σは次のような行列である．

ループ系での同定問題を扱うものでもない．

⎡ y⎤⎡u ⎤ ⎡⎢ M 11 M12 ⎤⎥⎥M=⎢x⎥⎢ z ⎥ = ⎢⎢ M⎣ ⎦⎣ ⎦ ⎣⎢ 21 M 22 ⎥⎦⎥⎡0 ⎤⎥⎢δ p1⎢⎥⎥xz =Σx= ⎢⎢%⎥ただし対象システムの左既約分解表現の構造は利用する．

⎢⎢
⎣⎢

u
N0

Gp (s)
u1 +y1+ z

Gpyz+ y2⎥δ pn ⎥⎦⎥0uyMxUz

δp(a) Uncertain System(b) LFT

Fig.1 Left co-prime factorization based
Fig.2 Pulling out Uncertain Parameter

description of Gp(s)-2-u

δp xxu2 +(5)yD0 -1RV(4)（Step2）次に，図2(b)を図1に書き換える．式(5)中の伝達関数を用いて，図1の伝達関数は、

N 0 = M 11M 21−1
D0 = M
U = M 21−1M 22−121

V =M12 − M11M 21−1M
R=−δ piGp (s)
u++D0 -1⎫⎪⎪⎪⎪⎪⎪⎬⎪
⎪
22 ⎪⎪
⎪⎪⎭

x+y

N0+R(6)
UzVのように与えられる．この書き換えができる

Fig.3 Right Co-prime Factorization Based

ためには，M21が可逆である必要がある．これ

Description of Gp(s)はD0(s)が可逆であることを意味する．
（Step3）書き換えられた図1においてN0(s)，研究では左既約を主に扱う．ここで本論から

D0(s)，U(s)，V(s)が安定プロパであればx，z見た左既約と右既約の違いをまとめておく．

が入出力u，yから計算できる．このとき，図3に閉ループ右既約分解表現を示す．
ここでx( s ) = U ( s ) y + V ( s )u ⎫⎬z ( s ) = D0 ( s ) y − N 0 ( s )u ⎭対象システムGp(s)の伝達関数は左既約(図1)(7)の場合と同じく式(1)で表される．R(s)の前後である．したがって対象システムの入出力u，の入出力x，zは，左既約の構造では式(7)の形yを観測し，式(7)によってx，zを算出し，xかで得られるが，右既約の構造の場合にはUy + Vu ⎫D0V + N 0U ⎪⎪⎬D0 y − N 0u ⎪z=D0V + N 0U ⎪⎭らzへの伝達特性として通常の方法によりx=R=diag(δp)を同定することができる．例えば
未知パラメータRの同定は次のような方法が考えられる．推定誤差の評価をとなり，Bezout方程式(3)が成立すれば左既約TJ = ∫ ( z (t ) − δ px(t )) 2 dt0(8)と右既約は等価であることが知られている．のように定義すると，最小化の必要条件
TT
dJ= 2δ p ∫ x 2 (t )dt − 2 ∫ z (t ) x(t )dt = 000dp

しかし本同定法において，未知パラメータの偏差は，R ( s ) に分離されるとともに単なるゲ(9)イン(係数)であり，同定するためには前後のより次式が得られる．

入出力 x ， z がわかればよく，左既約と右既約が等価である必要はないため，式(3)を満たT

δp=2 ∫ z (t ) x(t )dt0∫T0(11)

すことは考えない．対象システムを閉ループ(10)x 2 (t )dt

左既約の構造に帰着させる場合には式(7)を，これにより未知パラメータの偏差δpが一意に右既約の場合には式(11)を用いればよい．も算出することができる．

ちろん左も右も等価となるように式(3)を満たすような因子の設計も可能ではあるだろう

2.3 右既約の扱いについてが，複雑化を招く場合もあるだろう．これらのことから本論では数式構造の簡単な式(7)対象システムの既約化は，左既約と右既約を採用している．の２通りの方法が存在するが，前述の通り本-3-式のように表される．

3. 構造パラメータ推定法⎫⎪⎪(Ts + 1) 2 ⎪−1
D0 ( s ) =⎪m0 s 2 + c0 s + k0 ⎪⎬1⎪
U (s) =⎪
(Ts + 1) 2⎪
V ( s) = 0⎪⎪⎭
R ( s ) = −δ k
N 0 ( s) =(問題設定)次式で表される簡単な１自由度系

G p ( s) =1
ms + cs + k
2(12)の各パラメータのうち質量 m ，粘性 c の値はわかっているが，剛性 k 減衰が未知とする．
(Step.1)対象システムを左既約分解表現に書き換える．この場合の必要条件は，既約因子(17)

このように未知パラメータの偏差δkはRに単N0(s)，D0-1(s)が安定であることと，D0(s)が存なるゲインとして表されることになる．した在することである．つまりD0-1(s)は可逆かつ安定プロパとなるよう注意を払う必要があり，がって対象システムの入出力u，yを観測すれば，式(7)によってR(s)の前後の入出力x，zは例えば1(Ts + 1) 2算出され，偏差δkは式(10)により一意に算出(13)することができる．

を用いると対象システムの左既約は

G p ( s ) = N 0 D0−1⎫⎪⎪⎬(Ts + 1) 2 ⎪−1D0 = 2
ms + cs + k ⎪⎭N0 =1
(Ts + 1) 2

4. 既約分解表現を利用した(14)

パラメータ推定と故障診断(15)

4.1 加振実験装置
本同定法の有効性を検証するために１自由と表され，D0-1(s)の必要条件を満足すること度系を用い加振入力に対する応答実験を行っが明らかであるとともに，必要条件を満たすた．図４に本実験に用いた実験装置の概要を調整の段階で上記のようなフィルタリング設計を施すことが可能であるという特徴を持つ．
(Step.2)

1(Ts + 1) 2

示す．

次に開ループである対象システムyを，閉ループ左既約分解表現の構造に帰着させる．既知パラメータとして質量m0，粘性c0mとし，未知パラメータである剛性kを基準値k0と偏差δkを用いて次のように表す．

k = k0 + δ kc,k(16)次に式(16)を用いて閉ループ左既約分解表現に書き換える．本問題の場合，LFTを行わなくともブロック線図上で行うことができ，次Fig.4 1 Degree of Freedom System-4-uTable.1 Specification

401.775[kg]c0.157[Ns/m]

k730.6[N/m]

1st mode

20.3[rad/s]20Error [%]m

0-20-40-60

はりは本来分布定数系であるが１次モードの10みで扱う集中質量系とみなし，天井面の質量

1520
1/T [rad/sec]

2530質量m，はりが持つ粘性cおよび剛性kを用い

Fig.5 Effect of Cut Off Frequency on theた．各パラメータは，質量mについては直接

Estimated Calculation

測定を，粘性cおよび剛性kについては自由振本研究では提案する手法の有効性や高精度動波形の実験とMATLABシミュレーションな推定のための同定器設計の指針を探る．の重ね合せにより求めて真値とした．求めた各パラメータの真値を表１に示す．

4.2.1 導入したフィルタリングの影響

4.2 既約分解表現を用いたパラメ
ータ推定と故障診断

対象システムを左既約分解表現に書き換え

対象システムへの入力は，偏心質量を回転る際の必要条件を満たすために導入したローさせたときの水平方向加速度成分による正弦パスフィルタ式(13)の設計パラメータである波入力とし，一定時間の加振を行う．このとカットオフ周波数1/Tが推定結果にどのようき，入力信号は出力が大きくなるように，構に影響するか検討を行った．サンプリング周造物の固有振動数付近の正弦波を用いる．ま波数100[Hz]とし，同定窓を20[sec]として50た出力は質量の変位とし，レーザセンサ測定組の入出力波形を採取し，同定器の設計パラおよびA/D変換後にPCに取り込み，MATLABメータ1/Tを変化させて推定計算を行った結により式(7)と式(10)の計算処理を行った．

果を図５に示す．横軸はカットオフ周波数1/T

パラメータの基準値k0は真値とした．つまを，縦軸は真値からの誤差をパーセント化しり推定計算が良好なほど算出される偏差δkはたものである．また図中の点は50組の推定結０に近くなることになる．逆にもし推定精度果の平均値を表す．図より，採取した入出力が良好な場合に偏差が現れたならば，それは波形の良悪により推定誤差は幅広く帯状に分パラメータの変化を示すことにもなり，構造布してしまう．
しかし50組の平均値をみると，部材それぞれの異常を検知することができる．

システムの固有角周波数の手前で推定結果がまた，本手法は構造パラメータ一つ一つにつ良好なことがわかる．このことからカットオ

いて計算上でそれぞれ同定器を設計すること

フ周波数の値を固有振動数付近かやや手前にが可能であり，実際に必要な測定は１対の入設定して推定計算を行うことで，推定精度を出力波形だけであるため非常に簡易である．

向上させることができる．
-5-

4.2.2 サンプリング長の影響

周波数においても，同定窓が約5[sec]以降になると推定結果がほぼ一定になることがわかる．

実際に機械や構造物の動特性を測定する場

また，サンプリング推定結果が良好な範囲は合には，一時的に設置して行うにしろ，常時200[Hz]付近であり，高精度化を目指してサン設置して測定を連続的に行うにしろ，計測シプリング周波数を大きくすると，逆に推定誤ステムは廉価で，かつ測定が容易であること差が大きくなる結果となった．適したサンプが望ましい．計測を高精度に行おうとしてサリング周波数が存在すると考えることもでき

ンプリング周波数を高めれば，それだけ計測るが，さらなる熟慮が必要である．

システムの演算処理部が高価になる．またデータ長や採取回数が大きくなれば，記憶容量の増大を招く．
そこで本同定法を用いた場合，5. おわりに推定精度に及ぼすデータ長の影響について実1自由度系を対象として，加振入力に対する験的に検討を行った．前述の１自由度系につ構造物の応答実験を行い，提案する同定法のいて，測定サンプリング周波数とデータ長を変え，50組の測定を行った結果を図６に示す．

推定精度向上について検討を行った。得られた結果を以下に示す．

横軸は同定計算処理に用いたデータの時間長(1) 既約化のために導入したフィルタリン

を，縦軸は真値からの誤差の平均値を表す．

グの適した設計により，推定精度の向上

図中，各点は測定するサンプリング周波数ごを見込めることがわかった．との50組の平均値を表す．どのサンプリング(2) データ長は，ある一定時間長さえあれば推定精度が保てることがわかった．
506. 文献

1000[Hz]
500[Hz] 200[Hz]

Error [%]

0100[Hz]1)例えば，山本鎭男(編著)：ヘルスモニタリン

50[Hz]

グ，共立出版(1991)．

-502)大日方，黒沢，川合：日本機械学会論文集C20[Hz]編，70巻，691号，106/109，(2004)．

-100
0510153)K. Zhou, J. C. Doyle, K. Glover: Robust and20

Identification Window [sec]Optimal Control, Prentice-Hall, (1996).
4)S. Boyd, C. Barratt: Linear Controller Design -Fig.6 Effect of Sampling Length on theLimits of performance -, Prentice Hall, (1991).

Estimated Calculation-6-