固有直交分解と統計的推定法を応用した 高サンプリングレート

高サンプリングレート流速プローブと低サンプリングレート

計測データを用いた 高サンプリングレート流れ場オフライン推定

コードの開発

青木 良尚^*1

Development of the High Sampling Rate Flow Field Offline Estimation C++ Code that Applies Proper Orthogonal Decomposition and Stochastic Estimation using High

Sampling Rate Velocity Probe Data and Low Sampling Rate PIV Data*

Yoshihisa AOKI^*1

Abstract

Three step method is the offline estimation method of the flow field using control theory. C++ code of the three- step method which doesn’t depend on any operating systems and applications was developed so I open this code to public. Three-step method is offline estimation method for time-resolved velocity fields using time-resolved point measurements and non-time-resolved particle image velocimetry data. This code was validated with the data in Florida State University and I confirmed the result agreed with the result from MATLAB code.

Key Words: Proper Orthogonal Decomposition, System Identification, Fluid Dynamics, Kalman Smoother

概要

スリーステップ法と呼ばれる，制御理論を応用した流れ場のオフライン推定法がある．オペレーティングシステムやア プリケーションに依存せずに，スリーステップ法による流れ場のオフライン推定が可能となるように，スリーステップ法 のC++コードを開発したので，これを公開する．スリーステップ法は，固有直交分解と統計的推定法を応用して，高サン プリングレート流速プローブと低サンプリングレートPIV計測データから，高サンプリングレート流れ場の離散線形状態 方程式を同定し，このシステムに固定区間カルマンスムーザを適用して高サンプリングレート流れ場をオフライン推定す る方法である．このコードは，フロリダ州立大学が所有するデータを用いて検証され，MATLABで作成されたコードと同 じ推定結果が得られることを確認した．

１．はじめに¹⁾

流体制御は，流れの剥離や境界層遷移などの流体現象を意図的にコントロールしようとする技術であ る．航空機の高速化と航続距離の延長につながる巡航抵抗低減を目的とした境界層乱流遷移遅延制御の 研究の歴史は古く，１９３０年代からNACAにおいて受動制御（自然層流翼に代表される）と能動制御

（境界層吸い込み等を用いた層流境界層制御に代表される）の両面から研究が行われた．エンジン性能 の急速な進歩と，境界層乱流遷移遅延制御の実用化への技術的な困難さからこの研究は下火となったが，

* 平成25年12月20日受付（Received 20 December, 2013）

*1 航空本部 風洞技術開発センター (Wind Tunnel Technology Center, Institute of Aeronautical Technology)

オイルショックによる燃料費高騰をきっかけに，巡航抵抗低減が注目され再び研究が活発化した．自然 層流翼は，工作精度や保守に課題があった為，第二次世界大戦期の軍用機を除いて量産機に採用されな かったが，本田技研工業が量産機への採用に適する自然層流翼を開発し，ビジネスジェット機に採用し た²⁾．現在，この機体はFAAの型式証明取得を目指した試験が行われており，実用化が近い．剥離制御 としての能動的境界層制御は，日本においてもUS-2に採用され，既に実用化されている．一方で，能動 的な境界層乱流遷移遅延制御は，エネルギー収支，アクチュエータの複雑な構造や重く複雑な制御シス テムなど課題が多く，実用化されていない．こうした中，1990年代後半には，能動的な境界層乱流遷移 遅延制御よりも性能が劣るが課題を改善可能な，自然層流翼と能動制御を組み合わせるハイブリッド制 御がNASAによって提案された．ハイブリッド制御は，現在，最も実用化に近い能動的境界層制御であ る．

これまで，能動的な流体制御は流体の安定性理論の発達とともに発達していった．これに加えて，近 年のコンピュータの発達に伴い複雑な制御が可能になったことで，より高い制御効果を得るため，能動 的な流体制御に対して制御理論の導入が進められている．流体制御への制御理論の応用を目的とした研 究の中で，固有直交分解を応用した非定常流れ場推定法であるスリーステップ法³⁾の，MATLAB等のア プリケーションに依存しない C++コードの開発を行ったので，これを公開する．スリーステップ法は，

二次元カルマン渦列の実験から計測されたPIVによる低サンプリングレートの流れ場データと，ある位 置に設置されたプローブによる高サンプリングレートの非定常流速データを用いて，高サンプリングレ ートの流れ場を表現する離散線形状態方程式モデルを同定し，このシステムに対してカルマンスムーザ を適用して高サンプリングレートの流れ場情報を推定する方法である．

２．スリーステップ法とその周辺の数学的準備 ２．１ 固有直交分解（Proper Orthogonal Decomposition, POD）の基礎⁴⁾

PODは，データのコンパクトな表現を得るための，多変数の統計的手法である．この方法は，高次元 データの低次元化や，特徴抽出に用いられる．PODの考え方は，多数の相関のある変数を，元々の値に 対する誤差を出来るだけ小さくしながら，より少数の相関の無い変数に減らすことである．これは，デ ータの共分散行列の固有値を基礎とする直交変換によって，データをこの固有値に対応する固有ベクト ルの張る部分空間に展開する．この変換は，変数間の相関を無くすとともに，分散を最大化する．

PODの最も大きな特徴は，元々のデータと低次元化された線形表現の二乗平均誤差を最小化すること である．POD は，非線形問題にも適用されることがあるが，POD から得られるのは，データを表現す る仮想空間において最適に近似される線形多様体である．PODが線形であることは，線形作用素が利用 できるという利点があるが，データが非線形多様体に属する場合，これは主要な制限にもなる．

PODの数学的な表現について述べる．領域Ωに属するデータをθ

( )

( )

µ ^{と時間変動成分}ϑ

( )

( )

( ) ( )

θ = + （１）

ある時間t_kにおける瞬時場を，スナップショットθ^k

( )

(

)

( )

( )

けることと等価である．

(

)

Maximize ^k （２）

(

)

∫

( ) ( )

⋅ ^{：平均作用素}

(

)

=

⋅ ：ノルム

⋅ ^：絶対値

式（２）は，変動成分ϑ^を基底ϕに沿って展開する場合，他の基底関数に沿って展開するよりも大き い平均エネルギーを含むことを意味する．制約条件 ϕ ² =1は，解を一意とするために課される．この制 約条件下における式（２）の解は，ラグランジュの未定乗数法で計算できる．

[ ] ( )

(

)

ϑ （３）

式（３）は，微分値が０の時に最大値を取るので，式（２）の条件は，次の積分固有値問題と等価と なる．

( ) ( ) ( )

( )

k ϑ ϕ λϕ

ϑ ′ ′ ′=

∫

( ) ( )

k ϑ ′

ϑ ^{：平均自己相関関数}

式（２）の最適化問題の解は，PODモードや固有直交モード（Proper Orthogonal Modes, POMs）と 呼ばれる，式（４）の積分方程式の直交固有関数ϕ_i

( )

(

)

固有直交値（Proper Orthogonal Values, POVs）と呼ぶ．POMは，変動成分ϑ

( )

( ) ∑

( ) ( )

=

=

1

,

i ai t i x

t

x ϕ

ϑ （５）

( )

( ( ) ( )

)

a = ϑ _{, ,}ϕ_i

( ) ( )

i t a t

a =δ λ

最も大きなPOVに対するPODモードは，スナップショット全体を特徴付ける最適なベクトルとなる．

データに含まれるエネルギーε ^は，POVの和として定義される（ε ⁼

∑

∑

２．２ PODスナップショット法

POD スナップショット法は，離散化された POD である．ある時間t_m^におけるM^{個のスナップショ}

ットの集合を考える．

( )

( ) (

)

u ^m = , _m , =1,2,2, （６）

平均値を下式とする．

( ) ∑

( )

=

= ^M

m

m x

M u x u

0 1 1 （７）

データは，平均値と変動成分に分解できる．

(

)

( )

(

)

u , = ₀ + ′ , （８）

変動成分の自己相関の離散形式より，POMφ_i

( )

( )

( )

Cφ_i =λ_iφ_i （９）

( ) ( )

(

)

C= M1 ′ , ′ , ^{：離散自己相関行列}

この直交基底を用いて，変動成分は下式に分解される．

( ) ∑ ( ) ( )

=

′ = ^M

i i m i

m a t x

t x u

1

, φ （１０）

( )

(

(

) ( )

)

a_{i m} = ′ _{, m} φ_i

( ) ( )

i t a t

a =δ λ

次に，重み係数を導入したスナップショット法について述べる．

データ行列X ^{を下式とする．}

( ) ( ) ( )

[

]

X = ′ , ₀ ′ , ₁  ′ , （１１）

データ行列X ^{の内積の重み行列を}Wとすると，自己相関行列C^{は下式となる．}

X W X

C= ^T （１２）

こ こ で ，例え ば 二 次元流 れ 場 データ に 対 して， 重 み 行列W を 各 グ リッド の 面 積diag

(

n

)

ndy ,dx

 ^{とすると，}W のベクトルノルムは，運動エネルギーの積分量と解釈できる．

( )

(

( )

)

( )

(

( )

( )

)

u^′ , _j ²2 ^{= ′} , _j ^{T ′} , _j ⁼

∫∫

この時，POMφ_j

( )

( )

直交行列，Σ：対角行列）とすると，下式となる．

( ) ( )

[

]

=

Φ φ x  φ_m x XV （１４）

( ) ( )

( ) ( )

T m

n n

m

V t

a t

a

t a t

a

a ²¹

1

1 1

1

Σ

=

















=   

j i T i

j Wφ δ

φ =

２．３ PODと統計的推定法による低次元化推定（Linear Stochastic Estimation POD, LSE-POD）

まず初めに一般的な統計的推定法について述べる．統計的推定値を，統計学の知識を用いた条件付き 平均の近似値の期待値とする．事象P^{の下で発生する事象}Sの条件付き期待値は，下式の確率p^と等価

である．

( ) ( )

∫

=

= ds

p f

p s s f p P S

P

SP ,

| （１５）

⋅ ：期待値

fSP^：事象S^と事象P^{の結合確率密度関数} fP^：事象P^{単独の確率密度関数}

式（１５）の近似解法として，テイラー級数展開を用いた条件付き平均の統計的推定法が提案されて いる．

+ +

=

=

= _{i j j ij j ijk j k}

i S P p A p B p p

sˆ | （１６）

sˆi：条件付き平均の推定値

線形統計的推定（Linear Stochastic Estimation, LSE）では，式（１６）の線形項のみを使用する．

スリーステップ法では，非定常圧力センサの計測値p_j

( )

( )

流れ場のLSEによる推定値は，下式となる．

( )

( )

uˆ_i = _{ij j} （１７）

式（１７）の係数A_ij^{は，流れ場}u_i

( )

( )

式（１７）は，複数の非定常圧力センサを想定しているが，１つの非定常圧力センサの時系列計測値 を用いる事ができる．これを，時間遅れ線形統計的推定（delay time-LSE）と呼ぶ．

( )

(

)

i t A pt

uˆ = −t （１８）

典型的なPIVによる流れ場データは，多大な計測点数を持つ．ここで，全ての流れ場u_i

( )

PODを用いて低次元化した流れ場の POD係数a_i

( )

( )

( )

aˆ_i = _{ij j} （１９）

これもまた，1つの非定常圧力センサの時系列計測値を用いることが出来る．これを，時間遅れ低次元

化推定（delay time-LSE-POD, delay time-mLSE）と呼ぶ．

( )

(

)

i t A pt

aˆ = −t （２０）

２．４ システムの入力が無い離散線形システムのカルマンスムーザ⁵⁾

スリーステップ法では，システムの入力が無い場合の離散カルマンスムーザを用いる．まず，離散シ ステムを下式とする．

k k

k Fx d

x = ₋₁+ （２１）

k k k

k = H x +n

η

xk：離散システムの状態変数

ηk：離散システムの観測値

dk^{：プロセスノイズ}

nk^{：センサノイズ}

Hk：観測行列

こ こ で ， プ ロ セ ス ノ イ ズ と セ ン サ ノ イ ズ は 白 色 ノ イ ズ （

[ ]

id Q

d

E = δ− _，

[ ]

in R

n

E = δ− _，

[ ]

E _{）とする．}

まず，カルマンフィルタによって推定を行う．初期値を下式とする．

[ ]

0

ˆ , E x

x⁺f = （２２）

( )( )

[

]

f E x x x x

P⁺_,0 = ₀ − ˆ₀⁺ ₀ − ˆ₀⁺

ここで，Pは推定誤差の共分散である．次に，式（２２）を初期値として，下式を用いて状態推定を 行う．

1 1

,

, + −

−_k = _{f k}− ^T + _k

f FP F Q

P （２３）

(

)

, ,

− −

− +

= _{f k} ^T_{k k f k k}^T _k

k

f P H H P H R

K

+ −

−_, = ˆ _{, 1} ˆ_{f k} Fx_{f k} x

(

)

−

+_{f k} = x_fk +K_{f k k} −H_kx_{f k} xˆ , ˆ , , η ˆ ,

( )

+_k = − _{f k k f k}

f I K H P

P , , ,

次に，カルマンフィルタによる推定値を用いて，固定区間のカルマンスムーザによる推定を行う．ま ず，初期値を下式とする．

= +_{f Nt} Nt

s x

xˆ , ˆ , （２４）

= _f+_Nt Nt

s P

P, ,

式（２４）を初期値として，下式を用いて状態推定を行う．

(

)

1 ,

− − +

−

+ = _{f k}

k

f P

L （２５）

− +

= +_{, , 1}

, T f k

k f k

s P F L

K

(

)

s k f k

s P K P P K

P, , , −, 1 , +1 ,

+ − + −

=

(

)

+ + −

= _{, , , 1 , 1}

, ˆ ˆ ˆ

ˆ_sk x_{f k} K_sk x_sk x_{f k} x

以上より，入力が無い離散システムにおける，カルマンスムーザによる状態推定値が求まる．

２．５ スリーステップ法

スリーステップ法は，低サンプリングレートのPIV 計測により得られた流れ場データと，ある位置に 設置されたプローブによる高サンプリングレートの非定常定常流速データから，高サンプリングレート の流れ場を推定する方法である．スリーステップ法は，以下の3段階の手順からなる．

1) 非定常流れ場のPOD係数の統計的推定

delay time mLSEを用いて，低サンプリングレートの流れ場データと高サンプリングレートのある位

置の非定常流速データから，高サンプリングレートの非定常流れ場POD係数の統計的推定値を求める．

2) 非定常流れ場の離散線形システムの同定

LSE を用いて，非定常流れ場 POD 係数の統計的推定値から，流れ場の構造を考慮した，流れ場のダ イナミクスを表現する離散線形システムを同定する．

3) カルマンスムーザによるPOD係数の推定

流れ場の離散線形システムと，低サンプリングレートの流れ場データに対してカルマンスムーザを適 用し，高サンプリングレートの非定常流れ場の推定値を求める．

1段階目の非定常流れ場の初期推定の後に，カルマンスムーザによる推定を行うのは，初期推定値に含 まれるセンサノイズの影響を低減するためである．

３．スリーステップ法によるカルマン渦列の非定常流れ場推定

スリーステップ法を用いて，図１の後流の非定常流れ場推定を行う．ここで，一様流流速は 4.2m/s， 低サンプリングレートのPIVデータのサンプリングレートは32Hz，計測点数は2317点，高サンプリン グレートの流速データのサンプリングレートは800Hzである．

図１ 風洞試験概要図

まず最初のステップとして，低サンプリングレートのPIVによる流れ場データと，ある位置に設置さ れたプローブによる高サンプリングレートの非定常流速データを用いて，非定常流れ場のPOD係数の統 計的推定を行う．統計的推定値を求める準備として，PODを用いてPIVによる流れ場データの低次元化 を行う．ここでは，PODの次数を７次とした．この場合，図２の通り得られたPOD係数a_i

( ) (

)

によって推定される流れ場に含まれる運動エネルギーは全体の 89.4%であるが，カルマン渦列を再現す るためには十分である．

図２ 使用するPODモード数に対する含まれるエネルギーの割合

この低サンプリングレートの流れ場の POD 係数と高サンプリングレートのある点の非定常流速デー

タv

( )

( )

を求める．

( ) ∑ ( )

−

=

+

= ^t

j t ij

i n A v n j

aˆ （２５）

t ：時間遅れ

式（２５）の時間遅れは，図３のように推定されたPOD係数と計測されたPOD係数から計算可能な，

流れ場の運動エネルギーの推定誤差を最小化する時間遅れが存在するので，この最適値を用いる．ここ では，時間遅れt を5とした．

図３ delay time-mLSEの時間遅れに対する流れ場のエネルギー誤差

次のステップとして，式（２５）から得られた非定常流れ場POD係数の統計的推定値から，流れ場の構 造を考慮した，高サンプリングレートの流れ場のダイナミクスを表現する離散線形システムを同定する．

カルマン渦列の場合，事前知識として，1番目と2番目のPOD係数で表現される周期的振動モードと，

その他のモードから表現される事がわかっているので，LSE を用いて，下式のように，振動モードの部 分システムと，その他のモードの部分システムを分けて，それぞれ同定する．

( ) ( ) ( )

( )





−

 −



 



=



 





ˆ 1 ˆ 1 0

0 ˆ

ˆ

n a

n a F F

n a

n a

LSE OSC LSE OSC

LSE

OSC （２６）

( ) ( ) ( )



= n a

n n a

a_OSC

2 1

ˆ ˆ ˆ

( )

( ) ( ) ( ) ( )















= n a

n a

n a

n a n a_LSE

6 5 4 3

ˆ ˆ ˆ ˆ ˆ

次に，振動モードの部分システムF^OSCが不安定にならないように，F^OSCの固有値の絶対値が０．９ ９９となるようにF^OSCを定数倍する．以上により，カルマン渦列の高サンプリングレート流れ場の離散 線形システムが同定できた．

最後のステップとして，１段階目で得られた低サンプリングレートの流れ場データと高サンプリング レートのPOD係数の統計的推定値を観測値として，式（２６）の離散線形システムに式（２３），（２４）

のカルマンスムーザを適用し，高サンプリングレートの非定常流れ場の推定値を求める．観測値の共分 散の初期値P_f⁺_,0^{とノイズの共分散（}Q_k^，R_k）を下式とした．

I

Pf⁺_,0 =5 （２７）

(

)

diag Q_k =





×

= × ₋ I R_{k 10}⁴ I

10 1

10 2

I^：サイズ5×5の単位行列

以上により，低サンプリングレートの流れ場データと，高サンプリングレートのある点の流速から，

高サンプリングレートの流れ場を推定できた．

非定常流速データのノイズを下式とした時，PODの統計的推定値とカルマンスムーザによる推定値に おける流れ場の運動エネルギー誤差の割合は，図４となる．全てのノイズの値でカルマンスムーザによ る推定値の誤差は統計的推定値の誤差を下回る．従って，カルマンスムーザによる推定はノイズに対し てロバストであることがわかる．

2



 





′

= ′

RMS RMS

v

γ n （２８）

nRMS′ ：非定常流速データのノイズのRMS vRMS′ ：非定常流速データの変動成分のRMS

図４ 高サンプリングレート流速プローブのノイズに対する推定POD係数のエネルギー誤差

謝辞

この開発は，２０１２年１１月１日から２０１３年１０月３１日の宇宙航空研究開発機構の長期在外 研修の一部として，アメリカ合衆国フロリダ州立大学 Florida Center for Advanced Aero-Propulsion

（FCAAP）で実施した．フロリダ州立大学教授及び FCAAP Associate Director である，Dr. Louis

Cattafesta には，この在外研修において受入ホストを務めて頂いた．研修を進めるに当たり，航空本部

事業推進部の久野克子氏，竹本勇介氏には事務手続きについてご助力を頂いた．航空本部風洞技術開発 センター計画管理室の内藤朋子氏には，事務手続きの代理を務めて頂いた．航空本部風洞技術開発セン ターの渡辺重也センター長（現，航空本部チーフエンジニア），浜本滋計画管理チーフマネージャ（現，

センター長），満尾和徳超音速風洞/フラッタ風洞セクションリーダ（現，研究計画マネージャ事務代理）

には，研修計画についてご助言を頂いた．また，超音速風洞/フラッタ風洞セクション及び宇宙輸送ミッ ション本部宇宙輸送系システム技術研究開発センターの調布在勤の方々には，不在の間の仕事を引き継 いで頂いた．末筆ながら，ここに感謝の意を表する．

参考文献

1) Ronald D. Joslin, “Overview of Laminar Flow Control”, NASA/TP-1998-208705 2) Fujino. M, “Design and Development of the HondaJet”, AIAA-2003-2530

3) Jonathan H. Tu, Clarence W. Rowley, John Griffin, Louis Cattafesta, Adam Hart, Lawrence S.

Ukeiley, “Integration of non-time-resolved PIV and time-resolved velocity point sensors for dynamic estimation of time-resolved velocity fields”, AIAA-2012-033

4) GAETAN KERSCHEN, JEAN-CLAUDE GOLINVAL, ALEXANDER F. VAKAKIS and

LAWRENCE A, BERGMAN, “The Method of Proper Orthogonal Decomposition for Dynamical Characterization and Order Reduction of Mechanical Systems : An Overview, Nonlinear Dynamics”, 2005, 41, 147-169

5) 片山 徹, 「新版応用カルマンフィルタ」, 朝倉書店, 2000

付録. スリーステップ法によるカルマン渦列のカルマンスムーザによるオフライン推定C++コード スリーステップ法による，高サンプリングレート非定常流速プローブと，低サンプリングレート PIV デー タを用いた，高サンプリングレートのカルマン渦列の流れ場を推定する C++コードを記載する．C++コード

は，”main.cpp”，”matclass.cpp”，”matclass.hpp”の３つのファイルから成る．入力データは，第１引数に高

サンプリングレート非定常流速プローブデータファイル名（文字列），第２引数に非定常流速プローブのサン プリングレート(Hz，整数)，第３引数に低サンプリングレートPIVデータファイル名（文字列），第４引数に PIVデータのサンプリングレート（Hz，整数），第５引数にスナップショットに含まれるPIVデータ数（vx+vy， 整数），第６引数に計測時間（秒で指定，整数），第７引数にdelay time-mLSEの時間遅れ（整数）とする．

それぞれのデータは，下図のように同期して計測する必要が有る．

付図１ 流速プローブデータ及びPIVデータの同期計測

データフォーマットはテキスト形式である．流速プローブデータは，１列の行列とする．また，PIVデータ は，行方向を位置，列方向を時間として，カンマ区切りのテキスト形式とする．但し，流速プローブ及び各位 置における PIV データの平均値がゼロとなるように，それぞれのデータから平均値を引いておくこと．コー ドのデータ読み込み部分を修正すれば，どんなデータフォーマットにも対応可能である．推定された流れ場は，

ファイル名”result.csv”のファイルに，PIV と同一のフォーマットで保存される．推定された流れ場のサンプ リングレートは，流速プローブのサンプリングレートと同一となる．その他，一般的な注意点としては，POD スナップショットの重み係数を1に，使用するPODモード数を7に固定しているほか，コードの開発では実 行速度やメモリ使用量，エラー処理については深くは考慮していない．

例えば，GNU Compiler Collectionのg++で最も単純にプログラムを作成するためには，下記のコマンドを

打てばよい．

g++ main.cpp matclass.cpp matclass.hpp

OSがWindowsの場合，実行可能形式ファイル名が”a.exe”，流速プローブデータファイル名”probe.csv”，

流速プローブのサンプリングレート800Hz，PIVデータファイル名”piv.csv”，PIVサンプリングレート32Hz， データ数4634，計測時間1秒，delay time-mLSEの時間遅れt ^を5，データファイルがexeファイルと同一 のディレクトリにある場合，プログラムを実行するためには下記のコマンドを打てばよい．

a probe.csv 800 piv.csv 32 4634 1 5

// “main.cpp”

// This program is free software: you can redistribute it and/or modify // it under the terms of the GNU General Public License as published by // the Free Software Foundation, either version 3 of the License, or // (at your option) any later version.

//

// This program is distributed in the hope that it will be useful,

// but WITHOUT ANY WARRANTY; without even the implied warranty of // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the // GNU General Public License for more details.

//

// You should have received a copy of the GNU General Public License // along with this program. If not, see <http://www.gnu.org/licenses/>.

#include <iostream>

#include <cstdio>

#include <complex>

#include <cstdlib>

#include "matclass.hpp"

//SNUM is the maximum number of the characters

#define SNUM 20000

int main(int argc, char* argv[]) {

int i,j,k,m,n,sp,num,pnum,delay,rep,kk,sec,fpp,fpiv,dnum,dn;

double dx,dy,tempd,*vec1,*vec2,*pdata;

complex<double> *eigen_val;

matclass *drms,*amat,*bmat,*cmat,*fit_podc,*fit_coef,*dmodel,*kf_podc,*kfm_podc;

matclass **pm,**pp,*kg,*qmat,*rmat0,*rmat1,*umat,*ks_podc,*ps;

FILE *fp;

char str[SNUM],*tok;

printf("Probe data file name:%s¥n",argv[1]);

printf("Probe sampling rate:%s[Hz]¥n",argv[2]);

printf("PIV data file name:%s¥n",argv[3]);

printf("PIV sampling rate:%s[Hz]¥n",argv[4]);

printf("The Number of the data in one PIV snapshot:%s¥n",argv[5]);

printf("Measurement time:%s[sec]¥n",argv[6]);

printf("Delay time:%s¥n",argv[7]);

dx = dy = 1.;

delay = atoi(argv[7]);//delay time pnum = 7;// the number of POD sec = atoi(argv[6]);

//if this routine is modified, this code can be used to any import file format data.

//drms(i,j) = u_i(j) puts("Import Data");

fpp = atoi(argv[2]);

pdata = new double[sec*fpp];

fp = fopen(argv[1],"r");

if(fp == NULL) {

printf("Can't open %s¥n",argv[1]);

return 1;

}

for(i=0;i<sec*fpp;i++) {

memset(str,0,SNUM);

fgets(str,SNUM,fp);

pdata[i] = atof(str);

}

fclose(fp);

dnum = atoi(argv[5]);

fpiv = atoi(argv[4]);

dn = fpp / fpiv;

num = sec*fpiv-1;

drms = new matclass(dnum,num+1);

fp = fopen(argv[3],"r");

if(fp == NULL) {

printf("Can't open %s¥n",argv[3]);

return 1;

}

for(i=0;i<dnum;i++) {

memset(str,0,SNUM);

fgets(str,SNUM,fp);

tok = strtok(str,",¥n");

(*drms)(i,0) = atof(tok);

for(j=1;j<num+1;j++) {

tok = strtok(NULL,",¥n");

(*drms)(i,j) = atof(tok);

} }

fclose(fp);

drms->set_weight(dx*dy);

drms->pod(pnum);

//data import finishes.

////////////////////////

//dynamic model

puts("Construct Dynamic Model");

amat = new matclass(num,delay*2+1);

bmat = new matclass(num,pnum);

cmat = new matclass(delay*2+1,pnum);

for(j=0;j<num;j++) {

for(k=0;k<delay*2+1;k++)

(*amat)(j,k) = pdata[(j+1)*dn-delay+k];

}

for(j=0;j<num;j++) {

for(k=0;k<pnum;k++)

(*bmat)(j,k) = drms->podc(k,j+1);

}

tempd = amat->lsq(bmat,cmat);

cmat->transpose();

fit_coef = new matclass(*cmat);

delete amat;

delete bmat;

delete cmat;

fit_podc = new matclass(pnum,sec*fpp);

for(k=delay;k<sec*fpp-delay;k++) {

for(j=0;j<pnum;j++) {

for(int kk=-delay;kk<=delay;kk++)

(*fit_podc)(j,k) += (*fit_coef)(j,kk+delay) * pdata[k+kk];

} }

amat = new matclass(sec*fpp-2*delay,pnum-2);

bmat = new matclass(sec*fpp-2*delay,pnum-2);

cmat = new matclass(pnum-2,pnum-2);

for(k=delay;k<sec*fpp-delay-1;k++) {

for(j=2;j<pnum;j++) {

(*amat)(k-delay,j-2) = (*fit_podc)(j,k);

(*bmat)(k-delay,j-2) = (*fit_podc)(j,k+1);

} }

amat->lsq(bmat,cmat);

cmat->transpose();

dmodel = new matclass(pnum,pnum);

for(i=2;i<pnum;i++) {

for(j=2;j<pnum;j++)

(*dmodel)(i,j) = (*cmat)(i-2,j-2);

}

delete amat;

delete bmat;

delete cmat;

amat = new matclass(sec*fpp-2*delay,2);

bmat = new matclass(sec*fpp-2*delay,2);

cmat = new matclass(2,2);

for(k=delay;k<sec*fpp-delay-1;k++) {

for(j=0;j<2;j++) {

(*amat)(k-delay,j) = (*fit_podc)(j,k);

(*bmat)(k-delay,j) = (*fit_podc)(j,k+1);

} }

amat->lsq(bmat,cmat);

cmat->transpose();

eigen_val = new complex<double>[2];

tempd = ((*cmat)(0,0)-(*cmat)(1,1))*((*cmat)(0,0)-(*cmat)(1,1))+4.*(*cmat)(0,1)*(*cmat)(1,0);

if(tempd > 0.) {

eigen_val[0] = complex<double>(((*cmat)(0,0)+(*cmat)(1,1))/2.+sqrt(tempd)/2.,0.);

eigen_val[1] = complex<double>(((*cmat)(0,0)+(*cmat)(1,1))/2.-sqrt(tempd)/2.,0.);

} else {

eigen_val[0] = complex<double>(((*cmat)(0,0)+(*cmat)(1,1))/2.,sqrt(-tempd)/2.);

eigen_val[1] = complex<double>(((*cmat)(0,0)+(*cmat)(1,1))/2.,-sqrt(-tempd)/2.);

}

tempd = abs(eigen_val[0]);

delete amat;

delete bmat;

for(i=0;i<2;i++) {

for(j=0;j<2;j++)

(*dmodel)(i,j) = (*cmat)(i,j) / tempd * 0.999;

}

delete cmat;

delete[] eigen_val;

//kalman filter

puts("Kalman Filter Estimation");

kf_podc = new matclass(pnum,sec*fpp);

for(i=0;i<pnum;i++)

(*kf_podc)(i,delay) = (*fit_podc)(i,delay);

kfm_podc = new matclass(pnum,sec*fpp);

pp = new matclass*[sec*fpp];

pm = new matclass*[sec*fpp];

for(i=0;i<sec*fpp;i++) {

pp[i] = new matclass(pnum,pnum);

pm[i] = new matclass(pnum,pnum);

}

for(i=0;i<pnum;i++)

(*pp[delay])(i,i) = 5.;

amat = new matclass(pnum,pnum);

bmat = new matclass(pnum,pnum);

cmat = new matclass(pnum,pnum);

qmat = new matclass(pnum,pnum);

(*qmat)(0,0) = 1.;

(*qmat)(1,1) = 1.;

(*qmat)(2,2) = 0.004;

(*qmat)(3,3) = 0.5;

(*qmat)(4,4) = 0.5;

(*qmat)(5,5) = 0.5;

(*qmat)(6,6) = 0.5;

rmat0 = new matclass(pnum,pnum);

rmat1 = new matclass(pnum,pnum);

for(i=0;i<pnum;i++) {

(*rmat0)(i,i) = 20000.;

(*rmat1)(i,i) = 1e-10;

}

kg = new matclass(pnum,pnum);

umat = new matclass(pnum,pnum);

for(i=0;i<pnum;i++) (*umat)(i,i) = 1.;

for(i=delay+1;i<sec*fpp-delay;i++) {

amat->matrix_mult(*dmodel,*pp[i-1]);

bmat->matrix_mult(*amat,*dmodel,0,1);

pm[i]->matrix_add(*bmat,*qmat);

if(i%dn == 0)

amat->matrix_add(*pm[i],*rmat1);

else

amat->matrix_add(*pm[i],*rmat0);

kg->matrix_mult_inv(*pm[i],*amat,0,1);

for(j=0;j<pnum;j++) {

for(k=0;k<pnum;k++)

(*kfm_podc)(j,i) += (*dmodel)(j,k) * (*kf_podc)(k,i-1);

}

if(i%dn == 0) {

for(j=0;j<pnum;j++) {

(*kf_podc)(j,i) = (*kfm_podc)(j,i);

for(k=0;k<pnum;k++)

(*kf_podc)(j,i) += (*kg)(j,k)* (drms->podc(k,i/dn) - (*kfm_podc)(k,i));

} }

else {

for(j=0;j<pnum;j++) {

(*kf_podc)(j,i) = (*kfm_podc)(j,i);

for(k=0;k<pnum;k++)

(*kf_podc)(j,i) += (*kg)(j,k) * ((*fit_podc)(k,i) - (*kfm_podc)(k,i));

} }

amat->matrix_sub(*umat,*kg);

pp[i]->matrix_mult(*amat,*pm[i]);

}

//kalman smoother

puts("Kalman Smoother Estimation");

ks_podc = new matclass(pnum,sec*fpp);

for(i=0;i<pnum;i++)

(*ks_podc)(i,sec*fpp-1) = (*kf_podc)(i,sec*fpp-1);

ps = new matclass(*pp[sec*fpp-delay-1]);

for(i=sec*fpp-delay-1;i>=delay;i--) {

amat->matrix_mult(*pp[i],*dmodel,0,1);

kg->matrix_mult_inv(*amat,*pm[i+1],0,1);

amat->matrix_sub(*pm[i+1],*ps);

bmat->matrix_mult(*kg,*amat);

cmat->matrix_mult(*bmat,*kg,0,1);

ps->matrix_sub(*pp[i],*cmat);

for(j=0;j<pnum;j++) {

(*ks_podc)(j,i) = (*kf_podc)(j,i);

for(k=0;k<pnum;k++)

(*ks_podc)(j,i) += (*kg)(j,k) * ((*ks_podc)(k,i+1) - (*kfm_podc)(k,i+1));

} }

delete amat;

delete bmat;

delete cmat;

delete qmat;

delete kg;

delete rmat0;

delete rmat1;

delete umat;

delete ps;

//calcurate velocity from POD coefficients puts("Calcurate Estimated Velocity");

amat = new matclass(dnum,sec*fpp);

vec1 = new double[pnum];//POD coefficients vec2 = new double[dnum];//velocity

for(j=delay;j<sec*fpp-delay;j++) {

for(i=0;i<pnum;i++)

vec1[i] = (*ks_podc)(i,j);

drms->vel(vec1,vec2);

for(i=0;i<dnum;i++)

(*amat)(i,j) = vec2[i];

}

puts("Output Result");

fp = fopen("result.csv","w");

for(i=0;i<dnum;i++) {

for(j=0;j<sec*fpp;j++)

fprintf(fp,"%lf,",(*amat)(i,j));

fprintf(fp,"¥n");

}

fclose(fp);

delete amat;

delete[] vec1;

delete[] vec2;

delete drms;

delete fit_podc;

delete fit_coef;

delete dmodel;

delete kf_podc;

delete kfm_podc;

for(i=0;i<sec*fpp;i++) {

delete pp[i];

delete pm[i];

}

delete[] pp;

delete[] pm;

delete ks_podc;

delete pdata;

return 0;

}

//matclass.hpp

#ifndef _MATCLASS_DEFINED__

#define _MATCLASS_DEFINED__

#include <iostream>

#include <cmath>

#include <cstring>

#define TOR 1e-6

#define PI 3.14159265358979 using namespace std;

class matclass {

public:

//least square method:

// Ax=b,A is this matrix, b is matclass b, x is matclass c, returned values are stored in matclass c double lsq(matclass *b, matclass *c);

//return POD modes Phi_m,n double phi(int m, int n);

//return POD coefficients r_m,n double podc(int m, int n);

//calcurate velocity from POD coefficients void vel(double *podc, double *v);

//data(matclass data) is convert to POD coefficients(matclass podc) int podc(matclass *data, matclass *podc);

//matclass a - matclass b is stored in this matclass void matrix_sub(matclass a, matclass b);

//matclass a + matclass b is stored in this matclass void matrix_add(matclass a, matclass b);

//matclass a * matclass b is stored in this matclass //if matclass a is inversed, ia != 0

//if matclass b is inversed, ib != 0

void matrix_mult_inv(matclass a, matclass b, int ia, int ib);

//matclass a * matclass b is stored in this matclass

//if matclass a is transposed, ia != 0 //if matclass b is transposed, ib != 0 void matrix_mult(matclass a, matclass b);

//matclass a * matclass b is stored in this matrix void matrix_mult(matclass a, matclass b, int ta, int tb);

//return matrix size, m is row, n is column void getsize(int *m, int *n);

//calculate singular value decomposition of this matclass and store here int svd();

//calculate proper orthogonal decomposition of this matclass and store here //the number of the POD modes is npod

int pod(int npod);

//weight matrix for POD is stored as all the same value void set_weight(double val);

double& operator()(int m, int n) {

m_tempd = 0.;

if(m > m_m || 0 > m)

return m_tempd;

if(n > m_n || 0 > n)

return m_tempd;

return m_mat[m][n];

};

//transpose this matclass void transpose();

~matclass();

//allocated matrix size m*n matclass(int m, int n);

//allocated matrix size as same as matclass mat

//and initial values are the same values as matclass mat matclass(matclass& mat);

protected:

int svd(double tor, int snum, int flag);

int pod(int npod, double tor);

int pod_(int npod, double tor);

void delete_all();

matclass();

double **m_mat,**m_tmat,*m_weight,**m_lmat,**m_rmat,*m_singular,*m_eigen,**m_emat,**m_phi,**m_podc;

double m_tempd;

int m_m,m_n,m_snum,m_enum,m_pnum;

};

#endif

//matclass.cpp

#include "matclass.hpp"

int matclass::podc(matclass *data, matclass *podc) {

int md,nd,mp,np,i,j,k;

if(m_pnum <= 0) return 0;

data->getsize(&md,&nd);

podc->getsize(&mp,&np);

if(md != m_m) return 0;

if(nd != np) return 0;

if(mp != m_pnum) return 0;

for(i=0;i<m_pnum;i++) {

for(j=0;j<nd;j++) {

(*podc)(i,j) = 0.;

for(k=0;k<m_m;k++)

(*podc)(i,j) += m_phi[k][i] * m_weight[k] * (*data)(k,j);

} }

return 1;

}

void matclass::vel(double *podc, double *v) {

int i,j;

double ret = 0.;

for(i=0;i<m_m;i++) {

v[i] = 0.;

for(j=0;j<m_pnum;j++)

v[i] += m_phi[i][j] * podc[j];

} }

int matclass::svd(double tor, int snum, int flag) {

int i,j,k,num;