Nota
Esta página foi gerada, a partir de tutorials/optimization/7_examples_vehicle_routing.ipynb.
Execute interativamente no IBM Quantum lab.
Roteamento de Veículos¶
Introdução¶
A Logística é uma importante área industrial, pois - com base em algumas estimativas - foi avaliada globalmente em 8183 bilhões de dólares, no ano de 2015. E, a maioria dos prestadores de serviços opera vários veículos (por exemplo, caminhões e navios de contêineres), efetuam vários depósitos, onde os veículos passam a noite e servem vários locais de clientes, com cada veículo, durante cada dia. Existem muitos problemas de otimização e controle, que consideram tais parâmetros. Computacionalmente, o desafio chave é como projetar rotas, a partir de depósitos para vários locais de clientes e retornar aos mesmos, de modo a minimizar os quilômetros viajados pelos veículos, tempo gasto ou funções objetivas semelhantes. Neste artigo, formalizaremos uma versão idealizada do problema e mostraremos sua solução utilizando a abordagem de otimização aproximada quântica de Farhi, Goldstone e Gutman (2014).
O fluxo de trabalho geral, que demonstramos, compreende:
estabelecer os locais dos clientes. Normalmente, eles estariam disponíveis, antes do dia das entregas, a partir de um banco de dados. No nosso caso de uso, geramos so casos, aleatoriamente.
computar as distâncias dois a dois, tempos de viagem ou similares. No nosso caso, consideramos a distância euclidiana, em linha reta, que é, talvez, o mais simples possível.
computar as rotas reais. Esta etapa é executada duas vezes, na verdade. Primeiro, obtemos um valor de referência por uma execução de um solver clássico (IBM CPLEX) no computador clássico. Depois, executamos um algoritmo híbrido alternativo parcialmente no computador quântico.
visualização dos resultados. No nosso caso, este é novamente um plot simplista.
A seguir, explicamos primeiro o modelo, antes de prosseguir com a instalação dos pré-requisitos e o carregamento dos dados.
Modelo¶
Matematicamente falando, o problema de roteamento de veículos (VRP) é um problema combinatório, no qual as melhores rotas de um depósito para vários clientes e voltado para o depósito são procuradas, dado um certo número de veículos disponíveis. Há uma série de formulações possíveis, estendendo-se várias formulações do problema do caixeiro viajante [Applegate et al, 2006]. Aqui, apresentamos uma formulação conhecida como MTZ [Miller, Tucker, Zemlin, 1960].
Seja \(n\) o número de clientes (indexados como \(1,\dots,n\)) e \(K\) o número de veículos disponíveis. Seja \(x_{ij} = \{0,1\}\) a variável de decisão binária que, se for \(1\), ativa o segmento, a partir do nó \(i\) até o nó \(j\). O índice do nó vai de \(0\) até \(n\), sendo que \(0\) é (por convenção) o depósito. Há duas vezes mais variáveis de decisão distintas do que arestas. Por exemplo, em um gráfico completo, há \(n(n+1)\) variáveis de decisão binárias.
Se dois nós \(i\) e \(j\) possuem uma ligação de \(i\) para \(j\), escrevemos \(i \sim j\). Também denotamos por \(\delta(i)^+\) o conjunto de nós para os quais \(i\) possui uma ligação, ou seja, \(j \in \delta(i)^+\) se, e somente se, \(i \sim j\). Semelhantemente, denotamos por \(\delta(i)^-\) o conjunto de nós, que estão conectados a \(i\), no sentido que \(j \in \delta(i)^-\) se e somente se \(j \sim i\).
Além disto, consideramos variáveis contínuas para todos os nós \(i = 1,\dots, n\), denotadas \(u_i\). Estas variáveis são necessárias na formulação do MTZ do problema para eliminar subtours entre clientes.
O VRP pode ser formulado como:
sujeito à restrição de visita dos nós:
as restrições de visita dos depósitos:
e as restrições de eliminação de subtour:
Em particular, - A função de custo é linear nas funções de custo e pesa os diferentes arcos com base em um peso positivo \(w_{ij}>0\) (normalmente a distância entre o nó \(i\) e o nó \(j\)); - O primeiro conjunto de restrições impõe que de e para cada cliente, apenas um link é permitido; - O segundo conjunto de restrições impõe que de e para o depósito, exatamente \(K\) links são permitidos; - O terceiro conjunto de restrições impõe as restrições de eliminação de subtour e são limites em \(u_i\), com \(Q>q_j>0\) e \(Q,q_i \in \mathbb{R}\).
Solução clássica¶
Podemos resolver o VRP, classicamente, por exemplo, usando CPLEX. CPLEX usa um método branch-and-bound-and-cut para encontrar uma solução aproximada para o VRP, que nesta formulação é um programa linear inteiro misto (MILP). Para fins de notação, empacotamos as variáveis de decisão em um vetor como
onde \({\bf z} \in \{0,1\}^N\), com \(N = n (n+1)\). Portanto, a dimensão do problema escala, quadraticamente, com o número de nós. Vamos denotar a solução ideal por \({\bf z}^*\), e o custo ideal associado \(f^*\).
Solução quântica¶
Aqui, demonstraremos uma abordagem que combina etapas da computação clássica e quântica, seguindo a abordagem de otimização aproximada quântica de Farhi, Goldstone e Gutman (2014). Em particular, utilizaremos o algoritmo Variational Quantum Eigensolver (VQE). Salientamos que, dada a profundidade limitada dos circuitos quânticos empregados (formas variacionais), é difícil discutir o speedup do algoritmo, já que a solução obtida é heurística por natureza. Ao mesmo tempo, devido à natureza e à importância dos problemas alvo, vale a pena investigar abordagens heurísticas, que podem valer a pena para algumas classes de problemas.
Seguindo [5], o algoritmo pode ser resumido da seguinte forma: - Etapas de preparação: - Transformar o problema combinatório em um problema de otimização polinomial binária com restrições de igualdade apenas; - Mapear o problema resultante em um Hamiltoniano de Ising (\(H\)) para as variáveis \({\bf z}\) e base \(Z\), através de métodos de penalização se necessário; - Escolher a profundidade do circuito quântico \(m\). Observe que a profundidade pode ser modificada adaptativamente. - Escolher um conjunto de controles \(\theta\) e fazer uma função teste \(\big|\psi(\boldsymbol\theta)\rangle\), construída usando um circuito quântico feito de portas C-Fase e rotações Y de qubits únicos, parametrizada pelos componentes de \(\boldsymbol\theta\).
Etapas do algoritmo:
Avalie \(C(\boldsymbol\theta) = \langle\psi(\boldsymbol\theta)\big|H\big|\psi(\boldsymbol\theta)\rangle\) amostrando o resultado do circuito na base Z e adicionando os valores esperados dos termos de Ising individuais. Em geral, diferentes pontos de controle em torno de \(\boldsymbol\theta\) têm que ser estimados, dependendo do otimizador clássico escolhido.
Use um otimizador clássico para escolher um novo conjunto de controles.
Continue até \(C(\boldsymbol\theta)\) atingir um mínimo, próximo o suficiente da solução \(\boldsymbol\theta^*\).
Use o último \(\boldsymbol\theta\) para gerar um conjunto final de amostras da distribuição \(\Big|\langle z_i\big|\psi(\boldsymbol\theta)\rangle\Big|^2\;\forall i\) para obter a resposta.
Existem muitos parâmetros, notadamente a escolha da função da onda de teste. Abaixo, consideramos:
onde \(U_\mathrm{entangler} é uma coleção de portas C-Fase (portas totalmente entrelaçadas), e ::math:\), onde \(N\) é o número de qubits e \(m\) é a profundidade do circuito quântico.
Construa o Hamiltoniano de Ising¶
A partir do \(VRP\), pode-se construir uma otimização polinomial binária com restrições de igualdade, apenas considerando casos em que \(K=n-1\). Nestes casos, as restrições de eliminação de subtour não são necessárias e o problema é, apenas, na variável \({\bf z}\). Em particular, podemos escrever uma Lagrangiana aumentada como
onde \(A\) é um parâmetro grande o suficiente.
Do Hamiltoniano à formulação da QP¶
No vetor \({\bf z}\), e para um grafo completo (\(\delta(i)^+ = \delta(i)^- = \{0,1,\dots,i-1,i+1,\dots,n\}\)), \(H\) pode ser escrito da seguinte forma.
Isto é:
Onde: primeiro termo:
Segundo termo:
Terceiro termo:
A formulação QP do Hamiltoniano de Ising está pronta para uso do VQE. Vamos resolver o QP usando recursos de otimização disponíveis no Qiskit.
Referências¶
[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.
Inicialização¶
Primeiramente, carregamos todos os pacotes que precisamos: - Python 3.6 ou posterior é necessário; - CPLEX 12.8 ou posterior é necessário para os cálculos clássicos; - Qiskit mais atual é necessário para os cálculos quânticos.
[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
Inicializamos, então, as variáveis
[2]:
# Initialize the problem by defining the parameters
n = 3 # number of nodes + depot (n+1)
K = 2 # number of vehicles
Definimos uma classe inicializadora, que coloca os nós, aleatoriamente, em um plano 2-D e calcula a distância entre eles.
[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()
Solução clássica usando o IBM ILOG CPLEX¶
Para uma solução clássica, usamos o IBM ILOG CPLEX. O CPLEX é capaz de encontrar a solução exata deste problema. Definimos uma classe ClassicalOptimizer, que codifica o problema de uma forma que o CPLEX possa resolver e, em seguida, instanciamos a classe e a resolvemos.
[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')
Se você tem CPLEX, a solução mostra o depósito com uma estrela e as rotas selecionadas para os veículos com flechas.
Solução quântica do zero¶
Para a solução quântica, usamos o Qiskit.
Primeiro, derivamos a solução do zero, usando uma classe QuantumOptimizer que codifica a abordagem quântica para resolver o problema e, então, instanciamos e resolvemos. Definimos os seguintes métodos dentro da classe: - binary_representation : codifica o problema \((M)\) em termos do QP (que é basicamente álgebra linear); - construct_problem : constrói um problema de otimização QUBO como uma instância de QuadraticProgram; - solve_problem: resolve o problema \((M)\) construído na etapa anterior via MinimunEigenOptimizer, utilizando VQE com parâmetros padrões;
[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
Passo 1¶
Instancie a classe do otimizador quântico com parâmetros: - a instância; - o número de nós e veículos n e K;
[10]:
# Instantiate the quantum optimizer class with parameters:
quantum_optimizer = QuantumOptimizer(instance, n, K)
Passo 2¶
Codifique o problema como uma formulação binária (IH-QP).
Teste de sanidade: certifique-se de que a formulação binária no otimizador quântico esteja correta (ou seja, produza o mesmo custo dada a mesma solução).
[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
Passo 3¶
Codifique o problema como uma instância de QuadraticProgram.
[12]:
qp = quantum_optimizer.construct_problem(Q, g, c)
Passo 4¶
Resolva o problema, via MinimumEigenOptimizer do módulo de otimização. Nota: Dependendo do número de qubits, a simulação statevector pode demorar um pouco; por exemplo, com 12 qubits, leva mais de 12 horas. A criação de logs é útil para ver o que o programa está fazendo.
[13]:
quantum_solution, quantum_cost = quantum_optimizer.solve_problem(qp)
print(quantum_solution, quantum_cost)
[1. 1. 1. 0. 1. 0.] 132.11148115684045
Passo 5¶
Visualize a solução
[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')
Os gráficos apresentam o depósito com uma estrela e as rotas selecionadas para os veículos com flechas. Observe que, neste caso específico, podemos encontrar a solução ideal da formulação do QP, que por acaso coincide com a solução ideal do ILP.
Tenha em mente, no entanto, que o VQE é uma heurística trabalhando na formulação do QP do Hamiltoniano de Ising. Para escolhas adequadas de A, ótimos locais da formulação de QP serão soluções viáveis para o ILP. Enquanto que, para algumas pequenas instâncias, como acima, podemos encontrar soluções ideais da formulação do QP, que coincidem com a ótima do ILP, encontrar soluções ideais do ILP é mais difícil, do que encontrar ótimos locais da formulação do QP, em geral, o que por sua vez é mais difícil do que encontrar soluções factíveis do ILP. Mesmo dentro do VQE, pode-se fornecer garantias mais fortes, para formas variacionais específicas (funções de onda tentativa).
[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.
[ ]: