粒子フィルタによる圧力調整器シミュレーションモデルの
パラメータ同定
Parameter Identification of Pressure Regulator
by Using Particle Filter
石垣 司
1∗樋口 知之
1,2Tsukasa Ishigaki
1and Tomoyuki Higuchi
1,21
情報・システム研究機構 統計数理研究所 予測発見戦略研究センター
1Prediction and Knowledge Discovery Research Center,
The Institute of Statistical Mathematics,
Research Organization of Ingormation and Systems
2
情報・システム研究機構 統計数理研究所 モデリング研究系
2Department of Statistical Modeling,
The Institute of Statistical Mathematics,
Research Organization of Ingormation and Systems
Abstract: This paper describes a parameter identification method of the damping coefficient of pressure regulator with the nonlinear structure by using the particle filter on the nonlinear Gaussian state space model. In the framework of data assimilation, we give a method of the fusion of simulation model and observation data for nonlinear vibration system.
1
はじめに
対象とするシステムの動的特性を正確に再現可能な 物理・数学モデルを構築する作業（モデリング）は、シ ミュレーション物理、逆問題解析、制御系の設計、状 態の推定・予測などモデルベースの技術にとって根底 となる重要なものである。その中でも、直接観測する ことが不可能もしくは困難なモデル内のパラメータに ついてその値を決定する方法は、パラメータ同定問題 として数多くの研究がなされている [1, 2]。 従来、非線形モデルのパラメータ同定は、平衡点付 近での線形近似によるモデル自身の線形化や、各時刻 での状態量の周辺で線形近似を行い拡張 Kalman フィ ルタを適用する方法が主であった。しかしながら、こ れらの手法はモデルが強い非線形性をもつ場合、その 推定誤差は非常に大きくなることが知られている [3]。 近年、非線形な状態推定問題を取り扱うことが可能 な逐次 Bayes 推定の手法であるアンサンブル Kalman フィルタ [4] や粒子フィルタ [5, 6, 7] が発表され、非 線形モデルのための状態推定手法として注目されてい る。また、それらの手法は推定される状態の次元が高 ∗連絡先：情報・システム研究機構 統計数理研究所 住所 〒 106-8569 東京都港区南麻布 4-6-7 E-mail ishigaki@ism.ac.jp 次元の場合でも現実的に適用可能であるため、海洋学 や気象学などの分野において、逐次型データ同化手法 [8, 9, 10, 11, 12]の一つとして研究されている。データ 同化とはシミュレーションモデルと観測データとの融 合を目的とし、非線形モデルにおいて初期条件・境界 条件・モデル内パラメータなどの与え方などに起因す るシミュレーション結果と実現象との乖離を補正する 一つの方法論である。実際に観測・計測されたデータ を用いて各種条件・パラメータの値を修正することに より、より現実に適合するシミュレーションモデルの 構築を目指す。 一方、工学の分野においても、物理モデルによるシ ミュレーションは製品の設計や対象の理解において強 力な道具となることが認識されて久しい。自動車や飛 行機を例に取ると、順問題のシミュレーションモデルを 解くことによる再現実験と風洞実験等による実データ 計測の両者を行うものが常である。しかしながら、そ の多くはシミュレーションと実験の両者を並行して行 い、互いの結果を検証するに留まるものであり、両者 の融合という観点からの検証は乏しい。工学の分野に おいてもシミュレーションモデルと観測データの融合 によるモデル自身の改善や計測点の設計などは非常に 有効な手法であると考えられる。 人工知能学会研究会資料 SIG-DMSM-A703-08 (2/29)本発表では、非線形な物理特性をもつ振動系の一つで ある高圧ガス圧力調整器の物理モデルとある観測デー タを利用し、そのフィードバック機構のダンパ係数を 同定する。その手法として調整器の物理モデルから非 線形ガウス型状態空間モデルを構成し、粒子フィルタ により逐次 Bayes 推定を行う。 工学分野でのモデルと観測データ両者の融合を目 指す研究として、計測融合シミュレーション [13, 14] や ICEME (Interactive Computational-Experimental
Methodologies) [15]などが提案されている。前者は制 御工学におけるオブザーバを利用し、シミュレーション と観測データの差をシミュレーションモデルにフィー ドバックさせている。その目的はシミュレーション結果 を漸近的に観測データに収束させることであり、デー タ同化の目的とは異なる。後者の目的はデータ同化と 同様であり、手法として一般化最小二乗法を採用して いる。しかしながら、両者ともその構成には確率密度 の概念は含まれていないため、本発表の内容とは異な るものである。
2
シミュレーションモデル
2.1
高圧ガス圧力調整器
ガス圧力調整器には様々な型があるが，その原理と 基本機構は同じである．本論では図 1 に示す単段式圧 力調整器を例にとり検討を進める．図中の記号は次節 に示す物理モデルの各変数に対応する。 高圧ガス容器内の圧力は許容最大値で 1.56MPa，許 容最小値で 0.07MPa と高圧である．圧力調整器は，こ れらの高圧ガスを給湯器やコンロなどの家庭用ガス燃 焼機器で使用できる 2.8 ± 0.5kPa まで減圧・調整する機 器である．その機構はダイアフラム (Diaphragm) 大気 圧側の空気室 (Air chamber) と減圧室 (Decompressionchamber)内部間の差圧とバネ (Spring) により，復元力 をバランスさせることにより弁 (Rubber valve) の開度 を調整するものである．減圧室内の圧力が目標値より 高くなると，ダイアフラムが押し上げられ，その変位 がリンクを介してゴム弁とノズル間隔を狭めガスの流 量を絞る．一方，減圧室内の圧力が目標値より低くな ると，ダイアフラムが下がり，ゴム弁とノズル間隔が 広がり，ガス流量が増加する．圧力調整器はこのよう にフィードバック系として構成されている．
3
圧力調整器の物理モデル
圧力調整器の物理モデルは、(1) オリフィスへ流れ込 む管路内流体とベルヌーイの定理 [16]、(2) 理想気体の 状態方程式、(3) ノズル−フラッパ系の流体方程式、(4) ! " !# ! $ ! % & 図 1: 単段式高圧ガス用圧力調整器 質量−ダンパ−バネの二次振動系、を利用してモデル 化することができる。その詳細を以下に示す。3.1
変数とパラメータ
このモデルで使用する変数とパラメータを定義する． Qt(t):LPガス容器から調整器入口への流量 [kg/s]，Q1(t): ノズルから減圧室への流量 [kg/s]，Q2(t):減圧室出口か らゴムホース入口への流量 [kg/s]，Qa(t):空気室と大気 間の流量 [kg/s]，Pt:LPガス容器内の圧力 [Pa]，P1(t): 調整器入口での圧力 [Pa]，Pd(t):減圧室内の圧力 [Pa]，Patm:大気圧 [Pa]，Pa(t):空気室内の圧力 [Pa]，A1:調
整器入口の断面積 [m2]，A3:調整器出口の断面積 [m2]， Aa:通気口の断面積 [m2]，dn:ノズルの直径 [m]，ρa:空 気の密度 [kg/m3]，ρ g:ガスの密度 [kg/m3]，S:ダイアフ ラムの有効受圧面積 [m2]，V 1:調整器入口からノズルま での容積 [m3]，V 2:減圧室の初期容積 [m3]，Va:空気ダ ンパ用空気室の初期容積 [m3]，R:気体定数 [J/(kg·K)]， T :温度 [K]，m:ダイアフラムとばねの質量 [kg]，d:ダイ アフラムとばねのダンピング係数 [N·s/m]，k:ダイアフ ラムとばねのばね定数 [N/m]，x(t):ダイアフラムの変位 [m]，x0:ばねの締め上げ量 [m]，l:ダイアフラムに対し てのゴム弁の変動比，z(t):ノズルとダイアフラムの距離 [m]，g:重力加速度 [m/s2]，ci:各流量係数 (i=1,2,3,a)， cf:単位補正係数 [m3/mol]．
3.2
流量と圧力
流量 Qt，Q2はある管路からオリフィスのある管路 へ流れ込む流体としてモデル化できる．流量 Qt，Q2 はベルヌーイの定理よりそれぞれ式 (1)、式 (2) 式で与 えられる [16]． Qt(t) = c1A1sgn(Pt− P1(t)) × s 2 ρg |Pt− P1(t)| (1)Q2(t) = c3A3sgn(Pd(t)− Patm) × s 2 ρg |Pd(t)− Patm| (2) 管路内圧力 P1はそれぞれの管路に流れ込む流量と流 れ出る流量の差を用いることでモデル化できる．機器 内部の LP ガスは理想気体とみなせこの状態方程式よ り管路内圧力 P1は式 (3) であたえられる． P1(t) = cfRT V1 Z (Qt(t)− Q1(t))dt (3)
3.3
ノズルとゴム弁間の流量と減圧室圧力
ノズルとゴム弁間の流量はノズルフラッパ系とみな すことができ式 (4) のように与えられる. Q1(t) = πc2dnz(t)sgn(P1(t)− Pd(t)) × s 2 ρg |P1(t)− Pd(t)| (4) 減圧室の容積はダイアフラムの変動の影響を受ける．そ の影響を考慮すると，状態方程式より減圧室の圧力は 式 (5) のようになる. Pd(t) = cfRT V2− x(t)S Z (Q1(t)− Q2(t))dt (5)3.4
空気室の流量と圧力
空気室と大気との流量 Qaは上記と同様に式 (6) で 表される． Qa(t) = caAasgn(Patm− Pa(t)) × r 2 ρa|P atm− Pa(t)| (6) 空気室は大気と連通している．加えて，その容積はダ イアフラムの影響を受けるため，その圧力は式 (7) の ようになる． Pa(t) = cfRT Va+ x(t)S Z (Qa(t))dt +Patm Va Va+ x(t)S (7)3.5
ダイアフラムと弁のダイナミクス
ダイアフラムは減圧室の圧力と空気室の圧力との平 衡からその位置が決まる. ダイアフラムの平衡状態の位 置が仕様の範囲内におさまる用に x0の値を定める．ダ イアフラムの振動は質量−ダンパ−バネ系の振動から 誘発されフィードバック制御の不安定性として継続さ れる．x はダイアフラム下方向への変位を正とし，正 の値をとる．ばねの締め付け量と重力を考慮するとダ イアフラムの運動方程式は式 (8) のようになる． md 2x(t) dt2 + d dx(t) dt + kx(t) = S(Pa(t)− Pd(t)) +kx0+ mg (8) 変位 x により，z はリンク構造のレバー比に比例して 変動し，弁の開度を調整する． z(t) = l· x(t) (9) 圧力調整器のパラメータには以下の値を使用した． Pt=0.1713∼1.6613[MPa] (絶対圧)， Patm = 0.1013[M P a] (絶対圧)，ρa = 1.2[kg/m3]， ρg = 2.07[kg/m3]，R=8.31[J/(mol·K)]，T =300[K]， x0=6×10−2[m]，g=9.8[m/s2]，c1,c2,c3,c4,c5=1， cf=22.4×101 −3[m3/mol]，シミュレーションサンプリン グ時間間隔 ∆t=1×10−5[s]とした．上記以外のパラメー タの値は企業情報に含まれるため公表しない。 モデル内のパラメータは実際の圧力調整器とその部 品を測定することにより値を決定することができる。し かし、ダイアフラムとバネのダンピング係数 d を実測 することは、他のパラメータと比べ困難である。よっ て、本発表ではパラメータ d をデータ同化の手法を用 い、その値の同定を行う。4
パラメータ同定手法
本章では、パラメータ同定手法の枠組みについて述 べる。4.1
非線形状態空間モデル
パラメータ同定では、線形状態空間モデルと Kalman フィルタを用い各状態変数やパラメータを推定する手 法や、非線形モデルを一次近似し拡張 Kalman フィル タに同定を行う手法が存在する [17, 18]。しかしなが ら、前章に示したように圧力調整器の物理モデルは強 い非線形性を持つため、これらの手法ではその推定が 困難である。よって、ここでは非線形状態空間モデル の採用が適切である。また、調整器内部で発生する乱 流によるノイズはガウス分布で近似できる。そのため、 ここでは非線形ガウス型状態空間モデルを採用し、調 整器内部振動系のダンピング係数を同定する。 非線形ガウス型状態空間モデルは以下のように表現Q1 Qt Q2 Qa P1 z x Pa Pd Q1 Qt Q2 Qa P1 z x Pa Pd Q1 Qt Q2 Qa P1 z x Pa Pd
t
t-1
t+1
d
d
d
...
...
: Latent variable : Observed variable
図 2: パラメータ同定のためのグラフィカルモデル できる。 xt = ft(xt−1, vt), (10) yt = ht(xt, wt), (11) vt∼ N(0, Qt), wt∼ N(0, Rt) ここで、xtは時刻 t での状態ベクトル、ytは観測時系 列ベクトル、ftと htは非線形関数、Qtと Rtはそれ ぞれシステムノイズと観測ノイズの分散共分散行列、 N (0, V )は平均 0 分散共分散行列 V のガウス分布とす る。 本手法ではパラメータ d も時変の変数とみなし、状 態ベクトルの要素の一つとすることによりその値を推 定する。
4.2
粒子フィルタ
ここでは、粒子フィルタ [7] の概略について述べる。 粒子フィルタは、状態ベクトル x の予測分布とフィ ルタ分布を p(xt|y1:t)≅ 1 N N X i=1 δ(xt− x (i) t|t) (12) p(xt|y1:t−1)≅ 1 N N X i=1 δ(xt− x (i) t|t−1) (13) として、確率密度関数 p(xt|y1:·)から生成された N 個 のサンプル集合{x(i)t|·}N i=1を用い近似する。ただし、こ こで xt|jは p(xt|y1:j)の条件付平均を表す。また、δ(·) はデルタ関数である。この近似を用い、以下のアルゴ リズムにより逐次的にフィルタ粒子を求めることがで きる。 1) i = 1,· · · , N について、初期分布 p0(x)からの初 期粒子を生成する。 2) 各時刻 t = 1,· · · , T について以下のステップを実 行する。 a) 乱数{v(i)t }Ni=1を生成し{x (i) t|t−1} N i=1を x(i)t|t−1= ft(x (i) t−1|t−1, v (i) t )に従い計算する。 b) 尤度関数に基づきそれぞれの粒子の重み wt(i)= ht(yt|x (i) t|t−1)を計算する。 c) b)で計算された重みに比例する割合で {x(i) t|t−1} N i=1からの復元抽出を行い、 フィルタ粒子{x(i)t|t}N i=1を生成する。4.3
観測データ
ここでの観測データは圧力調整器の故障診断のため に計測された圧力変動データ [19] をパラメータ同定の ために用いる。本モデルでは Qaの値を観測しているこ とになる。その計測と計測器の状況を付録に記す。実 際の観測データは電圧値であるため、適宜、その値を 圧力に変換したものを用いる。また、このデータはサ ンプリングタイム 0.1ms で計測されている。よってシ ミュレーションモデルの時間更新の 10 回に 1 度の頻度 でデータが観測される。図 2 に本モデルのグラフィカ ルモデルを示す。0 100 200 300 400 500 600 700 800 900 1000 0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 Assimilation step Da mp ing co eff ici en t d 図 3: ダンピング係数の同定結果
5
パラメータ同定実験
5.1
実験の条件
ここでは粒子数を 100、シミュレーション時間を 0.1 秒とした。初期条件、システムノイズ、観測ノイズ、乱 数の振り方等の各種条件はは本発表時に詳述する。5.2
実験結果
図 3 に同定実験の結果を示す。実験の結果、ダンピ ング係数は約 3 から 4 の間と推定されている。これは ダイアフラムとバネの質量・バネ係数との兼ね合いか ら、減衰する 2 次振動のダンピング係数としては妥当 な値である。また、結果の詳細は本発表時に詳述する。6
まとめ
本発表では、非線形状態空間モデルと粒子フィルタ を用いた非線形構造をもつ圧力調整器の振動ダンピン グ係数の同定について述べた。その結果、未知のダン ピング係数の値は現実的に妥当な範囲の値として推定 された。これにより、未解決の課題は多いが、非線形 モデルのためのデータ同化手法の工学的応用の一つの 枠組みを提供できる可能性が示された。謝辞
高圧ガス圧力調整器の物理モデルの作成にあたり法 政大学工学部渡辺嘉二郎教授には多大なるご指導をい ただきました。厚く御礼申し上げます。A
付録：観測データの計測装置
図 4 に開発したシステムの概要を示す．圧力調整器 を分解できない．そのため図 4 に示す圧力調整器の空 気室に設けられている通気口 (Vent hole) からの空気流 の変化に伴う圧力変動をとらえる．空気室全体を被う 半密閉容器を調整器上部に置き，空気流により変化す る容器内の圧力を容器内に設置された安価で小型な低 周波広帯域マイクロフォンでとらえる．容器を密閉し ない理由は空気室と大気の連通を維持するためである． 実際上はこの容器を調整器に置くだけであるから密閉 にはならない．振動空気流成分は調整器の非線形性の ために波形が歪む．少し広めに 5Hz∼500Hz の範囲で 信号をとらえ，商用電源誘導ノイズを除去するノッチ フィルタ・増幅器を介してコンピュータに取り込む．参考文献
[1] 足立修一：MATLABによる制御のための上級システム 同定、東京電機大学出版局（2004） [2] 片山徹：システム同定入門、朝倉書店（1994）[3] S. Thrun, W. Burgard and D.Fox: Probabilistic
Robotics, MIT Press, (2005)
[4] —:The ensemble Kalman filter: Theoretical formula-tion and practical implementaformula-tion, Ocean Dynamics, vol.53, pp.343-367 (2003)
[5] G. Kitagawa:Monte Carlo filter and smoother for non-Gaussian nonlinear state space model, Journal of Computational and Graphical Statistics, vol.5 pp.1-25 (1996)
[6] A. Doucet, J. de Freitas and N. Gordon: Sequential Monte Carlo Methods in Practice, Springer-Verlag, (2001) [7] 樋口知之：粒子フィルタ、電子情報通信学会誌、Vol.88、 No.12、pp. 989-994、(2005)
Pressure regulator
Condenser
microphone
Semi-closed
vessel
Vent hole
Microcomputer 5Hz～500Hz band-pass filter Notch filter for electrical hum noise[8] G. Evensen: Data Assimilation: The Ensemble Kalman Filter, Springer-Verlag (2006)
[9] 中村和幸,上野玄太,樋口知之,データ同化：その概念と計 算アルゴリズム,統計数理, Vol.53、No.2、pp.211-229、 （2005）
[10] K. Nakamura, T. Higuchi, N. Hirose, Application of particle filter to identification of tsunami simulation model, Proceedings of SCIS & ISIS 2006, pp.1890-1895, (2006)
[11] G. Ueno, T. Higuchi, T. Kagimoto, and N.
Hi-rose: Application of the ensemble Kalman filter
and smoother to a coupled atmosphere-ocean model, Scientific Online Letters on the Atmosphere, Vol.3, pp.5-8 (2007)
[12] S. Nakano, G. Ueno, and T. Higuchi: Merging par-ticle filter for sequential data assimilation, Nonlinear Processes in Geophysics, Vol.14, pp.395-408 (2007)
[13] 早瀬敏幸：計測とシミュレーションの融合による流れの
実現象の再現、日本応用数理学会誌, Vol.16-1, pp.78-84 (2006)
[14] T. Hayase and S. Hayashi: State Estimator of
Flow as an Integrated Computational Method with the Feedback of Online Experimental Measurement, Journal of Fluids Engineering, Transactions of the ASME, Vol.119, pp.814-822 (1997)
[15] J. S. Szmyd, K. Suzuki, Z. S. Kollenda and J. A. C. Humphrey: A Study of Thermo-Fluid Phenom-ena with Uncertainties by making use of Interactive Computational-Experimental Methodology, Vol.35-II, No.4, pp. 599-607 (1992) [16] 増渕，川田：システムのモデリングと非線形制御，コロ ナ社(1996) [17] 西山清：最適フィルタリング、培風館(2001) [18] 村上、堀、登坂、鈴木：有限要素法・境界要素法による 逆問題解析、コロナ社（2002） [19] 石垣、樋口、渡辺：Kullback-Leiblerカーネルによる正 規化周波数スペクトル判別とその圧力調整器劣化診断へ の応用、電子情報通信学会論文誌D、Vol.90、No.10、 pp.2787-2797 (2007)
