コラム

第4回 Excelでもう一度

カルマンフィルタで逆解析/第4回担当:里 優(2020.06)


前回説明しました、カルマンフィルタを使った変形問題の逆解析を、ExcelのVBAを使って体験してみましょう。なお、観測は繰り返して行いますが、観測される値はいつも同じである場合を想定します。興味のある方は、ExcelのVisual Basicのエディタに下記のとおり(Sub kalman() から End Sub まで)入力し実行してみてください。


(Excelのメニュー[開発]が表示されていない方は、Alt+F11キーでVisual Basicエディタを表示。または[オプション]-[リボンのユーザー設定]-[リボンのユーザー設定]-[開発]にチェックを付ける→[開発]タブ-[Visual Basic]をクリック))



Sub kalman()
  Cells(3, 4) = "項目"
  Cells(3, 5) = "初期値"
  Cells(4, 4) = "観測値"
  Cells(5, 4) = "前期の状態 (の推定値)"
  Cells(6, 4) = "前期の状態の予測誤差の分散"
  Cells(7, 4) = "状態方程式のノイズの分散"
  Cells(8, 4) = "観測方程式のノイズの分散"
  Cells(9, 4) = "荷重×長さ÷断面積"

  If Cells(4, 5) = "" Then Cells(4, 5) = 1
  u = Cells(4, 5)    '観測値,一定=1
  If Cells(5, 5) = "" Then Cells(5, 5) = 10
  xPre = Cells(5, 5)    '前期の状態(の推定値),初期値=10
  If Cells(6, 5) = "" Then Cells(6, 5) = 10
  pPre = Cells(6, 5)    '前期の状態の予測誤差の分散=10
  If Cells(7, 5) = "" Then Cells(7, 5) = 100
  sigmaW = Cells(7, 5)    '状態方程式のノイズの分散=100
  If Cells(8, 5) = "" Then Cells(8, 5) = 1
  sigmaV = Cells(8, 5)    '観測方程式のノイズの分散=1
  If Cells(9, 5) = "" Then Cells(9, 5) = 50
  cc = Cells(9, 5)    '荷重×長さ÷断面積=50

  Cells(10 + 0, 5) = "繰返し回数"
  Cells(10 + 0, 6) = "ヤング率"
  Cells(10 + 0, 7) = "予測誤差の分散"
  For i = 1 To 20
    '状態の予測(予測値は,前回の値と同じ)
    xForecast = xPre
    
    'A=(∂f/∂x)k=k-1
    aa = -cc / xForecast ^ 2
    
    '状態の予測誤差の分散
    pForecast = pPre + sigmaW

    'カルマンゲイン
    kGain = pForecast * aa / (pForecast * aa * aa + sigmaV)

    'カルマンゲインを使って補正された状態
    xFiltered = xForecast + kGain * (u - cc / xForecast)

    '補正された状態の予測誤差の分散
    pFiltered = (1 - kGain * aa) * pForecast

    '値の格納
    Cells(10 + i, 5) = i
    Cells(10 + i, 6) = xFiltered
    Cells(10 + i, 7) = pFiltered

    '繰り返しのため値を更新
    xPre = xFiltered
    pPre = pFiltered
  Next

End Sub


結果は下図のようになり、観測値をもとにカルマンフィルタで補正を繰り返すと、「状態」(ヤング率)の推定値が正解の50に収束していきます。収束の速さはノイズの分散で左右されるので、分散の値を変えて試してみてください。


次回は、有限要素解析にカルマンフィルタを組み合わせます。


グラフ

※上記のExcelサンプルファイルはこちらからダウンロードすることができます