制御系解析基礎 - 東京工業大学
Download
Report
Transcript 制御系解析基礎 - 東京工業大学
パラメトリックな手法
東京工業大学 機械制御システム専攻
山北 昌毅
確定系のパラメータ推定(1)
ARXモデルの場合(一括最小自乗法)
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
2T N 1 N 1ˆ 2T N 1YN 0
ˆ
1
ˆ T N 1 N 1 T N 1YN ( ただし 、 T N 1 N 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)
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
確定系のパラメータ推定(3)
ARXモデルの逐次最小自乗法
(k )
(k 0)
ˆ(0), P(1) 0を 与える
確定系のパラメータ推定(3)’
逆行列の補題
( I XY )1 I X ( I YX )1 Y , ( A XY )1 ( I A1 XY )1 A1
( 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 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
[ 弱P E 性]
1
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 )}は観測可能
両辺の信号は等し いので、 ' 両辺をフィ ルタ リ ング’ し た
も のも 等し い. こ こ ではフィ ルタ ーと し て
F (s) :
a
(a 0)
sa
F (s)(my(t ) dy(t ) ky(t ) u(t )) F (s)u(t )
こ こ でsはラ プラ シアンではなく て微分作用素と し て考えている
物理モデルを用いたパラメータ同定(2)
実際の計算では、 例えば
F (s) y(t )の計算には
y f (t ) : ( F (s) y)(t ) ea(t ) ay( )d
t
0
[ea(t ) ay( )]t 0 a ea(t ) ay( )d ( 部分積分の公式)
t
0
ay(t ) aeat y(0) a ea (t ) ay( )d
t
0
a
) y(t )
0
sa
こ の式は加速度にフィ ルタ ーをかけた信号は速度信号よ り
近似的に計算でき る こ と を示し ている 。
ay(t ) a ea(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 )
sa
sa
sa
sa
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)
sa
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 )
sa
sa
sa
a
a
a 1
m cos( y(t )) y(t )
m( sin( y(t )) y(t )) y(t )
sa
sa
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
sa
推定量の性質
真のパラ メ ータ 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は発散する 。 し かし 、 その場合でも 誤差の期待値は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), , 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 )
確率系のパラメータ同定(2)
こ こ で、
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
確率系のパラメータ同定(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 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 )と 独立と なり 、 入力が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 (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 : [n1 , ,0 , n1 , , 0 ] R2n
モデルの構造
パラ メ ータ 将来の入力から 将来の出力が予測可能
過去の入出力
確率系の場合はどう か?
最小分散推定値
(xと は独立ではない)
推定したい
パラメータ x
観測量 y p( y x)
推定量
x̂
推定ルール
g (y)
評価関数 E x x̂
2
を最小にする推定値
条件付き期待値
(
xˆ Ex 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
EEf ( X , Y ) Y
2
E E g(Y ) X Y
2
2
E E g (Y ) EX Y EX Y X Y
E E g xy xy X T g xy xy X Y
E E g xy xy X 2 2g 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 ) EX 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)qiの場合は表現が簡単になる 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 )
出力の予測(1段先予測器)(3)
特に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 )である ので、 予測誤差は白色信号
と な っている 。 ( パラ メ ータ が正し い場合!)
確率系の逐次パラメータ同定(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)
n1
0
u(k n)
e(k 1, )
n1
0
u(k n) e(k 1, )
n1
e(k n, )
0
n1
n1
e(k n)
n1
: 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 dnd 1 p1
1
w(k )
D( p)
一般化最小自乗法
(GLS: Generalized Least Square)
d0 pnd
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 : [n1 , ,0 , n1 , , 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 : [dn1 , , 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)
相良ら:システム同定、コロナ社(1995)
足立:ユーザのためのシステム同定理論、SIC
E(1993)