知財求人 - 知財ポータルサイト「IP Force」
(19)【発行国】日本国特許庁(JP)
(12)【公報種別】特許公報(B2)
(11)【特許番号】
(24)【登録日】2023-08-01
(45)【発行日】2023-08-09
(54)【発明の名称】最適化装置および最適化方法
(51)【国際特許分類】
G06N 99/00 20190101AFI20230802BHJP
【ＦＩ】
G06N99/00 180
【請求項の数】
9
(21)【出願番号】P 2019112547
(22)【出願日】2019-06-18
(65)【公開番号】P2020204928
(43)【公開日】2020-12-24
【審査請求日】2022-03-08
(73)【特許権者】
【識別番号】000005223
【氏名又は名称】富士通株式会社
(74)【代理人】
【識別番号】110002918
【氏名又は名称】弁理士法人扶桑国際特許事務所
(72)【発明者】
【氏名】田村 泰孝
(72)【発明者】
【氏名】此島 真喜子
【審査官】今城 朋彬
(56)【参考文献】
【文献】特許第６４６５２３１（ＪＰ，Ｂ１）
【文献】小熊祐司、外１名，電力ピークシフトを目的とした生産設備スケジューリング問題とその近似解法，電気学会研究会資料，日本，一般社団法人電気学会，2018年12月15日，pp.17-24
(58)【調査した分野】(Int.Cl.，ＤＢ名)
Ｇ０６Ｎ ３／００－９９／００
(57)【特許請求の範囲】
【請求項1】
エネルギーを表す評価関数に含まれる複数の状態変数の値と状態変数の組毎の重み値とを保持する状態保持部と、前記複数の状態変数の何れかの値が変化する場合に、前記複数の状態変数の値と前記重み値とに基づいて、前記複数の状態変数の値のそれぞれを次の変化候補とする場合のエネルギーの変化値を計算するエネルギー変化計算部と、
前記複数の状態変数に対して計算された複数のエネルギーの変化値のそれぞれに、不等式制約における前記複数の状態変数のそれぞれの重みに応じた接続係数と閾値とに基づいて算出される、前記不等式制約に違反する超過量に応じたペナルティ値を加算して総エネルギーの変化値を算出するペナルティ加算部と、
設定された温度値と乱数値と複数の前記総エネルギーの変化値とに基づいて、前記状態保持部に保持される前記複数の状態変数の何れかの値を変化させる更新制御部と、
を有し、
前記状態保持部は、前記複数の状態変数と前記不等式制約に対応する不等式制約用変数との組ごとに前記接続係数を保持し、
前記接続係数は、前記不等式制約用変数の第１の局所場における前記複数の状態変数の重みを０以外の値とし、前記複数の状態変数の局所場における前記不等式制約用変数の重みを０とする非対称接続係数である、
最適化装置。
【請求項2】
前記ペナルティ加算部は、前記接続係数と前記閾値とに基づいて算出される前記超過量の変化値に前記不等式制約の重要度を示す重要度係数を乗じることで前記ペナルティ値を算出する、請求項１記載の最適化装置。
【請求項3】
前記不等式制約が上限制約値に対する不等式制約である場合、前記ペナルティ加算部は、前記複数のエネルギーの変化値のそれぞれに、第１の接続係数と第１の閾値とに基づいて算出される、前記第１の閾値に対する正方向の第１の超過量の変化値に応じた第１のペナルティ値を加算する、
請求項１または２記載の最適化装置。
【請求項4】
前記不等式制約が下限制約値に対する不等式制約である場合、前記ペナルティ加算部は、前記複数のエネルギーの変化値のそれぞれに、第２の接続係数と第２の閾値とに基づいて算出される、前記第２の閾値に対する負方向の第２の超過量の変化値に応じた第２のペナルティ値を加算する、
請求項１または２記載の最適化装置。
【請求項5】
前記不等式制約が上限制約値および下限制約値に対する不等式制約である場合、前記ペナルティ加算部は、前記複数のエネルギーの変化値のそれぞれに、第３の接続係数と第３の閾値とに基づいて算出される、前記第３の閾値に対する正方向の第３の超過量の変化値、および、前記第３の接続係数と前記第３の閾値よりも小さい第４の閾値とに基づいて算出される、前記第４の閾値に対する負方向の第４の超過量の変化値のうち、大きい方の変化値に応じた第３のペナルティ値を加算する、
請求項１または２記載の最適化装置。
【請求項6】
所定の制約値に対する等式制約の場合、前記ペナルティ加算部は、前記複数のエネルギーの変化値のそれぞれに、第４の接続係数と第５の閾値とに基づいて算出される、前記第５の閾値に対する正方向の第５の超過量の変化値、および、前記第５の閾値に対する負方向の第６の超過量の変化値のうち、大きい方の変化値に応じた第４のペナルティ値を加算する、
請求項１または２記載の最適化装置。
【請求項7】
前記更新制御部により選択された変化対象の状態変数の識別情報に基づいて、前記エネルギー変化計算部により出力される前記複数のエネルギーの変化値のうちの何れかの変化値を選択し、選択した変化値を累積することで前記状態保持部に保持される前記複数の状態変数の値に対応するエネルギーの値を計算するエネルギー計算部を更に有する、請求項１乃至６の何れか１項に記載の最適化装置。
【請求項8】
それぞれが前記状態保持部と前記エネルギー変化計算部と前記ペナルティ加算部と前記更新制御部とを備え、互いに異なる前記温度値が設定される複数の探索部と、前記複数の探索部のそれぞれによる基底状態探索の繰り返し回数到達または一定時間経過後に、前記複数の探索部間で前記温度値または前記複数の状態変数の値を入れ替える交換制御部と、
を更に有する請求項１乃至７の何れか１項に記載の最適化装置。
【請求項9】
エネルギー変化計算部が、状態保持部に保持される複数の状態変数の何れかの値が変化する場合に、前記複数の状態変数の値と前記状態保持部に保持される重み値とに基づいて、前記複数の状態変数の値のそれぞれを次の変化候補とする場合のエネルギーの変化値を計算し、ペナルティ加算部が、前記複数の状態変数に対して計算された複数のエネルギーの変化値のそれぞれに、不等式制約における前記複数の状態変数のそれぞれの重みに応じた接続係数と閾値とに基づいて算出される、前記不等式制約に違反する超過量に応じたペナルティ値を加算して総エネルギーの変化値を算出し、
更新制御部が、設定された温度値と乱数値と複数の前記総エネルギーの変化値とに基づいて、前記状態保持部に保持される前記複数の状態変数の何れかの値を変化させ、
前記状態保持部は、前記複数の状態変数と前記不等式制約に対応する不等式制約用変数との組ごとに前記接続係数を保持し、
前記接続係数は、前記不等式制約用変数の第１の局所場における前記複数の状態変数の重みを０以外の値とし、前記複数の状態変数の局所場における前記不等式制約用変数の重みを０とする非対称接続係数である、
最適化方法。
【発明の詳細な説明】
【技術分野】
【0001】
本発明は最適化装置および最適化方法に関する。
【背景技術】
【0002】
ノイマン型コンピュータが不得意とする多変数の最適化問題を解く方法として、イジング型のエネルギー関数（評価関数や目的関数とも呼ばれる）を用いた最適化装置（イジングマシンやボルツマンマシンとも呼ばれる）がある。最適化装置は、計算対象の問題を、磁性体のスピンの振る舞いを表すモデルであるイジングモデルに置き換えて計算する。
【0003】
最適化装置は、例えば、ニューラルネットワークを用いてモデル化することもできる。その場合、イジングモデルに含まれる複数のスピンに対応した複数の状態変数のそれぞれが、他の状態変数の値と、他の状態変数と自身の状態変数との相互作用の大きさを示す重み係数とに応じて０または１を出力するニューロンとして機能する。最適化装置は、例えば、シミュレーテッド・アニーリングなどの確率的探索法により、上記のようなエネルギー関数の値（以下、エネルギーと言う）の最小値が得られる各状態変数の値の組合せを、解として求める。
【0004】
ここで、例えば、デジタル回路を用いてシミュレーテッド・アニーリングを行うことでエネルギーが最小となる各状態変数の値の組合せを計算する最適化装置の提案がある（例えば、特許文献１参照）。また、イジングモデルにおける各スピン間の相互作用を模擬するイジングチップとイジングチップを制御する情報処理部とを有する情報処理装置の提案もある（例えば、特許文献２参照）。
【0005】
なお、二値二次計画問題を半正定計画問題に緩和して、半正定計画問題の解を導出し、導出された半正定計画問題の解を二値二次計画問題の解に変換することで最適解を導出する最適化システムの提案もある（例えば、特許文献３参照）。
【先行技術文献】
【特許文献】
【0006】
【文献】特開２０１８－４１３５１号公報
【文献】国際公開第２０１７／０１７８０７号
【文献】国際公開第２０１７／０５６３６６号
【発明の概要】
【発明が解決しようとする課題】
【0007】
上記の最適化装置は、状態変数の二次形式で表される評価関数を最小化する。一方、最適化問題では、ある量（例えば、積載量や予算など）に対して上限や下限を設けるといった、不等式で表される制約条件（不等式制約）が課されることがある。しかし、こうした不等式制約を状態変数の二次形式で定式化した問題を、上記の最適化装置で解くのが難しいことがある。
【0008】
１つの側面では、本発明は、不等式制約を含む問題の求解を効率化する最適化装置および最適化方法を提供することを目的とする。
【課題を解決するための手段】
【0009】
１つの態様では、最適化装置が提供される。最適化装置は、状態保持部とエネルギー変化計算部とペナルティ加算部と更新制御部とを有する。状態保持部は、エネルギーを表す評価関数に含まれる複数の状態変数の値と状態変数の組毎の重み値とを保持する。エネルギー変化計算部は、複数の状態変数の何れかの値が変化する場合に、複数の状態変数の値と重み値とに基づいて、複数の状態変数の値のそれぞれを次の変化候補とする場合のエネルギーの変化値を計算する。ペナルティ加算部は、複数の状態変数に対して計算された複数のエネルギーの変化値のそれぞれに、不等式制約における複数の状態変数のそれぞれの重みに応じた接続係数と閾値とに基づいて算出される、不等式制約に違反する超過量に応じたペナルティ値を加算して総エネルギーの変化値を算出する。更新制御部は、設定された温度値と乱数値と複数の総エネルギーの変化値とに基づいて、状態保持部に保持される複数の状態変数の何れかの値を変化させる。状態保持部は、複数の状態変数と不等式制約に対応する不等式制約用変数との組ごとに接続係数を保持する。接続係数は、不等式制約用変数の第１の局所場における複数の状態変数の重みを０以外の値とし、複数の状態変数の局所場における不等式制約用変数の重みを０とする非対称接続係数である。
【0010】
また、１つの態様では、最適化方法が提供される。
【発明の効果】
【0011】
１つの側面では、不等式制約を含む問題の求解を効率化できる。
【図面の簡単な説明】
【0012】
【図1】第１の実施の形態の最適化装置を示す図である。
【図2】関数Ｅ（ｘ）、Ｃ（ｘ）の例を示す図である。
【図3】変数間の関係の例を示す図である。
【図4】最適化装置の回路構成例を示す図である。
【図5】最適化装置の回路のブロック構成例を示す図である。
【図6】判定部の回路構成例を示す図である。
【図7】選択部の回路構成例を示す図である。
【図8】変数毎のΔＥおよびΔＣの計算例を示す図である。
【図9】関数Ｅ（ｘ）に対応するエネルギーを計算する回路構成例を示す図である。
【図10】最適化装置の演算例を示すフローチャートである。
【図11】求解結果の比較の例を示す図である。
【図12】最適化システムの例を示す図である。
【図13】最適化装置の他の回路構成例を示す図である。
【図14】第２の実施の形態の最適化装置の例を示す図である。
【発明を実施するための形態】
【0013】
以下、本実施の形態について図面を参照して説明する。
［第１の実施の形態］
第１の実施の形態を説明する。
［第１の実施の形態］
第１の実施の形態を説明する。
【0014】
図１は、第１の実施の形態の最適化装置を示す図である。
最適化装置１０は、計算対象の問題を変換したイジングモデルに含まれる複数のスピンに対応する複数の状態変数のそれぞれの値の組合せ（状態）のうち、エネルギー関数が最小値となるときの各状態変数の値（基底状態）を探索する。状態変数は、「バイナリ変数」や「ビット（スピンビット）」と呼ばれてもよい。
最適化装置１０は、計算対象の問題を変換したイジングモデルに含まれる複数のスピンに対応する複数の状態変数のそれぞれの値の組合せ（状態）のうち、エネルギー関数が最小値となるときの各状態変数の値（基底状態）を探索する。状態変数は、「バイナリ変数」や「ビット（スピンビット）」と呼ばれてもよい。
【0015】
最適化問題が不等式制約を含まない場合、最適化装置１０は、次のように、エネルギー関数Ｅ（ｘ）に基づいて基底状態の探索を行える。イジング型のエネルギー関数Ｅ（ｘ）は、例えば以下の式（１）で定義される。
【0016】
【数1】
【0017】
式（１）に基づく制約なしの状態変数（バイナリ変数）の二次形式の関数の最適化問題は、（ＱＵＢＯ：Quadratic Unconstraint Binary Optimization）と呼ばれることがある。また、この二次形式での定式化は、ＱＵＢＯ形式と呼ばれることがある。
【0018】
式（１）の右辺第１項は、全状態変数から選択可能な２つの状態変数の全組合せについて、漏れと重複なく、２つの状態変数の値と重み値との積を積算したものである。ｘｉは、ｉ番目の状態変数である。ｘｊは、ｊ番目の状態変数である。重み値Ｗｉｊは、ｉ番目の状態変数とｊ番目の状態変数との間の重み（例えば、結合の強さ）を示す。なお、行列Ｗ＝｛Ｗｉｊ｝について、Ｗｉｊ＝Ｗｊｉ、Ｗｉｉ＝０であることが多い。状態変数ｘｉなどの変数に付加される添え字ｉは、当該変数の識別情報であり、インデックスと呼ばれる。
【0019】
式（１）の右辺第２項は、全状態変数のそれぞれのバイアス値と状態変数の値との積の総和である。ｂｉは、ｉ番目の状態変数に対するバイアス値を示す。
例えば、イジングモデルにおけるスピンの「－１」は、状態変数の値「０」に対応する。イジングモデルにおけるスピンの「＋１」は、状態変数の値「１」に対応する。
例えば、イジングモデルにおけるスピンの「－１」は、状態変数の値「０」に対応する。イジングモデルにおけるスピンの「＋１」は、状態変数の値「１」に対応する。
【0020】
状態変数ｘｉの値が変化して１－ｘｉとなると、状態変数ｘｉの増加分は、δｘｉ＝（１－ｘｉ）－ｘｉ＝１－２ｘｉと表せる。したがって、エネルギー関数Ｅ（ｘ）に対して、状態変数ｘｉのスピン反転（値の変化）に伴うエネルギー変化ΔＥｉは、式（２）で表される。
【0021】
【数2】
【0022】
ｈｉは局所場（ローカルフィールド）と呼ばれ、式（３）で表される。
【0023】
【数3】
【0024】
状態変数ｘｊが変化（ビット反転）したときの局所場ｈｉの変化分δｈｉ
（ｊ）は、式（４）で表される。
【0025】
【数4】
【0026】
最適化装置１０は、局所場ｈｉを保存するレジスタ（図示を省略している）を有し、状態変数ｘｊの値が変化したときに変化分δｈｉ
（ｊ）を、ｈｉに加算することで、ビット反転後の状態に対応するｈｉを得る。
【0027】
最適化装置１０では、基底状態の探索において、エネルギー変化がΔＥｉとなる状態遷移（状態変数ｘｉの値の変化）を許容するか否かを決定するためにメトロポリス法やギブス法が用いられる。すなわち、最適化装置１０は、ある状態から当該状態よりもエネルギーの低い他の状態への遷移を探索する近傍探索において、エネルギーが下がる状態だけでなく、エネルギーが上がる状態への遷移を確率的に許容する。例えば、エネルギー変化ΔＥの状態変数の値の変化を受け入れる確率（遷移受入確率）Ａは、式（５）で表される。
【0028】
【数5】
【0029】
ここで、βは温度Ｔの逆数（β＝１／Ｔ）である。ｍｉｎ演算子は、引数のうちの最小値を取ることを示す。例えば、メトロポリス法を用いる場合、式（５）におけるＡ＝ｅｘｐ（－β・ΔＥ）に対して、両辺の自然対数を取って変形すると式（６）を得る。
【0030】
【数6】
【0031】
したがって、最適化装置１０は、一様乱数ｕ（０＜ｕ≦１）に対して、エネルギー変化ΔＥが式（７）を満たす場合に、該当の状態変数の値の変化を許容する。
【0032】
【数7】
【0033】
そして、最適化装置１０は、変化が許容された何れかの状態変数の値を変化させる。最適化装置１０は、温度Ｔを初期温度から最低の温度まで下げながら、各温度において状態変数の値を変化させる処理を繰り返し実行することで、最適化問題に対する解を求める。
【0034】
最適化問題には不等式制約が付加されることがある。そこで、最適化装置１０は、不等式制約に対する演算機能を提供する。最適化装置１０は、状態保持部１１、エネルギー変化計算部１２、ペナルティ加算部１３、更新制御部１４、制約超過量計算部１５および制約識別符号入力部１６を有する。最適化装置１０は、例えば、ＦＰＧＡ（Field Programmable Gate Array）などの半導体集積回路を用いて実現される。最適化装置１０は、半導体チップとして実現され得る。
【0035】
状態保持部１１は、エネルギーを表す評価関数（上記のエネルギー関数Ｅ（ｘ））に含まれる複数の状態変数の値と状態変数の組毎の重み値とを保持する。図１では、複数の状態変数の値を、状態ベクトルｘで表す。例えば、複数の状態変数の数は、ｎ個（ｎは２以上の整数）である（すなわち、状態ベクトルｘはｎビットである）。
【0036】
エネルギー変化計算部１２は、状態保持部１１に保持される複数の状態変数の何れかの値が変化されると、複数の状態変数の値と重み値とに基づいて、複数の状態変数の値のそれぞれを次の変化候補とする場合のエネルギーの変化値を計算する。例えば、エネルギー変化計算部１２は、変化候補を状態変数ｘｊとする場合のエネルギーの変化値ΔＥｊを、ｊ＝１～ｎのそれぞれに対して求める。
【0037】
ペナルティ加算部１３は、複数の状態変数に対して計算された複数のエネルギーの変化値ΔＥｊのそれぞれに、不等式制約に違反する超過量（制約超過量）に応じたペナルティ値ΔＣｊを加算して総エネルギーの変化値（ΔＥｊ＋ΔＣｊ）を算出する。制約超過量は、不等式制約における複数の状態変数のそれぞれの重みを示す接続係数と閾値とに基づいて算出される値である。例えば、変化候補の状態変数ｘｊに対応する制約超過量の変化値が、制約超過量計算部１５により計算され、ペナルティ加算部１３に供給される。ペナルティ加算部１３は、状態変数ｘｊに対応する制約超過量の変化値に基づいてペナルティ値ΔＣｊを求め、ΔＥｊに加算する。例えば、ペナルティ加算部１３は、当該制約超過量の変化値に、不等式制約に対する重要度を表す所定の係数を乗じることでペナルティ値ΔＣｊを算出することが考えられる。
【0038】
更新制御部１４は、設定された温度値と乱数値と複数の総エネルギーの変化値とに基づいて、状態保持部１１に保持される複数の状態変数の何れかの値を変化させる。更新制御部１４は、複数の総エネルギーの変化値と熱励起エネルギー（熱ノイズ）との相対関係によって複数の状態変数の何れかを受け入れるか否かを確率的に決定するともいえる。更新制御部１４は、遷移可否判定部１４ａおよび選択部１４ｂを有する。
【0039】
遷移可否判定部１４ａは、温度値Ｔと乱数値ｕと総エネルギーの変化値（ΔＥｊ＋ΔＣｊ）とに基づいて、当該総エネルギーの変化値（ΔＥｊ＋ΔＣｊ）に対応する状態遷移を許容するか否かを示すフラグＦｊを出力する。例えば、状態変数ｘｊの値の変化による状態遷移を許容する場合にＦｊ＝１であり、状態変数ｘｊの値の変化による状態遷移を許容しない場合にＦｊ＝０である。なお、温度値Ｔは、図示を省略している制御部（あるいは、制御回路）により、遷移可否判定部１４ａに入力される。ここで、遷移可否判定部１４ａによる判定は、式（７）のΔＥを、総エネルギーの変化値（ΔＥｊ＋ΔＣｊ）に置き換えた式（８）に基づいて行われる。
【0040】
【数8】
【0041】
選択部１４ｂは、遷移可否判定部１４ａにより出力されたＦｊ（ｊ＝１～ｎ）のうち、Ｆｊ＝１に対応する状態遷移の番号（インデックス）を１つ選択する。選択部１４ｂは、選択したインデックスを状態保持部１１に出力することで、状態保持部１１に保持される複数の状態変数の何れかの値を変化させる。例えば、状態保持部１１は、選択部１４ｂから受け付けたインデックスに対応する状態変数の値を変化（ビット反転）させる。こうして、状態保持部１１に入力されるクロック信号（ｃｌｋ）に同期して、状態保持部１１に保持される状態変数の値が更新される。
【0042】
制約超過量計算部１５は、不等式制約における閾値と、不等式制約に対する状態変数ｘｊの接続係数とに基づいて、状態変数ｘｊを変化候補とした場合の制約超過量の変化値を計算し、ペナルティ加算部１３に供給する。閾値は、不等式制約における上限値でもよいし、下限値でもよい。閾値は、上限値および下限値の両方でもよく、その場合、不等式制約に違反する制約超過量の変化値として、上限値に対する制約超過量の変化値および下限値に対する制約超過量の変化値のうちの大きい方を採用することが考えられる。
【0043】
制約識別符号入力部１６は、制約識別符号を制約超過量計算部１５に入力する。制約識別符号は、閾値として用いられる上限値や下限値を指定する情報である。制約識別符号は、前述の制御部により、制約識別符号入力部１６に入力される。
【0044】
また、最適化装置１０で保持可能な状態変数の最大数をＮ（Ｎはｎより大きい整数）として、Ｎ個の変数のうちｎ個を通常の状態変数として用い、Ｋ（Ｋは１以上の整数）個を不等式制約を表すための変数（不等式制約用変数）として用いる構成が考えられる。Ｎ＝ｎ＋Ｋである。Ｋは、不等式制約の数に対応し、不等式制約の数が１個であればＫ＝１、不等式制約の数が２個であればＫ＝２などとなる。
【0045】
その場合、制約識別符号入力部１６は、Ｎ個の変数のうち、状態ベクトルに属する通常の状態変数と、状態ベクトルに属さない不等式制約用変数とを、制約識別符号により指定する。そのため、制約識別符号入力部１６は、エネルギー変化計算部１２、ペナルティ加算部１３および遷移可否判定部１４ａに対して制約識別符号を入力する構成としてもよい（図１では当該構成を示している）。このようにすると、後述されるように、Ｋ個の不等式制約用変数のそれぞれに対応する局所場ｈｉを記憶するレジスタに格納された現在のｈｉの値を用いてペナルティ値ΔＣｉを計算することができる。
【0046】
Ｎ個の変数のうち、ｎ個を通常の状態変数として用い、Ｋ個を不等式制約用変数として用いる場合、不等式制約が付加された最適化問題を次のように定式化できる。ここで、通常の状態変数のインデックスの集合をΨと表す。集合Ψの要素数はｎ個である。不等式制約用変数のインデックスの集合をΩ（Ω＝｛ｋ１，ｋ２，…，ｋＫ｝）と表す。集合Ωの要素数はＫ個である。
【0047】
最小化するエネルギー関数は、ｎ個の状態変数の関数Ｅ（ｘ）として次の式（９）で表される。
【0048】
【数9】
【0049】
ｎ個の状態変数から構成されるＫ個の不等式制約は式（１０）で表される。
【0050】
【数10】
【0051】
接続係数ａｉｊは、インデックスｉに対応する不等式制約（不等式制約ｉと表記する）における状態変数ｘｊ（ｊ∈Ψ）の重みを示す。Ｃｌｉは、不等式制約ｉにおける下限値である。Ｃｕｉは、不等式制約ｉにおける上限値である。ただし、不等式制約として、下限値または上限値の何れか一方のみを設けてもよい。
【0052】
このとき、総エネルギーに対応する目的関数Ｅｔｏｔ（ｘ）を式（１１）で表す。
【0053】
【数11】
【0054】
関数Ｃ（ｘ）は、不等式制約に対する制約超過量と重要度係数λｉの積を、全ての不等式制約について和をとったものであり、式（１２）で表される。
【0055】
【数12】
【0056】
ここで、重要度係数λｉは、０以上の実数であり、不等式制約ｉの重要度を表す。また、ｍａｘ演算子は、引数のうちの最大値をとることを示す。関数Ｃ（ｘ）の変化値がペナルティ値に相当する。
【0057】
図２は、関数Ｅ（ｘ）、Ｃ（ｘ）の例を示す図である。
図２（Ａ）は、関数Ｅ（ｘ）の一例のグラフＧ１を示す。図２（Ｂ）は、関数Ｃ（ｘ）の一例のグラフＧ２を示す。グラフＧ１，Ｇ２の横軸は、状態ｘ（状態ベクトルで表される状態）を示す。グラフＧ１の縦軸はＥ（ｘ）である。グラフＧ２の縦軸はＣ（ｘ）である。関数Ｅ（ｘ）に関数Ｃ（ｘ）を加えた目的関数Ｅｔｏｔ（ｘ）に関して、例えば、式（８）（メトロポリス法）に基づいて状態を遷移させることで、関数Ｅ（ｘ）に対して不等式制約を充足する（あるいは不等式制約に関する制約超過を最小化する）解を得られる。
図２（Ａ）は、関数Ｅ（ｘ）の一例のグラフＧ１を示す。図２（Ｂ）は、関数Ｃ（ｘ）の一例のグラフＧ２を示す。グラフＧ１，Ｇ２の横軸は、状態ｘ（状態ベクトルで表される状態）を示す。グラフＧ１の縦軸はＥ（ｘ）である。グラフＧ２の縦軸はＣ（ｘ）である。関数Ｅ（ｘ）に関数Ｃ（ｘ）を加えた目的関数Ｅｔｏｔ（ｘ）に関して、例えば、式（８）（メトロポリス法）に基づいて状態を遷移させることで、関数Ｅ（ｘ）に対して不等式制約を充足する（あるいは不等式制約に関する制約超過を最小化する）解を得られる。
【0058】
図３は、変数間の関係の例を示す図である。
変数群ｇ１は、状態変数の集合｛ｘｊ｝（ｊ∈Ψ）を示す。変数群ｇ２は、不等式制約用変数の集合｛ｘｋ｝（ｋ∈Ω）を示す。
変数群ｇ１は、状態変数の集合｛ｘｊ｝（ｊ∈Ψ）を示す。変数群ｇ２は、不等式制約用変数の集合｛ｘｋ｝（ｋ∈Ω）を示す。
【0059】
状態変数のペアは、双方向のカップリングである。すなわち、当該ペアの何れか一方の状態変数の反転に対して、当該ペアの他方の状態変数の局所場が影響を受ける関係となる。
【0060】
状態変数と不等式制約用変数とのペアは、状態変数から不等式制約用変数へ向かう単方向のカップリングである。すなわち、状態変数の反転に対して、不等式制約用変数の局所場が影響を受けるが、不等式制約用変数の反転に対して、状態変数の局所場は影響を受けない関係となる。不等式制約ｋに対する状態変数ｘｊの影響の大きさ（重み）は、接続係数ａｋｊにより表される。接続係数ａｋｊは、前述の行列Ｗの要素として与えられる。また、不等式制約用変数のバイアスｂｋは、ｂｋ＝０が設定される。
【0061】
状態変数と不等式制約用変数との関係を実装するため、一例では、最適化装置１０は、不等式制約用変数ｘｋの値を０に設定し、不等式制約用変数の反転が生じないように制御することで、状態変数の局所場への影響を抑制する。他の例では、接続係数ａｋｊを非対称接続係数（すなわち、ａｋｊ≠０，ａｊｋ＝０）とすることで、不等式制約用変数ｘｋによる状態変数の局所場への影響を抑制し得る。更に他の例では、遷移可否判定部１４ａにより、制約識別符号に基づいて、不等式制約用変数ｘｋに対して常にフラグＦｋ＝０を出力するように制御することも考えられる。
【0062】
不等式制約用変数は他の不等式制約用変数とカップリングされない。すなわち、不等式制約用変数のペアに対しては、Ｗｉｊ，Ｗｊｉ（ｉ，ｊ∈Ω）を０とする。
不等式制約用変数ｘｋ（ｋ∈Ω）の局所場ｈｋは式（１３）で表される。
不等式制約用変数ｘｋ（ｋ∈Ω）の局所場ｈｋは式（１３）で表される。
【0063】
【数13】
【0064】
不等式制約用変数は、通常の状態変数で表される状態のエネルギーＥ（ｘ）には影響を与えないが、不等式制約用変数ｘｋに対応する局所場ｈｋの値は、通常の状態変数の反転に伴って、最新の値に更新される。そして、局所場ｈｋと、不等式制約の上限・下限とを比較することで、不等式制約の充足状況が分かる。
【0065】
式（１１）で表されるように、総エネルギーＥｔｏｔ（ｘ）は、ＱＵＢＯ項Ｅ（ｘ）と不等式制約項（ペナルティ項）Ｃ（ｘ）との和となる。通常の状態変数ｘｊ（ｊ∈Ψ）の反転に対する総エネルギーのエネルギー変化ΔＥｔｏｔ ｊ（ｘ）は、式（１４）で表される。
【0066】
【数14】
【0067】
ΔＥｊは、例えば、特開２０１８－４１３５１号公報に開示された既存の回路構成を用いて、ＱＵＢＯ項に基づき生成される。
ΔＣｊは、不等式制約用変数ｘｉ（ｉ∈Ω）の局所場ｈｉの現在の値を用いて、式（１５）に基づき生成される。
ΔＣｊは、不等式制約用変数ｘｉ（ｉ∈Ω）の局所場ｈｉの現在の値を用いて、式（１５）に基づき生成される。
【0068】
【数15】
【0069】
ここで、ΔＣｉｊは、式（１６）により表される。
【0070】
【数16】
【0071】
ただし、δｘｊは、状態変数ｘｊ反転時の状態変数ｘｊの変化分であり、式（１７）で表される。
【0072】
【数17】
【0073】
すなわち、ΔＣｉｊは、状態変数ｘｊの反転に伴うインデックスｉの不等式制約（不等式制約ｉ）に対する制約超過量の変化分の値（変化値）である。また、ペナルティ値ΔＣｊは、不等式制約毎の制約超過量の変化値の、重要度係数による重み付き総和である。なお、不等式制約ｉで上限値のみを与える場合、式（１５）、（１６）のｍａｘ演算において、「Ｃｌｉ－ａｉｊδｘｊ－ｈｉ」の引数や「Ｃｌｉ－ｈｉ」の引数を含まなくてよい。また、不等式制約ｉで下限値のみを与える場合、式（１５）、（１６）のｍａｘ演算において、「ｈｉ＋ａｉｊδｘｊ－Ｃｕｉ」の引数や「ｈｉ－Ｃｕｉ」の引数を含まなくてよい。更に、制約ｉとして等式制約を表す場合、式（１５）、（１６）において、Ｃｕｉ＝Ｃｌｉとすればよい。
【0074】
次に、最適化装置１０に対する上記原理の実装例を更に説明する。以下では、最適化装置１０で保持可能な状態変数の最大数をＮ個とする。また、Ωに属するインデックスについて、便宜的にｋ１を１、ｋ２を２，…，ｋＫをＫのように、ｋｉの添え字ｉで表記することがある。
【0075】
図４は、最適化装置の回路構成例を示す図である。
エネルギー変化計算部１２は、ΔＥ計算部１２ａ１～１２ａＮを有する。ΔＥ計算部１２ａ１～１２ａＮのそれぞれは、状態保持部１１から供給されるｊ番目の状態変数を反転候補とする場合のエネルギーの変化値ΔＥｊを、式（２）～（４）に基づいて計算し、ペナルティ加算部１３に出力する。ΔＥ計算部１２ａ１～１２ａＮのそれぞれには、局所場ｈ１～ｈＮを記憶するレジスタ（図示を省略している）から、局所場ｈ１～ｈＮの現在の値が供給される。
エネルギー変化計算部１２は、ΔＥ計算部１２ａ１～１２ａＮを有する。ΔＥ計算部１２ａ１～１２ａＮのそれぞれは、状態保持部１１から供給されるｊ番目の状態変数を反転候補とする場合のエネルギーの変化値ΔＥｊを、式（２）～（４）に基づいて計算し、ペナルティ加算部１３に出力する。ΔＥ計算部１２ａ１～１２ａＮのそれぞれには、局所場ｈ１～ｈＮを記憶するレジスタ（図示を省略している）から、局所場ｈ１～ｈＮの現在の値が供給される。
【0076】
なお、ΔＥ計算部１２ａ１～１２ａＮには、不等式制約用変数に対応するΔＥ計算部が含まれるが、後述されるように、不等式制約用変数に対応するΔＥ計算部は、制約識別符号に基づいて、比較的大きな値を出力することで、不等式制約用変数の反転を抑制する。
【0077】
ペナルティ加算部１３は、係数レジスタ１３ａ１～１３ａＫ、乗算器１３ｂ１～１３ｂＫ、加算器１３ｃ１１～１３ｃ１Ｎ，…，１３ｃＫ１～１３ｃＫＮを有する。ここで、係数レジスタ１３ａ１～１３ａＫ、乗算器１３ｂ１～１３ｂＫのように、不等式制約用変数のインデックスに対応するＫ個の係数レジスタおよび乗算器を図示しているが、当該係数レジスタおよび乗算器は、それぞれＮ個設けられる。すなわち、不等式制約用変数のインデックスに対応していない他の係数レジスタおよび他の乗算器の図示を省略している。同様に、加算器１３ｃ１１～１３ｃ１Ｎ，…，１３ｃＫ１～１３ｃＫＮのように、不等式制約用変数に対応するＫ群（あるいはＫ行）の加算器を図示しているが、当該加算器群は、Ｎ群（あるいはＮ行）設けられる。すなわち、不等式制約用変数のインデックスに対応していない他の加算器群の図示を省略している。
【0078】
係数レジスタ１３ａ１～１３ａＫは、式（１２）、（１５）における重要度係数λ１～λＫを保持し、係数レジスタ１３ａ１～１３ａＫに対応する乗算器１３ｂ１～１３ｂＫに供給する。
【0079】
乗算器１３ｂ１～１３ｂＫは、制約超過量計算部１５により計算されたΔＣ１ｊ～ΔＣＫｊ（ｊ∈Ψ）に、係数λ１～λＫを乗じて、加算器１３ｃ１１～１３ｃ１Ｎ，…，１３ｃＫ１～１３ｃＫＮに供給する。
【0080】
具体的には、乗算器１３ｂ１は、ΔＣ１１にλ１を乗じて、加算器１３ｃ１１に供給する。また、乗算器１３ｂ１は、ΔＣ１２にλ１を乗じて、加算器１３ｃ１２に供給する。以降、同様に、乗算器１３ｂ１は、ΔＣ１Ｎにλ１を乗じて、加算器１３ｃ１Ｎに供給する。
【0081】
更に、乗算器１３ｂＫは、ΔＣＫ１にλＫを乗じて、加算器１３ｃＫ１に供給する。また、乗算器１３ｂＫは、ΔＣＫ２にλＫを乗じて、加算器１３ｃＫ２に供給する。以降、同様に、乗算器１３ｂＫは、ΔＣＫＮにλＫを乗じて、加算器１３ｃＫＮに供給する。他の乗算器についても同様である。
【0082】
加算器１３ｃ１１～１３ｃ１Ｎ，…，１３ｃＫ１～１３ｃＫＮは、乗算器１３ｂ１～１３ｂＫから供給されるλ１ΔＣ１１～λ１ΔＣ１Ｎ，…，λＫΔＣＫ１～λＫΔＣＫＮを、エネルギー変化計算部１２で計算されたΔＥ１～ΔＥＮに加算する。
【0083】
具体的には、加算器１３ｃ１１は、ΔＥ計算部１２ａ１から供給されるΔＥ１にλ１ΔＣ１１を加算し、インデックス１に対応する次段の加算器（図示を省略している）に供給する。加算器１３ｃ１２は、ΔＥ計算部１２ａ２から供給されるΔＥ２にλ１ΔＣ１２を加算し、インデックス２に対応する次段の加算器に供給する。以降、同様に、加算器１３ｃ１Ｎは、ΔＥ計算部１２ａＮから供給されるΔＥＮにλ１ΔＣ１Ｎを加算し、インデックスＮに対応する次段の加算器に供給する。
【0084】
更に、加算器１３ｃＫ１は、インデックス１に対応する前段の加算器（図示を省略している）から供給される値にλＫΔＣＫ１を加算し、遷移可否判定部１４ａに出力する。加算器１３ｃＫ２は、インデックス２に対応する前段の加算器から供給される値にλＫΔＣＫ２を加算し、遷移可否判定部１４ａに出力する。以降、同様に、加算器１３ｃＫＮは、インデックスＮに対応する前段の加算器から供給される値にλＫΔＣＫＮを加算し、遷移可否判定部１４ａに出力する。
【0085】
遷移可否判定部１４ａは、判定部１４ａ１～１４ａＮを有する。判定部１４ａ１～１４ａＮのそれぞれは、ペナルティ加算部１３から供給されるΔＥ１＋ΔＣ１～ΔＥＮ＋ΔＣＮに基づいて、状態変数ｘｊの反転を許容するか否かを、式（８）により確率的に判定し、反転可否を示すフラグＦ１～ＦＮを選択部１４ｂに供給する。判定部１４ａ１～１４ａＮの回路構成の詳細は後述される。
【0086】
選択部１４ｂは、遷移可否判定部１４ａから供給されるフラグＦ１～ＦＮに基づいて、反転させる状態変数を選択し、当該状態変数のインデックスを状態保持部１１に供給する。選択部１４ｂの回路構成の詳細は後述される。
【0087】
制約超過量計算部１５は、ΔＣ計算部１５ａ１～１５ａＫを有する。ΔＣ計算部１５ａ１～１５ａＫは、式（１６）に基づいて、ΔＣ１ｊ～ΔＣＫｊ（ｊ∈Ψ）を計算し、ペナルティ加算部１３に供給する。
【0088】
具体的には、ΔＣ計算部１５ａ１は、ΔＣ１１～ΔＣ１Ｎを、乗算器１３ｂ１に供給する。以降、同様に、ΔＣ計算部１５ａＫは、ΔＣＫ１～ΔＣＫＮを乗算器１３ｂＫに供給する。
【0089】
上記のように、最適化装置１０の機能は、電子回路によって実現できる。例えば、状態保持部１１、エネルギー変化計算部１２、ペナルティ加算部１３、遷移可否判定部１４ａ、選択部１４ｂおよび制約超過量計算部１５を、それぞれ状態保持回路、エネルギー変化計算回路、ペナルティ加算回路、遷移可否判定回路、選択回路および制約超過量計算回路のように呼ぶこともできる。
【0090】
図５は、最適化装置の回路のブロック構成例を示す図である。
図４で例示した回路要素は、インデックスｉ＝１，２，…，Ｎのそれぞれに対応するブロックＢｉに分けて実装することができる。図５では、最適化装置１０のｊ番目のブロックＢｊの回路構成例が示されている。
図４で例示した回路要素は、インデックスｉ＝１，２，…，Ｎのそれぞれに対応するブロックＢｉに分けて実装することができる。図５では、最適化装置１０のｊ番目のブロックＢｊの回路構成例が示されている。
【0091】
ブロックＢｊは、ΔＥ計算部１２ａｊ、係数レジスタ１３ａ１ｊ～１３ａＫｊ、乗算器１３ｂ１ｊ～１３ｂＫｊ、加算器１３ｃ１ｊ～１３ｃＫｊ、判定部１４ａｊおよびΔＣ計算部１５ａ１ｊ～１５ａＫｊを有する。
【0092】
ΔＥ計算部１２ａｊは、インデックスｊに対応する状態変数ｘｊの値、および、局所場ｈｊの値に基づいて、ΔＥｊを計算し、加算器１３ｃ１ｊに供給する。
係数レジスタ１３ａ１ｊ～１３ａＫｊは、係数λ１～λＫを保持する。係数レジスタ１３ａ１ｊ～１３ａＫｊをブロック毎に設けるのではなく、各ブロックで共有のレジスタとしてもよい。
係数レジスタ１３ａ１ｊ～１３ａＫｊは、係数λ１～λＫを保持する。係数レジスタ１３ａ１ｊ～１３ａＫｊをブロック毎に設けるのではなく、各ブロックで共有のレジスタとしてもよい。
【0093】
乗算器１３ｂ１ｊ～１３ｂＫｊは、それぞれ係数レジスタ１３ａ１ｊ～１３ａＫｊに保持されるλ１～λＫを、ΔＣ計算部１５ａ１ｊ～１５ａＫｊから供給されるΔＣ１ｊ～ΔＣＫｊに乗算して、加算器１３ｃ１ｊ～１３ｃＫｊに供給する。乗算器１３ｂ１ｊ～１３ｂＫｊは、それぞれ乗算器１３ｂ１～１３ｂＫのうちのインデックスｊに対応する回路要素であると言える（例えば、乗算器１３ｂ１ｊは乗算器１３ｂ１のインデックスｊに対応する回路である）。
【0094】
加算器１３ｃ１ｊは、ΔＥ計算部１２ａｊから供給されるΔＥｊに、乗算器１３ｂ１ｊから供給されるλ１ΔＣ１ｊを加算して、次段の加算器（図示を省略している）に供給する。以降、同様にλｋＣｋｊ（ｋ∈Ω）が順次加算される。加算器１３ｃＫｊは、前段の加算器（図示を省略している）から供給される値にλＫΔＣＫｊを加算し、その結果（ΔＥｊ＋ΔＣｊに相当）を判定部１４ａｊに供給する。
【0095】
判定部１４ａｊは、総エネルギーの変化値ΔＥｊ＋ΔＣｊに基づいて、インデックスｊに対応する反転可否を判定し、判定結果に応じたフラグＦｊを選択部１４ｂに出力する。
ΔＣ計算部１５ａ１ｊ～１５ａＫｊは、それぞれΔＣ１ｊ～ΔＣＫｊを計算し、乗算器１３ｂ１ｊ～１３ｂＫｊに供給する。ΔＣ計算部１５ａ１ｊ～１５ａＫｊは、それぞれΔＣ計算部１５ａ１～１５ａＫのうちのインデックスｊに対応する回路要素であると言える（例えば、ΔＣ計算部１５ａ１ｊはΔＣ計算部１５ａ１のインデックスｊに対応する回路要素である）。
ΔＣ計算部１５ａ１ｊ～１５ａＫｊは、それぞれΔＣ１ｊ～ΔＣＫｊを計算し、乗算器１３ｂ１ｊ～１３ｂＫｊに供給する。ΔＣ計算部１５ａ１ｊ～１５ａＫｊは、それぞれΔＣ計算部１５ａ１～１５ａＫのうちのインデックスｊに対応する回路要素であると言える（例えば、ΔＣ計算部１５ａ１ｊはΔＣ計算部１５ａ１のインデックスｊに対応する回路要素である）。
【0096】
ここで、ブロックＢｊには、インデックスｊが通常の状態変数に対応するか、または、不等式制約用変数に対応するかを示す制約識別符号が、制約識別符号入力部１６により入力される。制約識別符号は、制約無し（不等式制約用変数に該当しない）、上限のみの制約、下限のみの制約、あるいは、上限および下限の両方の制約かを識別する識別コードを含む。
【0097】
例えば、識別コードｄｉは、２ビット（ｄｉ１，ｄｉ２）で表される。識別コードｄｉ＝００は、インデックスｉが通常の状態変数のインデックスであることを示す。識別コードｄｉ＝０１は、インデックスｉが上限のみの制約（あるいは上限のみの制約に対応する不等式制約用変数）のインデックスであることを示す。識別コードｄｉ＝１０は、インデックスｉが下限のみの制約（あるいは下限のみの制約に対応する不等式制約用変数）のインデックスであることを示す。識別コードｄｉ＝１１は、インデックスｉが上限および下限の両方の制約（あるいは上限および下限の両方の制約に対応する不等式制約用変数）のインデックスであることを示す。
【0098】
制約識別符号は、識別コードに対する上限値および下限値の指定も含む。例えば、制約識別符号は、（識別コード，上限値，下限値）の組となる。
等式制約の場合は、識別コードｄｉ＝１１を用い、前述のように、上限値および下限値を一致させる。ΔＥ計算部１２ａｊは、識別コードに応じて、インデックスｊに対するΔＥｊを計算する。ΔＣ計算部１５ａ１ｊ～１５ａＫｊは、識別コードに応じて、インデックスｊに対するΔＣ１ｊ～ΔＣＫｊを計算する。識別コードに応じたΔＥｊおよびΔＣ１ｊ～ΔＣＫｊの計算方法の例は後述される。
等式制約の場合は、識別コードｄｉ＝１１を用い、前述のように、上限値および下限値を一致させる。ΔＥ計算部１２ａｊは、識別コードに応じて、インデックスｊに対するΔＥｊを計算する。ΔＣ計算部１５ａ１ｊ～１５ａＫｊは、識別コードに応じて、インデックスｊに対するΔＣ１ｊ～ΔＣＫｊを計算する。識別コードに応じたΔＥｊおよびΔＣ１ｊ～ΔＣＫｊの計算方法の例は後述される。
【0099】
図６は、判定部の回路構成例を示す図である。
判定部１４ａｊは、オフセット値生成部２１、乱数生成部２２、ノイズ値生成部２３、符号反転回路２４、加算器２５，２６および比較器２７を有する。
判定部１４ａｊは、オフセット値生成部２１、乱数生成部２２、ノイズ値生成部２３、符号反転回路２４、加算器２５，２６および比較器２７を有する。
【0100】
オフセット値生成部２１は、選択部１４ｂから出力される遷移可否を示すフラグに基づいて、オフセット値Ｅｏｆｆ（Ｅｏｆｆ≧０）を生成し、加算器２５に供給する。具体的には、オフセット値生成部２１は、選択部１４ｂから出力されるフラグが遷移可を示す場合、オフセット値を０にリセットする。オフセット値生成部２１は、選択部１４ｂから出力されるフラグが遷移不可を示す場合、オフセット値に増分値ΔＥｏｆｆを加算する。当該フラグが連続して遷移不可を示す場合、オフセット値生成部２１は、ΔＥｏｆｆを積算することで、ＥｏｆｆをΔＥｏｆｆずつ増加させる。
【0101】
乱数生成部２２は、０＜ｕ≦１の一様乱数ｕを生成し、ノイズ値生成部２３に出力する。
ノイズ値生成部２３は、使用法則（例えば、メトロポリス法）に応じた所定の変換テーブルを保持する。ノイズ値生成部２３は、一様乱数ｕと、制御部（あるいは制御回路）により供給された温度Ｔを示す温度情報とを用いて、当該変換テーブルにより、式（８）に基づくノイズ値（熱ノイズ）に対応する－Ｔ・ｌｎ（ｕ）の値を生成する。ノイズ値生成部２３は、生成した－Ｔ・ｌｎ（ｕ）の値を加算器２６に出力する。
ノイズ値生成部２３は、使用法則（例えば、メトロポリス法）に応じた所定の変換テーブルを保持する。ノイズ値生成部２３は、一様乱数ｕと、制御部（あるいは制御回路）により供給された温度Ｔを示す温度情報とを用いて、当該変換テーブルにより、式（８）に基づくノイズ値（熱ノイズ）に対応する－Ｔ・ｌｎ（ｕ）の値を生成する。ノイズ値生成部２３は、生成した－Ｔ・ｌｎ（ｕ）の値を加算器２６に出力する。
【0102】
符号反転回路２４は、加算器１３ｃＫｊから供給される総エネルギーの変化値（ΔＥｊ＋ΔＣｊ）の符号を反転させ、加算器２５に供給する。
加算器２５は、符号反転回路２４から供給される－（ΔＥｊ＋ΔＣｊ）にオフセット値Ｅｏｆｆを加算し、加算器２６に供給する。
加算器２５は、符号反転回路２４から供給される－（ΔＥｊ＋ΔＣｊ）にオフセット値Ｅｏｆｆを加算し、加算器２６に供給する。
【0103】
加算器２６は、加算器２５から供給される－（ΔＥｊ＋ΔＣｊ）＋Ｅｏｆｆに、熱ノイズ－Ｔ・ｌｎ（ｕ）を加算し、比較器２７に供給する。
比較器２７は、加算器２６により出力された評価値（－（ΔＥｊ＋ΔＣｊ）＋Ｅｏｆｆ－Ｔ・ｌｎ（ｕ））を閾値（具体的には０）と比較することで、式（８）に基づく判定を行う。比較器２７は、評価値が０以上の場合、遷移可を示すフラグ（Ｆｊ＝１）を、選択部１４ｂに出力する。比較器２７は、評価値が０未満の場合、遷移不可を示すフラグ（Ｆｊ＝０）を、選択部１４ｂに出力する。
比較器２７は、加算器２６により出力された評価値（－（ΔＥｊ＋ΔＣｊ）＋Ｅｏｆｆ－Ｔ・ｌｎ（ｕ））を閾値（具体的には０）と比較することで、式（８）に基づく判定を行う。比較器２７は、評価値が０以上の場合、遷移可を示すフラグ（Ｆｊ＝１）を、選択部１４ｂに出力する。比較器２７は、評価値が０未満の場合、遷移不可を示すフラグ（Ｆｊ＝０）を、選択部１４ｂに出力する。
【0104】
ここで、選択部１４ｂから出力されるフラグが遷移不可を示す場合、現在の状態が局所解に陥っていると考えられる。オフセット値生成部２１による、－（ΔＥｊ＋ΔＣｊ）へのＥｏｆｆの加算やＥｏｆｆの漸増により、状態遷移が許容されやすくなり、現在の状態が局所解にある場合、その局所解からの脱出が促進される。
【0105】
次に、選択部１４ｂの回路構成例を説明する。
図７は、選択部の回路構成例を示す図である。
選択部１４ｂは、複数段にわたってツリー状に接続された複数のセレクタ部および乱数ビット生成部３２ａ，３２ｂ，…，３２ｒを有する。乱数ビット生成部３２ａ～３２ｒは、ツリー状に接続された複数のセレクタ部の段毎に設けられる。乱数ビット生成部３２ａ～３２ｒの各々は、０又は１の値をとる１ビット乱数を生成し、各段の選択回路に供給する。１ビット乱数は、入力されたフラグのペアのうちの何れか一方を選択するために用いられる。
図７は、選択部の回路構成例を示す図である。
選択部１４ｂは、複数段にわたってツリー状に接続された複数のセレクタ部および乱数ビット生成部３２ａ，３２ｂ，…，３２ｒを有する。乱数ビット生成部３２ａ～３２ｒは、ツリー状に接続された複数のセレクタ部の段毎に設けられる。乱数ビット生成部３２ａ～３２ｒの各々は、０又は１の値をとる１ビット乱数を生成し、各段の選択回路に供給する。１ビット乱数は、入力されたフラグのペアのうちの何れか一方を選択するために用いられる。
【0106】
初段のセレクタ部３１ａ１，３１ａ２，…，３１ａｐの各々には、判定部１４ａ１～１４ａＮの各々が出力する遷移可否を示すフラグの組が入力される。例えば、セレクタ部３１ａ１には、１番目の判定部１４ａ１が出力するフラグと、２番目の判定部１４ａ２が出力するフラグとのペアが入力される。また、セレクタ部３１ａ２には、３番目の判定部が出力するフラグと、４番目の判定部が出力するフラグとのペアが入力される。以降、同様にして、セレクタ部３１ａｐには、Ｎ－１番目の判定部が出力するフラグと、Ｎ番目の判定部１４ａＮが出力するフラグとのペアが入力される。このように、隣り合う判定部が出力するフラグのペアが初段のセレクタ部に入力される。初段のセレクタ部３１ａ１～３１ａｐの数は、Ｎ／２となる。以降、段を経る毎にセレクタ部の数は半分になる。
【0107】
セレクタ部３１ａ１～３１ａｐの各々は、入力されたフラグのペアと、乱数ビット生成部３２ａが出力する１ビット乱数に基づいて、入力されたフラグのペアのうちの一方を選択する。セレクタ部３１ａ１～３１ａｐの各々は、選択したフラグと選択したフラグに対応する１ビットの識別値とを含む状態信号を、２段目のセレクタ部３１ｂ１～３１ｂｑに出力する。例えば、セレクタ部３１ａ１が出力する状態信号と、セレクタ部３１ａ２が出力する状態信号とがセレクタ部３１ｂ１に入力される。同様に、セレクタ部３１ａ１～３１ａｐの隣り合うセレクタ部が出力する状態信号のペアが、２段目のセレクタ部に入力される。
【0108】
セレクタ部３１ｂ１～３１ｂｑの各々は、入力された状態信号のペアと、乱数ビット生成部３２ｂが出力する１ビット乱数に基づいて、入力された状態信号のうちの一方を選択する。セレクタ部３１ｂ１～３１ｂｑの各々は、選択した状態信号を、３段目のセレクタ部３１ｃ１，…に出力する。ここで、セレクタ部３１ｂ１～３１ｂｑは、選択した方の状態信号に含まれる状態信号について、何れの状態信号を選択したかを示すように１ビットを付加して更新し、選択した状態信号を出力する。
【0109】
３段目以降のセレクタ部においても同様の処理が行われ、各段のセレクタ部で１ビットずつ識別値のビット幅が増えていき、最終段のセレクタ部３１ｒから、選択部１４ｂの出力である状態信号が出力される。選択部１４ｂが出力する状態信号に含まれる識別値が、２進数で表されたインデックスに対応する。ただし、図７で示す回路構成例では、変数のインデックスを０始まりとする場合としている（１始まりのインデックスに対応させるには、選択部１４ｂが出力する識別値に１を加算した値をインデックスとすればよい）。
【0110】
例えば、図７では、セレクタ部３１ｂｑの回路構成例が示されている。２段目以降の他のセレクタ部も、セレクタ部３１ｂｑと同様の回路構成により実現される。
セレクタ部３１ｂｑの入力は、１つ目の状態信号（ｓｔａｔｕｓ＿１）と、２つ目の状態信号（ｓｔａｔｕｓ＿２）である。セレクタ部３１ｂ１の出力は、状態信号（ｓｔａｔｕｓ）である。セレクタ部３１ｂｑは、ＯＲ回路４１、ＮＡＮＤ回路４２およびセレクタ４３，４４を有する。
セレクタ部３１ｂｑの入力は、１つ目の状態信号（ｓｔａｔｕｓ＿１）と、２つ目の状態信号（ｓｔａｔｕｓ＿２）である。セレクタ部３１ｂ１の出力は、状態信号（ｓｔａｔｕｓ）である。セレクタ部３１ｂｑは、ＯＲ回路４１、ＮＡＮＤ回路４２およびセレクタ４３，４４を有する。
【0111】
ＯＲ回路４１には、状態信号（ｓｔａｔｕｓ＿１）に含まれるフラグ（ｆｌａｇ１）と、状態信号（ｓｔａｔｕｓ＿２）に含まれるフラグ（ｆｌａｇ２）とが入力される。例えば、状態信号（ｓｔａｔｕｓ＿１）が前段の２つのセレクタ部のうちの上位側（ｉｎｄｅｘの大きい側）の出力、状態信号（ｓｔａｔｕｓ＿２）が前段の２つのセレクタ部のうちの下位側（ｉｎｄｅｘの小さい側）の出力である。ＯＲ回路４１は、ｆｌａｇ１とｆｌａｇ２とのＯＲ演算結果（ｆｌａｇ）を出力する。
【0112】
ＮＡＮＤ回路４２には、ｆｌａｇ１とｆｌａｇ２とが入力される。ＮＡＮＤ回路４２は、ｆｌａｇ１とｆｌａｇ２とのＮＡＮＤ演算結果を、セレクタ４３の選択信号入力端子に出力する。
【0113】
セレクタ４３には、ｆｌａｇ１と１ビット乱数（ｒａｎｄ）とが入力される。セレクタ４３は、ＮＡＮＤ回路４２から入力されるＮＡＮＤ演算結果に基づいて、ｆｌａｇ１又はｒａｎｄの何れかを選択し、出力する。例えば、セレクタ４３は、ＮＡＮＤ回路４２のＮＡＮＤ演算結果が「１」の場合、ｆｌａｇ１を選択し、ＮＡＮＤ回路４２のＮＡＮＤ演算結果が「０」の場合、ｒａｎｄを選択する。
【0114】
セレクタ４４には、状態信号（ｓｔａｔｕｓ＿１）に含まれる識別値（ｉｎｄｅｘ１）と、状態信号（ｓｔａｔｕｓ＿２）に含まれる識別値（ｉｎｄｅｘ２）とが入力される。セレクタ４４の選択信号入力端子には、セレクタ４３の選択結果が入力される。セレクタ４４は、セレクタ４３の選択結果に基づいて、ｉｎｄｅｘ１又はｉｎｄｅｘ２の何れかを選択し、出力する。例えば、セレクタ４４は、セレクタ４３の選択結果が「１」の場合、ｉｎｄｅｘ１を選択し、セレクタ４３の選択結果が「０」の場合、ｉｎｄｅｘ２を選択する。
【0115】
ＯＲ回路４１およびセレクタ４３，４４の出力の組が、セレクタ部３１ｂｑが出力する状態信号（ｓｔａｔｕｓ）となる。
なお、初段のセレクタ部の入力は識別値を含まない。このため、初段のセレクタ部は、識別値（図中ではｉｎｄｅｘと表記）として、選択した方に対応するビット値（下位側の場合に「０」、上位側の場合に「１」）を追加して出力する回路となる。
なお、初段のセレクタ部の入力は識別値を含まない。このため、初段のセレクタ部は、識別値（図中ではｉｎｄｅｘと表記）として、選択した方に対応するビット値（下位側の場合に「０」、上位側の場合に「１」）を追加して出力する回路となる。
【0116】
このように、選択部１４ｂは、遷移可であるスピンビットから１つをトーナメント方式で選択する。トーナメントの各試合（すなわち、各選択回路での選択）では、勝った方（すなわち、選択された方）のエントリ番号（０又は１）をｉｎｄｅｘ ｗｏｒｄの上位ビットに付け加えていく。最終段のセレクタ部３１ｒが出力するインデックスが選ばれたスピンビットを示す。例えば、変数の数Ｎが１０２４個の場合、最終段のセレクタ部３１ｒが出力する状態信号は、遷移可否を示すフラグと、１０ビットで表されるインデックスとを含む。
【0117】
ただし、インデックスの出力方法は、上記のように選択部１４ｂで生成する方法以外の方法も考えられる。例えば、判定部１４ａ１～１４ａＮの各々から選択部１４ｂに各判定部に対応するインデックスを供給し、選択部１４ｂにより遷移可否を示すフラグとともに、当該フラグに対応するインデックスを選択してもよい。この場合、判定部１４ａ１～１４ａＮのそれぞれは、自身に対応するインデックスを格納するインデックスレジスタを更に有し、インデックスレジスタから選択部１４ｂにｉｎｄｅｘを供給する。
【0118】
図８は、変数毎のΔＥおよびΔＣの計算例を示す図である。
最適化装置１０は、全てのｉ（ｉ＝１～Ｎ）に関する（ｘ１，ｈ１），（ｘ２，ｈ２），…，（ｘＮ，ｈＮ）、および、識別コードｄｉに対して下記の手順を並列に実行する。
最適化装置１０は、全てのｉ（ｉ＝１～Ｎ）に関する（ｘ１，ｈ１），（ｘ２，ｈ２），…，（ｘＮ，ｈＮ）、および、識別コードｄｉに対して下記の手順を並列に実行する。
【0119】
（Ｓ１）最適化装置１０は、ｉ∈Ωであるか否かを判定する。ｉ∈Ωでない場合、ステップＳ２，Ｓ２ａに処理が進む。ｉ∈Ωである場合、ステップＳ３，Ｓ３ａに処理が進む。ｉ∈Ωであるか否かは、ｄｉ＝（ｄｉ１，ｄｉ２）により判定できる。すなわち、ｄｉ１およびｄｉ２の少なくとも一方が１（ｄｉ＝０１，１０，１１）ならば、ｉ∈Ωである。また、ｄｉ１およびｄｉ２の両方が０（ｄｉ＝００）ならば、ｉ∈Ωでない。
【0120】
（Ｓ２）最適化装置１０は、ΔＥｉを式（２）により計算される値（ΔＥｉ＝（２ｘｉ－１）ｈｉ）とする。
（Ｓ２ａ）最適化装置１０は、全てのインデックスｊにわたって、ΔＣｉｊ＝０とする。
（Ｓ２ａ）最適化装置１０は、全てのインデックスｊにわたって、ΔＣｉｊ＝０とする。
【0121】
（Ｓ３）最適化装置１０は、全てのインデックスｊにわたって、行列Ｗを保持するレジスタＲ１に記憶された係数ａｉｊに基づき、ΔＣｉｊを式（１６）により計算される値とする。ΔＣｉｊの計算に用いられる局所場ｈｉは、ｈｉを保持するレジスタから、不等式制約ｉに対するΔＣｉｊを計算する各ΔＣ計算部に供給される。
【0122】
（Ｓ３ａ）最適化装置１０は、ΔＥｉ＝Ｅｍａｘとする。ここで、Ｅｍａｘは、非常に大きな値であり、判定部１４ａｉにおける判定で常にＦｉ＝０となる値である。例えば、ＥｍａｘをΔＥｉとしてとれる値の最大値とすることが考えられる。
【0123】
また、最適化装置１０は、総エネルギーＥｔｏｔ（ｘ）のうち、ＱＵＢＯ項に相当するエネルギーＥ（ｘ）を次の回路によってモニタリングする。
図９は、関数Ｅ（ｘ）に対応するエネルギーを計算する回路構成例を示す図である。
図９は、関数Ｅ（ｘ）に対応するエネルギーを計算する回路構成例を示す図である。
【0124】
最適化装置１０は、更に、エネルギー計算部１７および制御部１８を有する。
エネルギー計算部１７は、総エネルギーＥｔｏｔ（ｘ）のうち、ＱＵＢＯ項に相当するエネルギーＥ（ｘ）を計算する。エネルギー計算部１７は、更新制御部１４により選択された変化対象の状態変数の識別情報に基づいて、エネルギー変化計算部１２により出力される複数のエネルギーの変化値のうちの何れかの変化値を選択する。エネルギー計算部１７は、選択した変化値を累積することで状態保持部１１に保持される複数の状態変数の値に対応するエネルギーの値Ｅ（ｘ）を計算する。具体的には、エネルギー計算部１７は、選択部１７ａおよび累算器１７ｂを有する。
エネルギー計算部１７は、総エネルギーＥｔｏｔ（ｘ）のうち、ＱＵＢＯ項に相当するエネルギーＥ（ｘ）を計算する。エネルギー計算部１７は、更新制御部１４により選択された変化対象の状態変数の識別情報に基づいて、エネルギー変化計算部１２により出力される複数のエネルギーの変化値のうちの何れかの変化値を選択する。エネルギー計算部１７は、選択した変化値を累積することで状態保持部１１に保持される複数の状態変数の値に対応するエネルギーの値Ｅ（ｘ）を計算する。具体的には、エネルギー計算部１７は、選択部１７ａおよび累算器１７ｂを有する。
【0125】
選択部１７ａは、反転させる状態変数のインデックスｉの選択結果を選択部１４ｂから受け付ける。選択部１７ａは、エネルギー変化計算部１２により出力されるΔＥｉ（ｉ＝１，…，ｎ）のうち、選択部１４ｂから供給されるインデックスｉに対応するΔＥｉを選択し、累算器１７ｂに出力する。
【0126】
累算器１７ｂは、エネルギーＥ（ｘ）の値を保持し、選択部１７ａから供給されるΔＥｉを、エネルギーＥ（ｘ）の値に累算することで、現在の状態変数の値に対応するエネルギーＥ（ｘ）を更新し、当該エネルギーＥ（ｘ）を出力する。なお、エネルギーＥ（ｘ）の初期値は、状態ベクトルの初期値に対して、演算開始時に累算器１７ｂに与えられる。
【0127】
制御部１８は、最適化装置１０における探索部Ｄ１による演算を制御する。ここで、探索部Ｄ１は、状態保持部１１、エネルギー変化計算部１２、ペナルティ加算部１３、更新制御部１４、制約超過量計算部１５、制約識別符号入力部１６、選択部１７ａおよび累算器１７ｂを含む回路ブロックである。制御部１８は、探索部Ｄ１における探索回数（状態変数の反転回数）や温度設定などを制御する。制御部１８は、電子回路により実現される。制御部１８は、制御回路とも呼べる。
【0128】
最適化装置１０は、図９の回路構成により、ＱＵＢＯ項に相当するエネルギーＥ（ｘ）をモニタリングできる。図９の回路構成によってエネルギーＥ（ｘ）をモニタリングすることで、不等式制約における不等式制約項における係数λｉの値が比較的大きい場合でも、エネルギーＥ（ｘ）の計算の桁あふれがないという利点がある。
【0129】
次に、最適化装置１０による演算の手順を説明する。
図１０は、最適化装置の演算例を示すフローチャートである。
（Ｓ１０）制御部１８は、探索部Ｄ１の初期化を行う。例えば、制御部１８は、温度や、温度毎の繰り返し回数Ｃ１および温度の更新回数Ｃ２などの外部からの設定を受け付け、探索部Ｄ１に初期温度を設定する。また、制御部１８は、制約識別符号を制約識別符号入力部１６に設定する。制御部１８は、初期状態を状態保持部１１に設定する。このとき、制御部１８は、制約識別符号により不等式制約用変数として指定された変数の値を０に設定する。更に、制御部１８は、重みＷ、係数λ、各変数に対応する局所場の初期値を所定のレジスタに設定するとともに、初期状態に対応するエネルギーＥ（ｘ）を累算器１７ｂに設定する。制御部１８は、初期化が完了すると探索部Ｄ１による演算を開始させる。
図１０は、最適化装置の演算例を示すフローチャートである。
（Ｓ１０）制御部１８は、探索部Ｄ１の初期化を行う。例えば、制御部１８は、温度や、温度毎の繰り返し回数Ｃ１および温度の更新回数Ｃ２などの外部からの設定を受け付け、探索部Ｄ１に初期温度を設定する。また、制御部１８は、制約識別符号を制約識別符号入力部１６に設定する。制御部１８は、初期状態を状態保持部１１に設定する。このとき、制御部１８は、制約識別符号により不等式制約用変数として指定された変数の値を０に設定する。更に、制御部１８は、重みＷ、係数λ、各変数に対応する局所場の初期値を所定のレジスタに設定するとともに、初期状態に対応するエネルギーＥ（ｘ）を累算器１７ｂに設定する。制御部１８は、初期化が完了すると探索部Ｄ１による演算を開始させる。
【0130】
下記のステップＳ１１，Ｓ１２についてブロックＢｊに着目して説明するが、ステップＳ１１，Ｓ１２は、ブロックＢ１～ＢＮにおいて、並列に実行される。
（Ｓ１１）ブロックＢｊは、ｈｊとｘｊとを読み出す。前述のようにｈｊは、所定のレジスタに保持される。ｘｊは、状態保持部１１に保持される。
（Ｓ１１）ブロックＢｊは、ｈｊとｘｊとを読み出す。前述のようにｈｊは、所定のレジスタに保持される。ｘｊは、状態保持部１１に保持される。
【0131】
（Ｓ１２）ブロックＢｊは、図８の手順に基づいてΔＣｉｊおよびΔＥｊの計算を行うとともに、ΔＣｉｊに基づいてΔＣｊを計算し、ΔＥｔｏｔ ｊ＝ΔＥｊ＋ΔＣｊを求める。ブロックＢｊは、－（ΔＥｊ＋ΔＣｊ）に対してＥｏｆｆを加算し、更に、熱ノイズ－Ｔ・ｌｎ（ｕ）を加算した評価値と閾値との判定に基づいて、フラグＦｊを出力する。
【0132】
（Ｓ１３）選択部１４ｂは、ブロックＢ１～ＢＮにより出力されるフラグＦ１～ＦＮに基づいて、反転対象のインデックスｉを選択し、選択したインデックスｉを状態保持部１１に供給する。
【0133】
（Ｓ１４）状態保持部１１は、選択部１４ｂから供給されるインデックスｉに対応する状態変数ｘｉの値を反転させる（変化させる）ことで、状態を更新するとともに、所定のレジスタに保持された各変数に対応する局所場ｈを更新する（式（４）に基づく更新）。なお、本ステップの説明では、変化対象のインデックスをｉと表記している。当該更新により、式（１３）に基づく不等式制約用変数に対する局所場の値も更新される。
【0134】
（Ｓ１５）制御部１８は、現温度において、規定回数Ｃ１だけ状態変数の更新を行ったか否かを判定する。規定回数Ｃ１だけ更新を行った場合、ステップＳ１６に処理が進む。規定回数Ｃ１だけ更新を行っていない場合、ステップＳ１１に処理が進む。
【0135】
（Ｓ１６）制御部１８は、温度の更新を行う。これにより、探索部Ｄ１に設定される温度Ｔが下げられる。制御部１８は、規定回数Ｃ１の更新のたびに、予め定められた方法により、探索部Ｄ１に設定する温度Ｔを下げる。
【0136】
（Ｓ１７）制御部１８は、規定回数Ｃ２だけ温度を更新したか否かを判定する。規定回数Ｃ２だけ温度を更新した場合、ステップＳ１８に処理が進む。規定回数Ｃ２だけ温度を更新していない場合、ステップＳ１１に処理が進む。
【0137】
（Ｓ１８）制御部１８は、累算器１７ｂにより出力された最低エネルギー（エネルギーＥ（ｘ）の最小値）に対応する状態を出力する。例えば、制御部１８は、探索部Ｄ１が最終的に到達した状態を、最低エネルギーに対応する状態として出力してもよい。あるいは、制御部１８は、上記の探索の過程で、探索部Ｄ１が到達した最低エネルギーに対応する状態を出力してもよい。
【0138】
ここで、ステップＳ１２において、ペナルティ加算部１３は、次のようにペナルティ値をエネルギーの変化値に加算することで、総エネルギーの変化値ΔＥｔｏｔ ｊ＝ΔＥｊ＋ΔＣｊ（ｊ∈Ψ）を算出する。
【0139】
例えば、ペナルティ加算部１３は、接続係数ａｉｊ（ｉ∈Ω）と閾値（ＣｕｉおよびＣｌｉ）に基づいて算出される制約超過量の変化値ΔＣｉｊに不等式制約ｉの重要度を示す重要度係数λｉを乗じることでペナルティ値ΔＣｊを算出し、ΔＥｊに加算する。ペナルティ加算部１３は、不等式制約が複数ある場合は、不等式制約毎のλｉΔＣｉｊをｉについて和をとって、ΔＣｊを求める。重要度係数により、ＱＵＢＯ項に対する不等式制約の重みや、不等式制約が複数ある場合の各不等式制約の重みを調整可能にして、求解を行える。
【0140】
また、不等式制約が上限制約値Ｃｕｉ（ｉ∈Ω）に対する不等式制約である場合、ペナルティ加算部１３は、複数のエネルギーの変化値ΔＥｊのそれぞれに対して、第１の接続係数ａｉｊと第１の閾値Ｃｕｉとに基づいて算出される、第１の閾値Ｃｕｉに対する正方向の第１の超過量の変化値に応じた第１のペナルティ値ΔＣｊを加算する。
【0141】
また、不等式制約が下限制約値Ｃｌｉに対する不等式制約である場合、ペナルティ加算部１３は、複数のエネルギーの変化値ΔＥｊのそれぞれに対して、第２の接続係数ａｉｊと第２の閾値Ｃｌｉとに基づいて算出される、第２の閾値Ｃｌｉに対する負方向の第２の超過量の変化値に応じた第２のペナルティ値ΔＣｊを加算する。
【0142】
また、不等式制約が上限制約値Ｃｕｉおよび下限制約値Ｃｌｉに対する不等式制約である場合、ペナルティ加算部１３は、複数のエネルギーの変化値ΔＥｊのそれぞれに対して、第３の接続係数ａｉｊと第３の閾値Ｃｕｉとに基づいて算出される、第３の閾値Ｃｕｉに対する正方向の第３の超過量の変化値、および、第３の接続係数ａｉｊと第３の閾値Ｃｕｉよりも小さい第４の閾値Ｃｌｉとに基づいて算出される、第４の閾値Ｃｌｉに対する負方向の第４の超過量の変化値のうち、大きい方の変化値に応じた第３のペナルティ値ΔＣｊを加算する。
【0143】
また、所定の制約値に対する等式制約の場合、ペナルティ加算部１３は、複数のエネルギーの変化値ΔＥｊのそれぞれに対して、第４の接続係数ａｉｊと第５の閾値Ｃｕｉ＝Ｃｌｉとに基づいて算出される、第５の閾値Ｃｕｉ＝Ｃｌｉに対する正方向の第５の超過量の変化値、および、第５の閾値Ｃｕｉ＝Ｃｌｉに対する負方向の第６の超過量の変化値のうち、大きい方の変化値に応じた第４のペナルティ値ΔＣｊを加算する。
【0144】
図１１は、求解結果の比較の例を示す図である。
図１１（Ａ）は、不等式制約を含む最適化問題に対する最適化装置１０を用いた求解結果の例を示す。図１１（Ｂ）は、ペナルティ加算部１３を有していない既存の最適化装置を用いた、不等式制約を状態変数の二次形式で表す場合（二乗制約の場合）の最適化問題に対する求解結果の例を示す。
図１１（Ａ）は、不等式制約を含む最適化問題に対する最適化装置１０を用いた求解結果の例を示す。図１１（Ｂ）は、ペナルティ加算部１３を有していない既存の最適化装置を用いた、不等式制約を状態変数の二次形式で表す場合（二乗制約の場合）の最適化問題に対する求解結果の例を示す。
【0145】
ここで、図１１では、最適化問題の一例として、ナップザック問題を考える。（ペナルティ加算部１３を有していない）既存の最適化装置を用いてナップザック問題を解く場合、次のように問題を定式化し得る。
【0146】
ここで考えるナップザック問題は、Ｎ個の品物（インデックスｉの品物の重さはｗｉ、値段はｖｉ）を、積載容量がＣのナップザックに入れ、ナップザック中の品物の値段の合計（総価値）を最大にする問題である。ナップザックに品物ｉを入れる場合に状態変数ｘｉ＝１とし、入れない場合に状態変数ｘｉ＝０とする。
【0147】
この場合、ナップザック中の品物の総価値Ｖは、式（１８）で表される。
【0148】
【数18】
【0149】
また、ナップザックにおける容量制約は、式（１９）で表される。
【0150】
【数19】
【0151】
式（１９）の容量制約は、スラック変数Ｓ（Ｓ≧０）を用いて、式（２０）のように表される。
【0152】
【数20】
【0153】
式（２０）に基づいて、ＱＵＢＯの目的関数Ｅを式（２１）で表せる。
【0154】
【数21】
【0155】
Ｓを式（２２）のようにバイナリ展開し、状態変数ｘｉ，ｙｉのＱＵＢＯ形式として解くことが考えられる。状態変数ｘｉ，ｙｉのＱＵＢＯ形式は式（２３）で表される。
【0156】
【数22】
【0157】
【数23】
【0158】
スラック変数Ｓは、最適化で比較的小さな値となるため、Ｓ＝０としても最適解が取得され得る。よって、式（２４）を得る。
【0159】
【数24】
【0160】
そこで、図１１の比較では、次のようなテストデータを作成して実験を行った。
ナップザック容量Ｃ＝４９５２９９１である。品物毎の重さＷ＝Ｄ１～ＤＮは、それぞれ区間［０，１００］の一様乱数に１０００００を加算した値である。品物毎の値段Ｐ＝Ｐ１～ＰＮは、それぞれ区間［１，１０００］の一様乱数として発生させた値である。
ナップザック容量Ｃ＝４９５２９９１である。品物毎の重さＷ＝Ｄ１～ＤＮは、それぞれ区間［０，１００］の一様乱数に１０００００を加算した値である。品物毎の値段Ｐ＝Ｐ１～ＰＮは、それぞれ区間［１，１０００］の一様乱数として発生させた値である。
【0161】
また、実験条件を次の通りとした。第１に、不等式制約に対する重要度係数λをλ＝１０＾ｚとして、ｚ＝－１０から１０まで、１刻みで網羅する。第２に、温度Ｔ＝５，１０，５０，１００，５００（初期温度）として、それぞれメトロポリス基準に従ったモンテカルロ法による演算を行う。第３に、乱数ｕの初期値を１００個用意し、この中で１つでも制約を満たす解が得られる初期値があれば、「解が得られる」ものとする。
【0162】
図１１の比較における最適化装置１０に対する不等式制約（絶対値制約）を含むエネルギー関数Ｅ（ｘ）は、式（２５）で表される。
【0163】
【数25】
【0164】
また、図１１の比較における二乗制約の場合のエネルギー関数Ｅ（ｘ）は、式（２６）で表される。
【0165】
【数26】
【0166】
ただし、式（２６）では、スラック変数は小さな値が最適であるため、スラック変数を０と置いた。
図１１（Ａ）および図１１（Ｂ）の各グラフの横軸はλ＝１０＾ｚにおけるｚである。図１１（Ａ）および図１１（Ｂ）の各グラフの縦軸は価値Ｐｉの合計値（最大）である。なお、制約を満たさない場合には、Ｐｉ＝０である。
図１１（Ａ）および図１１（Ｂ）の各グラフの横軸はλ＝１０＾ｚにおけるｚである。図１１（Ａ）および図１１（Ｂ）の各グラフの縦軸は価値Ｐｉの合計値（最大）である。なお、制約を満たさない場合には、Ｐｉ＝０である。
【0167】
上記の実験条件では、絶対値制約では何れの温度でもλが１０＾（－１）以上であれば、解を得られるのに対して、二乗制約では全てＰｉ＝０となり、制約を満たす解を得られなかった。そこで、λの範囲を広げ、温度ＴをＴ＝１０＾１０，１０＾１５，１０＾２０というように大きく設定したところ、二乗制約では限られた範囲のみで解が得られ、二乗制約では探索に比較的大きなコストがかかることが分かる。
【0168】
係数λが比較的小さい場合に解が得られない理由は、総価値に比べて不等式制約項の値が小さくなることで、制約が軽視され、制約を満たす解を得られる可能性が低下するためと考えられる。
【0169】
図１１（Ａ）において、系列ｐ１は、Ｔ＝１０＾１０の場合である。系列ｐ２は、Ｔ＝１０＾１５の場合である。系列ｐ３は、Ｔ＝１０＾２０の場合である。系列ｐ４は、Ｔ＝５の場合である。系列ｐ５は、Ｔ＝１０の場合である。系列ｐ６は、Ｔ＝５０の場合である。系列ｐ７は、Ｔ＝１００の場合である。系列ｐ８は、Ｔ＝５００の場合である。
【0170】
図１１（Ｂ）において、系列ｑ１は、Ｔ＝１０＾１０の場合である。系列ｑ２は、Ｔ＝１０＾１５の場合である。系列ｑ３は、Ｔ＝１０＾２０の場合である。図１１（Ｂ）の二乗制約の例では、温度Ｔ＝５，１０，５０，１００，５００に対応する系列は、解が得られなかったため、図示を省略している。
【0171】
二乗制約の場合において、このような結果となった原因は、二乗制約のエネルギー関数の式にあると考えられる。すなわち、ナップザックに入っている品物の重さの合計がナップザックの容量付近である場合、品物を取り除いても式（２６）に基づくエネルギーが増加してしまい、探索ができなくなるためである。特に、二乗制約の場合は、二次の項により、一次の項で表される絶対値制約の場合よりも、式（２６）のエネルギー関数Ｅ（ｘ）の値に対する不等式制約項の寄与が過大になり得る。
【0172】
ここで、最適化問題では様々な制約条件が入るため、二次形式だけでは問題を定式化するのか困難なことがある。上記の例以外にも、例えば、配車配送問題などの積載量上限、工場を建設するための予算の上限など、必要なリソース量が上限を超えない、あるいは下限を下回らないなど不等式で表される制約が頻出し得る。
【0173】
これらの不等式制約は、二次形式で表して既存の最適化装置で解くことが難しいことがあり、また、不等式制約を厳密に満たすように探索すると最適解に到達し難いという問題がある。これは、不等式制約を破るような中間状態を経由しないと最適解に到達できないことが多いためである。
【0174】
そこで、ペナルティ加算部１３を有する最適化装置１０を用いることで、不等式制約を二次形式に変形しなくても、元の一次式のままで容易に扱えるようになる。
図１１で例示されるように、最適化装置１０によれば、既存の最適化装置よりも、係数λや温度Ｔについて、広い範囲で解を得ることができる。このため、二乗制約の場合よりも、不等式制約を容易に扱うことができ、演算対象にできる問題の種類を拡張できる。すなわち、最適化装置１０の応用範囲を広げることができ、より多くの問題の求解に、最適化装置１０を利用可能になる。
図１１で例示されるように、最適化装置１０によれば、既存の最適化装置よりも、係数λや温度Ｔについて、広い範囲で解を得ることができる。このため、二乗制約の場合よりも、不等式制約を容易に扱うことができ、演算対象にできる問題の種類を拡張できる。すなわち、最適化装置１０の応用範囲を広げることができ、より多くの問題の求解に、最適化装置１０を利用可能になる。
【0175】
上記の例では、ｋ∈Ωに対して、ｘｋ＝０かつΔＥｋ＝Ｅｍａｘとすることで、変数ｘｋによるエネルギーＥ（ｘ）への影響を抑制するものとしたが、係数ａｋｊとして、非対称接続係数を用いて、変数ｘｋによるエネルギーＥ（ｘ）への影響を抑制してもよい。非対称接続係数は、前述のように、ａｋｊ≠０、ａｊｋ＝０と表される。なお、非対称接続係数は、変数間の重み値の記号Ｗを用いてＷｋｊ≠０、Ｗｊｋ＝０と表されてもよい。
【0176】
すなわち、状態保持部１１は、状態変数と不等式制約に対応する不等式制約用変数との組に対する接続係数ａｉｊを保持してもよい。接続係数ａｉｊは、不等式制約用変数ｘｋ（ｋ∈Ω）の局所場ｈｋにおける状態変数ｘｊの重み（ａｋｊ）を０以外の値とし、状態変数ｘｊの局所場ｈｊにおける不等式制約用変数ｘｋの重み（ａｊｋ）を０とする非対称接続係数であってもよい。
【0177】
非対称接続係数を用いる場合、選択部１４ｂにより不等式制約用変数が反転対象として選択されたとしても、エネルギーＥ（ｘ）に対する変数ｘｋの反転の影響を抑制できる。
次に、第１の実施の形態の他の構成例を説明する。まず、最適化装置１０を含む最適化システムの例を説明する。
次に、第１の実施の形態の他の構成例を説明する。まず、最適化装置１０を含む最適化システムの例を説明する。
【0178】
図１２は、最適化システムの例を示す図である。
最適化システム１００は、最適化装置１０および情報処理装置１１０を有する。
情報処理装置１１０は、ＣＰＵ（Central Processing Unit）１１１、メモリ１１２およびバス１１３を有する。
最適化システム１００は、最適化装置１０および情報処理装置１１０を有する。
情報処理装置１１０は、ＣＰＵ（Central Processing Unit）１１１、メモリ１１２およびバス１１３を有する。
【0179】
ＣＰＵ１１１は、情報処理装置１１０を制御するプロセッサである。ＣＰＵ１１１は、メモリ１１２に記憶されたプログラムを実行する。メモリ１１２は、例えば、ＲＡＭ（Random Access Memory）であり、ＣＰＵ１１１により実行されるプログラムやＣＰＵ１１１により処理されるデータを記憶する。バス１１３は、ＰＣＩｅ（Peripheral Component Interconnect Express）などの内部バスである。バス１１３には、ＣＰＵ１１１、メモリ１１２および最適化装置１０が接続される。
【0180】
ＣＰＵ１１１は、バス１１３を介して、最適化装置１０に対する温度Ｔの設定や規定回数Ｃ１，Ｃ２に関する情報の入力を行う。また、ＣＰＵ１１１は、制約識別符号（識別コード、上限値、下限値を含む）や不等式制約毎の係数λを、最適化装置１０に入力する。制御部１８は、こうして入力された情報に基づいて、探索部Ｄ１に対する設定を行う。例えば、ユーザは、情報処理装置１１０に接続された入力部（図示を省略している）を操作して、制約識別符号や係数λを、情報処理装置１１０に設定することができる。
【0181】
このように、最適化装置１０は、情報処理装置１１０に接続されて、不等式制約を含む組合せ最適化問題に対する演算をハードウェアにより高速に実行するアクセラレータとして利用される。
【0182】
次に、最適化装置の他の回路構成例を説明する。
図１３は、最適化装置の他の回路構成例を示す図である。
最適化装置１０ａは、探索部Ｄ２および制御部１８ａを有する。探索部Ｄ２は、状態保持部１１ａ、レジスタ１１ｂ１～１１ｂＮ、ｈ計算部１１ｃ１～１１ｃＮ、ΔＥ計算部１２ａ１～１２ａＮ、ΔＣ加算部１３ｄ１～１３ｄＮ、判定部１４ａ１～１４ａＮおよび選択部１４ｂを有する。
図１３は、最適化装置の他の回路構成例を示す図である。
最適化装置１０ａは、探索部Ｄ２および制御部１８ａを有する。探索部Ｄ２は、状態保持部１１ａ、レジスタ１１ｂ１～１１ｂＮ、ｈ計算部１１ｃ１～１１ｃＮ、ΔＥ計算部１２ａ１～１２ａＮ、ΔＣ加算部１３ｄ１～１３ｄＮ、判定部１４ａ１～１４ａＮおよび選択部１４ｂを有する。
【0183】
ここで、図１３では、ｈ計算部１１ｃ１～１１ｃＮ、ΔＥ計算部１２ａ１～１２ａＮおよびΔＣ加算部１３ｄ１～１３ｄＮについて、インデックスｉの状態変数に対応することが分かり易いように「ｈｉ」計算部などのように添え字ｉを付して名称を表記している。また、状態保持部１１ａ、レジスタ１１ｂ１～１１ｂＮおよびｈ計算部１１ｃ１～１１ｃＮは、前述の状態保持部１１に対応する。
【0184】
状態保持部１１ａは、状態ベクトルを保持する。状態保持部１１ａは、選択部１４ｂから供給されるｕｐｄａｔｅ信号に基づいて、状態ベクトルのうちの１つの状態変数の値を反転させることで、状態ベクトルを更新する。ｕｐｄａｔｅ信号は、反転させる状態変数のインデックスを示す。
【0185】
レジスタ１１ｂ１、ｈ計算部１１ｃ１、ΔＥ計算部１２ａ１、ΔＣ加算部１３ｄ１および判定部１４ａ１が１番目の状態変数に関する演算を行う。また、レジスタ１１ｂ２、ｈ計算部１１ｃ２、ΔＥ計算部１２ａ２、ΔＣ加算部１３ｄ２および判定部１４ａ２が２番目のスピンビットに関する演算を行う。同様に、「１１ｂ１」、「１２ａ１」などの符号の末尾の数値ｉがインデックスｉの状態変数に対応する演算を行うことを示す。すなわち、探索部Ｄ２は、レジスタ、ｈ計算部、ΔＥ計算部、ΔＣ加算部および判定部のセット（１スピンビットに関する演算を行う演算処理回路の一単位であり、「ニューロン」と呼ばれることもある）を、Ｎ個有する。Ｎ個のセットが並列に、各セットに対応する状態変数に関する演算を行う。
【0186】
なお、あるニューロンにおけるΔＣ加算部には、他のニューロンにおけるレジスタに保持される重み値およびｈ計算部により保持される局所場が入力されるが、当該ΔＣ加算部と他のニューロンのレジスタおよびｈ計算部との信号線の図示を省略している。
【0187】
例えば、ΔＣ加算部１３ｄ１の例でいえば、Ｗ１１，Ｗ２１，…，ＷＮ１のうち、不等式制約用変数ｘｋ（ｋ∈Ω）に対応する１以上の重み値（Ｗｋ１＝ａｋ１）がΔＣ加算部１３ｄ１に供給される。また、ｈ計算部１１ｃ１～１１ｃＮのうち、不等式制約用変数ｘｋに対応する１以上の局所場（ｈｋ）がΔＣ加算部１３ｄ１に供給される。なお、状態変数ｘｉ（ｉ∈Ψ）に対応するＷｉ１やｈｉがΔＣ加算部１３ｄ１に供給されても、図８で示した識別コードに基づく制御によりΔＣ加算部１３ｄ１は、ΔＣｉ１＝０とする。
【0188】
また、ΔＣ加算部１３ｄ２の例でいえば、Ｗ１２，Ｗ２２，…，ＷＮ２のうち、不等式制約用変数ｘｋ（ｋ∈Ω）に対応する１以上の重み値（Ｗｋ２＝ａｋ２）がΔＣ加算部１３ｄ２に供給される。また、ｈ計算部１１ｃ１～１１ｃＮのうち、不等式制約用変数ｘｋに対応する１以上の局所場（ｈｋ）がΔＣ加算部１３ｄ２に供給される。なお、状態変数ｘｉ（ｉ∈Ψ）に対応するＷｉ２やｈｉがΔＣ加算部１３ｄ２に供給されても、識別コードに基づく制御によりΔＣ加算部１３ｄ２は、ΔＣｉ２＝０とする。
【0189】
以下では、主に、レジスタ１１ｂ１、ｈ計算部１１ｃ１、ΔＥ計算部１２ａ１、ΔＣ加算部１３ｄ１および判定部１４ａ１を例示して説明する。同名の構成であるレジスタ１１ｂ２～１１ｂＮ、ｈ計算部１１ｃ２～１１ｃＮ、ΔＥ計算部１２ａ２～１２ａＮ、ΔＣ加算部１３ｄ２～１３ｄＮおよび判定部１４ａ２～１４ａＮも同様の機能である。
【0190】
レジスタ１１ｂ１は、状態変数ｘ１とそれ以外の他の状態変数との間の重み値Ｗ１ｊ（ｊ＝１～Ｎ）を記憶する記憶部である。ここで、状態変数の数Ｎに対して、重み値の総数はＮ＾２である。レジスタ１１ｂ１には、Ｎ個の重み値が格納される。
【0191】
レジスタ１１ｂ１は、インデックス１に対して、Ｎ個の重み値Ｗ１１，Ｗ１２，…，Ｗ１Ｎを記憶する。なお、Ｗｉｉ＝Ｗ１１＝０である。レジスタ１１ａ１は、選択部１４ｂにより供給されるｉｎｄｅｘ＝ｊに対応する重み値Ｗ１ｊをｈ計算部１１ｃ１に出力する。
【0192】
ｈ計算部１１ｃ１は、レジスタ１１ｂ１から供給される重み係数Ｗ１ｊを用いて、式（３）、（４）に基づく局所場ｈ１を計算する。ｈ計算部１１ｃ１は、前回計算された局所場ｈ１を保持するレジスタを有し、ｉｎｄｅｘ＝ｊで示される状態変数の反転方向に応じたδｈ１
（ｊ）を、ｈ１に積算することで、当該レジスタに格納されるｈ１を更新する。ｈ１の初期値は、問題に応じて、ｈ計算部１１ｃ１のレジスタに予め設定される。また、ｂ１の値も、問題に応じて、ｈ計算部１１ｃ１のレジスタに予め設定される。ｈ計算部１１ｃ１は、計算した局所場ｈ１をΔＥ計算部１２ａ１に出力する。
【0193】
ΔＥ計算部１２ａ１は、ｈ計算部１１ｃ１から供給される局所場ｈ１に基づいて、式（２）によりΔＥ１を計算し、ΔＣ加算部１３ｄ１に出力する。ここで、ΔＥ計算部１２ａ１は、状態変数ｘ１の反転方向に応じてｈ１の符号を決めることでΔＥ１を計算する。
【0194】
ΔＣ加算部１３ｄ１は、式（１５）に基づいて不等式制約に対する超過量の変化値に応じたペナルティ値ΔＣ１を計算し、ΔＥ計算部１２ａ１から供給されるΔＥ１に対して、ΔＣ１を加算する。ΔＣ加算部１３ｄ１は、ΔＥ１＋ΔＣ１を判定部１４ａ１に出力する。
【0195】
判定部１４ａ１は、ΔＣ加算部１３ｄ１から供給されるΔＥ１＋ΔＣ１に基づいて、反転可否を示すフラグＦ１を選択部１４ｂに出力する。
選択部１４ｂは、判定部１４ａ１～１４ａＮから供給されるフラグＦ１～ＦＮに基づいて、反転対象の状態変数のインデックスｊを選択し、状態保持部１１ａおよびレジスタ１１ｂ１～１１ｂＮに供給する。なお、選択部１４ｂは、判定部１４ａ１～１４ａＮに対して、反転対象の状態変数のインデックスｊを選択したか否かを示すフラグを供給することで、オフセット値の生成を促すが、当該フラグを供給する信号線の図示を省略している。
選択部１４ｂは、判定部１４ａ１～１４ａＮから供給されるフラグＦ１～ＦＮに基づいて、反転対象の状態変数のインデックスｊを選択し、状態保持部１１ａおよびレジスタ１１ｂ１～１１ｂＮに供給する。なお、選択部１４ｂは、判定部１４ａ１～１４ａＮに対して、反転対象の状態変数のインデックスｊを選択したか否かを示すフラグを供給することで、オフセット値の生成を促すが、当該フラグを供給する信号線の図示を省略している。
【0196】
制御部１８ａは、探索部Ｄ２における判定部１４ａ１～１４ａＮに対する温度Ｔの設定、状態変数の更新回数の制御、ΔＥ計算部１２ａ１～１２ａＮおよびΔＣ加算部１３ｄ１～１３ｄＮに対する制約識別符号の入力を行う。ΔＥ計算部１２ａ１～１２ａＮのΔＥ計算およびΔＣ加算部１３ｄ１～１３ｄＮのΔＣ計算は、図８で示されるように、制約識別符号により制御される。
【0197】
最適化装置１０ａも、最適化装置１０と同様に、不等式制約を含む最適化問題の求解を効率的に実行できる。最適化装置１０ａを用いることで、不等式制約を容易に扱うことができ、演算対象にできる問題の種類を拡張できる。すなわち、最適化装置１０ａの応用範囲を広げることができ、より多くの問題の求解に、最適化装置１０ａを利用可能になる。
【0198】
［第２の実施の形態］
次に、第２の実施の形態を説明する。前述の第１の実施の形態と相違する事項を主に説明し、共通する事項の説明を省略する。
次に、第２の実施の形態を説明する。前述の第１の実施の形態と相違する事項を主に説明し、共通する事項の説明を省略する。
【0199】
第１の実施の形態では、最適化装置１０または最適化装置１０ａが単一の探索部を有する構成を例示した。一方、複数の探索部を有する最適化装置において、レプリカ交換法を用いて、最適化問題の求解を行う構成も考えられる。
【0200】
図１４は、第２の実施の形態の最適化装置の例を示す図である。
最適化装置５０は、探索部５１ａ１～５１ａＭ（Ｍ（Ｍは２以上の整数）個の探索部）および交換制御部５２を有する。
最適化装置５０は、探索部５１ａ１～５１ａＭ（Ｍ（Ｍは２以上の整数）個の探索部）および交換制御部５２を有する。
【0201】
探索部５１ａ１～５１ａＭは、第１の実施の形態の探索部Ｄ１または探索部Ｄ２と同様の回路構成を有する。図１４では、探索部５１ａｊの例が示されている。
探索部５１ａｊは、状態保持部１１、エネルギー変化計算部１２、ペナルティ加算部１３、更新制御部１４、制約超過量計算部１５および制約識別符号入力部１６を有する。状態保持部１１、エネルギー変化計算部１２、ペナルティ加算部１３、更新制御部１４、制約超過量計算部１５および制約識別符号入力部１６のそれぞれの機能は、第１の実施の形態における同名の要素と同様である。また、探索部５１ａｊは、エネルギー計算部１７を更に有する（図示を省略している）。
探索部５１ａｊは、状態保持部１１、エネルギー変化計算部１２、ペナルティ加算部１３、更新制御部１４、制約超過量計算部１５および制約識別符号入力部１６を有する。状態保持部１１、エネルギー変化計算部１２、ペナルティ加算部１３、更新制御部１４、制約超過量計算部１５および制約識別符号入力部１６のそれぞれの機能は、第１の実施の形態における同名の要素と同様である。また、探索部５１ａｊは、エネルギー計算部１７を更に有する（図示を省略している）。
【0202】
交換制御部５２は、探索部５１ａ１～５１ａＭのそれぞれに互いに異なる温度値を設定し、探索部５１ａ１～５１ａＭを用いたレプリカ交換法による基底状態探索を制御する。交換制御部５２は、探索部５１ａ１～５１ａＭのそれぞれによる基底状態探索の繰り返し回数到達または一定時間経過後に、探索部間で温度値または各探索部に保持される状態変数の値を入れ替える。例えば、交換制御部５２は、温度が隣接する２つの探索部のペアで、所定の確率（交換確率）により温度を交換することが考えられる。交換確率としては、例えば、メトロポリス法に基づく確率が用いられる。あるいは、交換制御部５２は、温度の代わりに、温度が隣接する２つの探索部のペアで、当該２つの探索部が保持する状態（複数の状態変数の値）や各状態変数に対応する局所場を交換してもよい。
【0203】
このように、最適化装置５０によれば、レプリカ交換法を用いる場合でも、不等式制約を含む最適化問題の求解を効率的に実行できる。最適化装置５０によれば、不等式制約を容易に扱うことができ、演算対象にできる問題の種類を拡張できる。最適化装置５０の応用範囲を広げることができ、より多くの問題の求解に、最適化装置５０を利用可能になる。
【符号の説明】
【0204】
１０ 最適化装置
１１ 状態保持部
１２ エネルギー変化計算部
１３ ペナルティ加算部
１４ 更新制御部
１４ａ 遷移可否判定部
１４ｂ 選択部
１５ 制約超過量計算部
１６ 制約識別符号入力部
１１ 状態保持部
１２ エネルギー変化計算部
１３ ペナルティ加算部
１４ 更新制御部
１４ａ 遷移可否判定部
１４ｂ 選択部
１５ 制約超過量計算部
１６ 制約識別符号入力部
【図1】
【図2】
【図3】
【図4】
【図5】
【図6】
【図7】
【図8】
【図9】
【図10】
【図11】
【図12】
【図13】
【図14】