この記事の概要
Pilon et al., 2021, "Aircraft Loading Optimization - QUBO models under multiple constraints"を読み、その理解を深めるために内容をまとめたものです。
Aircraft Loading Optimization
航空機に搭載できる荷物には限りがあります。荷物コンテナをどのように航空機内に配置すれば、安全かつ運べる荷物重量が最大となる輸送が行えるかを最適化するのがこの最適化問題の目的です。
モデル化
モデル化にあたり以下のようなバイナリ変数を用います。
添字は利用可能コンテナの番号、は荷物を置くことができる位置を表します。
目的関数
コンテナは重量を持ちます。またコンテナサイズとしてT1(中くらいのサイズで1の広さを占有), T2(小さいサイズで1/2の広さを占有), T3(大きいサイズで2の広さを占有)の3種類を考えます。今、最大化したいのは運ぶ荷物の総重量としましょう。すると、目的関数はこれらを用いて
のように表現することができます。ここで
です。例えばあるT3コンテナの重さが2Mでこれを1つ搭載することは、重さMのT1コンテナを2つ搭載することと等価です。よってT3をこのように区別しています。
制約
Payload Limits (荷重制限)
- No-overlaps (コンテナにおける位置に対する制約)
1つの場所におけるコンテナは中/大サイズのコンテナを1つ、または小サイズのコンテナなら2つまでです。この制約を数理モデルで表現すると
ここで
は位置の占有を表現するための係数です。この制約をQUBOで定式化すると以下のようになります。
ここではスラック変数、そしてはその係数です。もしという同じ位置に不適切なコンテナ量が割り当てられたら、これが正の値となります。
- No-duplication (1つのコンテナは1つの位置にしか置けない制約)
小/中サイズのコンテナは1つの位置にしか置けません。大サイズのコンテナは隣り合う2つの位置に置けますが、これもある1つの位置に限られます。これを数理モデルで表現すると
となります。よってこの制約をQUBOで定式化すると以下のようになります。
- Contiguity for big containers (大サイズコンテナに対する制約)
(7)式において、大サイズのコンテナに関しては連続する位置を占有させる必要があります。よって(7)式を拡張する必要があります。もしコンテナが大サイズであるとすると
が成り立ちます。このコンテナを置くことが決まっているのであれば右辺は1、どこにも置かない場合には0です。よってこれをQUBOで定式化すると以下のようになります。
この制約式の形はすでに2次式になっているので、これ以上の変形の仕様がありません。しかしとの値を適切に取ることで、の値を負にならないようにすることができます。のときはとなるのは明かなので、ここからはのときにとなる場合を考えましょう。
のとき、 が最も小さくなる最悪の場合として考えられるのはです。またが最も小さくなるときを考えると
これが正の値となるには
であれば良いことがわかります。はのとき最大値2を取ることがわかるので、
であれば常にとなります。
- 最大積載量以下でなければならない制約
飛行機に搭載できる荷物の最大総重量をと書くと、この数理モデルは
となります。これをQUBOで定式化すると
です。
Center of Gravity Limits (重心制約)
荷物コンテナの配置が悪いと、航空機は姿勢制御に燃料を費やし最適な輸送となりません。よってここでは搭載した荷物の重心を計算し、それに対する制約を課します。ここではコンテナを配置する位置を
- ターゲットとする重心位置との一致
重心がこの位置に一致していることが望ましいとする値をとすると
ここでは荷物を搭載していないときの航空機の重量、そしては航空機の重心です。これをQUBOで定式化すると
となります。
- 上限と下限の制約
先ほどのターゲットにピタリと一致していなくとも、荷物を搭載したときの重心が]の範囲に収まってほしいというのを制約として課します。下限に対する制約は
そして上限に対する制約は
です。はよりも大きくなければいけません。(16), (17)式は(15)式に比べて強い制約でなければならないからです。この論文ではとしています。
Shear Limits (シアに対する制約)
重心位置が同じ2通りの置き方があるとして、一方は重心の周りに重たいものを配置する置き方、もう一方は重心から離れた位置に重たいものを配置する置き方だとします。後者の方法では航空機体に大きなシアが発生し、安全な航行ではなくなる可能性があります。よってここではシアの大きさに対する制約を考えましょう。
線形・対称な場合を考えるとシアの最大値は
と書かれます。ここでは適当な定数、は航空機の長さです。そしてにおける位置は
そしての位置に働く、左側からのシアは
同様にの位置に働く、右側からのシアは
わかりやすくするため、ここではが偶数の場合を考えます。
- 左側のシア
これをQUBOで定式化すると
- 右側のシア
これをQUBOで定式化すると
Total QUBO
これまでの制約と目的関数をまとめると、この最適化問題を解くための全ハミルトニアンは以下のようになります。
不等式制約の計算法
例えばのような不等式制約では、Log encodingスラック変数を用いて
のように等式制約に直してから最適化計算を行います。同様にのような不等式制約では
のようにして等式制約にします。
最適化計算の実装
この論文では以下の3つの手法を用いて最適化計算を実装し、数値実験を行っています。
- QBsolv: D-Waveが公開している、Tabu-Searchに基づいたメタヒューリスティック最適化アルゴリズム。古典ソルバーとして使用。
- D-Wave 2000Q: D-Waveが開発した量子アニーリングハードウェアD-Wave 2000Q_5のキメラグラフに問題を埋め込み、最適化計算を行う手法。量子ソルバーとして使用。
- Exact solver: 全ての状態を計算しその中から一番エネルギーの低い状態を探すという、力まかせ(brute-force)ソルバーによる計算手法。
数値実験結果
論文では3つの場合に対してExact solver, QBSolv, D-Wave 2000Q手法を用いて計算を行い、それぞれのベンチマークを取得しました。
- Case: PL (Only Payload Limits)
- Case: PL + CL (Payload + Center of Gravity Limits)
- Case: PL + CL + SL (Payload + Center of Gravity + Shear Limits)
それではベンチマーク結果を見ていきましょう。
古典ソルバー
Average Time-to-Solution
QBsolvでこの最適化問題を解いたときのAverage Time-to-Solution (TTS)を示します。制約が増えるに連れてスラック変数の数が増えるため、計算時間が増えていると解釈できます。
成功確率
以下は解の成功確率を示したものです。CaseがPL+CLの場合、PLが満たされている解が出た確率は98.0%、CLが満たされている解が出た確率は100%というように結果を見ます。
下図はPLとPL+CLのときの、重心の出現分布です。PLのみ(左パネル)では重心の位置が(黄色線)よりも小さな値になっている解が出現していますが、PL+CL(右パネル)ではそれが解消されています。
シア制約
以下の表はシアのエラー割合です。0 errorsは全ての位置がシア制約を満たした解が出現した割合、1 errorsはシア制約を満たせなかった位置が1つあるような解が出現した割合などです。
シア制約(SL)を導入することでシアの大きさが規定範囲内に収まっているような解が出現する割合が若干増えているように見えます。しかし、PL+CLとPL+CL+SLのerrors=1列を見ると、PL+CL+SLの方が値が大きくなっています。これは重心制約とシア制約が拮抗しているためと考えられますが、さらなる調査が必要としています。
量子ソルバー
以下の表は厳密計算(exact), 古典ソルバー(qbsolv), 量子ソルバー(dwave)でのベンチマーク比較です。
qbsolvは100%の確率で実行可能解を算出し、そのうち33%は最適解を導いています。dwaveは80%で実行可能解を算出し、最適解はわずか4%です。これはD-Wave 2000Qにノイズが多く存在するためと考えられます。
さらに下の表は、得られた解から算出された荷物の総重量の最大値と平均値です。qbsolvに比べてdwaveが少ない重量の解が現れていることがわかります。
結言
航空機に積む荷物の重量や重心・シアに関しての制約の定式化の方法と、それら制約を満たしながらも積載量を最大化する問題を解く方法をご紹介いたしました。
またこの問題において最適解を見つける速度はD-Wave量子アニーラー上で量子ソルバーを用いたときの方が高速ですが、安定して良い解を出せるのは古典ソルバーであることもわかりました。