グラフ処理と動的シュミレーション

2000年5月19日 小林 茂雄

2001年6月2日 加筆(Ruby拡張ライブラリ作成)

目次

はじめに

シュミレーションにおいて第一に遭遇する困難は数式 モデルの構築である。 ここでは「数式モデル」を、あるモデルを表現する ために必要な要素(変数)と、それを計算するための 数式の集合という意味でもちいる。 数式モデルを構築するときの問題点と、求められる 解決手段を整理すると大体以下のようになるだろう。
1.どれくらい詳細な数式モデルを定義するべきか?

例えば、化学プロセスの装置モデルを定義する場合等 では、反応がある場合と無い場合では数式表現に大 変な差が出てくる。最初から反応を無視したモデルを 構築するのは簡単だが、それでは将来も反応は扱えない。 かといって複雑な数式モデルを定義すると、反応の無い 簡単な場合でも余計な計算や、ときとして有害な計算 をしてしまう可能性がある。

明らかに、複雑な数式モデルを与えても、条件に応じて 適切な数式を自動的に選ぶと同時に、不必要な計算は しないような機構が求められる。

2.変数条件等の変更に柔軟に対応できるか?

たとえば y = f(x) と数式を定義したとする。 これは x にある値を入れたときに y が計算 されるということである。通常は、これで何も 問題はない。ところが逆に、 y がある値(a としよう) になるときの x の値を計算したい場合がある。 数式処理システムなら、自動的に x = g(y) と並び替えが できるかもしれないが、そうでなければ f(x) - a => 0.0 を Newton 法等で解くことになる。

ここでは、与えた条件から自動的に連立(一般に非線形)方程式 を組み立てる機構が求められる。
本稿では、上記1.や2.に対してグラフ処理が 有効であることを説明する。

数式モデルのグラフ表現

数式モデルを組み立てるということを、グラフ処理の 表現でいえば「個々の変数の因果関係を有効グラフで 表現する」ことである。 グラフ処理では因果関係のみを扱う。 従って y = f(x) という数式は、グラフでは y <--- x というように x から出た矢印が y に入るというだけの意味になる。 f(x) の具体的な内容や、数値的に計算できるか どうかには立ち入らない。 数式処理システムでは、式を並び替えて(つまり、因果関係を変更して)y <--- x を x <--- y にできるが、本稿では「一旦定義した因果関係は普遍」 という前提で話を進める。

動的シュミレーションモデル

シュミレーションには時間によらない「静的モデル」 と、外乱の伝播状況などをシュミレートする「動的モデル」 がある。一般に、静的モデルは代数方程式系に帰着される。 一方、動的モデルは、数式表現に時間に依存する微分方程式 と代数方程式を含んだ系に帰着される。

本稿では微分方程式は連立常微分方程式を想定する。

連立常微分方程式の計算方法として、 いちばん簡単な Euler 法を例に考えてみよう。

連立常微分方程式系なので、時間に関係する 微係数(微分変数)と積分変数が連立常微分方程式の元数分だけ存在する。

シュミレーション開始である時刻ゼロの時間断面においては、 当然、境界条件から微係数や積分変数の値が与えられる。 そして、これらの値を元に依存関係から、その他の変数の 値が代数方程式を解くことにより計算される。 つまり、ある特定の時間平面上の全ての変数は「静的モデル」と同様に 計算できる。次に積分の刻み時間(dt 時間)後の積分変数の値が 連立常微分方程式を解くことで計算される。

Euler 法では
< dt 時刻後の積分変数の値 > = < 今の積分変数の値 > + (<今の微係数の値>)*dt
である。
そして、その他の変数の dt 時間後の値は 変化した積分変数の値をもとに代数方程式をとくことで計算される。 後は、同じ繰り返しである。

Euler 法や Runge-Kutta法では、ある時間断面での 情報のみから微係数と積分変数を計算できるので、 以上のようなイメージで計算が進行すると考えてよい。 ただし、過去の履歴から次の積分の刻みや微係数を 計算する多段解法ではこうはいかない。もちろん 過去の変数も含めて因果関係をグラフとして構築することで 同様に扱うことができるが、因果関係は複雑になる。

グラフ処理の適用項目

詳細は後述するとして、静的モデルを含む動的モデル にグラフ処理を適用する順序と項目は以下のようになる。
  1. 全ての変数間の因果関係を有向グラフとして構築
  2. 個々の変数に与えた条件に応じて不必要な計算(変数) の除去と、連立代数方程式の構築
  3. 連立代数方程式の解の存在性をチェック
  4. ループの検出
  5. 大型連立代数方程式を、連続した小連立方程式 に分解(ブロック三角化)
  6. 代数方程式系の計算順序の決定
  7. 時間に関する計算量の最適化
  8. 微分方程式を代数方程式に変換して定常状態の計算

グラフ処理の詳細

以下に、具体的な処理を順序に従って説明する。
全ての変数間の因果関係を有向グラフとして構築
変数相互の因果関係を有向グラフとして構築することは、 自明のこととして多くを説明する必要はないだろう。 x = f(y,z) という関数関係は y と z から出る矢印が x に入るという関係を内部的に定義するだけである。 この時点では x、y、z の各変数の「属性」は未定義である。

注意
便宜上、本稿で採用する関数表現は x = f(y,z) のように 式の左辺は常に変数が一個の場合を想定する。

変数の属性設定と不必要な計算の除去
変数相互の因果関係が有向グラフとして構築されたら、 個々の変数に属性(計算条件)を与える。 例えば、温度を計算する式が定義されていたとしても、 何らかの方法で温度を一定にすることが可能であり、 かつ、一定にすること自体には興味が無い場合がある。 この場合、温度を保持する変数を定数にしてから 入ってくる矢印を削除すればよい。このような「計算条件」 を列挙すると以下のようなものが考えられる。

S(Set)型指定
先の温度の例のように、計算する必要の無い変数に 与える属性である。S型指定された変数は定数として 扱われるので、入る矢印は属性指定した時点で全て削除される。
G(Given)型指定
これは、計算の結果として所望の値(定数)になってほしい場合に 設定する属性である。例えば x <-- y <-- z というグラフ 関係で y にG型指定を設定したとする。x から見れば、 y の値は常に一定なので y はS型と同じである。 一方、z は y を所望の値にするために自由に値を変動させる ことができなければならない(z にはG型やS型等の固定された 値を持つ属性を指定できない)。同時に y をG型指定するということは、 指定した値を a として a - f(z) ==> 0.0 を z を独立変数として (Newton 法などを用いて)解くことを指定していることに他ならない。 グラフ的には、y をG型指定した時点で2つ(y1,y2)に分割する。 従って x <-- y1 , y2 <-- z となる。そして y1 はS型として 扱う。x <-- y <-- z なる一つの接続関係が、 x <-- y1 と y2 <-- z の(独立した)二つの接続関係(Connected components)に分かれる ことに注意していただきたい。
A(Assumable)型
前のG型変数の例では z がA型変数となる。 A型変数は特に指定するものではない。 入る矢印が無く(つまり、式の左辺には現れない)、かつ、 値を自由に変動させることのできる変数は自動的に A型とみなされる。明らかにG型変数の数とA型変数の数は 一致しないと連立方程式は 解くことができない。さらに、個々のA型変数からG型変数へ 独立したルートが存在しないと、やはり連立方程式は 解くことができない。
R(Required)型
計算は目に見える形で出力しなければ意味が無い。 大規模なシミュレーションでは、数多ある変数を 全て印刷するのは得策ではない。値を印刷したい変数にのみ R型指定をする(R型以外は印刷されない)。 R型はS型やG型等と異なる意味合いの属性であり、 A型で同時にR型である変数もあり得る。 R型指定は「不必要な計算をしない」という意味で非常に重要である。 一般に、変数の因果関係は(S型やG型指定の結果も含めて)独立した 接続関係(Connected components) を有する別々の複数ブロックに分割される。 そして、R型変数が存在しない計算ブロックは、 ブロック全体を計算から除外できる。 さらに、同一ブロック中でも R型変数の計算に寄与しない変数は削除してかまわないのである。
I(Integral)型
常微分方程式を解くことで計算される積分変数である。 仮に、積分関数を Integral()、微分関数を Differentiate() とする(積分と微分は、互いに逆の 関係になるので、関数関係を表現するときはどちらか 一つの表現を用意すればよい)。 y = Integral(z) と表現すれば y は z を積分した ものと表現することになる。z = Differentiate(y) と表現することも可能である。本稿では y = Integral(z) の表現を採用する。 Euler 法をイメージした表現を用いれば 「次の時刻の y の値は、今の z を積分したもの」と 表現できる。 このとき y が I型変数とみなされる。 y = Integral(z) と表現しても、グラフ処理上は y <--- z を意味しないので、積分(微分)を扱う 関数関係はやや特殊である。 これは、ある特定の時間断面における代数方程式を 解くためにグラフ処理を適用することに起因する。 ある特定の時刻では y の値は変化しないので S型と同じように扱うことになる。定常状態が存在する 物理モデルでは、逆に y から z に至る何らかの代数関係 があるはずである(「定常状態の直接計算」を参照)。

上記以外で計算途中(ルート途中)の変数をX型と呼ぶ。 その他D(Divisible)型やDD(Divided)型があるが、 それらは「ループの検出」で説明する。



連立代数方程式の解の存在性
G型変数の説明でもふれたが、連立代数方程式が 解を持つには
  • A型変数の数とG型変数の数が等しい
  • A型変数からG型変数へは独立したルートが必要
という条件が満たされなければならない。 後者の条件をグラフ的に表現すれば「それぞれのA型変数とG型変数 が独立したルートを辿って一対一に対応付けできる」 というマッチング問題に帰着される。 A型変数が多い場合は、設定すべき条件が足りないし、 逆は設定し過ぎであることを意味する。

もちろん、上記二つの条件が満足されたからといって 「数値的に解が存在」することを保証するわけではない。

ループの検出
シュミレーションの世界ではループの存在は自然である。 例えば、原料や製品の一部を還流させるようなパイプライン モデルでは計算ループは避けがたい。 ループを関数的に表現すると y = f(y) となる。 ここで y はループ上に存在する変数である。 ループを検出した場合は y - f(y) ==> 0.0 という条件を 生成して、A型とG型の連立方程式に追加する。このとき y はA型であると同時にG型でもあるとして扱われる。 ところで、ループ上に複数存在する変数のどれを y として 選択すればよいかという問題がまだ残っている。 通常は、より多くのループに関係する(含まれる) y を選ぶ(選ばれた y を DD(Divided)型と呼ぶ)のが 自然である。しかし、数値計算上、解を得るには適切な 初期値を与えなければならない。このことを考慮すると、 初期値を与えやすい変数をなるべくDD型としたいわけで、 D(Divisible)型指定された変数を優先的にDD型に選定するような ロジックが必要である。


大型連立代数方程式を連続した小連立方程式に分解する方法
A型とG型変数が非常に多い大規模シミュレーションでは、当然、 大型の非線型連立代数方程式を解くことになる。しかし、多くの場合 (Newton 法での)微係数行列(Jacobian 行列)の要素は大半がゼロである スパース行列になることが多い。このことは Jacobian 行列の行と列を 並び替えることで小さな連続した小連立方程式に分解できる 可能性があることを意味する。 実際に Jacobian 行列を並べ替えなくても、以下のようにグラフ処理を 適用することで同じことが達成できる。

まず、全てのA型とG型変数を一対一に対応付ける(Jacobian 行列の 対角要素を非ゼロにすることに相当する)。 次に、全てのG型変数に対して、G型変数から出て自身に一対一対応付けされたA型変数に 入る仮想的な矢印を一本追加する。これにより、一対一対応付けされた全てのA-Gペア上に 小さなループがそれぞれ生成される。一つのA型変数から複数のG型変数へ至るルートが存在するのが 普通なので、この小さなループから出発してA-Gのルートを辿ることで次第により大きい ループへとマージすることが可能になる。 このような処理の結果、順序を持ったループが幾つか残ることになる。 そして、これらの順序を持ったループが小連立方程式になるのである。

代数方程式系の計算順序の決定
大型連立代数方程式を連続した小連立方程式に分解することで、 全ての変数の計算順序が決定する。ここで、計算順序を決定するロジックを 改めて説明することはしない。

時間に関する計算量の最適化
代数計算における計算順序が決定したとして、動的モデルの シミュレーションでは更に以下のような最適化が考えられる。

積分開始前に一度だけ計算する変数の抽出
S型変数のみからに直接または間接的に計算される変数は 最初に一度計算すれば十分である。
積分刻みに応じて毎回計算しなければならない変数の抽出
時刻、積分変数、微分変数と、それらの上流に位置する変数である。
印刷指定時刻にのみ計算しなければならない変数の抽出
時間に影響される変数を上流に持ち、それ自身は 時刻、積分変数、微分変数等を計算することに寄与しないもの。
積分終了時に一度計算するばよい変数の抽出
時間に全く影響されない(かつ影響もしない)もの。 つまり、いつ計算してもかまわないもの。


微分方程式を代数方程式に変換して定常状態の計算
動的モデルのシミュレーション(時間に関する積分)を続けると、 最終的に「定常状態」に達するのが普通である。しかし、定常状態 に達するまでの時間が長い場合は、一気に定常状態を代数的に 計算することが可能である。定常状態では微係数がゼロになるので z = Differentiate(y) の関係から微係数 z がゼロになるときの 積分変数 y を計算するわけである。このとき z は値ゼロのG型であり y は対応するA型として扱われる。代数計算では、y の初期値が重要である。 代数的に定常状態計算する前の y は時間と共に定常状態に向かって 変動しているので、しばらく積分してから代数計算に移行する のが良い。

終わりに

本稿の内容は、既に退官された東京大学の 伊理正夫先生のグループと、かつて私の上司であった 恒川純吉氏らがまとめ上げた成果を概説したものです。 恒川氏らの成果は市販のソフトウェアに組み込まれ 大きな実績をあげています。ただ、学会等で アナウンスされたことはあっても、今日まで余り一般的に 知られていないようなので、ここにまとめて 概説した次第です。

文献

理論的な裏付については、

Kazuo MUROTA: Structural Solvability and Controllability of Systems
Research Memorandum RMI 83-01,Department of Mathematical Engineering and Instrumentation Physics,Faculty of Engineering,University of Tokyo

等に詳細に記述されています。


小林 茂雄 (E-Mail:<shigeo@tinyforest.gr.jp>)