同時推定の考え方

線形離散時間モデルのパラメータと状態を同時に推定する方法を解説する。 尤度方程式を解くとパラメータに関する線型方程式と状態に関する線形方程式が得られる。この方法は2つの方程式を交互に繰り返して推定する。

  • 本方法は、逐次解法であるので、オンラインの適応制御に有用である。
  • 1入力1出力系で示しているが、多入力多出力の高次系にも、拡張できる。

ここで提案した方法は漸近的に有効な推定値を与える。

  • 漸近有効性とは、ある母数θ に対し、漸近正規推定量のうち分散nが最小のものを、漸近有効推定量と呼び、この良好な性質をもつ推定方法である。

線形離散時間モデル

入力utと出力ytの観測データに基づいて、離散時間モデルを逐次推定する方法を紹介する。

線形動的システムのモデルとして、次の最も簡単な1次モデルをまず考えよう。

xt+1=fxt+gut+vt     状態推移式
yt=xt+wt        観測式
  • xtが状態であり、vtはシステムノイズ、wtは観測ノイズである。
  • 未知パラメータは(f,g)で定数である。
  • ノイズは、それぞれ独立な正規確率変数であり、互いに無相関であると仮定する。
    E(vt)=0,E(vtvj)=σv^2・δtj
    E(wt)=0,E(wtwj)=σw^2・δtj
    E(vt,wj)=0,E(x0,
  • 初期状態x0は、正規分布N(μx0,σ0)に従い、ノイズと無相関とする。
     E(x0)=μx0 E((x0-μx0)^2)=σ0

尤度関数

パラメータが所与のもとで、n期までの観測データの確率密度関数を次のように表記し、この確率密度関数を最大とする推定値を求める。パラメータベクトルをΘ=(f,g)'とする。

Ln(Θ)=P(YDn,XDn|Θ)
パラメータベクトルをΘ=(f,g)'とする。
YDnはn期までの出力データ:YDn={yt}t=1,~n
XDnはn期までの入力データ:xDn={xt}t=1,~n

ここで、n期までのデータに基づく推定値を次のように表記する。

f*n,g*n:n期までのデータに基づくfとgの推定値。
Θn*=(f*n,g*n):n期までのデータに基づくパラメータの推定値
x*k|n :n期までのデータに基づく事後的な状態推定値。
  • 状態推定値x*k|nは、k<nのとき内挿値と呼び、k=nのときフィルタリング値と呼ぶ。

簡単な試算

t期の状態変数が、初期状態とt期までの入力とシステムノイズから、正規確率密度関数であらわされることを試算してみよう。

x1=fx0+gu0+v0
x2=fx1+gu1+v1=f^2x0+fgu0+gu1+fv0+v1
x3=fx2+gu2+v2=f^3x0+f^2gu0+fgu1+gu2+f^2v0+fv1+v2
  =f^3x0+(f^2gu0+fgu1+gu2)+(f^2v0+fv1+v2)
.....

より

xt=f^tx0+ Σf^(t-1-i)ui+ Σf^(t-i)vi    t=0,1,2,...
Σ は i=0からt-1までの合計

で表される。 これは、確率変数x0,viの加法和であり、互いに無相関で、それぞれ独立な正規分布であった。 uiは既知の入力系列である。そこで、t期の状態変数が正規確率密度関数で表されることがわかる。

また、観測地も観測方程式より、正規確率密度関数で表されることがわかる。

尤度関数の計算

逐次推定を行うので、尤度関数を書き直して、正規分布の積の形に表せることを示す。 ベイズの定理より、尤度関数の漸化式がえられる。

Ln(Θ)=P(YDn,XDn|Θ)
     =p(yn,xn|Θ,XDn-1,YDn-1)P(XDn-1,YDn-1|Θ)
     =P(yn|xn)P(xn|Θ,xn-1)Ln-1(Θ)

尤度関数のn=0での同時確率は、次式の通り。

L0(Θ)=p(yo,x0|Θ)=P(y0|x0)p(x0)

尤度関数の第1項の計算は、観測式yt=xt+wtより、xnが所与のときのynの確率密度関数は

P(yn|xn)=N(xn,σw^2) n=0,1,2,...

で表される。 尤度関数の第2項の計算は,

P(xn|Θ,xn-1)=N(fxn-1+gun-1,σv^2)

これらを用いて、尤度関数の漸化式から、最終的な尤度関数は、正規分布の積で表され

Ln(Θ)=P(yn|xn)P(xn|Θ,xn-1)Ln-1(Θ)
     =N(xn,σw^2)N(fxn-1+gun-1,σv^2)Ln-1(Θ)

これより、求める尤度関数は以下の正規分布である。

Ln(Θ)=(2πσw)^(-n-1)・σv^(-n)・σ0^(-1)
     EXP{-Σ[(yi-xi)^2/(2σw^2)]-Σ[(xi-fxi-gui)^2/(2σv^2)]-(xo-μx0)^2/(2σ0^2)}
ただし、最初のΣはi=0~nの合計で、次のΣはi=1~nの合計

パラメータと状態の同時推定

上記の尤度関数を最大とする状態xiとパラメータΘを求めればよい。

そこで、尤度関数の対数をとり、これをパラメータと状態に関して一回微分をとりゼロと置く。

  • 初期状態で微分
    (y0-x0)/σw^2-(x0-μx0)/σ0^2 + f(x1-fx0-gu0)/σv^2=0
  • xk(k=1~n-1)で微分
    (yk-xk)/σw^2-(xk-fxk-1-guk-1)/σv^2 + f(xk+1-fxk-guk)/σv^2=0
  • xnで微分
    (yn-xn)/σw^2-(xn-fxn-1-gun-1)/σv^2=0
  • パラメータf で微分
    Σ {xi-1(xi-fxi-1-gui-1)/σv^2}=0  Σはi=1~nの合計
  • パラメータg で微分
    Σ {ui-1(xi-fxi-1-gui-1)/σv^2}=0 Σはi=1~nの合計

上記の尤度方程式を連立して解くことで、最尤推定できる。 しかしながら、全体としては非線形であり収束計算が必要となる。

  • 状態で微分した上側の3式については、x0,x1,...xnに関して線形である。
  • パラメータf,gで微分した下側の2式については、f,gに関して線形である。 そこで、次の2つの問題に分けて、相互依存的な線形問題として解くことを考える。

問題1:状態推定問題:パラメータ所与

パラメータ所与のもとでの状態の推定(平滑:スムージング)問題である。

  • 状態で微分した上側の3式をxtの線形関数として整理してみよう。
  • 第1の式は
     (y0-x0)/σw^2-(x0-μx0)/σ0^2 + f(x1-fx0-gu0)/σv^2=0
    より  box0+b2x1=q0
    ただし
    b0=1+(σw^2/σ0^2)+f^2(σw^2/σv^2)
    b2=-f(σw^2/σv^2)
    q0=y0+μx0(σw^2/σ0^2)+fgu0(σw^2/σv^2)
    となる。
  • 第2の式は  b2xk-1+b1xk+b2xk+1=qk
    ただし
    b2:同じ
    b1=1+(1+f^2)(σw^2/σv^2)
    qk=yk+(uk-1-fuk-1)g(σw^2/σv^2)
  • 第3の式は  b2xn-1+b3xn=qn
    b2:同じ
    b3=1+(1+f^2)(σw^2/σv^2)
    qn=yn+un-1g(σw^2/σv^2)

以上より、''パラメータからb0,b1,b2を与え、q0,qk(k=1~n-1),qnを観測データから計算しておけば 逐次内挿値Xk,k=1~nを求めることができる。''

問題2:パラメータ推定問題:状態推定値所与

最尤方程式の4式と5式を使って、パラメータを求めよう。 この場合内挿値xkが問題1を解いて与えられているものとする。 n個の状態内挿値の内積を定義する。

x'x=Σxi^2 Σはi=0~n-1の合計
x'u=Σxiui Σはi=0~n-1の合計
u'u=Σui^2 Σはi=0~n-1の合計

これから、下記の2つの式の形に、5式,6式が整理できる。

Σ(xi^2)f+Σ(xiui)g=Σ(xi+1xi)
Σ(xiui)f+Σ(ui^2)g=Σ(ui+1xi)

f,gは、上記の2個の線形連立方程式を解けばよい。

逐次推定

すべての過去の状態の内挿値が、各期で必要になることが、難点であろう。 そこで、次の繰り返し計算を提案する。

  • 1.パラメータΘ=(f,g)の初期推定値を適当に与え、問題1で状態の推定を行う。
  • 2.求められた状態推定値から、パラメータの再推定を行う。
  • 3.上記を繰り返して、次の収束判定を行う。
    ||Θ*m-Θ*m-1||<ε εは十分小さい値。(例えば、0.001)

良好な収束性を示したとの報告が、研究論文にある。 パラメータが、十分に収束した時点では、通常のカルマンフィルタ を用いた推定値で、問題1の最適内挿値の代用をしても良いであろう。

参考

  • y.BAR-SHAROM:"Optimal Simulations State Estimation and Parameter Identification in Linear Discrete-Time Systems",IEEE AC-17,no.3,JUNE(1972)

トップ   差分 バックアップ リロード   一覧 単語検索 最終更新   ヘルプ   最終更新のRSS
Last-modified: 2010-09-07 (火) 00:05:00 (3333d)