JP7051771B2 - 計算装置、計算方法およびプログラム - Google Patents

Description

本発明の実施形態は、計算装置、計算方法およびプログラムに関する。

変数が離散値（２値）で表された最適化問題を解くための様々なアルゴリズムおよびハードウェアが提案されている。最適化問題を解くためのアルゴリズムおよびハードウェアは、簡易な構成および処理であることが望ましい。

また、変数が離散値および連続値の両方で表される最適化問題もある。例えば、金融工学の最も基本的な問題として、２次計画問題により定式化されるポートフォリオ最適化問題が知られている。２次計画問題は、変数が連続値で表された連続最適化問題である。

特許第５８６５４５６号公報 特開２０１８－００５５４１号公報 特開２０１８－０１０４７４号公報 特許第５８６８０４８号公報

Hayato Goto, Kosuke Tatsumura, Alexander R. Dixon, "Combinatorial optimization by simulating adiabatic bifurcations in nonlinear Hamiltonian systems",Science Advances, Vol. 5, no. 4, eaav2372, 19 Apr. 2019 Brendan O’Donoghue and Chris J. Maddison, "Hamiltonian descent for composite objectives", arXiv:1906.02608v1, Submitted on 6 Jun 2019 Chris J. Maddison et al.,"Hamiltonian Descent Methods", arXiv:1809.05042v1, Submitted on 13 Sep 2018 T. Inagaki et al.,"A coherent Ising machine for 2000-node optimization problems", Science 354, 603 (2016). H. Goto,"Bifurcation-based adiabatic quantum computation with a nonlinear oscillator network", Sci. Rep. 6, 21686 (2016). Y. Haribara et al., "Performance evaluation of coherent Ising machines against classical neural networks", Quantum Sci. Technol. 2, 044002 (2017).

発明が解決しようとする課題は、変数が連続値で表された最適化問題を簡易に解くことである。

実施形態に係る計算装置は、変数算出部と、時間発展部と、を備える。前記変数算出部は、第１時刻におけるＮ個（Ｎは２以上の整数）の第１変数を相互作用させたＮ個の第１中間変数と、前記第１時刻におけるＮ個の第２変数とに基づき、前記第１時刻からサンプリング期間を経過した第２時刻における前記Ｎ個の第２変数を算出する。前記時間発展部は、時間発展処理を実行して、前記第２時刻における前記Ｎ個の第１変数を生成する。前記時間発展部は、前記第２時刻における前記Ｎ個の第２変数と、前記第１時刻における前記Ｎ個の第１変数とに基づき、前記第２時刻における前記Ｎ個の第１変数を算出する。前記時間発展部は、前記第２時刻における前記Ｎ個の第１変数が予め定められた制約条件を満たしていない場合、前記第２時刻における前記Ｎ個の第２変数を前記制約条件を満たす方向に変化させる。

第１実施形態に係る計算装置の構成図。 第１実施形態に係る演算部の構成図。 第１実施形態に係る演算部の処理の流れを示すフローチャート。 制約条件を満たす方向に第１変数を変化させる処理を示す図。 第１実施形態に係る後処理部の処理の流れを示すフローチャート。 第２実施形態に係る時間発展部の構成図。 第２実施形態に係る演算部の処理の流れを示すフローチャート。 第３実施形態に係るｍ段目の発展部の構成図。 第３実施形態に係る演算部の処理の流れを示すフローチャート。 第４実施形態に係るｍ段目の発展部の構成図。 第４実施形態に係る演算部の処理の流れを示すフローチャート。 第４実施形態に係る後処理部の処理の流れを示すフローチャート。 計算装置のハードウェア構成図。

以下、図面を参照しながら実施形態に係る計算装置１０について詳細に説明する。実施形態に係る計算装置１０は、変数が連続値で表された最適化問題を解くことを目的とする。

（最適化問題）
まず、計算装置１０により解かれる最適化問題について説明する。計算装置１０により解かれる最適化問題は、Ｋ個（Ｋは、１以上の整数）の不等式による制約を含む２次計画問題である。

計算装置１０により解かれる最適化問題の目的関数（Ｅ（ｘ））は、式（１）により表される。

ｘは、変数を表す。Ｎは、変数の個数を表し、２以上の整数である。ｉおよびｊは、変数の識別番号を表し、１以上Ｎ以下の任意の整数である。

Ａ_ｉｊは、ｉ番目の変数ｘ_ｉと、ｊ番目の変数ｘ_ｊとの間の相互作用を表す。Ａ_ｉｊは、Ｎ行Ｎ列の行列である係数行列Ａに含まれるｉ行ｊ列の係数である。係数行列Ａは、半正定値対称行列である。

さらに、計算装置１０により解かれる最適化問題は、Ｋ個の不等式による制約を含む。計算装置１０は、Ｋ個の不等式を満たしながら最適化問題を解く。

Ｋ個の不等式は、一例として、Ｋ個の不等式を含む連立一次不等式により表される。例えば、Ｋ個の不等式のうちのｋ番目（ｋは、１以上Ｋ以下の整数）の不等式は、式（２）のように表される。

Ｂ_ｉｋは、Ｎ行Ｋ列の制約行列Ｂに含まれるｉ行ｋ列の係数である。Ｂ_ｉｋは、ｋ番目の不等式に含まれるｉ番目の変数ｘ_ｉに対する係数を表す。

さらに、計算装置１０により解かれる最適化問題は、Ｎ個の変数ｘ_ｉのそれぞれを０以上とする制約を含んでもよい。この場合、計算装置１０は、１からＮまでの全てのｉについて、式（３）により表される式を満たすように、最適化問題を解く。

さらに、計算装置１０による解かれる最適化問題は、Ｎ個の変数ｘ_ｉの合計を１に規格化する制約を含んでもよい。この場合、計算装置１０は、式（４）により表される式を満たすように、最適化問題を解く。

（ハミルトニアンＨ）
本発明者は、上述の最適化問題を、Ｎ次元空間を運動する粒子のエネルギーを表すハミルトニアンＨに組み込むことを考えた。

式（５）は、最適化問題を組み込んだハミルトニアンＨの一例を示す。

式（５）のハミルトニアンＨは、式（１）の目的関数、および、式（２）のＫ個の不等式による制約が、ポテンシャルエネルギーとして組み込まれている。

式（５）において、右辺の第１項は、粒子の運動エネルギーを表す。式（５）において、ｙ_ｉは、粒子の運動量の第ｉ成分を表す。

式（５）において、右辺の第２項は、粒子のポテンシャルエネルギーを表す。式（５）において、ｘ_ｊは、粒子の位置の第ｊ成分を表す。

右辺の第２項のポテンシャルエネルギーは、目的関数により表され、下に凸の２次関数であり、ｘ_ｉ＝０で極小値を取る。なお、第２項の係数である１／２は、係数行列Ａに含まれるＮ×Ｎ個の係数に、予め乗じられていてもよい。

式（５）において、右辺の第３項は、Ｋ個の不等式の制約を取り込んだポテンシャルエネルギー表す。式（５）において、ｔは、時刻を表す。式（５）において、ｔは、開始時刻ｔ_０から終了時刻ｔ_Ｔまで増加する。

α（ｔ）は、ｔを変数とした非減少関数であり、０以上であり、少なくとも終了時刻ｔ_Ｔにおいて１となる。例えば、α（ｔ）は、開始時刻ｔ_０において０となり、終了時刻ｔ_Ｔにおいて１となる単調増加関数であってもよい。また、α（ｔ）は、ｔに関わらず、常時１であってもよい。

Φ（ξ）は、式（６）により表される関数である。式（６）において、φは、予め定められた定数である。

さらに、θ（ξ）は、式（７）により表される関数である。

右辺の第３項は、－ξ＝｛－Σ^Ｎ _ｉ＝１（（Ｂ_ｉｋｘ_ｉ）－α（ｔ））｝が０より大きい場合、すなわち、Σ^Ｎ _ｉ＝１（（Ｂ_ｉｋｘ_ｉ）－α（ｔ））が０より小さい場合、Σ^Ｎ _ｉ＝１（（Ｂ_ｉｋｘ_ｉ）－α（ｔ））を正の方向に変化させるように、ｘ_ｉを変化させるポテンシャルエネルギーを表す。なお、右辺の第３項は、－ξ＝｛－Σ^Ｎ _ｉ＝１（（Ｂ_ｉｋｘ_ｉ）－α（ｔ））｝が０以下である場合、すなわち、Σ^Ｎ _ｉ＝１（（Ｂ_ｉｋｘ_ｉ）－α（ｔ））が０以上である場合に、最小となる。

ここで、α（ｔ）は、終了時刻ｔ_Ｔにおいて１となる関数である。従って、右辺の第３項は、ハミルトニアンＨを時間発展させた場合、終了時刻ｔ_Ｔにおいて、Ｋ個の不等式による制約を満たすように、ｘ_ｉを変化させることができる。従って、式（５）に示すハミルトニアンＨは、時間発展した場合、式（２）のＫ個の不等式による制約を満たしながら、式（１）の目的関数を最小にするように、ｘ_ｉを変化させることができる。

また、式（８）は、最適化問題を組み込んだハミルトニアンＨの他の一例を示す。

式（８）のハミルトニアンＨは、式（１）の目的関数、式（２）のＫ個の不等式による制約、式（３）のＮ個の変数ｘ_ｉのそれぞれを０以上とする制約、および、式（４）のＮ個の変数ｘ_ｉの合計を１に規格化する制約が、ポテンシャルエネルギーとして組み込まれている。

式（８）において、右辺の第１項、右辺の第２項および右辺の第３項は、式（５）と同一である。

式（８）において、右辺の第４項は、Ｎ個の変数ｘ_ｉのそれぞれを０以上とする制約を取り込んだポテンシャルエネルギーを表す。

Ψ（ξ）は、式（９）により表される関数である。式（９）において、ψは、予め定められた定数である。

右辺の第４項は、－ξ＝（－ｘ_ｉ）が０より大きい場合、すなわち、ｘ_ｉが０より小さい場合、ｘ_ｉを正の方向に変化させるポテンシャルエネルギーを表す。なお、右辺の第４項は、－ξ＝（－ｘ_ｉ）が０以下である場合、すなわち、ｘ_ｉが０以上である場合に、最小となる。

式（８）において、右辺の第５項は、Ｎ個の変数ｘ_ｉの合計を１に規格化する制約を取り込んだポテンシャルエネルギー表す。βは、予め定められた定数を表す。

右辺の第５項は、｛Σ^Ｎ _ｉ＝１（ｘ_ｉ）－α（ｔ）｝^２を最小にさせるように、ｘ_ｉを変化させるポテンシャルエネルギーを表す。なお、右辺の第５項は、｛Σ^Ｎ _ｉ＝１（ｘ_ｉ）－α（ｔ）｝＝０である場合に、最小となる。

ここで、α（ｔ）は、終了時刻ｔ_Ｔにおいて１となる関数である。従って、右辺の第５項は、ハミルトニアンＨを時間発展させた場合、終了時刻ｔ_Ｔにおいて、Ｎ個の変数ｘ_ｉの合計を１に規格化する制約を満たすように、ｘ_ｉを変化させることができる。以上のような式（８）に示すハミルトニアンＨは、時間発展した場合、式（２）、式（３）および式（４）の制約を満たしながら、式（１）の目的関数を最小にするように、ｘ_ｉを変化させることができる。

計算装置１０は、以上のようなハミルトニアンＨを開始時刻ｔ_０から終了時刻ｔ_Ｔまで時間発展させる。そして、計算装置１０は、終了時刻ｔ_Ｔにおける第１から第ＮまでのＮ個の位置変数ｘ_ｉの値を取得し、最適化問題の解として出力する。

計算装置１０は、例えば、ハミルトニアンＨの時間発展を演算処理によりシミュレーションすることにより、終了時刻ｔ_ＴにおけるＮ個のポテンシャルｘ_ｉを算出する。計算装置１０は、ハミルトニアンＨを設定可能な他の物理系システムを時間発展させることにより、終了時刻ｔ_ＴにおけるＮ個のポテンシャルｘ_ｉを取得してもよい。

（第１実施形態）
第１実施形態に係る計算装置１０について説明する。

第１実施形態に係る計算装置１０は、式（５）に示すハミルトニアンＨを運動量および位置のそれぞれにより偏微分することにより得られる連立常微分方程式を解く。

式（５）に示すハミルトニアンＨを偏微分することにより得られる運動方程式は、式（１０）および式（１１）により表される。

第１実施形態において、Ｎは、変数の個数を表し、２以上の整数である。第１実施形態において、ｉおよびｊは、変数の識別番号を表し、１以上Ｎ以下の任意の整数である。

第１実施形態において、Ｋは、１以上の整数である。また、第１実施形態において、ｋは、１以上Ｋ以下の任意の整数である。

第１実施形態において、ｘ_ｉは、粒子の位置の第ｉ成分を表す変数に相当し、第１変数と呼ぶ。第１実施形態において、ｙ_ｉは、粒子の運動量の第ｉ成分を表す変数に相当し、第２変数と呼ぶ。

第１実施形態において、Ａ_ｉｊは、ｉ番目の第１変数ｘ_ｉと、ｊ番目の第１変数ｘ_ｊとの間の相互作用を表す。Ａ_ｉｊは、Ｎ行Ｎ列の係数行列Ａに含まれるｉ行ｊ列の係数である。係数行列Ａは、密な半正定値対象対称行列である。第１実施形態において、係数行列Ａに含まれるＮ×Ｎ個の係数は、予め設定されている。

第１実施形態において、Ｂ_ｉｋは、Ｋ個の不等式のうちのｋ番目の不等式に含まれるｉ番目の第１変数ｘ_ｉに対する係数を表す。Ｂ_ｉｋは、Ｎ行Ｋ列の制約行列Ｂに含まれるｉ行ｋ列の係数である。第１実施形態において、制約行列Ｂに含まれるＫ×Ｎ個の係数は、予め設定されている。

第１実施形態において、ｒ_ｋは、式（１２）により表される。

α（ｔ）は、ｔを変数とした非減少関数であり、０以上であり、少なくとも終了時刻ｔ_Ｔにおいて１となる。第１実施形態において、α（ｔ）は、ｔに関わらず、１である。Ｂ_ｊｋは、制約行列Ｂに含まれるｊ行ｋ列の係数である。

第１実施形態において、θ（ξ）は、式（１３）により表される。

従って、式（１１）において、θ（－ｒ_ｋ）は、ｒ_ｋが０より小さい場合、１となり、ｒ_ｋが０以上である場合、０となる。つまり、式（１１）は、任意のｉにおいて、ｒ_ｋが０より小さい場合、第２項が０以外の値となり、ｒ_ｋが０以上である場合、第２項は、０となる。

第１実施形態に係る計算装置１０は、式（１０）および式（１１）により表される連立常微分方程式を離散解法により解く。より具体的には、計算装置１０は、ｔを開始時刻ｔ_０からサンプリング期間ずつ増加させて、ｔ毎に、式（１０）および式（１１）を用いてｘ_ｉおよびｙ_ｉを更新する。そして、計算装置１０は、ｔを十分に大きくした終了時刻ｔ_Ｔにおけるｘ_ｉの値を、解として出力する。これにより、計算装置１０は、Ｎ個のｘ_ｉおよびＮ個のｙ_ｉを並行して更新しながら、解を算出することができる。

式（１０）および式（１１）において、最も計算量が大きい係数行列Ａに対する行列乗算は、式（１１）には含まれるが、式（１０）には含まれない。従って、計算装置１０は、最も計算量が大きい係数行列Ａに対する行列乗算をｙ_ｉの更新のためにのみ実行すればよく、ｘ_ｉの更新のためには実行しなくてよい。

また、ｘ_ｉの時間微分値（ｄｘ_ｉ／ｄｔ）を算出するための式（１０）は、ｙ_ｉを含むが、ｘ_ｉを含まない。また、ｙ_ｉの時間微分値（ｄｙ_ｉ／ｄｔ）を算出するための式（１１）は、ｘ_ｉを含むが、ｙ_ｉを含まない。

つまり、式（１０）および式（１１）を用いる場合、ｘ_ｉとｙ_ｉとは、互いに分離されている。従って、計算装置１０は、計算量が小さく安定な離散解法を適用して、ｘ_ｉおよびｙ_ｉを更新することが可能となる。例えば、計算装置１０は、シンプレクティック・オイラー法等を適用して、ｘ_ｉおよびｙ_ｉを更新する。従って、計算装置１０は、簡易な演算および簡易な構成で、最適化問題の解を算出することができる。

図１は、第１実施形態に係る計算装置１０の構成を示す図である。計算装置１０は、演算部２０と、入力部２２と、出力部２４と、設定部２６とを備える。

演算部２０は、ＣＰＵ（Central Processing Unit）等の１または複数のプロセッシング回路がプログラムを実行することにより実現される。また、演算部２０は、ＦＰＧＡ（Field Programmable Gate Array）、ゲートアレイまたは特定用途向け集積回路（ＡＳＩＣ）等であってもよく、さらに、プロセッサ回路を含んでもよい。

演算部２０は、時刻を表すパラメータであるｔを、開始時刻ｔ_０（例えば０）からサンプリング期間ずつ増加させる。演算部２０は、それぞれの時刻毎に、Ｎ個の第１変数ｘ_ｉおよびＮ個の第２変数ｙ_ｉを更新する。そして、演算部２０は、予め定められた時刻である終了時刻ｔ_ＴにおけるＮ個の第１変数ｘ_ｉを出力する。

入力部２２は、演算部２０による演算処理に先だって、開始時刻ｔ_０におけるＮ個の第１変数ｘ_ｉおよび開始時刻ｔ_０におけるＮ個の第２変数ｙ_ｉを取得して、演算部２０に与える。開始時刻ｔ_０におけるＮ個の第１変数ｘ_ｉおよびＮ個の第２変数ｙ_ｉは、予め設定された値であってもよいし、例えば、乱数により発生された値であってもよい。

出力部２４は、演算部２０による演算処理の終了後に、終了時刻ｔ_ＴにおけるＮ個の第１変数ｘ_ｉを取得する。そして、出力部２４は、終了時刻ｔ_ＴにおけるＮ個の第１変数ｘ_ｉの値を最適解として出力する。

設定部２６は、演算部２０による演算処理に先だって、各パラメータを演算部２０に対して設定する。

図２は、第１実施形態に係る演算部２０の構成を示す図である。演算部２０は、開始時刻ｔ_０（例えば、ｔ＝０）から終了時刻ｔ_Ｔ（例えば、ｔ＝Ｔ）までサンプリング期間の間隔で順次に時刻ｔを増加させ、それぞれの時刻ｔに対するＮ個の第１変数ｘ_ｉ（ｔ）およびＮ個の第２変数ｙ_ｉ（ｔ）を算出する。そして、演算部２０は、終了時刻ｔ_ＴにおけるＮ個の第１変数ｘ_ｉ（ｔ_Ｔ）を出力する。

演算部２０は、相互作用部３２と、第１加算部３４と、時間発展部３６と、第１変数メモリ３８と、第２変数メモリ４０と、後処理部４２と、管理部４４とを有する。

相互作用部３２は、第１変数メモリ３８から、任意の時刻を表す第１時刻ｔ_１におけるＮ個の第１変数ｘ_ｉ（ｔ_１）を取得する。相互作用部３２は、第１時刻ｔ_１におけるＮ個の第１変数ｘ_ｉ（ｔ_１）と、係数行列Ａとを行列演算することにより算出されたＮ個の値のそれぞれに、第１時間定数ｄｔの正負を反転させた値を乗じたＮ個の第１中間変数ａ_ｉを生成する。第１時間定数ｄｔは、サンプリング期間を表す値である。

相互作用部３２は、例えば、係数行列メモリ５２と、行列演算部５４と、第１乗算部５６（第１中間変数算出部）とを有する。係数行列メモリ５２は、係数行列Ａを記憶する。係数行列Ａは、ユーザ等により予め設定されている。

行列演算部５４は、第１変数メモリ３８から、第１時刻ｔ_１におけるＮ個の第１変数ｘ_ｉ（ｔ_１）を取得する。行列演算部５４は、第１時刻ｔ_１におけるＮ個の第１変数ｘ_ｉ（ｔ_１）と、係数行列メモリ５２に記憶された係数行列Ａとを行列演算することにより、Ｎ個の値ａ_ｉ´を算出する。

第１乗算部５６は、第１時間定数ｄｔの正負を反転させた値が予め設定されている。第１乗算部５６は、行列演算部５４からＮ個の値ａ_ｉ´を取得する。そして、第１乗算部５６は、Ｎ個のａ_ｉ´のそれぞれに、第１時間定数ｄｔの正負を反転させた値を乗じたＮ個の第１中間変数ａ_ｉを算出する。

第１加算部３４は、第２変数メモリ４０から、第１時刻ｔ_１におけるＮ個の第２変数ｙ_ｉ（ｔ_１）を取得する。第１加算部３４は、相互作用部３２から、Ｎ個の第１中間変数ａ_ｉを取得する。そして、第１加算部３４は、第１時刻ｔ_１におけるＮ個の第２変数ｙ_ｉ（ｔ_１）のそれぞれに、Ｎ個の第１中間変数ａ_ｉのうちの対応する第１中間変数ａ_ｉを加算することにより、第１時刻ｔ_１からサンプリング期間を経過した第２時刻ｔ_２におけるＮ個の第２変数ｙ_ｉ（ｔ_２）を算出する。

時間発展部３６は、第１加算部３４から、第２時刻ｔ_２におけるＮ個の第２変数ｙ_ｉ（ｔ_２）を取得する。時間発展部３６は、第１変数メモリ３８から、第１時刻ｔ_１におけるＮ個の第１変数ｘ_ｉ（ｔ_１）を取得する。時間発展部３６は、第２時刻ｔ_２におけるＮ個の第２変数ｙ_ｉ（ｔ_２）に基づいて、第１時刻ｔ_１におけるＮ個の第１変数ｘ_ｉ（ｔ_１）を、第２時刻ｔ_２におけるＮ個の第１変数ｘ_ｉ（ｔ_２）に時間発展させる。

さらに、時間発展部３６は、第２時刻ｔ_２におけるＮ個の第１変数ｘ_ｉ（ｔ_２）が予め設定された制約条件を満たしていない場合、第２時刻ｔ_２におけるＮ個の第１変数ｘ_ｉ（ｔ_２）を制約条件を満たす方向に変化させる。また、時間発展部３６は、第２時刻ｔ_２におけるＮ個の第１変数ｘ_ｉ（ｔ_２）が制約条件を満たしていない場合、第２時刻ｔ_２におけるＮ個の第２変数ｙ_ｉ（ｔ_２）を制約条件を満たす方向に変化させる。

そして、時間発展部３６は、制約条件を満たした、または、制約条件を満たす方向に変化させた後の第２時刻ｔ_２におけるＮ個の第１変数ｘ_ｉ（ｔ_２）を第１変数メモリ３８に記憶させる。また、時間発展部３６は、制約条件を満たした、または、制約条件を満たす方向に変化させた後の第２時刻ｔ_２におけるＮ個の第２変数ｙ_ｉ（ｔ_２）を第２変数メモリ４０に記憶させる。

時間発展部３６は、例えば、第２乗算部６２（第２中間変数算出部）と、第２加算部６４と、第１制約部６６と、第２制約部６８とを有する。

第２乗算部６２は、第１時間定数ｄｔが予め設定されている。第２乗算部６２は、第１加算部３４から、第２時刻ｔ_２におけるＮ個の第２変数ｙ_ｉ（ｔ_２）を取得する。第２乗算部６２は、第２時刻ｔ_２におけるＮ個の第２変数ｙ_ｉ（ｔ_２）のそれぞれに第１時間定数ｄｔを乗じたＮ個の第２中間変数ｂ_ｉを算出する。

第２加算部６４は、第１変数メモリ３８から、第１時刻ｔ_１におけるＮ個の第１変数ｘ_ｉ（ｔ_１）を取得する。また、第２加算部６４は、第２乗算部６２から、Ｎ個の第２中間変数ｂ_ｉを取得する。そして、第２加算部６４は、第１時刻ｔ_１におけるＮ個の第１変数ｘ_ｉ（ｔ_１）のそれぞれに、Ｎ個の第２中間変数ｂ_ｉのうちの対応する第２中間変数ｂ_ｉを加算することにより、第２時刻ｔ_２におけるＮ個の第１変数ｘ_ｉ（ｔ_２）を算出する。

第１制約部６６は、第２加算部６４から、第２時刻ｔ_２におけるＮ個の第１変数ｘ_ｉ（ｔ_２）を取得する。第１制約部６６は、第２時刻ｔ_２におけるＮ個の第１変数ｘ_ｉ（ｔ_２）が、制約条件を満たしているか否かを判断する。第１制約部６６は、第２時刻ｔ_２におけるＮ個の第１変数ｘ_ｉ（ｔ_２）が制約条件を満たしている場合、第２加算部６４から取得した第２時刻ｔ_２におけるＮ個の第１変数ｘ_ｉ（ｔ_２）を、そのまま第１変数メモリ３８に記憶させる。

第１制約部６６は、第２時刻ｔ_２におけるＮ個の第１変数ｘ_ｉ（ｔ_２）が制約条件を満たしていない場合、第２加算部６４から取得した第２時刻ｔ_２におけるＮ個の第１変数ｘ_ｉ（ｔ_２）を制約条件を満たす方向に変化させる。第１制約部６６は、第２時刻ｔ_２におけるＮ個の第１変数ｘ_ｉ（ｔ_２）が制約条件を満たしていない場合、制約条件を満たす方向に変化させた後の第２時刻ｔ_２におけるＮ個の第１変数ｘ_ｉ（ｔ_２）を第１変数メモリ３８に記憶させる。なお、第１制約部６６は、少なくとも制約条件を満たす方向に変化させれば、変化後の第２時刻ｔ_２におけるＮ個の第１変数ｘ_ｉ（ｔ_２）が制約条件を満たしていなくてもよい。

第２制約部６８は、第１加算部３４から、第２時刻ｔ_２におけるＮ個の第２変数ｙ_ｉ（ｔ_２）を取得する。第２制約部６８は、第２時刻ｔ_２におけるＮ個の第１変数ｘ_ｉ（ｔ_２）が、制約条件を満たしているか否かを判断する。なお、第２制約部６８は、第１制約部６６による判断結果に基づき、第２時刻ｔ_２におけるＮ個の第１変数ｘ_ｉ（ｔ_２）が制約条件を満たしているか否かを判断してもよい。第２制約部６８は、第２時刻ｔ_２におけるＮ個の第１変数ｘ_ｉ（ｔ_２）が制約条件を満たしている場合、第１加算部３４から取得した第２時刻ｔ_２におけるＮ個の第２変数ｙ_ｉ（ｔ_２）を、そのまま第２変数メモリ４０に記憶させる。

第２制約部６８は、第２時刻ｔ_２におけるＮ個の第１変数ｘ_ｉ（ｔ_２）が制約条件を満たしていない場合、第１加算部３４から取得した第２時刻ｔ_２におけるＮ個の第２変数ｙ_ｉ（ｔ_２）を制約条件を満たす方向に変化させる。第２制約部６８は、第２時刻ｔ_２におけるＮ個の第１変数ｘ_ｉ（ｔ_２）が制約条件を満たしていない場合、制約条件を満たす方向に変化させた後の第２時刻ｔ_２におけるＮ個の第２変数ｙ_ｉ（ｔ_２）を第２変数メモリ４０に記憶させる。

ここで、制約条件は、Ｋ個の不等式を含む連立一次不等式で表される。Ｋ個の不等式のそれぞれは、Ｎ個の第１変数ｘ_ｉのうちの１以上Ｎ以下の複数の第１変数ｘ_ｉを、変数として含む。

第１制約部６６および第２制約部６８は、Ｋ個の不等式のそれぞれ毎に、不等式を満たしているか否かを判断する。そして、Ｋ個の不等式のそれぞれ毎に、第１制約部６６は、不等式を満たしていない場合、第２時刻ｔ_２におけるＮ個の第１変数ｘ_ｉ（ｔ_２）のうちの不等式に含まれる複数の第１変数ｘ_ｉを、不等式を満たす方向に変化させる。また、Ｋ個の不等式のそれぞれ毎に、第２制約部６８は、不等式を満たしていない場合、第２時刻ｔ_２におけるＮ個の第２変数ｙ_ｉ（ｔ_２）のうちの不等式に含まれる複数の第１変数ｘ_ｉに対応する複数の第２変数ｙ_ｉを、不等式を満たす方向に変化させる。

第１変数メモリ３８は、開始時刻ｔ_０から終了時刻ｔ_Ｔまでの各時刻におけるＮ個の第１変数ｘ_ｉ（ｔ）を記憶する。なお、各時刻におけるＮ個の第１変数ｘ_ｉ（ｔ）は、次時刻におけるＮ個の第１変数ｘ_ｉ（ｔ＋ｄｔ）が算出された後、計算には用いられない。従って、第１変数メモリ３８は、各時刻におけるＮ個の第１変数ｘ_ｉ（ｔ）を記憶した後、その時刻以前のＮ個の第１変数ｘ_ｉ（ｔ）が消去されてもよい。

第２変数メモリ４０は、開始時刻ｔ_０から終了時刻ｔ_Ｔまでの各時刻におけるＮ個の第２変数ｙ_ｉ（ｔ）を記憶する。なお、各時刻におけるＮ個の第２変数ｙ_ｉ（ｔ）は、次時刻におけるＮ個の第２変数ｙ_ｉ（ｔ＋ｄｔ）が算出された後、計算には用いられない。従って、第２変数メモリ４０は、各時刻におけるＮ個の第２変数ｙ_ｉ（ｔ）を記憶した後、その時刻以前のＮ個の第２変数ｙ_ｉ（ｔ）が消去されてもよい。

後処理部４２は、終了時刻ｔ_ＴにおけるＮ個の第１変数ｘ_ｉが算出された後、終了時刻ｔ_ＴにおけるＮ個の第１変数ｘ_ｉを第１変数メモリ３８から取得する。後処理部４２は、終了時刻ｔ_ＴにおけるＮ個の第１変数ｘ_ｉが制約条件を満たしているか否かを判断する。

後処理部４２は、終了時刻ｔ_ＴにおけるＮ個の第１変数ｘ_ｉが制約条件を満たしている場合、第１変数メモリ３８から取得した終了時刻ｔ_ＴにおけるＮ個の第１変数ｘ_ｉ（ｔ_Ｔ）を、そのまま出力部２４に出力する。

後処理部４２は、終了時刻ｔ_ＴにおけるＮ個の第１変数ｘ_ｉが制約条件を満たしていない場合、第１変数メモリ３８から取得した終了時刻ｔ_ＴにおけるＮ個の第１変数ｘ_ｉ（ｔ_Ｔ）が制約条件を満たすまで、終了時刻ｔ_ＴにおけるＮ個の第１変数ｘ_ｉ（ｔ_Ｔ）を制約条件を満たす方向に変化させる。そして、後処理部４２は、制約条件を満たした終了時刻ｔ_ＴにおけるＮ個の第１変数ｘ_ｉ（ｔ_Ｔ）を出力部２４に出力する。

例えば、後処理部４２は、Ｎ個の第１変数ｘ_ｉを制約条件を満たす方向に変化させる修正処理を、Ｎ個の第１変数ｘ_ｉが制約条件を満たすまで複数回繰り返す。制約条件は、Ｋ個の不等式で表される。従って、後処理部４２は、それぞれの修正処理において、Ｋ個の不等式のそれぞれ毎に、不等式を満たしているか否かを判断し、不等式を満たしていない場合、終了時刻ｔ_ＴにおけるＮ個の第１変数ｘ_ｉ（ｔ_Ｔ）のうちの不等式に含まれる複数の第１変数ｘ_ｉを、不等式を満たす方向に変化させる。

管理部４４は、相互作用部３２、時間発展部３６および後処理部４２による処理のタイミングを管理する。管理部４４は、開始時刻ｔ_０から終了時刻ｔ_Ｔまでサンプリング期間の間隔で順次に時刻を増加させる。そして、管理部４４は、それぞれの時刻に対するＮ個の第１変数ｘ_ｉおよびＮ個の第２変数ｙ_ｉを算出させるように、相互作用部３２、時間発展部３６および後処理部４２による処理を管理する。

以上のような構成の演算部２０は、後処理部４２から出力された終了時刻ｔ_ＴにおけるＮ個の第１変数ｘ_ｉ（ｔ_Ｔ）を出力部２４に出力する。出力部２４は、演算部２０から受け取った終了時刻ｔ_ＴにおけるＮ個の第１変数ｘ_ｉ（ｔ_Ｔ）を外部装置へと出力する。

図３は、第１実施形態に係る演算部２０の処理の流れを示すフローチャートである。第１実施形態に係る演算部２０は、図３に示す流れで処理を実行する。なお、本例においては、演算部２０は、α（ｔ）を、時刻ｔに関わらず１に固定した処理を実行している。

まず、Ｓ１１において、演算部２０は、時刻ｔ、Ｎ個の第１変数ｘ_ｉおよびＮ個の第２変数ｙ_ｉを初期化する。例えば、演算部２０は、時刻ｔを０にする。また、演算部２０は、Ｎ個の第１変数ｘ_ｉのそれぞれを１／Ｎにする。また、演算部２０は、Ｎ個の第２変数ｙ_ｉのそれぞれを０にする。なお、演算部２０は、Ｎ個の第１変数ｘ_ｉのそれぞれおよびＮ個の第２変数ｙ_ｉのそれぞれを、入力部２２が外部から取得した値に初期化してもよい。

続いて、演算部２０は、時刻ｔが予め設定されたＴとなるまで、Ｓ１３からＳ２３までの処理を繰り返す（Ｓ１２とＳ２４との間のループ処理）。Ｔは、終了時刻ｔ_Ｔを表す値である。

Ｓ１３において、相互作用部３２は、相互作用演算を実行することにより、Ｎ個の第１中間変数ａ_ｉを生成する。相互作用演算として、演算部２０は、Ｎ個の第１変数ｘ_ｉと係数行列Ａとを行列演算することにより算出されたＮ個の値のそれぞれに、第１時間定数ｄｔの正負を反転させた値を乗じる。例えば、Ｓ１３において、相互作用部３２は、ｉ＝１～Ｎのそれぞれについて、式（１４）の演算を実行する。

続いて、Ｓ１４において、第１加算部３４は、Ｎ個の第２変数ｙ_ｉのそれぞれに、Ｎ個の第１中間変数ａ_ｉのうちの対応する第１中間変数ａ_ｉを加算することにより、Ｎ個の第２変数ｙ_ｉを算出する。例えば、Ｓ１４において、第１加算部３４は、ｉ＝１～Ｎのそれぞれについて、式（１５）の演算を実行する。

続いて、Ｓ１５において、時間発展部３６は、Ｎ個の第２変数ｙ_ｉのそれぞれに第１時間定数ｄｔを乗じたＮ個の第２中間変数ｂ_ｉを生成する。例えば、Ｓ１５において、時間発展部３６は、ｉ＝１～Ｎのそれぞれについて、式（１６）の演算を実行する。

続いて、Ｓ１６において、時間発展部３６は、Ｎ個の第１変数ｘ_ｉのそれぞれに、Ｎ個の第２中間変数ｂ_ｉのうちの対応する第２中間変数ｂ_ｉを加算する。例えば、Ｓ１６において、時間発展部３６は、ｉ＝１～Ｎのそれぞれについて、式（１７）の演算を実行する。

続いて、時間発展部３６は、Ｓ１８からＳ２１までの処理を、Ｋ回繰り返して実行する（Ｓ１７とＳ２２との間のループ処理）。これにより、時間発展部３６は、Ｋ個の不等式のそれぞれ毎に、不等式の判断処理および制約処理を実行することができる。なお、本フローにおいて、ｋは、ループ内の繰り返し回数を表す。

Ｓ１８において、時間発展部３６は、Ｋ個の不等式のうちのｋ番目の不等式を満たしているか否かを判断するための条件値ｒ_ｋを算出する。例えば、Ｓ１８において、時間発展部３６は、式（１８）の演算を実行する。

続いて、Ｓ１９において、時間発展部３６は、ｋ番目の不等式を満たしているか否かを判断する。より具体的には、時間発展部３６は、条件値ｒ_ｋが０より小さいか否かを判断する。ｋ番目の不等式を満たしていない場合、すなわち、条件値ｒ_ｋが０より小さい場合（Ｓ１９のＹｅｓ）、時間発展部３６は、Ｓ２０およびＳ２１の処理を実行する。

ｋ番目の制約を満たしている場合、すなわち、条件値ｒ_ｋが０より小さくない場合（Ｓ１９のＮｏ）、時間発展部３６は、Ｓ２０およびＳ２１の処理を実行しない。

Ｓ２０において、時間発展部３６は、第１制約処理を実行する。すなわち、時間発展部３６は、Ｎ個の第１変数ｘ_ｉを制約条件を満たす方向に変化させる。例えば、Ｓ２０において、時間発展部３６は、ｉ＝１～Ｎのそれぞれについて、式（１９）の演算を実行する。

Ｓ２１において、時間発展部３６は、第２制約処理を実行する。すなわち、時間発展部３６は、Ｎ個の第２変数ｙ_ｉを制約条件を満たす方向に変化させる。例えば、Ｓ２１において、時間発展部３６は、ｉ＝１～Ｎのそれぞれについて、式（２０）の演算を実行する。

時間発展部３６は、Ｓ１８からＳ２１までの処理をＫ回繰り返した場合、Ｓ１７とＳ２２との間のループ処理を抜けて、処理をＳ２３に進める。

Ｓ２３において、演算部２０は、時刻ｔを更新する。例えば、演算部２０は、時刻ｔに第１時間定数ｄｔを加算する。

Ｓ２４において、演算部２０は、時刻ｔがＴ以上であるか否かを判断する。演算部２０は、時刻ｔがＴ以上となった場合、Ｓ１２とＳ２４との間のループ処理を抜けて、処理をＳ２５に進める。

Ｓ２５において、演算部２０は、後処理を実行する。演算部２０は、後処理において、Ｎ個の第１変数ｘ_ｉが制約条件を満たしていない場合、制約条件を満たすようにＮ個の第１変数ｘ_ｉを修正する。例えば、演算部２０は、最急降下法を用いて、Ｎ個の第１変数ｘ_ｉを制約条件を満たす方向に変化させる修正処理を実行する。そして、演算部２０は、このような修正処理を、Ｎ個の第１変数ｘ_ｉが制約条件を満たすまで繰り返す。

Ｓ２５の処理を終えると、演算部２０は、本フローを終了する。

図４は、制約条件を満たす方向にＮ個の第１変数ｘ_ｉおよびＮ個の第２変数ｙ_ｉを変化させる処理について説明する図である。時間発展部３６は、Ｓ２０において、Ｎ個の第１変数ｘ_ｉを制約条件を満たす方向に変化させる。また、時間発展部３６は、Ｓ２１において、Ｎ個の第２変数ｙ_ｉを制約条件を満たす方向に変化させる。

例えば、Ｎ＝２の場合、制約条件は、例えば１次不等式で表される。この場合、制約条件を満たす領域と制約条件を満たさない領域との境界は、図３に示すような一次方程式で表される。時間発展部３６は、制約条件を満たさない場合（不等式を満たさない場合）、不等式により表される境界に対して垂直に近づける方向に、Ｎ個の第１変数ｘ_ｉにより表される位置を変化させる。また、時間発展部３６は、制約条件を満たさない場合（不等式を満たさない場合）、Ｎ個の第１変数ｘ_ｉにより表される位置が不等式により表される境界に対して近づく方向に、Ｎ個の第２変数ｙ_ｉを変化させる。

また、制約条件がＫ個の不等式を含む連立一次不等式で表される場合、時間発展部３６は、Ｋ個の不等式のそれぞれ毎に、Ｓ２０およびＳ２１の処理を実行する。すなわち、時間発展部３６は、不等式を満たしていない場合、第２時刻におけるＮ個の第１変数ｘ_ｉのうちの不等式に含まれる複数の第１変数により表される位置を、不等式により表される境界に対して垂直に近づける方向に、第１変数ｘ_ｉを変化させ、不等式に含まれる複数の第１変数に対応する複数の第２変数ｙ_ｉを変化させる。

また、後処理において、制約条件を満たすようにＮ個の第１変数ｘ_ｉを修正する場合、演算部２０は、同様の処理を実行してもよい。

図５は、第１実施形態に係る後処理部４２の処理の流れを示すフローチャートである。第１実施形態に係る後処理部４２は、図５に示す流れで処理を実行する。

まず、Ｓ３１において、後処理部４２は、Ｎ個の基準変数ｘ_ｉ´の全てを０に初期化する。Ｎ個の基準変数ｘ_ｉ´は、直前の修正処理によって変化する前のＮ個の第１変数ｘ_ｉを表す。

続いて、Ｓ３２において、後処理部４２は、直前の修正処理によりＮ個の第１変数ｘ_ｉが変化したか否かを判断する。変化した場合（Ｓ３２のＹｅｓ）、後処理部４２は、処理をＳ３３に進めて、再度、修正処理を実行する。変化しなかった場合（Ｓ３２のＮｏ）、後処理部４２は、Ｎ個の第１変数ｘ_ｉが制約条件を満たしたとして、本フローを終了する。

例えば、Ｓ３２において、後処理部４２は、式（２１）に示す判定処理を実行する。

式（２１）において、後処理部４２は、Ｎ個の基準変数ｘ_ｉ´と、修正処理の後のＮ個の第１変数ｘ_ｉとの二乗誤差を算出する。そして、後処理部４２は、二乗誤差の平方根が許容誤差εより大きくない場合（Ｓ３２のＮｏ）、Ｎ個の第１変数ｘ_ｉが制約条件を満たしたとして、本フローを終了する。また、後処理部４２は、二乗誤差の平方根が許容誤差εより大きい場合には、処理をＳ３３に進め、再度、修正処理を実行する。なお、許容誤差εは、例えば計算機イプシロンであり、１より大きく計算機で取り扱える最小の値と、１との差である。

Ｓ３３において、後処理部４２は、Ｎ個の基準変数ｘ_ｉ´に、Ｎ個の第１変数ｘ_ｉを代入する。

続いて、後処理部４２は、Ｓ３５からＳ３７までの処理を、Ｋ回繰り返して実行する（Ｓ３４とＳ３８との間のループ処理）。これにより、後処理部４２は、Ｋ個の不等式のそれぞれ毎に、不等式の判断処理および制約処理を実行することができる。

Ｓ３５において、後処理部４２は、Ｋ個の不等式のうちのｋ番目の不等式を満たしているか否かを判断するための条件値ｒ_ｋを算出する。例えば、Ｓ３５において、後処理部４２は、式（２２）の演算を実行する。

続いて、Ｓ３６において、後処理部４２は、ｋ番目の不等式を満たしているか否かを判断する。より具体的には、後処理部４２は、条件値ｒ_ｋが０より小さいか否かを判断する。ｋ番目の不等式を満たしていない場合、すなわち、条件値ｒ_ｋが０より小さい場合（Ｓ３６のＹｅｓ）、後処理部４２は、Ｓ３７の処理を実行する。

ｋ番目の制約を満たしている場合、すなわち、条件値ｒ_ｋが０より小さくない場合（Ｓ３６のＮｏ）、後処理部４２は、Ｓ３７の処理を実行しない。

Ｓ３７において、後処理部４２は、修正処理を実行する。すなわち、後処理部４２は、Ｎ個の第１変数ｘ_ｉを制約条件を満たす方向に変化させる。例えば、Ｓ３７において、後処理部４２は、ｉ＝１～Ｎのそれぞれについて、式（２３）の演算を実行する。なお、式（２３）において、γは、予め定められた定数である。

後処理部４２は、Ｓ３５からＳ３７までの処理をＫ回繰り返した場合、Ｓ３４とＳ３８との間のループ処理を抜けて、処理をＳ３２に戻す。

以上のような第１実施形態に係る計算装置１０は、Ｋ個の不等式による制約を含む２次計画問題をシミュレーションにより解くことができる。

（第２実施形態）
第２実施形態に係る計算装置１０について説明する。第２実施形態に係る計算装置１０は、第１実施形態に係る計算装置１０と略同一の構成を有するので、略同一の構成を有するブロックには、同一の符号を付けて、相違点を除き詳細な説明を省略する。

図６は、第２実施形態に係る時間発展部３６の構成を示す図である。第２実施形態に係る計算装置１０は、第１実施形態と比較して、時間発展部３６の構成が異なる。

第２実施形態に係る時間発展部３６は、第１時刻ｔ_１におけるＮ個の第１変数ｘ_ｉ（ｔ_１）に対してＭ回（Ｍは２以上の整数）の時間発展処理を実行して、第２時刻ｔ_２におけるＮ個の第１変数ｘ_ｉ（ｔ_２）を生成する。第２実施形態に係る時間発展部３６は、直列に接続されたＭ段の発展部７２－１～７２－Ｍを有する。

１段目の発展部７２－１は、第１変数メモリ３８から、第１時刻ｔ_１におけるＮ個の第１変数ｘ_ｉ（ｔ_１）を取得する。また、１段目の発展部７２－１は、第１加算部３４から、第２時刻ｔ_２におけるＮ個の第２変数ｙ_ｉ（ｔ_２）を取得する。そして、１段目の発展部７２－１は、１回目の時間発展処理を実行して、第１時刻ｔ_１から、サンプリング期間の１／Ｍの微小期間を経過した時刻（ｔ_１＋ｄｔ´）におけるＮ個の第１変数ｘ_ｉ（ｔ_１＋ｄｔ´）を出力する。また、１段目の発展部７２－１は、１回目の時間発展処理を実行して、第２時刻ｔ_２におけるＮ個の第２変数ｙ_ｉ（ｔ_２）を出力する。

ｆ段目（ｆは、２以上、Ｍ－１以下の整数）の発展部７２－ｆは、（ｆ－１）段目の発展部７２－（ｆ－１）から、第１時刻ｔ_１に、微小期間の（ｆ－１）倍の時間を経過した時刻（ｔ_１＋（ｆ－１）×ｄｔ´）おけるＮ個の第１変数ｘ_ｉ（ｔ_１＋（ｆ－１）×ｄｔ´）を取得する。また、ｆ段目の発展部７２－ｆは、（ｆ－１）段目の発展部７２－（ｆ－１）から、第２時刻ｔ_２におけるＮ個の第２変数ｙ_ｉ（ｔ_２）を取得する。そして、ｆ段目の発展部７２－１は、ｆ回目の時間発展処理を実行して、第１時刻ｔ_１から、微小期間のｆ倍の時間を経過した時刻（ｔ_１＋ｆ×ｄｔ´）におけるＮ個の第１変数ｘ_ｉ（ｔ_１＋ｆ×ｄｔ´）を出力する。また、ｆ段目の発展部７２－ｆは、ｆ回目の時間発展処理を実行して、第２時刻ｔ_２におけるＮ個の第２変数ｙ_ｉ（ｔ_２）を出力する。

Ｍ段目の発展部７２－Ｍは、（Ｍ－１）段目の発展部７２－（Ｍ－１）から、第１時刻ｔ_１に、微小期間の（Ｍ－１）倍の時間を経過した時刻（ｔ_１＋（Ｍ－１）×ｄｔ´）おけるＮ個の第１変数ｘ_ｉ（ｔ_１＋（Ｍ－１）×ｄｔ´）を取得する。また、Ｍ段目の発展部７２－Ｍは、（Ｍ－１）段目の発展部７２－（Ｍ－１）から、第２時刻ｔ_２におけるＮ個の第２変数ｙ_ｉ（ｔ_２）を取得する。そして、Ｍ段目の発展部７２－１は、Ｍ回目の時間発展処理を実行して、第２時刻ｔ_２におけるＮ個の第１変数ｘ_ｉ（ｔ_２）を出力する。また、Ｍ段目の発展部７２－Ｍは、Ｍ回目の時間発展処理を実行して、第２時刻ｔ_２におけるＮ個の第２変数ｙ_ｉ（ｔ_２）を出力する。

Ｍ段の発展部７２－１～７２－Ｍのそれぞれは、第１実施形態に係る時間発展部３６と同一の構成を含む。すなわち、Ｍ段の発展部７２－１～７２－Ｍのそれぞれは、第２乗算部６２と、第２加算部６４と、第１制約部６６と、第２制約部６８とを含む。ただし、Ｍ段の発展部７２－１～７２－Ｍのそれぞれは、第１実施形態における第１時間定数ｄｔを、微小期間を表す第２時間定数ｄｔ´に代えて処理をする。

従って、第ｍ段目（ｍは、１以上、Ｍ以下の整数）の発展部７２－ｍは、Ｍ回の時間発展処理のうちのｍ回目の時間発展処理において、次のような処理を実行する。

すなわち、第ｍ段目の発展部７２－ｍは、第（ｍ－１）回目の時間発展処理で出力した第２時刻ｔ_２におけるＮ個の第２変数ｙ_ｉ（ｔ_２）のそれぞれに、第２時間定数ｄｔ´を乗じたＮ個の第２中間変数ｂ_ｉを算出する。なお、ｍ＝１の場合、第ｍ段目の発展部７２－ｍは、第１加算部３４から取得した、第２時刻ｔ_２におけるＮ個の第２変数ｙ_ｉ（ｔ_２）のそれぞれに、第２時間定数ｄｔ´を乗じたＮ個の第２中間変数ｂ_ｉを算出する。

第ｍ段目の発展部７２－ｍは、第１時刻ｔ_１から微小期間の（ｍ－１）倍の時間を経過した時刻（ｔ_１＋（ｍ－１）×ｄｔ´）におけるＮ個の第１変数ｘ_ｉ（ｔ_１＋（ｍ－１）×ｄｔ´）のそれぞれに、Ｎ個の第２中間変数ｂ_ｉのうちの対応する第２中間変数ｂ_ｉを加算することにより、第１時刻ｔ_１から微小期間のｍ倍の時間を経過した第ｍ経過時刻（ｔ_１＋ｍ×ｄｔ´）におけるＮ個の第１変数ｘ_ｉ（ｔ_１＋ｍ×ｄｔ´）を算出する。

第ｍ段目の発展部７２－ｍは、第ｍ経過時刻（ｔ_１＋ｍ×ｄｔ´）におけるＮ個の第１変数ｘ_ｉ（ｔ_１＋ｍ×ｄｔ´）が、制約条件を満たしているか否かを判断する。第ｍ段目の発展部７２－ｍは、制約条件を満たしている場合、第ｍ経過時刻（ｔ_１＋ｍ×ｄｔ´）におけるＮ個の第１変数ｘ_ｉ（ｔ_１＋ｍ×ｄｔ´）をそのまま出力する。また、第ｍ段目の発展部７２－ｍは、制約条件を満たしている場合、取得した第２時刻ｔ_２におけるＮ個の第２変数ｙ_ｉ（ｔ_２）を、そのまま出力する。

第ｍ段目の発展部７２－ｍは、制約条件を満たしていない場合、第ｍ経過時刻（ｔ_１＋ｍ×ｄｔ´）におけるＮ個の第１変数ｘ_ｉ（ｔ_１＋ｍ×ｄｔ´）を、制約条件を満たす方向に変化させて、出力する。また、第ｍ段目の発展部７２－ｍは、制約条件を満たしていない場合、取得した第２時刻ｔ_２におけるＮ個の第２変数ｙ_ｉ（ｔ_２）を、制約条件を満たす方向に変化させて、出力する。

第ｍ段目の発展部７２－ｍは、制約条件を満たす方向に変化させる処理を、第１実施形態に係る時間発展部３６と同様に処理する。ただし、第ｍ段目の発展部７２－ｍは、第１実施形態における第１時間定数ｄｔを第２時間定数ｄｔ´に代えて処理をする。

図７は、第２実施形態に係る演算部２０の処理の流れを示すフローチャートである。第２実施形態に係る演算部２０は、図７に示す流れで処理を実行する。図７のフローチャートは、図３のフローチャートと比較して、Ｓ１４とＳ１５との間にＭ回のループの開始処理（Ｓ４１）が追加され、Ｓ２２とＳ２３との間に、Ｍ回のループの終了処理（Ｓ４２）が追加されている。

第２実施形態に係る時間発展部３６は、Ｓ１４の処理の後において、Ｓ１５からＳ２２までの処理を、Ｍ回繰り返して実行する（Ｓ４１とＳ４２との間のループ処理）。これにより、時間発展部３６は、時間発展処理をＭ回繰り返して実行することができる。

ただし、Ｓ１５からＳ２２までの処理において、時間発展部３６は、第１実施形態では第１時間定数ｄｔを用いた処理を、第２時間定数ｄｔ´を用いて行う。

より詳しくは、Ｓ１５において、時間発展部３６は、Ｎ個の第２変数ｙ_ｉのそれぞれに第２時間定数ｄｔ´を乗じたＮ個の第２中間変数ｂ_ｉを生成する。例えば、Ｓ１５において、時間発展部３６は、ｉ＝１～Ｎのそれぞれについて、式（２４）の演算を実行する。

また、例えば、Ｓ２１において、時間発展部３６は、ｉ＝１～Ｎのそれぞれについて、式（２５）の演算を実行する。

以上のような第２実施形態に係る計算装置１０は、Ｋ個の不等式による制約を含む２次計画問題をシミュレーションにより解くことができる。

（第３実施形態）
第３実施形態に係る計算装置１０について説明する。第３実施形態に係る計算装置１０は、第２実施形態に係る計算装置１０と略同一の構成を有するので、略同一の構成を有するブロックには、同一の符号を付けて、相違点を除き詳細な説明を省略する。

第３実施形態に係る計算装置１０は、式（２６）および式（２７）に示される連立常微分方程式を解く。

μは、０より大きい緩和定数である。式（２７）における、－μｙ_ｉは、運動量の緩和項を表す。

第３実施形態に係る計算装置１０は、式（２６）および式（２７）により表される連立常微分方程式を離散解法により解く。より具体的には、計算装置１０は、ｔを開始時刻ｔ_０からサンプリング期間ずつ増加させて、ｔ毎に、式（２６）および式（２７）を用いてｘ_ｉおよびｙ_ｉを更新する。そして、計算装置１０は、ｔを十分に大きくした終了時刻ｔ_Ｔにおけるｘ_ｉの値を、解として出力する。これにより、計算装置１０は、Ｎ個のｘ_ｉおよびＮ個のｙ_ｉを並行して更新しながら、解を算出することができる。

図８は、第３実施形態に係るｍ段目の発展部７２－ｍの構成を示す図である。Ｍ段の発展部７２－１～７２－Ｍのそれぞれは、緩和部７４をさらに含む。

ｍ段目の発展部７２－ｍの緩和部７４は、ｍ－１段目の発展部７２－（ｍ－１）から、第２時刻ｔ_２におけるＮ個の第２変数ｙ_ｉ（ｔ_２）を取得する。なお、ｍ＝１の場合、ｍ段目の発展部７２－ｍの緩和部７４は、第１加算部３４から、第２時刻ｔ_２におけるＮ個の第２変数ｙ_ｉ（ｔ_２）を取得する。

緩和部７４は、取得した第２時刻ｔ_２におけるＮ個の第２変数ｙ_ｉ（ｔ_２）のそれぞれに、緩和項を加算する。緩和項は、第２時間定数ｄｔ´と、予め定められた緩和定数μの正負を反転させた値と、取得した第２時刻ｔ_２におけるＮ個の第２変数ｙ_ｉ（ｔ_２）のうちの対応する第２変数ｙ_ｉ（ｔ_２）とを乗じた値である。

緩和部７４は、緩和項を加算した第２時刻ｔ_２におけるＮ個の第２変数ｙ_ｉ（ｔ_２）を、第２制約部６８に与える。第２制約部６８は、緩和項を加算した第２時刻ｔ_２におけるＮ個の第２変数ｙ_ｉ（ｔ_２）に対して制約処理を実行する。

図９は、第３実施形態に係る演算部２０の処理の流れを示すフローチャートである。第３実施形態に係る演算部２０は、図９に示す流れで処理を実行する。図９のフローチャートは、図７のフローチャートと比較して、Ｓ１６とＳ１７との間に、緩和処理（Ｓ５１）が追加されている。

第３実施形態に係る時間発展部３６は、Ｓ１６の処理の後において、Ｎ個の第２変数ｙ_ｉ（ｔ_２）のそれぞれに、緩和項を加算する緩和処理を実行する。例えば、Ｓ５１において、時間発展部３６は、ｉ＝１～Ｎのそれぞれについて、式（２８）の演算を実行する。

なお、第１実施形態に係る時間発展部３６が緩和項を加算する緩和処理を実行する構成であってもよい。この場合、緩和項は、第１時間定数ｄｔと、緩和定数μの正負を反転させた値と、取得した第２時刻ｔ_２におけるＮ個の第２変数ｙ_ｉ（ｔ_２）のうちの対応する第２変数ｙ_ｉ（ｔ_２）とを乗じた値となる。すなわち、この場合、第１実施形態に係る時間発展部３６は、第２時間定数ｄｔ´を、第１時間定数ｄｔに置き換えて処理を実行する。

以上のような第３実施形態に係る計算装置１０は、Ｋ個の不等式による制約を含む２次計画問題をシミュレーションにより解くことができる。

（第４実施形態）
第４実施形態に係る計算装置１０について説明する。第４実施形態に係る計算装置１０は、第３実施形態に係る計算装置１０と略同一の構成を有するので、略同一の構成を有するブロックには、同一の符号を付けて、相違点を除き詳細な説明を省略する。

第４実施形態に係る計算装置１０は、式（８）に示すハミルトニアンＨを運動量および位置のそれぞれにより偏微分することにより得られる連立常微分方程式を解く。

式（８）に示すハミルトニアンＨを偏微分することにより得られる運動方程式は、式（２９）および式（３０）により表される。

式（３０）において、βおよびψは、予め定められた定数である。

第４実施形態に係る計算装置１０は、式（２９）および式（３０）により表される連立常微分方程式を離散解法により解く。より具体的には、計算装置１０は、ｔを開始時刻ｔ_０からサンプリング期間ずつ増加させて、ｔ毎に、式（２９）および式（３０）を用いてｘ_ｉおよびｙ_ｉを更新する。そして、計算装置１０は、ｔを十分に大きくした終了時刻ｔ_Ｔにおけるｘ_ｉの値を、解として出力する。これにより、計算装置１０は、Ｎ個のｘ_ｉおよびＮ個のｙ_ｉを並行して更新しながら、解を算出することができる。

図１０は、第４実施形態に係るｍ段目の発展部７２－ｍの構成を示す図である。Ｍ段の発展部７２－１～７２－Ｍのそれぞれは、第３制約部７６と、第４制約部７８と、第５制約部８０とをさらに含む。

ｍ段目の発展部７２－ｍに含まれる第３制約部７６、第４制約部７８および第５制約部８０は、次のような動作をする。

第３制約部７６は、緩和部７４から、第２時刻ｔ_２におけるＮ個の第２変数ｙ_ｉ（ｔ_２）を取得する。第３制約部７６は、第２加算部６４から、第ｍ経過時刻（ｔ_１＋ｍ×ｄｔ´）におけるＮ個の第１変数ｘ_ｉ（ｔ_１＋ｍ×ｄｔ´）を取得する。

第３制約部７６は、第ｍ経過時刻（ｔ_１＋ｍ×ｄｔ´）におけるＮ個の第１変数ｘ_ｉ（ｔ_１＋ｍ×ｄｔ´）の合計値が予め定められた規格値となる方向に、取得した第２時刻ｔ_２におけるＮ個の第２変数ｙ_ｉ（ｔ_２）のそれぞれを変化させる。本例において、予め定められた規格値は、１である。なお、予め定められた規格値は、α（ｔ）により表される値であってもよい。

第４制約部７８は、第２加算部６４から、第ｍ経過時刻（ｔ_１＋ｍ×ｄｔ´）におけるＮ個の第１変数ｘ_ｉ（ｔ_１＋ｍ×ｄｔ´）を取得する。第４制約部７８は、第ｍ経過時刻（ｔ_１＋ｍ×ｄｔ´）におけるＮ個の第１変数ｘ_ｉ（ｔ_１＋ｍ×ｄｔ´）のそれぞれが、予め定められた範囲から外れているか否かを判断する。例えば、第４制約部７８は、第ｍ経過時刻（ｔ_１＋ｍ×ｄｔ´）におけるＮ個の第１変数ｘ_ｉ（ｔ_１＋ｍ×ｄｔ´）のそれぞれが、０より小さいか否かを判断する。

第４制約部７８は、第ｍ経過時刻（ｔ_１＋ｍ×ｄｔ´）におけるＮ個の第１変数ｘ_ｉ（ｔ_１＋ｍ×ｄｔ´）のうち予め定められた範囲から外れている１または複数の第１変数ｘ_ｉ（ｔ_１＋ｍ×ｄｔ´）のそれぞれを、予め定められた範囲の値とする。例えば、第４制約部７８は、第ｍ経過時刻（ｔ_１＋ｍ×ｄｔ´）におけるＮ個の第１変数ｘ_ｉ（ｔ_１＋ｍ×ｄｔ´）のうち０より小さい１または複数の第１変数ｘ_ｉ（ｔ_１＋ｍ×ｄｔ´）のそれぞれを、０とする。なお、第４制約部７８は、第ｍ経過時刻（ｔ_１＋ｍ×ｄｔ´）におけるＮ個の第１変数ｘ_ｉ（ｔ_１＋ｍ×ｄｔ´）のうち予め定められた範囲内の１または複数の第１変数ｘ_ｉ（ｔ_１＋ｍ×ｄｔ´）のそれぞれについては、そのまま出力する。

第４制約部７８は、第ｍ経過時刻（ｔ_１＋ｍ×ｄｔ´）におけるＮ個の第１変数ｘ_ｉ（ｔ_１＋ｍ×ｄｔ´）を、第１制約部６６に与える。そして、第１制約部６６は、第４制約部７８から取得した第ｍ経過時刻（ｔ_１＋ｍ×ｄｔ´）におけるＮ個の第１変数ｘ_ｉ（ｔ_１＋ｍ×ｄｔ´）に対して制約処理を実行する。

第５制約部８０は、第３制約部７６から、第２時刻ｔ_２におけるＮ個の第２変数ｙ_ｉ（ｔ_２）を取得する。第５制約部８０は、第２加算部６４から、第ｍ経過時刻（ｔ_１＋ｍ×ｄｔ´）におけるＮ個の第１変数ｘ_ｉ（ｔ_１＋ｍ×ｄｔ´）を取得する。

第５制約部８０は、第ｍ経過時刻（ｔ_１＋ｍ×ｄｔ´）におけるＮ個の第１変数ｘ_ｉ（ｔ_１＋ｍ×ｄｔ´）のそれぞれが、予め定められた範囲から外れているか否かを判断する。なお、第５制約部８０は、第４制約部７８による判断結果に基づき、第ｍ経過時刻（ｔ_１＋ｍ×ｄｔ´）におけるＮ個の第１変数ｘ_ｉ（ｔ_１＋ｍ×ｄｔ´）のそれぞれが、予め定められた範囲から外れているか否かを判断してもよい。例えば、第５制約部８０は、第ｍ経過時刻（ｔ_１＋ｍ×ｄｔ´）におけるＮ個の第１変数ｘ_ｉ（ｔ_１＋ｍ×ｄｔ´）のそれぞれが、０より小さいか否かを判断する。

第５制約部８０は、取得した第２時刻ｔ_２におけるＮ個の第２変数ｙ_ｉ（ｔ_２）のうちの、予め定められた範囲から外れている１または複数の第１変数ｘ_ｉ（ｔ_１＋ｍ×ｄｔ´）に対応する１または複数の第２変数ｙ_ｉ（ｔ_２）のそれぞれを、対応する第１変数ｘ_ｉ（ｔ_１＋ｍ×ｄｔ´）が予め定められた範囲内となる方向に変化させる。

例えば、第５制約部８０は、予め定められた範囲から外れている１または複数の第１変数ｘ_ｉ（ｔ_１＋ｍ×ｄｔ´）に対応する１または複数の第２変数ｙ_ｉ（ｔ_２）のそれぞれを、対応する第１変数ｘ_ｉ（ｔ_１＋ｍ×ｄｔ´）が０となる方向に変化させる。なお、第５制約部８０は、第ｍ経過時刻（ｔ_１＋ｍ×ｄｔ´）におけるＮ個の第１変数ｘ_ｉ（ｔ_１＋ｍ×ｄｔ´）のうち予め定められた範囲内の１または複数の第１変数ｘ_ｉ（ｔ_１＋ｍ×ｄｔ´）に対応する１または複数の第２変数ｙ_ｉ（ｔ_２）については、そのまま出力する。

第５制約部８０は、第２時刻ｔ_２におけるＮ個の第２変数ｙ_ｉ（ｔ_２）を、第２制約部６８に与える。第２制約部６８は、第５制約部８０から取得した第２時刻ｔ_２におけるＮ個の第２変数ｙ_ｉ（ｔ_２）に対して制約処理を実行する。

図１１は、第４実施形態に係る演算部２０の処理の流れを示すフローチャートである。第４実施形態に係る演算部２０は、図１１に示す流れで処理を実行する。図１１のフローチャートは、図９のフローチャートと比較して、Ｓ５１とＳ１７との間に、Ｓ６１～Ｓ６４の処理が追加されている。従って、第４実施形態に係る時間発展部３６は、Ｍ回のループの処理（Ｓ４１とＳ４２との間の処理）内において、緩和処理（Ｓ５１）の後に、Ｓ６１～Ｓ６４の処理を実行する。

まず、Ｓ６１において、時間発展部３６は、第３制約処理を実行する。時間発展部３６は、第３制約処理において、Ｎ個の第１変数ｘ_ｉの合計値が予め定められた規格値となる方向に、Ｎ個の第２変数ｙ_ｉ（ｔ_２）のそれぞれを変化させる。本例において、予め定められた規格値は、１である。例えば、Ｓ６１において、時間発展部３６は、ｉ＝１～Ｎのそれぞれについて、式（３１）の演算を実行する。

なお、式（３１）において、βは、予め定められた定数である。

続いて、Ｓ６２において、時間発展部３６は、Ｎ個の第１変数ｘ_ｉのそれぞれが、０より小さいか否かを判断する。時間発展部３６は、Ｎ個の第１変数ｘ_ｉのうち、０より小さい第１変数ｘ_ｉについては（Ｓ６２のＹｅｓ）、Ｓ６３およびＳ６４の処理を実行する。時間発展部３６は、Ｎ個の第１変数ｘ_ｉのうち、０以上である第１変数ｘ_ｉについては（Ｓ６２のＮｏ）、Ｓ６３およびＳ６４の処理を実行せずに、処理をＳ１７に進める。

Ｓ６３において、時間発展部３６は、第５制約処理を実行する。時間発展部３６は、第５制約処理において、Ｎ個の第１変数ｘ_ｉのうち０より小さい第１変数ｘ_ｉに対応する第２変数ｙ_ｉに対して、式（３２）の演算を実行する。

なお、式（３２）において、ψは、予め定められた定数である。これにより、時間発展部３６は、予め定められた範囲から外れている第１変数ｘ_ｉに対応する第２変数ｙ_ｉを、対応する第１変数ｘ_ｉが予め定められた範囲内となる方向に変化させることができる。

Ｓ６４において、時間発展部３６は、第４制約処理を実行する。時間発展部３６は、第４制約処理において、Ｎ個の第１変数ｘ_ｉのうち０より小さい第１変数ｘ_ｉを０とする。これにより、時間発展部３６は、予め定められた範囲から外れている第１変数ｘ_ｉを予め定められた範囲の値とすることができる。Ｓ６４の処理を終えると、時間発展部３６は、処理をＳ１７に進める。

なお、第１実施形態に係る時間発展部３６が、第３制約処理、第４制約処理および第５制約処理を実行する構成であってもよい。この場合、式（３１）および式（３２）の演算における第２時間定数ｄｔ´は、第１時間定数ｄｔに置き換えられる。

図１２は、第４実施形態に係る後処理部４２の処理の流れを示すフローチャートである。第４実施形態に係る後処理部４２は、図１２に示す流れで処理を実行する。図１２のフローチャートは、図５のフローチャートと比較して、Ｓ３３とＳ３４との間に、Ｓ７１からＳ７３の処理が追加されている。従って、第４実施形態に係る後処理部４２は、Ｓ３３の後に、Ｓ７１の処理を実行する。

Ｓ７１において、後処理部４２は、Ｎ個の第１変数ｘ_ｉの合計値を予め定められた規格値とするように、Ｎ個の第１変数ｘ_ｉのそれぞれの値を修正する。本例において、予め定められた規格値は、１である。例えば、Ｓ７１において、後処理部４２は、ｉ＝１～Ｎのそれぞれについて、式（３３）の演算を実行する。

続いて、Ｓ７２において、後処理部４２は、Ｎ個の第１変数ｘ_ｉのそれぞれが、０より小さいか否かを判断する。後処理部４２は、Ｎ個の第１変数ｘ_ｉのうち、０より小さい第１変数ｘ_ｉについては（Ｓ７２のＹｅｓ）、Ｓ７３の処理を実行する。後処理部４２は、Ｎ個の第１変数ｘ_ｉのうち、０以上である第１変数ｘ_ｉについては（Ｓ７２のＮｏ）、Ｓ７３の処理を実行せずに、処理をＳ３４に進める。

Ｓ７３において、後処理部４２は、Ｎ個の第１変数ｘ_ｉのうち０より小さい第１変数ｘ_ｉを０とする。これにより、後処理部４２は、予め定められた範囲から外れている第１変数ｘ_ｉを予め定められた範囲の値とすることができる。後処理部４２は、Ｓ７３の処理を終えると、処理をＳ３４に進める。

以上のような第４実施形態に係る計算装置１０は、Ｋ個の不等式による制約、Ｎ個の変数のそれぞれを予め定められた範囲内とする制約、および、Ｎ個の変数の合計を予め定められた規格値とする制約を含む２次計画問題をシミュレーションにより解くことができる。

図１３は、計算装置１０のハードウェアブロック図である。計算装置１０は、一例として、図１３に示すような一般のコンピュータ（情報処理装置）と同様のハードウェア構成により実現される。計算装置１０は、図１３に示すような１つのコンピュータにより実現されてもよいし、協働して動作する複数のコンピュータにより実現されてもよい。また、計算装置１０は、一部に専用のハードウェア回路を備える構成であってもよい。

計算装置１０は、メモリ２０４と、１または複数のハードウェアプロセッサ２０６と、記憶装置２０８と、操作装置２１０と、表示装置２１２と、通信装置２１４とを備える。各部は、バスにより接続される。

メモリ２０４は、例えば、ＲＯＭ２２２と、ＲＡＭ２２４とを含む。ＲＯＭ２２２は、計算装置１０の制御に用いられるプログラムおよび各種設定情報等を書き換え不可能に記憶する。ＲＡＭ２２４は、ＳＤＲＡＭ（Synchronous Dynamic Random Access Memory）等の揮発性の記憶媒体である。ＲＡＭ２２４は、１または複数のハードウェアプロセッサ２０６の作業領域として機能する。

１または複数のハードウェアプロセッサ２０６は、メモリ２０４（ＲＯＭ２２２およびＲＡＭ２２４）にバスを介して接続される。１または複数のハードウェアプロセッサ２０６のそれぞれは、例えば、ＣＰＵ（Central Processing Unit）であってもよいし、演算用のハードウェア回路であってもよい。

１または複数のハードウェアプロセッサ２０６は、ＲＡＭ２２４の所定領域を作業領域としてＲＯＭ２２２または記憶装置２０８に予め記憶された各種プログラムとの協働により各種処理を実行し、計算装置１０を構成する各部の動作を統括的に制御する。また、１または複数のハードウェアプロセッサ２０６は、ＲＯＭ２２２または記憶装置２０８に予め記憶されたプログラムとの協働により、操作装置２１０、表示装置２１２および通信装置２１４等を制御する。

記憶装置２０８は、フラッシュメモリ等の半導体による記憶媒体、磁気的または光学的に記録可能な記憶媒体等の書き換え可能な記録装置である。記憶装置２０８は、計算装置１０の制御に用いられるプログラムおよび各種設定情報等を記憶する。

操作装置２１０は、マウスおよびキーボード等の入力デバイスである。操作装置２１０は、ユーザから操作入力された情報を受け付け、受け付けた情報を１または複数のハードウェアプロセッサ２０６に出力する。

表示装置２１２は、情報をユーザに表示する。表示装置２１２は、１または複数のハードウェアプロセッサ２０６から情報等を受け取り、受け取った情報を表示する。なお、通信装置２１４または記憶装置２０８等に情報を出力する場合、計算装置１０は、表示装置２１２を備えなくてもよい。通信装置２１４は、外部の機器と通信して、ネットワーク等を介して情報を送受信する。

本実施形態の計算装置１０で実行されるプログラムは、インストール可能な形式または実行可能な形式のファイルでＣＤ－ＲＯＭ、フレキシブルディスク（ＦＤ）、ＣＤ－Ｒ、ＤＶＤ（Digital Versatile Disk）等のコンピュータで読み取り可能な記録媒体に記録されて提供される。

また、本実施形態の計算装置１０で実行されるプログラムを、インターネット等のネットワークに接続されたコンピュータ上に格納し、ネットワーク経由でダウンロードさせることにより提供するように構成してもよい。また、本実施形態の計算装置１０で実行されるプログラムをインターネット等のネットワーク経由で提供または配布するように構成してもよい。また、本実施形態の計算装置１０で実行されるプログラムを、ＲＯＭ等に予め組み込んで提供するように構成してもよい。

情報処理装置を計算装置１０として機能させるためのプログラムは、例えば、相互作用モジュール（行列演算モジュールおよび第１乗算モジュール）、第１加算モジュール、時間発展モジュール（第２乗算モジュール、第２加算モジュール、第１制約モジュール、第２制約モジュール、緩和モジュール、第３制約モジュール、第４制約モジュールおよび第５制約モジュール）、および、後処理モジュールを含むモジュール構成となっている。このプログラムは、１または複数のハードウェアプロセッサ２０６により実行されることにより各モジュールがメモリ２０４のＲＡＭ２２４にロードされ、１または複数のハードウェアプロセッサ２０６を、相互作用部３２（行列演算部５４および第１乗算部５６）、第１加算部３４、時間発展部３６（第２乗算部６２、第２加算部６４、第１制約部６６、第２制約部６８、緩和部７４、第３制約部７６、第４制約部７８および第５制約部８０）、および、後処理部４２として機能させる。なお、これらの構成は、一部または全部がハードウェアにより構成されていてもよい。

本発明のいくつかの実施形態を説明したが、これらの実施形態は、例として提示したものであり、発明の範囲を限定することは意図していない。これら新規な実施形態は、その他の様々な形態で実施されることが可能であり、発明の要旨を逸脱しない範囲で、種々の省略、置き換え、変更を行うことができる。これら実施形態やその変形は、発明の範囲や要旨に含まれるとともに、請求の範囲に記載された発明とその均等の範囲に含まれる。

（技術案）
なお、上記の実施形態を、以下の技術案にまとめることができる.

技術案１
第１時刻におけるＮ個（Ｎは２以上の整数）の第１変数を相互作用させたＮ個の第１中間変数と、前記第１時刻におけるＮ個の第２変数とに基づき、前記第１時刻からサンプリング期間を経過した第２時刻における前記Ｎ個の第２変数を算出する変数算出部と、
時間発展処理を実行して、前記第２時刻における前記Ｎ個の第１変数を生成する時間発展部と、
を備え、
前記時間発展部は、
前記第２時刻における前記Ｎ個の第２変数と、前記第１時刻における前記Ｎ個の第１変数とに基づき、前記第２時刻における前記Ｎ個の第１変数を算出し、
前記第２時刻における前記Ｎ個の第１変数が予め定められた制約条件を満たしていない場合、前記第２時刻における前記Ｎ個の第２変数を前記制約条件を満たす方向に変化させる
計算装置。

技術案２
前記変数算出部は、
前記第１時刻における前記Ｎ個の第１変数と、予め設定されたＮ行×Ｎ列の係数を含む係数行列とを行列演算することにより算出されたＮ個の値のそれぞれに、前記サンプリング期間を表す第１時間定数の正負を反転させた値を乗じたＮ個の第１中間変数を生成する相互作用部と、
前記第１時刻における前記Ｎ個の第２変数のそれぞれに、前記Ｎ個の第１中間変数のうちの対応する第１中間変数を加算することにより、前記第２時刻における前記Ｎ個の第２変数を算出する第１加算部と、
を有し、
前記時間発展部は、
前記第２時刻における前記Ｎ個の第２変数のそれぞれに前記第１時間定数を乗じたＮ個の第２中間変数を算出し、
前記第１時刻における前記Ｎ個の第１変数のそれぞれに、前記Ｎ個の第２中間変数のうちの対応する第２中間変数を加算することにより、前記第２時刻における前記Ｎ個の第１変数を算出する
技術案１に記載の計算装置。

技術案３
開始時刻から終了時刻まで前記サンプリング期間の間隔で順次に時刻を増加させ、それぞれの時刻に対する前記Ｎ個の第１変数および前記Ｎ個の第２変数を算出させるように処理を管理する管理部と、
前記終了時刻における前記Ｎ個の第１変数を出力する出力部と、
をさらに備える技術案２に記載の計算装置。

技術案４
前記時間発展部は、前記時間発展処理において、前記第２時刻における前記Ｎ個の第１変数が前記制約条件を満たしていない場合、前記第２時刻における前記Ｎ個の第１変数を前記制約条件を満たす方向に変化させる
技術案３に記載の計算装置。

技術案５
前記終了時刻における前記Ｎ個の第１変数が前記制約条件を満たしていない場合、前記終了時刻における前記Ｎ個の第１変数が前記制約条件を満たすまで、前記終了時刻における前記Ｎ個の第１変数を前記制約条件を満たす方向に変化させる修正処理を繰り返す後処理部をさらに備え、
前記出力部は、前記制約条件を満たした後の前記終了時刻における前記Ｎ個の第１変数を出力する
技術案４に記載の計算装置。

技術案６
前記時間発展部は、前記第１時刻における前記Ｎ個の第２変数に対してＭ回（Ｍは１以上の整数）の時間発展処理を実行して、前記第２時刻における前記Ｎ個の第１変数を生成し、
前記Ｍ回の時間発展処理のうちのｍ回目（ｍは、１からＭまでの整数）の時間発展処理において、前記時間発展部は、
前記第２時刻における前記Ｎ個の第２変数のそれぞれに、前記サンプリング期間の１／Ｍの微小期間を表す第２時間定数を乗じた前記Ｎ個の第２中間変数を算出し、
前記第１時刻から前記微小期間の（ｍ－１）倍の時間を経過した時刻における前記Ｎ個の第１変数のそれぞれに、前記Ｎ個の第２中間変数のうちの対応する第２中間変数を加算することにより、前記第１時刻から前記微小期間のｍ倍の時間を経過した第ｍ経過時刻における前記Ｎ個の第１変数を算出し、
前記第ｍ経過時刻における前記Ｎ個の第１変数が前記制約条件を満たしていない場合、前記第ｍ経過時刻における前記Ｎ個の第１変数、および、前記Ｎ個の第２変数を前記制約条件を満たす方向に変化させる
技術案５に記載の計算装置。

技術案７
前記制約条件は、不等式で表される
技術案１から６の何れか１項に記載の計算装置。

技術案８
前記時間発展部は、
前記不等式を満たしていない場合、前記第２時刻における前記Ｎ個の第１変数のうちの前記不等式に含まれる複数の第１変数により表される位置を、前記不等式により表される境界に対して垂直に近づける方向に、前記複数の第１変数、および、前記不等式に含まれる前記複数の第１変数に対応する複数の第２変数を変化させる
技術案７に記載の計算装置。

技術案９
前記制約条件は、Ｋ個（Ｋは１以上の整数）の不等式を含む連立一次不等式で表され、
前記Ｋ個の不等式のそれぞれは、前記Ｎ個の第１変数のうちの１以上Ｎ以下の第１変数を変数として含む
技術案６に記載の計算装置。

技術案１０
前記ｍ回目の時間発展処理において、前記時間発展部は、
前記Ｋ個の不等式のそれぞれ毎に、
不等式を満たしていない場合、前記第２時刻における前記Ｎ個の第１変数のうちの前記不等式に含まれる複数の第１変数を前記不等式を満たす方向に変化させ、および、前記不等式に含まれる前記複数の第１変数に対応する複数の第２変数を、前記不等式を満たす方向に変化させる
技術案９に記載の計算装置。

技術案１１
前記ｍ回目の時間発展処理において、前記時間発展部は、前記第２時刻における前記Ｎ個の第２変数のそれぞれに、緩和項を加算し、
前記緩和項は、前記第２時間定数と、予め定められた緩和定数の正負を反転させた値と、前記第２時刻における前記Ｎ個の第２変数のうちの対応する第２変数とを乗じた値である
技術案１０に記載の計算装置。

技術案１２
前記ｍ回目の時間発展処理において、前記時間発展部は、
前記第ｍ経過時刻における前記Ｎ個の第１変数の合計値が予め定められた規格値となる方向に、前記第２時刻における前記Ｎ個の第２変数のそれぞれを変化させる
技術案１０または１１に記載の計算装置。

技術案１３
前記ｍ回目の時間発展処理において、前記時間発展部は、
前記第ｍ経過時刻における前記Ｎ個の第１変数のうちの予め定められた範囲から外れている１または複数の第１変数のそれぞれを、前記予め定められた範囲の値とし、
前記第２時刻における前記Ｎ個の第２変数のうちの、前記１または複数の第１変数に対応する１または複数の第２変数のそれぞれを、対応する第１変数が前記予め定められた範囲となる方向に変化させる
技術案１０から１２の何れか１項に記載の計算装置。

技術案１４
前記相互作用部は、下記の式（１０１）の演算をｉ＝１からＮまで実行し、

前記第１加算部は、下記の式（１０２）の演算をｉ＝１からＮまで実行し、

ｉ，ｊは、１以上Ｎ以下の整数であり、
ａ_ｉは、前記Ｎ個の第１中間変数のうちのｉ番目の第１中間変数であり、
ｄｔは、前記第１時間定数であり、
Ａ_ｉｊは、前記係数行列におけるｉ行ｊ列の値であり、
ｘ_ｊは、前記Ｎ個の第１変数のうちのｊ番目の第１変数であり、
ｙ_ｉは、前記Ｎ個の第２変数のうちのｉ番目の第２変数である
技術案１０に記載の計算装置。

技術案１５
前記時間発展部は、前記Ｋ個の不等式のそれぞれ毎に、ｒ_ｋが０より小さいか否かを判断し、
ｒ_ｋは、前記Ｋ個の不等式のうちのｋ番目（ｋは、１以上Ｋ以下の整数）の不等式の判断において、下記の式（１０３）により算出される値であり、

Ｂ_ｊｋは、前記ｋ番目の不等式における、前記ｊ番目の第１変数に対する係数である
技術案１４に記載の計算装置。

技術案１６
前記時間発展部は、前記ｋ番目の不等式の判断において、ｒ_ｋが０より小さい場合、下記の式（１０４）および式（１０５）の演算をｉ＝１からＮまで実行し、

ｘ_ｉは、前記Ｎ個の第１変数のうちのｉ番目の第１変数であり、
ｄｔ´は、前記第２時間定数であり、
φは、予め設定された定数であり、
Ｂ_ｉｋは、前記ｋ番目の不等式における、前記ｉ番目の第１変数に対する係数である
技術案１５に記載の計算装置。

技術案１７
前記ｍ回目の時間発展処理において、前記時間発展部は、下記の式（１０６）の演算をｉ＝１からＮまで実行する

μは、予め定められた緩和定数である
技術案１６に記載の計算装置。

技術案１８
前記ｍ回目の時間発展処理において、前記時間発展部は、下記の式（１０７）の演算をｉ＝１からＮまで実行する

βは、予め定められた定数である
技術案１６または１７に記載の計算装置。

技術案１９
前記ｍ回目の時間発展処理において、前記時間発展部は、
前記Ｎ個の第２変数のうち、前記０より小さい第１変数に対応する第２変数に対して、式（１０８）の演算を実行し、

前記Ｎ個の第１変数のうち、０より小さい第１変数を０とし、
ψは、予め定められた定数である
技術案１６から１８の何れか１項に記載の計算装置。

技術案２０
前記修正処理において、前記後処理部は、
下記の式（１０９）の演算をｉ＝１からＮまで実行し、

前記Ｎ個の第１変数のうち、０より小さい第１変数を０とし、
前記Ｋ個の不等式のそれぞれ毎に、
ｒ_ｋが０より小さいか否かを判断し、
ｒ_ｋが０より小さい場合、下記の式（１１０）の演算をｉ＝１からＮまで実行する

γは、予め定められた定数である
技術案１６から１９の何れか１項に記載の計算装置。

技術案２１
情報処理装置が計算をするための計算方法であって、
前記情報処理装置が、第１時刻におけるＮ個（Ｎは２以上の整数）の第１変数を相互作用させたＮ個の第１中間変数と、前記第１時刻におけるＮ個の第２変数とに基づき、前記第１時刻からサンプリング期間を経過した第２時刻における前記Ｎ個の第２変数を算出し
前記情報処理装置が、時間発展処理を実行して、前記第２時刻における前記Ｎ個の第１変数を生成し、
前記情報処理装置が、前記時間発展処理において、
前記第２時刻における前記Ｎ個の第２変数と、前記第１時刻における前記Ｎ個の第１変数とに基づき、前記第２時刻における前記Ｎ個の第１変数を算出し、
前記第２時刻における前記Ｎ個の第１変数が予め定められた制約条件を満たしていない場合、前記第２時刻における前記Ｎ個の第２変数を前記制約条件を満たす方向に変化させる
計算方法。

技術案２２
情報処理装置に、計算装置として機能させるためのプログラムであって、
前記情報処理装置を、
第１時刻におけるＮ個（Ｎは２以上の整数）の第１変数を相互作用させたＮ個の第１中間変数と、前記第１時刻におけるＮ個の第２変数とに基づき、前記第１時刻からサンプリング期間を経過した第２時刻における前記Ｎ個の第２変数を算出する変数算出部と、
時間発展処理を実行して、前記第２時刻における前記Ｎ個の第１変数を生成する時間発展部と、
して機能させ、
前記時間発展部は、
前記第２時刻における前記Ｎ個の第２変数と、前記第１時刻における前記Ｎ個の第１変数とに基づき、前記第２時刻における前記Ｎ個の第１変数を算出し、
前記第２時刻における前記Ｎ個の第１変数が予め定められた制約条件を満たしていない場合、前記第２時刻における前記Ｎ個の第２変数を前記制約条件を満たす方向に変化させる
プログラム。

技術案２３
ハミルトニアンＨを設定可能な計算装置であって、
前記ハミルトニアンＨを終了時刻まで時間発展させ、
前記終了時刻におけるｘ_ｉの値を取得して出力し、
前記ハミルトニアンＨは、式（１１１）により表され、

ｘ_ｉ，ｘ_ｊは、位置を表し、
ｙ_ｉは、運動量を表し、
Ｎは、２以上の整数であり、
ｉおよびｊは、１以上Ｎ以下の任意の整数であり、
Ａ_ｉｊは、ｉ番目のｘ_ｉと、ｊ番目のｘ_ｊとの間の相互作用を表し、
Ｂ_ｉｋは、Ｎ行Ｋ列の制約行列Ｂに含まれるｉ行ｋ列の係数であり、
ｔは、時刻を表し、
α（ｔ）は、ｔを変数とした非減少関数であり、少なくとも前記終了時刻において１となり、
Φ（ξ）は、式（１１２）により表される関数であり、

φは、予め定められた定数であり、
θ（ξ）は、式（１１３）により表される関数である

計算装置。

１０ 計算装置
２０ 演算部
２２ 入力部
２４ 出力部
２６ 設定部
３２ 相互作用部
３４ 第１加算部
３６ 時間発展部
３８ 第１変数メモリ
４０ 第２変数メモリ
４２ 後処理部
４４ 管理部
５２ 係数行列メモリ
５４ 行列演算部
５６ 第１乗算部
６２ 第２乗算部
６４ 第２加算部
６６ 第１制約部
６８ 第２制約部
７４ 緩和部
７６ 第３制約部
７８ 第４制約部
８０ 第５制約部

Claims (13)

1. 第１時刻におけるＮ個（Ｎは２以上の整数）の第１変数を相互作用させたＮ個の第１中間変数と、前記第１時刻におけるＮ個の第２変数とに基づき、前記第１時刻からサンプリング期間を経過した第２時刻における前記Ｎ個の第２変数を算出する変数算出部と、
   時間発展処理を実行して、前記第２時刻における前記Ｎ個の第１変数を生成する時間発展部と、
   を備え、
   前記時間発展部は、
   前記第２時刻における前記Ｎ個の第２変数と、前記第１時刻における前記Ｎ個の第１変数とに基づき、前記第２時刻における前記Ｎ個の第１変数を算出し、
   前記第２時刻における前記Ｎ個の第１変数が予め定められた制約条件を満たしていない場合、前記第２時刻における前記Ｎ個の第２変数を前記制約条件を満たす方向に変化させる
   計算装置。
2. 前記変数算出部は、
   前記第１時刻における前記Ｎ個の第１変数と、予め設定されたＮ行×Ｎ列の係数を含む係数行列とを行列演算することにより算出されたＮ個の値のそれぞれに、前記サンプリング期間を表す第１時間定数の正負を反転させた値を乗じたＮ個の第１中間変数を生成する相互作用部と、
   前記第１時刻における前記Ｎ個の第２変数のそれぞれに、前記Ｎ個の第１中間変数のうちの対応する第１中間変数を加算することにより、前記第２時刻における前記Ｎ個の第２変数を算出する第１加算部と、
   を有し、
   前記時間発展部は、
   前記第２時刻における前記Ｎ個の第２変数のそれぞれに前記第１時間定数を乗じたＮ個の第２中間変数を算出し、
   前記第１時刻における前記Ｎ個の第１変数のそれぞれに、前記Ｎ個の第２中間変数のうちの対応する第２中間変数を加算することにより、前記第２時刻における前記Ｎ個の第１変数を算出する
   請求項１に記載の計算装置。
3. 開始時刻から終了時刻まで前記サンプリング期間の間隔で順次に時刻を増加させ、それぞれの時刻に対する前記Ｎ個の第１変数および前記Ｎ個の第２変数を算出させるように処理を管理する管理部と、
   前記終了時刻における前記Ｎ個の第１変数を出力する出力部と、
   をさらに備える請求項２に記載の計算装置。
4. 前記時間発展部は、前記時間発展処理において、前記第２時刻における前記Ｎ個の第１変数が前記制約条件を満たしていない場合、前記第２時刻における前記Ｎ個の第１変数を前記制約条件を満たす方向に変化させる
   請求項３に記載の計算装置。
5. 前記終了時刻における前記Ｎ個の第１変数が前記制約条件を満たしていない場合、前記終了時刻における前記Ｎ個の第１変数が前記制約条件を満たすまで、前記終了時刻における前記Ｎ個の第１変数を前記制約条件を満たす方向に変化させる修正処理を繰り返す後処理部をさらに備え、
   前記出力部は、前記制約条件を満たした後の前記終了時刻における前記Ｎ個の第１変数を出力する
   請求項４に記載の計算装置。
6. 前記時間発展部は、前記第１時刻における前記Ｎ個の第２変数に対してＭ回（Ｍは１以上の整数）の時間発展処理を実行して、前記第２時刻における前記Ｎ個の第１変数を生成し、
   前記Ｍ回の時間発展処理のうちのｍ回目（ｍは、１からＭまでの整数）の時間発展処理において、前記時間発展部は、
   前記第２時刻における前記Ｎ個の第２変数のそれぞれに、前記サンプリング期間の１／Ｍの微小期間を表す第２時間定数を乗じた前記Ｎ個の第２中間変数を算出し、
   前記第１時刻から前記微小期間の（ｍ－１）倍の時間を経過した時刻における前記Ｎ個の第１変数のそれぞれに、前記Ｎ個の第２中間変数のうちの対応する第２中間変数を加算することにより、前記第１時刻から前記微小期間のｍ倍の時間を経過した第ｍ経過時刻における前記Ｎ個の第１変数を算出し、
   前記第ｍ経過時刻における前記Ｎ個の第１変数が前記制約条件を満たしていない場合、前記第ｍ経過時刻における前記Ｎ個の第１変数、および、前記Ｎ個の第２変数を前記制約条件を満たす方向に変化させる
   請求項５に記載の計算装置。
7. 前記制約条件は、不等式で表される
   請求項１から６の何れか１項に記載の計算装置。
8. 前記時間発展部は、
   前記不等式を満たしていない場合、前記第２時刻における前記Ｎ個の第１変数のうちの前記不等式に含まれる複数の第１変数により表される位置を、前記不等式により表される境界に対して垂直に近づける方向に、前記複数の第１変数、および、前記不等式に含まれる前記複数の第１変数に対応する複数の第２変数を変化させる
   請求項７に記載の計算装置。
9. 前記制約条件は、Ｋ個（Ｋは１以上の整数）の不等式を含む連立一次不等式で表され、
   前記Ｋ個の不等式のそれぞれは、前記Ｎ個の第１変数のうちの１以上Ｎ以下の第１変数を変数として含む
   請求項６に記載の計算装置。
10. 前記時間発展部は、
    前記Ｎ個の第１変数の合計値が予め定められた規格値となる方向に、前記第２時刻における前記Ｎ個の第２変数のそれぞれを変化させる
    請求項１から９の何れか１項に記載の計算装置。
11. 前記時間発展部は、
    前記Ｎ個の第１変数のうちの予め定められた範囲から外れている１または複数の第１変数のそれぞれを、前記予め定められた範囲の値とし、
    前記第２時刻における前記Ｎ個の第２変数のうちの、前記１または複数の第１変数に対応する１または複数の第２変数のそれぞれを、対応する第１変数が前記予め定められた範囲となる方向に変化させる
    請求項１から１０の何れか１項に記載の計算装置。
12. 情報処理装置が計算をするための計算方法であって、
    前記情報処理装置が、第１時刻におけるＮ個（Ｎは２以上の整数）の第１変数を相互作用させたＮ個の第１中間変数と、前記第１時刻におけるＮ個の第２変数とに基づき、前記第１時刻からサンプリング期間を経過した第２時刻における前記Ｎ個の第２変数を算出し
    前記情報処理装置が、時間発展処理を実行して、前記第２時刻における前記Ｎ個の第１変数を生成し、
    前記情報処理装置が、前記時間発展処理において、
    前記第２時刻における前記Ｎ個の第２変数と、前記第１時刻における前記Ｎ個の第１変数とに基づき、前記第２時刻における前記Ｎ個の第１変数を算出し、
    前記第２時刻における前記Ｎ個の第１変数が予め定められた制約条件を満たしていない場合、前記第２時刻における前記Ｎ個の第２変数を前記制約条件を満たす方向に変化させる
    計算方法。
13. 情報処理装置に、計算装置として機能させるためのプログラムであって、
    前記情報処理装置を、
    第１時刻におけるＮ個（Ｎは２以上の整数）の第１変数を相互作用させたＮ個の第１中間変数と、前記第１時刻におけるＮ個の第２変数とに基づき、前記第１時刻からサンプリング期間を経過した第２時刻における前記Ｎ個の第２変数を算出する変数算出部と、
    時間発展処理を実行して、前記第２時刻における前記Ｎ個の第１変数を生成する時間発展部と、
    して機能させ、
    前記時間発展部は、
    前記第２時刻における前記Ｎ個の第２変数と、前記第１時刻における前記Ｎ個の第１変数とに基づき、前記第２時刻における前記Ｎ個の第１変数を算出し、
    前記第２時刻における前記Ｎ個の第１変数が予め定められた制約条件を満たしていない場合、前記第２時刻における前記Ｎ個の第２変数を前記制約条件を満たす方向に変化させる
    プログラム。

Images