(57)【特許請求の範囲】
【請求項１】
シミュレーションを実行する方法であって、コンピュータが、
複数の第１のシミュレーションを並列に走らせるステップであって、前記複数の第１のシミュレーションが、異なるパラメータによって特徴づけられる動作を有し且つ、前記動作が所定の観測時点から開始する、前記走らせるステップと、
前記所定の観測時点後の所定のタイミングでの実際の事象の実測値と、前記所定のタイミングでの前記複数の第１のシミュレーションそれぞれの計算結果との間の距離に従う離散分布を生成するステップと、
前記離散分布を用いて、前記複数の第１のシミュレーションから複数の第２のシミュレーションを生成するステップであって、前記複数の第２のシミュレーションの生成が、前記離散分布を用いて、前記所定のタイミングでの前記実測値に近い状態値を持つ複数のシミュレーション・サンプルを生成し、前記複数の第１のシミュレーションの一つと前記複数のシミュレーション・サンプルそれぞれとの間の重み付けの離散分布を生成することを含む、前記第２のシミュレーションを生成するステップと
を実行することを含む、前記方法。
【請求項２】
前記複数の第２のシミュレーションを生成するステップが、
前記重み付けの離散分布を移行尤度とし、当該移行尤度を利用して前記複数のシミュレーション・サンプルから前記複数の第２のシミュレーションを選択するステップと
を含む、請求項１に記載の方法。
【請求項３】
前記複数の第２のシミュレーションを生成するステップが、
前記複数のシミュレーション・サンプルから、状態値と前記パラメータを受け取るステップ
を含む、請求項１又は２に記載の方法。
【請求項４】
前記コンピュータが、
前記生成した複数の第２のシミュレーションを並列に走らせるステップ
をさらに実行することを含む、請求項１〜３のいずれか一項に記載の方法。
【請求項５】
前記コンピュータが、
所定の観測時点から開始して、実際の事象の実測値を観測して観測値を得るステップ
をさらに実行することを含む、請求項１〜４のいずれか一項に記載の方法。
【請求項６】
前記複数の第１のシミュレーションを並列に走らせるステップが、
前記複数の第１のシミュレーションそれぞれに、スレッドを割り当てるステップ
を含む、請求項１〜５のいずれか一項に記載の方法。
【請求項７】
前記所定のタイミングが、
前記所定の観測時点から所定時間経過後、又は、
前記実測値が得られたタイミング
である、請求項１〜６のいずれか一項に記載の方法。
【請求項８】
前記第１のシミュレーションが、交通シミュレーション、金融工学のリスク・シミュレーション、商品価格若しくは株価の変動シミュレーション、又は天気予報である、請求項１〜７のいずれか一項に記載の方法。
【請求項９】
前記実測値が、プローブカーデータ又は定点観測データである、請求項１〜８のいずれか一項に記載の方法。
【請求項１０】
シミュレーションを実行するシステムであって、
複数の第１のシミュレーションを並列に走らせる手段であって、前記複数の第１のシミュレーションが、異なるパラメータによって特徴づけられる動作を有し且つ、前記動作が所定の観測時点から開始する、前記走らせる手段と、
前記所定の観測時点後の所定のタイミングでの実際の事象の実測値と、前記所定のタイミングでの前記複数の第１のシミュレーションそれぞれの計算結果との間の距離に従う離散分布を生成する手段と、
前記離散分布を用いて、前記複数の第１のシミュレーションから複数の第２のシミュレーションを生成する手段であって、前記複数の第２のシミュレーションの生成が、前記離散分布を用いて、前記所定のタイミングでの前記実測値に近い状態値を持つ複数のシミュレーション・サンプルを生成し、前記複数のシミュレーション・サンプルの一つと前記複数のシミュレーション・サンプルそれぞれとの間の重み付けの離散分布を生成することを含む、前記複数の第２のシミュレーションを生成する手段と
を備えている、前記システム。
【請求項１１】
前記複数の第２のシミュレーションを生成する手段が、
前記重み付けの離散分布を移行尤度とし、当該移行尤度を利用して前記複数のシミュレーション・サンプルから前記複数の第２のシミュレーションを選択する、
請求項１０に記載のシステム。
【請求項１２】
前記複数の第２のシミュレーションを生成する手段が、
前記複数のシミュレーション・サンプルから、状態値と前記パラメータを受け取る処理を行う、
請求項１０又は１１に記載のシステム。
【請求項１３】
前記生成した複数の第２のシミュレーションを並列に走らせる手段をさらに備えている、請求項１０〜１２のいずれか一項に記載のシステム。
【請求項１４】
所定の観測時点から開始して、実際の事象の実測値を観測して観測値を得る手段をさらに備えている、請求項１０〜１３のいずれか一項に記載のシステム。
【請求項１５】
前記複数の第１のシミュレーションを並列に走らせる手段が、
前記複数の第１のシミュレーションそれぞれに、スレッドを割り当てる処理を行う、
請求項１０〜１４のいずれか一項に記載のシステム。
【請求項１６】
前記所定のタイミングが、
前記所定の観測時点から所定時間経過後、又は、
前記実測値が得られたタイミング
である、請求項１０〜１５のいずれか一項に記載のシステム。
【請求項１７】
前記第１のシミュレーションが、交通シミュレーション、金融工学のリスク・シミュレーション、商品価格若しくは株価の変動シミュレーション、又は天気予報であり、
前記実測値が、プローブカーデータ又は定点観測データである、
請求項１０〜１６のいずれか一項に記載のシステム。
【請求項１８】
シミュレーションを実行するコンピュータ・プログラムであって、コンピュータに、請求項１〜９のいずれか一項に記載の方法の各ステップを実行させる、前記コンピュータ・プログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
この発明は、コンピュータ・ベースのシミュレーション技法に関し、より詳しくは、交通シミュレーションなどの短期予測に適合するシミュレーション技法に関する。
【背景技術】
【０００２】
近年、都市計画の一環として、いわゆる、高度道路交通システム（Intelligent Transportation System）が、特に１９９０年代後半以降、活発に研究されている。
【０００３】
最近特に要望されてきているのは、自動車流入調整などの道路交通管制業務支援、混雑予測などのカーナビへの情報提供などの、リアルタイムの道路情報処理技術である。このような道路情報処理のベースとなるデータとして、例えば、実際に走行している自動車に取り付けられたセンサから得られたデータである、プローブカーデータがある。プローブカーデータが得られるようになると、それを用いてどのようにしてリアルタイムに、交通量、道路の混み具合などの交通情報を適切に処理するかが課題になりつつある。
【０００４】
シミュレーション一般の背景技術として、実測データに基づきシミュレーションを行う従来技術として、下記のものがある。
【０００５】
特開２０００−１６３４５１号公報は、シミュレーション計算を行う演算装置と、シミュレーション計算範囲の所定位置の観測データを求める観測機器を有する環境予測システムを開示する。このシステムにおいて、演算装置は、シミュレーション計算範囲を二次元座標マトリクスに区切り、それぞれの座標位置に対して、シミュレーション計算のための初期値を設定され、且つ境界条件が入力され、前記設定された初期値及び、入力された境界条件に基づき、所定期間経過後の各座標位置における観察項目の値をシミュレーション計算により求める。観測機器は、所定期間経過後の所定の座標位置に実際の観測データを求め、前記演算装置は、更に所定の座標位置のシミュレーション計算により求めた値と、観測機器による実際の観測データを比較し、差分があれば、前記初期値及び又は、境界条件を補正して再度シミュレーション計算を行うことを開示する。
【０００６】
特開２００２−２８６８６０号公報は、地下水の流れの解析を行う地下水流れ解析処理と、この地下水流れ解析処理にて得られる流れの条件下で移流分散遅延モデルにより流れ方向分散係数と遅延係数を物質輸送パラメータとして地下水中の物質移動解析とを行うシミュレーション装置において、前記移流分散遅延モデルの要素である地下水中の対象物質の地下水流れ方向に対する分散を表す流れ方向分散係数ａLおよび地下水の流れに対する地下水中の対象物質の移動の遅れを表す遅延係数Ｒdのうち、前記遅延係数Ｒdは所望の固定値ｒ0を用いると共に、流れ方向分散係数ａLは複数種を用い、これら複数種それぞれにおける前記物質の濃度分布または特定地点の濃度を経過時間対応に算出する物質輸送解析手段と、所要の経過時間ｍをおいて解析対象地での実測により得られた物質の濃度分布を実測値情報として用い、この実測値情報と前記物質輸送解析手段の算出結果とを比較して、前記算出した結果のうち、前記実測値情報に対する類似度が高い算出結果を選定し、当該算出結果の取得にあたって用いた流れ方向分散係数を以て流れ方向分散係数として決定し、かつ、この算出結果の該当する経過時間ｎと、前記実測での経過時間ｍから遅延係数をr0×ｍ／ｎとして決定する遅延係数決定手段と、を具備することを開示する。
【０００７】
また、交通状況監視に関する従来技術として、下記のものがある。
【０００８】
特開平９−１２８６７７号公報は、走行軌跡を計測する車載機を搭載した車両を走行させ、それによって得た走行軌跡データを用いて交通状況を監視し、また、信号タイミングデータの実績値および、車両走行モデルを用いて走行軌跡を実測の軌跡に近似させ、そのモデルを用いてシミュレーションを行い、信号タイミングを最適化することを開示する。
【０００９】
特開２００４−１１８７３５号公報は、現在から比較的長い将来にわたった交通状況を高い精度で予測するために、交通状況計算部で、道路ネットワークデータ、道路ネットワーク上の信号データ、及びＯＤ交通量を用いた交通流シミュレーションを行なって、交通量の計算値を計算し、将来予測処理部で、過去の交通量の観測値と現在の交通量の観測値とに基づいて、将来の交通量の予測値を予測し、ＯＤ交通量推定部で計算値が現在の交通量の観測値及び予測値に一致するように交通状況計算部で使用するＯＤ交通量を修正することを開示する。
【先行技術文献】
【特許文献】
【００１０】
【特許文献１】特開２０００−１６３４５１号公報
【特許文献２】特開２００２−２８６８６０号公報
【特許文献３】特開平９−１２８６７７号公報
【特許文献４】特開２００４−１１８７３５号公報
【発明の概要】
【発明が解決しようとする課題】
【００１１】
上記従来技術は、観測結果の実測値に基づき、実測値を近似するようにシミュレータを補正しつつ、シミュレーション動作を行うことは開示するが、シミュレータの精度を高めることについては、特に示唆すものではない。
【００１２】
従って、本発明の目的は、実測値に対する近似精度を高めることを可能ならしめるシミュレーション技法を提供することにある。
【００１３】
本発明の更なる目的は、コンピュータ・システムのマルチスレッド環境を利用して、高速且つ、実測値に対して高精度でシミュレーションを実行する技法を提供することにある。
【課題を解決するための手段】
【００１４】
本発明によれば、上記目的は、コンピュータ・システムにより、モンテカルロ・シミュレーションの原理を利用して、実測値に基づき、乱数的に複数のシミュレーションのスレッドを生成し、後の時点で別の実測値が得られたとき、その実測値から近いシミュレーションを複数残して、そこから更に複数のシミュレーションを再スタートさせるようにすることによって達成される。
【００１５】
本発明によれば、最初に、実測値に基づき、少しずつ動作の異なる複数のシミュレーションが並列に起動される。
【００１６】
シミュレーションの計算は、一般的に実世界の動きよりも高速なので、時刻T0で実測値が得られたとき、シミュレーションの計算の計算結果は、時刻T > T0まで進んでいる。ここで時刻Tは、コンピュータ上のシミュレーション時刻である。そこで、本発明によれば、既に時刻Tまで進んでいるシミュレーションの遡った値である、時刻T0での計算結果が、時刻T0で実測値と比較される。
【００１７】
これを可能とするため、生成された各シミュレーションにおいて、各時刻毎の計算結果が、後に参照可能にログとして保存される。
【００１８】
次に、再スタートされる前の複数のシミュレーションから、好適には重複を許して、実測値から近いほど大きい値をとるような第１の重み付けの離散分布が得られる。
【００１９】
次に、第１の重み付けの離散分布に従い複数のシミュレーションが生成され、そのシミュレーションの計算値相互間の近さを勘案して、近いほど大きい値をとるような第２の重み付けの離散分布、すなわち、シミュレーション間の移行尤度が得られる。
【００２０】
そして、この移行尤度に従い、再スタートされる前の複数のシミュレーションから個別に、次のシミュレーションが生成されて実行される。
【００２１】
本発明の１つの側面によれば、生成される各シミュレーションは、好適には、ある１つのシミュレーション動作関数で、それぞれ異なる引数またはパラメータを付与ことでバリエーションを与えられる。こうしておくと、シミュレーションを再スタートさせるとき、パラメータをコピーすることで、シミュレーションの動作を、再スタートするシミュレーションに継承することができる。
【００２２】
本発明の別の側面によれば、再スタートされるシミュレーションの数は、再スタートされる前のシミュレーションの数と同数で、発生された乱数と、上記第１及び第２の重み付け係数の組み合せにより、重複を許して、可能なシミュレーションから選ばれる。この際、重み付け係数の存在により、実測値に近く且つ移行尤度の高い傾向を示すシミュレーションがより選ばれやすくなる。
【００２３】
本発明の更に別の側面によれば、コンピュータ・システムは、マルチスレッド環境をもち、各々のシミュレーションは、個別のスレッドとして生成される。その際、好適には、コンピュータ・システムは、マルチコアまたはマルチプロセッサのシステムとして実装され、個別のコアまたはプロセッサに、シミュレーションのスレッドを割り当てることにより、シミュレーションがより高速に実行される。
【発明の効果】
【００２４】
この発明によれば、乱数的に複数生成されたシミュレーションのスレッドから、実測値に近いものがより高い尤度になるように重みをつけてシミュレーションを選び実行されることで、より高い精度のシミュレーションを実行することが可能となる。
【００２５】
また、マルチコアまたはマルチプロセッサのコンピュータ・システム上で、個別のコアまたはプロセッサに、シミュレーションのスレッドを割り当てることにより、シミュレーションがより高速に実行される。このことは、非常に複雑なシミュレーション動作を、実時間の計測よりも高速に実行させることを可能ならしめる。
【図面の簡単な説明】
【００２６】
【図１】シミュレーションを実施するための一例のハードウェアのブロック図である。
【図２】シミュレーションを実施するための一例の機能ブロック図である。
【図３】シミュレーション全体の処理のフローチャートを示す図である。
【図４】初期化の処理のフローチャートを示す図である。
【図５】シミュレーションの再配置処理のフローチャートを示す図である。
【図６】離散分布に従い、シミュレーションの複数のサンプルを生成する処理のフローチャートを示す図である。
【図７】シミュレーション処理を模式的に説明する図である。
【図８】シミュレーション処理を模式的に説明する図である。
【発明を実施するための形態】
【００２７】
以下、図面を参照して、本発明の一実施例の構成及び処理を説明する。以下の記述では、特に断わらない限り、図面に亘って、同一の要素は同一の符号で参照されるものとする。なお、ここで説明する構成と処理は、一実施例として説明するものであり、本発明の技術的範囲をこの実施例に限定して解釈する意図はないことを理解されたい。
【００２８】
先ず、図１を参照して、本発明を実施するために使用されるコンピュータのハードウェアについて説明する。図１において、ホスト・バス１０２には、複数のＣＰＵ１  １０４ａ、ＣＰＵ２  １０４ｂ、ＣＰＵ３  １０４ｃ、・・・ＣＰＵｎ  １０４ｎが接続されている。ホスト・バス１０２にはさらに、ＣＰＵ１  １０４ａ、ＣＰＵ２  １０４ｂ、ＣＰＵ３  １０４ｃ、・・・ＣＰＵｎ  １０４ｎの演算処理のためのメイン・メモリ１０６が接続されている。
【００２９】
一方、Ｉ／Ｏバス１０８には、キーボード１１０、マウス１１２、ディスプレイ１１４及びハードティスク・ドライブ１１６が接続されている。Ｉ／Ｏバス１０８は、Ｉ／Ｏブリッジ１１８を介して、ホスト・バス１０２に接続されている。キーボード１１０及びマウス１１２は、オペレータが、コマンドを打ち込んだり、メニューをクリックするなどして、操作するために使用される。ディスプレイ１１４は、必要に応じて、後述する本発明に係るプログラムをＧＵＩで操作するためのメニューを表示するために使用される。
【００３０】
この目的のために使用される好適なコンピュータ・システムのハードウェアとして、ＩＢＭ（Ｒ）Ｓｙｓｔｅｍ  Ｘがある。その際、ＣＰＵ１  １０４ａ、ＣＰＵ２  １０４ｂ、ＣＰＵ３  １０４ｃ、・・・ＣＰＵｎ  １０４ｎは、例えば、インテル（Ｒ）Ｘｅｏｎ（Ｒ）であり、オペレーティング・システムは、Ｗｉｎｄｏｗｓ（商標）Ｓｅｒｖｅｒ  ２００３である。オペレーティング・システムは、ハードティスク・ドライブ１１６に格納され、コンピュータ・システムの起動時に、ハードティスク・ドライブ１１６からメイン・メモリ１０６に読み込まれる。
【００３１】
本発明を実施するためには、必須ではないが、マルチプロセッサ・システムを用いることが望ましい。十分な数のコアまたはプロセッサがあれば、１つのコアまたはプロセッサに１つのシミュレーションのスレッドを割り当てて、高速に実行できるからである。ここでマルチプロセッサ・システムとは、一般に、独立に演算処理し得るプロセッサ機能のコアを複数もつプロセッサを用いるシステムを意図しており、従って、マルチコア・シングルプロセッサ・システム、シングルコア・マルチプロセッサ・システム、及びマルチコア・マルチプロセッサ・システムのどれかでよいことを理解されたい。
【００３２】
なお、本発明を実施するために使用可能なコンピュータ・システムのハードウェアは、ＩＢＭ（Ｒ）Ｓｙｓｔｅｍ  Ｘに限定されず、本発明のシミュレーション・プログラムを走らせることができるものであれば、通常のパーソナル・コンピュータなどの任意のコンピュータ・システムを使用することができる。オペレーティング・システムも、Ｗｉｎｄｏｗｓ（Ｒ）に限定されず、Ｌｉｎｕｘ（Ｒ）、Ｍａｃ  ＯＳ（Ｒ）など、任意のオペレーティング・システムを使用することができる。さらに、シミュレーション・プログラムを高速で動作させるために、ＰＯＷＥＲ（商標）６ベースで、オペレーティング・システムがＡＩＸ（商標）のＩＢＭ（Ｒ）Ｓｙｓｔｅｍ  Ｐなどのコンピュータ・システムを使用してもよい。
【００３３】
ハードティスク・ドライブ１１６にはさらに、図２で後述する本発明に係るシミュレータ２０４、重み計算モジュール２０６、シミュレータ選択モジュール２０８、シミュレーションを実行するためのパラメータ２１０などが格納されており、オペレータのキーボードやマウス操作に応答して、メイン・メモリ１０６にロードされて実行される。
【００３４】
Ｉ／Ｏバス１０８にはさらに、通信インターフェース１２０により、インターネット１２２を介して、観測装置１２４に接続されている。観測装置１２４は、この実施例では、プローブカー・データの観測装置であり、図１のシステムは、通信インターフェース１２０を介して、一定周期、あるいは不定周期にプローブカー・データからデータを取得することができる。尚ここで、プローブカー・データでなく、定点観測データを用いるようにしてもよい。
【００３５】
通信インターフェース１２０は、例えば、イーサネット・カードであり、この場合、図示しないが、特定のウェブ・サーバを介して、観測装置１２４に有線接続される。
【００３６】
あるいは通信インターフェース１２０は、データ通信カードを用いて、サーバにダイアルアップ接続するタイプのものでもよく、さらには公衆無線ＬＡＮによるタイプのものでもよい。本発明において、観測装置１２４に対する接続形態に制限はなく、任意の接続構成を用いることができることを理解されたい。
【００３７】
次に、図２の機能ブロックを参照して、この実施例で使用されるプログラム・モジュールとその機能について説明する。
【００３８】
図２において、制御モジュール２０２は、ハードティスク・ドライブ１１６に保存されており、キーボード１１０やマウス１１２などの操作に応答して、オペレーティング・システムの働きにより、メイン・メモリ１０６にロードされて、実行される。
【００３９】
制御モジュール２０２は、シミュレーション全体の動作を管理する役割を果たすものであり、これには限定されないが、複数のＣＰＵ１  １０４ａ、ＣＰＵ２  １０４ｂ、ＣＰＵ３  １０４ｃ、・・・ＣＰＵｎ  １０４ｎのうちの１つに割り当てられる管理スレッドとして実装されてもよい。
【００４０】
本発明によれば、一般的に、多くのシミュレーションを並列に走らせるほど、精度が高まる。並列に走らせることができるシミュレーションの最大数は、システムの可用なリソース、オペレーティング・システムの機能、システムのハードウェア・アーキテクチャ、必要なシミュレーションの動作速度などで決定される。
【００４１】
シミュレータ２０４は、この実施例では、道路上の自動車の動きをシミュレートするものであり、関数Fで以って、下記のように記述される。
xn+1 = xn + Δt
yn+1 = F(xn,yn,θi)
【００４２】
この式で、n → n+1 を繰り返して、xnが所定の時刻に到達するまで、シミュレーションのスレッドが走ることになる。
ここで、xnはシミュレーション時間であり、Δtはクロック1ステップ分の時間刻みである。
また、ynは、一般的に、複数次元をもつベクトル量であり、シミュレーションの状態量をあらわす。θiは、シミュレータ２０４のパラメータであり、このθiとして異なる値を与えることにより、シミュレータ２０４は、異なる挙動を示す。θiは一般的に、複数の要素からなるベクトル量である。パラメータθiの要素には、シミュレータ２０４で使用される乱数の種も含まれてもよい。
【００４３】
シミュレータ２０４は、動作中に、シミュレーション時刻T = xnと、
計算結果ynの特徴量S(T)を、例えば、ハードディスク・ドライブ１１６の所定領域に記録する。すると、制御モジュール２０２は、任意の時点で、所与の時刻のシミュレータ２０４の計算結果を取得できる。なお、実測値が得られる時間が予め分かっていれば、計算結果ynの特徴量S(T)は、その時間付近だけ記録するようにしてもよい。
【００４４】
制御モジュール２０２は、起動した任意のシミュレーションのスレッドを停止、あるいは消滅させることができる。このような機能は、Ｗｉｎ３２ＡＰＩ、Ｌｉｎｕｘ  ＡＰＩなど、通常のオペレーティング・システムのＡＰＩとして標準で備わっている機能である。
【００４５】
制御モジュール２０２は、メイン・メモリ１０６などのリソースが許す限り、任意の数のシミュレータ２０４を起動して、個別のスレッドとして並列に走らせることができる。
【００４６】
シミュレーションのためのパラメータ２１０がハードディスク・ドライブ１１６に保存されて、そこには、シミュレータ２０４に引数として与えるためのパラメータθi(i = 1,...,k)が含まれる。
【００４７】
制御モジュール２０２が参照しやすくするため、パラメータθi(i = 1,...,k)は、システムの起動時に、メイン・メモリ１０６の所定の領域にロードしておいてもよい。
【００４８】
制御モジュール２０２は、シミュレータ２０４を複数個、個別のスレッドとして起動する際、パラメータの引数として、異なるθi(i = 1,...,k)を順次与えて、少しずつ挙動の異なるシミュレーションを走らせることができる。
【００４９】
制御モジュール２０２がシミュレータ２０４を複数起動する際、利用可能なプロセッサまたはコアが、起動するシミュレータ２０４の数より多い場合、制御モジュール２０２は、個々のプロセッサまたはコアに、起動された個々のシミュレータ２０４のスレッドが１つずつアサインされる。
【００５０】
重み計算モジュール２０６は、制御モジュール２０２からコールされて、観測値と起動されたシミュレータ２０４の計算結果に基づき、後述するシミュレーションの再配置のために、重みの係数を計算する機能をもつ。
【００５１】
シミュレータ選択モジュール２０８は、制御モジュール２０２からコールされて、重み計算モジュール２０６から得られた重み係数と、起動されたシミュレータ２０４の計算結果に基づき、重複を許して、ｋ個のシミュレータ２０４を生成する機能をもつ。
【００５２】
重み計算モジュール２０６とシミュレータ選択モジュール２０８のより詳しい機能は、後で、図３以下のフローチャートを参照して説明する。
【００５３】
表示モジュール２１２は、制御モジュール２０２により制御されて、シミュレーションの途中経過を、ディスプレイ１１４に表示する機能を有する。
【００５４】
データ観測モジュール２１４は、制御モジュール２０２により制御されて、観測装置１２４から、観測データと、その観測時刻を取得し、メイン・メモリ１０６の所定領域に配置する機能をもつ。
【００５５】
尚、図２に示すモジュールは、使用されているオペレーティング・システム及びシステム環境に適合する限り、Ｃ、Ｃ＋＋、Ｃ＃、Ｊａｖａ(R)などの既知の任意のプログラミング言語により作成することができる。
【００５６】
次に、図３〜図６のフローチャートを参照して、この実施例のシミュレーション実行処理を説明する。図３のステップ３０２で、制御モジュール２０２は、シミュレーション開始に当たっての初期化を行う。
【００５７】
図４のフローチャートは、初期化の処理をより詳細に示す。図４において、制御モジュール２０２は、ステップ４０２で、データ観測モジュール２１４の機能を利用して、実状態s(0)を観測する。これは、時間T = 0での観測値に対応する。s(0)は、シミュレータ２０４の関数Fの計算結果に対応する次元と状態量をもつ。s(0)は、観測値自体であってもよいし、観測値に基づき計算された特徴量であってもよい。この実施例では、s(0)は、観測値に基づき計算された特徴量であるとする。
【００５８】
ステップ４０４では、制御モジュール２０２は、パラメータθi(i = 1,...,k)を、メイン・メモリ１０６または、ハードディスク・ドライブ１１６の所定の領域から読み込む。
【００５９】
図３に戻って、制御モジュール２０２は、観測した初期値と、パラメータθi(i = 1,...,k)により、下記のようにk個のシミュレーションを並列に走らせる。
ここで、y(i)0 = s(0) (i = 1,...,k)とする。
x(1)n+1 = x(1)n + Δt
y(1)n+1 = F(x(1)n,y(1)n,θ1)
---------------------------------------------------------
x(2)n+1 = x(2)n + Δt
y(2)n+1 = F(x(2)n,y(2)n,θ2)
---------------------------------------------------------
x(3)n+1 = x(3)n + Δt
y(3)n+1 = F(x(3)n,y(3)n,θ3)
---------------------------------------------------------
・・・・・・・・・・・・・・・・・・・・・・・・・・
・・・・・・・・・・・・・・・・・・・・・・・・・・
---------------------------------------------------------
x(k)n+1 = x(k)n + Δt
y(k)n+1 = F(x(k)n,y(k)n,θk)
【００６０】
ステップ３０６では、制御モジュール２０２は、所定時間（例えば、T0)が経過したかどうか、判断する。もしそうでなければ、ステップ３０８で、k個のシミュレーションは、そのまま並列実行される。
【００６１】
ステップ３０６で、所定時間が経過したと判断されると、ステップ３１０で、シミュレーションの再配置が行われ、ステップ３１２では、再配置された結果のk個のシミュレーションが並列実行される。尚、ステップ３１０の再配置処理は、図５のフローチャートを参照して、後でより詳細に説明する。
【００６２】
ステップ３０６の判断は、例えば、カウンタとしての所定の変数に値をセットして、デクリメントし、その値がゼロに達したことに応答して行う。シミュレーションの再配置が行われた後は、カウンタとしての所定の変数に値をセットして、デクリメントする処理が再開される。こうして、シミュレーションの再配置は、所定時間毎に行われることになる。
【００６３】
ステップ３０６の判断の段階で、所定時間が経過しても、何かの事情で、データ観測モジュールで実測値が得られていない、という場合がありえる。その場合は、所定時間が経過してもステップ３０６の判断を直ちに真にせず、後で実測値が更新されていたことに応答して、ステップ３０６の判断を真にしてもよい。
【００６４】
ステップ３１４では、制御モジュール２０２は、終了条件に達したかどうかを判断する。ここで終了条件とは例えば、シミュレーション開始後、所定の時間が経過したかどうか、あるいは、観測結果が所定の値になったかどうか、などである。
【００６５】
さらに、観測結果がある程度定期的に得られるなら、ステップ３０６の判断を所定時間の経過とせず、観測結果が入来したタイミングとすることもできる。
【００６６】
次に、図５のフローチャートを参照して、ステップ３１０の再配置処理を、より詳細に説明する。
【００６７】
ステップ５０２では、制御モジュール２０２は、実状態sを観測した時間T0と、現在のシミュレーション時間Tを読み取る。実状態sを観測した時間T0は、データ観測モジュール２１４が、観測値とその観測時間を、メモリ上の所定の記憶領域に配置するので、読み取ることができる。また、シミュレーション時間Tは、システムのクロックから読み取ることができる。このとき、シミュレーションは、実世界よりも高速なので、T > T0である。
【００６８】
ステップ５０４では、制御モジュール２０２が、重み計算モジュール２０６を呼び出して、重み値w(1)(i)を、i = 1,...,kに亘って、下記の式により、計算する。
【００６９】
すなわち、k個のシミュレーションの時間T0における値をそれぞれ、
S1(T0),S2(T0),...,Sk(T0)とする。ここで、Si(T0)というのは、上記の式で、y(i)nと同じものである。T0がステップ数nに対応する。
また、時間T0における実測値の特徴量をs(T0)とする。すなわち、シミュレーションの計算結果は、特徴量に対応している。
そこで、ある定義された距離関数d1を以って、
d1(Si(T0),s(T0))が大きいほど、w(1)(i)が小さくなるようにw(1)(i)が選ばれる。ここで、d1は、これには限定されないが、例えば、Si(T0),s(T0)をベクトルとしたとき、ユークリッド距離でよい。あるいは、用途に応じて、マンハッタン距離なども使用することができる。
また、これには限定されないが、例えば、 w(1)(i) 〜  exp(-d1(Si(T0),s(T0)))のように選んでよい。w(1)(i) は、実測値とシミュレーションの計算値の単調減少関数になっていればよい。
こうしておいて、下記のように、w(1)(i)をスケーリングする。
【数１】
0005633961.tif 000002
重み値w(1)(i)は、実測値と近い状態値をもつシミュレーションが選ばれやすくなるように分布される。
【００７０】
ステップ５０６では、制御モジュール２０２が、P(i) = w(1)(i) i = 1, ...,kの離散分布に従って、シミュレータ選択モジュール２０８を呼び出し、k個のシミュレーション・サンプル{j}を得る。尚、シミュレータ選択モジュール２０８の詳細な処理は、図６のフローチャートを参照して、後で説明する。
【００７１】
ステップ５０８では、制御モジュール２０２が、ある定義された距離関数d2を以って、ステップ５０６で生成されたk個のシミュレーション・サンプル{j}の各要素jにおいて、
d2((Si(T),θi),(Sj(T),θj))が大きいほど、w(2)j(i)が小さくなるようにw(2)j(i)を選ぶ。ここで、d2は、これには限定されないが、例えば、(Si(T),θi),(Sj(T),θj)をベクトルとしたとき、ユークリッド距離でよい。d2の式は、d1と同様でもよく、異なっていてもよい。さらに、(Si(T),θi)という記法は、ベクトルSi(T)と、ベクトルθiを並べたベクトルを意味する。前述のように、Si(T)は特徴量なので、パラメータθiも一種の特徴量と解釈してもよい。
また、これには限定されないが、例えば、 w(2)j(i) 〜  exp(-d2((Si(T),θi),(Sj(T),θj))))のように選んでよい。w(2)j(i)も、実測値とシミュレーションの計算値の単調減少関数になっていればよい。
こうしておいて、下記のように、w(2)j(i)をスケーリングする。
【数２】
0005633961.tif 000003
重み値w(2)j(i)は、シミュレーションjからシミュレーションiへの移行尤度を示す。
【００７２】
ステップ５１０では、制御モジュール２０２が、各jについて、Pj(i) = w(2)j(i) i = 1, ...,kの離散分布に従って、シミュレータ選択モジュール２０８を呼び出し、k個のシミュレーション・サンプル{i*}を得る。すなわち、このステップでは、各jについて、
Pj(i) = w(2)j(i) i = 1, ...,kの離散分布に従って、その都度シミュレータ選択モジュール２０８を呼び出し、１つずつシミュレーション・サンプルが選ばれ、j = 1,..,kまで進んだとき、k個のシミュレーション・サンプル{i*}が選ばれている、という次第である。
【００７３】
ステップ５１２では、制御モジュール２０２が、ステップ５１０で生成されたk個のシミュレーション・サンプル{i*}を新たなシミュレータ集合とするように、{i*}のシミュレーション時刻Tでのシミュレーション状態S*i(T)と、パラメータθ*iを、元のシミュレータ集合に下記のように順次コピーする。
for (i = 1;i <= k; i++ )  {
Si(T) = S*i(T);
θi = θ*i;
}
ここでパラメータθ*iとあるのは、実際は、θi ( i = 1,...,k)から選ばれたものであることに留意されたい。
【００７４】
次に図６のフローチャートを参照して、シミュレータ選択モジュール２０８の処理を説明する。この処理は、離散分布P(i) = w(i) i = 1,...,kと、選択するシミュレータの数Nを引数とする。
【００７５】
ステップ６０２では、シミュレータ選択モジュール２０８は、ローカル変数m = 1とおく。
【００７６】
ステップ６０４では、シミュレータ選択モジュール２０８は、既知のアルゴリズムを用いて、範囲(0,1)から一様乱数を生成し、それをxとおく。
【００７７】
ステップ６０４では、シミュレータ選択モジュール２０８は、下記の式をみたすjを選択する。
【数３】
0005633961.tif 000004
但し、
【数４】
0005633961.tif 000005
とする。
【００７８】
これによってjが決まるので、ステップ６０８で、m番目のシミュレータとしてシミュレータjが選択される。これは言い換えると、m番目のシミュレータは、パラメータとして、θjを使用することを意味する。
【００７９】
ステップ６１０では、mが1増分され、ステップ６１２でm > Nかどうかが判断される。そうでなければ、処理はステップ６０４に戻り、そうであるなら、所定の数のシミュレータが選ばれたので、処理を終って、シミュレータ選択モジュール２０８は、もとの処理に戻る。
【００８０】
尚、xの乱数性により、ステップ６０８で、同一のシミュレータが重複して選ばれることがあることを理解されたい。
【００８１】
図７は、この実施例のシミュレータの動作を模式的に示す図である。ここで、縦軸は、シミュレーションの状態値、横軸は時間である。また、説明の便宜上、並列に走るシミュレータは３個だけだとする。
【００８２】
この図は、実測値がt = T0で得られたとき、シミュレーション時刻は、T > T0まで進んでいることを示す。ここで、図３のステップ３０６により、シミュレーションの再配置のステップ３１０が呼び出される。
【００８３】
すると、本発明によれば、t = T0での実測値とシミュレーションの計算値の距離に従い、ステップ５０４で、重みw(1)(i)が計算される。
【００８４】
ステップ５０６で、シミュレータ選択モジュール２０８が呼び出され、重みw(1)(i)に従い、k個のシミュレーション・サンプルが選ばれる。
【００８５】
そのようにして選ばれたシミュレーション・サンプルに基づき、ステップ５０８において、重みw(2)j(i)が計算され、それに基づき、各j毎にシミュレータ選択モジュール２０８が呼ばれて、改めて結果的にk個のシミュレーション・サンプルが選ばれる。
【００８６】
そのk個のシミュレーション・サンプルの状態値とパラメータが、ステップ５１２で、元のシミュレータ集合の引数としてコピーされ、ステップ３１２で、その結果再配置されたシミュレーションが実行される。
【００８７】
図８は、その結果選ばれたシミュレータS'1, S'2及びS'3がシミュレーション時刻Tから始まって、先に進んでいる様子を示す。
【００８８】
このようなシミュレーションの挙動は、表示モジュール２１２の機能により、ディスプレイ１１４に表示される。
【００８９】
このようにシミュレーション・スレッドが複数あると、それらの状態値を複数のシミュレーションで平均した結果を、最終的にシミュレーションされた結果とみなすことができる。
【００９０】
上記実施例では、第１の重み係数としてのw(1)(i)と、第２の重み係数としてのw(2)j(i)の両方を用いて、シミュレーションの離散分布を生成するようにした。この実施例は、精度の点で好ましいが、第１の重み係数のみを用いてシミュレーションの離散分布を生成しても、シミュレーションを実現することができることを理解されたい。このとき、生成したj番目のシミュレーションがもとのシミュレーションのi番目に対応するなら、生成したj番目のシミュレーションに、もとのi番目のシミュレーションの計算結果とパラメータがコピーされる。
【００９１】
以上、本発明を、プローブカーデータまたは定点観測データの測定結果に基づく、交通シミュレーションを例にとって説明したが、本発明は、これには限定されず、金融工学のリスク・シミュレーション、商品価格や株価の変動シミュレーション、天気予報など、実測値に基づき離散分布に従う複数のシミュレーションを確率的に発生されることが有効であるさまざまな分野に適用可能であることを理解されたい。
【００９２】
また、例示されているハードウェア構成は一例であって、本発明は、マルチタスクあるいははマルチスレッドをサポートする環境であるなら、特定のハードウェア、オペレーティングシステムに限定されるものではない。
【符号の説明】
【００９３】
１０４ａ・・・１０４ｎ  ＣＰＵ
１０６  メイン・メモリ
１０８  バス
１１０  キーボード
１１２  マウス
１１４  ディスプレイ
１１６  ハードティスク・ドライブ
１２４  観測装置
２０２  制御モジュール
２０４  シミュレータ
２０６  計算モジュール
２０８  シミュレータ選択モジュール
２１０  パラメータ
２１２  表示モジュール
２１４  データ観測モジュール
【図１】
0005633961.tif 000006
【図２】
0005633961.tif 000007
【図３】
0005633961.tif 000008
【図４】
0005633961.tif 000009
【図５】
0005633961.tif 000010
【図６】
0005633961.tif 000011
【図７】
0005633961.tif 000012
【図８】
0005633961.tif 000013