微分代数方程式とINDEXの低減

Transcription

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

2 . 準備体操 2. 微分代数方程式とは? 3. INDEX の概念 4. INDEX の低減 5. ベンチマーク まとめ 発表内容 2009 CYBERNET SYSTEMS CO.,LTD. All Rights Reserved. 2

3 . 準備体操 : 初期値問題と境界値問題 y 2 点境界値問題 :FEM などの構造解析 y y 微分方程式の一般解 : 無数に存在 解の唯一性 : 特殊解を求める 多点境界値問題 : 工学的には特殊? y 2009 CYBERNET SYSTEMS CO.,LTD. All Rights Reserved. 3 t 初期値問題 : 時間軸シミュレータ

4 . 準備体操 : スティフ / 陽解法と陰解法 いわゆる スティフ ( 硬い ) な微分方程式とは, 微分方程式を陽解法で数値的に解く場合に重要な概念である. 陽解法と陰解法の具体例 微分方程式 : 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) より CYBERNET SYSTEMS CO.,LTD. All Rights Reserved. 4

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

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

7 Simulink によるシミュレーション Simulink モデル 条件数 ( 約 ) 2009 CYBERNET SYSTEMS CO.,LTD. All Rights Reserved. 7

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

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

10 代数方程式の解法 :Newton-Raphson 法 y 初期値 初期値によって見つかる解が変わる! 2009 CYBERNET SYSTEMS CO.,LTD. All Rights Reserved. 0

11 Simulink によるシミュレーション 2 y + y = とyの関係 Simulink によるモデル 初期値 + からスタート 2009 CYBERNET SYSTEMS CO.,LTD. All Rights Reserved. 初期値 - からスタート

12 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

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

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

15 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 INDEX は p となる CYBERNET SYSTEMS CO.,LTD. All Rights Reserved. 5

16 簡単な例題 制約のある微分方程式 y y 2 ( ) = q t = y y y q t = = () 2 拘束方程式微分方程式 2 ( ) y = y = q t q (t) を 2 回微分 INDEX CYBERNET SYSTEMS CO.,LTD. All Rights Reserved. 6

17 4. INDEX の低減 出典 : Dymola で使用されている INDEX 低減アルゴリズム :Pantelides アルゴリズム 参考文献 5) より 2009 CYBERNET SYSTEMS CO.,LTD. All Rights Reserved. 7

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

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

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

21 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 bの微分を加える. 2 2 E A b z z A = + 0 b 2 2 New E : 正則になるまで繰り返す. 繰り返し回数 :INDEX

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

23 直接ブロック線図で表現できないか? = λ 因果的モデリングシステムでは独立した微分方程式を関連付けることができない =? 2 = λ g CYBERNET SYSTEMS CO.,LTD. All Rights Reserved. 23

24 5. ベンチマーク 5. 固い微分方程式 :Robertson DAE 5.2 高 INDEX 問題 2009 CYBERNET SYSTEMS CO.,LTD. All Rights Reserved. 24

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

26 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 = = = ODE 問題 2009 CYBERNET SYSTEMS CO.,LTD. All Rights Reserved. 26 = = = IINEX- DAE 問題

27 3e7 Gain5 s Integrator s Integrator s Integrator2 ODE Out Gain Out2 e4 Gain u 2 Math Function Out3 3e7 Gain4 DAE vs ODE: Simulation Results, e4* 2, Robertson DAE problem converting to ODE Time (s) 2 3 error: 0= order Constrain error Error divergence 0= Time (s) Robertson DAE problem with a Conservation Law 0.9 誤差のオーダ 0-6 order 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, Time (s) 2 3 error: 0= = Time (s) 2009 CYBERNET SYSTEMS CO.,LTD. All Rights Reserved. 27

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

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

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

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