WO2022220077A1

WO2022220077A1 - 常微分方程式の数値計算装置、計算装置において常微分方程式を解く演算の実行方法、及びプログラムを記憶した記憶媒体

Info

Publication number: WO2022220077A1
Application number: PCT/JP2022/014497
Authority: WO
Inventors: 仁城所
Original assignee: ＷｈｉｔｅＲｏｏｋ合同会社
Priority date: 2021-04-13
Filing date: 2022-03-25
Publication date: 2022-10-20
Also published as: EP4325381A1; JP6905298B1; JP2022162801A

Abstract

状態を表す状態変数ベクトル項及び入力を表す入力ベクトル項に夫々第１及び第２の係数行列を乗算した項を含む常微分方程式であって、所定の計算ステップで定まる離散時間に関する差分方程式の形式を有し、第１及び第２の係数行列が夫々、行列指数関数を近似する、所定の計算ステップと上記定数行列との積の冪乗項により構成される有限な次数の第１の行列冪級数を含む、状態方程式が設定される。事前計算手段であるＹｕｍｍｙ行列計算手段５０２は、計算ステップを時間的に不変な時間不変計算ステップｈ_ｅｘｐと時間的に変化し得る時間可変計算ステップｈとに分離し、時間不変計算ステップｈ_ｅｘｐと定数行列Ａとの積の冪乗項により構成される有限な次数の第２の行列冪級数であるＹｕｍｍｙ行列と、Ｙｕｍｍｙ行列に基づく第１の定数係数行列及び第２の定数係数行列を、状態方程式の解法演算に先立って事前に計算する。状態方程式計算手段５０９は、Ｙｕｍｍｙ行列計算手段５０２が計算した第１及び第２の定数係数行列の夫々と、パルス信号計測手段５０８が順次算出する時間可変計算ステップとの間の簡単な乗算演算により夫々第１の係数行列及び第２の係数行列を算出しながら、差分方程式に対する解法演算を実行する。行列指数近似の次数Ｎがどんなに長くなっても、その冪級数部分は予めＹｕｍｍｙ行列として計算しておくことができるため、状態方程式の解法演算を精度高く短い計算時間で実行可能となる。

Description

常微分方程式の数値計算装置、計算装置において常微分方程式を解く演算の実行方法、及びプログラムを記憶した記憶媒体

　本発明は、物理現象を数学的にモデル化した常微分方程式を数値計算により解く演算を実行する数値計算装置、そのような計算装置において常微分方程式を解く演算の実行方法、及びプログラムを記憶した記憶媒体に関する。

　物理現象を数学的モデルとして表現する運動方程式や電子回路は、常微分方程式となる事が多い。しかしながら、多くの常微分方程式は解析的に求めることは困難である。このため、古くは１９世紀のＣａｒｌ　Ｒｕｎｇｅ、Ｍａｒｔｉｎ　Ｋｕｔｔａの時代より、様々な数値計算手法が提案されてきている。
　今日では、制御対象を接続することにより、その制御対象に対してフィードバック制御を行う制御装置が無数に存在する。その制御対象は、機械装置である場合もあり、或いは電子回路である場合もある。例えば、自動車には、ＥＣＵ（Ｅｌｅｃｔｒｏｎｉｃｓ　Ｃｏｎｔｒｏｌ　Ｕｎｉｔ：電子制御ユニット）と呼ばれる制御装置が多く搭載されている。夫々のＥＣＵは、例えばエンジン、トランスミッション、又はサスペンションなどの機械装置である制御対象に対して、フィードバック制御を行う。また、近年の車両はいくつかの異なる電圧を電子機器に使用しているため、電圧コンバータ装置を制御するＥＣＵはスイッチ素子の断続により所定の電圧あるいは電流となるよう、スイッチ素子の断続状態に対しフィードバック制御を行う。これらの制御装置は、制御対象を接続して初めて機能するという特性を持つ。従って、制御装置の開発においては、かつては、実際に制御対象を接続して実験／検証を行う必要があった。しかし、開発段階の制御装置を実際の制御対象に接続して検証を行った場合、制御装置による制御が不完全であることにより、制御対象が破壊されてしまうという事故が起こり得た。また、例えば月着陸船の制御装置のように、実際の使用環境での実験／検証が不可能な制御装置の開発においては、制御対象を実際に接続してその制御装置を開発することができなかった。或いは、開発段階において制御対象の実機が存在しない場合においても、制御対象を実際に接続してその制御装置を開発することができなかった。

　このような問題を解決するために、従来、ＨＩＬＳ（Ｈａｒｄｗａｒｅ　Ｉｎ　ｔｈｅ　Ｌｏｏｐ　Ｓｉｍｕｌａｔｏｒ）と呼ばれるテスト装置が知られている（例えば特許文献１）。ＨＩＬＳでは、機械装置又は電子回路である制御対象が計算装置におけるソフトウェア処理によって模擬される。この場合、制御対象の動作は、計算機内のソフトウェアによって実現される数学モデルで表現される。また、計算装置は、電圧又は電流等の実際の信号を入出力する物理的入出力手段を装備する。そして、計算装置のこれらの物理的入出力手段に、制御装置の実際の物理的入出力手段が接続される。計算装置は、常微分方程式を解く演算を実時間で実行することにより制御対象の動作を模擬する。そして、計算装置は実際の信号を制御装置に対して実時間で入出力することにより、ハードウェアである制御装置を実際に動作させてテストする。即ち、ＨＩＬＳは、リアルタイムシミュレータとして動作する。ＨＩＬＳにより、制御対象の実機を用意することなく、ＥＣＵなどの制御装置の設計及び検証を支援することが可能となる。

　更に、ＳＩＬＳ（Ｓｏｆｔｗａｒｅ　Ｉｎ　ｔｈｅ　Ｌｏｏｐ　Ｓｉｍｕｌａｔｏｒ）と呼ばれるテスト装置も知られている（例えば特許文献１）。ＳＩＬＳでは、制御装置自体、及びその物理的入出力手段も計算装置によって模擬される。そして、ＳＩＬＳは、制御装置と制御対象の全体の動作を、計算装置が実行するソフトウェア処理によって模擬する。

　ＨＩＬＳ又はＳＩＬＳで用いられる数学モデルは、制御対象である、機械装置の動作を表現する運動方程式や、電子回路の動作を表す回路方程式である。この方程式は、一般的に常微分方程式である状態方程式として表現することができる。特に、状態変数ベクトルに第１の係数行列を乗算し、入力ベクトルに第２の係数行列を乗算した各乗算結果の和が、状態変数ベクトルの１階の導関数に等しい線形状態方程式は、制御対象の挙動を表す数学モデルとして有用であり、非線形の場合であっても上記形式に近似して用いられることが多い。ここで、「状態変数ベクトル」とは、常微分方程式の解をいう。

　ＨＩＬＳにおいては、実時間で動作する制御装置に合わせて、線形状態方程式を実時間で解きながら、その結果導出される信号を制御装置に対して実時間で入出力する必要がある。従って、状態方程式の数値解法演算を上記実時間処理に間に合うように計算装置で高速に実行する技術は、極めて重要である。

　連続系の状態方程式を解くための計算方法としては、その状態方程式を離散信号形式で表現することにより差分方程式として計算する数値計算手法が良く知られている。この解法の代表例としては、Ｒｕｎｇｅ－Ｋｕｔｔａ法が知られている。

　ところで、線形状態方程式における係数行列の要素が一定時間において不変な定数であれば、その方程式は、所謂ＬＴＩ（Ｌｉｎｅａｒ　Ｔｉｍｅ　Ｉｎｖａｒｉａｎｔ：線形時間不変）システムの方程式となる。ＨＩＬＳ又はＳＩＬＳで用いられる物理モデルは、ＬＴＩシステムを仮定しても十分であることが多い。このようなＬＴＩシステムの状態方程式においては、係数行列を対角要素が互いに異なる固有値である対角行列に変換することができる場合がある。

　このような対角行列に含まれる固有値は、一般に複素数である。そして、複数の固有値の大きさが著しく異なる状態方程式は、解が周波数の比較的高い成分と低い成分を含む、いわゆるＳｔｉｆｆな方程式と呼ばれる（例えば特許文献２）。

特許第５６９２７３９号公報特開２００５－０２５６５１号公報

　上述のようなＬＴＩシステムを仮定できる状態方程式（常微分方程式）を、前述したＲｕｎｇｅ－Ｋｕｔｔａ法などによる差分方程式に変換して演算する場合について考える。このケースでは、差分方程式における離散時間の計算ステップと対角行列の固有値の関係が適切でない場合、計算がプラス無限大（∞）又はマイナス無限大（－∞）に発散することが知られている。特に、前述したＳｔｉｆｆな方程式を解く場合に、絶対値の大きな固有値でも計算が発散しないようにするためには、計算ステップｈを小さくしなければならない。しかしながら、計算ステップｈを小さくした場合、方程式の応答を支配する絶対値の小さな固有値に対しては、その計算ステップは過剰に小さくなりすぎる。この結果、計算ステップを小さくした場合、計算誤差が蓄積すると共に、それ以上に、例えばＨＩＬＳの場合には、計算ステップｈが意味する物理的時間以内で計算が間に合わなくなる、という事態を招いてしまうという問題があった。

　このような問題を解決するために、更に従来、複素平面のいわゆる左半平面全体を計算安定領域とする、後退Ｅｕｌｅｒ法、台形法（Ｔｕｓｔｉｎ法）、又はＧｅａｒ法と呼ばれる後退微分方程式を用いた解法が知られている。しかしながら、これらの解法は、陰解法と呼ばれ、計算が容易ではない。即ち、常微分方程式が前述したようなＬＴＩシステムの場合であっても逆行列を計算する必要があり、計算負荷が極めて高くなる。また、常微分方程式が非線形方程式であれば、Ｎｅｗｔｏｎ法などの反復計算が必要となり、この場合も計算負荷が極めて高くなる。このため、後退Ｅｕｌｅｒ法や台形法やＧｅａｒ法は、実時間での計算を要求されるＨＩＬＳに適した計算法ではなく、ＳＩＬＳに対しても制御装置のための検証時間が長くなってしまうという問題点があった。またＳＩＬＳで検証したのち、ＨＩＬＳを適用する場合、計算アルゴリズムが異なることよる検証精度の劣化を招くという問題点もあった。

　更に、計算ステップｈに対して、後退Ｅｕｌｅｒ法は１次の近似精度であるため、計算精度が低い。台形法は２次の近似精度となるが、原点から離れた固有値に対しては、本来は解が一定値に収束する場合にあっても、解が振動し一定値に収束しないという問題がある。一方、Ｇｅａｒ法においては、次数を上げると安定領域が狭まり、連続系常微分方程式の安定領域である左半平面内において、不安定な領域が残る。このため、Ｓｔｉｆｆな方程式を満足に解くことのできる計算アルゴリズムは、皆無であると言わざるを得ない。

　以上の考察からわかるように、計算装置において絶対値が小さな固有値と絶対値が大きな固有値を共に有するＳｔｉｆｆな方程式を適切に解く演算を実行することは、従来困難であった。

　そこで、本発明は、計算装置において、Ｓｔｉｆｆな常微分方程式を安定して解く演算を、計算負荷が低くかつ高精度に、実行可能とすることを目的とする。

　態様の一例の常微分方程式の数値計算装置は、制御対象の動作を数学的にモデル化した常微分方程式である状態方程式を解く演算を実行する計算装置であって、状態を表す状態変数ベクトル項に第１の係数行列を乗算した項と入力を表す入力ベクトル項に第２の係数行列を乗算した項を含む常微分方程式について、所定の計算ステップによって定まる離散時間間の差分関係を示す差分方程式の形式を有し、一定時間において不変な定数とみなせる定数行列を指数とする行列指数関数を近似する行列冪級数であって、所定の計算ステップと定数行列との積の冪乗項により構成される有限な次数の第１の行列冪級数を夫々含む第１の係数行列及び第２の係数行列が夫々定義される、状態方程式を設定する方程式設定手段と、計算ステップを時間的に不変な時間不変計算ステップと時間的に可変な時間可変計算ステップとに分離し、少なくとも、時間不変計算ステップと定数行列との積の冪乗項により構成される有限な次数の第２の行列冪級数を状態方程式の解法演算に先立って事前に計算する事前計算手段と、事前に計算された第２の行列冪級数に基づいて夫々第１の係数行列及び第２の係数行列を算出しながら、順次算出される前記時間可変計算ステップに基づいて、差分方程式に対する解法演算を実行する状態方程式計算手段と、を備える。

　本発明によれば、計算装置において、Ｓｔｉｆｆな常微分方程式を安定して解く演算を、計算負荷が低くかつ高精度に、実行することが可能となる。これにより、例えばＨＩＬＳが適用できる制御装置の開発対象の範囲を広げることが可能となる。その他、ＳＩＬＳはもとより、一般的な数値計算の分野においても、計算時間の短縮に貢献し、ＳＩＬＳとＨＩＬＳの検証結果の整合性を高め、様々な解析の効率化に供することが可能となる。

ｈ＝１の場合においてＮ次の行列指数近似法の安定領域を説明する図である。Ｒｕｎｇｅ－Ｋｕｔｔａ法において安定に計算できる固有値λと計算ステップhの積ｈλの複素平面領域を説明する図である。固定ステップ法によるコンバータ回路の模擬動作を説明する図である。部分的可変ステップ法によるコンバータ回路の模擬動作を説明する図である。計算装置の一実施形態を示すブロック図である。昇圧コンバータ装置の回路図の例を説明する図である。パルス信号計測手段の処理例を示すフローチャートである。状態方程式計算手段の処理例を示すフローチャートである。Ｙｕｍｍｙ法の安定領域を説明する図（その１）である。Ｙｕｍｍｙ法の安定領域を説明する図（その２）である。Ｙｕｍｍｙ法の安定領域を説明する図（その３）である。Ｙｕｍｍｙ法の時系列の計算結果を説明する図（その１）である。Ｙｕｍｍｙ法の時系列の計算結果を説明する図（その２）である。計算装置の他の実施形態を示すブロック図である。他の実施形態において、Ｙｕｍｍｙ行列記憶手段が記憶する各配列データの値を示すテーブル例を示す図である。他の実施形態において、Ｙｕｍｍｙ行列選択処理の例を示すフローチャートである。

　以下、本発明を実施するための形態（以下「本実施形態」と記載）について図面を参照しながら詳細に説明する。以下の説明においては、本実施形態について詳細に説明する前に、まず常微分方程式の数値計算による解法に関する一般的な技術背景について説明する。そして、その一般技術の問題点を具体的に明らかにした上で、本実施形態について詳細に説明する。

＜常微分方程式＞
　前述したように、ＨＩＬＳ又はＳＩＬＳで用いられる数学モデルは、制御対象である、機械装置の動作を表現する運動方程式や、電子回路の動作を表す回路方程式である。この方程式は、任意数の要素を持つ状態変数ベクトルＸ、時間を表す変数ｔ、及びベクトル関数Ｆ（ｔ，Ｘ）を用いて、下記（１）式の常微分方程式として表現されることがある。

　上記（１）式の常微分方程式は、公知の技術である動作点による区分線形化手法などにより、Ａ、Ｂ行列と、外力を表す入力ベクトルＵ（ｔ）とを用いた、線形状態方程式と呼ばれる下記（２）式の形式で表現できることが知られている。

　このように、線形状態方程式は、状態変数ベクトルＸ（ｔ）に係数行列Ａを乗算し、入力ベクトルＵ（ｔ）に係数ベクトルＢを乗算した各乗算結果ＡＸ（ｔ）とＢＵ（ｔ）の和が、状態変数ベクトルの１階の導関数

に等しい、という形を有する。

　ＨＩＬＳまたはＳＩＬＳにおいて（２）式の形式の状態方程式を解く演算を効率的に実行することは重要である。特に、ＨＩＬＳにおいては、実時間で動作する制御装置に合わせて、（２）式の状態方程式を実時間で解きながら、その結果導出される信号を制御装置に対して実時間で供給する必要がある。従って、状態方程式の数値解法演算を上記実時間処理に間に合うように計算装置で高速に実行する技術は、極めて重要である。

＜行列指数関数法＞
　ところで、上記（２）式の各係数行列Ａ、Ｂは必ずしも定数行列である必要はない。しかしながら、各係数行列Ａ、Ｂの要素が一定時間において不変な定数であれば、（２）式は、所謂ＬＴＩ（Ｌｉｎｅａｒ　Ｔｉｍｅ　Ｉｎｖａｒｉａｎｔ：線形時間不変）システムの線形状態方程式となる。ＨＩＬＳ又はＳＩＬＳで用いられる物理モデルは、ＬＴＩシステムを仮定しても十分であることが多い。前述した（２）式がＬＴＩシステムの状態方程式となっており、かつ入力Ｕ（ｔ）が一定値Ｕであるとみなせる場合、行列指数関数

を用いて、下記（３）式の一般解を得ることができる。

　　ここで、ｅ^Atは式（４）で与えられる。

　しかし、上記（４）式において無限大（∞）の項までの計算は事実上不可能なため、実用上は有限の整数Ｎで打ち切る下記（５）式を用いるか、パデ近似などの近似的手法を用いる。

　Ｎでの打切りにしろ、パデ近似にせよ、これらは近似値である。そこで、例えば上記（５）式の計算誤差を少なくするために、所定の時間間隔ｈを導入し、下記（６）式のように、ｈ間隔の差分方程式で、Ｘ_１、Ｘ_２、・・・、Ｘ_ｎというように離散計算を行うことは広く行われている。ここで、Ｕ（ｔ）は時間間隔ｈの間一定とし、Ｕ_ｎと表現する。

　本実施形態で、上記（６）式で示される差分方程式の表現法を、Ｎ次の行列指数近似法と呼ぶ。

＜計算安定性＞
　係数行列Ａ、Ｂが定数であるＬＴＩシステムでは、前述した（２）式は下記（７）、（８）式に示されるように固有値分解できる場合がある。

　上述の（８）式の固有値λ_ｓ（１≦ｓ≦ｎ）は、一般に複素数である。そして、前述したように、複数の固有値λ_ｓの大きさが著しく異なる状態方程式は、いわゆるＳｔｉｆｆな方程式と呼ばれる。

　また上記（８）式のように、Λは対角要素が固有値である対角行列である。従って、（２）式は下記（９）式で示されるスカラー連立方程式として扱うことができる。

　上記（９）式を数値計算する場合は、所定の計算ステップｈを導入し、差分方程式で近似して計算することができる。従来知られている簡便な前進Ｅｕｌｅｒ法を用いて具体例を示せば、下記（１０）式及び（１１）式が得られる。

　上記（１０）式及び（１１）式では、近似がなされているから、近似誤差が発生する。この近似誤差の存在により、元の常微分方程式と差分方程式の関係によっては、計算結果が発散することがある。上記（１０）式及び（１１）式の前進Ｅｕｌｅｒ法では、

の条件が満たされると、

となり、ｎが増えるに従って値は単調に増加し、最終的には∞／-∞となり、発散と呼ばれる状態となる。
　従って、（２）式は（９）式のスカラー連立方程式と等価であるから、（２）式のＡ行列の固有値が一つでも（１２）式の条件を満たせば、（２）式の解法計算は発散することとなる。常微分方程式が前述したＳｔｉｆｆな方程式である場合に、（１２）式の条件が満たされ易くなる。

　図１は、前述した（６）式で定義した本実施形態でいうところのＮ次の行列指数近似法のｈ＝１において、安定に計算できる固有値λの複素平面領域を説明する図である。
　図１において、横軸はλの値の実部を示す実部軸、縦軸はλの値の虚数部を示す虚数部軸である。
　図１の安定領域曲線１０１～１０８は夫々、前述の（６）式においてＮ＝１、２、３、４、５、６、７、及び８とした場合に計算したときの、発散せずに計算できる安定境界上でｈ＝１の場合の固有値λの値をプロットしたものである。
　すなわち、計算する係数行列Ａのすべての固有値λが、あるＮが与える上記安定領域として示した曲線の内側にあれば、このＮを用いた（６）式を用いて、発散せずに計算することが可能である。
　なお、λは（６）式の係数行列Ａと正則な行列Ｐを用いて前述の（７）式によって計算される。

　ここで、行列指数近似法の次数とは前述した（６）式におけるＮの値を指すものである。
　また、前述の（１０）式及び（１１）式の前進Ｅｕｌｅｒ法は、上記（６）式でＮ＝１とした場合であり、これを次数が１であるところの１次の行列指数近似法と呼称する。

　図１から理解されるように、行列指数近似法における次数Ｎが大きくなれば、安定領域も広がる。Ｎ＝∞であれば行列指数関数の厳密解となって、安定領域は左半平面となり、連続系の常微分方程式の安定領域と一致する。

＜Ｓｔｉｆｆな方程式＞
　前述した（７）、（８）式に示されるＬＴＩシステムを仮定できる状態方程式（常微分方程式）を、例えばＲｕｎｇｅ－Ｋｕｔｔａ法などの差分方程式に変換して演算する場合について考える。このケースでは、上記前進Ｅｕｌｅｒ法の場合と同様に、差分方程式における離散時間の計算ステップｈと、Ａ行列の固有値λの値とが適切でない場合、計算がプラス無限大（∞）又はマイナス無限大（－∞）に発散することが知られている。逆に発散せず安定に状態方程式を解くことができる条件は、計算アルゴリズムごとに異なる。

　図２は、Ｒｕｎｇｅ－Ｋｕｔｔａ法で安定に計算できる固有値λと計算ステップｈの積ｈλの複素平面領域を説明する図である。図２で、図１の場合と同様に、横軸は積ｈλの値の実部を示す実部軸、縦軸は積ｈλの値の虚数部を示す虚数部軸である。図２の安定領域曲線２０１、２０２、及び２０３は夫々、前述の（２）式の状態方程式を代表的な数値解法である１次、２次、及び４次のＲｕｎｇｅ－Ｋｕｔｔａ法で解いたときの、発散せずに計算できる安定境界上で積ｈλの値をプロットしたものである。積ｈλの計算で、ｈは上記Ｒｕｎｇｅ－Ｋｕｔｔａ法が実行されるときの計算ステップであり、λは（２）式の係数行列Ａと正則な行列Ｐを用いて前述の（７）式により計算される。

　例えば、４次のＲｕｎｇｅ－Ｋｕｔｔａ法は、下記（１４）式で示される。

　図２において、破線領域２０４は、Ｓｔｉｆｆな固有値に基づく積ｈλの存在領域である。また、破線領域２０５は、応答を支配する固有値に基づく積ｈλの存在領域である。前述したＳｔｉｆｆな方程式を解く場合に、絶対値が大きな固有値でも計算が発散しないようにするためには、計算ステップｈを小さくして、積ｈλの存在領域２０４が安定領域曲線２０１、２０２、又は２０３の内側に入るようにする必要がある。計算ステップｈを小さくすれば、図２において、積ｈλの存在領域が領域２０４から破線矢印２０６で示される方向に移動して安定領域曲線２０１、２０２、又は２０３の内側に入るため、発散しない計算が可能である。しかしながら、計算ステップｈを小さくした場合、方程式の応答を支配する絶対値の小さな固有値に基づく積ｈλの存在領域２０５付近では、その計算ステップｈは過剰に小さくなりすぎる。

　この結果、計算ステップｈを小さくした場合、計算誤差が蓄積すると共に、それ以上に、例えばＨＩＬＳの場合には、計算ステップｈが意味する物理的時間以内で計算が間に合わなくなる、という事態を招いてしまうという問題があった。

　なお、以上は（１）式のような一般的な常微分方程式に適用できるＲｕｎｇｅ－Ｋｕｔｔａ法を例に説明したが、常微分方程式が（２）式のように線形状態方程式の場合、（１４）式で説明した４次のＲｕｎｇｅ－Ｋｕｔｔａ法は（６）式で表した行列指数近似法のＮ＝４の場合と完全に一致し、行列指数近似法の場合においても上記説明したＲｕｎｇｅ－Ｋｕｔｔａの場合と同様な計算ステップｈと安定領域の関係が成立する。

＜Ｓｔｉｆｆな方程式に対応する数値計算アルゴリズム＞
　以上説明したＲｕｎｇｅ－Ｋｕｔｔａ法及び行列指数近似法は、広く使われる数値計算手法である。しかし、これらの数値計算法は、安定領域に限りがある。そして、Ｓｔｉｆｆな方程式の場合、前述したように、不安定になり計算が発散するか、過剰に小さな計算ステップｈを用いなければならないという問題がある。

　このような問題を解決するために、複素平面のいわゆる左半平面全体を計算安定領域として含む、下記（１５）式で示される後退Ｅｕｌｅｒ法、下記（１６）式で示される台形法、或いはＧｅａｒ法（数式は省略）として知られる、後退微分方程式の数値計算解法が一般的に用いられている。

　しかしながら、前述したように、これらの解法は、陰解法と呼ばれる、Ｘ_ｎ＋１を得るのにＸ_ｎ＋１を必要とする解法であり、計算が容易ではない。
　即ち、常微分方程式が（２）式に示されるような線形方程式であれば、逆行列を計算する必要がある。一般に逆行列の計算負荷は、行列が“ｎ×ｎ”要素のサイズを有する場合、ｎの３乗に比例すると言われている。従って、ｎが大きなシステムの場合には、計算負荷が極めて高くなる。
　また、常微分方程式が（１５）式又は（１６）式に示されるような非線形を含む一般的な関数で表される方程式であれば、前述したように、Ｎｅｗｔｏｎ法などの反復計算が必要となり、この場合も計算負荷が極めて高くなる。

　このため、これらの解法は、実時間での計算を要求されるＨＩＬＳに適した計算法ではなく、ＳＩＬＳに対しても制御装置のための検証時間が長くなってしまうという問題点があった。
　常微分方程式の解法をＨＩＬＳなどに適用する場合、ＨＩＬＳは別名リアルタイムシミュレータであるから、計算は所定の時間以内で終了することが要請される。従って、Ｅｕｌｅｒ法や台形法やＧｅａｒ法などの計算負荷の重いアルゴリズムは、採用しづらい。また、必ず所定の時間以内に完了することが求められているため、Ｎｅｗｔｏｎ法による反復計算など、計算終了時間が一定でないアルゴリズムを採用することは困難である。
　ＳＩＬＳであっても、シミュレーション時間が長くなるため、開発効率の低下を招くという問題が生じる。また、ＳＩＬＳで検証したのちにＨＩＬＳを適用する場合、計算アルゴリズムが異なることよる検証精度の劣化を招くという問題も生じる。

　更に、前述したように、計算ステップｈに対して、後退Ｅｕｌｅｒ法は１次の近似精度であるため、計算精度が低い。台形法は２次の近似精度であるが、原点から離れた固有値に対しては、本来は解が一定値に収束する場合にあっても、解が振動し一定値に収束しないという問題がある。更に、Ｇｅａｒ法においては、次数を上げると安定領域が狭まり、連続系常微分方程式の安定領域である左半平面内において、不安定な領域が残る。このため、Ｓｔｉｆｆな方程式を満足に解くことのできる計算アルゴリズムは、皆無であると言わざるを得ない。

　そこで、以下に説明する本実施形態による計算アルゴリズムは、複素平面上に極めて広い計算安定領域を持ちながら、最も計算負荷の低い前進Ｅｕｌｅｒ法と同等の計算負荷を持つという特徴を有する。本実施形態の詳細について、以下に説明する。

＜Ｙｕｍｍｙ法＞
　まず、本実施形態の動作原理について説明する。Ｎ次の行列指数近似法において、次数Ｎが大きい場合は行列指数関数を計算する計算負荷が高く、例えばパデ近似などの他の手法を用いることが多い。しかし、（２）式の係数行列Ａが定数で不変であり、かつ計算ステップｈが定数の場合には、前述した（５）式の級数部分も定数となる。この場合には、行列指数近似法で使用する級数部分は、Ａ行列が不変である一定期間の間に１度だけ計算すればよいこととなる。
　ＳＩＬＳにおいてもこの特性は有用であるが、特にＨＩＬＳ等への応用の場合は、上記行列指数関数は、システムが実際にリアルタイムに稼働する前の準備段階で計算しておけばよい。この場合には、（５）式の級数部分の計算にたとえ数秒程度の時間を要しても実用上問題にはならない。

　従って、前述した（６）式において、Ｎ＝２０やＮ＝１００という、計算負荷が極めて大きな数としても、実用性が損なわれることはなく、行列指数近似法は有用な計算手法である。
　しかしながら、計算ステップｈが定数であることは、応用範囲を狭めるものである。例えば、本願の発明者が、下記論文文献１又は２にて提案した「部分可変ステップ法」として知られる従来の数値計算手法においては、スイッチング回路のスイッチングタイミングに応じた回路挙動を正確に計算するため、計算ステップｈ_ｆは実時間稼働時点での計測値となっている。

　　＜論文文献１＞：城所　仁　中原　正俊　“DC/DCコンバータのリアルタイム・シミュレーションを実現する陽解法に基づく高速計算手法”電子情報通信学会論文誌B　Vol.J97-B,No.6,pp.-,Jun. 2014.
　　＜論文文献２＞：城所　仁　中原　正俊　“FPGAを用いたスイッチングコンバータのリアルタイム・シミュレーション”電子情報通信学会　研究会資料（技術研究報告）EE2013-55

　本実施形態による数値計算手法では、所定の定数ｈ_ｅｘｐと（上記論文文献ではｈに相当）、計測値に基づいて変化する変数であるｈ（上記論文文献ではｈ_ｆに相当）を分離し、下記（１７）式、（１８）式、又は（１９）式に基づく計算を実行する。以下、本実施形態によるこの数値計算手法を「Ｙｕｍｍｙ法」と呼び、（１７）式により計算されるＡ_ｙｍ行列を「Ｙｕｍｍｙ行列」と呼ぶ。

　上記（１７）式のＹｕｍｍｙ行列Ａ_ｙｍは予め計算できる定数であるから、ＳＩＬＳにおいては最初の１回だけ、ＨＩＬＳにおいては制御装置を接続して実時間で動作する前に、計算して記憶しておけばよい。実時間動作時に計算するのは（１８）式である。更に、（１９）式のように、定数係数行列（Ａ_ｙｍＡ）と、係数行列Ｂが定数と見なせるのであれば、更に定数係数行列（Ａ_ｙｍＢ）も、予め計算しておけば、この形式の計算規模は、前述した（１０）式及び（１１）式で示される前進Ｅｕｌｅｒ法と同等となる。従って、（１９）式の形式は、現在知られている計算負荷が最も低い数値計算手法である。
　なお、連続系との対応は（２０）式で示され、ｈが計算ステップ幅である。

　上記動作原理に基づく本発明の実施形態の詳細について、以下に説明する。本計算装置の実施形態は、パルス信号のＯＮ（オン）／ＯＦＦ（オフ）時間により所定の電圧を生成するコンバータ装置を制御対象として模擬するものである。

　まず、パルス信号のＯＮ／ＯＦＦ時間を扱う場合の基本的な問題について説明する。図３は、固定ステップ法によるコンバータ回路の模擬動作を説明する図、図４は、部分的可変ステップ法によるコンバータ回路の模擬動作を説明する図である。図３及び図４共に、前述の論文文献１又は２を基礎としている。

　通常、ＨＩＬＳでは、図３に示されるように、リアルタイム性の制約により、計算時間ｈを一定の値とする固定ステップ法を用いる。しかし、本実施形態のように、パルス時間のＯＮ／ＯＦＦ時間により動作するコンバータ装置を模擬する場合、次のような問題がある。即ち、この場合、パルス時間の分解能＝計算ステップ時間＝ｈとなり、パルス時間の分解能、即ちコンバータ装置のパルス時間のＯＮ／ＯＦＦ時間の時間分解能が著しく悪化する。
　一方、ＳＩＬＳでは、リアルタイム性の制約はないが、ＨＩＬＳとの整合を図るためにＨＩＬＳの場合と同様な計算を行うことが一般的である。このため、ＳＩＬＳにおいても、上記と同様の問題が生じる。

　これに対し、本願の発明者は、前述した論文文献１又は２において、部分可変ステップ法と呼ぶ、パルスＯＮ／ＯＦＦ時間の分解能を著しく向上させることができる手法を提案している。この手法では、図４に示されるように、基本的には一定の計算時間ｈで計算するものの、パルス信号のＯＮ／ＯＦＦ変化が生じた場合、その時間＝ｈ_ｆを計測し、この値を用いた計算を実行する。

＜具体的な実施形態＞
　上述したような計算方法を用いる計算装置の一実施形態について、図５、図６、図７、及び図８を用いて具体的に説明する。図５は本発明の一実施形態である計算装置５００の実施形態を示すブロック図、図６はシミュレーション対象であるコンバータ装置の回路図の例を説明する図、図７は図５の計算装置５００内のパルス信号計測手段５０８の処理例を示すフローチャート、図８は図５の計算装置５００内の状態方程式計算手段５０９の処理例を示すフローチャートである。

　図５は、本発明の一実施形態である、計算装置５００を示すブロック図である。計算装置５００は例えば、物理的な信号入出力手段を有する、いわゆるＨＩＬＳの実施形態であり、更に前述したように例えば、特には図示しない制御装置からのパルス信号により所定の電圧を生成するコンバータ装置を制御対象として模擬するものである。以降、この計算装置を、ＨＩＬＳ計算装置５００と記載する。

　ＨＩＬＳ計算装置５００は、特には図示しないが、所謂組込みコンピュータのハードウェアとして、ＣＰＵ（中央演算処理装置）、ＲＯＭ（リードオンリーメモリ）、ＲＡＭ（ランダムアクセスメモリ）、フラッシュメモリ装着装置などを有する。

　図５に示されるＨＩＬＳ計算装置５００において、通信手段５０１は、所謂イーサネットあるいはＣＡＮ通信、またはこれに準ずるものである。

　Ｙｕｍｍｙ行列計算手段５０２、Ｙｕｍｍｙ行列記憶手段５０３、ベクトル算出手段５０６、パルス信号計測手段５０８、及び状態方程式計算手段５０９は、上記ＣＰＵが上記フラッシュメモリ装置に装着されたメモリカードに記憶されたＨＩＬＳ制御プログラムを上記ＲＡＭに読み出して実行する機能ブロック群である。

　基準時間生成手段５０７は、ＨＩＬＳ計算装置５００のコンピュータが内蔵するタイマ装置のハードウェアである。

　デジタル信号入力手段５０４、電圧入力手段５０５、及び電圧出力手段５１０は、上記入出力インタフェース装置によって構成される。

　通信手段５０１には、ローカルエリアネットワークやインターネット等のネットワークまたは専用通信線を介してホストコンピュータ５２０が通信接続される。

　デジタル信号入力手段５０４及び電圧入力手段５０５は、計算装置側物理的入力手段として動作し、特には図示しない制御装置の制御装置側物理的出力手段のデジタル信号出力手段及び電圧出力手段に接続される。

　デジタル信号入力手段５０４は、制御装置からの複数の信号線の各電圧を監視し、夫々所定の閾値を超えた場合に“１”、超えない場合に“０”の各デジタル入力信号をベクトル算出手段５０６及びパルス信号計測手段５０８に出力する。

　電圧入力手段５０５は、例えば、複数のアナログ／デジタル変換器を内蔵し、制御装置からの複数のアナログ入力電圧値をデジタル入力電圧値に変換し、各デジタル入力電圧値をベクトル算出手段５０６に出力する。

　ベクトル算出手段５０６は、上記複数のデジタル入力信号の“１”又は“０”値及び複数のデジタル入力電圧値をベクトル要素とする入力ベクトルＵ（前述の（２）式のＵ（ｔ）に対応）を算出する。なお、入力ベクトルＵの第１要素としては値“１”を設定する。この第１要素は、状態方程式が定数項（前述した（２）式のＢＵ（ｔ）＝定数）を必要とした場合に、係数行列Ｂの値により任意の定数を生成するために使用する。

　一方、ホストコンピュータ５２０は、一般的なパーソナルコンピュータであってよく、例えばＨＩＬＳのコンバータ装置を模擬する数学モデルを生成するための回路シミュレーションソフトウェアのプログラムを実行することができる。ユーザは、例えばコンバータ装置を制御対象とし、その数学モデルを生成したい場合には、上記回路シミュレーションソフトウェアにより、例えば、図６に示される昇圧コンバータ装置の回路図をホストコンピュータ５２０のキーボード、マウス、及びディスプレイを使って作成する。

　図６の回路は、ハイサイドスイッチＱ_ｈとローサイドスイッチＱ_ｌを具備する。これらのスイッチは、交互にＯＮ／ＯＦＦ動作を行う。例えば外部より入力されたＰＷＭ（パルス波形変調）信号が、ハイ（Ｈｉ）レベルの時は、Ｑ_ｈがＯＮ状態、Ｑ_ｌがＯＦＦ状態であり、ロー（Ｌｏ）レベルの時は、Ｑ_ｈがＯＦＦ状態、Ｑ_ｌがＯＮ状態である。

　上記回路シミュレーションソフトウェアは、図６で示す回路において、スイッチ素子Ｑ_ｈがＯＮの時の抵抗値をｒ_ｓｈ、Ｑ_ｌがＯＦＦの時の抵抗値をｒ_ｓｌに設定する。そして、このソフトウェアは、（２）式で示した状態方程式のＡ、Ｂ行列と、下記（２５）式で示される出力方程式を構成するＣ、Ｄ行列とを、前述した論文文献１又は２の記述に準じて、例えば下記（２１）式から（２７）式の演算によって生成する。なお、（２）式は、簡略化して（２１）式として再掲する。

　前述した回路シミュレーションソフトウェアは、外部信号がＨｉレベルの場合に、図６で示す回路において、スイッチ素子Ｑ_ｈがＯＮの時の抵抗値をｒ_ｓｈ、Ｑ_ｌがＯＦＦの時の抵抗値をｒ_ｓｌに設定する。そして、上記ソフトウェアは、Ａ、Ｂ、Ｃ、及びＤ行列を、上記（２１）式から（２７）式の演算により生成し、これらの行列をＡ_H、Ｂ_H、Ｃ_H、及びＤ_H行列として出力する。

　前述した回路シミュレーションソフトウェアは、外部信号がＬｏレベルの場合は、図６で示す回路において、スイッチ素子Ｑ_ｈがＯＦＦの時の抵抗値をｒ_ｓｈ、Ｑ_ｌがＯＮの時の抵抗値をｒ_ｓｌに設定する。そして、上記ソフトウェアは、Ａ、Ｂ、Ｃ、及びＤ行列を、上記（２１）式から（２７）式の演算により生成し、これらの行列をＡ_Ｌ、Ｂ_Ｌ、Ｃ_Ｌ、及びＤ_Ｌ行列として出力する。

　更に、上記回路シミュレーションソフトウェアは、上記Ａ_Ｈ及びＡ_Ｌ行列を夫々固有値分解し、得られた各固有値の値に基づいて予め定められた手続きに従い、定数計算ステップｈ_ｅｘｐを算出する。

　以上の動作の結果、ホストコンピュータ５２０は例えば、デジタル信号入力手段５０４に入力されるパルス信号がハイレベル（ＯＮ）のときに（２）式の状態方程式を構成するＯＮ時係数行列Ａ_Ｈ、Ｂ_Ｈ、Ｃ_Ｈ、及びＤ_Ｈと、パルス信号がローレベル（ＯＦＦ）のときに（２）式の状態方程式を構成するＯＦＦ時係数行列Ａ_Ｌ、Ｂ_Ｌ、Ｃ_Ｌ、及びＤ_Ｌと、所定の定数である定数計算ステップｈ_expとを出力する。これらのデータは、ＨＩＬＳ計算装置５００の方程式設定手段である通信手段５０１で受信され、記憶される。

　図５のＨＩＬＳ計算装置５００において、事前計算手段であるＹｕｍｍｙ行列計算手段５０２は、通信手段５０１が受信した上記係数行列Ａ_Ｈ及びＡ_Ｌと定数計算ステップｈ_ｅｘｐとを受け取り、これらのデータに基づいて前述した（１７）式の演算により、Ｙｕｍｍｙ行列Ａ_ｙｍＨ及びＡ_ｙｍＬを計算し、それらのＹｕｍｍｙ行列をＹｕｍｍｙ行列記憶手段５０３に記憶する。Ｙｕｍｍｙ行列記憶手段５０３は、記憶しているＹｕｍｍｙ行列Ａ_ｙｍＨ及びＡ_ｙｍＬを状態方程式計算手段５０９に出力する。

　ＨＩＬＳ計算装置５００により例えばＨＩＬＳの模擬動作が開始されたら、デジタル信号入力手段５０４が特には図示しない制御装置からパルス信号を取得する。デジタル信号入力手段５０４は、入力された信号を値「０」又は「１」に変換し、変数ｐｕｌｓｅの値として、パルス信号計測手段５０８及びベクトル算出手段５０６に送る。ベクトル算出手段５０６は、上述のように与えられた数値を各要素とするＵベクトルを作成する。なお、Ｕベクトルの第１要素は定数項のため、必ず値「１」を設定する。

　次に、図５のパルス信号計測手段５０８の動作について説明する。図７は、パルス信号計測手段５０８の処理例を示すフローチャートである。まず、変数Ｃｏｕｎｔの値を「０」にリセットする（図７のステップ０１）。なお、以降の説明において、「変数」は、計算装置５００を構成する特には図示しないＲＡＭ等に記憶される。

　次に、ＨＩＬＳ計算装置５００内の基準時間生成手段５０７からの信号を調べ、これが変化するまで待機する（ステップ０２の判定がＮＯ（図中「Ｎ」と記載）の繰返し）。

　基準時間生成手段５０７からの信号が変化し、ステップ０２の判定がＹＥＳ（図７中「Ｙ」と記載）になったら、変数Ｃｏｕｎｔの値を＋１インクリメントする（図中「＋＋」と記載）（図７のステップ０３）。

　次に変数ｏｌｄ＿ｐｕｌｓｅの値が「０」であるか否かを判定する（図中「＝＝」は等しいことを示す演算記号）（図７のステップ０４）。

　ステップ０４の判定がＹＥＳならば、デジタル信号入力手段５０４より与えられる変数ｐｕｌｓｅの値が「１」であるか否かを判定する（図７のステップ０５）。

　ステップ０５の判定がＹＥＳであれば、ＰＷＭ信号が変化したとして、下記（２８）式で示される設定処理を実行する（図７のステップ０７）。

　ステップ０５の判定がＮＯであれば、ＰＷＭ信号はＬｏレベルで変化していないとして、下記（２９）式で示される設定処理を実行する（図７のステップ０８）。

　ステップ０４の判定がＮＯの場合は、デジタル信号入力手段５０４より与えられる変数ｐｕｌｓｅの値が「０」であるか否かを判定する（図７のステップ０６）。

　ステップ０６の判定がＹＥＳであれば、ＰＷＭ信号が変化したとして、下記（３０）式で示される設定処理を実行する（図７のステップ０９）。ここで、「ｃｏｕｎｔ」は変数ｃｏｕｎｔに設定されている値を意味する。

　ステップ０６の判定がＮＯであれば、ＰＷＭ信号はＨｉレベルで変化していないとして、下記（３１）式で示される設定処理を実行する（図７のステップ１０）。

　ステップ０７、０８、０９、又は１０の後、変数Ｃｏｕｎｔの値がｈ_ｅｘｐの値以上になったか否かを判定する（図７のステップ１１）。

　ステップ１１の判定がＮＯならば、ステップ０２に戻り、上記と同様の処理を繰り返し実行する。

　ステップ１１の判定がＹＥＳになったら、状態方程式計算手段５０９にタイミングパルスを発行する（図７のステップ１２）。
　Ｃｏｕｎｔ＝ｈ_ｅｘｐであればステップ１１の判定はＹＥＳであるので、タイミングパルスはｈ_ｅｘｐの時間間隔で発行される。

　ステップ１２の後、変数ｏｌｄ＿ｐｕｌｓｅに変数ｐｕｌｓｅの内容を設定する（ステップ１３）。

　ステップ１３の後、ステップ０１に戻って変数Ｃｏｕｎｔの値を「０」にリセットした後、上記と同様の処理を繰り返し実行する。

　以上の処理により、ｈ_Ｈ≦ｈ_ｅｘｐ、ｈ_Ｌ≦ｈ_ｅｘｐであることが保証される。

　ステップ１３で変数ｏｌｄ＿ｐｕｌｓｅに変数ｐｕｌｓｅの値を設定しているので、時間間隔ｈ_ｅｘｐの間にＰＷＭ信号が所謂チャタリングを起こすことにより、短時間にＨｉレベルとＬｏレベルを繰り返すような動作は無視される。時間間隔ｈ_ｅｘｐを超えた時間、ＨｉレベルあるいはＬｏレベルが継続した場合のみ、そのＰＷＭ信号が変化したと認識され、その変化がｈ_Ｈ及びｈ_Ｌの値に反映される。

　続いて、図５の状態方程式計算手段５０９の動作について説明する。図８は、図５の状態方程式計算手段５０９の処理例を示すフローチャートである。まず、状態変数ベクトルＸに初期値Ｘ_０を与え、変数ｏｌｄ＿ｐｕｌｓｅに値「０」を設定する（図８のステップ２０）。

　次に、パルス信号計測手段５０８より与えられるタイミングパルスを待機する（図８のステップ２１の判定がＮＯの繰返し）。

　タイミングパルスが入力しステップ２１の判定がＹＥＳになったら、変数ｏｌｄ＿ｐｕｌｓｅの値が「０」であるか否かを判定する（図８のステップ２２）。

　ステップ２２の判定がＹＥＳ（値「０」）であれば、前回処理におけるパルスはＬｏレベルであるので、出力変数ベクトルＹを、下記（３２）式の演算により算出する（図８のステップ２３）。

　ステップ２２の判定がＮＯ（値「１」）であれば、前回処理におけるパルスはＨｉレベルであるので、出力変数ベクトルＹを、下記（３３）式の演算により算出する（図８のステップ２４）。

　次に、変数ｐｕｌｓｅの値が「０」であるか否かを判定する（図８のステップ２５）。

　ステップ２５の判定がＹＥＳ（値「０」）であれば、下記（３４）式及び（３５）式で示される演算を実行する（ステップ２６）。
　なお、ステップ２１でタイミングパルスを待ち、タイミングパルスは図７のステップ１２でｈ_ｅｘｐの時間間隔で発行されるのであるから、以下の処理はｈ_ｅｘｐの時間間隔で行われる。

　ステップ２５の判定がＮＯ（値「１」）であれば、下記（３６）式及び（３７）式で示される演算を実行する（ステップ２７）。

　このとき、ＰＷＭ信号がＨｉからＬｏに変化していれば、ステップ２６の処理となり、まずＡ_ｙｍＨ及びＡ_Ｈを用いた計算を行い、次にＡ_ｙｍＬ及びＡ_Ｌを用いた計算となる。
　ステップ２６の処理は、ステップ２１においてｈ_ｅｘｐの時間間隔で発行されるタイミングパルスと同期させているため、ステップ２６の処理は、ｈ_ｅｘｐの時間の時間間隔において状態変数ベクトルＸを更新する処理は２回実行される。
　これに対し、ステップ２３、ステップ２４において行われる出力変数ベクトルの計算は１度のみである。すなわちステップ２６、ステップ２７において計算される状態変数ベクトルＸ_ｔｍｐの値に対しては、出力変数ベクトルＹは計算されず、後述するＨＩＬＳ計算装置５００内の電圧出力手段５１０より出力電圧として出力されない。
　なお、ＰＷＭ信号がＬｏからＨｉに変化せずＬｏレベルのままであれば、ｈ_Ｈ＝０であり、（３５）式のみを計算した場合と等価となるが、計算負荷を一定にするため（３４）式も計算する。
　ステップ２７における計算は、上述のＨｉレベルとＬｏレベルが逆転した関係となる。

　ステップ２６又は２７の後、変数ｏｌｄ＿ｐｕｌｓｅに変数ｐｕｌｓｅの値を設定する。これにより、状態変数ベクトルＸが計算されたときのＰＷＭパルスがＨｉレベルかＬｏレベルかを記憶し、ステップ２２の判定の根拠とする。

　図５において、電圧出力手段５１０は、状態方程式計算手段５０９が出力する出力変数ベクトルＹ_ｎに相当する電圧値を出力する。なお、出力変数ベクトルＹが更新される時間間隔、すなわちｈ_ｅｘｐの時間間隔において、この電圧値は更新される。図５には図示していないが、電圧出力手段５１０内のデジタル／アナログ変換器によりデジタル電圧信号から変換されて出力されたアナログ電圧値は、制御装置の制御装置側物理的入出力手段に入力される。
　制御装置は、このアナログ電圧値を、制御対象に設置されたセンサ装置からの信号と理解して、制御対象の挙動を判定し、所定の制御出力を出力する。
　この制御出力は、ＨＩＬＳ計算装置５００のデジタル信号入力手段５０４又は電圧入力手段５０５に入力される。このようにして、制御装置と制御対象であるＨＩＬＳ計算装置５００との間に、フィードバックのクローズドループ（閉ループ）が形成される。この結果、図５のＨＩＬＳ計算装置５００は、図６の回路図で示したコンバータ装置を模擬し、制御装置に制御対象を与えることができる。

＜Ｙｕｍｍｙ法の安定領域＞
　Ｙｕｍｍｙ法は前述した（１７）式及び（１８）乃至（１９）式により計算するが、前述した（７）式から（１１）式までの関係と同様に安定領域が定義される。この時、ｈ＝ｈ_ｅｘｐ＝１であれば、図１で示した領域となるが、ｈ≠ｈ_ｅｘｐの場合は異なる領域となる。

　図９、図１０、及び図１１は、本実施形態によるＹｕｍｍｙ法の安定領域を説明する図である。図９～図１１において、図１の場合と同様に、横軸は固有値λの値の実部を示す実部軸、縦軸はλの値の虚数部を示す虚数部軸である。図９において、９０１は、前進Ｅｕｌｅｒ法（ｈ＝０．５）の安定領域曲線である。９０２は、行列指数近似法（Ｎ＝４、ｈ＝１）の安定領域曲線である。９０３は、行列指数近似法（Ｎ＝８、ｈ＝１）の安定領域曲線である。９０４は、本実施形態によるＹｕｍｍｙ法（Ｎ＝８、ｈ_ｅｘｐ＝１、ｈ＝０．５）の安定領域曲線である。図１０において、１００１は、行列指数近似法（Ｎ＝８、ｈ＝１）の安定領域曲線である。１００２は、本実施形態によるＹｕｍｍｙ法（Ｎ＝８、ｈ_ｅｘｐ＝１、ｈ＝０．８）の安定領域曲線である。図１１において、１１０１は、前進Ｅｕｌｅｒ法（ｈ＝０．１）の安定領域曲線である。１１０２は、行列指数近似法（Ｎ＝８、ｈ＝１）の安定領域曲線である。１１０３は、本実施形態によるＹｕｍｍｙ法（Ｎ＝８、ｈ_ｅｘｐ＝１、ｈ＝０．１）の安定領域曲線である。前述した（２）式のＡ行列の全ての固有値が上記各安定領域曲線の内側の安定領域にあれば、状態方程式計算手段５０９での状態方程式の計算は安定である。

　図９、図１０、図１１に示すように、ｈ＜ｈ_ｅｘｐであるならば、安定領域は（６）式として示したＮ次の行列指数近似法において計算ステップサイズをｈ_ｅｘｐとした場合に対し、必ず広くなる。従って、図５の実施形態においては、ホストコンピュータ５２０が係数行列Ａ_Ｈ、Ａ_Ｌの固有値を計算し、これらの固有値の全てがＮ次打ち切りの行列指数近似法の安定領域に入るように、ｈ_ｅｘｐ及び次数Ｎを決定し、ＨＩＬＳ計算装置５００に送る。図５のパルス信号計測手段５０８は、ｈ_Ｈ若しくはｈ_Ｌがｈ_ｅｘｐより大きい場合はｈ_ｅｘｐとすることができる。以上により、状態方程式計算手段５０９における状態方程式の解法演算が不安定になることはない。

＜Ｙｕｍｍｙ法の安定性の証明＞
　本実施形態によるＹｕｍｍｙ法の安定性、即ち、ｈ＜ｈ_ｅｘｐであるならば、安定領域は計算ステップサイズをｈ_ｅｘｐとしたＮ次の行列指数近似法の安定領域より広くなることについて、以下に示す。
　（７）～（１１）式のように（１７）式はスカラーの連立式と等価であるから、以下にスカラー式で示す。また外部入力値ｕは安定性解析には無関係なため省略する。

　すなわち、ｈによる安定領域は、ｈｅｘｐによる安定領域を含有するため、ｈｅｘｐのものより広いということができる。
　なお、Ω（λ）＝０となるλのときは、

であり、Ω（λ）＝０となるλにおいて、ｈとｈ_ｅｘｐの安定領域は一致する。

　本実施形態によるＹｕｍｍｙ法の計算精度について説明する。前述の（１８）式及び（１９）式のＹｕｍｍｙ法の状態方程式（常微分方程式）の計算精度は、ｈ_ｅｘｐがｈに近い値であれば、Ｎ次の行列指数近似法に漸近するが、ｈ_ｅｘｐとｈが離れた値、例えばｈ＝０．２ｈ_ｅｘｐのような場合は、計算精度は劣化する。

　例えば、下記（４２）及び（４３）式で示される２つの運動方程式について考察する。

　上記（４２）式及び（４３）式の運動方程式を、下記の諸元のもとで、８次のＹｕｍｍｙ法で計算した結果を、図１２Ａ及び図１２Ｂに示す。なおこの場合、前進Ｅｕｌｅｒ法では計算が発散するため、比較対象としては、前述した（１５）式として示した後退Ｅｕｌｅｒ法を用いている。
　　ｍ_１＝３００　ｍ_２＝１
　　ｋｓ＝０．９　ｋ＝１０　ｃ＝１０
　図１２Ａは、前述の（４２）式及び（４３）式の横軸：時間（秒）に対する縦軸：「ｘ_１：質量体ｍ_１の変位［ｍ］」の運動方程式の解法結果を示している。１２０１（「◇」マーク）は実施形態によるＹｕｍｍｙ法（ｈ_ｅｘｐ＝１、ｈ＝０．２）の計算結果であり、１２０２（実線）は厳密解、１２０３（「＊」マーク）は後退Ｅｕｌｅｒ法（ｈ＝０．２）、１２０４（破線）は後退Ｅｕｌｅｒ法（ｈ＝１）である。
　図１２Ｂは、前述の（４２）式および（４３）式の運動方程式の、横軸：時間（秒）に対する縦軸：「ｘ_２：質量体ｍ_２の変位［ｍ］」の解法結果を示している。１２１１（「◇」マーク）は、本実施形態によるＹｕｍｍｙ法（ｈ_ｅｘｐ＝１、ｈ＝０．２）の計算結果であり、１２１２（実線）は、厳密解、１２１３（「＊」マーク）は後退Ｅｕｌｅｒ法（ｈ＝０．２）、１２１４（破線）は後退Ｅｕｌｅｒ法（ｈ＝１）である。

　図１２に示したように、ｈ_ｅｘｐ＝１、ｈ＝０．２の場合においては、（１５）式として示した後退Ｅｕｌｅｒ法でｈ＝０．２として計算したほうが厳密解に近いという結果となっている。

　このような計算精度の低下を補うためには、図１３の例に示されるように、ｈ_ｅｘｐを複数用意し、これに対応するＡ_ｙｍＨ、Ａ_ｙｍＬを複数計算し、これらすべてを記憶し、計測されたｈ_Ｈまたはｈ_Ｌの値により、上記複数の計算結果から１つを選択して計算するようにすればよい。

　図１３は、ｈ＝０．２においてｈ_ｅｘｐ＝０．２５を選択した場合の結果を示している。図１３Ａは、前述の（４２）式および（４３）式の運動方程式の、横軸：時間（秒）に対する縦軸：「ｘ_１：質量体ｍ_１の変位［ｍ］」の解法結果を示している。１３０１（「◇」マーク）は本実施形態によるＹｕｍｍｙ法（ｈ_ｅｘｐ＝０．２５、ｈ＝０．２）の計算結果であり、１３０２（実線）は厳密解、１３０３（「＊」マーク）は後退Ｅｕｌｅｒ法（ｈ＝０．２）、１３０４（破線）は後退Ｅｕｌｅｒ法（ｈ＝１）である。
　図１３Ｂは、横軸：時間（秒）に対する縦軸：「ｘ_２：質量体ｍ_２の変位［ｍ］」の、前述の（４２）式および（４３）式の運動方程式の、横軸：時間（秒）に対する縦軸：「ｘ_２：質量体ｍ_２の変位［ｍ］」の解法結果を示している。１３１１（「◇」マーク）は本実施形態によるＹｕｍｍｙ法（ｈ_ｅｘｐ＝０．２５、ｈ＝０．２）の計算結果であり、１３１２（実線）は厳密解、１３１３（「＊」マーク）は後退Ｅｕｌｅｒ法（ｈ＝０．２）、１３１４（破線）は後退Ｅｕｌｅｒ法（ｈ＝１）である。

　図１３より、ｈ＝ｈ_ｅｘｐであれば、その精度はＮ次の行列指数関数近似法と一致し、Ｎを大きくとれば、極めて高い精度を得ることができる。ｈ≠ｈ_ｅｘｐであっても、近似なｈ_ｅｘｐによるＡ_ｙｍＨ、Ａ_ｙｍＬを使用することで、高い精度を極めて高速な計算で得ることができることがわかる。

＜ＳＩＬＳ計算装置＞
　図１４は、計算装置の他の実施形態を示すブロック図である。図１４は、計算装置１４００がＳＩＬＳを構成する例である。以降、この計算装置をＳＩＬＳ計算装置１４００と記載する。図１４では、図５の通信手段５０１の代わりに、回路シミュレーション実行手段１４０２を有し、状態方程式において用いられる前述したＡ，Ｂ，Ｃ，Ｄ行列が設定される。図１４において、図５の場合と同じ参照番号を付したブロックは、図５の場合と同じ機能を有する。

　また、ＳＩＬＳ計算装置１４００は、制御装置制御則計算手段１４０１を備えている。図５のＨＩＬＳ計算装置５００では、制御装置は、デジタル信号入力手段５０４、電圧入力手段５０５、又は電圧出力手段５１０を介してＨＩＬＳとして動作する計算装置５００の外部に接続されていた。これに対して、図１４では、上記制御装置の動作を、ＳＩＬＳ計算装置１４００内の制御装置制御則計算手段１４０１によって模擬することができる。

　図１４のＳＩＬＳ計算装置１４００によれば、それ単体で、内蔵される制御装置制御則計算手段１４０１によって制御装置の制御則の動作を検証することが可能である。状態方程式計算手段５０９の出力は、制御装置制御則計算手段１４０１にフィードバックされると共に、ディスプレイ装置１４０４に表示させて確認することができる。制御装置制御則計算手段１４０１から状態方程式計算手段５０９への入力ベクトルＵ（ｔ）についても、ディスプレイ装置１４０４で確認できるようにしてもよい。

　ＳＩＬＳ計算装置１４００は、回路シミュレーション実行手段１４０２を内蔵し、図５のホストコンピュータ５２０が実行した回路シミュレーションソフトウェアを実行できる。回路シミュレーション実行手段１４０２は、例えばユーザによるキーボード装置１４０３の操作により入力された回路図（図６を参照）から、状態方程式を導き、図５のＨＩＬＳ計算装置５００の場合と同様に、Ａ_Ｈ、Ｂ_Ｈ、Ｃ_Ｈ、Ｄ_Ｈ行列及びＡ_Ｌ、Ｂ_Ｌ、Ｃ_Ｌ、Ｄ_Ｌ行列を生成する。

　また、回路シミュレーション実行手段１４０２は、状態方程式を解析し、安定に計算を行える、計算ステップｈ_ｅｘｐ０からｈ_ｅｘｐＭまでの、所定の定数であるＭ＋１個の計算ステップを算出し、これらをＹｕｍｍｙ行列計算手段５０２に送る。なお、これらの計算ステップの値は、ｈ_ｅｘｐＭ＜・・・＜ｈ_ｅｘｐ１＜ｈ_ｅｘｐ０が成立する順とする。

　Ｙｕｍｍｙ行列計算手段５０２は、ｈ_ｅｘｐ０からｈ_ｅｘｐＭの各値に基づき、Ａ_ｙｍ０Ｈ～Ａ_ｙｍＭＨ行列と、Ａ_ｙｍ０Ｌ～Ａ_ｙｍＭＬ行列を計算し、これらをＹｕｍｍｙ行列記憶手段５０３に記憶させる。Ｙｕｍｍｙ行列記憶手段５０３は、これらのデータを、図１６のテーブルに示される対応関係で、各配列のデータｈ_ｅｘｐ［ｉ］、Ａ_ｙｍＬ［ｉ］、Ａ_ｙｍＨ［ｉ］（０≦ｉ≦３）に記憶させ、また各データセット毎にＹｕｍｍｙ法における次数Ｎも設定する。

　模擬動作が開始されると、Ｙｕｍｍｙ行列記憶手段５０３は、パルス信号計測手段５０８から送られてくる計算ステップｈ_Ｈ及びｈ_Ｌの値により、Ａ_ｙｍＨ、Ａ_ｙｍＬ行列を選択し、これらを状態方程式計算手段５０９に送る。
　図１６は、Ｙｕｍｍｙ行列記憶手段５０３が実行するＹｕｍｍｙ行列選択処理の例を示すフローチャートである。

　まず図１６のテーブルに例示されるように、各配列のデータｈ_ｅｘｐ［ｉ］、Ａ_ｙｍＬ［ｉ］、Ａ_ｙｍＨ［ｉ］（０≦ｉ≦３）を設定し、各データセット毎にＹｕｍｍｙ法における次数Ｎも設定する。そして、繰り返しの変数ｉを値「０」にリセットする（図１６のステップ３０）。

　次に、図１６において、ステップ３２で変数ｉの値を＋１ずつインクリメントしながら、ステップ３１で配列データｈ_ｅｘｐ［ｉ］の値が計算ステップｈ_Ｌの値未満となったと判定するまで、各配列データｈ_ｅｘｐ［ｉ］の値をチェックする動作を繰り返す（ステップ３１の判定がＮＯ→ステップ３２→ステップ３１の判定の繰返し）。

　配列データｈ_ｅｘｐ［ｉ］の値が計算ステップｈ_Ｌの値未満となってステップ３１の判定がＹＥＳになると、そのときの変数ｉの値に対応する配列データＡ_ｙｍＬ［ｉ］の値をＡ_ｙｍＬ行列としてセットする（図１６のステップ３３）。なお、図１５のテーブルに示されるように、配列データｈ_ｅｘｐ［ｉ］の値は、ｉが増加するに従って値が最大値（＝ｈ_ｅｘｐ［０］）から減ってゆくように設定されており、最後の値は「－１」のため、ステップ３１の判定は何れかのｉの値で必ずＹＥＳになる。

　次に、図１６において、ステップ３４で変数ｉの値を「０」にリセットした後、ステップ３６で変数ｉの値を＋１ずつインクリメントしながら、ステップ３５で配列データｈ_ｅｘｐ［ｉ］の値が計算ステップｈ_Ｈの値未満となったと判定するまで、各配列データｈ_ｅｘｐ［ｉ］の値をチェックする動作を繰り返す（ステップ３１の判定がＮＯ→ステップ３２→ステップ３１の判定の繰返し）。

　配列データｈ_ｅｘｐ［ｉ］の値が計算ステップｈ_Ｈの値未満となってステップ３５の判定がＹＥＳになると、そのときの変数ｉの値に対応する配列データＡ_ｙｍＨ［ｉ］の値をＡ_ｙｍＨ行列としてセットする（図１６のステップ３７）。このときも、配列データｈ_ｅｘｐ［ｉ］の値は、ｉが増加するに従って値が最大値（＝ｈ_ｅｘｐ［０］）から減ってゆくように設定されており、最後の値は「－１」のため、ステップ３７の判定は何れかのｉの値で必ずＹＥＳになる。

　図１５のテーブルから理解できるように、例えば

である。従って、Ａ_ｙｍＬ、Ａ_ｙｍＨの計算に使用するＮは、使用するｈ_ｅｘｐの値が小さくなれば計算精度が向上すること、及び計算安定領域も拡大することに鑑みて、図１５のテーブルに示されるように、ｈ_ｅｘｐの値が小さくなれば、Ｎの値も小さくすることにより、Ａ_ｙｍＬ、Ａ_ｙｍＨを計算する時間を短縮することができる。
　前述したように、ｈ_ｅｘｐ３＜ｈ_ｅｘｐ２＜ｈ_ｅｘｐ１＜ｈ_ｅｘｐ０の関係がある。このため、図１６に示されるフローチャートの制御処理により、例えばｈ_ｅｘｐ０≦ｈ_Ｈ＜ｈ_ｅｘｐ１であるｈ_Ｈであれば、ｈ_ｅｘｐ０で計算したＡ_ｙｍＬ行列が選択される結果となる。この結果、図１４の状態方程式計算手段５０９は、精度の高い計算が安定に行えるＡ_ｙｍＬ行列を用いて、計算を行えることになる。

　このように、図１４の構成を有するＳＩＬＳ計算装置１４００においては、図５のＨＩＬＳ計算装置５００におけるデジタル信号入力手段５０４や電圧出力手段５１０を備える必要なく、制御装置制御則計算手段１４０１により、テストを行う制御則を実行し、模擬された制御対象が適切に制御しうるかを実験／検証することが可能となる。

　以上説明したように、本実施形態による計算装置は、Ｓｔｉｆｆな常微分方程式を安定して解く演算を、計算負荷が低くかつ高精度に、実行することが可能となる。これにより、例えばＨＩＬＳが適用できる制御装置の開発対象の範囲を広げることが可能となる。その他、ＳＩＬＳはもとより、一般的な数値計算の分野においても、計算時間の短縮に貢献し、様々な解析の効率化に供することが可能となる。

　本実施形態で説明したＹｕｍｍｙ法は、複素平面上に極めて広い計算安定領域を持ちながら、最も計算負荷の低い前進Ｅｕｌｅｒ法と同等の計算負荷を持つ、という特徴を有する画期的な計算アルゴリズムである。
　従って、ＨＩＬＳなど、計算時間に制約がありながら、Ｓｔｉｆｆな方程式に対応しなければならない用途に適用することで、従来では困難であった、図６で例示したような機構や電子回路などの模擬、即ち、リアルタイムのシミュレーションを可能とするものである。
　加えて、ＨＩＬＳ以外の、物理的信号を介さず計算機内で制御機器の制御アルゴリズムを検証するＳＩＬＳ、更には一般的な数値計算の分野においても、計算時間の短縮に貢献し、ＳＩＬＳとＨＩＬＳの検証結果の整合性を高め、様々な解析の効率化に供するものである。

　１０１、１０２、１０３、１０４、１０５、１０６、１０７、１０８、２０１、２０２、２０３、９０１、９０２，９０３，９０４、１００１、１００２、１１０１、１１０２、１１０３　安定領域曲線
　５００　ＨＩＬＳ計算装置
　５０１　通信手段
　５０２　Ｙｕｍｍｙ行列計算手段
　５０３　Ｙｕｍｍｙ行列記憶手段
　５０４　デジタル信号入力手段
　５０５　電圧入力手段
　５０６　ベクトル算出手段
　５０７　基準時間生成手段
　５０８　パルス信号計測手段
　５０９　状態方程式計算手段
　５１０　電圧出力手段
　５２０　ホストコンピュータ
　１４００　ＳＩＬＳ計算装置
　１４０２　回路シミュレーション実行手段
　１４０３　キーボード装置
　１４０４　ディスプレイ装置

Claims

　制御対象の動作を数学的にモデル化した常微分方程式である状態方程式を解く演算を実行する計算装置であって、
　状態を表す状態変数ベクトル項に第１の係数行列を乗算した項と入力を表す入力ベクトル項に第２の係数行列を乗算した項を含む常微分方程式について、所定の計算ステップによって定まる離散時間間の差分関係を示す差分方程式の形式を有し、一定時間において不変な定数とみなせる定数行列を指数とする行列指数関数を近似する行列冪級数であって、所定の計算ステップと前記定数行列との積の冪乗項により構成される有限な次数の第１の行列冪級数を夫々含む前記第１の係数行列及び前記第２の係数行列が夫々定義される、状態方程式を設定する方程式設定手段と、
　前記計算ステップを時間的に不変な時間不変計算ステップと時間的に可変な時間可変計算ステップとに分離し、少なくとも、前記時間不変計算ステップと前記定数行列との積の冪乗項により構成される前記有限な次数の第２の行列冪級数を前記状態方程式の解法演算に先立って事前に計算する事前計算手段と、
　前記事前に計算された前記第２の行列冪級数に基づいて夫々前記第１の係数行列及び前記第２の係数行列を算出しながら、順次算出される時間可変計算ステップに基づいて、前記差分方程式に対する解法演算を実行する状態方程式計算手段と、
　を備える常微分方程式の数値計算装置。
　前記定数行列をＡ、前記計算ステップをｈ、前記有限な次数をＮとして、前記第１の行列冪級数は、

と定義され、
　前記（１）式として定義される前記第１の行列冪級数を用いて、

の演算式により前記第１の係数行列が定義され、
　前記定数行列Ａの逆行列をＡ^－１、第２の定数行列をＢ、単位行列をＩとして、前記（１）式として定義される前記第１の行列冪級数を用いて、

の演算式により前記第２の係数行列が定義され、
　前記方程式設定手段は、
　前記計算ステップｈによって定まる前記離散時間をｎ、Ｘ_ｎ及びＸ_ｎ＋１を夫々前記離散時間ｎ及びｎ＋１における前記状態変数ベクトル項、Ｕ_ｎを前記離散時間ｎにおける前記入力ベクトル項として、前記（２）式の演算及び前記（３）式の演算に基づいて、前記状態方程式である前記差分方程式を、

として設定し、
　前記事前計算手段は、
　前記計算ステップｈを時間的に不変な時間不変計算ステップｈ_ｅｘｐと時間的に可変な時間可変計算ステップｈとに分離し、前記第２の行列冪級数を、

として事前に計算し、
　更に前記（５）式に基づく第１の定数係数行列（Ａ_ｙｍＡ）及び第２の定数係数行列（Ａ_ｙｍＢ）を夫々事前に計算し、
　前記状態方程式計算手段は、前記第１の定数係数行列（Ａ_ｙｍＡ）、前記第２の定数係数行列（Ａ_ｙｍＢ）、及び前記順次算出される前記時間可変計算ステップｈに基づいて、

で示される形式の前記（４）式に対応する差分方程式に対して前記解法演算を実行する、
　請求項１に記載の常微分方程式の数値計算装置。
　前記事前に算出される前記第２の行列冪級数を記憶する記憶手段を更に備え、
　前記状態方程式計算手段は、前記記憶手段に記憶されている前記第２の行列冪級数を用いて、前記差分方程式に対して前記解法演算を実行する、
　請求項１又は２に記載の常微分方程式の数値計算装置。
　前記時間不変計算ステップは、前記状態方程式の固有値と前記時間不変計算ステップとの積が前記状態方程式の計算結果が安定となる安定領域内に存在するように設定され、
　前記時間可変計算ステップの算出において、算出された前記時間可変計算ステップの値が前記時間不変計算ステップの値よりも大きくなった場合に、前記時間可変計算ステップの値を前記時間不変計算ステップの値として算出する、
　請求項１乃至３の何れか一項に記載の常微分方程式の数値計算装置。
　前記数値計算装置は、制御装置の開発において、前記制御装置が制御を行う制御対象を模擬する装置であって、
　Ｘ_ｎを離散時間ｎにおける前記状態変数ベクトル項、Ｕ_ｎを前記離散時間ｎにおける前記入力ベクトル項とし、Ｃ及びＤを所定の定数行列として、出力変数ベクトルＹ_ｎを、以下の（７）式で示される演算により算出し、

　前記出力変数ベクトルＹ_ｎに応じた物理的な信号を前記制御対象に供するにおいて、前記時間不変計算ステップの時間間隔で前記（７）式を算出する出力変数算出手段を更に備える、
　請求項１乃至４の何れか一項に記載の常微分方程式の数値計算装置。
　前記計算装置は、制御装置の開発において、前記制御装置が制御を行う制御対象を模擬する装置であって、
　前記制御装置から得られる信号により正の数で前記時間可変計算ステップを算出するパルス信号計測手段を更に備える、
　請求項１乃至５の何れか一項に記載の常微分方程式の数値計算装置。
　前記事前計算手段は、
　前記時間不変計算ステップを複数設定し、
　前記複数の時間不変計算ステップの夫々に対して複数の前記第２の行列冪級数の夫々を、前記状態方程式の解法演算に先立って事前に計算し、
　前記状態方程式計算手段は、前記事前計算手段が事前に計算した前記複数の第２の行列冪級数から１つを選択して、前記差分方程式に対する解法演算に使用する、
　請求項１乃至６の何れか一項に記載の常微分方程式の数値計算装置。
　前記事前計算手段は、前記複数の時間不変計算ステップの夫々に対して複数の前記第２の行列冪級数の夫々を前記状態方程式の解法演算に先立って事前に計算するときに、前記複数設定した時間不変計算ステップの夫々に対応した次数で計算を実行する、
　請求項７に記載の常微分方程式の数値計算装置。
　制御対象の動作を数学的にモデル化した常微分方程式である状態方程式を解く演算を実行する計算機に、
　状態を表す状態変数ベクトル項に第１の係数行列を乗算した項と入力を表す入力ベクトル項に第２の係数行列を乗算した項を含む常微分方程式について、所定の計算ステップによって定まる離散時間間の差分関係を示す差分方程式の形式を有し、一定時間において不変な定数とみなせる定数行列を指数とする行列指数関数を近似する行列冪級数であって所定の計算ステップと前記定数行列との積の冪乗項により構成される有限な次数の第１の行列冪級数を夫々含む前記第１の係数行列及び前記第２の係数行列が夫々定義される、状態方程式を設定する処理と、
　前記計算ステップを時間的に不変な時間不変計算ステップと時間的に可変な時間可変計算ステップとに分離し、少なくとも、前記時間不変計算ステップと前記定数行列との積の冪乗項により構成される前記有限な次数の第２の行列冪級数を前記状態方程式の解法演算に先立って事前に計算する処理と、
　前記事前に計算された前記第２の行列冪級数に基づいて夫々前記第１の係数行列及び前記第２の係数行列を算出しながら、順次算出される時間可変計算ステップに基づいて、前記差分方程式に対する解法演算を実行する処理と、
　を実行させる、計算装置において常微分方程式を解く演算の実行方法。
　制御対象の動作を数学的にモデル化した常微分方程式である状態方程式を解く演算を実行する計算機に、
　状態を表す状態変数ベクトル項に第１の係数行列を乗算した項と入力を表す入力ベクトル項に第２の係数行列を乗算した項を含む常微分方程式について、所定の計算ステップによって定まる離散時間間の差分関係を示す差分方程式の形式を有し、一定時間において不変な定数とみなせる定数行列を指数とする行列指数関数を近似する行列冪級数であって所定の計算ステップと前記定数行列との積の冪乗項により構成される有限な次数の第１の行列冪級数を夫々含む前記第１の係数行列及び前記第２の係数行列が夫々定義される、状態方程式を設定する処理と、
　前記計算ステップを時間的に不変な時間不変計算ステップと時間的に可変な時間可変計算ステップとに分離し、少なくとも、前記時間不変計算ステップと前記定数行列との積の冪乗項により構成される前記有限な次数の第２の行列冪級数を前記状態方程式の解法演算に先立って事前に計算する処理と、
　前記事前に計算された前記第２の行列冪級数に基づいて夫々前記第１の係数行列及び前記第２の係数行列を算出しながら、順次算出される時間可変計算ステップに基づいて、前記差分方程式に対する解法演算を実行する処理と、
　を実行させるためのプログラムを記憶した記憶媒体。