特開 2019-160180 シミュレーションモデルの時系列シミュレーション値の推定のためのシステム、方法、プログラム、プログラムを記憶した記憶媒体

(19)【発行国】日本国特許庁(JP)

(12)【公報種別】公開特許公報(A)

(11)【公開番号】特開2019-160180(P2019-160180A)

(43)【公開日】2019年9月19日

(54)【発明の名称】シミュレーションモデルの時系列シミュレーション値の推定のためのシステム、方法、プログラム、プログラムを記憶した記憶媒体

(51)【国際特許分類】

   G16Z 99/00 20190101AFI20190823BHJP

   B64G 7/00 20060101ALI20190823BHJP

【ＦＩ】

   G06F19/00 110

   B64G7/00 B

【審査請求】未請求

【請求項の数】6

【出願形態】ＯＬ

【全頁数】19

(21)【出願番号】特願2018-49582(P2018-49582)

(22)【出願日】2018年3月16日

(71)【出願人】

【識別番号】503361400

【氏名又は名称】国立研究開発法人宇宙航空研究開発機構

(74)【代理人】

【識別番号】100094569

【弁理士】

【氏名又は名称】田中 伸一郎

(74)【代理人】

【識別番号】100088694

【弁理士】

【氏名又は名称】弟子丸 健

(74)【代理人】

【識別番号】100103610

【弁理士】

【氏名又は名称】▲吉▼田 和彦

(74)【代理人】

【識別番号】100067013

【弁理士】

【氏名又は名称】大塚 文昭

(74)【代理人】

【識別番号】100086771

【弁理士】

【氏名又は名称】西島 孝喜

(74)【代理人】

【識別番号】100109070

【弁理士】

【氏名又は名称】須田 洋之

(74)【代理人】

【識別番号】100109335

【弁理士】

【氏名又は名称】上杉 浩

(74)【代理人】

【識別番号】100120525

【弁理士】

【氏名又は名称】近藤 直樹

(74)【代理人】

【識別番号】100151987

【弁理士】

【氏名又は名称】谷口 信行

(72)【発明者】

【氏名】加藤 博司

(72)【発明者】

【氏名】安藤 麻紀子

(72)【発明者】

【氏名】福添 森康

【テーマコード（参考）】

5L049

【Ｆターム（参考）】

5L049DD02

(57)【要約】      （修正有）

【課題】予測精度が高く、不確定因子の網羅的検討を実施することができる、時系列シミュレーション値の推定のためのシステム及び方法を提供する。
【解決手段】時系列シミュレーション値を推定するシステム１０は、初期設定部と、サンプル点生成部と、複数の初期サンプル点の各々に対応するシミュレーション値を計算するシミュレーション部と、シミュレーション値から、所定の複数の時刻の各々に対する各時刻エミュレーションモデルを生成する各時刻エミュレーションモデル生成部と、複数のサンプル点の値を各時刻エミュレーションモデルに入力し、所定の複数の時刻の各々に対する、複数のサンプル点の各々に対応するエミュレーション値を計算し、所定の複数の時刻の各々に対する、複数のサンプル点の各々に対応するエミュレーション値に基づいて、複数のサンプル点の各々に対する時系列シミュレーション値を推定するエミュレーション部とを備える。
【選択図】図１

【特許請求の範囲】

【請求項1】

不確実性を有する複数の説明変数を含むシミュレーションモデルの時系列シミュレーション値を推定するシステムであって、
前記複数の説明変数の各々の不確実性を設定する初期設定部と、
前記複数の説明変数の各々の前記不確実性にしたがって、前記複数の説明変数の、複数の初期サンプル点及び複数のサンプル点を生成するサンプル点生成部と、
前記複数の初期サンプル点の値を、前記シミュレーションモデルに入力し、所定の複数の時刻の各々に対する、前記複数の初期サンプル点の各々に対応するシミュレーション値を計算するシミュレーション部と、
前記複数の初期サンプル点の各々と、前記所定の複数の時刻の各々に対応する、前記複数の初期サンプル点の各々に対応するシミュレーション値から、前記所定の複数の時刻の各々に対する各時刻エミュレーションモデルを生成する各時刻エミュレーションモデル生成部と、
前記複数のサンプル点の値を前記各時刻エミュレーションモデルに入力し、前記所定の複数の時刻の各々に対する、前記複数のサンプル点の各々に対応するエミュレーション値を計算し、前記所定の複数の時刻の各々に対する、前記複数のサンプル点の各々に対応するエミュレーション値に基づいて、前記複数のサンプル点の各々に対する前記時系列シミュレーション値を推定するエミュレーション部と、
を備えるシステム。

【請求項2】

前記各時刻エミュレーションモデルはノンパラメトリック回帰モデルである請求項１に記載のシステム。

【請求項3】

前記ノンパラメトリック回帰モデルは、ガウス過程回帰モデルである請求項２に記載のシステム。

【請求項4】

不確実性を有する複数の説明変数を含むシミュレーションモデルの時系列シミュレーション値を推定する、コンピュータにより実行される方法であって、
前記複数の説明変数の各々の不確実性を設定するステップと、
前記複数の説明変数の各々の前記不確実性にしたがって、前記複数の説明変数の、複数の初期サンプル点を生成するステップと、
前記複数の初期サンプル点の値を、前記シミュレーションモデルに入力し、所定の複数の時刻の各々に対する、前記複数の初期サンプル点の各々に対応するシミュレーション値を計算するステップと、
前記複数の初期サンプル点の各々と、前記所定の複数の時刻の各々に対する、前記複数の初期サンプル点の各々に対応するシミュレーション値から、前記所定の複数の時刻の各々に対する各時刻エミュレーションモデルを生成するステップと、
前記複数の説明変数の各々の前記不確実性及び／又は条件にしたがって、前記複数の説明変数の、複数のサンプル点を生成し、
前記複数のサンプル点の値を前記各時刻エミュレーションモデルに入力し、前記所定の複数の時刻の各々に対する、前記複数のサンプル点の各々に対応するエミュレーション値を計算するステップと、
前記所定の複数の時刻の各々に対する、前記複数のサンプル点の各々に対応するエミュレーション値に基づいて、前記複数のサンプル点の各々に対する前記時系列シミュレーション値を推定するステップと、
を含む方法。

【請求項5】

請求項４に記載の方法をコンピュータに実行させるためのプログラム。

【請求項6】

請求項５に記載のプログラムを記憶したコンピュータ読み取り可能な記憶媒体。

【発明の詳細な説明】

【技術分野】

【0001】

本発明は、シミュレーションモデルの時系列シミュレーション値の推定のためのシステム、方法、プログラム、プログラムを記憶した記憶媒体に関する。より詳細には、不確実性を有する複数の説明変数を用いるシミュレーションモデルの時系列シミュレーション値の推定のためのシステム、方法、プログラム、プログラムを記憶した記憶媒体に関する。

【背景技術】

【0002】

システムの動的評価を行うためには、対象とするシステムの挙動を記述するモデル（以下、「システムモデル」という。）を時系列で解く必要がある。時系列データに対してシステムモデルの入力−出力間の関係性を近似モデル化して時系列評価を実施する方法としては、ＡＲ (Autoregressive)モデルや、ＲＮＮ（Recurrent Neural Network）モデル等がある（下記特許文献１、２参照）。

【先行技術文献】

【特許文献】

【0003】

【特許文献1】特開２０１７−１２４０１４号公報

【特許文献2】特開２０１６−２１２２７３号公報

【発明の概要】

【発明が解決しようとする課題】

【0004】

しかしながら、システムモデルの時系列計算は、一般的に、定常計算と比較して計算コストが高く、また、システムモデル自体の計算負荷が高い場合は、限られた計算回数での評価しか実施できない。そのため、外部環境条件や、システムモデル内に不確定性が存在する場合には、不確定性因子の網羅的検討が実施できずに、得られた計算結果に対して安全率やマージン等を考慮し、システムの安全性等を評価する必要がある。しかし、過剰な安全率やマージンは、システムの肥大化を招き、システムの性能向上や開発コスト低減を目指す上では、適切に不確定性因子の影響を考慮していく必要がある。

【0005】

ＡＲモデルやＲＮＮモデルは、時系列で統一的なモデルを作成することで、入力条件に対する時系列情報を予測しようとするアプローチであるが、時系列で統一的なモデルを作成しようとした時、大量の学習データが必要であり、学習データを計算から作成するような場合、計算コストが高くなる。そのため、少ないデータ数しか事前収集できない、もしくは、限られた計算しか実施できない場合は、ＡＲモデルやＲＮＮモデルでは十分な予測精度が出ない。

【0006】

そこで、本発明は、少ないデータ数しか事前収集できない、もしくは、限られた計算しか実施できない場合でも、予測精度が高く、不確定因子の網羅的検討を実施することができる、シミュレーションモデルの時系列シミュレーション値の推定のためのシステム、方法を提供することを目的の１つとする。

【課題を解決するための手段】

【0007】

本発明の１つの態様は、不確実性を有する複数の説明変数を含むシミュレーションモデルの時系列シミュレーション値を推定するシステムであって、前記複数の説明変数の各々の不確実性を設定する初期設定部と、前記複数の説明変数の各々の前記不確実性にしたがって、前記複数の説明変数の、複数の初期サンプル点及び複数のサンプル点を生成するサンプル点生成部と、前記複数の初期サンプル点の値を、前記シミュレーションモデルに入力し、所定の複数の時刻の各々に対する、前記複数の初期サンプル点の各々に対応するシミュレーション値を計算するシミュレーション部と、前記複数の初期サンプル点の各々と、前記所定の複数の時刻の各々に対応する、前記複数の初期サンプル点の各々に対応するシミュレーション値から、前記所定の複数の時刻の各々に対する各時刻エミュレーションモデルを生成する各時刻エミュレーションモデル生成部と、前記複数のサンプル点の値を前記各時刻エミュレーションモデルに入力し、前記所定の複数の時刻の各々に対する、前記複数のサンプル点の各々に対応するエミュレーション値を計算し、前記所定の複数の時刻の各々に対する、前記複数のサンプル点の各々に対応するエミュレーション値に基づいて、前記複数のサンプル点の各々に対する前記時系列シミュレーション値を推定するエミュレーション部とを備えるシステムを提供するものである。

【0008】

前記各時刻エミュレーションモデルはノンパラメトリック回帰モデルであるものとすることができる。

【0009】

前記ノンパラメトリック回帰モデルは、ガウス過程回帰モデルであるものとすることができる。

【0010】

本発明の１つの態様は、不確実性を有する複数の説明変数を含むシミュレーションモデルの時系列シミュレーション値を推定する、コンピュータにより実行される方法であって、前記複数の説明変数の各々の不確実性を設定するステップと、前記複数の説明変数の各々の前記不確実性にしたがって、前記複数の説明変数の、複数の初期サンプル点を生成するステップと、前記複数の初期サンプル点の値を、前記シミュレーションモデルに入力し、所定の複数の時刻の各々に対する、前記複数の初期サンプル点の各々に対応するシミュレーション値を計算するステップと、前記複数の初期サンプル点の各々と、前記所定の複数の時刻の各々に対する、前記複数の初期サンプル点の各々に対応するシミュレーション値から、前記所定の複数の時刻の各々に対する各時刻エミュレーションモデルを生成するステップと、前記複数の説明変数の各々の前記不確実性及び／又は条件にしたがって、前記複数の説明変数の、複数のサンプル点を生成し、前記複数のサンプル点の値を前記各時刻エミュレーションモデルに入力し、前記所定の複数の時刻の各々に対する、前記複数のサンプル点の各々に対応するエミュレーション値を計算するステップと、前記所定の複数の時刻の各々に対する、前記複数のサンプル点の各々に対応するエミュレーション値に基づいて、前記複数のサンプル点の各々に対する前記時系列シミュレーション値を推定するステップとを含む方法を提供するものである。

【0011】

本発明の１つの態様は、前記方法をコンピュータに実行させるためのプログラムを提供するものである。

【0012】

本発明の１つの態様は、前記プログラムを記録したコンピュータ読み取り可能な記録媒体を提供するものである。

【発明の効果】

【0013】

上記構成を有する本発明によれば、少ないデータ数しか事前収集できない、もしくは、限られた計算しか実施できない場合でも、予測精度が高く、不確定因子の網羅的検討を実施することができる、シミュレーションモデルの時系列シミュレーション値の推定のためのシステム、方法を提供することができる。

【図面の簡単な説明】

【0014】

【図1】本発明の１つの実施形態に係る時系列シミュレーション値推定システムの全体構成図である。

【図2】本発明の１つの実施形態に係る時系列シミュレーション値推定システムのハードウエア構成の例を示す図である。

【図3A】本発明の１つの実施形態に係る時系列シミュレーション値推定処理の一例のフローチャートである。

【図3B】本発明の１つの実施形態に係る時系列シミュレーション値推定処理の一例のフローチャートである。

【図4】人工衛星モデルの斜視図である。

【図5】人工衛星モデルの分解斜視図である。

【図6A】各説明変数が公称値で不確実性がなく、発熱機器Ａを起動し続けた場合の発熱機器Ａのシミュレーションによる時系列温度の一例を示す図である。

【図6B】各説明変数が公称値で不確実性がなく、発熱機器Ｄを起動し続けた場合の発熱機器Ｄのシミュレーションによる時系列温度の一例を示す図である。

【図7A】発熱機器Ａを起動し続けた場合の説明変数の不確実性を考慮した発熱機器Ａのシミュレーションによる時系列温度の一例を示す図である。

【図7B】発熱機器Ｄを起動し続けた場合の説明変数の不確実性を考慮した発熱機器Ｄのシミュレーションによる時系列温度の一例を示す図である。

【図8A】同一のサンプル点に対する時刻３０００秒時の発熱機器Ａのシミュレーション値とエミュレーション値を示す図である。

【図8B】同一のサンプル点に対する時刻６０００秒時の発熱機器Ａのシミュレーション値とエミュレーション値を示す図である。

【図8C】同一のサンプル点に対する時刻３０００秒時の発熱機器Ｄのシミュレーション値とエミュレーション値を示す図である。

【図8D】同一のサンプル点に対する時刻６０００秒時の発熱機器Ｄのシミュレーション値とエミュレーション値を示す図である。

【図9A】同一のサンプル点に対する時刻３０００秒時の発熱機器Ａのシミュレーション値とエミュレーション値の差のヒストグラムを示す図である。

【図9B】同一のサンプル点に対する時刻６０００秒時の発熱機器Ａのシミュレーション値とエミュレーション値の差のヒストグラムを示す図である。

【図9C】同一のサンプル点に対する時刻３０００秒時の発熱機器Ｄのシミュレーション値とエミュレーション値の差のヒストグラムを示す図である。

【図9D】同一のサンプル点に対する時刻６０００秒時の発熱機器Ｄのシミュレーション値とエミュレーション値の差のヒストグラムを示す図である。

【図10A】各発熱機器の発熱時間の合計を最大とする発熱パターンにおける発熱機器Ａの時系列温度のエミュレーション値とシミュレーション値を示す図である。

【図10B】各発熱機器の発熱時間の合計を最大とする発熱パターンにおける発熱機器Ｄの時系列温度のエミュレーション値とシミュレーション値を示す図である。

【発明を実施するための形態】

【0015】

以下、本発明の実施形態について図面を参照して説明する。本実施形態では、不確実性を有する複数の説明変数を含む宇宙機の熱数学モデルの時系列シミュレーション値の推定値であるエミュレーション値を計算する場合を例として説明する。

【0016】

図１は、本発明の１つの実施形態に係る時系列シミュレーション値推定システム１０の全体構成図である。時系列シミュレーション値推定システム１０は、モデル設定部１０１、初期設定部１０３、サンプル点生成部１０５、シミュレーション部１０７、各時刻エミュレーションモデル生成部１１１、エミュレーション部１１３を備える。

【0017】

モデル設定部１０１は、ユーザからの入力等によって、目的変数と、不確実性を有する複数の説明変数を含むシミュレーションモデルである熱数学モデルを設定する。

【0018】

初期設定部１０３は、モデル設定部１０１によって設定された熱数学モデルの、不確実性を有する複数の説明変数の各々の不確実性を設定する。また、熱数学モデルの不確実性を有さない説明変数の値を設定する。

【0019】

サンプル点生成部１０５は、１つ以上の説明変数から構成されるサンプル点を生成する。すなわち、初期設定部１０３によって設定された複数の説明変数の各々の不確実性にしたがって、複数の説明変数の、複数の初期サンプル点や、後述の生成された各時刻エミュレーションモデルに基づく時系列温度のシミュレーション値を推定するための複数のサンプル点を生成する。

【0020】

シミュレーション部１０７は、サンプル点生成部１０５によって生成された複数の初期サンプル点の値を、シミュレーションモデルである宇宙機の熱数学モデルに入力し、所定の複数の時刻の各々に対する、複数の初期サンプル点の各々に対応するシミュレーション値を計算する。

【0021】

各時刻エミュレーションモデル生成部１１１は、複数の初期サンプル点の各々と、所定の複数の時刻の各々に対する複数の初期サンプル点の各々に対応するシミュレーション値から、所定の複数の時刻の各々に対する各時刻エミュレーションモデルを生成する。

【0022】

エミュレーション部１１３は、サンプル点生成部１０５によって生成された複数のサンプル点を、各時刻エミュレーションモデル生成部１１１によって生成された各時刻エミュレーションモデルに入力し、所定の複数の時刻の各々に対する、該複数のサンプル点の各々に対応するエミュレーション値を計算する。そして、所定の複数の時刻の各々に対する、複数のサンプル点の各々に対応するエミュレーション値に基づいて、複数のサンプル点の各々に対する時系列シミュレーション値を推定する。

【0023】

図２は、本実施形態に係る時系列シミュレーション値推定システム１０のハードウエア構成の例を示す図である。時系列シミュレーション値推定システム１０は、ＣＰＵ１３０ａ、ＲＡＭ１３０ｂ、ＲＯＭ１３０ｃ、外部メモリ１３０ｄ、入力部１３０ｅ、出力部１３０ｆ、通信部１３０ｇを含む。ＲＡＭ１３０ｂ、ＲＯＭ１３０ｃ、外部メモリ１３０ｄ、入力部１３０ｅ、出力部１３０ｆ、通信部１３０ｇは、システムバス１３０ｈを介して、ＣＰＵ１３０ａに接続されている。

【0024】

図１に示される時系列シミュレーション値推定システム１０の各部は、ＲＯＭ１３０ｃや外部メモリ１３０ｄに記憶された各種プログラムが、ＣＰＵ１３０ａ、ＲＡＭ１３０ｂ、ＲＯＭ１３０ｃ、外部メモリ１３０ｄ、入力部１３０ｅ、出力部１３０ｆ、通信部１３０ｇ等を資源として使用することで実現される。

【0025】

以上のシステム構成を前提に、本発明の第１つの実施形態に係る時系列シミュレーション値推定処理の例を、図１、３Ａ、３Ｂ等を参照して、以下に説明する。図３Ａ、３Ｂは、本実施形態に係る時系列シミュレーション値推定システムの時系列シミュレーション値推定処理のフローチャートである。

【0026】

モデル設定部１０１は、ユーザからの入力等によって、目的変数と不確実性を有する複数（ｎ個）の説明変数を含むシミュレーションモデルである宇宙機の熱数学モデルを設定する（Ｓ１０１）。

【0027】

初期設定部１０３は、複数の説明変数の各々の不確実性を設定する（Ｓ１０３）。本実施形態において、説明変数の不確実性は、説明変数の公称値に対する上限・下限値である。説明変数の不確実性はこれに限定されるものではなく、公称値に対する標準偏差等の説明変数の確率分布等他の任意の適切なものとすることができる。

【0028】

サンプル点生成部１０５は、設定された各説明変数の不確実性にしたがって、高次元不確実性空間中を均一に選択されるように、LHS(Latin Hypercube Sampling)法により、Ｌ個の初期サンプル点を生成し、記憶部１１９に記憶する（Ｓ１０５）。初期サンプル点を生成する手法は、これに限定されるものではなく、一様乱数、正規乱数等の他の任意の適切な手法を用いることができる。

【0029】

シミュレーション部１０７は、生成されたＬ個の初期サンプル点を、シミュレーションモデルである宇宙機の熱数学モデルに入力して、Ｌ個の目的変数の時系列シミュレーション値ｆ₁(ｔ)，ｆ₂(ｔ)，・・・，ｆ_L(ｔ)を算出し、記憶部１１９に記憶する（Ｓ１０７）。具体的には、例えば、時刻ｔ₀のシミュレーション値をまず計算し、そのシミュレーション値に基づいて、所定時間Δｔ後のシミュレーション値を計算し、これを所定時間間隔Δｔごとに繰り返し、ｆ₁(ｔ₀＋Δｔ)、・・・、ｆ₁(ｔ₀＋ｉΔｔ)（ｉは整数）を計算する。そして、時刻ｔ₀＋ｊΔｔ（ｊは１以上（ｉ−１）以下の整数）と時刻ｔ₀＋(ｊ＋１)Δｔの間の時刻のｆ₁(ｔ)の値は、ｆ₁(ｔ₀＋ｊΔｔ)の値とｆ₁(ｔ₀＋(ｊ＋１)Δｔ)の値から補間により計算する。これらの計算により得られた時系列シミュレーション値を記憶部１１９に記憶する。

【0030】

次に、各時刻エミュレーションモデル生成部１１１は、記憶部１１９からＬ個の初期サンプル点とそれに対応する時刻ｔ₀のシミュレーション値を読み込む（Ｓ１０９）。

【0031】

各時刻エミュレーションモデル生成部１１１は、読み込まれたＬ個の初期サンプル点、すなわちＬ個のｎ次元の入力ベクトルとそれに対応する時刻ｔ₀のシミュレーション値から、時刻ｔ₀のエミュレーションモデル

を生成する（Ｓ１１１）。エミュレーションモデルの生成に用いる入力ベクトルの数や次元はＬ個やｎ次元でなくても、Ｌ個未満やｎ未満の次元の適切な数や次元としてもよい。本実施形態では、エミュレーションモデル

は、式（１）で表されるガウス過程回帰モデルであるが、これに限定されるものではなく、単項式回帰モデル、多項式回帰モデル、サポートベクトル回帰モデル、ニューラルネットワーク回帰モデル等の他の任意の適切な回帰モデルとすることができる。

ここで、

は、ｎ個の説明変数を成分とする初期サンプル点、すなわちｎ次元の入力ベクトル

である。

は、応答出力でスカラー量である。μは、定数を仮定した大域モデルである。

は、

におけるμからの偏差に相当する局所モデルである。

は、以下の条件でガウス過程に従うと仮定する。

ここで、

は、

の平均である。

は、

の分散である。

は、

と

の共分散である。

は、任意の２つの位置

と

の間の関係性を示すカーネル関数である。

【0032】

本実施形態においては、カーネル関数について、定常過程を仮定する。各時刻エミュレーションモデル生成部１１１は、読み込まれたＬ個のｎ次元の入力ベクトルとそれに対応するシミュレーション値から、線形関数カーネル、指数関数カーネル、動径基底関数（ＲＢＦ）カーネル、Matern32カーネル、Matern52等のカーネル関数の各々について、カーネル関数のハイパーパラメータ（定数、分散、長さスケール等）を、準ニュートン法や、勾配法であるＢＦＧＳ法を用いて最適化する。次いで、各時刻エミュレーションモデル生成部１１１は、最適化された各カーネル関数の中から、k-fold交差検定により選択する。候補となるカーネル関数は、上述の５つのカーネル関数に限定されるものではなく、シグモイドカーネル等の他の適切な任意のカーネル関数を用いることができる。また、カーネル関数の選択手法は、k-fold交差検定に限定されるものではなく、leave-one-out等の他の交差検定や、各カーネル関数を当てはめたときに得られる尤度の比較等の他の適切な任意の手法を用いることができる。

【0033】

今、Ｍ個のサンプル点を

、これに対応する出力を

としたとき、

におけるμとσ_z²の値は式（５）、（６）を用いて予測される。

ここで、

は、ｉ成分がｋ(ｘ⁽ⁱ⁾)であるＭ次元のベクトルである。

は、(ｉ,ｊ)成分がｋ(ｘ⁽ⁱ⁾, ｘ^(j))であるＭ×Ｍ行列である。ｋ_**は、(ｉ,ｊ)成分がｋ(ｘ_i, ｘ_j)であるｎ×ｎ行列である。雑音の分散σ_n²は、ゼロと仮定する。

【0034】

このようにして、各時刻エミュレーションモデル生成部１１１によって、時刻ｔ₀のガウス過程回帰モデルが生成される。

【0035】

続いて、各時刻エミュレーションモデル生成部１１１は、同様にして、記憶部１１９からＬ個の初期サンプル点とそれに対応する所望の時刻ｔ₁のシミュレーション値を読み込み、読み込まれたＬ個の初期サンプル点とそれに対応する時刻ｔ₁のシミュレーション値から、時刻ｔ₁のエミュレーションモデル

を生成する。そして、同様にして、所望の時刻ｔ₂、・・・、ｔ_m（ｍは整数）についてのエミュレーションモデル

、・・・、

を生成する（Ｓ１１３）。

【0036】

次に、サンプル点生成部１０５は、初期サンプル点の生成と同様にして、生成された各時刻エミュレーションモデルに基づく時系列温度のシミュレーション値を推定するための、初期サンプル点以外の、説明変数から構成されるサンプル点を多数生成する（Ｓ１１５）。

【0037】

次に、エミュレーション部１１３は、生成された多数のサンプル点の入力ベクトルを、生成された時刻ｔ₁、ｔ₂、・・・、ｔ_mについてのガウス過程回帰モデル

、

、・・・、

に代入し、多数のサンプル点の各々に対応する、時刻ｔ₁、ｔ₂、・・・、ｔ_mについての目的変数の値である応答出力ｙ₁、ｙ₂、・・・、ｙ_mの値（エミュレーション値）を計算することによって、時刻ｔ₁、ｔ₂、・・・、ｔ_mに対するシミュレーション値を推定するまた、時刻ｔ_kとｔ_(k+1)（ｋは、１以上ｍ以下の整数）の間の時刻についてのエミュレーション値（シミュレーション値の推定値）は、時刻ｔ₁、ｔ₂、・・・、ｔ_mについてのエミュレーション値で補間を行うことにより計算する（Ｓ１１７）。ここで、時刻ｔ₁、ｔ₂、・・・、ｔ_mは、シミュレーション値が補間によらずに計算された時刻とすると、精度が高くなる。

【0038】

＜実施例＞
本発明の実施例として、目的変数を人工衛星に搭載された発熱機器の温度とする人工衛星の熱数学モデルについて行った計算実験結果について説明する。

【0039】

図４は、実施例の熱数学モデルで用いた人工衛星モデルの斜視図である。また、図５は、人工衛星モデルの分解斜視図である。人工衛星２０は、１辺が５００ｍｍの立方体の形状を有し、立方体の各面が、１辺が５００ｍｍの正方形のアルミ製のハニカムパネルである＋Ｘパネル２０１、−Ｘパネル２０３、＋Ｙパネル２０５、−Ｙパネル２０７、＋Ｚパネル２０９、−Ｚパネル２１１から構成されていた。また、人工衛星２０の内部には、３枚の機器搭載用のアルミ製の矩形のデッキである＋Ｘデッキ２１３、−Ｘデッキ２１５、＋Ｚデッキ２１７が設けられていた。各デッキは、一方の長辺を互いに共有するように接続されていた。また、＋Ｘデッキ２１３、−Ｘデッキ２１５の長辺の他方は、−Ｚパネル２１１の対向するＹ方向に延びる辺に接するように接続されていた。＋Ｚデッキ２１７の長辺の他方は、＋Ｚパネル２０９のＸ方向の中央部分に接続されていた。また、＋Ｘパネル２０１、−Ｘパネル２０３、＋Ｙパネル２０５、−Ｙパネル２０７、−Ｚパネル２１１の外側表面には、ソーラーパネルが、＋Ｚパネル２０９の外側表面には、多層断熱材（ＭＬＩ）が取り付けられていた。各デッキの外側及び内側表面と、全てのパネルの内側表面は、黒色塗装が施されていた。

【0040】

−Ｘパネル２０３、＋Ｙパネル２０５、−Ｙパネル２０７、＋Ｘパネル２０１の裏面には、発熱機器Ａ２２１、発熱機器Ｂ２２３、発熱機器Ｃ２２５、発熱機器Ｄ２２７がそれぞれ設けられていた。各発熱機器の形状は、１５０ｍｍ×２００ｍｍ×１００ｍｍの直方体形状であった。発熱機器Ａ２２１、発熱機器Ｂ２２３、発熱機器Ｃ２２５、発熱機器Ｄ２２７の発熱量は、それぞれ、スイッチＯＮ時に５０、６０、７０、４０Ｗであった。各発熱機器と各パネル間の接触熱コンダクタンスは３０Ｗ／Ｋ（１０００Ｗ／ｍ²・Ｋ）であった。

【0041】

人工衛星２０の重さは、５０ｋｇであった。人工衛星２０は、地上から７９８ｋｍの太陽同期軌道上で、軌道に対して約９８度傾いていた。

【0042】

人工衛星の熱数学モデルは、人工衛星を、構成部品である、構体パネルや搭載機器などをいくつかの要素に分割し、各要素単位に熱特性（温度、比熱、熱伝導係数、輻射特性など）を代表する節点を設けることで構築された。太陽輻射、アルベド、地球赤外放射などの外部からの熱入力源や搭載機器からの発熱などによる内部熱入力もそれぞれ節点として考えることができ、これら節点間の熱交換を記述することで支配方程式が求められた。シミュレーションモデルである熱数学モデルは、式（７）で表された。

ここで、Ｃ_i、Ｔ_i、Ｑ_iは節点ｉの熱容量[J/K]、温度[K]、内外の熱入力[W] であった。Ｃ_ijは節点i、ｊ間の熱コンダクタンス[W/K]、Ｒ_ijは輻射係数[m²]、σはStefan-Boltzmann係数（5.669×10^-8[W/m²/K⁴]）であった。Ｎは総節点数であり、Ｎ個の支配方程式を連立させて解くことで各節点での温度を求めることができた。本計算実験においては、接点は、各パネル及び各デッキで、総節点数Ｎは９であった。熱コンダクタンスは節点ｉ、ｊが同一物体内の場合は物体の熱伝導率で表された。一方、節点ｉ、ｊが異種物体である場合は、接触熱伝達率で表された。

【0043】

人工衛星２０の不確実性を有する説明変数である軌道・熱環境条件は、人工衛星２０の軌道面と太陽方向とのなす角度であるβ角、太陽光強度、地球全体の反射率であるアルベド、地球表面温度からの輻射エネルギーである地球赤外輻射強度の計４個であった。各説明変数の公称値、下限値、上限値は表１の通りであった。

【0044】

表１

【0045】

上述のように、目的変数は、人工衛星に搭載された発熱機器の温度としたが、図６Ａ、図６Ｂは、各説明変数が公称値で不確実性がなく、発熱機器Ａ、Ｄを起動し続けた場合の発熱機器Ａ、Ｄのシミュレーションによる時系列温度の一例を示す図である。図７Ａ、図７Ｄは、発熱機器Ａ、Ｄを起動し続けた場合の説明変数の不確実性を考慮した発熱機器Ａ、Ｄのシミュレーションによる時系列温度の一例を示す図である。ここで、発熱機器Ｂ、Ｃについては、発熱機器Ａ、Ｄと同様の挙動を示しているので図示を省略する。

【0046】

シミュレーションは、C&R Technologies Company製のソフトウエア「SINDA/FLUINT」により行った。

【0047】

各時刻エミュレーションモデルを生成するための初期サンプル点は、４８点生成した。

【0048】

各時刻エミュレーションモデルの生成に適用されたカーネル関数は、線形関数カーネル、指数関数カーネル、動径基底関数（ＲＢＦ）カーネル、Matern32カーネル、Matern52の５つで、それぞれのカーネル関数内に存在するハイパーパラメータは、勾配法の１つの前処理付き共役勾配法（Scaled Conjugate Gradient method）により最適化された。ハイパーパラメータのうちの長さスケールについては、説明変数の次元毎で異なる値を割り当てた。カーネル関数は、ｋ＝１０のk-fold交差検定により２乗平均平方根誤差値が最小となるカーネル関数を選択した。

【0049】

時系列シミュレーション値は、シミュレーションモデルに基づいて１０秒間隔で６０００秒時まで計算され、それらの間は線形補間により計算された。また、時系列エミュレーション値は、各時刻エミュレーションモデルが３００秒間隔で６０００秒時まで生成され、生成された各時刻エミュレーションモデルに基づいて計算され、それらの間は線形補間により計算された。

【0050】

（実施例１）
生成された各時刻エミュレーションモデルに基づく時系列温度のシミュレーション値を推定するためのサンプル点を、表１に示される不確実性の範囲から、ＬＨＳ法を用いて１０００点生成した。

【0051】

図８Ａ〜図８Ｄは、生成された各時刻エミュレーションモデルに基づく時系列温度のエミュレーション値を推定するために生成したサンプル点１０００点についての、時刻が３０００秒時と６０００秒時の、同一のサンプル点に対する発熱機器Ａ、Ｄの時系列温度のシミュレーション値とエミュレーション値を示す図であり、図９Ａ〜図９Ｄは、同一のサンプル点に対する発熱機器Ａ、Ｄの時系列温度のシミュレーション値とエミュレーション値の差のヒストグラムを示す図である。また、表２に、同一のサンプル点に対する発熱機器Ａ、Ｄの時系列温度のシミュレーション値とエミュレーション値の差の平均と標準偏差を示す。

【0052】

表２

【0053】

図８Ａ〜図９Ｄ、表２から、本実施形態による各時刻エミュレーションモデルに基づくエミュレーション値は、９５％信頼区間（平均±２σ）で考えると、最大でも±４℃程度の誤差でシミュレーションモデルによるシミュレーション値を推定できていることが確認でき、本実施形態による各時刻エミュレーションモデルに基づくエミュレーション値の信頼性が確認できた。

【0054】

（実施例２）
各発熱機器の温度が許容温度上限以下となるという条件の下で、各発熱機器の発熱時間の合計を最大とする、発熱機器の発熱パターン（以下、「最適発熱パターン」という。）を抽出した。

【0055】

具体的には、発熱機器の発熱時間のパターンを１０００通り生成した。また、軌道・熱環境条件の不確実性として１００通りの軌道・熱環境条件を考慮するために、生成された各時刻エミュレーションモデルに基づく時系列温度のエミュレーション値を推定するためのサンプル点を、表１に示される不確実性の範囲から、ＬＨＳ法を用いて１００点生成した。そして、１０００通りの発熱パターンに対して、軌道・熱環境条件の不確実性として１００通りの軌道・熱環境条件を考慮し、合計で、１０００００通りの計算条件を設定した。

【0056】

この１０００００通りの計算条件について、実施例１で初期サンプル点４８点に基づいて生成された各時刻エミュレーションモデルに基づいて各発熱機器の時系列温度のエミュレーション値を計算した。

【0057】

各発熱機器の許容温度上限を５０℃として、各発熱機器の時系列温度のエミュレーション値が５０℃以下となる、発熱機器の発熱パターンを抽出し、抽出された発熱機器の発熱パターンのうちから、最適発熱パターンを抽出した。表３に、各発熱機器の発熱時間の合計が最も大きかったパターン（最適発熱パターン）から５番目に大きかったパターンを示す。

【0058】

表３

【0059】

また、最適発熱パターンにおける発熱機器Ａについての時系列温度のエミュレーション値と、同じ発熱機器の発熱パターンでのシミュレーションモデル（熱数学モデル）のシミュレーションを行った結果得られた時系列温度のシミュレーション値を図１０Ａに示す。また、発熱機器Ｄについての時系列温度のエミュレーション値とシミュレーション値を図１０Ｂに示す。

【0060】

図１０Ａ、図１０Ｂから、各発熱機器の最適発熱パターンでの時系列温度の評価は、シミュレーションモデルに基づく評価、各時刻エミュレーションモデルに基づく評価共にトレンドは一致しており、また、発熱機器の許容温度上限温度とした５０℃を逸脱していなかったことも確認できた。つまり、各時刻エミュレーションモデルに基づく最適発熱パターンの抽出の妥当性が確認できた。

【0061】

また、最適発熱パターンの探索時間についてみると、シミュレーションモデル（熱数学モデル）に基づく探索時間が約１２日であったのに対して、各時刻エミュレーションモデルに基づく探索時間は約１．５時間と、各時刻エミュレーションモデルに基づく探索がシミュレーションモデルに基づく探索よりも格段に短い時間で可能であることが分かった。

【0062】

本実施形態によれば、データ間の時系列情報を捨てて、時刻毎のデータごとにエミュレーションモデルを作成し、そのエミュレーションモデルのエミュレーション値をつなぎ合わせることで、時系列情報を再構成するという従来にない新規な発想によって、シミュレーションモデルの時系列シミュレーション値の推定に係る計算コストや事前に必要なデータ量を大幅に減少させることができ、シミュレーションモデルの時系列シミュレーション値を精度よく推定することができる。

【0063】

また、本実施形態によれば、システムの安全性評価を、不確定性因子の網羅的な検討で実現できることになり、過剰なマージン設定の防止等システムの性能を向上させることができる。

【0064】

上記実施形態においては、宇宙機の熱数学モデルを例として説明したが、宇宙機の熱数学モデルに限定されることなく、他の設計モデルや、設計モデル以外の例えば気象モデル等の一般のモデルにも本発明は適用可能である。

【0065】

以上、本発明について、例示のためにいくつかの実施形態に関して説明してきたが、本発明はこれに限定されるものでなく、本発明の範囲及び精神から逸脱することなく、形態及び詳細について、様々な変形及び修正を行うことができることは、当業者に明らかであろう。

【符号の説明】

【0066】

１０ 時系列シミュレーション値推定システム
１０１ モデル設定部
１０３ 初期設定部
１０５ サンプル点生成部
１０７ シミュレーション部
１１１ 各時刻エミュレーションモデル生成部
１１３ エミュレーション部
１１９ 記憶部
２０ 人工衛星
２０１ ＋Ｘパネル
２０３ −Ｘパネル
２０５ ＋Ｙパネル
２０７ −Ｙパネル
２０９ ＋Ｚパネル
２１１ −Ｚパネル
２１３ ＋Ｘデッキ
２１５ −Ｘデッキ
２１７ ＋Ｚデッキ
２２１ 発熱機器Ａ
２２３ 発熱機器Ｂ
２２５ 発熱機器Ｃ
２２７ 発熱機器Ｄ

【図1】

【図2】

【図3A】

【図3B】

【図4】

【図5】

【図6A】

【図6B】

【図7A】

【図7B】

【図8A】

【図8B】

【図8C】

【図8D】

【図9A】

【図9B】

【図9C】

【図9D】

【図10A】

【図10B】