Note
Cette page a été générée à partir de tutorials/optimization/7_examples_vehicle_routing.ipynb.
Exécuter en mode interactif dans le ` IBM Quantum lab <https://quantum-computing.ibm.com/jupyter/tutorial/optimization/7_examples_vehicle_routing.ipynb>` _.
** Routage de véhicules * *¶
L’introduction¶
Logistics is a major industry, with some estimates valuing it at USD 8183 billion globally in 2015. Most service providers operate a number of vehicles (e.g., trucks and container ships), a number of depots, where the vehicles are based overnight, and serve a number of client locations with each vehicle during each day. There are many optimization and control problems that consider these parameters. Computationally, the key challenge is how to design routes from depots to a number of client locations and back to the depot, so as to minimize vehicle-miles traveled, time spent, or similar objective functions. In this notebook we formalize an idealized version of the problem and showcase its solution using the quantum approximate optimization approach of Farhi, Goldstone, and Gutmann (2014).
Le workflow global que nous allons détailler comprend:
établir les emplacements des clients à livrer : normalement ceux-ci sont disponibles avant le jour des livraisons dans une base de données. Dans notre cas d’utilisation, nous les générons aléatoirement.
calculer les distances point à point, les temps de trajet (ou équivalent). Dans notre cas, nous considérons la distance euclidienne, « à vol d’oiseau », ce qui est peut-être le plus simple possible.
calculer les routes réelles. Cette étape est exécutée deux fois, en fait. Tout d’abord, nous obtenons une valeur de référence par une exécution d’un solveur classique (IBM CPLEX) sur un ordinateur classique. Ensuite, nous avons une alternative, un algorithme hybride en partie sur l’ordinateur quantique.
visualisation des résultats, dans notre cas, il s’agit à nouveau d’un graphique simple.
Dans ce qui suit, nous expliquons d’abord le modèle utilisé ci-dessus, avant de procéder à l’installation des pré-requis et au chargement des données.
Le Modèle¶
En termes mathématiques, le problème du routage des véhicules (VRP - vehicle routing problem) est un problème de combinatoire où l’on recherche les meilleurs itinéraires d’un dépôt à un certain nombre de clients et retour, étant donné un certain nombre de véhicules disponibles. Il existe plusieurs formulations possibles qui étendent un certain nombre de formulations du problème du voyageur de commerce [Applegate et al, 2006]. Ici, nous présentons une formulation connue sous le nom de MTZ [Miller, Tucker, Zemlin, 1960].
Soit \(n\) s le nombre de clients (indexés par \(1, \dots, n\)), et \(K\) le nombre de véhicules disponibles. Soit \(x_{ij} = \ { 0,1 \ }\) la variable de décision binaire qui, si elle vaut \(1\), active le segment du noeud \(i\) au noeud \(j\). L’index du noeud passe de \(0\) à \(n\), où \(0\) est (par convention) le dépôt. Il existe deux fois plus de variables de décision distinctes que d’arêtes. Par exemple, dans un graphe entièrement connecté, il y a :math:”n(n+1) variables de décision binaires.
Si deux noeuds \(i\) et \(j\) sont liés de \(i\) à \(j\), nous écrivons \(i \sim j\). Nous notons également \(\delta (i) ^ +\) l’ensemble de noeuds avec lesquels le noeud \(i\) a un lien, c’est-à-dire \(j \in \delta (i) ^ +\) si et seulement si \(i \sim j\). De même, nous notons :math:` delta (i) ^-` l’ensemble de noeuds qui sont connectés à \(i\), dans le sens où \(j \in \delta (i) ^-\) si et seulement si \(j \sim i\).
En outre, nous considérons que les variables continues, pour tous les noeuds \(i = 1, \dots, n\), sont notées \(u_i\). Ces variables sont nécessaires dans la formulation MTZ du problème pour éliminer les sous-circuits entre clients.
Le VRP peut être formulé ainsi :
sous réserve de la contrainte de visite des dépots :
les contraintes de visite des dépôts :
et les contraintes d’élimination des sous-circuits :
En particulier, -La fonction de coût est une fonction linéaire des fonctions de coût et prend en compte les différents arcs avec un poids positif \(w_{ij}> 0\) (généralement il s’agit de la distance entre le noeud \(i\) et le noeud \(j\)) ; -Le premier ensemble de contraintes impose qu’un seul lien aille sur chaque client ; -Le deuxième ensemble de contraintes impose que, depuis et vers le dépôt, exactement \(K\) liens sont autorisés ; -Le troisième ensemble de contraintes applique les contraintes d’élimination des cous-circuits et appliquent les limites sur \(u_i\), avec \(Q > q_j> 0\), et \(Q, q_i \in \mathbb{R}\).
Solution classique¶
Nous pouvons résoudre le VRP classiquement, par exemple en utilisant CPLEX. Le CPLEX utilise une méthode de branchement et de coupe pour trouver une solution approximative de la VRP, qui, dans cette formulation, est un programme linéaire mixte entier (MILP). Pour des raisons de notation, nous comptons les variables de décision dans un vecteur comme
où \({\bf z} \in \{0,1\}^N\), avec \(N = n (n+1)\). Ainsi, la dimension du problème varie de façon quadratique avec le nombre de nœuds. Désignons la solution optimale de \({\bf z}^*\), et le coût optimal associé \(f^*\).
Solution Quantique¶
Here, we demonstrate an approach that combines classical and quantum computing steps, following the quantum approximate optimization approach of Farhi, Goldstone, and Gutmann (2014). In particular, we use the variational quantum eigensolver (VQE). We stress that given the use of limited depth of the quantum circuits employed (variational forms), it is hard to discuss the speed-up of the algorithm, as the solution obtained is heuristic in nature. At the same time, due to the nature and importance of the target problems, it is worth investigating heuristic approaches, which may be worthwhile for some problem classes.
D’après [5], l’algorithme peut être résumé comme suit: -Etapes de préparation: - Transformation du problème combinatoire en un problème d’optimisation polynomiale binaire avec des contraintes d’égalité uniquement ;- Mapping du problème résultant dans un hamiltonien Ising (\(H\)) pour les variables \({ \bf z }\) et la base \(Z\), via des méthodes de pénalité si nécessaire ; -Choix de la profondeur du circuit quantique \(m\). Notez que la profondeur peut être modifiée de façon adaptative. - Choix d’un ensemble de commandes \(\theta\) et construction d’une fonction d’essai \(\big | \psi (\boldsymbol\theta)\rangle\), à l’aide d’un circuit quantique fait de portes C-Phase et de rotations sur un seul qubit selon Y, paramétrées par les composants de :math:` boldsymboltheta `.
Etapes de l’algorithme:
Évaluer \(C( \boldsymbol\theta) = \langle\psi (\boldsymbol\theta) \big | H \big | \psi (\boldsymbol\theta)\rangle\) en échantillonant le résultat du circuit selon la base Z et en additionnant les valeurs d’espérance des différent termes d’Ising. En général, différents points de contrôle autour de :math:` boldsymboltheta doivent être estimés, en fonction de l’optimiseur classique choisi.
Utiliser un optimiseur classique pour choisir un nouvel ensemble de contrôles.
Continuer jusqu’à ce que \(C(\boldsymbol\theta)\) atteigne un minimum, assez proche de la solution \(\boldsymbol\theta^*\).
Utiliser le dernier \(\boldsymbol\thet\) pour générer un ensemble final d’échantillons à partir de la distribution \(\Big|\langle z_i\big|\psi(\boldsymbol\theta)\rangle\Big|^2\; a forall i\) afin d’obtenir la réponse.
Il existe de nombreux paramètres, notamment le choix de la fonction d’onde de l’essai. Ci-dessous, nous considérons:
où \(U_\mathrm{entangler}\) est un ensemble de portes C-Phase (portes d’intrication maximale), \(U_\mathrm{single}(\theta) = \prod_{i=1}^N Y(\theta_{i})\), où \(N\) est le nombre de qubits et \(m\) est la profondeur du circuit quantique.
Construction de l’hamiltonien Ising¶
A partir de :math:` VRP ` on peut construire une optimisation polynomiale binaire avec des contraintes d’égalité en considérant seulement les cas pour lesquels \(K=n-1\). Dans ces cas, les contraintes d’élimination de sous-circuits ne sont pas nécessaires et le problème ne dépend que de la variable :math:` { bf z } `. En particulier, nous pouvons écrire un lagrangien augmenté comme ceci :
où \(A\) est un paramètre suffisamment grand.
Du hamiltonien à la formulation QP¶
Par rapport au vecteur \({ \bf z }\), et pour un graphe complet (:math:` delta (i) ^ + = delta (i) ^-= { 0 ,1, dots, i-1,i + 1, dots, n} ), :math:`H peut être écrit comme suit.
Il s’agit de:
Où: premier terme:
Deuxième terme:
Troisième terme:
La formulation QP du Hamiltonien Ising est prête pour l’utilisation de VQE. Nous allons résoudre le QP en utilisant la pile d’optimisation disponible dans Qiskit.
Références¶
[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.
Initialisation¶
Tout d’abord, nous chargeons tous les packages dont nous avons besoin : - Python 3.6 ou supérieur; - CPLEX 12.8 ou supérieur sont requis pour les calculs classiques; - Le dernier niveau de Qiskit est requis pour les calculs quantiques.
[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
Nous initialisons ensuite les variables
[2]:
# Initialize the problem by defining the parameters
n = 3 # number of nodes + depot (n+1)
K = 2 # number of vehicles
Nous définissons une classe d’initialisation qui place aléatoirement les nœuds dans un plan en dimension 2 et calcule la distance entre eux.
[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()
Solution classique utilisant IBM ILOG CPLEX¶
Pour une solution classique, nous utilisons IBM ILOG CPLEX. CPLEX est capable de trouver la solution exacte pour ce problème. Nous définissons d’abord une classe ClassicalOptimizer qui code le problème d’une manière que CPLEX peut résoudre, puis nous instancions la classe et nous résolvons le problème.
[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')
Si vous avez CPLEX, la solution marque le dépôt avec une étoile et les routes sélectionnées pour les véhicules avec des flèches.
Solution quantique depuis le début¶
Pour la solution quantique, nous utilisons Qiskit.
Tout d’abord, nous construisons la solution en utilisant une classe QuantumOptimizer qui encode l’approche quantique pour résoudre le problème et ensuite nous l’instancions et le résolvons. Nous définissons les méthodes suivantes dans la classe : - binary_representation : pour encoder le problème \((M)\) en termes de QP (il s’agit essentiellement d’algèbre linéaire); - construct_problem : construit un problème d’optimisation QUBO en tant qu’instance de QuadraticProgram ; - solve_problem : résout le problème \((M)\) qui a été construit à l’étape précédente via MinimunEigenOptimizer en utilisant VQE avec des paramètres par défaut ;
[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
Étape 1¶
Instancions la classe d’optimiseur quantique avec les paramètres suivants : - l’instance ; - le nombre de nœuds et de véhicules n et K ;
[10]:
# Instantiate the quantum optimizer class with parameters:
quantum_optimizer = QuantumOptimizer(instance, n, K)
Étape 2¶
Encodons le problème en tant que formulation binaire (IH-QP).
Vérification : assurons-nous que la formulation binaire dans l’optimiseur quantique est correcte (c’est-à-dire qu’elle génère le même coût compte tenu de la même solution).
[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
Étape 3¶
Encodons le problème sous la forme d’une instance de QuadraticProgram.
[12]:
qp = quantum_optimizer.construct_problem(Q, g, c)
Étape 4¶
Résolvons le problème via MinimumEigenOptimizer depuis la pile d’optimisation. N.B. Selon le nombre de qubits, la simulation par vecteur d’état peut prendre un certain temps ; par exemple avec 12 qubits, elle prend plus de 12 heures. La journalisation est donc utile pour voir ce que fait le programme.
[13]:
quantum_solution, quantum_cost = quantum_optimizer.solve_problem(qp)
print(quantum_solution, quantum_cost)
[1. 1. 1. 0. 1. 0.] 132.11148115684045
Étape 5¶
Visualisons la solution
[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')
Les graphique présente le dépôt avec une étoile et les routes sélectionnées pour les véhicules avec des flèches. Notez que dans ce cas particulier, nous pouvons trouver la solution optimale de la formulation QP, qui se trouve coïncider avec la solution optimale de l’ILP.
Gardez à l’esprit que VQE est un travail heuristique sur la formulation QP du hamiltonien Ising. Pour les choix appropriés de A, l’optimisation locale de la formulation QP sera une solution réalisable pour l’ILP. Bien que pour des instances de petite taille, comme ci-dessus, nous pouvons trouver des solutions optimales de la formulation QP qui coïncident avec optimum de l’ILP, il est plus difficile de trouver des solutions optimales de l’ILP que de trouver des optimum locaux de la formulation QP, en général, ce qui à son tour est plus difficile que de trouver des solutions réalisables de l’ILP. Ainsi, dans le cadre du VQE, on peut fournir des garanties plus fortes, pour des formes variationnelles spécifiques (fonctions d’onde d’essai).
[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.
[ ]: