計算機科学特論 I レポート課題1

線形システムの同定
はじめに

1

システムの解析，予測，制御，シミュレーション及び異常診断のため，システムの動特性を知る必要がある．システムの入出力信号の観測値に基づいて，ある評価規範のもとでそのシステムをもっともよく記述する数学モデルを決定することをシステム同定という．このような統計学に基づいた同定法は大量かつ複雑なデータ処理を必要とするが，ディジタルコンピュータの発達はこのようなデータ処理を可能とした．したがって，この 20 年間に，システムの解析，予測，制御などの基礎となるシステム同定理論の進歩はめざましく，成熟しつつある．現在では，最小二乗法で代表されるシステムのパラメータ推定手法は実システムに使えるようになり，工学の基礎技術の一つに数えられるようになった．
本演習は，2 次の離散時間線形システムのパラメータ推定問題をシミュレーションによって検討する．

2

最小二乗法による線形システムのパラメータ推定

2.1

最小二乗法の導出

Fig. 1 に示すような線形システムの入出力関係を表す数式モデルの記述について考える．

Figure 1: 線形システムの同定のためのモデル.

線形システムの次数 n を既知と仮定すると，その入出力関係は次式のような差分方程式で表される．
nn∑∑xt = −
ai xt−i +bi ut−i(1)i=1i=11

この資料は，学部 3 年以上のレベルの学生が，講義を直接受けなくても，読んでわかるように配慮して作成したものです。最後まで落ち着いて読んでから始めてください。

1ただし，ut と xt はそれぞれシステムの入出力信号の離散時刻 t におけるサンプル値であり，aiと bi はシステムモデルのパラメータである。
ここで，出力信号 xt には，平均値ゼロの離散値確率雑音 vt が加わると仮定する．

y t = xt + v t(2)一般に，観測雑音 vt はシステムの入力信号 ut とは，無相関である．すなわち，

E[ut vt ] = 0(3)ただし，E[·] は期待値を意味する．
我々の目的は，システムの入出力の観測値 ut と yt からシステムのパラメータを推定することである．
ここで，(2) 式を (1) 式に代入すると，

yt = −n∑ai yt−i +i=1n∑bi ut−i + rt(4)i=1∑nが得られる．ただし，rt = i=0 ai vt (a0 = 1) は観測雑音による式誤差である．
上式はまた次のようなベクトル形式で表現できる．
yt = z Tt θ + rtz Tt = [−yt−1 , · · · , −yt−n , ut−1 , · · · , ut−n ](5)θ T = [a1 , · · · , an , b1 , · · · , bn ]システムのパラメータを推定するため，次のような二乗誤差評価を導入する．

J=N∑

ρN −t rt2 =t=1N∑

ρN −t (yt − z Tt θ)2(6)t=1

ただし，0 < ρ ≤ 1 は忘却係数である．N を現在時刻とみれば，ρN −t は過去値になるほど，二乗誤差に対する重みを小さくして推定値への寄与を小さくする役割を果たしており，一種の忘却機能を持っている．この忘却の効果のため，ゆっくり変動する時変パラメータの推定にも適用できる．
J を最小にする必要条件は ∂J/∂θ = 0 である．すなわち，
NN∑∑∂J= −2
ρN −t z t yt + 2
ρN −t z t z Tt θ = 0∂θ
t=1t=1(7)より，θ の最小二乗推定値が得られる．

b=
θ( N∑

ρN −t z t z Tt)−1 N∑

t=1

上式が成り立つためには，

t=1( N∑)−1
ρN −t z t z Ttt=12

ρN −t z t yt(8)が存在しなければならない．一般に，入力信号が十分多くの周波数成分を含んでいれば，上記の逆行列が存在する．(8) 式が J を最小化しているのを示すのは容易である．
J = θ −(N∑ρN −t z t z Tt)−1 N
∑

t=1+N
∑
t=1

ρN −t yt2 −(ρN −t z t yt N∑)
(N)−1 N∑∑
ρN −t z t z Tt θ −
ρN −t z t z Tt
ρN −t z t yt 

t=1t=1(N∑)T ( N∑ρN −t z t yt

t=1

ρN −t z t z Tt

(9)b は確かに上式を最小にしている．と書けるので，θ

2.2

ρN −t z t yt

t=1t=1t=1)−1 N∑t=1

推定値の漸近的性質

データ数が多くなるにつれて，確率極限

b=θplim θ(10)N →∞b は一致推定量であるという．
が成立するとき，θ
ここでは，忘却係数 ρ = 1，すなわち過去のデータを忘却しない場合について考察する．(8)式の推定値を (5) 式に代入して，確率極限をとると，()−1 ()NN∑∑{ T}
11Tb =plim
plim θztzt
plimz t z t θ + rt
N →∞ N t=1N →∞N →∞ N t=1(11)()−1NN∑
1 ∑1= θ + plimz t z Tt
plimz t rtNN →∞N →∞ Nt=1t=1b が一致推定量であるためには，ベクトル z t と式誤差 rt は無相関が得られる．したがって，θでなければならない，すなわち，
1 ∑z t rt = 0
N →∞ Nt=1Nplim(12)という条件が成り立たなければならない．一般に上式が成立することが期待できないので，観測雑音のレベルが高い場合，最小二乗推定は一致推定量でなく，偏り（バイアス）を持つ．

2.3

補助変数法によるパラメータの一致推定ベクトル z t と同サイズの補助変数ベクトル mt を導入し，次のように補助変数法による推定値を得る．
( N)−1 N∑∑bIV =(13)θρN −t mt z Tt
ρN −t mt yt
t=1t=1

そこで，忘却係数 ρ = 1，すなわち過去のデータを忘却しない場合における補助変数推定値の漸近的性質について考察する．上式の推定値を (5) 式に代入して，確率極限をとると，())−1 (NN∑∑}{
11b =plim θplimmt z Ttplimmt z Tt θ + rt
N →∞N →∞ N t=1
N →∞ N t=1(14)()−1NN∑∑
11= θ + plim
plim
mt z Tt
mt r tN →∞ NN →∞ Nt=1t=13が得られる．
以上より明らかに，補助変数ベクトル mt が次のような条件を満足すれば，パラメータ推定値は一致推定量である．
（１）逆行列)−1(N1 ∑T
mt z tplim
N →∞ N t=1が存在する．
（２）補助変数ベクトル mt と式誤差 rt は無相関である．すなわち，

1 ∑mt r t = 0
N →∞ Nt=1N

plim(15)が成り立つ．
補助変数ベクトルの幾つかの構成：（１）遅延された入力を利用する．

mt = [ut−d , · · · , ut−d−n , ut−1 , · · · , ut−n ](16)

ただし，d > n．
（２）遅延された出力を利用する．

mt = [yt−1−d , · · · , yt−n−d , ut−1 , · · · , ut−n ](17)

ただし，d の選び方は，明確ではない．
（３）推定された出力を利用する．

mt = [−bxt−1 , · · · , −bxt−n , ut−1 , · · · , ut−n ]ただし，

xbt = −n∑bai xbt−i +i=1n∑

bbi ut−i(18)

(19)

i=1

ここで，bai 及び bbi は何らかの方法で得たパラメータの推定値である（一致推定量でなくてもいい）．
本演習は (3) 番目の構成法を利用する．詳細は後述する．

3逐次推定アルゴリズム

前節で導出した最小二乗法および補助変数法は，データを蓄えておき，データ行列の逆行列を計算して推定値を得るバッチ処理である．このような方法は繰り返し計算を必要としないことから one-shot 法またはオフライン法とも呼ばれている．しかし，システムの特性を実時間での監視や，異常の検出と診断には不向きである．
bN とする．そこで，θ
bN を過去の推定値及び時いま，N 個のデータから得られる推定値を θ
刻 N でのデータより求めることを考える．新しくデータが得られる度に直前の推定値を修正していく方式（逐次計算式）が実現されれば，オンライン推定や実時間推定への適用が可能である．
具体的導出過程は紙面の制限のため，参考文献に委ねたいが，最小二乗法と補助変数法に関する逐次推定アルゴリズムは以下のようになる．

4

bN = θbN −1 + LN εN
θbN −1εN = yN − φTN θP N −1 ψNTρN +[ φN P N −1 ψN]
1P N −1 ψN φTN P N −1=
P N −1 −ρNρN + φTN P N −1 ψN

LN =PN(20)ただし，P N は 2n × 2n の行列で，LN は 2n × 1 のベクトルである．
上式において，逐次最小二乗法は次のように与えられる．

φN = z N ,ψN = z N(21)

また，逐次補助変数法は次のように与えられる．

φN = z N ,ψN = m N(22)

ここで，下記の点に留意しておくべきである．
bN は，1 時刻前の推定値 θbN −1 に，式誤差の推定値 yN − φT θb（１）推定値 θN N −1 に比例した修正項を加えることによって生成される．したがって，バッチ処理に比べて記憶容量が少なくすむ．また，逐次計算では逆行列の計算を必要としない．
（２）上述の逐次アルゴリズムを実行するためには初期値の設定が必要である．通常，次のように選べばよい．
b0 = Oθ(23)b 0 = αIP
ただし，α = 103 ∼ 105 ．
（３）推定された出力を使って補助変数を構成する場合，推定されたシステムの安定性をチェックしなければならない．ここで，補助変数ベクトルは次のように構成されるとする．

mN = [−bxN −1 , · · · , −bxN −n , uN −1 , · · · , uN −n ]ただし，xbN = −n∑āi,N −1 xbt−i +i=1n∑bbi,N −1 ut−i

(24)

(25)

i=1

ここで，bbi,N −1 は (20) 式の逐次推定アルゴリズムから得られた時刻 N − 1 における推定値である．āi,N −1 は，安定性を考慮して，(20) 式の逐次推定アルゴリズムから得られた時刻 N − 1 における推定値 bai,N −1 を修正したものである．時刻 N において，推定出力 xbN を計算するとき，上記のシステムの極が安定でなければ，値が発散してしまう．したがって，システムの極，すなわち，z n + ā1,N −1 z n−1 + · · · + ān,N −1 = 0(26)の解が単位円の内部に入れば（複素数の半径が 1 より小さい），システムが安定である．
時刻 N において，推定出力 xbN を計算するとき，bai,N −1 (i = 1, 2, · · · , n) で構成されたシステムの極が不安定であった場合，対策として，bai,N −1 (i = 1, 2, · · · , n) を安定な領域に射影することによって修正する．すなわち，bai,N −1 (i = 1, 2, · · · , n) が安定な極を構成していれば，āi,N −1 = bai,N −1 (i = 1, 2, · · · , n) とする．そうでなければ，bai,N −1 (i = 1, 2, · · · , n) を安定な領域へ修正する．アルゴリズムは下記のとおりである．

5

Step a: If bai,N −1 (i = 1, 2, · · · , n) stable thenāi,N −1 = bai,N −1 (i = 1, 2, · · · , n)Stop
Elseāi,N −1 = āi,N −2 + 0.5 ∗ (bai,N −1 − āi,N −2 )EndifStep b: If āi,N −1 (i = 1, 2, · · · , n) stable thenStopElseāi,N −1 = āi,N −2 + 0.5 ∗ (āi,N −1 − āi,N −2 )EndifStep c: Go to Step b.
注意：上述の修正は，推定出力の計算だけのために行われることに注意されたい．(20) 式の逐次推定アルゴリズムに対しては，推定されたシステムの安定性の如何に関わらず，このような修正をしない．
このように，システムの出力とパラメータを交互に推定することは Bootstrap（靴のひもかけ）法という．
（４）忘却係数 ρN も (20) 式の推定アルゴリズムの性能を大きく左右する．一般に，忘却係数が小さいほど，忘却機能が強く働くので，推定速度が速く，パラメータの変動にも対応できるが，システムに関する情報を速く忘却するので，観測雑音に弱く，推定値が振動的になる．通常，時不変システムのパラメータ推定において，ρN = 1 としても，差し支えないが，忘却係数を時間の経過とともに 1 に近づくように設定した方が，パラメータの推定値が速く真値へ収束する．
ρN = (1 − 0.01)ρN −1 + 0.01, ρ0 = 0.95(27)すなわち，最初の内に忘却係数を小さくすることによって推定速度を上げるが，推定値がある程度その真値に近づいたら，雑音の影響を軽減するために，忘却係数を 1 とすることによって推定速度を落とす．
一方，時変システムのパラメータ推定において，ρN = 0.97 ∼ 0.999 と，1 よりやや小さい定数にすればよい．忘却係数が小さいほど，時変パラメータへの追従速度が速くなるが，雑音の影響を受けやすくなる．いずれにしても，雑音のレベルが高い場合，時変システムのパラメータ推定は困難である．

4

シミュレーション実験

次のような 2 次のシステムを考える．

xt = −a1 xt−1 − a2 xt−2 + b1 ut−1 + b2 ut−2
y t = xt + v t

(28)

ただし，

a1 = −1.5,

a2 = 0.7,

b1 = 1.0,

b2 = 0.5(29)

ここで，入力信号を [0, 1] の範囲で一様分布する乱数とする．観測雑音はゼロ平均（確認すること）の乱数で，入力信号と無相関（確認すること）である．
雑音のレベルは一般に雑音信号比（Noise Signal Ratio, NSR）で評価される．ここで，NSRは雑音 vt の標準偏差と真の出力信号 xt の標準偏差との比である．

6

上記の要領でシステムのシミュレータを作成すること．ただし、信号の初期値をすべてゼロとする．特に，雑音に課する仮定条件を確認すること．同じ実験をもう一回確認する場合，乱数の種数を前回のと同じでなければ，前回の結果が再現できないことに注意すること．

4.1

逐次最小二乗法によるパラメータ推定

同じ入力信号に対して，出力の雑音信号比 NSR をそれぞれ 0%，10%，20% とし，入出力信号の観測値 {ut , yt }(t = 1, 2, · · · , N ) を用いて，パラメータを逐次最小二乗法で推定せよ．ただし，N は十分大きな整数（2000 ∼ 3000）とする．時刻を横軸として，各推定パラメータの挙動を，パラメータの真値とあわせて図示せよ．観測雑音によるパラメータ推定値への影響を考察せよ．

4.2

逐次最小二乗法による時変系のパラメータ推定

次のような 2 次の時変システムを考える．

xt = −a1,t xt−1 − a2,t xt−2 + b1,t ut−1 + b2,t ut−2
y t = xt + v t
ただし，

(30)

a1,t = −1.0 + 0.5 sin(2πt/2000)
a2,t = 0.7(31)

b1,t = 1.0b2,t = 0.5 + 0.4 cos(2πt/2000)

同じ入力信号に対して，出力の雑音信号比 NSR をそれぞれ 0%，5%，10% とした場合，逐次最小二乗法で上記のシステムの時変パラメータを推定せよ．時刻を横軸として，各推定パラメータの挙動を，パラメータの真値とあわせて図示せよ．忘却係数及び観測雑音によるパラメータ推定値への影響を考察せよ．

4.3

逐次補助変数法によるパラメータ推定

逐次最小二乗法によるパラメータ推定の実験と同じ実験条件で，パラメータを逐次補助変数法で推定せよ．時刻を横軸として，各推定パラメータの挙動を，パラメータの真値とあわせて図示せよ．補助変数法と最小二乗法による結果を比較し，検討せよ．

参考文献
相良，秋月，中溝，片山 (1981)：システム同定（計測自動制御学会）.
中溝高好 (1988)：信号解析とシステム同定（コロナ社）．
相良，和田，中野 (1992)：ディジタル制御の基礎（コロナ社）.
足立修一 (1993)：ユーザーのためのシステム同定理論（計測自動制御学会）.
片山徹 (1994)：システム同定入門（朝倉書店）.

7