状態とパラメータの逐次決定

同時推定の考え方 †

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

• 本方法は、逐次解法であるので、オンラインの適応制御に有用である。
• 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期までのデータに基づくｆと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に関して線形である。 そこで、次の２つの問題に分けて、相互依存的な線形問題として解くことを考える。

問題１：状態推定問題：パラメータ所与 †

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

• 状態で微分した上側の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を求めることができる。''

問題２：パラメータ推定問題：状態推定値所与 †

最尤方程式の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)