2009-12-29[n年前へ]
■「カルマンフィルタ」と「エクセルで解く2次元非定常熱伝導問題」 
正月に、（自分用の）汎用「カルマンフィルタ」ライブラリをRubyとCとExcelで書いてみることにした。たとえば、さまざまなデータ、たとえば、信頼性が低く、誤差の大きなセンサデータや、安定性に欠ける実験データから、現実に近い状態量を推定するツールを作ってみることにした。そして、何か（解析式による）モデル計算や各種シミュレーション計算と比較をしてみたり、それらの計算改善へのフィードバック例を作ってみよう、と考えた。
そこで、扱う題材を考えつつ、実際に上記のようなことを行っている例を探してみた。すると、たとえば、
といったものがある。これらの記事が（下に張り付けた動画でその一端がわかると思うが）実にわかりやすく・面白くて楽しく・役に立ちそうに見える。何というか、つまるところ、魅力を持つに必要な三拍子がすべて備わっている。
先日、「２次元非定常熱伝導問題を解く」エクセル・シート、しかも、そのシートに、センサ機能／フィードバック機能なども付けてみた。そんな素材・材料が揃ってきたこともあるので、まずは、PID制御で（疑似三次元空間における）温度制御を行う例をいくつか作り、その後は、上記記事を参考にしつつ「（誤差を付加した）センサ→カルマンフィルタ→制御量最適化」という例でも作ってみることにしよう。
2009-12-31[n年前へ]
■Rubyで単純なカルマンフィルタを書いてみた
An Introduction to the Kalman Filterの一定出力プラント例に対するカルマンフィルタをPythonで実装した
Cookbook / KalmanFilteringを参考に、Rubyで単純なカルマンフィルタを適当に書いてみました。
動かし方は、
ruby kalman.rb 100というようにでも起動すると、「出力」→「観測」→「予測」のカルマン・フィルタの動作を99回ほど繰り返す、という具合です。出力結果例は、また後ほど示すころにしますが、開始後およそ5回程度からは、「真の出力値」に近い予測を行うことができています。
（最後に配列を使ってグラフ化などをするわけでもないことと、一回前のデータ以外は不必要だということがわかりにくくなってしまいそうだったので）本来は不要な配列を使いたくありませんでした。ただ、変数が（一見）増えたように見えるのもどうかと思い、Python例と同じように配列を使った実装にしてしまいました。
一応、実際に動作する際に「発生すること」「処理すること」の「手順・ループ」を意識しながら、コメント文を書いてみました（正式な用語とは異なるものも多いかもしれませんが）。
ただし、今回の「一定出力」の例では、「状態方程式」「出力方程式」が登場しないので、そういう拡張をすることを考えた場合には、変数名が適切でないような気がします。そういった辺りのことは、また次の機会に整理し直すことにします。さて、これで、今年を終わりにしたいと思います。
というわけで、下記がRubyソースになります。
# 簡単なカルマンフィルタRuby実装例
def gaussian(n) # 正規分布作成用適当関数
sum = 0.0
((n*12).to_i).times{ sum += rand() }
return sum-(12/2*n)
end
operationNum=ARGV[0].to_i # Operation Num
x=[] # 出力の真の値
z=[] # 測定値
xhat=[] # 事後出力推定値
xhatMinus=[] # 事前出力推定値
p=[] # 事後誤差推定値
pMinus=[] # 事前誤差推定値
kG=[] # カルマンゲイン
xVall=55.5 # 出力の真の値
q=1.0 # 出力自身が持つノイズの分散
r=2.0 # 観測における誤差の分散
xhat＜＜20 # 初期出力推定値
p＜＜0 # 初期誤差推定値
puts 'k,x,z,xhat,kG'
1.upto(operationNum) do |k|
# DoProcess
x[k]=xVall+gaussian(q) # 出力値生成
z[k]=x[k]+gaussian(r) # 測定値生成
# Time Updade ("Predict")
xhatMinus[k]=xhat[k-1] # 状態予測を進める
pMinus[k]=p[k-1]+q # 誤差共分散計算を進める
# Measurement Update("Correct")
kG[k] = pMinus[k]/( pMinus[k]+r) # カルマンゲイン算出
xhat[k] =xhatMinus[k]+ #測定値を用い事後出力推定値算出
kG[k]*(z[k]-xhatMinus[k])
p[k] = (1-kG[k])*pMinus[k] # 事後誤差算出
# この回のオペレーションでの結果を表示
puts [k,x[k],z[k],xhat[k],kG[k]].join(',') if k＞1
end
2010-01-01[n年前へ]
■Rubyで書いた単純なカルマンフィルタの出力グラフ例
「Rubyで単純なカルマンフィルタを書いてみた」ので、そのスクリプトの動作させた場合の出力結果をグラフにしてみました(RubyソースはsimpleKalman.rbとしてここに置いておきました)。
Ruby simpleKalman.rb 50 > simpleKalman.csvという風にCSVファイルにして、Excelで開き、結果を示したのが下のグラフです。観測ノイズがある中で（観測値が朱線で示したzです）、水色点線で示したxhat（この単純例では出力量の推定値）が、水色で示したx(この単純例では出力量)をそこそこ推定できているようです。
さて、次は「C言語で整数演算だけを使いこの単純な例を実装するときの苦労をする」か「エクセルで説明用の実装をする」か、はたまた、もう少し面白そうなモデルベース予測あるいはシミュレーション計算予測との組み合わせ（たとえば、アクティブ制御のブラ）実現にでも挑戦してみるようか悩んでいるところです。
」など。