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

パラメトリックな手法
東京工業大学 機械制御システム専攻
山北 昌毅
確定系のパラメータ推定（１）
ARＸモデルの場合（一括最小自乗法）
y (k )   T (k  1) 0
 T (k  1) : [ y (k  1), y (k  2), , y (k  n), u (k  1), , u (k  n)]
 0T : [ n 1 , ,  0 ,  n 1 , ,  0 ]  R 2 n
 (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)  :  ˆ  Y

eN : 
N 1
N


 




  T


e( N )   ( N  1) 
 y( N )
と し て、 eN の l2ノ ルム を 最小にする こ と を 考える 。 つま り 、
J :|| eN ||2
を 最小化する ˆを 求める 。

J   N 1ˆ  YN
 
T

ˆ  YN  ˆT T N 1 N 1ˆ  2Y T N  N 1ˆ  Y T N YN
N 1
J
 2T N 1 N 1ˆ  2T N 1YN  0
ˆ
1
ˆ   T N 1 N 1  T N 1YN （ ただし 、 T N 1 N 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)
a (k 1)
a2 T (k 1) (k 1)
T
e
(
k
)

e
(
k
)
e(k )
2
c+ T (k 1) (k 1)
c   T (k 1) (k 1)
a T (k 1) (k 1)
||  (k 1) || (2 
) || e(k ) ||2
T
c+ (k 1) (k 1)
c  0, a  0
2
||  (k 1) ||2 (2  a) || e(k ) ||2
a 1
0


確定系のパラメータ推定（３）
ARＸモデルの逐次最小自乗法
(k )
(k  0)
ˆ(0), P(1)  0を 与える
確定系のパラメータ推定（３）’
逆行列の補題
( I  XY )1  I  X ( I  YX )1 Y , ( 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 N 1
P ( N  1) :   (k ) T (k )
N k 0
1 t tN 1
P 1 ( N  1, t ) :
 (k ) T (k )

N k t
[ 弱Ｐ Ｅ 性]
1
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 )}は観測可能
両辺の信号は等し いので、 ' 両辺をフィ ルタ リ ング’ し た
も のも 等し い. こ こ ではフィ ルタ ーと し て
F (s) :
a
(a  0)
sa
F (s)(my(t )  dy(t )  ky(t )  u(t ))  F (s)u(t )
こ こ でsはラ プラ シアンではなく て微分作用素と し て考えている
物理モデルを用いたパラメータ同定（２）
実際の計算では、 例えば
F (s) y(t )の計算には
y f (t ) : ( F (s) y)(t )   ea(t  ) ay( )d
t
0
 [ea(t  ) ay( )]t 0  a ea(t  ) ay( )d （ 部分積分の公式）
t
0
 ay(t )  aeat y(0)  a  ea (t  ) ay( )d
t
0
a
) y(t )
0
sa
こ の式は加速度にフィ ルタ ーをかけた信号は速度信号よ り
近似的に計算でき る こ と を示し ている 。
 ay(t )  a ea(t  ) ay( )d  a(1 
t
ma 1  F (s)  y(t )  dF (s) y(t )  kF (s) y(t )  F (s)u(t )
こ れは形式的に次のよ う にも 解釈でき る
a
a
as
a
y(t ) 
sy(t ) 
y(t )  a(1 
) y(t )
sa
sa
sa
sa

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 cos( y(t )) y(t ) 
m( sin( y(t )) y(t )) y(t )

sa
 sa

a 

x
(
t
)

a
1

 s  a  m cos( y(t )) y(t )
 f 1


 m( x f 1 (t )  x f 2 (t )),

 x (t )  a sin( y(t ) y 2 (t )
f2
sa

推定量の性質
真のパラ メ ータ pと その推定値pˆの誤差を e : p  pˆと する
[ 不偏性]
E{e}  0のと き 、 pˆを pの不偏推定量と いう
[ 一致性]
eNを N 個のデータ から 推定し た際の誤差である と する 。 こ のと き
plim eN  0
N 
と なる と き 、 pˆ Nを pの（ 弱い） 一致推定量と いう 。
[ 有効推定]
不偏推定値の中でも っと も 分散が小さ い推定値
推定値の良さ を評価する には、 不偏性だけでは意味がない！
簡単なシステムの不偏推定値
推定値の誤差e(k )が次の漸化式に従う と する 。 ただし 、 w(k )はホワ イ ト ノ イ ズと する 。
e(k  1)  ae(k )  w(k ),
E{e(0)}  0

a  R ただし 、 E{e2 (0)}  p(0)
E{w2 (k )}   2

k
k
i 1
i 1
e(k )  ak e(0)   ai 1w(k  i)　 E{e(k )}  a k E{e(0)}   ai 1E{w(k  i)}  0
よ って全ての時刻kで, aの値に依ら ずゼロ はe(k )の不偏推定値. 各時刻の分散は
p(k  1) : E{e2 (k  1)}  E{a2e2 (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), , y( k  n), u (k  1), , u (k  n)]
 0T : [ n 1 , ,  0 ,  n 1 , ,  0 ]  R 2 n
T
u (k )と w(k )が独立であ る と する と 、  (k  1)の全ての要素と w(k )は独立と なる .
（ y (n)(n  k )と w(k )には相関があ る ） [ yを 観測し て uを 決定する よ う なフ ィ ード バッ ク 系を 構成し ていないと いう こ と ]
こ の系に最小自乗法を 適用する こ と を 考える
1
ˆ   T N 1 N 1  T N 1YN （ ただし 、 T N 1 N 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 .
 w 
N N 1
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) 
u (1  n) 





u ( N  n) 
T
 w(1) 
 w(2) 


 . 


 . 
 . 


 w( N ) 
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  N N 1
であ れば、 推定値は一致推定量と なる 。 ま た、 期待値と 確率収束の交換性に
よ り 、 漸近的に不偏推定量と も なる 。
補助変数（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), , y (k  n), u (k  1), , u (k  n)]
 0T : [ n 1 , ,  0 ,  n 1 , ,  0 ]  R 2 n
こ れに対し て、 外乱がなかっ たと き に観測さ れる であ ろ う 信号を x(k )と し て
以下のシス テム を 考える 。
x(k )   T (k  1) 0
 T (k  1) : [ x(k  1), x(k  2), , x(k  n), u (k  1), , u (k  n)]
 0T : [ n 1 , ,  0 ,  n 1 , ,  0 ]  R 2 n
こ こ で 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 (k ), u (k )}よ り 補助変数を 用いて同定する 。
p
3. ˆ* ( p)が収束し ていれば終了。 そう でなければ4.へ
4. p  p  1と する 。
5.x (k )を ˆ* ( p  1)と u (k )を 用いて生成し 、 2. へ
p
こ の繰り 返し が一致推定を 与える かど う かは理論的にはわから ない。
出力の予測
確定系の場合
y(k )   T (k  1)0
 T (k  1) : [ y(k  1), y(k  2), , y(k  n), u(k  1), , u(k  n)]
0T : [n1 , ,0 , n1 , , 0 ]  R2n
モデルの構造

 パラ メ ータ  将来の入力から 将来の出力が予測可能
過去の入出力

確率系の場合はどう か？
最小分散推定値
(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 g (Y )  X
 EEf ( X , Y ) Y 
2
 E E g(Y )  X Y 
2
2


 E E g (Y )  EX Y  EX Y  X Y  

 
 E E g  xy  xy  X T g  xy  xy  X Y



 E E g  xy  xy  X 2  2g  xy T xy  X Y
2
期待値をとると
0
23

続き

E g (Y )  X
2
 E E g(Y )  x
2
y
 xy  X
2
Y


2


 E E g (Y )  xy Y   E xy  X

 

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 ) 
1
v(k )
H ( p)
vˆ(k ) 
 H ( p)  1 v(k ) 
H ( p)

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)  y(k )  G( p)u(k ) 


 H 1 ( p)G( p)u(k )  1  H 1 ( p) y(k )
出力の予測（１段先予測器）（３）
特にARMAXモデルでは
y (k )   n 1 y (k  1)   n 2 y (k  2)    0 y (k  n)
  n 1u (k  1)    0 u (k  n)  w(k )   n 1 w(k  1) 
A(q) y(k )  B(q)u(k )  C(q)w(k )
(1   n 1 q 1    0 q  n ) y (k )  (  n 1q 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 ) 
 A(q) 
B(q)
u(k )  1 
 y (k )
C(q)
 C(q) 

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 1q 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  dnd 1 p1
1
w(k )
D( p)
一般化最小自乗法
（GLS: Generalized Least Square）
 d0 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), , u(k  n)]

y(k )
u(k )
0T : [n1 , ,0 , n1 , , 0 ]  R2n

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 , )  dnd 1e(k  1, ) 
e(k  i, )  e(k  i, (k  1))
 d0 e(k  nd , )  w(k ) : d (k  1)d  w(k )  d T (k  1) : [e(k  1, ), , e(k  nd , )]

d T : [dn1 , , d0 ]

参考文献
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（１９９３）