微分代数方程式とINDEXの低減 - PDF 無料ダウンロード

SICE プラントモデリング研究会微分代数方程式と INDEX の低減 2009/09/25 モデルベース開発推進室石塚真一

. 準備体操 : スティフ / 陽解法と陰解法いわゆるスティフ ( 硬い ) な微分方程式とは, 微分方程式を陽解法で数値的に解く場合に重要な概念である. 陽解法と陰解法の具体例微分方程式 : y = f( y,t) 陽解法 ( 前進 Euler 法 ): 陰解法 ( 後退 Euler 法 ): y y f y ( ) (, t ) = +h n+ n n, t n y = y + hf y n+ n n+ n+ * 時間軸シミュレータでは, 現時点の情報から未来を計算できる陽解法の方が都合が良いが, 安定性を保証できず, 刻み幅の設定が重要となる. スティフな微分方程式の定義ある区間 [0,b] において, 前進 Euler 法の安定性を保つための刻み幅が, 解の精度を満たすために要する刻み幅よりはるかに小さい場合, 初期値問題はこの区間でスティフである. * 参考文献 3) より. 2009 CYBERNET SYSTEMS CO.,LTD. All Rights Reserved. 4

どんな系がスティフか? 線形システムでは : 系の固有値によって決まり, 固有値の原点からの距離が大きく離れている場合. 物理的には : 系の共振周波数に大きな差がある場合. 系の減衰に大きな差がある場合. 具体的には : メカトロニクスで, 機械系の時定数と電気系の時定数に大きな差がある. 金属とゴムなど剛性の大きく異なる複合材料. 条件数 ( 固有値の原点からの距離の差を評価する尺度 ) 固有値の絶対値で評価 : * A: システム行列 ma λ min λ ( A) ( A) ma 特異値で比較 (MATLAB): min λ λ T ( A A) T ( A A) 2009 CYBERNET SYSTEMS CO.,LTD. All Rights Reserved. 5

スティフな微分方程式の例題 m 2 y 2 k 2 c 2 y, f m k c 運動方程式 : パラメータ : 2 自由度振動モデルの初期値応答 ( ) ( ) ( ) ( ) m = c k c k 2 2 2 2 m = c + k 2 2 2 2 2 2 m = 2( kg), c = 5 0, k = 5 0 m = ( kg), c = 0, k = 0 0 8 2 2 2 3 2009 CYBERNET SYSTEMS CO.,LTD. All Rights Reserved. 6

固定ステップシミュレーション a) 刻み幅 : 5 0 5 b) 刻み幅 : 2 0 5 Dormand-Prince 法による固定ステップシミュレーション Dormand-Prince 法は,Runge-Kutta 系列に属する 5 次のソルバ. 陽解法であるため, 安定性が保証されず, 刻み幅の選択によっては解が発散することがある. 2009 CYBERNET SYSTEMS CO.,LTD. All Rights Reserved. 8

可変ステップシミュレーション a) ode45 (Dormand-Prince) b) ode23 (Bogacki-Shampine) c) ode5s (Stiff Solver) d) ode23s (Stiff / 修正 Rosenbrock 法 ) 可変刻みシミュレーションソルバにより解が大きく異なり, ソルバの選択を慎重に行う必要がある. 低次のソルバ, あるいはスティフなソルバは, 減衰が高めに見積もられる傾向にある. 多くの場合, 計算時間は掛かるが,ode45 が精度が高い. 2009 CYBERNET SYSTEMS CO.,LTD. All Rights Reserved. 9

2. 微分代数方程式とは? 微分代数方程式とは, ひとことで言えば : 微分方程式と代数方程式の連立方程式である参考文献 3) に従って, 微分代数方程式の要点を説明する. 微分代数方程式 (DAE: Differential-Algebraic Equation, 以下 DAE) の最も一般的な形は, F ( y y ) t,, = 0 と表され, この形は陰的微分方程式とも呼ばれる. 特別な場合,(2) 式のような制約のある常微分方程式となる. ( t,, ) ( t,, ) = f z 0 = g z これは,(t) に関する常微分方程式 (2a) は, 追加した代数変数 z(t) によって決まり, 解は代数方程式の条件 (2b) も満たさなくてはならず, 微分代数方程式の半陽的多次元方程式となる. () (2a) (2b) 2009 CYBERNET SYSTEMS CO.,LTD. All Rights Reserved. 2

DAE の実問題 DAE(Differential Algebraic Equation) とは? ダイナミクスを表す微分方程式と, 拘束を与える代数方程式が連立した問題で, マルチボディダイナミクス, 油圧機器, 電気回路など, モデリング時に頻出する問題. 古典的な DAE 数値解法アプローチでは? ステップごとに常微分方程式ソルバ (Runge-Kutta, etc) と, 代数方程式ソルバ (Newton-Laphson, etc) を交互に解く. 問題はあるか? シミュレーション速度が低下する ( 毎ステップごとの代数方程式収束計算 ). 収束計算が入るのでリアルタイムシステムには工夫が必要! スライダークランク機構サーボ弁 2009 CYBERNET SYSTEMS CO.,LTD. All Rights Reserved. 3 電気回路

簡単な例題 : 単振り子右図の振り子の問題を, 平行座標系 (, 2 ) で表す. λ をLagrange 乗数とすると,Newtonの運動方程式は,(3) となる. = λ = λ g (3) 2 2 棒による長さの制約は (4) となり, 微分代数方程式となる. 2 2 + = (4) 2 階の微分方程式に書き改めると (5) となり, (2) の形となる. = 3 = 2 4 (5) = λ 3 = λ g 4 2 2 2 + = 2 2009 CYBERNET SYSTEMS CO.,LTD. All Rights Reserved. 4 o 2 θ l mg * 質量を無視できる長さの剛体の棒に大きさを無視できる質量の質点が付いた振り子. θ

3. INDEX の概念 INDEX は DAE と ODE の距離のようなもの. INDEXはy を y とt に関して解くときに必要な最小の微分回数のことである. F df dt ( t, y, y ) = ( t, yy,, y ) = p d F ( ( p ) t, yy +,,, y ) = dt 0 0 0 INDEX は p となる. 2009 CYBERNET SYSTEMS CO.,LTD. All Rights Reserved. 5

具体例 : 振り子 /2 右図の振り子の問題を, 平行座標系 (, 2 ) で表す. λ をLagrange 乗数とすると,Newtonの運動方程式は,(3) となる. = λ = λ g (3) 2 2 棒による長さの制約は (4) となり, 微分代数方程式となる. 2 2 + = (4) 2 階の微分方程式に書き改めると (5) となり, (2) の形となる. = 3 = 2 4 (5) = λ 3 = λ g 4 2 2 2 + = 2 2009 CYBERNET SYSTEMS CO.,LTD. All Rights Reserved. 8 o 2 θ l mg * 質量を無視できる長さの剛体の棒に大きさを無視できる質量の質点が付いた振り子. θ

微分と代入 ( 回目 ) 微分と代入 (2 回目 ) 微分と代入 (3 回目 ) = 3 = 2 4 = λ 3 = λ g 4 2 2 2 + = 2 2 + 2 = 0 2 2 2 + 2 = 0 3 2 4 具体例 : 振り子 2/2 ODE 化 ( ) ( ) 2 + + 2 + = 0 3 3 2 4 2 4 2 2 2 2 ( ) λ( ) + + g = 3 4 2 2 2 2 ( ) λ ( ) 2 + 2 λ 2 + 2 + g = 0 3 3 4 4 2 2 2 2 λ = 4λ( + 3 2 4) 3g4 INDEX-3 問題 2009 CYBERNET SYSTEMS CO.,LTD. All Rights Reserved. 9 0 陽的常微分方程式 3 2 4 d λ = 3 dt λ g 4 2 λ 4λ ( + 3 2 4) 3g 4 高 Inde DAE を通常の常微分方程式化するためには, 解析的微分と代入が必要! 数式処理アプローチ.

初期条件微分と代入 ( 回目 ) 微分と代入 (2 回目 ) 微分と代入 (3 回目 ) = 3 = 2 4 = λ 3 = λ g 4 2 2 2 + = 2 2 + 2 = 0 2 2 2 + 2 = 0 3 2 4 ( ) ( ) 2 + + 2 + = 0 3 3 2 4 2 4 2 2 2 2 ( ) λ( ) + + g = 3 4 2 2 ( ) 2 2 ( ) λ ( ) 2 + 2 λ 2 + 2 + g = 0 3 3 4 4 2 2 2 2 λ = 4λ + 3g 3 2 4 4 0 同時に満足する必要がある 2009 CYBERNET SYSTEMS CO.,LTD. All Rights Reserved. 20

Linear DAE の INDEX 低減 Linear DAE はディスクリプタ形式で表現できる. Ez = Az + b Non-singular な P と Q を見つける. PEQ 2009 CYBERNET SYSTEMS CO.,LTD. All Rights Reserved. 2 E 0 = y=q として, 左から P を掛ける ( 途中省略 ). Az E A b z z 0 = + A b 2 2 + bの微分を加える. 2 2 E A b z z A = + 0 b 2 2 New E : 正則になるまで繰り返す. 繰り返し回数 :INDEX

そもそも DAE 化する必要があるのか? 普通だったら o 点回りの極座標系で運動方程式を立てればすむのではないか? o I θ = T ml 2 = mg sinθ θ l 長さ l に対する拘束が陽に無い円弧状を動くことを案に仮定 sin 関数で表現 sin 関数は正確に計算できない. 2 θ 通常は許容誤差範囲内に入る場合が多いが, スライダークランク機構などの, 閉ループ機構などでは無視できなくなる場合がある. mg * 質量を無視できる長さの剛体の棒に大きさを無視できる質量の質点が付いた振り子. 2009 CYBERNET SYSTEMS CO.,LTD. All Rights Reserved. 22

5. 固い微分方程式 :Robertson DAE 線形拘束方程式 = 0.04 + 0 4 2 3 = 0.04 0 3 0 4 7 2 2 2 3 2 0= + + 2 3 拘束方程式を回微分して ODE 化 INDEX- DAE = 0.04 + 0 4 2 3 = 0.04 0 3 0 4 7 2 2 2 3 2 = 3 0 7 2 3 2 通常の陽的常微分方程式積分器だけで初期値問題として解ける! 2009 CYBERNET SYSTEMS CO.,LTD. All Rights Reserved. 25

ODE vs DAE with Simulinik Out s 0.04 Integrator Gain s Integrator 2 Out2 e4 Gain y_ 0.04 /s y_ y_ u 2 Math Function 3e7 Gain4 y_2 u 2 e4 3e7 y_'>y_ /s y_2'>y_2 2 y_2 y_2 f(z) Solve z f(z) = 0 Algebraic Constraint 3 y_3 3e7 Gain5 s Integrator2 3 Out3 Math Function y_ 3 = 0.04 + 0 4 2 3 = 0.04 0 3 0 4 7 2 2 2 3 2 = 3 0 7 2 3 2 ODE 問題 2009 CYBERNET SYSTEMS CO.,LTD. All Rights Reserved. 26 = 0.04 + 0 4 2 3 = 0.04 0 3 0 4 7 2 2 2 3 2 0= + + 2 3 IINEX- DAE 問題

3e7 Gain5 s Integrator s Integrator s Integrator2 ODE Out 0.04 2 3 Gain Out2 e4 Gain u 2 Math Function Out3 3e7 Gain4 DAE vs ODE: Simulation Results, e4* 2, 3 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0. 0 Robertson DAE problem converting to ODE 0 0 0 5 Time (s) 2 3 error: 0= + 2 + 3-0 - -2-3 0-5 order 2 0-5 Constrain error Error divergence 0= + + 2 3-4 0-3 0-0 0 3 0 5 0 7 Time (s) Robertson DAE problem with a Conservation Law 0.9 誤差のオーダ 0-6 order 0 0-6 Constrain error 0.04 y_ e4 y_2 u 2 Math Function 3e7 y_ 3 y_ y_ /s y_'>y_ 2 y_2 /s y_2 y_2'>y_2 DAE f(z) Solve z f(z) = 0 Algebraic Constraint 3 y_3, e4* 2, 3 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0. 0 0 0 0 5 Time (s) 2 3 error: 0= + 2 + 3 - -0.2 0= + + -0.4 2 3-0.6-0.8-0 -3 0-0 0 3 0 5 0 7 Time (s) 2009 CYBERNET SYSTEMS CO.,LTD. All Rights Reserved. 27

5.2 高 INDEX 問題 ODE 形式 2 ml = mg sin θ ODE 数値ソルバ s Integrator s Integrator Scope -9.8 Gain sin Trigonometric Function 高 INDEX DAE 形式 = λ = λ g 2 2 2 2 + = 2 数式操作による INDEX 低減と INDEX- ODE 数値ソルバシステム方程式の自動生成と INDEX 低減 2009 CYBERNET SYSTEMS CO.,LTD. All Rights Reserved. 28

ODE vs DAE: Simulation Results -position y-position position, y 0.5 0-0.5-0 2 3 4 5 Time (s) 0~5(s) ODE の結果 0~0 4 (s) -position y-position 0.5 0 0.5 0 2 3 4 5 Time (s) 0~5(s) DAE の結果 0~0 4 (s) 2009 CYBERNET SYSTEMS CO.,LTD. All Rights Reserved. 29

まとめ物理システムは自然に DAE 形式で表現される. DAE 問題では ODE 化すると, 拘束条件を満足しなくなることがあるので注意が必要. DAE 問題として解くことが精度的には有効. 多くの工学, 科学問題は高 INDEX DAE 問題となる. 高 INDEX DAE 問題における INDEX 低減には数式処理が必須. スティフなシステムを ODE として解く場合は注意が必要! 複合領域問題はスティフになりやすい. 物理システムのモデリング / シミュレーションには数式処理 + 数値解析技術が重要と考える. 2009 CYBERNET SYSTEMS CO.,LTD. All Rights Reserved. 30

参考文献 ) 一松信著, 数値解析, 朝倉書店,998 2) 日本機械学会編, 数値積分法の基礎と応用, コロナ社,2003 3) U.M. アッシャー,L.R. ペツォルド共著, 中森眞理雄監訳, 常微分方程式と微分代数方程式の数値解法, 培風館,2006 4) 戸川隼人, 微分方程式の数値計算, オーム社,99 5) Michael M. Tilller 著, 古田勝久監訳,Modelica による物理モデリング入門, オーム社,2003 2009 CYBERNET SYSTEMS CO.,LTD. All Rights Reserved. 3