(57)【特許請求の範囲】
【請求項１】
不規則外乱を受ける物理系の状態量ｘの時間変化を、関数と乱数によって記述し、時間刻みｈ毎に、状態量ｘの初期値から、状態量ｘの近似値を陽的に求め、シミュレーション結果を出力するシミュレーション装置であって、
関数ｆ（ｘ）およびｇ（ｘ）と時刻ｔ０における状態量ｘの
0003735128.tif 000042
を入力する初期設定入力手段と、
0003735128.tif 000043
を格納する第１のパラメータ格納手段と、
0003735128.tif 000044
を格納する第２のパラメータ格納手段と、
時刻ｔｎにおける状態量ｘの
0003735128.tif 000045
を格納する近似値格納手段と、
ｎ番目の時刻ｔｎ＝ｔ０＋ｎｈにおける第１の乱数ξ、時刻ｔｎにおける
0003735128.tif 000046
、時間刻みｈ、正の実数ν、関数ｆ（ｘ）、関数ｇ（ｘ）、および前記ａパラメータを用いて、状態量ｘの
0003735128.tif 000047
から
0003735128.tif 000048
を計算するための
0003735128.tif 000049
というｋパラメータ計算式を格納する第１の計算式格納手段と、
前記ｋパラメータ、時間刻みｈ、正の実数ν、ΔＷ、
0003735128.tif 000050
、および前記ｂパラメータを用いて状態量ｘの時間増分Ψを計算するための
0003735128.tif 000051
という時間増分計算式を格納する第２の計算式格納手段と、
時間増分Ψを用いて状態量ｘの
0003735128.tif 000052
から時刻ｔｎ＋１における状態量ｘの
0003735128.tif 000053
を計算するための
0003735128.tif 000054
という近似値計算式を格納する第３の計算式格納手段と、
第１の乱数ξ、
0003735128.tif 000055
、時間刻みｈ、正の実数ν、関数ｆ（ｘ）、関数ｇ（ｘ）、前記ａパラメータ、および前記ｋパラメータ計算式を用いて、前記ｋパラメータを生成するｋパラメータ生成手段と、
時間増分Ψの値を格納する時間増分格納手段と、
前記時間増分格納手段に０を設定した後、時間刻みｈ、ｂパラメータｂｉ、ｋパラメータｋｉ、および前記時間増分計算式を用いて該時間増分計算式の第１項の値を生成して該時間増分格納手段に加算し、ΔＷ、
0003735128.tif 000056
、および該時間増分計算式を用いて該時間増分計算式の第２項の値を生成して該時間増分格納手段に加算し、
0003735128.tif 000057
、および該時間増分計算式を用いて該時間増分計算式の第３項の値を生成して該時間増分格納手段に加算し、時間刻みｈ、正の実数ν、
0003735128.tif 000058
、および該時間増分計算式を用いて該時間増分計算式の第４項の値を生成して該時間増分格納手段に加算することで、時間増分Ψを生成する増分生成手段と、
生成された時間増分Ψと前記近似値計算式を用いて状態量ｘの
0003735128.tif 000059
を生成し、前記近似値格納手段の
0003735128.tif 000060
を
0003735128.tif 000061
に更新する近似値更新手段と、
前記近似値更新手段から与えられる時間刻みｈ毎の状態量ｘの各近似値を、前記物理系に対する前記シミュレーション結果として出力する出力手段と
を有することを特徴とするシミュレーション装置。
【請求項２】
前記近似値更新手段は更新された近似値を前記増分生成手段に入力し、
前記増分生成手段は前記近似値更新手段より与えられる前記更新された近似値から、前記時間増分を生成することを特徴とする請求項１記載のシミュレーション装置。
【請求項３】
前記増分生成手段は前記第１の乱数を発生する第１の乱数発生手段と、前記第１の乱数とは独立な前記第２の乱数を発生する第２の乱数発生手段を有することを特徴とする請求項１または２記載のシミュレーション装置。
【請求項４】
前記第１の乱数と前記第２の乱数は、互いに独立な、平均０、分散１の正規乱数の系列に属することを特徴とする請求項３記載のシミュレーション装置。
【請求項５】
前記第１のパラメータ格納手段および第２のパラメータ格納手段は、前記ａパラメータおよびｂパラメータとして、前記時間増分を前記時間刻みの巾について展開した時の、前記時間刻みについて１．５次までの各項の係数が、前記物理系を記述する前記状態量に関する微分方程式の解を前記時間刻みの巾について展開した時の、前記時間刻みについて１．５次までの各項の係数に一致するような値を格納することを特徴とする請求項１または２記載のシミュレーション装置。
【請求項６】
前記第１のパラメータ格納手段および第２のパラメータ格納手段は、前記ａパラメータおよびｂパラメータとして、
0003735128.tif 000062
0003735128.tif 000063
0003735128.tif 000064
0003735128.tif 000065
0003735128.tif 000066
という条件を満たす値を格納することを特徴とする請求項５記載のシミュレーション装置。
【請求項７】
不規則外乱を受ける物理系の状態量ｘの時間変化を、関数と乱数によって記述し、時間刻みｈ毎に、状態量ｘの初期値から、状態量ｘの近似値を陽的に求め、シミュレーション結果を出力するシミュレーション方法であって、
初期設定入力手段が、関数ｆ（ｘ）およびｇ（ｘ）と時刻ｔ０における状態量ｘの
0003735128.tif 000067
を入力し、
ｋパラメータ生成手段が、ｎ番目の時刻ｔｎ＝ｔ０＋ｎｈにおける第１の乱数ξ、時刻ｔｎにおける
0003735128.tif 000068
、時間刻みｈ、正の実数ν、関数ｆ（ｘ）、関数ｇ（ｘ）、第１のパラメータ格納手段に格納された
0003735128.tif 000069
、および第１の計算式格納手段に格納された、
0003735128.tif 000070
というｋパラメータ計算式を用いて、近似値格納手段に格納された時刻ｔｎにおける状態量ｘの
0003735128.tif 000071
から
0003735128.tif 000072
を生成し、
増分生成手段が、状態量ｘの時間増分Ψの値を格納する時間増分格納手段に０を設定した後、時間刻みｈ、第２のパラメータ格納手段に格納されたｂパラメータｂｉ（ｉ＝１，．．．，ｓ）、生成されたｋパラメータｋｉ、および第２の計算式格納手段に格納された、
0003735128.tif 000073
という時間増分計算式を用いて、該時間増分計算式の第１項の値を生成して該時間増分格納手段に加算し、ΔＷ、前記第２のパラメータ格納手段に格納された
0003735128.tif 000074
、生成された
0003735128.tif 000075
、および該時間増分計算式を用いて該時間増分計算式の第２項の値を生成して該時間増分格納手段に加算し、
0003735128.tif 000076
、前記第２のパラメータ格納手段に格納された
0003735128.tif 000077
、生成された
0003735128.tif 000078
、および該時間増分計算式を用いて該時間増分計算式の第３項の値を生成して該時間増分格納手段に加算し、時間刻みｈ、正の実数ν、前記第２のパラメータ格納手段に格納された
0003735128.tif 000079
、生成された
0003735128.tif 000080
、および該時間増分計算式を用いて該時間増分計算式の第４項の値を生成して該時間増分格納手段に加算することで、時間増分Ψを生成し、
近似値更新手段が、生成された時間増分Ψと第３の計算式格納手段に格納された、
0003735128.tif 000081
という近似値計算式を用いて、状態量ｘの
0003735128.tif 000082
から時刻ｔｎ＋１における状態量ｘの
0003735128.tif 000083
を生成して、前記近似値格納手段の
0003735128.tif 000084
を
0003735128.tif 000085
に更新し、
出力手段が、前記近似値更新手段から与えられる時間刻みｈ毎の状態量ｘの各近似値を、前記物理系に対する前記シミュレーション結果として出力する
ことを特徴とするシミュレーション方法。
【発明の詳細な説明】
【０００１】
【産業上の利用分野】
本発明は不規則外乱を受ける系の数値シミュレーションを行うシミュレーション装置に関する。
【０００２】
【従来の技術】
自然科学の分野において、一般に物理系を記述する数学モデルの解を解析的に求めるのは難しく、その数学モデルの初期値問題は数値処理により、近似的に解かれる場合が多い。近年の情報処理装置の発達を背景にして、こうした数値処理を行う種々のシミュレーション装置が考案されているが、シミュレーション精度を上げるためには膨大な計算量が必要となり、汎用性が損なわれることがある。そこでシミュレーションの信頼性と汎用性を考慮しつつ、いかにシミュレーション装置を構築するかが重要な問題である。
【０００３】
物理系として、工業プラントやロボットなどの挙動を扱う工学の分野においては、センサからの観測データをもとに物理系の挙動を推定し、この物理系を制御することが重要な課題である。この推定処理は一般にシミュレーション装置を用いた数値シミュレーションにより行われる。
【０００４】
シミュレーション装置の性能は数値シミュレーションの結果の信頼性を左右し、この結果を用いて物理系を制御する制御装置の性能を決定付けることになる。このため、数値シミュレーションの処理対象である物理系や、制御装置に要求される信頼性に応じて、多種多様なシミュレーション装置が開発されている。
【０００５】
処理対象となる物理系のモデルは自然法則に基づいて導出されるが、その際に、現実の多くの物理系に混在している雑音などの不規則外乱の影響が無視され、理想化あるいは単純化が施されている場合が多い。しかし対象とする物理系によっては、不規則外乱の影響が無視できないものもある。
【０００６】
例えば人工衛星の制御装置を設計する際には、太陽の放射圧、重力、宇宙空間に浮遊する他の衛星の残骸との衝突などの不規則外乱を無視することができない。また工業プラントやロボットなどの制御装置においては、制御対象を観測するセンサからの信号に混入する雑音などの不規則外乱が無視できないことがある。
【０００７】
このように、不規則外乱の影響を受ける物理系の解析や制御を目的とするシミュレーションにおいては、不規則外乱を反映した物理モデルを与え、これを数値処理する必要がある。
【０００８】
一般に任意の物理系のモデルは次のような常微分方程式により表現できる。
【０００９】
【数１】
0003735128.tif 000002
【００１０】
ここで、ｘ＝ｘ（ｔ）は物理系における変位や温度等の、シミュレーションが対象とする物理量であり、時刻ｔに依存して変動する。  外１  はｘのｔに関する
【００１１】
【外１】
0003735128.tif 000003
【００１２】
微分を表す。ｆ（ｘ，ｔ）はｘとｔの関数であり、ｘ0 は時刻ｔ0 において与えられるｘの初期値である。
（１）式の物理モデルに対する数値シミュレーションにおいては、時刻ｔn ＝ｔ0 ＋ｎｈにおける物理量ｘ（ｔn ）の近似値をｘn とし、次式により次の時刻ｔn+1 における物理量ｘ（ｔn+1 ）の近似値ｘn+1 を計算する。
【００１３】
【数２】
0003735128.tif 000004
【００１４】
ただし、ｘ0 は近似値ではなく、時刻ｔ0 におけるｘの初期値である。ここでｈはシミュレーションの時間ステップの刻み幅であり、Ψ（ｘn ，ｔn ，ｈ）は時刻ｔn における近似値ｘn の時間増分を与えている。このΨ（ｘn ，ｔn ，ｈ）をｘn 、ｔn 、ｈを用いてどのように構築するかが重要であり、数値シミュレーションの結果の信頼性を左右することになる。
【００１５】
（２）式のΨ（ｘn ，ｔn ，ｈ）の構築方法はテーラ級数項の方法とルンゲ−クッタ法の２種類に大別される。このうちテーラ級数項の方法は原理的には単純で、Ψ（ｘn ，ｔn ，ｈ）をシミュレーションに必要な次数までのテーラ級数項により与える方法である。しかしこの方法では一般に、ｐ次の項の係数に関数ｆ（ｘ，ｔ）のｐ次の導関数∂p ｆ（ｘ，ｔ）／∂ｘp が含まれる。関数ｆ（ｘ，ｔ）がｘの線型関数等の簡単な関数形で与えられる場合は、導関数は簡単に求められるが、ｆ（ｘ，ｔ）が複雑な非線型関数の場合、その高階の導関数をｘ、ｔの関数として与えることは困難である。
【００１６】
またこのテーラ級数項の方法ではオペレータがｆ（ｘ，ｔ）の導関数を求めて、シミュレーション装置に入力しなければならず、さらに導関数の値の計算はシミュレーション装置にとって多大な負担となるので汎用性に乏しい。
【００１７】
そこで数値シミュレーションでは、もう一つの方法として、次式によりΨ（ｘn ，ｔn ，ｈ）を構築するルンゲ−クッタ法が多く用いられている。
【００１８】
【数３】
0003735128.tif 000005
【００１９】
ここでｓは自然数であり、ｋパラメータｋi は適切に選ばれたｘ、ｔに対応する関数ｆ（ｘ，ｔ）の値で、次式のような計算処理により生成される。
【００２０】
【数４】
0003735128.tif 000006
【００２１】
（３）、（４）式でｂパラメータｂi やａパラメータａijはｋパラメータに対する重みを表すパラメータであり、シミュレーション装置内のメモリに予め格納されている。ｂパラメータとａパラメータの選び方によって近似値ｘn に対する増分Ψ（ｘn ，ｔn ，ｈ）の値が異なり、近似値ｘn+1 の精度が左右されるため、結果としてシミュレーション装置の性能が決められる。
【００２２】
このｋパラメータを用いた方法は、（２）式の数値処理にあたって係数関数ｆ（ｘ，ｔ）の微分形を与える必要がないため、テーラ級数の方法よりも遙かに汎用的である。
【００２３】
時刻ｔn における（１）式の解をｘ（ｔn ）とすると、数値シミュレーションによる近似値ｘn のｘ（ｔn ）に対する誤差はｘ（ｔn ）−ｘn で与えられる。シミュレーション結果の信頼性は誤差ｘ（ｔn ）−ｘn が小さいほど高くなる。この誤差の大きさを評価する尺度として、刻み幅ｈが用いられる。誤差がｈr+1 程度であることが保証されるとき、そのシミュレーション装置あるいはシミュレーションの結果はｒ次の近似精度を持つといわれ、ｒの値が大きいほど、信頼性の高いシミュレーション結果が得られる。
【００２４】
対象とする物理系が雑音等の不規則外乱の影響を受けている場合は、この系の物理モデルは、（１）式の代わりに次のような伊藤型確率微分方程式を用いて記述される。
【００２５】
【数５】
0003735128.tif 000007
【００２６】
ここで、ｘ＝ｘ（ｔ）は対象とする物理系における状態量であり、ｆ（ｘ）、ｇ（ｘ）はｘの関数である。またｗ（ｔ）は雑音モデルを表し、その値はブラウン運動モデルの確率過程に依存する。このモデルは（１）式の右辺に、不規則外乱による第２項を付加した形になっている。ｘ0 は時刻ｔ0 において与えられるｘの初期値である。
【００２７】
不規則外乱を受ける物理モデルを表す（５）式は、（１）式と異なって不規則外乱による付加項があるため、その数値シミュレーションを行うためには、（３）、（４）式のΨ（ｘn ，ｔn ，ｈ）をそのまま用いることはできない。この場合には、確率過程に依存する不規則外乱の項を扱うために、（３）、（４）式のΨ（ｘn ，ｔn ，ｈ）を確率微分方程式の場合に拡張する必要がある。
【００２８】
この拡張された時間増分Ψ（ｘn ，ｔn ，ｈ）をｋパラメータを用いて構築する従来の方法として、確率的Heun法、１．５次精度を保証する３段陽的解法、漸近的に有効な解法がある。
【００２９】
確率的Heun法（RUMELIN, SIAM J. Vol.19, No.3, June 1982）は（５）式の物理モデルに対するΨ（ｘn ，ｔn ，ｈ）として、（３）式の右辺に確率過程に依存する項を付加したものを用いている。この方法は確率過程を扱う不規則外乱を受ける物理モデルに対する（３）式の素直な拡張であり、１次の近似精度を保証している。また（３）式と異なってｋパラメータが関数ｇ（ｘ）とｇ（ｘ）の１階の導関数の値を含むことを特徴とする。
【００３０】
１．５次精度を保証する３段陽的解法（斉藤、三井  日本応用数理学会論文誌 Vol.2, No.1, 1992）は２種類の確率過程を用いて、確率的Heun法のΨ（ｘn ，ｔn ，ｈ）をさらに拡張している。この方法は一般に１．５次の高い近似精度を保証しており、ｋパラメータは関数ｆ（ｘ）とｇ（ｘ）のそれぞれの１階の導関数ならびに、ｇ（ｘ）の２階の導関数の値を含んでいる。
【００３１】
漸近的に有効な解法（NEWTON, SIAM J. Vol.51, No.2, April 1991）は確率的Heun法のΨ（ｘn ，ｔn ，ｈ）にさらに（ｈ）1/2 に関する項を付加して、ｋパラメータの計算処理を簡単化している。この方法では、ｋパラメータが関数ｇ（ｘ）の１階の導関数の値を含まず、関数ｆ（ｘ）とｇ（ｘ）の値を用いて計算できるので、汎用性に優れたシミュレーション方法を与えている。この方法により保証される近似精度は確率的Heun法と同様に１次である。
【００３２】
【発明が解決しようとする課題】
しかしながら上述のような従来の時間増分の構築方法による数値シミュレーションには以下のような問題がある。
【００３３】
確率的Heun法では、近似値ｘn に対する時間増分Ψ（ｘn ，ｔn ，ｈ）を計算する時、既に計算されているｘn の値とｔn 、ｈを用いて２組のｋパラメータを計算し、ｋパラメータに適当な重みを付けて加算する。このとき一方の組のｋパラメータの計算式が、不規則外乱を記述している関数ｇ（ｘ）の１階の導関数を含むため、数値シミュレーションを行う際に、関数ｇ（ｘ）の１階の導関数の値を計算しなければならない。このためテーラ級数項の方法と同様に計算量が増大して、シミュレーション装置の負荷が大きくなるという欠点をもち、汎用的でない。また保証されるシミュレーションの近似精度は１次に留まっており、必ずしもその信頼性は高くない。
【００３４】
１．５次精度を保証する３段陽的解法は、確率的Heun法よりもシミュレーション結果の信頼性は高いが、ｋパラメータの計算式が関数ｆ（ｘ）とｇ（ｘ）の１階の導関数ならびに、ｇ（ｘ）の２階の導関数を含むため、これらの導関数の関数形を与える必要がある。関数ｇ（ｘ）だけでなく、関数ｆ（ｘ）の導関数の値が必要な上、さらに関数ｇ（ｘ）の２階の導関数の値も必要とするため、確率的Heun法よりも著しく計算量が増大し、実際にはシミュレーション装置として実用化するのは困難である。
【００３５】
漸近的に有効な解法は、確率的Heun法と１．５次精度を保証する３段陽的解法と比べて、数値シミュレーションを行う際に、関数ｇ（ｘ）等の導関数を与える必要がないという利点を有する。この方法の場合、ｋパラメータの計算式は関数ｆ（ｘ）とｇ（ｘ）により記述されており、数値処理が容易で、汎用性に優れている。しかし保証されるシミュレーションの近似精度は１次に留まっており、その信頼性は低いという問題がある。
【００３６】
不規則外乱の影響を受ける系のシミュレーションを行うには、（５）式に示される確率微分方程式の解を数値的に求めるシミュレーション装置の開発が不可欠である。それにも係わらず、従来のシミュレーション方法が信頼性と汎用性の面で上述のような問題を抱えているため、実用的なシミュレーション装置はほとんど開発されていない。
【００３７】
本発明は、不規則外乱を受ける系の数値シミュレーションにおいて、計算量が少なくて汎用性に優れ、かつ、高い近似精度を持つシミュレーション結果が得られるシミュレーション装置を提供することを目的とする。
【００３８】
【課題を解決するための手段】
図１は本発明の原理図である。
本発明は、不規則外乱を受ける物理系を、状態量の時間変化を表す関数と乱数によって記述し、時間刻み毎に、状態量の初期値から、状態量の近似値を陽的に求め、シミュレーション結果を出力する情報処理装置におけるシミュレーション装置１０である。
【００３９】
本シミュレーション装置１０は、初期設定入力手段１１と、増分生成手段１２、近似値更新手段１３、出力手段１４を有する。
初期設定入力手段１１は、初期データである上記関数と初期値を増分生成手段１２に入力する。
【００４０】
増分生成手段１２は、上記乱数に対応する第１の乱数と第２の乱数を発生し、これらの乱数と時間刻みを用いて、初期設定入力手段１１より与えられる関数の関数値として、３種類以上の異なるｋパラメータを状態量の初期値から生成する。そして、時間刻みの平方根に比例し、かつ、第１の乱数と第２の乱数に依存しない項を用いて、生成したｋパラメータからルンゲ−クッタ型の状態量の時間増分を生成する。
【００４１】
近似値更新手段１３は、増分生成手段１２が生成した時間増分を用いて状態量の近似値を更新し、更新した近似値を出力手段１４と増分生成手段１２に与える。
【００４２】
増分生成手段１２は、上記の処理に従い、さらに近似値更新手段１３から与えられる更新された近似値に対する時間増分を生成し、近似値更新手段１３に与える。
【００４３】
出力手段１４は、近似値更新手段１３から与えられる時間刻み毎の状態量の各近似値を、シミュレーション結果として外部に出力する。
上記増分生成手段１２は、さらに数値処理を簡単化するために、不図示のパラメータ生成手段と時間増分生成手段、第１のパラメータ格納手段、第２のパラメータ格納手段を備えるものとする。
【００４４】
第１のパラメータ格納手段は、パラメータ生成手段が用いる第１のパラメータを格納し、第２のパラメータ格納手段は、時間増分生成手段が用いる第２のパラメータを格納する。
【００４５】
パラメータ生成手段は、時間刻み、第１の乱数、第２の乱数のうちの一つと、第１のパラメータおよび状態量の近似値とを用いて、状態量に対応する標本状態量を生成し、与えられた関数を用いて、生成した標本状態量から上記ｋパラメータを生成し、時間増分生成手段に与える。
【００４６】
時間増分生成手段は、パラメータ生成手段から与えられるｋパラメータと、第１の乱数、第２の乱数、時間刻みおよび第２のパラメータとを用いて、状態量の時間増分を生成する。
【００４７】
【作用】
増分生成手段１２が、２種類の乱数を用いて時間増分を生成するので、時間増分を時間刻みの巾について展開した時に、１種類の乱数のみを用いて時間増分を生成するシミュレーション装置よりも、高い次数の項を含むことが可能になる。
【００４８】
また時間増分が、時間刻みの平方根に比例し、かつ、上記２種類の乱数の値に依存しない項を含むため、（５）式における関数ｆ（ｘ）やｇ（ｘ）の導関数の値を計算する必要がなくなり、このような項を用いないシミュレーション装置よりも数値処理が容易になる。
【００４９】
また状態量の近似値から標本状態量を生成し、生成した標本状態量からｋパラメータを生成するパラメータ生成手段と、ｋパラメータから時間増分を生成する時間増分生成手段を備えることにより、数値処理を段階的に行うことができる。
【００５０】
【実施例】
以下、図面を参照しながら、本発明の実施例を説明する。
本実施例においては、不規則外乱の影響を受ける系として、（５）式の形で記述される物理モデルのシミュレーションを行う場合を考える。ここで状態量ｘ（ｔ）は一般にｄ次元の状態量であり、（５）式はｄ次元の確率微分方程式を記述する数学モデルに対応する。
【００５１】
本実施例のシミュレーション装置では、時刻ｔn ＝ｔ0 ＋ｎｈにおける状態量ｘ（ｔ）の近似値  外２  （以後ｘn のバーと記す。）と、時刻ｔn+1 ＝ｔ0 ＋
【００５２】
【外２】
0003735128.tif 000008
【００５３】
（ｎ＋１）ｈにおける状態量ｘ（ｔ）の近似値  外３  （以後ｘn+1 のバーと記
【００５４】
【外３】
0003735128.tif 000009
【００５５】
す。）の関係を次式により記述する。
【００５６】
【数６】
0003735128.tif 000010
【００５７】
Ψ（ｘn のバー，ｔn ，ｈ）は時刻ｔn における状態量ｘ（ｔ）の近似値ｘn のバーに対する時間増分を表している。近似値ｘn のバー、ｘn+1 のバー、時間増分Ψ（ｘn のバー，ｔn ，ｈ）もまたｄ次元の値を持つ。
【００５８】
図２は本実施例のシミュレーション装置の構成図である。図２において初期設定入力部２１は図１の初期設定入力手段１１に相当し、外部から与えられる初期データである（５）式における関数ｆ（ｘ）とｇ（ｘ）の具体的な関数形と、時刻ｔ0 における状態量ｘの初期値ｘ0 ＝ｘ（ｔ0 ）とを初期設定としてｋパラメータ生成装置２２に入力する。
【００５９】
ｋパラメータ生成装置２２、時間増分生成装置２４、メモリ２３─１〜２３─４、メモリ２５─１〜２５─４は、それぞれパラメータ生成手段、時間増分生成手段、第１のパラメータ格納手段、第２のパラメータ格納手段に相当する。そして、これらは乱数発生装置２６─１、２６─２と共に、図１の増分生成手段１２を形成している。
【００６０】
乱数発生装置２６─１、２６─２はそれぞれ、第１の乱数、第２の乱数を発生する。
ｋパラメータ生成装置２２と時間増分生成装置２４は、初期設定入力手段１１または近似値更新装置２７からの入力に基づいて、時刻ｔn における状態量ｘ（ｔ）の近似値ｘn のバーに対する時間増分Ψ（ｘn のバー，ｔn ，ｈ）を生成し、近似値更新装置２７に入力する。
【００６１】
近似値更新装置２７は図１の近似値更新手段１３に相当し、時間増分生成装置２４から入力される時間増分Ψ（ｘn のバー，ｔn ，ｈ）を近似値ｘn のバーに加算してこれをｘn+1 のバーとし、ｋパラメータ生成装置２２と出力部２８に与える。
【００６２】
出力部２８は図１の出力手段１４に相当し、近似値更新装置２７から逐次与えられる状態量の近似値をまとめてシミュレーション結果とし、外部に出力する。図２のシミュレーション装置では、（５）式の右辺第２項の不規則外乱を記述する項ｇ（ｘ）ｗ（ｔ）を数値処理して、近似値ｘn のバーに反映させるため、４種類のｋパラメータを導入し、時間増分Ψ（ｘn のバー，ｔn ，ｈ）を次式に基づいて生成する。
【００６３】
【数７】
0003735128.tif 000011
【００６４】
（７）式の右辺第１項は常微分方程式に関するルンゲ−クッタ法の場合に類似した項であり、第３、第４項は（５）式の解を１．５次以上の近似精度で求めるために必要な項である。ここでパラメータΔＷn+1 、  外４  （以後ΔＷn+1 の
【００６５】
【外４】
0003735128.tif 000012
【００６６】
波と記す。）はそれぞれ、時刻ｔn における平均０、分散１の正規乱数のサンプルξn+1 、  外５  （以後ξn+1 の波と記す。）を用いてΔＷn+1 ＝ξn+1 （ｈ
【００６７】
【外５】
0003735128.tif 000013
【００６８】
）1/2 、ΔＷn+1 の波＝ξn+1 の波（ｈ）1/2 のように与えられる。サンプルξn+1 、ξn+1 の波は、それぞれ第１の乱数、第２の乱数に相当し、近似値ｘn のバーを更新する度に、乱数発生装置２６─１、２６─２より独立に生成される。
【００６９】
（７）式の右辺第４項は、サンプルξn+1 、ξn+1 の波の値に依存せず、（ｈ）1/2 に比例する。第３、第４項があるために、関数ｆやｇの導関数の値を用いなくても、高い近似精度を実現できる。従ってこれらの項は、数値処理を容易にし、シミュレーション装置の汎用性を高めるために必要である。ここでνは適切な正の実数である。
【００７０】
（７）式においてｂパラメータｂi 、  外６  （以後ｂi のバーと記す。）、
【００７１】
【外６】
0003735128.tif 000014
【００７２】
外７  （以後ｂi の波と記す。）、  外８  （以後ｂi の山と記す。）、（ｉ
【００７３】
【外７】
0003735128.tif 000015
【００７４】
【外８】
0003735128.tif 000016
【００７５】
＝１，・・・，ｓ）はそれぞれ、時間増分を生成する際のｋパラメータｋi 、  外９  （以後ｋi のバーと記す。）、  外１０  （以後ｋi の波と記す。）、
【００７６】
【外９】
0003735128.tif 000017
【００７７】
【外１０】
0003735128.tif 000018
【００７８】
外１１  （以後ｋi の山と記す。）、（ｉ＝１，・・・，ｓ）に対する重み付け
【００７９】
【外１１】
0003735128.tif 000019
【００８０】
のパラメータで、予めメモリ２５─１、２５─２、２５─３、２５─４、に格納されている。これらのｂパラメータは第２のパラメータに相当する。
ｋパラメータｋi （ｉ＝１，・・・，ｓ）は近似値ｘn のバーに対して、（８）式により与えられる関数ｆ（ｘ）の値で、ｋパラメータｋi のバー、ｋi の波、ｋi の山（ｉ＝１，・・・，ｓ）はそれぞれ、（９）、（１０）、（１１）式により与えられる関数ｇ（ｘ）の値である。これらの４種類のｋパラメータはｋパラメータ生成装置２２により生成される。
【００８１】
（８）、（９）、（１０）、（１１）式におけるａパラメータａij、  外１２
【００８２】
【外１２】
0003735128.tif 000020
【００８３】
（以後ａijのバーと記す。）、  外１３  （以後ａijの波と記す。）、  外１
【００８４】
【外１３】
0003735128.tif 000021
【００８５】
４  （以後ａijの山と記す。）、（ｉ＝１，・・・，ｓ）はそれぞれ、ｋパラメ
【００８６】
【外１４】
0003735128.tif 000022
【００８７】
ータを生成する際のｋパラメータｋi 、ｋi のバー、ｋi の波、ｋi の山（ｉ＝１，・・・，ｓ）に対する重み付けのパラメータで、予めメモリ２３─１、２３─２、２３─３、２３─４、に格納されている。これらのａパラメータは第１のパラメータに相当する。
【００８８】
次にこれらのａパラメータとｂパラメータの決定方法について説明する。
一般にルンゲ−クッタ法による微分方程式の近似解ｘn のバーに対して、ｒ次の近似精度を保証するためには、時刻ｔn+1 における解ｘ（ｔn+1 ）を、時刻ｔn における解ｘ（ｔn ）のまわりでｈに関してテーラ展開し、ｈのｒ次の次数の項までが一致するように時間増分Ψ（ｘn のバー，ｔn ，ｈ）を決める必要がある。
【００８９】
本シミュレーション装置においては、（５）式の近似解ｘn のバーに対して、１．５次の近似精度を保証するために、ｈに関して１．５次の項までのテーラ展開を用いることにする。時刻ｔn+1 における（５）式の解ｘ（ｔn+1 ）のｘ（ｔn ）のまわりでのテーラ展開は、次式で記述されることが知られている（Ito, Applied Math. and Optimization, Vol.1, 1975）。
【００９０】
【数８】
0003735128.tif 000023
【００９１】
（１２）式におけるＦとＧはｄ次元の微分作用素で、次式により定義される。
【００９２】
【数９】
0003735128.tif 000024
【００９３】
また（１２）式のΔＷn+1 、ΔＷn+1 の波はそれぞれ独立なブラウン運動モデルの増分を表し、（７）式のΔＷn+1 、ΔＷn+1 の波と同様に正規乱数のサンプルξn+1 、ξn+1 の波を用いてΔＷn+1 ＝ξn+1 （ｈ）1/2 、ΔＷn+1 の波＝ξn+1 の波（ｈ）1/2 のように定められる。Ｏ（ｈ2 ）はｈに関して２次以上の次数の項を表す。
【００９４】
ここで時刻ｔn+1 における近似値ｘn+1 のバーを与える（６）式の右辺を、近似値ｘn のバーのまわりでｈに関してテーラ展開し、ｈに関して１．５次の項までの各項の係数を（１２）の各項の係数と比較すれば、ａパラメータとｂパラメータに対する拘束条件を得ることができる。このとき、（７）式右辺の第３、第４項があるために、関数ｆ、ｇの導関数を用いることなく、（１２）式右辺のＯ（ｈ2 ）以外の項を得ることができる。こうして求められた拘束条件を次に示す。
【００９５】
【数１０】
0003735128.tif 000025
【００９６】
【数１１】
0003735128.tif 000026
【００９７】
ただし、ｃi 、  外１５  （以後ｃi のバーと記す。）、  外１６  （以後ｃi 
【００９８】
【外１５】
0003735128.tif 000027
【００９９】
【外１６】
0003735128.tif 000028
【０１００】
の波と記す。）、  外１７  （以後ｃi の山と記す。）は各ａパラメータを用い
【０１０１】
【外１７】
0003735128.tif 000029
【０１０２】
て次式により定義する。
【０１０３】
【数１２】
0003735128.tif 000030
【０１０４】
ａパラメータとｂパラメータに対する拘束条件が（１４．１）〜（１４．２２）式の２２個であることを考慮して、ｓ＝４と選べば、シミュレーションに必要なａパラメータとｂパラメータの総数は４０個になる。ここで、簡単のためｂi ＝ｂi のバー  、ａij＝ａijのバー、ｃi ＝ｃi の波、ｂ3 の波＝ｂ3 の山、ｂ4 の波＝ｂ4 の山と置くと、パラメータの総数は２５個になる。またこのとき、（１４．１）、（１４．４）の２式が同値となり、（１４．２）、（１４．９）、（１４．１０）、（１４．１１）の４式が同値となり、（１４．１９）、（１４．２１）の２式が同値となる。この結果２２個の拘束条件は１７個に縮退する。従ってａパラメータとｂパラメータを決定するにあたっては、８個のパラメータを自由パラメータとして選び、さらにνの値を適当に与える必要がある。
【０１０５】
一例としてａ31の波＝ａ41の波＝ａ42の波＝０、ａ32＝１／４、ａ43の山＝１／６、ｃ2 ＝１／２、ｃ4 ＝１、ｃ4 の山＝０、ν＝３と選び、これらを１７個の拘束条件を表す代数方程式に代入して、他の独立パラメータを求めると、各メモリに格納される４０個のａパラメータ、ｂパラメータが次のように決定される。
【０１０６】
【数１３】
0003735128.tif 000031
【０１０７】
【数１４】
0003735128.tif 000032
【０１０８】
（１６．１）式のｂパラメータの値およびν＝３を（７）式に代入すると、時間増分を与える式は次式のようになる。
【０１０９】
【数１５】
0003735128.tif 000033
【０１１０】
また（１６．２）式のａパラメータの値およびν＝３を（８）〜（１１）式に代入すると、ｋパラメータを与える式は次式のようになる。
【０１１１】
【数１６】
0003735128.tif 000034
【０１１２】
【数１７】
0003735128.tif 000035
【０１１３】
（６）、（１７）、（１８．１）、（１８．２）式により生成される近似値ｘn+1 のバーは、ｈに関して１．５次の次数の項まで、（１２）式の右辺のテーラ展開と一致するため、１．５次の近似精度が保証されている。また関数ｆやｇの導関数を必要としないため、オペレータがこれらの導関数の関数形を求める必要がなく、近似値の数値計算における計算量も少なくて済む。
【０１１４】
図３は図２のシミュレーション装置の動作の概要を示したフローチャートである。
シミュレーションの考察対象の物理系が与えられると、オペレータは前処理Ｐ１において対象とする物理系を記述する物理モデルとその初期状態を決める。この物理モデルと初期状態は（５）式に対応している。そしてここで決めた関数ｆ（ｘ）、ｇ（ｘ）の関数形と、状態量ｘの初期値ｘ0 ＝ｘ（ｔ0 ）を初期データとしてシミュレーション装置に入力する。
【０１１５】
次にシミュレーション装置の内部では、初期設定入力部２１が与えられた初期データをｋパラメータ生成装置２２に入力する。ただし、このときｘ0 のバー＝ｘ0 と設定する。
【０１１６】
ステップＱ１でｋパラメータ生成装置２２は、関数ｆ（ｘ）、ｇ（ｘ）の関数形を格納した後、メモリ２３─１〜２３─４に格納されたａパラメータと乱数発生装置２６─１、２６─２が発生する正規乱数のサンプルξn+1 、ξn+1 の波を用いて、初期値ｘ0 より各ｋパラメータを逐次生成する。このｋパラメータ生成処理はｋパラメータ生成装置２２内に格納された（８）〜（１１）式を用いて行う。ただし、このときｎ＝０、ｓ＝４とする。
【０１１７】
ステップＱ２で時間増分生成装置２４は、メモリ２５─１〜２５─４に格納されたｂパラメータと乱数発生装置２６─１、２６─２が発生する正規乱数のサンプルξn+1 、ξn+1 の波を用いて、生成されたｋパラメータから時間増分Ψ（ｘn のバー，ｔn ，ｈ）を生成する。この時間増分生成処理は時間増分生成装置２４内に格納された（７）式を用いて行う。ただし、このときｎ＝０、ｓ＝４、ν＝３とする。
【０１１８】
ステップＱ３で近似値更新装置２７は、格納している（６）式を用いて、生成された時間増分Ψ（ｘn のバー，ｔn ，ｈ）とｘn のバーから、次のシミュレーションの参照時刻ｔn+1 における近似値ｘn+1 のバーを計算し、これをｋパラメータ生成装置２２と出力部２８に出力する。ただし、このときｎ＝０とする。
【０１１９】
ステップＱ４で近似値更新装置２７は、ｎがシミュレーションの参照終了時刻に対応するかどうかを判定する。
まだ参照終了点の近似値が得られていない場合はｎ＝ｎ＋１とし、ステップＱ１に戻って新たにｋパラメータを生成し、ステップＱ２で時間増分Ψ（ｘn のバー，ｔn ，ｈ）を生成し、ステップＳ３で次の参照点の近似値ｘn+1 のバーを求める。以後ステップＱ１〜Ｑ４を繰り返し、ステップＱ４で参照終了時刻に対応した最終参照点に到達したことを確認すると処理を終える。
【０１２０】
図４〜７は図３のステップＱ１のｋパラメータ生成装置２２の詳細な動作を示すフローチャートである。図５〜７におけるｘn (m) 、  外１８  （以後ｘn (m
【０１２１】
【外１８】
0003735128.tif 000036
【０１２２】
) のバーと記す。）、  外１９  （以後ｘn (m) の波と記す。）、  外２０  （
【０１２３】
【外１９】
0003735128.tif 000037
【０１２４】
【外２０】
0003735128.tif 000038
【０１２５】
以後ｘn (m) の山と記す。）、（ｍ≧１）はそれぞれ、ｉ≧２の場合の（８）、（９）、（１０）、（１１）式において、関数ｆまたはｇによりｋパラメータに変換される（  ）内の数値を表し、標本状態量に相当する。ただし、ｍ＝ｉ−１とする。ｉ＝１の場合はｘn のバーが標本状態量に相当する。
【０１２６】
図４のステップＳ１においてまず、初期設定入力部２１または近似値更新装置２７から与えられたｘn のバーを関数ｆと関数ｇを用いて変換し、パラメータｋ1 、ｋ1 のバー、ｋ1 の波、ｋ1 の山を生成する。ただしｋ1 のバー、ｋ1 の波、ｋ1 の山は同じ値なので、これらの値の生成処理は一回でよい。
【０１２７】
次に図５のステップＳ２において、ステップＳ１で生成されたｋ1 、ｋ1 のバー、ｋ1 の波、ｋ1 の山の値から、パラメータｋ2 、ｋ2 のバー、ｋ2 の波、ｋ2 の山を生成する。
【０１２８】
パラメータｋ2 を生成するには、まずｋ1 、ｋ1 のバー、ｋ1 の波の値にメモリ２３─１、２３─２、２３─３に格納されているａ11、ａ11のバー、ａ11の波の値を乗算する。さらに乱数発生装置２６─１、２６─２が発生する正規乱数のサンプルを用いて（８）式の（  ）内の第２、第３、第４項を生成し、それらをｘn のバーに順次加算して、ｘn (1) を生成する。次に生成されたｘn (1) の値を関数ｆを用いて変換し、ｋ2 の値を得る。
【０１２９】
パラメータｋ2 のバーを生成するには、まずｋ1 、ｋ1 のバー、ｋ1 の山の値にメモリ２３─１、２３─２、２３─４に格納されているａ11、ａ11のバー、ａ11の山の値を乗算する。さらに乱数発生装置２６─１が発生する正規乱数のサンプルを用いて（９）式の（  ）内の第２、第３、第４項を生成し、それらをｘn のバーに順次加算して、ｘn (1) のバーを生成する。次に生成されたｘn (1) のバーの値を関数ｇを用いて変換し、ｋ2 のバーの値を得る。
【０１３０】
パラメータｋ2 の波を生成するには、まずｋ1 、ｋ1 の山の値にメモリ２３─１、２３─４に格納されているａ11、ａ11の山の値を乗算して（１０）式の（  ）内の第２、第３項を生成し、それらをｘn のバーに順次加算して、ｘn (1) の波を生成する。次に生成されたｘn (1) の波の値を関数ｇを用いて変換し、ｋ2 の波の値を得る。
【０１３１】
パラメータｋ2 の山を生成するには、まずｋ1 の山の値にメモリ２３─４に格納されているａ11の山の値を乗算して（１１）式の（  ）内の第２項を生成し、ｘn のバーに加算して、ｘn (1) の山を生成する。次に生成されたｘn (1) の山の値を関数ｇを用いて変換し、ｋ2 の山の値を得る。
【０１３２】
図６、図７のパラメータｋ3 、ｋ3 のバー、ｋ3 の波、ｋ3 の山およびｋ4 、ｋ4 のバー、ｋ4 の波、ｋ4 の山の生成処理についても、基本的に図５の処理と同様である。
【０１３３】
図６のステップＳ３では、ステップＳ１で生成されたｋ1 、ｋ1 のバー、ｋ1 の波、ｋ1 の山の値と、ステップＳ２で生成されたｋ2 、ｋ2 のバー、ｋ2 の波、ｋ2 の山の値から、ｘn (2) 、ｘn (2) のバー、ｘn (2) の波、ｘn (2) の山の値を求め、パラメータｋ3 、ｋ3 のバー、ｋ3 の波、ｋ3 の山を生成する。
【０１３４】
図７のステップＳ４では、ステップＳ１で生成されたｋ1 、ｋ1 のバー、ｋ1 の波、ｋ1 の山の値と、ステップＳ２で生成されたｋ2 、ｋ2 のバー、ｋ2 の波、ｋ2 の山の値と、ステップＳ３で生成されたｋ3 、ｋ3 のバー、ｋ3 の波、ｋ3 の山の値とから、ｘn (3) 、ｘn (3) のバー、ｘn (3) の波、ｘn (3) の山の値を求め、パラメータｋ4 、ｋ4 のバー、ｋ4 の波、ｋ4 の山を生成する。
【０１３５】
図８は時間増分生成装置２４の詳細な動作を示すフローチャートである。
時間増分生成装置２４は、まず時間増分Ψ（ｘn のバー，ｔn ，ｈ）の値をゼロに設定する。
【０１３６】
次にｋパラメータ生成装置２２から与えられるｋ1 、ｋ2 、ｋ3 、ｋ4 の値にメモリ２５─１に格納されたｂパラメータｂ1 、ｂ2 、ｂ3 、ｂ4 の値をそれぞれ乗算して、（７）式右辺の第１項の値を生成する。そして生成された第１項の値をΨ（ｘn のバー，ｔn ，ｈ）に加算する。
【０１３７】
次にｋパラメータ生成装置２２から与えられるｋ1 のバー、ｋ2 のバー、ｋ3 のバー、ｋ4 のバーの値にメモリ２５─２に格納されたｂパラメータｂ1 のバー、ｂ2 のバー、ｂ3 のバー、ｂ4 のバーの値をそれぞれ乗算する。さらに乱数発生装置２６─１が発生する正規乱数のサンプルを用いて（７）式右辺の第２項の値を生成する。そして生成された第２項の値をΨ（ｘn のバー，ｔn ，ｈ）に加算する。
【０１３８】
次にｋパラメータ生成装置２２から与えられるｋ1 の波、ｋ2 の波、ｋ3 の波、ｋ4 の波の値にメモリ２５─３に格納されたｂパラメータｂ1 の波、ｂ2 の波、ｂ3 の波、ｂ4 の波の値をそれぞれ乗算する。さらに乱数発生装置２６─２が発生する正規乱数のサンプルを用いて（７）式右辺の第３項の値を生成する。そして生成された第３項の値をΨ（ｘn のバー，ｔn ，ｈ）に加算する。
【０１３９】
次にｋパラメータ生成装置２２から与えられるｋ1 の山、ｋ2 の山、ｋ3 の山、ｋ4 の山の値にメモリ２５─４に格納されたｂパラメータｂ1 の山、ｂ2 の山、ｂ3 の山、ｂ4 の山の値をそれぞれ乗算して、（７）式右辺の第４項の値を生成する。そして生成された第３項の値をΨ（ｘn のバー，ｔn ，ｈ）に加算する。
【０１４０】
こうして得られたΨ（ｘn のバー，ｔn ，ｈ）の値が時刻ｔn における近似値ｘn のバーに対する時間増分に相当する。
近似値更新装置２７は時間増分生成装置２４から与えられる時間増分Ψ（ｘn のバー，ｔn ，ｈ）の値をｘn のバーの値に加算し、次の参照時刻ｔn+1 における状態量ｘ（ｔn+1 ）の近似値ｘn+1 のバーを生成する。そして生成されたｘn+1 のバーの値をｋパラメータ生成装置２２と出力部２８に与える。
【０１４１】
出力部２８は近似値更新装置２７から与えられる各参照時刻の近似値を格納しておき、近似値の更新が終了すると、これらの近似値をシミュレーション結果として、外部のディスプレイ装置や、シミュレーション対象の物理系を制御する制御装置等に出力する。
【０１４２】
次に本実施例のシミュレーション装置を具体的に、不規則外乱を受ける物理系に適用した結果を説明する。
図９は地球を周回する人工衛星の衛星軌道を示している。衛星軌道は、地球の大気密度の急速な変動と上層大気中の外乱等の影響を受けて変化し、これらが不規則外乱として作用する。
【０１４３】
図９において、状態量ｘは衛星軌道を含む平面内で、地球の中心から人工衛星３０に向かう放射方向における、与えられた衛星軌道からの摂動を表す。ここでｘ1 ＝ｘ、ｘ2 ＝ｄｘ／ｄｔと置くと、図９の衛星軌道は次のような２次元のモデルにより記述される。
【０１４４】
【数１８】
0003735128.tif 000039
【０１４５】
ここで、Ｘは（５）式の状態量に相当する２次元の状態量を表し、関数ｆとｇは（５）式の右辺の関数に相当し、２次元の関数値をもつ。またａ、ｂ、ｃは定数である。
【０１４６】
図１０は（１９）式で記述される衛星軌道のモデルを対象にして、図２のシミュレーション装置により数値シミュレーションを行った結果を示す。図１０のシミュレーションでは、ａ＝０．７５、ｂ＝１．５、ｃ＝０．５とした。またシミュレーションの時間ステップの刻み幅ｈはｈ＝２-8（秒）とし、時刻ｔ0 ＝０における初期条件として、ｘ1 （０）＝１．５（ｍ）、ｘ2 （０）＝−０．９（ｍ／ｓ）を与えた。またメモリ２３─１〜２３─４、２５─１〜２５─４に格納するａパラメータ、ｂパラメータとしては、（１６）式の値を用いた。
【０１４７】
図１１は（１９）式で不規則外乱の影響を無視した理想的な衛星軌道のモデルに対するシミュレーション結果を示す。図１０を図１１と比べると、大気密度の変動等の不規則外乱が人工衛星の軌道に与える影響が無視できないことがわかる。
【０１４８】
次に本シミュレーション装置の精度を他のシミュレーション装置と比較するために、次式で与えられる１次元のモデルを考える。
【０１４９】
【数１９】
0003735128.tif 000040
【０１５０】
これは、（５）式においてｘを１次元の状態量とし、関数ｆ、ｇの関数形をｆ（ｘ）＝ａｘ、ｇ（ｘ）＝ｂｘと置いたモデルに相当する。ただし、ａ、ｂは定数であり、ｘ0 は時刻ｔ0 におけるｘの初期値である。
【０１５１】
（２０）式の確率微分方程式は関数ｆ、ｇの関数形が簡単なため、その厳密解を解析的に求めることができる。この厳密解は次式で与えられる。
【０１５２】
【数２０】
0003735128.tif 000041
【０１５３】
ここでＷ（ｔ）はブラウン運動モデルの確率過程で、雑音モデルｗ（ｔ）に対応している。またｅｘｐの中の−（１／２）ｂ2 ｔの項は、確率積分に現れる特徴的な項である。
【０１５４】
図１２は（２１）式でａ＝１、ｂ＝１、ｘ0 ＝１と置き、２-9秒の時間刻みでサンプリングした結果を示す。また図１３、図１４はそれぞれ０．５次の近似精度のオイラ法を用いたシミュレーション装置、本実施例のシミュレーション装置によるシミュレーション結果を示す。図１３、図１４のシミュレーションにおいてもａ＝１、ｂ＝１、ｘ0 ＝１と置き、シミュレーションの時間刻みはｈ＝２-6（秒）とした。
【０１５５】
図１２、図１３、図１４を比較すると、０．５次の近似精度のオイラ法による結果よりも、１．５の近似精度を保証する本実施例のシミュレーション結果の方が厳密解のサンプリング結果に近いことがわかる。
【０１５６】
本実施例では、ａパラメータとｂパラメータの値の一例として（１６．１）、（１６．２）式の値を採用したが、本発明はこれに限定されるものではなく、（１４．１）〜（１４．２２）および（１５）式の拘束条件を満足する他の値も用いることができる。
【０１５７】
また（７）〜（１１）式の形の時間増分を用いれば、１．５次より高い近似精度を保証するシミュレーション装置の実現も可能である。
【０１５８】
【発明の効果】
本発明のシミュレーション装置によれば、不規則外乱の影響を受ける系に対して、１．５次の高精度を保証するシミュレーション結果を得ることができる。
【０１５９】
また数学モデルで与えられる関数の導関数を求める必要がなく、与えられた関数形のみを入力すればよいので、オペレータの負担が軽減される。さらに導関数の値を計算しなくて済むので、数値処理における計算量が少なく、汎用性に富むシミュレーション装置を提供する。
【０１６０】
このように、不規則外乱の影響を受ける系に対して、高い信頼性を持つシミュレーション結果を容易に得ることができるので、対象とする系の挙動の精密な解析と、その系を制御する高性能な制御装置の実現が可能になる。
【図面の簡単な説明】
【図１】本発明の原理図である。
【図２】本発明の一実施例のシミュレーション装置の構成図である。
【図３】一実施例のシミュレーション装置の動作を示すフローチャートである。
【図４】一実施例のｋパラメータ生成装置の動作を示すフローチャートである。
【図５】一実施例のｋパラメータ生成装置の動作を示すフローチャートである。
【図６】一実施例のｋパラメータ生成装置の動作を示すフローチャートである。
【図７】一実施例のｋパラメータ生成装置の動作を示すフローチャートである。
【図８】一実施例の時間増分生成装置の動作を示すフローチャートである。
【図９】人工衛星の衛星軌道を示す図である。
【図１０】図９の衛星軌道のモデルに対する本実施例のシミュレーション装置によるシミュレーション結果を示す図である。
【図１１】不規則外乱を受けない衛星軌道のモデルに対するシミュレーション結果を示す図である。
【図１２】一具体例のモデルの厳密解を示す図である。
【図１３】図１２のモデルに対する０．５次のオイラ法によるシミュレーション結果を示す図である。
【図１４】図１２のモデルに対する本実施例のシミュレーション装置によるシミュレーション結果を示す図である。
【符号の説明】
１０  シミュレーション装置
１１  初期設定入力手段
１２  増分生成手段
１３  近似値更新手段
１４  出力手段
２１  初期設定入力部
２２  ｋパラメータ生成装置
２３─１、２、３、４  メモリ
２５─１、２、３、４  メモリ
２４  時間増分生成装置
２６─１、２  乱数発生装置
２７  近似値更新手段
２８  出力部
３０  人工衛星
【図１】
0003735128.tif 000086
【図２】
0003735128.tif 000087
【図３】
0003735128.tif 000088
【図４】
0003735128.tif 000089
【図５】
0003735128.tif 000090
【図６】
0003735128.tif 000091
【図７】
0003735128.tif 000092
【図８】
0003735128.tif 000093
【図９】
0003735128.tif 000094
【図１０】
0003735128.tif 000095
【図１１】
0003735128.tif 000096
【図１２】
0003735128.tif 000097
【図１３】
0003735128.tif 000098
【図１４】
0003735128.tif 000099