パラメトリックな手法
東京工業大学 機械制御システム専攻
山北 昌毅
確定系のパラメータ推定(1)
ARXモデルの場合(一括最小自乗法)
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
確定系のパラメータ推定(2)
ARXモデルの場合(最も簡単な逐次推定法)
一括最小自乗法の場合、 データ 数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 )
確定系のパラメータ推定(3)
ARXモデルの逐次最小自乗法
(k )
(k  0)
ˆ(0), P(1)  0を 与える
確定系のパラメータ推定(3)’
逆行列の補題
( 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
確定系のパラメータ推定(4)
確定系のパラメータ推定(5)
(
持続励振(PE)条件
(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
[ 弱P E 性]
N  0に対し て、 P 1 ( N  1)が正則であ る と き 入力( 列) は弱P E 性があ る と いう
[ 強P E 性]
あ る 大き さ N  0を 固定し た時、 任意の初期時間tに対し て P 1 ( N  1, t )が正則のと き 、
入力( 列) は強P E 性があ る と いう
物理モデルを用いたパラメータ同定(1)
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はラ プラ シア ン ではなく て微分作用素と し て考えている
物理モデルを用いたパラメータ同定(2)
実際の計算では、 例え ば
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 
物理モデルを用いたパラメータ同定(3)
微分方程式は線形でなく と も 良い。 例えば、
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は発散する 。 し かし 、 その場合でも 誤差の期待値は0 であ る !
確率変数の収束に関する性質
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 
確率系のパラメータ同定(1)
外乱のあるARXモデルの場合(一括最小自乗法)
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 ) 
確率系のパラメータ同定(2)
こ こ で、
 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
確率系のパラメータ同定(3)
また
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 )と 独立と なり 、 入力がP E 条件を
満たせば補助変数の条件を 満たす。 こ の式は 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
出力の予測(1段先予測器)(1)
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 )
出力の予測(1段先予測器)(2)

特に 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 ) 
出力の予測(1段先予測器)(3)
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 )であ る ので、 予測誤差は白色信号
と な っ ている 。 ( パラ メ ータ が正し い 場合!)
確率系の逐次パラメータ同定(1)
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 )
こ の式はA R X モデルの場合と 同様の式を し ている 。 こ の式に、 繰り 返し 最小自乗法を 適用する 。
確率系の逐次パラメータ同定(2)
(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
こ と が知ら れている が、 こ こ では説明を 省略する
確率系の逐次パラメータ同定(2)
ノイズが特殊なモデルで表現できる場合(ARARXモデル)
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)
相良ら:システム同定、コロナ社(1995)
足立:ユーザのためのシステム同定理論、SIC
E(1993)
ダウンロード

制御系解析基礎 - 東京工業大学