注釈
当ページは tutorials/optimization/7_examples_vehicle_routing.ipynb から生成されました。
IBM Quantum lab でインタラクティブに実行します。
配車ルーティング¶
はじめに¶
ロジスティクスは主要産業であり、2015年には世界全体で8兆1830億米ドルになるという推計もあります。ほとんどのサービスプロバイダーは、多数の車両（トラックやコンテナ船など）、車両が夜間に拠点とする多数のデポを運営しています。 車両は夜間にはデポを拠点とし、毎日、各車両で多くのクライアントの場所にサービスを提供します。 これらのパラメータを考慮する最適化および制御の問題がたくさんあります。 計算上の重要な課題は、デポから多数のクライアントの場所に行ってまたデポに戻るルートをどのように設計すれば、車両の走行距離、費やした時間、またはそれに類する目的関数を最小にすることができるか、ということです。 このノートブックでは、理想化された形の問題を定式化し、Farhi, Goldstone, and Gutman（2014）の量子近似最適化アプローチを使用してその解決策を紹介します。
私たちが示す全体的なワークフローは次のとおりです。
クライアントの場所を確立します。 通常、これらはデータベースからの配信日の前に利用可能です。 私たちのユースケースでは、これらをランダムに生成します。
ペアワイズ距離、移動時間などを計算します。 私たちの場合、「カラスが飛ぶように」ユークリッド距離を考慮します。これはおそらく最も単純な距離です。
実際のルートを計算します。 このステップは、実際には2回実行されます。 最初に、古典的なコンピューターで古典的なソルバー（IBM CPLEX）を実行して参照値を取得します。 次に、量子コンピューターで部分的に代替のハイブリッドアルゴリズムを実行します。
結果の視覚化。 私たちの場合、これも単純なプロットです。
以下では、前提条件のインストールとデータのロードに進む前に、まずモデルについて説明します。
モデル¶
数学的に言えば、配車ルート問題（vehicle routing problem、VRP）は組み合わせ問題であり、利用可能な車両の数を考慮して、デポから多数のクライアントに戻り、デポに戻る最適なルートが求められます。 巡回セールスマン問題の多くの定式化を拡張して、可能な多くの定式化があります [Applegate et al, 2006] 。 ここでは、MTZ [Miller, Tucker, Zemlin, 1960] として知られる定式化を紹介します。
\(n\) をクライアントの数（ \(1,\dots,n\) としてインデックス付け）、 \(K\) を利用可能な車両の数とします。 \(x_{ij} = \{0,1\}\) をバイナリ決定変数とします。これが \(1\) の場合、ノード \(i\) からノード \(j\) へのセグメントをアクティブにします。 ノードインデックスは \(0\) から \(n\) まで実行されます。ここで、 \(0\) は（慣例により）デポです。 エッジの2倍の異なる決定変数があります。 たとえば、完全接続されたグラフには、 \(n(n+1)\) 個のバイナリ決定変数があります。
2つのノード \(i\) と \(j\) に \(i\) から \(j\) へのリンクがある場合、 \(i \sim j\) と記述します。 また、 \(i \sim j\) の場合に限り、 \(i\) がリンクしているノードのセットを \(\delta(i)^+\) で表します。つまり、\(j \in \delta(i)^+\) です。 同様に、 \(j \sim i\) の場合に限り、 \(j \in \delta(i)^-\) という意味で、iに接続されているノードのセットをδ（i）-で表します。
さらに、 \(u_i\) で表されるすべてのノード \(i = 1,\dots, n\) について、連続変数を検討します。 これらの変数は、クライアント間のサブツアーを排除するために、問題のMTZ定式化で必要です。
VRPは以下のように定式化できます:
ノード訪問制約の対象：
デポ訪問の制約：
およびサブツアー除去の制約：
特に、 - コスト関数はコスト関数で線形であり、正の重み \(w_{ij}>0\) （通常はノード \(i\) とノード \(j\) の間の距離）に基づいてさまざまなアーチの重みを付けます。 - 最初の一連の制約により、すべてのクライアントとの間で、1つのリンクのみが許可されます。 - 制約の2番目のセットは、デポとの間で正確に \(K\) 個のリンクが許可されることを強制します。 - 制約の3番目のセットは、サブツアー除去制約を適用し、 \(Q>q_j>0\) 、および \(Q,q_i \in \mathbb{R}\) の境界です。
古典的なソリューション¶
CPLEXを使用するなどして、VRPを古典的に解決できます。 CPLEXは、branch-and-bound-and-cut 法を使用して、VRPの近似解を見つけます。これは、この定式化では、混合整数線形計画（mixed-integer linear program、MILP）です。 表記のために、決定変数を1つのベクトルに次のようにパックします。
ここでは、 \({\bf z} \in \{0,1\}^N\) 、 \(N = n (n+1)\) です。 したがって、問題の次元はノードの数に比例して二次関数的にスケーリングします。 最適解を \({\bf z}^*\) で表し、関連する最適コストを \(f^*\) で表します。
量子ソリューション¶
ここでは、Farhi, Goldstone, and Gutman（2014）の量子近似最適化アプローチに従って、古典的な計算ステップと量子的な計算ステップを組み合わせたアプローチを示します。 特に、変分量子固有値ソルバー（variational quantum eigensolver、VQE）を使用します。 強調しておきたいのは、使用できる量子回路の深さが限られていること（変分形式）を考えると、アルゴリズムの高速化について議論するのは難しいということです。というのは、得られる解は本質的にヒューリスティックなものとなるからです。 同時に、対象となる問題の性質と重要性のため、ヒューリスティックなアプローチについて調べるのは価値のあることです。いくつかの問題クラスではこれが有用でしょう。
[5] に従い、アルゴリズムは次のように要約できます。 - 準備手順： - 組み合わせ問題を等式制約のみのバイナリ多項式最適化問題に変換します。 - 必要に応じてペナルティ法を使用して、結果の問題を変数 \({\bf z}\) と基底 \(Z\) のイジングハミルトニアン（ \(H\) ）にマッピングします。 - 量子回路の深さ \(m\) を選択します。 深さは適応的に変更できることに注意してください。 - コントロール \(\theta\) のセットを選択し、 \(\boldsymbol\theta\) の成分によってパラメータ化された制御位相ゲートと単一量子ビットY回転で構成される量子回路を使用して構築された試行関数 \(\big|\psi(\boldsymbol\theta)\rangle\) を作成します。
アルゴリズムの手順：
Zベースで回路の結果をサンプリングし、個々のIsing項の期待値を合計して、 \(C(\boldsymbol\theta) = \langle\psi(\boldsymbol\theta)\big|H\big|\psi(\boldsymbol\theta)\rangle\) を評価します。 一般に、選択した古典的なオプティマイザに応じて、 \(\boldsymbol\theta\) 周辺のさまざまな制御点を推定する必要があります。
従来のオプティマイザを使用して、新しいコントロールのセットを選択します。
\(C(\boldsymbol\theta)\) が最小値に達するまで続け、解 \(\boldsymbol\theta^*\) に十分に近づけます。
最後の \(\boldsymbol\theta\) を使用して、分布 \(\Big|\langle z_i\big|\psi(\boldsymbol\theta)\rangle\Big|^2\;\forall i\) からサンプルの最終セットを生成し、答えを取得します。
全体を通して多くのパラメータがあり、特に試行波動関数の選択があります。 以下では、次を検討します。
ここで、 \(U_\mathrm{entangler}\) は制御位相ゲート（完全にエンタングルしたゲート）の集合であり、 \(U_\mathrm{single}(\theta) = \prod_{i=1}^N Y(\theta_{i})\) です。ここで、 \(N\) は量子ビットの数、 \(m\) は量子回路の深さです。
イジングハミルトニアンを構築する¶
\(VRP\) から、 \(K=n-1\) の場合を考慮することによってのみ、等式制約を使用して二値多項式最適化を構築できます。 これらの場合、サブツアー除去制約は必要なく、問題は変数 \({\bf z}\) にのみあります。 特に、拡張ラグランジュは次のように書くことができます。
ここで、 \(A\) は十分に大きいパラメーターです。
ハミルトニアンからQP定式化へ¶
ベクトル \({\bf z}\) で、完全グラフ（ \(\delta(i)^+ = \delta(i)^- = \{0,1,\dots,i-1,i+1,\dots,n\}\) ）の場合、 \(H\) は次のように記述できます。 。
それは：
ここで：最初の項：
第二の項：
第三の項：
イジングハミルトニアンのQP定式化はVQEの使用の準備ができています。 Qiskitで利用可能な最適化スタックを使用してQPを解決します。
参考文献¶
[1] E. Farhi, J. Goldstone, S. Gutmann e-print arXiv 1411.4028, 2014
[3] C. E. Miller, E. W. Tucker, and R. A. Zemlin (1960). “Integer Programming Formulations and Travelling Salesman Problems”. J. ACM. 7: 326–329. doi:10.1145/321043.321046.
[4] D. L. Applegate, R. M. Bixby, V. Chvátal, and W. J. Cook (2006). The Traveling Salesman Problem. Princeton University Press, ISBN 978-0-691-12993-8.
初期化¶
まず、必要なすべてのパッケージをロードします。 - Python3.6以降が必要です。 - 従来の計算にはCPLEX12.8以上が必要です。 - 量子計算には最新のQiskitが必要です。
[1]:
# Load the packages that are required
import numpy as np
import operator
import matplotlib.pyplot as plt
import sys
if sys.version_info < (3, 6):
raise Exception('Please use Python version 3.6 or greater.')
try:
import cplex
from cplex.exceptions import CplexError
except:
print("Warning: Cplex not found.")
import math
# Qiskit packages
from qiskit import BasicAer
from qiskit.quantum_info import Pauli
from qiskit.aqua import QuantumInstance, aqua_globals
from qiskit.aqua.algorithms import NumPyMinimumEigensolver, VQE
from qiskit.circuit.library import TwoLocal
from qiskit.aqua.components.optimizers import SPSA
from qiskit.aqua.operators import WeightedPauliOperator
# setup aqua logging
import logging
from qiskit.aqua._logging import set_logging_config, build_logging_config
#set_logging_config(build_logging_config(logging.DEBUG)) # choose INFO, DEBUG to see the log
次に、変数を初期化します
[2]:
# Initialize the problem by defining the parameters
n = 3 # number of nodes + depot (n+1)
K = 2 # number of vehicles
ノードを2次元平面にランダムに配置し、ノード間の距離を計算する初期化子クラスを定義します。
[3]:
# Get the data
class Initializer():
def __init__(self, n):
self.n = n
def generate_instance(self):
n = self.n
# np.random.seed(33)
np.random.seed(1543)
xc = (np.random.rand(n) - 0.5) * 10
yc = (np.random.rand(n) - 0.5) * 10
instance = np.zeros([n, n])
for ii in range(0, n):
for jj in range(ii + 1, n):
instance[ii, jj] = (xc[ii] - xc[jj]) ** 2 + (yc[ii] - yc[jj]) ** 2
instance[jj, ii] = instance[ii, jj]
return xc, yc, instance
[4]:
# Initialize the problem by randomly generating the instance
initializer = Initializer(n)
xc,yc,instance = initializer.generate_instance()
IBM ILOGCPLEX を使用した従来のソリューション¶
従来のソリューションでは、IBM ILOGCPLEX を使用します。 CPLEX は、この問題の正確な解決策を見つけることができます。 まず、CPLEX が解決できる方法で問題をエンコードする ClassicalOptimizer クラスを定義し、次にクラスをインスタンス化して解決します。
[5]:
class ClassicalOptimizer:
def __init__(self, instance,n,K):
self.instance = instance
self.n = n # number of nodes
self.K = K # number of vehicles
def compute_allowed_combinations(self):
f = math.factorial
return f(self.n) / f(self.K) / f(self.n-self.K)
def cplex_solution(self):
# refactoring
instance = self.instance
n = self.n
K = self.K
my_obj = list(instance.reshape(1, n**2)[0])+[0. for x in range(0,n-1)]
my_ub = [1 for x in range(0,n**2+n-1)]
my_lb = [0 for x in range(0,n**2)] + [0.1 for x in range(0,n-1)]
my_ctype = "".join(['I' for x in range(0,n**2)]) + "".join(['C' for x in range(0,n-1)])
my_rhs = 2*([K] + [1 for x in range(0,n-1)]) + [1-0.1 for x in range(0,(n-1)**2-(n-1))] + [0 for x in range(0,n)]
my_sense = "".join(['E' for x in range(0,2*n)]) + "".join(['L' for x in range(0,(n-1)**2-(n-1))])+"".join(['E' for x in range(0,n)])
try:
my_prob = cplex.Cplex()
self.populatebyrow(my_prob,my_obj,my_ub,my_lb,my_ctype,my_sense,my_rhs)
my_prob.solve()
except CplexError as exc:
print(exc)
return
x = my_prob.solution.get_values()
x = np.array(x)
cost = my_prob.solution.get_objective_value()
return x,cost
def populatebyrow(self,prob,my_obj,my_ub,my_lb,my_ctype,my_sense,my_rhs):
n = self.n
prob.objective.set_sense(prob.objective.sense.minimize)
prob.variables.add(obj = my_obj, lb = my_lb, ub = my_ub, types = my_ctype)
prob.set_log_stream(None)
prob.set_error_stream(None)
prob.set_warning_stream(None)
prob.set_results_stream(None)
rows = []
for ii in range(0,n):
col = [x for x in range(0+n*ii,n+n*ii)]
coef = [1 for x in range(0,n)]
rows.append([col, coef])
for ii in range(0,n):
col = [x for x in range(0+ii,n**2,n)]
coef = [1 for x in range(0,n)]
rows.append([col, coef])
# Sub-tour elimination constraints:
for ii in range(0, n):
for jj in range(0,n):
if (ii != jj)and(ii*jj>0):
col = [ii+(jj*n), n**2+ii-1, n**2+jj-1]
coef = [1, 1, -1]
rows.append([col, coef])
for ii in range(0,n):
col = [(ii)*(n+1)]
coef = [1]
rows.append([col, coef])
prob.linear_constraints.add(lin_expr=rows, senses=my_sense, rhs=my_rhs)
[6]:
# Instantiate the classical optimizer class
classical_optimizer = ClassicalOptimizer(instance,n,K)
# Print number of feasible solutions
print('Number of feasible solutions = ' + str(classical_optimizer.compute_allowed_combinations()))
Number of feasible solutions = 3.0
[7]:
# Solve the problem in a classical fashion via CPLEX
x = None
z = None
try:
x,classical_cost = classical_optimizer.cplex_solution()
# Put the solution in the z variable
z = [x[ii] for ii in range(n**2) if ii//n != ii%n]
# Print the solution
print(z)
except:
print("CPLEX may be missing.")
[1.0, 1.0, 1.0, 0.0, 1.0, 0.0]
[8]:
# Visualize the solution
def visualize_solution(xc, yc, x, C, n, K, title_str):
plt.figure()
plt.scatter(xc, yc, s=200)
for i in range(len(xc)):
plt.annotate(i, (xc[i] + 0.15, yc[i]), size=16, color='r')
plt.plot(xc[0], yc[0], 'r*', ms=20)
plt.grid()
for ii in range(0, n ** 2):
if x[ii] > 0:
ix = ii // n
iy = ii % n
plt.arrow(xc[ix], yc[ix], xc[iy] - xc[ix], yc[iy] - yc[ix], length_includes_head=True, head_width=.25)
plt.title(title_str+' cost = ' + str(int(C * 100) / 100.))
plt.show()
if x is not None:
visualize_solution(xc, yc, x, classical_cost, n, K, 'Classical')
CPLEXを使用している場合、ソリューションには星付きのデポと矢印付きの車両用に選択されたルートが表示されます。
ゼロからの量子ソリューション¶
量子ソリューションには、Qiskitを使用します。
まず、問題を解決するために量子アプローチをエンコードするクラスQuantumOptimizerを使用して、ゼロからソリューションを導き出し、次にそれをインスタンス化して解決します。 クラス内で次のメソッドを定義します。– binary_representation ：問題 \((M)\) をQP項 (基本的に線形代数) にエンコードします。 – construct_problem ： QuadraticProgram のインスタンスとしてQUBO最適化問題を構築します。 – solve_problem ：デフォルトのパラメーターでVQEを使用することにより、 MinimunEigenOptimizer を介して前のステップで構築された問題 \((M)\) を解決します。
[9]:
from qiskit.optimization import QuadraticProgram
from qiskit.optimization.algorithms import MinimumEigenOptimizer
class QuantumOptimizer:
def __init__(self, instance, n, K):
self.instance = instance
self.n = n
self.K = K
def binary_representation(self,x_sol=0):
instance = self.instance
n = self.n
K = self.K
A = np.max(instance) * 100 # A parameter of cost function
# Determine the weights w
instance_vec = instance.reshape(n ** 2)
w_list = [instance_vec[x] for x in range(n ** 2) if instance_vec[x] > 0]
w = np.zeros(n * (n - 1))
for ii in range(len(w_list)):
w[ii] = w_list[ii]
# Some variables I will use
Id_n = np.eye(n)
Im_n_1 = np.ones([n - 1, n - 1])
Iv_n_1 = np.ones(n)
Iv_n_1[0] = 0
Iv_n = np.ones(n-1)
neg_Iv_n_1 = np.ones(n) - Iv_n_1
v = np.zeros([n, n*(n-1)])
for ii in range(n):
count = ii-1
for jj in range(n*(n-1)):
if jj//(n-1) == ii:
count = ii
if jj//(n-1) != ii and jj%(n-1) == count:
v[ii][jj] = 1.
vn = np.sum(v[1:], axis=0)
# Q defines the interactions between variables
Q = A*(np.kron(Id_n, Im_n_1) + np.dot(v.T, v))
# g defines the contribution from the individual variables
g = w - 2 * A * (np.kron(Iv_n_1,Iv_n) + vn.T) - \
2 * A * K * (np.kron(neg_Iv_n_1, Iv_n) + v[0].T)
# c is the constant offset
c = 2 * A * (n-1) + 2 * A * (K ** 2)
try:
max(x_sol)
# Evaluates the cost distance from a binary representation of a path
fun = lambda x: np.dot(np.around(x), np.dot(Q, np.around(x))) + np.dot(g, np.around(x)) + c
cost = fun(x_sol)
except:
cost = 0
return Q, g, c, cost
def construct_problem(self, Q, g, c) -> QuadraticProgram:
qp = QuadraticProgram()
for i in range(n * (n - 1)):
qp.binary_var(str(i))
qp.objective.quadratic = Q
qp.objective.linear = g
qp.objective.constant = c
return qp
def solve_problem(self, qp):
aqua_globals.random_seed = 10598
quantum_instance = QuantumInstance(BasicAer.get_backend('qasm_simulator'),
seed_simulator=aqua_globals.random_seed,
seed_transpiler=aqua_globals.random_seed)
vqe = VQE(quantum_instance=quantum_instance)
optimizer = MinimumEigenOptimizer(min_eigen_solver=vqe)
result = optimizer.solve(qp)
# compute cost of the obtained result
_,_,_,level = self.binary_representation(x_sol=result.x)
return result.x, level
ステップ 1¶
パラメータを使用して量子オプティマイザクラスをインスタンス化します。 - インスタンス。 - ノードと車両の数 n と K ;
[10]:
# Instantiate the quantum optimizer class with parameters:
quantum_optimizer = QuantumOptimizer(instance, n, K)
ステップ２¶
問題をバイナリ定式化（IH-QP）としてエンコードします。
サニティチェック：量子オプティマイザのバイナリ定式化が正しいことを確認します（つまり、同じソリューションで同じコストが得られることを確認します）。
[11]:
# Check if the binary representation is correct
try:
if z is not None:
Q, g, c, binary_cost = quantum_optimizer.binary_representation(x_sol = z)
print("Binary cost:", binary_cost, "classical cost:", classical_cost)
if np.abs(binary_cost - classical_cost) < 0.01:
print('Binary formulation is correct')
else: print('Error in the binary formulation')
else:
print('Could not verify the correctness, due to CPLEX solution being unavailable.')
Q, g, c, binary_cost = quantum_optimizer.binary_representation()
print("Binary cost:", binary_cost)
except NameError as e:
print("Warning: Please run the cells above first.")
print(e)
Binary cost: 132.11148115684045 classical cost: 132.1114811568365
Binary formulation is correct
ステップ３¶
QuadraticProgram のインスタンスとして問題をエンコードします。
[12]:
qp = quantum_optimizer.construct_problem(Q, g, c)
ステップ 4¶
最適化スタックから MinimumEigenOptimizer を介して問題を解決します。 N.B. 量子ビットの数によっては、状態ベクトルシミュレーションに時間がかかる場合があります。 たとえば、12量子ビットの場合、12時間以上かかります。 ロギングは、プログラムが何をしているかを確認するのに役立ちます。
[13]:
quantum_solution, quantum_cost = quantum_optimizer.solve_problem(qp)
print(quantum_solution, quantum_cost)
[1. 1. 1. 0. 1. 0.] 132.11148115684045
ステップ5¶
ソリューションを視覚化する
[14]:
# Put the solution in a way that is compatible with the classical variables
x_quantum = np.zeros(n**2)
kk = 0
for ii in range(n ** 2):
if ii // n != ii % n:
x_quantum[ii] = quantum_solution[kk]
kk += 1
# visualize the solution
visualize_solution(xc, yc, x_quantum, quantum_cost, n, K, 'Quantum')
# and visualize the classical for comparison
if x is not None:
visualize_solution(xc, yc, x, classical_cost, n, K, 'Classical')
プロットには、デポに星が表示され、車両用に選択されたルートが矢印で表示されます。 この特定のケースでは、QP定式化の最適解を見つけることができます。これはたまたまILPの最適解と一致します。
ただし、VQEはイジングハミルトニアンのQP定式化に取り組んでいるヒューリスティックであることに注意してください。 Aの適切な選択については、QP定式化の局所最適点がILPの実行可能解になります。 いくつかの小さな例では、上記のように、ILPの最適化と一致するQP公式の最適解を見つけることができますが、ILPの最適解を見つけることは、一般にQP公式の局所最適を見つけることよりも困難です。 ILPの実行可能な解決策を見つけるよりも難しい。 VQE内でも、特定の変分形式（試行波動関数）に対してより強力な保証を提供する場合があります。
[15]:
import qiskit.tools.jupyter
%qiskit_version_table
%qiskit_copyright
Version Information
| Qiskit Software | Version |
|---|---|
| Qiskit | 0.19.6 |
| Terra | 0.15.0 |
| Aer | 0.5.2 |
| Ignis | 0.4.0 |
| Aqua | 0.7.5 |
| IBM Q Provider | 0.7.2 |
| System information | |
| Python | 3.7.7 (default, May 6 2020, 11:45:54) [MSC v.1916 64 bit (AMD64)] |
| OS | Windows |
| CPUs | 4 |
| Memory (Gb) | 31.83771514892578 |
| Thu Aug 13 11:01:07 2020 GMT Daylight Time | |
This code is a part of Qiskit
© Copyright IBM 2017, 2020.
This code is licensed under the Apache License, Version 2.0. You may
obtain a copy of this license in the LICENSE.txt file in the root directory
of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.
Any modifications or derivative works of this code must retain this
copyright notice, and modified files need to carry a notice indicating
that they have been altered from the originals.
[ ]: