Jij Tech Blog

Jij inc.の開発日記です

航空機荷物搭載最適化問題

Aircraft Loading Optimization

航空機に搭載できる荷物には限りがあります。荷物コンテナをどのように航空機内に配置すれば、安全かつ運べる荷物重量が最大となる輸送が行えるかを最適化するのがこの最適化問題の目的です。

モデル化

モデル化にあたり以下のようなバイナリ変数を用います。

 \displaystyle{
p_{i, j} = \left\{ \begin{array}{ll} 
 1 & (コンテナiを位置jに配置するとき) \\
 0 & (それ以外)
\end{array} \right. \tag{1}
}

添字 i=1, \dots, nは利用可能コンテナの番号、 j=1, \dots, Nは荷物を置くことができる位置を表します。

目的関数

コンテナ iは重量 m_i(>0)を持ちます。またコンテナサイズとしてT1(中くらいのサイズで1の広さを占有), T2(小さいサイズで1/2の広さを占有), T3(大きいサイズで2の広さを占有)の3種類を考えます。今、最大化したいのは運ぶ荷物の総重量としましょう。すると、目的関数はこれらを用いて

 \displaystyle{
f^\mathrm{obj}(\mathbf{p}) = - \sum_{i=1}^{n} \sum_{j=1}^{N} t_i m_i p_{i, j} \tag{2}
}

のように表現することができます。ここで

 \displaystyle{
t_i = \left\{ \begin{array}{ll}
1 & (\mathrm{T1}, \mathrm{T2}のとき) \\
\frac{1}{2} & (\mathrm{T3}のとき)
\end{array} \right. \tag{3}
}

です。例えばあるT3コンテナの重さが2Mでこれを1つ搭載することは、重さMのT1コンテナを2つ搭載することと等価です。よってT3をこのように区別しています。

制約

Payload Limits (荷重制限)
  • No-overlaps (コンテナにおける位置に対する制約)

1つの場所におけるコンテナは中/大サイズのコンテナを1つ、または小サイズのコンテナなら2つまでです。この制約を数理モデルで表現すると

 \displaystyle{
\sum_{i=1}^n d_i p_{i, j} \leq 1 \tag{4}
}

ここで

 \displaystyle{
d_i = \left\{ \begin{array}{ll}
1 & (\mathrm{T1, T3}のとき) \\
\frac{1}{2} & (\mathrm{T2}のとき)
\end{array} \right. \tag{5}
}

は位置の占有を表現するための係数です。この制約をQUBOで定式化すると以下のようになります。

 \displaystyle{
f^{\bar{O}}_j = P^{\bar{O}} \left(\sum_{i=1}^n d_i p_{i, j} + \sum_k c_k v_k -1 \right)^2 \tag{6}
}

ここで v_kはスラック変数、そして c_kはその係数です。もし jという同じ位置に不適切なコンテナ量が割り当てられたら、これが正の値となります。

  • No-duplication (1つのコンテナは1つの位置にしか置けない制約)

小/中サイズのコンテナは1つの位置にしか置けません。大サイズのコンテナは隣り合う2つの位置に置けますが、これもある1つの位置に限られます。これを数理モデルで表現すると

 \displaystyle{
t_i \sum_{j=1}^N p_{i, j} \leq 1 \tag{7}
}

となります。よってこの制約をQUBOで定式化すると以下のようになります。

 \displaystyle{
f_i^{\bar{D}} = P^{\bar{D}} \left(t_i \sum_{j=1}^N p_{i, j} + \sum_k c_k v_k -1 \right)^2 \tag{8}
}

  • Contiguity for big containers (大サイズコンテナに対する制約)

(7)式において、大サイズのコンテナに関しては連続する位置を占有させる必要があります。よって(7)式を拡張する必要があります。もしコンテナ iが大サイズであるとすると

 \displaystyle{
\sum_{j=1}^{N-1} p_{i,j } p_{i, j+1} = \frac{1}{2}\sum_{j=1}^N p_{i, j} \tag{9}
}

が成り立ちます。このコンテナを置くことが決まっているのであれば右辺は1、どこにも置かない場合には0です。よってこれをQUBOで定式化すると以下のようになります。

 \displaystyle{
f_i^{C} = P^{C} \left( \frac{1}{2} \sum_{j=1}^N p_{i, j} - \sum_{j=1}^{N-1} p_{i, j} p_{i, j+1} \right) \tag{10}
}

この制約式の形はすでに2次式になっているので、これ以上の変形の仕様がありません。しかし P^C P^{\bar{D}}の値を適切に取ることで、 f^C_i + f^{\bar{D}}_iの値を負にならないようにすることができます。 \sum_{j=1}^N p_{i, j} = 1のときは f^C_i > 0となるのは明かなので、ここからは \sum_{j=1}^N p_{i, j} >2のときに f_i^C>0となる場合を考えましょう。

 \sum_{j=1}^N p_{i, j} = mのとき、 f_i^C が最も小さくなる最悪の場合として考えられるのは \sum_{j=1}^{N-1} p_{i, j} p_{i, j+1} = m-1です。また f_i^{\bar{D}}が最も小さくなるときを考えると

 \displaystyle{
f_i^C + f_i^{\bar{D}} 
= P^C \left( \frac{1}{2} m - (m-1) \right) + P^{\bar{D}} \left(\frac{1}{2} m-1 \right)^2 
= \frac{1}{2} P^C (2-m) + \frac{1}{4} P^{\bar{D}} (2-m)^2 
}

これが正の値となるには

 \displaystyle{ 
P^{\bar{D}} > \frac{2(m-2)}{(m-2)^2} P^{C}
}

であれば良いことがわかります。 \frac{2(m-2)}{(m-2)^2} \ (m>2) m=3のとき最大値2を取ることがわかるので、

 \displaystyle{
P^{\bar{D}} > 2 P^C \tag{11}
}

であれば常に f^C_i + f^{\bar{D}}_i > 0となります。

  • 最大積載量以下でなければならない制約

飛行機に搭載できる荷物の最大総重量を W_pと書くと、この数理モデルは

 \displaystyle{ 
\sum_{i=1}^n \sum_{j=1}^N  t_i m_i p_{i, j} \leq W_p \tag{12}
}

となります。これをQUBOで定式化すると

 \displaystyle{
f^W = P^W \left( \sum_{i=1}^n \sum_{j=1}^N t_i m_i p_{i, j} + \sum_k c_k v_k -W_p\right)^2 \tag{13}
}

です。

Center of Gravity Limits (重心制約)

荷物コンテナの配置が悪いと、航空機は姿勢制御に燃料を費やし最適な輸送となりません。よってここでは搭載した荷物の重心を計算し、それに対する制約を課します。ここではコンテナを配置する位置を

 \displaystyle{
x_j = \frac{L}{N} \left( j - \frac{N}{2} \right) - \frac{L}{2N} \tag{14}
}

  • ターゲットとする重心位置との一致

重心がこの位置に一致していることが望ましいとする値を x_{cg}^tとすると

 \displaystyle{
\hat{x}_{cg} = \frac{\bar{x}_{cg}}{\underline{x}_{cg}} = \frac{\sum_{i=1}^n \sum_{j=1}^N t_i m_i p_{i, j} x_j + W_e x_{cg}^e}{\sum_{i=1}^n \sum_{j=1}^N t_i m_i p_{i, j} + W_e x_{cg}^e} = x_{cg}^t \ \Longrightarrow \ \bar{x}_{cg} = x_{cg}^t \underline{x}_{cg}
}

ここで W_eは荷物を搭載していないときの航空機の重量、そして x_{cg}^eは航空機の重心です。これをQUBOで定式化すると

 \displaystyle{
f^{C_t} = P^{C_t} \left( \bar{x}_{cg} - x_{cg}^t \underline{x}_{cg} \right)^2 \tag{15}
}

となります。

  • 上限と下限の制約

先ほどのターゲットにピタリと一致していなくとも、荷物を搭載したときの重心が [x_{cg}^{min}, x_{cg}^{max}]の範囲に収まってほしいというのを制約として課します。下限に対する制約は

 \displaystyle{
f^{C_l} = P^{C_l} \left( \bar{x}_{cg} - x_{cg}^{min} \underline{x}_{cg} - \sum_{k} c_k v_k \right)^2 \tag{16}
}

そして上限に対する制約は

 \displaystyle{
f^{C_u} = P^{C_u} \left( \bar{x}_{cg} - x_{cg}^{max} \underline{x}_{cg} - \sum_{k} c_k v_k \right)^2 \tag{17}
}

です。 P^{C_l}, P^{C_u} P^{C_t}よりも大きくなければいけません。(16), (17)式は(15)式に比べて強い制約でなければならないからです。この論文では P^{C_l} = P^{C_u} \simeq 10 P^{C_t}としています。

Shear Limits (シアに対する制約)

重心位置が同じ2通りの置き方があるとして、一方は重心の周りに重たいものを配置する置き方、もう一方は重心から離れた位置に重たいものを配置する置き方だとします。後者の方法では航空機体に大きなシアが発生し、安全な航行ではなくなる可能性があります。よってここではシアの大きさに対する制約を考えましょう。

線形・対称な場合を考えるとシアの最大値は

 \displaystyle{
S^{max}(x) = \left\{ \begin{array}{ll} 
S_0^{max} \frac{L + 2x}{L} & (x<0) \\
S_0^{max} \frac{L - 2x}{L} & (x>0)
\end{array} \right. \tag{18}
}

と書かれます。ここで S_0^{max}は適当な定数、 L > 0は航空機の長さです。そして u = 1, \dots, Nにおける位置は

 \displaystyle{
x_u = \frac{L}{N} \left( u- \frac{N}{2} \right) \tag{19}
}

そして uの位置に働く、左側からのシアは

 \displaystyle{
\hat{S}^l (u) = \sum_{i=1}^n \sum_{j=1}^u t_i m_i p_{i, j} \quad (u: x_u < 0) \tag{20}
}

同様に uの位置に働く、右側からのシアは

 \displaystyle{
\hat{S}^r (u) = \sum_{i=1}^n \sum_{j=u+1}^N t_i m_i p_{i, j} \quad (u: x_u > 0) \tag{21}
}

わかりやすくするため、ここでは Nが偶数の場合を考えます。

  • 左側のシア

 \displaystyle{
\hat{S}^l (u) \leq S^{max} (x_u) \quad (u = 1, \dots, N/2) \tag{22}
}

これをQUBOで定式化すると

 \displaystyle{
f_{x_u}^{S_l} = P^{S_l} \left( \hat{S}^l (u)+ \sum_k c_k v_k - S^{max} (x_u) \right)^2 \tag{23}
}

  • 右側のシア

 \displaystyle{
\hat{S}^r (u) \leq S^{max} (x_u) \quad (u = N/2, \dots, N-1) \tag{22}
}

これをQUBOで定式化すると

 \displaystyle{
f_{x_u}^{S_r} = P^{S_r} \left( \hat{S}^r (u)+ \sum_k c_k v_k - S^{max} (x_u) \right)^2 \tag{23}
}

Total QUBO

これまでの制約と目的関数をまとめると、この最適化問題を解くための全ハミルトニアンは以下のようになります。

 \displaystyle{
H^{tot} = f^{\mathrm{obj}} + f^W + \sum_{i=1}^n f_i^{\bar{D}} + \sum_{i:iがT3の場合}f_i^{C} \sum_{j=1}^N f_j^{\bar{O}} \\
 + f^{C_t} + f^{C_l} + f^{C_u} \\
 + \sum_{u=1}^{N/2} f_{x_u}^{S_l} + \sum_{u=N/2}^{N-1} f_{x_u}^{S_r} \tag{24}
}

不等式制約の計算法

例えば (\mathrm{left \ side}) \leq uのような不等式制約では、Log encodingスラック変数を用いて

 \displaystyle{
(\mathrm{left \ side}) + \sum_{k=0}^{r_u} 2^k v_k = u
}

のように等式制約に直してから最適化計算を行います。同様に (\mathrm{left \ side}) \geq lのような不等式制約では

 \displaystyle{
(\mathrm{left \ side}) - \sum_{k=0}^{r_l} 2^k v_k = l
}

のようにして等式制約にします。

最適化計算の実装

この論文では以下の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)を示します。制約が増えるに連れてスラック変数の数が増えるため、計算時間が増えていると解釈できます。

f:id:jijinc_nakamura:20210225041846p:plain
500回の試行における平均TTS.
成功確率

以下は解の成功確率を示したものです。CaseがPL+CLの場合、PLが満たされている解が出た確率は98.0%、CLが満たされている解が出た確率は100%というように結果を見ます。

f:id:jijinc_nakamura:20210225042121p:plain
モデル別、制約を満たしている成功確率。

下図はPLとPL+CLのときの、重心の出現分布です。PLのみ(左パネル)では重心の位置が x_{cg}^{min}(黄色線)よりも小さな値になっている解が出現していますが、PL+CL(右パネル)ではそれが解消されています。

f:id:jijinc_nakamura:20210225042329p:plain
重心位置の分布
シア制約

以下の表はシアのエラー割合です。0 errorsは全ての位置がシア制約を満たした解が出現した割合、1 errorsはシア制約を満たせなかった位置が1つあるような解が出現した割合などです。

f:id:jijinc_nakamura:20210225042500p:plain
シア制約に終えるエラー割合

シア制約(SL)を導入することでシアの大きさが規定範囲内に収まっているような解が出現する割合が若干増えているように見えます。しかし、PL+CLとPL+CL+SLのerrors=1列を見ると、PL+CL+SLの方が値が大きくなっています。これは重心制約とシア制約が拮抗しているためと考えられますが、さらなる調査が必要としています。

量子ソルバー

以下の表は厳密計算(exact), 古典ソルバー(qbsolv), 量子ソルバー(dwave)でのベンチマーク比較です。

f:id:jijinc_nakamura:20210225042639p:plain
3つのソルバーにおける、平均TTS, 実行可能解の出現率, 最適解の出現率の比較

qbsolvは100%の確率で実行可能解を算出し、そのうち33%は最適解を導いています。dwaveは80%で実行可能解を算出し、最適解はわずか4%です。これはD-Wave 2000Qにノイズが多く存在するためと考えられます。

さらに下の表は、得られた解から算出された荷物の総重量の最大値と平均値です。qbsolvに比べてdwaveが少ない重量の解が現れていることがわかります。

f:id:jijinc_nakamura:20210225042829p:plain
複数回試行したときの最大重量と平均重量の比較

結言

航空機に積む荷物の重量や重心・シアに関しての制約の定式化の方法と、それら制約を満たしながらも積載量を最大化する問題を解く方法をご紹介いたしました。
またこの問題において最適解を見つける速度はD-Wave量子アニーラー上で量子ソルバーを用いたときの方が高速ですが、安定して良い解を出せるのは古典ソルバーであることもわかりました。

文責

中村翔、株式会社 Jij
Sho K. NAKAMURA, Jij Inc.
憂いの篩 -Pensieve-

※ 疑問点や間違いなどございましたら、お気軽にブログ上でご質問ください。