制御系解析基礎 - [yamakita-lab.ctrl.titech]Home

パラメトリックな手法
東京工業大学機械制御システム専攻
山北昌毅
確定系のパラメータ推定（１）
ARＸモデルの場合（一括最小自乗法）
y (k )   T (k  1) 0
 T (k  1) : [ y (k  1), y (k  2),
 0T : [ n 1 ,
,  0 ,  n 1 ,
, y ( k  n), u ( k  1),
, u ( k  n)]
, 0 ]  R2n
 (k  1)を回帰ベクトルと呼ぶ
ˆ k) は
{y( t ) , u( t ) }( t  k  1),  の推定値を ˆとする。この時の、 k時刻での予測出力y(
yˆ(k )   T (k  1)ˆ
となる。ここで、 yを時刻t =0,
, N( Nは十分大きく）まで観測したとすると , e(k ) : yˆ (k )  y (k )として
 e(1)    T (0) 
 y (1) 

 e(2)   T


    (1)  ˆ   y (2)  :  N 1ˆ  YN
eN : 


 




  T


e( N )   ( N  1) 
 y( N )
として、 eN の l2ノルムを最小にすることを考える。つまり、
J :|| eN ||2
を最小化する ˆを求める。

J   N 1ˆ  YN
J
ˆ
 
T

ˆ  YN  ˆT T N 1 N 1ˆ  2Y T N  N 1ˆ  Y T N YN
N 1
 2T N 1 N 1ˆ  2T N 1YN  0
ˆ   T N 1 N 1  T N 1YN （ただし、 T N 1 N 1は正則と仮定した）
1
確定系のパラメータ推定（２）
ARＸモデルの場合（最も簡単な逐次推定法）
一括最小自乗法の場合、データ数Nが大きくなると大きな行列の逆行列を計算しなければならない。
また、適応制御系のように、パラメータの変化に対して、逐次対応するような制御系設計では、オンライン
でパラメータを推定する必要がある。そのため、各時刻毎に予測誤差に基づいてパラメータを更新
[射影法]
[性質]
略証
a (k-1)
[ y (k )   T (k  1)ˆ(k  1)]
c+ (k  1)(k  1)
a (k-1)
ˆ(k )    ˆ(k  1)   
e(k ), e(k ) : y (k )   (k  1)ˆ(k  1)
T
c+ (k  1)(k  1)
a (k  1)
 (k )   (k  1) 
e( k )
T
c+ (k  1)(k  1)
ˆ(k )  ˆ(k  1) 
T
||  (k ) ||2 ||  ( k  1) ||2 2 T ( k  1)
||  (k  1) ||2 (2 
a (k  1) (k  1)
) || e(k ) ||2
T
c+ (k  1) (k  1)
T
c  0, a  0
||  (k  1) ||2 (2  a) || e( k ) ||2
a 1
0
a (k  1)
a 2 T (k  1) (k  1)
T
e
(
k
)

e
(
k
)
c+ T (k  1) (k  1)
c   T (k  1) (k  1)


2
e( k )
確定系のパラメータ推定（３）
ARＸモデルの逐次最小自乗法
(k )
(k  0)
ˆ(0), P(1)  0を与える
確定系のパラメータ推定（３）’
逆行列の補題
( I  XY )1  I  X ( I  YX )1Y , ( A  XY )1  ( I  A1 XY )1 A1
( I  XY ) 1 ( I  XY )  ( I  X ( I  YX ) 1 Y )( I  XY )
 I  XY  X ( I  YX ) 1 Y  X ( I  YX ) 1 YXY
 I  XY  X ( I  YX ) 1 ( I  YX )Y
I
確定系のパラメータ推定（４）
確定系のパラメータ推定（５）
(
持続励振（ＰＥ）条件
（Persistently Excitation）
1


1
 1

ˆ( N )    T N 1 N 1  P0 1   
P0 1ˆ(0)   N 1T YN 
N
 N

1
P ( N  1) :
N
1
P 1 ( N  1, t ) :
N 1
  (k )
T
(k )
k 0
1
N
t  tN 1
  (k )
T
(k )
k t
[ 弱ＰＥ性]
N  0に対して、 P 1 ( N  1)が正則であるとき入力（列）は弱ＰＥ性があるという
[ 強ＰＥ性]
ある大きさ N  0を固定した時、任意の初期時間tに対して P 1 ( N  1, t )が正則のとき、
入力（列）は強ＰＥ性があるという
物理モデルを用いたパラメータ同定（１）
my (t )  dy (t )  ky (t )  u (t )
仮定：{ y (t ), y (t ), u (t )}は観測可能
両辺の信号は等しいので、 ' 両辺をフィルタリング ’ した
ものも等しい. ここではフィルターとして
a
F ( s ) :
(a  0)
sa
F ( s )(my (t )  dy (t )  ky (t )  u (t ))  F ( s )u (t )
ここで sはラプラシアンではなくて微分作用素として考えている
物理モデルを用いたパラメータ同定（２）
実際の計算では、例えば
F ( s ) y (t )の計算には
t
y f (t ) : ( F ( s ) y )(t )   e  a ( t  ) ay ( )d
0
t
 [e  a ( t  ) ay ( )]t 0  a  e  a ( t  ) ay ( )d （部分積分の公式）
0
t
 ay (t )  ae  at y (0)  a  e  a ( t  ) ay ( )d
0
a
) y (t )
0
sa
この式は加速度にフィルターをかけた信号は速度信号より
t
 ay (t )  a  e  a ( t  ) ay ( )d  a (1 
近似的に計算できることを示している。
これは形式的に次のようにも解釈できる
a
a
as
a
y (t ) 
sy (t ) 
y (t )  a (1 
) y (t )
sa
sa
sa
sa
ma 1  F ( s )  y (t )  dF ( s ) y (t )  kF ( s ) y (t )  F ( s )u (t )

m
 a 1  F ( s )  y (t ) F ( s ) y (t ) F ( s ) y (t )   d   F ( s )u (t )
 
 k 
物理モデルを用いたパラメータ同定（３）
微分方程式は線形でなくとも良い。例えば、
m cos( y (t )) y (t )  dy (t )  ky (t )  u (t ))  F ( s )u (t )
F ( s ) :
a
(a  0)
sa
F ( s )(m cos( y (t )) y (t )  dy (t )  ky (t )  u (t ))  F ( s )u (t )
a
as
a
m cos( y (t )) sy (t )) 
m cos( y (t )) y (t ) 
m( s cos( y (t )) y (t )
sa
sa
sa
a 
a

 a 1 
m(  sin( y (t )) y (t )) y (t )
 m cos( y (t )) y (t ) 
sa
 sa

a 

x
(
t
)

a
1

f
1

 m cos( y (t )) y (t )

s

a


 m( x f 1 (t )  x f 2 (t )),

 x (t )  a sin( y (t ) y 2 (t )
 f 2
sa
推定量の性質
真のパラメータ pとその推定値pˆの誤差を e : p  pˆとする
[ 不偏性]
E{e}  0のとき、 pˆを pの不偏推定量という
[ 一致性]
eNを N 個のデータから推定した際の誤差であるとする。このとき
p lim eN  0
N 
となるとき、 pˆ Nを pの（弱い）一致推定量という。
[ 有効推定]
不偏推定値の中でもっとも分散が小さい推定値
推定値の良さを評価するには、不偏性だけでは意味がない！
簡単なシステムの不偏推定値
推定値の誤差e(k )が次の漸化式に従うとする。ただし、 w(k )はホワイトノイズとする。
e(k  1)  ae(k )  w(k ),
a  R ただし、
 E{e(0)}  0

2
 E{e (0)}  p (0)
 E{w2 (k )}   2

k
k
i 1
i 1
e(k )  a k e(0)   a i 1w(k  i )　 E{e(k )}  a k E{e(0)}   a i 1E{w(k  i )}  0
よって全ての時刻kで, aの値に依らずゼロは e(k )の不偏推定値. 各時刻の分散は
p (k  1) : E{e 2 (k  1)}  E{a 2e 2 (k )  2ae(k ) w(k )  w2 (k )}  a 2 p (k )   2
p (k )有界は一定値pに収束したとすると
p  a2 p   2

p
2
1  a2
pはゼロ以上であるので、
1  a2  0

| a | 1
| a | 1である場合は pは発散する。しかし、その場合でも誤差の期待値は０である！
確率変数の収束に関する性質
1.連続関数f ()に関して確率収束と概収束の順番は交換可能
plim f ( xn )  f (plim xn )
 n
n 

f ( xn )  f (a.s.lim xn )
a.s.lim
n 
n 
行列A, Bのサイズがnに関して不変の場合（変化する場合は一般に成り立たない）
plim An Bn  plim An plim Bn
n 
n 
 n
1



1
plim A n   plim An 
 n

 n
2.確率収束, 概収束, 自乗平均収束と期待値演算は交換可能
lim E{xn }  E{plim xn }
n 
 n

E{xn }  E{a.s.lim xn }
lim
n 
n 

E{xn }  E{l.i.m. xn }
lim
n 
n 
確率系のパラメータ同定（１）
外乱のあるARＸモデルの場合（一括最小自乗法）
y (k )   T (k  1) 0  w(k ),
w(k )は平均値ゼロの白色ノイズ
 (k  1) : [ y (k  1), y ( k  2),
T
 0T : [ n 1 ,
,  0 ,  n 1 ,
, y ( k  n), u ( k  1),
, u ( k  n)]
, 0 ]  R 2n
u (k )と w(k )が独立であるとすると、  (k  1)の全ての要素と w(k )は独立となる .
（ y (n)(n  k )と w(k )には相関がある） [ yを観測して uを決定するようなフィードバック系を構成していないということ ]
この系に最小自乗法を適用することを考える
ˆ   T N 1 N 1  T N 1YN （ただし、 T N 1 N 1は正則と仮定した）
1
 y (1) 
 y (2) 
 であるが、今の場合次のように書き表される
YN  




 y( N )
 w(1) 
1
 w(2) 
 :  N 1 0  w　 ˆ   0   T N 1 N 1 1 N 1 T N 1w   0   1 T N 1 N 1  1 T N 1w
YN   N 1 0  




N
N
 N


 w( N ) 
確率系のパラメータ同定（２）
ここで、
 y (0)
 y (1)

.
1 T
1
 N 1w  
.
N
N

.

 y ( N  1)
y (1)
.
 y (0)


1  y (n) y (1  n)
 
u (1)
N  u (0)


 u (n) u (1  n)
y (  n)
y (1  n)
u (0)
u (1)
y ( N  n) u ( N  1)
.
.
u ( n)   w(1) 
u (1  n)   w(2) 
  . 
 

  . 
  . 
 

u ( N  n)   w( N ) 
T
y ( N  1)   w(1) 
  w(2) 


y ( N  n)   . 


u ( N  1)   . 
 . 


u ( N  n)   w( N ) 
上記の演算で y (i ) w( j )に関してはi  jの演算しか含まず、 w()と u ()は独立であるとしているので、
u (k ), y (k )の期待値が有界であれば
1

1

E  T N 1w  E  T N 1  E{w}  0
N

N

1
となり、 T N 1
N
が収束すれば( 他の変数と独立になり）、 ˆは不偏推定量となる。
N 1
確率系のパラメータ同定（３）
また
1
1
1

plim ˆ   0  plim  T N 1 N 1  plim T N 1w
N
N 
N   N

N 
であるので、
1

1 T

plim   N 1 N 1  が存在する
 N   N


plim 1 T w  0
N 1
 N  N
であれば、推定値は一致推定量となる。また、期待値と確率収束の交換性に
より、漸近的に不偏推定量ともなる。
補助変数（Instrumental Variable)法(1)
w(k )が白色でなく、相関がある場合について考える。最小自乗法では推定値を以下で決定した。
ˆ   T N 1 N 1  T N 1YN
1
補助変数法では、 T N 1を行列のサイズが同じ V Tに置き換える。つまり、
ˆ*  V T  N 1  V T YN
1
このようにした場合でも、以下の条件がなり立てば一致推定量となる。
1

1 T

plim  V  N 1  が存在する
 n  N


plim 1 V T w  0
 n N
補助変数（Instrumental Variable)法(2)
雑音に汚されていない信号を用いる方法
実際のシステムは外乱w(k )に汚されて以下のように表される
y (k )   T (k  1) 0  w(k )
 T (k  1) : [ y (k  1), y (k  2),
 0 : [ n 1 ,
T
,  0 ,  n 1 ,
, y (k  n), u (k  1),
, 0 ]  R
, u ( k  n)]
2n
これに対して、外乱がなかったときに観測されるであろう信号を x(k )として
以下のシステムを考える。
x(k )   T (k  1) 0
 T (k  1) : [ x(k  1), x(k  2),
 0 : [ n 1 ,
T
,  0 ,  n 1 ,
, x( k  n), u ( k  1),
, 0 ]  R
, u ( k  n)]
2n
ここで u (k )が w(k )と独立であるとすると、明らかに x(k )はw(k )と独立となり、入力がＰＥ条件を
満たせば補助変数の条件を満たす。この式は x(k )の漸化式であるので、 x(0), x( 1),
, x( n)が与え
られれば、 u (k )は既知であるので計算可能である。
一番簡単な方法では、 ˆを求めることと、 x(k )を求めることを交互に繰り返す。つまり、
1.繰り返しの回数p  1とする。 y (k )を x(k )にセットする。これを x p (k )と表す。
2.ˆ* ( p )を{x p (k ), u (k )}より補助変数を用いて同定する。
3. ˆ* ( p )が収束していれば終了。そうでなければ4.へ
4. p  p  1とする。
5.x p (k )を ˆ* ( p  1)と u (k )を用いて生成し、 2. へ
この繰り返しが一致推定を与えるかどうかは理論的にはわからない。
出力の予測
確定系の場合
y (k )   T (k  1)0
 T (k  1) : [ y(k  1), y (k  2),
0T : [ n 1 ,
,  0 ,  n 1 ,
, y ( k  n), u (k  1),
, u ( k  n)]
, 0 ]  R 2n
モデルの構造

 パラメータ  将来の入力から将来の出力が予測可能
 過去の入出力

確率系の場合はどうか？
最小分散推定値
( xとは独立ではない)
推定したい
パラメータ x
観測量 y p( y x)
推定量 xˆ
推定ルール
g ( y)

評価関数 E x  xˆ
2
 を最小にする推定値
条件付き期待値
xˆ  Ex y と等価
（ x の分布の種類によらず）
22
•証明
E f ( X , Y )   f ( x, y ) p( x, y )dxdy
  f ( x, y ) p ( x y )dx p ( y )dy
 EEf ( X , Y ) Y 

E g (Y )  X
2
  E E g(Y )  X Y 
2
2


 E  E  g (Y )  EX Y  EX Y  X Y  

 

 E E g  x y  x y  X  g  x y  x y  X Y

 E E g  xy
T
2

 x y  X  2g  x y  x y  X Y
2
T
期待値をとると
0
23

続き

E g (Y )  X
2
  E  E  g (Y )  x
 E 

y
 xy  X
2
Y


2

E  g (Y )  x y Y   E x y  X


2

Y 

g (Y )と無関係
これが最小になるのは
g (Y )  xy
2
の時。
つまり、
g (Y )  EX Y  が最小分散推定値となる
24
出力の予測（１段先予測器）（１）
y (k )  G ( p )u (k )  H ( p ) w(k )  G ( p )u (k )  v (k )
であるシステムを考える。

v(k ) : H ( p ) w(k )  w(k )  H ( p )v(k ) :  h (i )v(k  i ),
1
i 0


 | h (i) | 
i 0
v(i )(i  k  1)が既知とする。 w(k )  v(k )h (0)   h (i )v(k  i )  v(k ) 
i 1
w(k )  h (i )

v (k  i )
h (0) i 1 h (0)
従って、 y (k )の v(i )(i  k  1)の条件付き期待値yˆ (k )は

E{w(k )}  h (i )
h (i )
yˆ (k )  G ( p )u (k ) 

v (k  i )  G ( p )u (k )  
v (k  i ) : G ( p )u (k )  vˆ(k )
h (0)
i 1 h (0)
i 1 h (0)
v(k ) 
w(k )
 vˆ(k )
h (0)
vˆ(k )   w(k ) / h (0)  v(k )   w(k ) / h (0)  H ( p ) w(k )  ( 1 / h (0)  H ( p )) w(k )
出力の予測（１段先予測器）（２）

特に H ( p)  1   h(i )q  iの場合は表現が簡単になる  h (0)  1
i 1
vˆ(k )   H ( p )  1 w(k )
w(k ) 
vˆ(k ) 
1
v (k )
H ( p)
 H ( p)  1
H ( p)


v(k )  1  H 1 ( p ) v(k )



yˆ (k )  G ( p)u (k )  1  H 1 ( p) v(k )  G ( p)u ( k )  1  H 1 ( p)


 H 1 ( p)G ( p)u (k )  1  H 1 ( p) y (k )
  y(k )  G( p)u(k ) 
出力の予測（１段先予測器）（３）
y (k )   n 1 y (k  1)   n  2 y (k  2) 
特に ARMAXモデルでは
  n 1u (k  1) 
A(q ) y (k )  B(q )u (k )  C (q) w(k )
y (k ) 
B(q )
C (q )
u (k ) 
w(k )
A(q )
A(q)
yˆ (k ) 

B(q )
A(q ) 
u (k )  1 
 y (k )
C (q )
 C (q ) 

(1   n 1 q 1 
  0 y ( k  n)
  0 u (k  n)  w(k )   n 1 w(k  1) 
  0 q  n ) y (k )  (  n 1 q 1 
A(q) y (k )  B(q)u (k )  C ( q) w(k )
G ( q) 
B(q)
,
A( q)
H ( q) 
C (q)
A( q)
yˆ (k )  B(q )u ( k )   C (q )  1  1  A( q)  y ( k )   C ( q )  1 yˆ ( k )
 B (q )u (k )  1  A(q)  y (k )   C (q )  1 y (k )  yˆ ( k ) 
出力の予測にはパラメータが必要（モデル構造は既知としても）
yˆ (k )に関する漸化式（フィルター）になっている
  0 w(k  n)
  0 q  n )u (k )  (1   n 1 q 1 
  0 q  n ) w(k )
予測誤差
y (k )  G ( p)u (k )  H ( p ) w(k )  G ( p )u (k )  v(k )  G ( p)u (k )  H ( p ) w(k )
yˆ (k )  G ( p)u (k )  vˆ(k )  G ( p )u (k )  ( H ( p )  1) w( k )
y (k )  G ( p)u (k )  ( H ( p)  1) w(k )  w(k )
e(k ,  ) : y (k )  yˆ (k )  w(k )
e(k ,  )を予測誤差と呼ぶ。最小分散推定値との予測誤差はw(k )であるので、予測誤差は白色信号
となっている。（パラメータが正しい場合！）
確率系の逐次パラメータ同定（１）
ARMAXモデルのパラメータ同定（拡張最小自乗法）
y (k )  G ( p )u (k )  H ( p ) w(k )  G ( p )u (k )  v(k )
であるシステムを考える。 v(k )は有色であるので、 v(k )を白色信号であるとして
推定するだけでは一般に一致推定量を得ることができない。しかし、予測モデルを
用いることにより、誤差を白色化することができる。
yˆ (k )  B(q )u (k )  1  A(q)  y (k )   C (q )  1 y (k )  yˆ (k ) 
   y (k  1)
   y (k  1)
 y (k  n) u (k  1)
 y (k  n) u (k  1)
  n 1 




 0 
u ( k  n)  
  e(k  1,  )
  n 1 




 0 
u (k  n) e( k  1,  )
 n 1 

e(k  n,  )  

  0 
 n 1 




  n 1 
e( k  n) 



  n 1 




:  T (k  1)
y (k )  yˆ (k )  w(k )   T (k  1)  w(k )
この式はＡＲＸモデルの場合と同様の式をしている。この式に、繰り返し最小自乗法を適用する。
確率系の逐次パラメータ同定（２）
(k )
ただし、
e(k ) : y (k )   T (k  1)ˆ(k  1)

T
e(k  i,  )  y (k  i )   (k  i  1)ˆ(k  1)
1
収束性はC ( z )  が強正実( St r i ct l y Posi t i ve Real : SPR) であれば一致推定になる
2
ことが知られているが、ここでは説明を省略する
確率系の逐次パラメータ同定（２）
ノイズが特殊なモデルで表現できる場合（ＡＲＡＲＸモデル）
A( p ) y (k )  B( p )u (k ) 
D( p )  1  d nd 1 p 1
1
w(k )
D( p)
一般化最小自乗法
（GLS: Generalized Least Square）
 d 0 p  nd
y (k )   T (k  1) 0  w(k )
A( p ) D( p ) y (k )  B ( p ) D( p )u (k )  w( k )   T (k  1) : [ y (k  1), y (k  2), , y ( k  n), u ( k  1),

y(k )
u (k )
 0T : [ n 1 , ,  0 ,  n 1 , ,  0 ]  R 2 n

, u ( k  n)]
A( p ), B ( p )が既知とすれば
1
1
A( p ) y (k )  B( p )u (k ) 
w(k ) 
w(k )  A( p ) y (k )  B ( p )u (k ) : e(k ,  )
D( p)
D( p)
D( p )e(k ,  )  w(k )
e(k ,  )  d T (k  1) d  w(k )
e(k ,  )  d nd 1e(k  1,  ) 
e(k  i,  )  e(k  i,  (k  1))
 d 0 e(k  nd , )  w(k ) : d (k  1) d  w(k )  d T (k  1) : [e(k  1,  ), , e( k  nd ,  )]

 d T : [d n 1 , , d 0 ]

参考文献
1.
2.
3.
4.
L.Lung:System Identification,Prentice Hall
PTR(1987)
G.C.Coodwin,K.S.Sin: Adaptive Filtering
Prediction and Control, Prentice-Hall(1984)
相良ら：システム同定、コロナ社（１９９５）
足立：ユーザのためのシステム同定理論、ＳＩＣ
E（１９９３）

Download Report