JP2020194273A - 最適化装置、最適化方法及び最適化プログラム - Google Patents

Abstract

【課題】組合せ最適化問題の解の精度を向上する。【解決手段】探索部１２が、組合せ最適化問題を変換したイジングモデルについての基底状態の探索を行うことで、組合せ最適化問題の仮最適解を算出し、実行部１３が、仮最適解を用いたシミュレーションを実行し、シミュレーション結果の評価基準を示す評価基準値に基づいてシミュレーション結果を評価し、シミュレーション結果が評価基準を満たす場合、仮最適解を最適解として出力し、更新部１４が、シミュレーション結果が評価基準を満たさない場合、イジングモデルにシミュレーション結果に基づいた第１の制約項を追加した更新イジングモデルを生成し、更新イジングモデルについての基底状態の探索を探索部１２に実行させる。【選択図】図１

Description

本発明は、最適化装置、最適化方法及び最適化プログラムに関する。

現実社会のモデル化を行う手法として、シミュレータに様々な入力パラメータの値の組合せを与えて、多数回の試行（シミュレーション）を実行させる手法がある。この手法では、多数回の試行のそれぞれに対するシミュレーション結果の中で最もよい結果が得られたときの入力パラメータの値の組合せが最適解として決定される。このような手法の例として、ＰＩＤＯ（Process Integration and Design Optimization）というジャンルのソフトウェアを用いた手法が知られている。しかし、上記のような手法により現実社会に存在する複雑な組合せ最適化問題を解こうとする場合、１回の試行にある程度の時間がかかるため、実用的な時間で最適解に近い解を得ることが難しい。

実用的な時間で組合せ最適化問題の最適解に近い解を得る手法として、たとえば、シミュレーテッド・アニーリング法や、レプリカ交換法（交換モンテカルロ法とも呼ばれる）などのマルコフ連鎖モンテカルロ法が知られている。また、量子アニーリング法を用いる手法も知られている。これらの手法では、組合せ最適化問題が磁性体のスピンの振る舞いを表すモデルであるイジングモデルに置き換えられて計算される。イジングモデルを利用して組合せ最適化問題を解くハードウェアとして、デジタル回路を用いてシミュレーテッド・アニーリングを行うことでイジング型の評価関数の値が最小となる状態変数の値の組合せを計算するものがある。また、超伝導回路を用いて量子アニーリングを行うことで同様の計算を行うものなどがある。

なお、従来、物流の最適化を行うシミュレーション手法として、シミュレーション時間が増大することを抑制するために、制約の範囲内において物流パラメータの最適化を行う手法があった（たとえば、特許文献１参照）。

また、制約条件を考慮して船の荷役を最適化する際に、シミュレーテッド・アニーリング法を用いる手法があった（たとえば、特許文献２参照）。

特開２００２−２６９１９２号公報 特開平１１−１９９０５９号公報

現実社会に存在する複雑な組合せ最適化問題をイジングモデルにより忠実に定式化しようとする場合、パラメータ数が膨大になる。たとえば、時間経過を考慮する問題（動的な問題）をイジングモデルにより定式化する場合、時系列を表現するために、パラメータ数が時系列を表現しない場合に比べて激増する。このため、イジングモデルを用いた従来の手法ではハードウェア的な制約などにより解くことが困難な動的な問題を、静的な問題へと変換して計算することが考えられるが、精度のよい（最適解に近い）解を得ることができない。

１つの側面では、本発明は、組合せ最適化問題の解の精度を向上可能な最適化装置、最適化方法及び最適化プログラムを提供することを目的とする。

１つの実施態様では、組合せ最適化問題を変換したイジングモデルについての基底状態の探索を行うことで、前記組合せ最適化問題の仮最適解を算出する探索部と、前記仮最適解を用いたシミュレーションを実行し、前記シミュレーションの結果の評価基準を示す評価基準値に基づいて前記結果を評価し、前記結果が前記評価基準を満たす場合、前記仮最適解を最適解として出力する実行部と、前記結果が前記評価基準を満たさない場合、前記イジングモデルに前記結果に基づいた第１の制約項を追加した更新イジングモデルを生成し、前記更新イジングモデルについての基底状態の探索を前記探索部に実行させる更新部と、を有する最適化装置が提供される。

また、１つの実施態様では、最適化方法が提供される。
また、１つの実施態様では、最適化プログラムが提供される。

１つの側面では、組合せ最適化問題の解の精度を向上できる。

第１の実施の形態の最適化装置の一例を示す図である。 第１の実施の形態の最適化方法の一例の流れを示すフローチャートである。 第２の実施の形態の最適化装置のハードウェアの一例を示す図である。 第２の実施の形態の最適化装置のハードウェアの他の例を示す図である。 第２の実施の形態の最適化装置の機能例を示すブロック図である。 物流センター内における作業例を説明する図である。 あるトラックがある入荷品置き場に荷物を置いたときのフォークリフトの移動経路の例を示す図である。 あるトラックが別の入荷品置き場に荷物を置いたときのフォークリフトの移動経路の例を示す図である。 問題の制約条件を説明する図である。 第２の実施の形態の最適化装置による処理の一例の流れを示すフローチャートである。 ３つ目の制約条件を説明する図である。 シミュレーションの一例の流れを示すフローチャートである。 待ち時間が発生する例を示す図である。 待ち時間が発生する他の例を示す図である。 シミュレーションによる増加時間の観測結果を度数分布として集計した例を示す図である。

以下、発明を実施するための形態を、図面を参照しつつ説明する。
（第１の実施の形態）
図１は、第１の実施の形態の最適化装置の一例を示す図である。

第１の実施の形態の最適化装置１０は、生成部１１、探索部１２、実行部１３、更新部１４を有する。
生成部１１は、計算対象の組合せ最適化問題の問題データに基づいてイジングモデルを生成する（組合せ最適化問題をイジングモデルに変換する）。イジングモデルは、イジング型の評価関数（目的関数とも呼ばれる）にて表される。イジング型の評価関数は、ＱＵＢＯ（Quadratic Unconstrained Binary Optimization）と呼ばれる問題を解く際に用いられる２次関数で表すことができる。

計算対象の組合せ最適化問題の一例として、物流センターにおける作業時間を最短化するために、複数の入荷便トラック（以下単にトラックという）のそれぞれを複数の入荷品置き場の何れに割当てるか（停車させるか）を決定する問題がある。なお、計算対象の組合せ最適化問題としては、この例に限定されず、自動車の渋滞発生を抑制するために、交通量を最適化する問題など、様々な組合せ最適化問題への適用が可能である。

評価関数は、たとえば、図１に示すように、Ｈ＝コスト項ｃｓｔ＋制約項ｐ１と表すことができる。コスト項ｃｓｔは、最小化したい値であり、たとえば、物流センター内において荷物を複数の入荷品置き場のそれぞれから、複数の出荷品置き場の何れかに荷物を運ぶフォークリフトの総移動距離をコスト項ｃｓｔとすることができる。総移動距離が短くなれば、作業時間も短くなるためである。また、制約項ｐ１は、組合せ最適化問題が満たすべき制約条件を表し、制約条件を満たさない場合に値が大きくなる。たとえば、上記物流センターに関する問題では、物流センターに到着したトラックはどこか１箇所の入荷品置き場に入ること、という制約条件がある。制約条件が複数ある場合には、上記制約項ｐ１が複数、評価関数に含まれる。物流センターに関する問題では、他に、複数のトラックの停車時刻に重なりがある場合、それら複数のトラックは同じ入荷品置き場には入れない、という制約条件などがある。

なお、最適化装置１０は、外部の装置によって生成されたイジングモデルを取得し、図示しない記憶部に記憶しておいてもよい。その場合、生成部１１はなくてもよい。
探索部１２は、イジングモデルについての基底状態の探索を行うことで、組合せ最適化問題の仮最適解を算出する。探索部１２は、シミュレーテッド・アニーリング法やレプリカ交換法などのマルコフ連鎖モンテカルロ法により基底状態の探索を行ってもよいし、量子アニーリング法により基底状態の探索を行ってもよい。なお、探索結果として出力される仮最適解である状態（評価関数の全状態変数の値の組合せ）は、たとえば、所定時間内に多数回更新された状態のうち、評価関数の値が最小となる状態（基底状態とは限らない）である。上記物流センターに関する問題の例では、仮最適解は、フォークリフトの総移動距離をコスト項ｃｓｔとした評価関数Ｈの値が最小となる、複数のトラックの複数の入荷品置き場への割当て方を示す。この場合、各状態変数は、各トラックが各入荷品置き場に入るか否かを示す２値の変数である。

なお、組合せ最適化問題の種類によっては、制約項ｐ１に対応する制約条件として、必ず満たさなければいけない必須の制約条件と、できるだけ満たしたほうがよい任意の制約条件の、２種類をもつものがある。必須の制約条件が成立しない仮最適解は、現実には実行不可能な解を意味している。前述の物流センターの事例において、たとえば２台のトラックがある瞬間に同じ荷物置き場に入っている状態では、必須の制約条件が成立していない。このような種類の組合せ最適化問題では、必須の制約条件を満たさない仮最適解が棄却され、必須の制約条件を満たす解が得られるまで、探索が繰り返される。

実行部１３は、探索部１２が算出した仮最適解を用いたシミュレーションを実行し、シミュレーション結果の評価基準を示す評価基準値に基づいて、シミュレーション結果を評価する。そして、実行部１３は、シミュレーション結果が評価基準を満たす場合、仮最適解を最適解として出力する。まず実行部１３は、仮最適解を反映したシミュレーション入力データを生成し、その入力データに基づいて、シミュレーションを実行する。

たとえば、上記物流センターに関する問題では図１に示すように、複数のトラックによって入荷品置き場に置かれた全ての荷物が、フォークリフトによって複数の出荷品置き場の何れかに運ばれるまでの作業時間がシミュレートされる。その場合、たとえば、図１に示すようにあるフォークリフト１５ａが通過するまで、フォークリフト１５ｂは待機（または速度を減少）していなければならないという状況が生じる可能性がある。シミュレーションでは、このような状況が生じる場合のフォークリフト１５ｂの待ち時間などについても作業時間にカウントされる。

シミュレーション結果は、実際に最適化したい対象（最適化対象）の数値データを含むものであり、たとえば、上記物流センターに関する問題では物流センター内における荷物運びの作業時間である。その場合、実行部１３は、シミュレーション結果として得られた作業時間が、外部から与えられる評価基準値である所定時間より短ければ、評価基準を満たすものとして、仮最適解を最適解として出力する。

なお、最適化対象は、上記のコスト項ｃｓｔが示す対象とは必ずしも一致しない。たとえば、上記のようにコスト項ｃｓｔが、フォークリフトの総移動距離であるのに対して、最適化対象は作業時間である。

また、上記「最適解」は、シミュレーション結果が最初に評価基準を満たしたときの仮最適解（探索部１２による探索結果）を示し、厳密な意味での最適解とは必ずしも一致しない。

更新部１４は、シミュレーション結果が評価基準を満たさない場合、イジングモデルにシミュレーション結果に基づいた制約項（これは前述の任意の制約条件に分類される）を追加した更新イジングモデルを生成し、更新イジングモデルについての基底状態の探索を探索部１２に実行させる。制約項は、シミュレーションの実行中に観測された事象を記録した数値データ（以下イベントデータという）に基づいて生成される。図１には、上記物流センターに関する問題における、イベントデータの一例が示されている。イベントデータには、たとえば、イベント番号と日時、イベントの内容、フォークリフト情報、ロスタイムが含まれる。イベント番号＝１のイベントは、トラックＴ１０が入荷品置き場Ｐ２に置いた荷物を運搬中のフォークリフト１５ａが交差点を通過するのを、トラックＴ１５が入荷品置き場Ｐ３に置いた荷物を運搬中のフォークリフト１５ｂが７秒待ったことを示す。

このような解を悪化させる要素（待ち時間など）が生じる、状態変数の値の組合せについては評価関数の値が大きくなるように、制約項が生成される。図１では、更新イジングモデルを表す評価関数Ｈａに、制約項ｐ２が追加された例が示されている。

生成部１１、探索部１２、実行部１３及び更新部１４は、たとえば、ＣＰＵ（Central Processing Unit）やＤＳＰ（Digital Signal Processor）などのプロセッサが実行するプログラムモジュールを用いて実装できる。なお、探索部１２は、デジタル回路を用いてシミュレーテッド・アニーリング法やレプリカ交換法などを実行する専用ハードウェアであってもよいし、量子アニーリングを行う専用ハードウェアであってもよい。

以下第１の実施の形態の最適化装置１０の動作（最適化方法）の例を説明する。
図２は、第１の実施の形態の最適化方法の一例の流れを示すフローチャートである。
生成部１１は、計算対象の組合せ最適化問題の問題データに基づいてイジングモデルを生成し（ステップＳ１）、探索部１２は、イジングモデルについての基底状態の探索を行うことで、組合せ最適化問題の仮最適解を算出する（ステップＳ２）。

その後、実行部１３は、探索部１２が算出した仮最適解を用いたシミュレーションを実行し（ステップＳ３）、評価基準値に基づいて、シミュレーション結果が評価基準を満たすか否かを判定する（ステップＳ４）。

シミュレーション結果が評価基準を満たさない場合、更新部１４は、イジングモデルに制約項を追加した更新イジングモデルを生成する（イジングモデルを更新する）（ステップＳ５）。ステップＳ５の処理後、生成された更新イジングモデルを用いて、ステップＳ２からの処理が繰り返される。

シミュレーション結果が評価基準を満たす場合、実行部１３は、仮最適解を最適解として出力し（ステップＳ６）、最適化装置１０の動作が終了する。
以上のように、最適化装置１０は、イジングモデルに変換した組合せ最適化問題の仮最適解を用いてシミュレーションを行い、その結果に基づいた制約項を加えてイジングモデルを更新し、仮最適解の探索を繰り返す。これにより、イジングモデルが、実際の組合せ最適化問題の特性（動的な特性など）をより反映させたものになり、解の精度を向上できる。

また、イジングモデルを用いて計算を行うため、シミュレーションだけで最適解を探索する場合よりもシミュレーションの実行回数を削減でき、最適解または最適解に近い解を短時間で求めることが可能となる。

（第２の実施の形態）
図３は、第２の実施の形態の最適化装置のハードウェアの一例を示す図である。
第２の実施の形態の最適化装置２０は、たとえばコンピュータであり、ＣＰＵ２１、ＲＡＭ（Random Access Memory）２２、ＨＤＤ（Hard Disk Drive）２３、画像信号処理部２４、入力信号処理部２５、媒体リーダ２６及び通信インタフェース２７を有する。上記ユニットは、バスに接続されている。

ＣＰＵ２１は、プログラムの命令を実行する演算回路を含むプロセッサである。ＣＰＵ２１は、ＨＤＤ２３に記憶されたプログラムやデータの少なくとも一部をＲＡＭ２２にロードし、プログラムを実行する。なお、ＣＰＵ２１は複数のプロセッサコアを備えてもよく、最適化装置２０は複数のプロセッサを備えてもよく、以下で説明する処理を複数のプロセッサまたはプロセッサコアを用いて並列に実行してもよい。また、複数のプロセッサの集合（マルチプロセッサ）を「プロセッサ」と呼んでもよい。

ＲＡＭ２２は、ＣＰＵ２１が実行するプログラムやＣＰＵ２１が演算に用いるデータを一時的に記憶する揮発性の半導体メモリである。なお、最適化装置２０は、ＲＡＭ以外の種類のメモリを備えてもよく、複数個のメモリを備えてもよい。

ＨＤＤ２３は、ＯＳ（Operating System）やミドルウェアやアプリケーションソフトウェアなどのソフトウェアのプログラム、及び、データを記憶する不揮発性の記憶装置である。プログラムには、たとえば、組合せ最適化問題を計算する最適化プログラムが含まれる。なお、最適化装置２０は、フラッシュメモリやＳＳＤ（Solid State Drive）などの他の種類の記憶装置を備えてもよく、複数の不揮発性の記憶装置を備えてもよい。

画像信号処理部２４は、ＣＰＵ２１からの命令にしたがって、最適化装置２０に接続されたディスプレイ２４ａに画像を出力する。ディスプレイ２４ａとしては、ＣＲＴ（Cathode Ray Tube）ディスプレイ、液晶ディスプレイ（ＬＣＤ：Liquid Crystal Display）、プラズマディスプレイ（ＰＤＰ：Plasma Display Panel）、有機ＥＬ（ＯＥＬ：Organic Electro-Luminescence）ディスプレイなどを用いることができる。

入力信号処理部２５は、最適化装置２０に接続された入力デバイス２５ａから入力信号を取得し、ＣＰＵ２１に出力する。入力デバイス２５ａとしては、マウスやタッチパネルやタッチパッドやトラックボールなどのポインティングデバイス、キーボード、リモートコントローラ、ボタンスイッチなどを用いることができる。また、最適化装置２０に、複数の種類の入力デバイスが接続されていてもよい。

媒体リーダ２６は、記録媒体２６ａに記録されたプログラムやデータを読み取る読み取り装置である。記録媒体２６ａとして、たとえば、磁気ディスク、光ディスク、光磁気ディスク（ＭＯ：Magneto-Optical disk）、半導体メモリなどを使用できる。磁気ディスクには、フレキシブルディスク（ＦＤ：Flexible Disk）やＨＤＤが含まれる。光ディスクには、ＣＤ（Compact Disc）やＤＶＤ（Digital Versatile Disc）が含まれる。

媒体リーダ２６は、たとえば、記録媒体２６ａから読み取ったプログラムやデータを、ＲＡＭ２２やＨＤＤ２３などの他の記録媒体にコピーする。読み取られたプログラムは、たとえば、ＣＰＵ２１によって実行される。なお、記録媒体２６ａは、可搬型記録媒体であってもよく、プログラムやデータの配布に用いられることがある。また、記録媒体２６ａやＨＤＤ２３を、コンピュータ読み取り可能な記録媒体ということがある。

通信インタフェース２７は、ネットワーク２７ａに接続され、ネットワーク２７ａを介して他の情報処理装置と通信を行うインタフェースである。通信インタフェース２７は、スイッチなどの通信装置とケーブルで接続される有線通信インタフェースでもよいし、基地局と無線リンクで接続される無線通信インタフェースでもよい。

図４は、第２の実施の形態の最適化装置のハードウェアの他の例を示す図である。図４において、図３に示した要素と同じ要素については同一符号が付されている。
最適化装置３０は、情報処理装置２０ａとイジングマシン２８ａを有する。情報処理装置２０ａは、インタフェース２８を有する。インタフェース２８は、イジングマシン２８ａに接続され、ＣＰＵ２１とイジングマシン２８ａとの間でデータの送受信を行う。インタフェース２８は、たとえば、ＰＣＩ（Peripheral Component Interconnect） Ｅｘｐｒｅｓｓなどの有線通信インタフェースでもよいし、無線通信インタフェースでもよい。

イジングマシン２８ａは、デジタル回路を用いてシミュレーテッド・アニーリング法やレプリカ交換法などを実行する専用ハードウェアであってもよいし、量子アニーリングを行う専用ハードウェアであってもよい。

次に、最適化装置２０，３０の機能及び処理手順を説明する。
図５は、第２の実施の形態の最適化装置の機能例を示すブロック図である。
なお、以下では、図３に示した最適化装置２０の機能例を示すが、図４に示した最適化装置３０についても同様の機能を有する。

最適化装置２０は、入力部３１、記憶部３２、生成部３３、探索部３４、シミュレーション部３５、制御部３６、更新部３７、出力部３８を有する。入力部３１、生成部３３、探索部３４、シミュレーション部３５、制御部３６、更新部３７、出力部３８は、たとえば、ＣＰＵ２１が実行するプログラムモジュールを用いて実装できる。なお、図４に示した最適化装置３０が用いられる場合、イジングマシン２８ａが探索部３４（またはその一部）として機能する。記憶部３２は、たとえば、ＲＡＭ２２またはＨＤＤ２３に確保した記憶領域を用いて実装できる。

入力部３１は、たとえば、入力デバイス２５ａによって入力されるデータ（入力データ）を取得する。入力データには、イジングモデルを生成するために用いられる問題データ、シミュレーションを実行する際に用いられるデータ、シミュレーション結果の評価基準を示す評価基準値が含まれる。入力部３１が取得した入力データは、記憶部３２に記憶される。

記憶部３２は、上記の入力データを記憶する。また、記憶部３２は、探索部３４によって得られる仮最適解、シミュレーション部３５によって得られるシミュレーション結果を記憶してもよい。

生成部３３は、記憶部３２に記憶された入力データに基づいて、イジングモデルを生成する。以下では、イジングモデルは、ＱＵＢＯの２次関数で表されるものとする。生成したイジングモデルの情報は、記憶部３２に記憶される。

探索部３４は、イジングモデルの情報を生成部３３（または記憶部３２）から受け、イジングモデルについての基底状態の探索を行うことで、組合せ最適化問題の仮最適解を算出する。探索部３４は、シミュレーテッド・アニーリング法やレプリカ交換法などのマルコフ連鎖モンテカルロ法により基底状態の探索を行ってもよいし、量子アニーリング法により基底状態の探索を行ってもよい。

シミュレーション部３５は、探索部３４が算出した仮最適解を用いたシミュレーションを実行し、シミュレーション結果を出力する。シミュレーション部３５は、制御部３６が仮最適解に基づいて生成したシミュレーションのための入力データと、記憶部３２に記憶されているシミュレーションのための別の入力データと、に基づいてシミュレーションを実行する。シミュレーション結果には、最適化対象の数値データのほかに、シミュレーション実行中に観測された事象を記録したイベントデータが含まれる。

制御部３６は、最適化装置２０の各部を制御するとともに、探索部３４が算出した仮最適解に基づいて、シミュレーションのための入力データを生成し、シミュレーション部３５に供給する。さらに、制御部３６は、記憶部３２に記憶されている評価基準値に基づいて、シミュレーション部３５が出力するシミュレーション結果を評価する。そして、制御部３６は、シミュレーション結果が評価基準を満たす場合、仮最適解を最適解として出力部３８に出力させる。このように、図５に示す第２の実施の形態の最適化装置２０の機能例では、シミュレーション部３５、制御部３６、出力部３８によって、第１の実施の形態の最適化装置１０の実行部１３の機能が行われる。

また、制御部３６は、シミュレーション結果が評価基準を満たさない場合、更新部３７に、イジングモデルの更新を行わせる。
更新部３７は、シミュレーション結果が評価基準を満たさない場合、シミュレーション結果に含まれるイベントデータに基づいてイジングモデルに追加する制約項を生成する。そして、更新部３７は、記憶部３２に記憶されている現在のイジングモデルの情報に基づいて、現在のイジングモデルに生成した制約項を追加することで、イジングモデルを更新する。更新部３７は、更新したイジングモデルの情報を探索部３４に供給し、更新したイジングモデルについての基底状態の探索を、探索部３４に実行させる。

出力部３８は、シミュレーション結果が評価基準値を満たした場合の、仮最適解を最適解として、たとえば、ディスプレイ２４ａに出力し表示させる。出力部３８は、最適解を、記憶部３２に記憶させてもよい。

次に、第２の実施の形態の最適化装置２０，３０の動作（最適化方法）の例を説明する。
なお、以下では、計算対象の組合せ最適化問題の一例として、荷物が運び込まれるある施設内における作業時間を最短化するために、荷物を運び込む複数の車両のそれぞれを複数の入荷品置き場の何れに割当てるかを決定する問題を例に挙げて説明する。施設内では、移動体により荷物を移動させる作業があるものとする。

なお、以下では荷物が運び込まれる施設を物流センター、荷物を運び込む車両をトラック、荷物を移動させる移動体をフォークリフトと呼ぶが、これらに限定されるわけではなく。たとえば、移動体は人であってもよい。

図６は、物流センター内における作業例を説明する図である。
複数のトラックのそれぞれは、複数の入荷品置き場（図６の例では１〜９の９つ）の何れかに荷物（入荷品）を置く。なお、トラックが発着する場所のことを、ピットと呼ぶ。以下では、ピットと入荷品置き場とは１対１で対応しているものとして、ピット番号（図６の例では１〜９）により入荷品置き場を識別するものとする。

各荷物には、たとえば、実際の荷物のあて先に対応した出荷品置き場のあて先（図６の例ではＡ〜Ｋ）が割当てられている。人または運搬装置（以下、フォークリフトとする）は、各入荷品置き場から荷物を１つ持ち出して、あて先に対応する出荷品置き場に移動する。この作業を、各入荷品置き場に荷物がある限り繰り返す。

なお、トラックが到着する時刻、出発する時刻、及び運び込む荷物の個数とあて先は、トラックごとに予め決められているものとする。また、決められた時刻に到着して、出荷品置き場から荷物を持ち出す出荷便トラックも存在するが、入荷便トラックとは別のピットに到着し、入荷と出荷が干渉しあうことはないものとする。

このような問題では、生成部３３が生成するイジングモデルは、フォークリフトの総移動距離をコスト項として含むＱＵＢＯの２次関数となる。そして、各トラックが、どの入荷品置き場に荷物を下ろすのかについての組合せのうち、ＱＵＢＯの２次関数の値が最小になる組合せがイジングモデルの基底状態に相当する。

図７は、あるトラックがある入荷品置き場に荷物を置いたときのフォークリフトの移動経路の例を示す図である。なお、図では移動経路は直線で示されているが、直線には限らない。

図７の例では、トラック番号＝ｉのトラック（図では“トラックｉ”と表記されている（以降の図も同じ））が、ピット番号＝１の入荷品置き場に、５つの荷物を置いた例が示されている。５つの荷物のあて先は、Ａ，Ｂ，Ｃ，Ｋ，Ｋである。この場合の、フォークリフトの移動距離の合計は、Ｃ_ｉ，１＝ｄｉｓｔ（１，Ａ）＋ｄｉｓｔ（１，Ｂ）＋ｄｉｓｔ（１，Ｃ）＋ｄｉｓｔ（１，Ｋ）×２と表せる。ｄｉｓｔ（１，Ａ）は、フォークリフトが、ピット番号＝１の入荷品置き場からあて先＝Ａに対応した出荷品置き場に荷物を運ぶ場合の移動距離を示す。ｄｉｓｔ（１，Ｂ）は、フォークリフトが、ピット番号＝１の入荷品置き場からあて先＝Ｂに対応した出荷品置き場に荷物を運ぶ場合の移動距離を示す。ｄｉｓｔ（１，Ｃ）は、フォークリフトが、ピット番号＝１の入荷品置き場からあて先＝Ｃに対応した出荷品置き場に荷物を運ぶ場合の移動距離を示す。ｄｉｓｔ（１，Ｋ）は、フォークリフトが、ピット番号＝１の入荷品置き場からあて先＝Ｋに対応した出荷品置き場に荷物を運ぶ場合の移動距離を示す。

これらの各入荷品置き場と各あて先の出荷品置き場との間における移動距離は、入力データに含まれ、記憶部３２に記憶されている。なお、上記のようなＣ_ｉ，１は、生成部３３が生成するイジングモデルの情報に含まれる重み係数を示す。

図８は、あるトラックが別の入荷品置き場に荷物を置いたときのフォークリフトの移動経路の例を示す図である。
図８の例では、トラック番号＝ｉのトラックが、ピット番号＝９の入荷品置き場に、図７の例と同じあて先の５つの荷物を置いた例が示されている。この場合の、フォークリフトの移動距離の合計は、Ｃ_ｉ，９＝ｄｉｓｔ（９，Ａ）＋ｄｉｓｔ（９，Ｂ）＋ｄｉｓｔ（９，Ｃ）＋ｄｉｓｔ（９，Ｋ）×２と表せる。もし、Ｃ_ｉ，９＜Ｃ_ｉ，１であるならば、トラック番号＝ｉのトラックは、ピット番号＝１の入荷品置き場に荷物を置くよりも、ピット番号＝９の入荷品置き場に荷物を置いたほうが、総移動距離を短くできるため、作業時間を短縮できる可能性が高い。

ただし、実際は複数のフォークリフトによる干渉によって待ち時間が発生し、作業時間が長くなる可能性もある。これらの点がシミュレーションによって検証される。そして、シミュレーション結果に応じて、大きな待ち時間が発生するような各トラックとピットの割当ての組合せについては、ＱＵＢＯの２次関数の値が大きくなるような制約項がイジングモデルに追加され、探索部３４による探索が繰り返される。

図９は、問題の制約条件を説明する図である。
図９では、各ピット（入荷品置き場）へのトラックの割当て例が示されている。図９において、横軸は時刻を表し、縦軸はピット番号を表している。

前述のように、各トラックが物流センターに到着する時刻及び出発する時刻は予め決まっている。この問題における制約条件として以下のようなものがある。１つ目は、物流センターに到着したトラックは必ずどこかのピットに入る、という制約条件である。２つ目は、１つのピットに停車時刻に重なりがあるトラックが同時に２台入らない、という制約条件である。たとえば、図９のように、ピット番号＝１の入荷品置き場（ピット）にトラック番号＝１のトラックが存在するときに、そのトラックと停車時刻が重なる他のトラックは、そのピットに入ることができない。図９の例では、トラック番号＝１のトラックが出発した後に、トラック番号＝１０のトラックが、ピット番号＝１の入荷品置き場に入っている。

探索部３４が、演算可能なビット数を考慮して、所定時間のタイムスロット内における最適化を行う場合、以下のような制約条件も加わる。すなわち、直前のタイムスロットから処理対象のタイムスロットに渡ってトラックが入っているピットには、停車時刻が重なるトラックは入れないという制約条件である。

以上のような問題に対する、最適化装置２０，３０の処理の流れを説明する。
図１０は、第２の実施の形態の最適化装置による処理の一例の流れを示すフローチャートである。

まず、入力部３１は、入力データを取得する（ステップＳ１１）。入力部３１が取得した入力データは、記憶部３２に記憶される。入力データには、イジングモデルを生成するためのデータが含まれる。イジングモデルを生成するためのデータとして、各トラックの到着時刻、出発時刻、及び各トラックが運び込む荷物の個数とあて先、各入荷品置き場と各あて先の出荷品置き場との間における移動距離の情報がある。さらに、イジングモデルを生成するためのデータとして、制約条件に関する情報、探索部３４が演算可能なビット数の情報がある。また、入力データには、シミュレーションを実行する際に用いられるデータとして、フォークリフトの数、フォークリフトの標準速度、各入荷品置き場から各あて先の出荷品置き場までのフォークリフトの移動経路を示すマップ情報などが含まれる。標準速度は、フォークリフト同士が干渉することがないとき、フォークリフトが自由に物流センター内を走行できるときの平均的な速度を示す。さらに、入力データには、シミュレーション結果の評価基準を示す評価基準値として、最適とみなせる作業時間が含まれる。

その後、生成部３３は入力データに基づいて、イジングモデルを生成する（ステップＳ１２）。イジングモデルはＱＵＢＯの２次元関数で定式化される。
ＱＵＢＯの２次関数に含まれる状態変数は、ｘ_ｉ，ｐと表せる。ｉは、トラック番号であり、ｐはピット番号である。ｘ_ｉ，ｐ＝１の場合、ピット番号＝ｐの入荷品置き場に、トラック番号＝ｉのトラックが入ることを表し、ｘ_ｉ，ｐ＝０の場合、ピット番号＝ｐの入荷品置き場に、トラック番号＝ｉのトラックが入らないことを表す。２次関数に含まれる状態変数の数は、トラックの台数×入荷品置き場の数（ピットの数）である。

たとえば、物流センターにおいて、１日あたりｎ台（たとえば、数百台）のトラックが到着し、それらが、２３箇所（図７の例では９箇所）のピットのどこかに入るものとする。その場合、各トラックが、どの入荷品置き場に荷物を下ろすのかについての組合せの数は、２３^ｎ通りある。探索部３４が演算可能なビット数（状態変数の数に相当する）が、たとえば、１０２４ビットである場合、トラック数が数百台ともなると、探索部３４は一度に演算ができない。

しかし、トラックが物流センターに停車する時間は、１０〜５０分程度であるので、ｎ台全てを一度に考慮しなくてもよく、たとえば、２４時間を２時間単位のタイムスロットに区切って解いても、解はそれほど悪くはならないと思われる。入荷品置き場の数が２３であり、演算可能なビット数が１０２４ビットである場合、１０２４／２３＝４４．５…となるため、トラックの台数が４４台以下になるようにタイムスロットを調整すればよい。

なお、以下では、生成部３３は、物流センターへの到着時刻が早い順を第１キーにして発車時刻が早い順を第２キーにしてトラックをソートして、探索部３４が演算可能なビット数に収まる台数のトラックを到着時刻順に抽出し、トラック番号を割当てて、定式化を行うものとする。

ＱＵＢＯの２次関数は、以下の式（１）で表せる。

式（１）において、右辺の１項目は、コスト項であり、フォークリフトの総移動距離を表す。Ｃ_ｉ，ｐは、前述のような重み係数であり、トラック番号＝ｉのトラックが、ピット番号＝ｐの入荷品置き場に荷物を置いた場合の、フォークリフトの移動距離の合計を表す。

右辺の２項目は、物流センターに到着したトラックは必ずどこかのピットに入る、という制約条件に対応した制約項である。たとえば、トラック番号＝ｉのトラックが、どのピットにも入らない場合（ｘ_ｉ，ｐ＝１となるピット番号＝ｐがない場合）、この制約項は０よりも大きい値をもつことになる。

右辺の３項目は、１つのピットに停車時刻に重なりがあるトラックが２台入らない、という制約条件に対応した制約項である。３項目のｉ，ｊは、停車時刻に重なりがある２つのトラックのトラック番号を示す。たとえば、トラック番号＝ｉ，ｊの２つのトラックが、同じピット番号＝ｐのピットに入る場合（ｘ_ｉ，ｐ＝ｘ_ｊ，ｐ＝１の場合）、この制約項は０よりも大きい値をもつことになる。

右辺の４項目は、直前のタイムスロットから処理対象のタイムスロットに渡ってトラックが入っているピットには、停車時刻が重なるトラックは入れないという制約条件に対応した制約項である。

図１１は、３つ目の制約条件を説明する図である。
たとえば、図１１は、図９にタイムスロット４０，４１を追加したものである。タイムスロット４０，４１は、ＱＵＢＯの２次関数の状態変数の数が、探索部３４が演算可能なビット数以内になるように決定されている。タイムスロット４０において、トラックが入るピットの割当てを決定する場合、対象となるトラックは、トラック番号＝１１，１２，１４のトラックである。なお、直前のタイムスロット４１における処理により、トラック番号＝１３，１５，１８，１９，２０，２１のトラックについては、ピットの割当てが決定している。

たとえば、トラック番号＝１１のトラックは、トラック番号＝１３，１５，１８，１９，２０，２１のトラックと停車時刻が重なるため、それらがいるピット番号＝２，３，６，７，８，９のピットには入れない。トラック番号＝１２のトラックと停車時刻が重なるのは、トラック番号＝２１のトラックだけであるので、ピット番号＝９のピット以外に入れる。

式（１）の右辺の４項目において、“Ｒｅｍａｉｎ”は、直前のタイムスロットにおける処理によりトラックが割当てられたピットであり、その占有時間が、トラック番号＝ｉのトラックの停車時刻と重なっているピットのピット番号の集合である。図１１の例では、トラック番号＝１１に対する“Ｒｅｍａｉｎ”は、ピット番号＝２，３，６，７，８，９であり、トラック番号＝１２に対する“Ｒｅｍａｉｎ”は、ピット番号＝９である。

図１１の例では、たとえば、ｘ_１２，９＝１の場合、この制約項が０よりも大きくなる。
式（１）において、α，β，γは各制約項の係数であり、入力データに含まれており、比較的大きな値が設定される。

以上のような定式化が終わると、探索部３４は、定式化されたイジングモデルについての基底状態の探索を行うことで、組合せ最適化問題の仮最適解を算出する（ステップＳ１３）。

探索部３４の機能を、デジタル回路を用いてシミュレーテッド・アニーリングを行うイジングマシン２８ａが行う場合、たとえば、以下のような処理が行われる。イジングマシン２８ａは、ある状態変数が変化した場合にその変化に伴って各状態変数が変化することによるエネルギー変化（式（１）のＥ（ｘ）の変化分）をそれぞれ計算する。そして、イジングマシン２８ａは、温度パラメータと乱数値に基づいて生成される閾値と各状態変数が変化することによるエネルギー変化とを比較する。さらに、イジングマシン２８ａは、その比較結果に基づいて、Ｅ（ｘ）が小さくなるようなエネルギー変化を引き起こす状態遷移が優先的に選択されるように変化させる状態変数を決定する。イジングマシン２８ａは、このような処理を温度パラメータの値を徐々に小さくしながら多数回繰り返すことで、基底状態の探索を行う。

なお、ステップＳ１３では、仮最適解が必須の制約条件が成立しているかチェックされる。必須の制約条件は、式（１）の２項目、３項目、４項目に対応しており、これらの項に、仮最適解における状態変数の値を代入して、値が０にならないときは、必須の制約条件が成立していない、と判断できる。仮最適解が必須の制約条件を満たしていないときは、基底状態の探索が継続され、必須の制約条件が成立する別の仮最適解が見つかるまで、探索が繰り返される。ここで、繰り返し回数の上限を設定しておいて、上限を超えたときは、エラー終了させてもよい。なお、必須の制約条件が成立しづらいときの原因として、入力として与えられる式（１）のα，β，γが不適切な値となっている可能性がある。

ステップＳ１３の処理後、制御部３６は、探索部３４が算出した仮最適解に基づいて、シミュレーションのための入力データを生成する（ステップＳ１４）。たとえば、制御部３６は、仮最適解において値が１である状態変数に基づいて、各入荷品置き場に、トラックの荷物を対応付けた入力データを生成する。

たとえば、ｘ_ｉ，１＝１の場合、制御部３６は、図７に示した例のように、ピット番号＝１の入荷品置き場に、トラック番号＝ｉのトラックの５つの荷物を対応付ける。ｘ_ｉ，９＝１の場合、制御部３６は、図８に示した例のように、ピット番号＝９の入荷品置き場に、トラック番号＝ｉのトラックの５つの荷物を対応付ける。

ステップＳ１４の処理後、シミュレーション部３５は、制御部３６が生成した入力データと、記憶部３２に記憶されている別の入力データ（フォークリフトの数、標準速度、マップ情報など）と、に基づいてシミュレーションを実行する（ステップＳ１５）。

図１２は、シミュレーションの一例の流れを示すフローチャートである。
なお、シミュレーションにおいて、フォークリフトは複数台あってもよく、フォークリフトごとに、入荷品置き場から出荷品置き場への荷物を運搬する役割が与えられているものとする。各フォークリフトは自律的に動作する。また、各フォークリフトは走行ルールを守る義務もある。走行ルールとしては、制限速度を守る、一時停止をする、走行可能な場所のみを走行する、他のフォークリフトとの接触を避けるように速度をコントロールする、入出荷品置き場に進入できるのは１台まで、などがある。このフォークリフトの動作と位置などが、たとえば０．１秒刻みなどの微小時間ごとに、プログラムによって繰り返し計算されることによって、シミュレーションを行うことが可能である。

シミュレーション部３５は、シミュレーションを開始すると、複数台のフォークリフトによる干渉が発生するたびに、待ち時間や速度の低下に関する情報を、その干渉を起こす状態変数の組とともに、たとえば、記憶部３２に記憶させる（ステップＳ２０）。そして、シミュレーション部３５は、入荷品置き場から全荷物が出荷品置き場に運搬が完了したか否かを判定する（ステップＳ２１）。入荷品置き場から全荷物が出荷品置き場に運搬が完了していない場合には、ステップＳ２０からの処理が繰り返される。入荷品置き場から全荷物が出荷品置き場に運搬が完了した場合、シミュレーション部３５は作業時間と記録した待ち時間や速度低下情報を含むシミュレーション結果を出力し（ステップＳ２２）、シミュレーションを終了する。なお、シミュレーションは、１つの仮最適解に基づいて生成された入力データに基づいて、複数のシミュレーション条件で行われるようにしてもよい。複数のシミュレーション条件は、たとえば、フォークリフトの台数を増減させることで、作成することが可能である。

図１３は、待ち時間が発生する例を示す図である。
たとえば、仮最適解に含まれるｘ_１０，２＝ｘ_１５，３＝１の場合、トラック番号＝１０のトラックが、ピット番号＝２のピットに置いた荷物を運搬するフォークリフト５０ａが存在する。また、トラック番号＝１５のトラックが、ピット番号＝３のピットに置いた荷物を運搬するフォークリフト５０ｂが存在する。図１３の例では、交差点においてフォークリフト５０ａが、フォークリフト５０ｂの移動経路をふさいでいるため、フォークリフト５０ｂは、フォークリフト５０ａが交差点を追加するまで停止している。シミュレーション部３５は、このときの待ち時間（図の例ではｔ秒）とｘ_１０，２，ｘ_１５，３とを対応付けて、たとえば、記憶部３２に記憶させる。

図１４は、待ち時間が発生する他の例を示す図である。
たとえば、仮最適解に含まれるｘ_１１，４＝ｘ_１６，５＝１の場合、トラック番号＝１１のトラックが、ピット番号＝４のピットに置いた荷物を運搬するフォークリフト５１ａが存在する。また、トラック番号＝１６のトラックが、ピット番号＝５のピットに置いた荷物を運搬するフォークリフト５１ｂが存在する。図１４の例では、フォークリフト５１ａが、出荷品置き場で荷物を降ろしているときに、同じ出荷品置き場に入るフォークリフト５１ｂがきた場合、フォークリフト５１ｂは、フォークリフト５１ａが出荷品置き場から出るまで停止している。シミュレーション部３５は、このときの待ち時間（図の例ではｔａ秒）とｘ_１１，４，ｘ_１６，５とを対応付けて、たとえば、記憶部３２に記憶させる。

以上のようなシミュレーションの実行が終了すると、制御部３６は、記憶部３２に記憶されている評価基準値に基づいて、シミュレーション結果が評価基準を満たすか否かを判定する（ステップＳ１６）。たとえば、制御部３６は、シミュレーション結果に含まれる作業時間が、評価基準である作業時間よりも短い場合、評価基準を満たすと判定する。

シミュレーション結果が評価基準を満たさない場合、更新部３７は、イジングモデルの更新を行う（ステップＳ１７）。
なお、評価基準値は、適宜変更されるようにしてもよい。たとえば、所定回数のイジングモデルの更新が行われても、シミュレーション結果の作業時間が評価基準を満たさない（評価基準値より大きい）場合、制御部３６は、ディスプレイ２４ａに評価基準値を大きくすることを促すようなメッセージを表示させてもよい。また、シミュレーション結果の作業時間が評価基準を比較的少ない更新回数で満たす場合、より短い作業時間となる解がある可能性があるため、制御部３６は、ディスプレイ２４ａに評価基準値を小さくすることを促すようなメッセージを表示させてもよい。

図１５は、シミュレーションによる増加時間の観測結果を度数分布として集計した例を示す図である。図１５において、横軸は増加時間（フォークリフト同士の干渉が１回発生たときの、［シミュレーションで得られた総作業時間］−［待ち時間などがないときの作業時間］）を表し、縦軸は度数（その増加時間が観測された回数）を表す。なお、フォークリフト同士の干渉とは、図１３、図１４のような状態のことであり、このような状態が発生するたびに、１回と数える。

更新部３７は、たとえば、図１５の増加時間が大きなサンプルにおいて、待ち時間や速度低下の情報と、その待ち時間や速度低下を引き起こす状態変数の組を、記憶部３２から抽出する。そして、更新部３７は、待ち時間や速度の低下に関する情報を距離の次元の値に変換して、抽出した状態変数の組を用いて新たな制約項を生成し、式（１）に示した２次関数に追加することでイジングモデルの更新を行う。図１３に示した例では、ｘ_１０，２＝ｘ_１５，３＝１の場合に、待ち時間ｔが生じるため、制約項として、ｘ_１０，２×ｘ_１５，３×ｔ×標準速度を、２次関数に追加する。速度の低下が発生するときについても同様に、低下した速度とその速度で走行した時間との積に基づいて、制約項が生成される。なお、ここで追加される制約項は、できるだけ避けたい現象（フォークリフトの干渉）を意味しているので、必須の制約条件ではなく任意の制約条件である。また、制約項では標準速度の代わりに、式（１）のα，β，γと同様の、重みを加味した値を用いてもかまわない。

ステップＳ１７の処理後、更新されたイジングモデルを用いてステップＳ１３からの処理が繰り返される。
ステップＳ１６においてシミュレーション結果が評価基準を満たすと判定された場合、出力部３８は、その時点での仮最適解を最適解として出力し（ステップＳ１８）、最適化処理を終える。

ステップＳ１７の処理により更新されたイジングモデルは、シミュレーションで検証された問題の動的な特性を反映しているため、探索部３４によって得られる新たな仮最適解も、問題の動的な特性を反映した解となる。上記のような処理を繰り返すことで、イジングモデルが、実際の組合せ最適化問題をより反映させたものにすることが可能となり、解の精度を向上できる。

なお、上記では物流センターに関する問題を例にして説明したが、この問題に限定されず他の問題にも適用可能である。
なお、前述のように、上記の処理内容は、最適化装置２０，３０にプログラム（最適化プログラム）を実行させることで実現できる。

プログラムは、コンピュータ読み取り可能な記録媒体（たとえば、記録媒体２６ａ）に記録しておくことができる。記録媒体として、たとえば、磁気ディスク、光ディスク、光磁気ディスク、半導体メモリなどを使用できる。磁気ディスクには、ＦＤ及びＨＤＤが含まれる。光ディスクには、ＣＤ、ＣＤ−Ｒ（Recordable）／ＲＷ（Rewritable）、ＤＶＤ及びＤＶＤ−Ｒ／ＲＷが含まれる。プログラムは、可搬型の記録媒体に記録されて配布されることがある。その場合、可搬型の記録媒体から他の記録媒体（たとえば、ＨＤＤ２３）にプログラムをコピーして実行してもよい。

以上、実施の形態に基づき、本発明の最適化装置、最適化方法及び最適化プログラムの一観点について説明してきたが、これらは一例にすぎず、上記の記載に限定されるものではない。

１０ 最適化装置
１１ 生成部
１２ 探索部
１３ 実行部
１４ 更新部
１５ａ，１５ｂ フォークリフト

Claims (6)

1. 組合せ最適化問題を変換したイジングモデルについての基底状態の探索を行うことで、前記組合せ最適化問題の仮最適解を算出する探索部と、
   前記仮最適解を用いたシミュレーションを実行し、前記シミュレーションの結果の評価基準を示す評価基準値に基づいて前記結果を評価し、前記結果が前記評価基準を満たす場合、前記仮最適解を最適解として出力する実行部と、
   前記結果が前記評価基準を満たさない場合、前記イジングモデルに前記結果に基づいた第１の制約項を追加した更新イジングモデルを生成し、前記更新イジングモデルについての基底状態の探索を前記探索部に実行させる更新部と、
   を有する最適化装置。
2. 前記第１の制約項は、前記結果に含まれる、前記シミュレーションの実行中に観測された事象を記録した数値データに基づいて生成される、請求項１に記載の最適化装置。
3. 入力される問題データと前記探索部が演算可能なビット数とに基づいて、前記ビット数以内の数の状態変数を用いた前記イジングモデルを生成する生成部、をさらに有する請求項１または２に記載の最適化装置。
4. 前記組合せ最適化問題は、複数の荷物が運び込まれる施設において、前記複数の荷物を運び込む複数の車両のそれぞれを複数の入荷品置き場の何れに割当てるかを決定する問題であり、
   前記探索部は、前記複数の車両のそれぞれを前記複数の入荷品置き場の何れかに割当てることによって生じる、前記施設内における前記複数の荷物を移動させる複数の移動体の総移動距離を表すコスト項と、前記問題の制約条件を表す第２の制約項とによって表される前記イジングモデルについての基底状態の探索を行い、
   前記実行部は、前記結果に含まれる前記複数の移動体による作業時間が、前記評価基準を満たす場合、前記仮最適解を前記最適解として出力する、
   請求項１乃至３の何れか一項に記載の最適化装置。
5. 最適化装置が有する探索部が、組合せ最適化問題を変換したイジングモデルについての基底状態の探索を行うことで、前記組合せ最適化問題の仮最適解を算出し、
   前記最適化装置が有する実行部が、前記仮最適解を用いたシミュレーションを実行し、前記シミュレーションの結果の評価基準を示す評価基準値に基づいて前記結果を評価し、前記結果が前記評価基準を満たす場合、前記仮最適解を最適解として出力し、
   前記最適化装置が有する更新部が、前記結果が前記評価基準を満たさない場合、前記イジングモデルに前記結果に基づいた第１の制約項を追加した更新イジングモデルを生成し、前記更新イジングモデルについての基底状態の探索を前記探索部に実行させる、
   最適化方法。
6. 組合せ最適化問題を変換したイジングモデルについての基底状態の探索を行うことで、前記組合せ最適化問題の仮最適解を算出し、
   前記仮最適解を用いたシミュレーションを実行し、
   前記シミュレーションの結果の評価基準を示す評価基準値に基づいて前記結果を評価し、
   前記結果が前記評価基準を満たす場合、前記仮最適解を最適解として出力し、
   前記結果が前記評価基準を満たさない場合、前記イジングモデルに前記結果に基づいた第１の制約項を追加した更新イジングモデルを生成し、前記更新イジングモデルについての基底状態の探索を実行する、
   処理をコンピュータに実行させる最適化プログラム。

Images