{ "cells": [ { "cell_type": "markdown", "metadata": { "collapsed": true }, "source": [ "# Vehicle Routing\n", "\n", "## The Introduction\n", "\n", "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). \n", "\n", "The overall workflow we demonstrate comprises:\n", "\n", "1. establish the client locations. Normally, these would be available ahead of the day of deliveries from a database. In our use case, we generate these randomly.\n", "\n", "2. compute the pair-wise distances, travel times, or similar. In our case, we consider the Euclidean distance, \"as the crow flies\", which is perhaps the simplest possible.\n", "\n", "3. compute the actual routes. This step is run twice, actually. First, we obtain a reference value by a run of a classical solver (IBM CPLEX) on the classical computer. Second, we run an alternative, hybrid algorithm partly on the quantum computer.\n", "\n", "4. visualization of the results. In our case, this is again a simplistic plot.\n", "\n", "In the following, we first explain the model, before we proceed with the installation of the pre-requisites and the data loading.\n", "\n", "## The Model \n", "\n", "Mathematically speaking, the vehicle routing problem (VRP) is a combinatorial problem, wherein the best routes from a depot to a number of clients and back to the depot are sought, given a number of available vehicles. There are a number of formulations possible, extending a number of formulations of the traveling salesman problem [Applegate et al, 2006]. Here, we present a formulation known as MTZ [Miller, Tucker, Zemlin, 1960]. \n", "\n", "Let $n$ be the number of clients (indexed as $1,\\dots,n$), and $K$ be the number of available vehicles. Let $x_{ij} = \\{0,1\\}$ be the binary decision variable which, if it is $1$, activates the segment from node $i$ to node $j$. The node index runs from $0$ to $n$, where $0$ is (by convention) the depot. There are twice as many distinct decision variables as edges. For example, in a fully connected graph, there are $n(n+1)$ binary decision variables. \n", "\n", "If two nodes $i$ and $j$ have a link from $i$ to $j$, we write $i \\sim j$. We also denote with $\\delta(i)^+$ the set of nodes to which $i$ has a link, i.e., $j \\in \\delta(i)^+$ if and only if $i \\sim j$. Similarly, we denote with \n", "$\\delta(i)^-$ the set of nodes which are connected to $i$, in the sense that $j \\in \\delta(i)^-$ if and only if $j \\sim i$. \n", "\n", "In addition, we consider continuous variables, for all nodes $i = 1,\\dots, n$, denoted $u_i$. These variables are needed in the MTZ formulation of the problem to eliminate sub-tours between clients. \n", "\n", "The VRP can be formulated as:\n", "\n", "$$\n", "(VRP) \\quad f = \\min_{\\{x_{ij}\\}_{i\\sim j}\\in \\{0,1\\}, \\{u_i\\}_{i=1,\\dots,n}\\in \\mathbb{R}} \\quad \\sum_{i \\sim j} w_{ij} x_{ij}\n", "$$\n", "\n", "subject to the node-visiting constraint:\n", "\n", "$$\n", "\\sum_{j \\in \\delta(i)^+} x_{ij} = 1, \\,\\sum_{j \\in \\delta(i)^-} x_{ji} = 1,\\, \\forall i \\in \\{1,\\dots,n\\},\n", "$$\n", "\n", "the depot-visiting constraints:\n", "\n", "$$\n", "\\sum_{i \\in \\delta(0)^+} x_{0i} = K, \\, \\sum_{j \\in \\delta(0)^+} x_{j0} = K,\n", "$$\n", "\n", "and the sub-tour elimination constraints:\n", "\n", "$$\n", "u_i - u_j + Q x_{ij} \\leq Q-q_j, \\, \\forall i \\sim j, \\,i ,j \\neq 0, \\quad q_i \\leq u_i \\leq Q,\\, \\forall i, i \\neq 0.\n", "$$\n", "\n", "In particular, \n", "\n", "- The cost function is linear in the cost functions and weighs the different arches based on a positive weight $w_{ij}>0$ (typically the distance between node $i$ and node $j$);\n", "- The first set of constraints enforce that from and to every client, only one link is allowed;\n", "- The second set of constraints enforce that from and to the depot, exactly $K$ links are allowed;\n", "- The third set of constraints enforce the sub-tour elimination constraints and are bounds on $u_i$, with $Q>q_j>0$, and $Q,q_i \\in \\mathbb{R}$.\n", "\n", "\n", "## Classical solution\n", "\n", "We can solve the VRP classically, e.g., by using CPLEX. CPLEX uses a branch-and-bound-and-cut method to find an approximate solution of the VRP, which, in this formulation, is a mixed-integer linear program (MILP). For the sake of notation, we pack the decision variables in one vector as\n", "\n", "$$\n", "{\\bf z} = [x_{01},x_{02},\\ldots,x_{10}, x_{12},\\ldots,x_{n(n-1)}]^T,\n", "$$\n", "\n", "wherein ${\\bf z} \\in \\{0,1\\}^N$, with $N = n (n+1)$. So the dimension of the problem scales quadratically with the number of nodes. Let us denote the optimal solution by ${\\bf z}^*$, and the associated optimal cost $f^*$. \n", "\n", "\n", "## Quantum solution\n", "\n", "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. \n", "\n", "Following [5], the algorithm can be summarized as follows:\n", "\n", "- Preparation steps: \n", "\t- Transform the combinatorial problem into a binary polynomial optimization problem with equality constraints only;\n", "\t- Map the resulting problem into an Ising Hamiltonian ($H$) for variables ${\\bf z}$ and basis $Z$, via penalty methods if necessary;\n", "\t- Choose the depth of the quantum circuit $m$. Note that the depth can be modified adaptively.\n", "\t- Choose a set of controls $\\theta$ and make a trial function $\\big|\\psi(\\boldsymbol\\theta)\\rangle$, built using a quantum circuit made of C-Phase gates and single-qubit Y rotations, parameterized by the components of $\\boldsymbol\\theta$.\n", "\n", "\n", "- Algorithm steps: \n", "\t- Evaluate $C(\\boldsymbol\\theta) = \\langle\\psi(\\boldsymbol\\theta)\\big|H\\big|\\psi(\\boldsymbol\\theta)\\rangle$ by sampling the outcome of the circuit in the Z-basis and adding the expectation values of the individual Ising terms together. In general, different control points around $\\boldsymbol\\theta$ have to be estimated, depending on the classical optimizer chosen.\n", "\t- Use a classical optimizer to choose a new set of controls.\n", "\t- Continue until $C(\\boldsymbol\\theta)$ reaches a minimum, close enough to the solution $\\boldsymbol\\theta^*$.\n", "\t- Use the last $\\boldsymbol\\theta$ to generate a final set of samples from the distribution $\\Big|\\langle z_i\\big|\\psi(\\boldsymbol\\theta)\\rangle\\Big|^2\\;\\forall i$ to obtain the answer.\n", "\n", "\n", "There are many parameters throughout, notably the choice of the trial wavefunction. Below, we consider:\n", "\n", "$$\n", "\\big|\\psi(\\theta)\\rangle = [U_\\mathrm{single}(\\boldsymbol\\theta) U_\\mathrm{entangler}]^m \\big|+\\rangle\n", "$$\n", "\n", "where $U_\\mathrm{entangler}$ is a collection of C-Phase gates (fully-entangling gates), and $U_\\mathrm{single}(\\theta) = \\prod_{i=1}^N Y(\\theta_{i})$, where $N$ is the number of qubits and $m$ is the depth of the quantum circuit. \n", "\n", "\n", "### Construct the Ising Hamiltonian\n", "\n", "From $VRP$ one can construct a binary polynomial optimization with equality constraints only by considering cases in which $K=n-1$. In these cases the sub-tour elimination constraints are not necessary and the problem is only on the variable ${\\bf z}$. In particular, we can write an augmented Lagrangian as\n", "\n", "$$\n", "(IH) \\quad H = \\sum_{i \\sim j} w_{ij} x_{ij} + A \\sum_{i \\in \\{1,\\dots,n\\}} \\Big(\\sum_{j \\in \\delta(i)^+} x_{ij} - 1\\Big)^2 + A \\sum_{i \\in \\{1,\\dots,n\\}}\\Big(\\sum_{j \\in \\delta(i)^-} x_{ji} - 1\\Big)^2 +A \\Big(\\sum_{i \\in \\delta(0)^+} x_{0i} - K\\Big)^2 + A\\Big(\\sum_{j \\in \\delta(0)^+} x_{j0} - K\\Big)^2\n", "$$\n", "\n", "where $A$ is a big enough parameter. \n", "\n", "### From Hamiltonian to QP formulation \n", "\n", "In the vector ${\\bf z}$, and for a complete graph ($\\delta(i)^+ = \\delta(i)^- = \\{0,1,\\dots,i-1,i+1,\\dots,n\\}$), $H$ can be written as follows.\n", "\n", "$$\n", "\\min_{{\\bf z}\\in \\{0,1\\}^{n(n+1)}} {\\bf w}^T {\\bf z} + A \\sum_{i \\in \\{1,\\dots,n\\}} \\Big({\\bf e}_i \\otimes {\\bf 1}_n^T {\\bf z} - 1\\Big)^2 + A \\sum_{i \\in \\{1,\\dots,n\\}}\\Big({\\bf v}_i^T {\\bf z} - 1\\Big)^2 + A \\Big(({\\bf e}_0 \\otimes {\\bf 1}_n)^T{\\bf z} - K\\Big)^2 + A\\Big({\\bf v}_0^T{\\bf z} - K\\Big)^2.\n", "$$\n", "\n", "That is:\n", "\n", "$$\n", "\\min_{\\bf z\\in \\{0,1\\}^{n(n+1)}} \\bf z^T {\\bf Q} \\bf z + {\\bf g}^T \\bf z + c,\n", "$$\n", "\n", "Where: first term:\n", "\n", "$$\n", "{\\bf Q} = A \\sum_{i \\in \\{0,1,\\dots,n\\}} \\Big[({\\bf e}_i \\otimes {\\bf 1}_n)({\\bf e}_i \\otimes {\\bf 1}_n)^T + {\\bf v}_i{\\bf v}_i^T \\Big] \n", "$$\n", "\n", "Second term:\n", "\n", "$$\n", "{\\bf g} = {\\bf w} -2 A \\sum_{i \\in \\{1,\\dots,n\\}} \\Big[({\\bf e}_i \\otimes {\\bf 1}_n) + {\\bf v}_i \\Big] -2 A K \\Big[({\\bf e}_0 \\otimes {\\bf 1}_n) + {\\bf v}_0 \\Big]\n", "$$\n", "\n", "Third term:\n", "\n", "$$\n", "c = 2An +2AK^2.\n", "$$\n", "\n", "The QP formulation of the Ising Hamiltonian is ready for the use of VQE. We will solve the QP using optimization stack available in Qiskit.\n", "\n", "\n", "\n", "## References\n", "\n", "[1] E. Farhi, J. Goldstone, S. Gutmann e-print arXiv 1411.4028, 2014\n", "\n", "[2] [Max-Cut and Traveling Salesman Problem](./06_examples_max_cut_and_tsp.ipynb)\n", "\n", "[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.\n", "\n", "[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." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Initialization\n", "\n", "First of all we load all the packages that we need: \n", " - Python 3.6 or greater is required;\n", " - CPLEX 12.8 or greater is required for the classical computations;\n", " - Latest Qiskit is required for the quantum computations." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Warning: Cplex not found.\n" ] } ], "source": [ "# Load the packages that are required\n", "import numpy as np\n", "import operator\n", "import matplotlib.pyplot as plt\n", "\n", "import sys\n", "\n", "if sys.version_info < (3, 6):\n", " raise Exception(\"Please use Python version 3.6 or greater.\")\n", "\n", "try:\n", " import cplex\n", " from cplex.exceptions import CplexError\n", "except:\n", " print(\"Warning: Cplex not found.\")\n", "import math\n", "\n", "# Qiskit packages\n", "from qiskit import BasicAer\n", "from qiskit.quantum_info import Pauli\n", "from qiskit.utils import QuantumInstance, algorithm_globals\n", "from qiskit.algorithms import NumPyMinimumEigensolver, VQE\n", "from qiskit.circuit.library import TwoLocal\n", "from qiskit.algorithms.optimizers import SPSA" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We then initialize the variables" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "# Initialize the problem by defining the parameters\n", "n = 3 # number of nodes + depot (n+1)\n", "K = 2 # number of vehicles" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We define an initializer class that randomly places the nodes in a 2-D plane and computes the distance between them. " ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "# Get the data\n", "class Initializer:\n", " def __init__(self, n):\n", " self.n = n\n", "\n", " def generate_instance(self):\n", "\n", " n = self.n\n", "\n", " # np.random.seed(33)\n", " np.random.seed(1543)\n", "\n", " xc = (np.random.rand(n) - 0.5) * 10\n", " yc = (np.random.rand(n) - 0.5) * 10\n", "\n", " instance = np.zeros([n, n])\n", " for ii in range(0, n):\n", " for jj in range(ii + 1, n):\n", " instance[ii, jj] = (xc[ii] - xc[jj]) ** 2 + (yc[ii] - yc[jj]) ** 2\n", " instance[jj, ii] = instance[ii, jj]\n", "\n", " return xc, yc, instance" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "# Initialize the problem by randomly generating the instance\n", "initializer = Initializer(n)\n", "xc, yc, instance = initializer.generate_instance()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Classical solution using IBM ILOG CPLEX\n", "\n", "For a classical solution, we use IBM ILOG CPLEX. CPLEX is able to find the exact solution of this problem. We first define a ClassicalOptimizer class that encodes the problem in a way that CPLEX can solve, and then instantiate the class and solve it. \n" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [], "source": [ "class ClassicalOptimizer:\n", " def __init__(self, instance, n, K):\n", "\n", " self.instance = instance\n", " self.n = n # number of nodes\n", " self.K = K # number of vehicles\n", "\n", " def compute_allowed_combinations(self):\n", " f = math.factorial\n", " return f(self.n) / f(self.K) / f(self.n - self.K)\n", "\n", " def cplex_solution(self):\n", "\n", " # refactoring\n", " instance = self.instance\n", " n = self.n\n", " K = self.K\n", "\n", " my_obj = list(instance.reshape(1, n ** 2)[0]) + [0.0 for x in range(0, n - 1)]\n", " my_ub = [1 for x in range(0, n ** 2 + n - 1)]\n", " my_lb = [0 for x in range(0, n ** 2)] + [0.1 for x in range(0, n - 1)]\n", " my_ctype = \"\".join([\"I\" for x in range(0, n ** 2)]) + \"\".join(\n", " [\"C\" for x in range(0, n - 1)]\n", " )\n", "\n", " my_rhs = (\n", " 2 * ([K] + [1 for x in range(0, n - 1)])\n", " + [1 - 0.1 for x in range(0, (n - 1) ** 2 - (n - 1))]\n", " + [0 for x in range(0, n)]\n", " )\n", " my_sense = (\n", " \"\".join([\"E\" for x in range(0, 2 * n)])\n", " + \"\".join([\"L\" for x in range(0, (n - 1) ** 2 - (n - 1))])\n", " + \"\".join([\"E\" for x in range(0, n)])\n", " )\n", "\n", " try:\n", " my_prob = cplex.Cplex()\n", " self.populatebyrow(my_prob, my_obj, my_ub, my_lb, my_ctype, my_sense, my_rhs)\n", "\n", " my_prob.solve()\n", "\n", " except CplexError as exc:\n", " print(exc)\n", " return\n", "\n", " x = my_prob.solution.get_values()\n", " x = np.array(x)\n", " cost = my_prob.solution.get_objective_value()\n", "\n", " return x, cost\n", "\n", " def populatebyrow(self, prob, my_obj, my_ub, my_lb, my_ctype, my_sense, my_rhs):\n", "\n", " n = self.n\n", "\n", " prob.objective.set_sense(prob.objective.sense.minimize)\n", " prob.variables.add(obj=my_obj, lb=my_lb, ub=my_ub, types=my_ctype)\n", "\n", " prob.set_log_stream(None)\n", " prob.set_error_stream(None)\n", " prob.set_warning_stream(None)\n", " prob.set_results_stream(None)\n", "\n", " rows = []\n", " for ii in range(0, n):\n", " col = [x for x in range(0 + n * ii, n + n * ii)]\n", " coef = [1 for x in range(0, n)]\n", " rows.append([col, coef])\n", "\n", " for ii in range(0, n):\n", " col = [x for x in range(0 + ii, n ** 2, n)]\n", " coef = [1 for x in range(0, n)]\n", "\n", " rows.append([col, coef])\n", "\n", " # Sub-tour elimination constraints:\n", " for ii in range(0, n):\n", " for jj in range(0, n):\n", " if (ii != jj) and (ii * jj > 0):\n", "\n", " col = [ii + (jj * n), n ** 2 + ii - 1, n ** 2 + jj - 1]\n", " coef = [1, 1, -1]\n", "\n", " rows.append([col, coef])\n", "\n", " for ii in range(0, n):\n", " col = [(ii) * (n + 1)]\n", " coef = [1]\n", " rows.append([col, coef])\n", "\n", " prob.linear_constraints.add(lin_expr=rows, senses=my_sense, rhs=my_rhs)" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Number of feasible solutions = 3.0\n" ] } ], "source": [ "# Instantiate the classical optimizer class\n", "classical_optimizer = ClassicalOptimizer(instance, n, K)\n", "\n", "# Print number of feasible solutions\n", "print(\"Number of feasible solutions = \" + str(classical_optimizer.compute_allowed_combinations()))" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "CPLEX may be missing.\n" ] } ], "source": [ "# Solve the problem in a classical fashion via CPLEX\n", "x = None\n", "z = None\n", "try:\n", " x, classical_cost = classical_optimizer.cplex_solution()\n", " # Put the solution in the z variable\n", " z = [x[ii] for ii in range(n ** 2) if ii // n != ii % n]\n", " # Print the solution\n", " print(z)\n", "except:\n", " print(\"CPLEX may be missing.\")" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [], "source": [ "# Visualize the solution\n", "def visualize_solution(xc, yc, x, C, n, K, title_str):\n", " plt.figure()\n", " plt.scatter(xc, yc, s=200)\n", " for i in range(len(xc)):\n", " plt.annotate(i, (xc[i] + 0.15, yc[i]), size=16, color=\"r\")\n", " plt.plot(xc[0], yc[0], \"r*\", ms=20)\n", "\n", " plt.grid()\n", "\n", " for ii in range(0, n ** 2):\n", "\n", " if x[ii] > 0:\n", " ix = ii // n\n", " iy = ii % n\n", " plt.arrow(\n", " xc[ix],\n", " yc[ix],\n", " xc[iy] - xc[ix],\n", " yc[iy] - yc[ix],\n", " length_includes_head=True,\n", " head_width=0.25,\n", " )\n", "\n", " plt.title(title_str + \" cost = \" + str(int(C * 100) / 100.0))\n", " plt.show()\n", "\n", "\n", "if x is not None:\n", " visualize_solution(xc, yc, x, classical_cost, n, K, \"Classical\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "If you have CPLEX, the solution shows the depot with a star and the selected routes for the vehicles with arrows. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Quantum solution from the ground up\n", "\n", "For the quantum solution, we use Qiskit. \n", "\n", "First, we derive the solution from the ground up, using a class QuantumOptimizer that encodes the quantum approach to solve the problem and then we instantiate it and solve it. We define the following methods inside the class:\n", "\n", "- binary_representation : encodes the problem $(M)$ into a QP terms (that's basically linear algebra);\n", "- construct_problem : constructs a QUBO optimization problem as an instance of QuadraticProgram;\n", "- solve_problem: solves the problem $(M)$ constructed at the previous step via MinimunEigenOptimizer by using VQE with default parameters;" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [], "source": [ "from qiskit_optimization import QuadraticProgram\n", "from qiskit_optimization.algorithms import MinimumEigenOptimizer\n", "\n", "\n", "class QuantumOptimizer:\n", " def __init__(self, instance, n, K):\n", "\n", " self.instance = instance\n", " self.n = n\n", " self.K = K\n", "\n", " def binary_representation(self, x_sol=0):\n", "\n", " instance = self.instance\n", " n = self.n\n", " K = self.K\n", "\n", " A = np.max(instance) * 100 # A parameter of cost function\n", "\n", " # Determine the weights w\n", " instance_vec = instance.reshape(n ** 2)\n", " w_list = [instance_vec[x] for x in range(n ** 2) if instance_vec[x] > 0]\n", " w = np.zeros(n * (n - 1))\n", " for ii in range(len(w_list)):\n", " w[ii] = w_list[ii]\n", "\n", " # Some variables I will use\n", " Id_n = np.eye(n)\n", " Im_n_1 = np.ones([n - 1, n - 1])\n", " Iv_n_1 = np.ones(n)\n", " Iv_n_1[0] = 0\n", " Iv_n = np.ones(n - 1)\n", " neg_Iv_n_1 = np.ones(n) - Iv_n_1\n", "\n", " v = np.zeros([n, n * (n - 1)])\n", " for ii in range(n):\n", " count = ii - 1\n", " for jj in range(n * (n - 1)):\n", "\n", " if jj // (n - 1) == ii:\n", " count = ii\n", "\n", " if jj // (n - 1) != ii and jj % (n - 1) == count:\n", " v[ii][jj] = 1.0\n", "\n", " vn = np.sum(v[1:], axis=0)\n", "\n", " # Q defines the interactions between variables\n", " Q = A * (np.kron(Id_n, Im_n_1) + np.dot(v.T, v))\n", "\n", " # g defines the contribution from the individual variables\n", " g = (\n", " w\n", " - 2 * A * (np.kron(Iv_n_1, Iv_n) + vn.T)\n", " - 2 * A * K * (np.kron(neg_Iv_n_1, Iv_n) + v[0].T)\n", " )\n", "\n", " # c is the constant offset\n", " c = 2 * A * (n - 1) + 2 * A * (K ** 2)\n", "\n", " try:\n", " max(x_sol)\n", " # Evaluates the cost distance from a binary representation of a path\n", " fun = (\n", " lambda x: np.dot(np.around(x), np.dot(Q, np.around(x)))\n", " + np.dot(g, np.around(x))\n", " + c\n", " )\n", " cost = fun(x_sol)\n", " except:\n", " cost = 0\n", "\n", " return Q, g, c, cost\n", "\n", " def construct_problem(self, Q, g, c) -> QuadraticProgram:\n", " qp = QuadraticProgram()\n", " for i in range(n * (n - 1)):\n", " qp.binary_var(str(i))\n", " qp.objective.quadratic = Q\n", " qp.objective.linear = g\n", " qp.objective.constant = c\n", " return qp\n", "\n", " def solve_problem(self, qp):\n", " algorithm_globals.random_seed = 10598\n", " quantum_instance = QuantumInstance(\n", " BasicAer.get_backend(\"qasm_simulator\"),\n", " seed_simulator=algorithm_globals.random_seed,\n", " seed_transpiler=algorithm_globals.random_seed,\n", " )\n", "\n", " vqe = VQE(quantum_instance=quantum_instance)\n", " optimizer = MinimumEigenOptimizer(min_eigen_solver=vqe)\n", " result = optimizer.solve(qp)\n", " # compute cost of the obtained result\n", " _, _, _, level = self.binary_representation(x_sol=result.x)\n", " return result.x, level" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Step 1\n", "\n", "Instantiate the quantum optimizer class with parameters: \n", "\n", "- the instance;\n", "- the number of nodes and vehicles n and K;" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [], "source": [ "# Instantiate the quantum optimizer class with parameters:\n", "quantum_optimizer = QuantumOptimizer(instance, n, K)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Step 2\n", "\n", "Encode the problem as a binary formulation (IH-QP).\n", "\n", "Sanity check: make sure that the binary formulation in the quantum optimizer is correct (i.e., yields the same cost given the same solution)." ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Could not verify the correctness, due to CPLEX solution being unavailable.\n", "Binary cost: 0\n" ] } ], "source": [ "# Check if the binary representation is correct\n", "try:\n", " if z is not None:\n", " Q, g, c, binary_cost = quantum_optimizer.binary_representation(x_sol=z)\n", " print(\"Binary cost:\", binary_cost, \"classical cost:\", classical_cost)\n", " if np.abs(binary_cost - classical_cost) < 0.01:\n", " print(\"Binary formulation is correct\")\n", " else:\n", " print(\"Error in the binary formulation\")\n", " else:\n", " print(\"Could not verify the correctness, due to CPLEX solution being unavailable.\")\n", " Q, g, c, binary_cost = quantum_optimizer.binary_representation()\n", " print(\"Binary cost:\", binary_cost)\n", "except NameError as e:\n", " print(\"Warning: Please run the cells above first.\")\n", " print(e)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Step 3\n", "\n", "Encode the problem as an instance of QuadraticProgram." ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [], "source": [ "qp = quantum_optimizer.construct_problem(Q, g, c)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Step 4\n", "\n", "Solve the problem via MinimumEigenOptimizer from the optimization stack. N.B. Depending on the number of qubits, the state-vector simulation can take a while; for example with 12 qubits, it takes more than 12 hours. Logging is useful to see what the program is doing." ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[1. 1. 1. 0. 1. 0.] 132.11148115684045\n" ] } ], "source": [ "quantum_solution, quantum_cost = quantum_optimizer.solve_problem(qp)\n", "\n", "print(quantum_solution, quantum_cost)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Step 5\n", "Visualize the solution" ] }, { "cell_type": "code", "execution_count": 14, "metadata": { "scrolled": true, "tags": [ "nbsphinx-thumbnail" ] }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAEICAYAAABRSj9aAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAy80lEQVR4nO3deXhU5fn/8fc9k5UsQAgkCCGRrbIpGIooVEBcgEpxgYpFfqhFKqho9VuqdqEuFS222KpIpbVoXSJSUbHgylZQZN+Rfd+3AAkhmczcvz9mwAgJScgkZya5X9d1rszMOfOczxyGe84855lzRFUxxhhTfbmcDmCMMaZyWaE3xphqzgq9McZUc1bojTGmmrNCb4wx1ZwVemOMqeas0BtjTDVnhd4YY6o5K/Q1lIjcKSKrROSkiOwTkfEiUruK1j1bRIZWxboqg4hMEpGnK6ntV0VkvYj4ROTOs+YNDMw7JiIHROR1EUkMzIsWkX+KyHYROSEiy0Wk93nW01ZEPhWRQyJyzq8mReR+EVksIvkiMinYr9NULSv0NZCIPAI8B/wKqA10BjKAz0Qk0sFoBlYAI4ClxcybD3RR1dpAUyACOP2BEwHsBLrh/zf9LTBZRDJKWI8HmAz8vIT5ewJtv1b+l2BCjqraVIMmIBHIAX561uPxwEFgSOD+JODpIvO7A7uK3H8U2AycANYCNxeZdycwD3geOApsBXoH5v0R8AKnAjlewv8ho0BEkTZmA0OLtDcfGAdkA1uAqwKP7wQOnM5dwmtOAv6Fv3gdBT4oMu8eYBNwBPgIuCjwuATWdwA4DqwC2gLD8BfJgkD+aZX07zQPuPM88+OBN4Dp51lmJXBrKetp7i8DJc5/Gpjk9PvWpopNtkdf81wFxADvF31QVXOA6cD1ZWxnM/Aj/HuPTwBvikjDIvOvANYDycCfgH+KiKjqb4D/Aferaryq3l/G9V2Bv3DVA94GsoAf4i9UdwAviUh8Cc/9N1ALaAM0wF/AEZFrgDHAT4GGwPZAu+DfDlcDLQOv8afAYVV9FXgL+FMgf9/iVigiK0Uku4RpfBlfc3HtdhWRY/g/YG8FXihhuZRA9jUXui5TfVihr3mSgUOqWljMvL1A/bI0oqrvqeoeVfWp6rvARqBTkUW2q+pEVfUCr+MvpCkVyL1VVf8VaO9dIA14UlXzVfUz/HvYzc9+UuDDpzdwr6oeVVWPqs4JzB4EvKaqS1U1H3gMuDLQ3eEBEoBLAFHVdaq6t6xhVfVSVa1TwjTiQjeCqs5Tf9dNY2AssK2Y1xyJ/8PodVX99kLXZaoPK/Q1zyEgWUQiipnXMDC/VCLy/wIH/LJFJBt/t0ZykUX2nb6hqicDN0va4y6L/UVu5wXaPfux4tpPA46o6tFi5l2Efy/+dM4c4DDQSFVn4u9Wehk4EDhImliB/EGlqruBT/juGwgAIuLC/w2mACjrtyVTzVmhr3m+BvKBW4o+GOj26I2/bxwgF393x2mpRZZNBybiLyT1VLUOsBp/v3ZZnD3KIzfwt9j1VdBOIElE6hQzbw+QfvqOiMTh7xraDaCqf1PVTKA1/m6QXwUWLfXc3iKyRkRySpgmVOwlnREBNCuyTgH+if+b062q6gnSekyYs0Jfw6jqMfx96i+KSC8RiQx0VUzGvzf/VmDR5UAfEUkSkVTgoSLNxOEvdgcBROQu/Hv0ZbUf/6iR05kO4i+ud4iIW0TupkgBq4hAd8sMYLyI1A283qsDs98B7hKR9iISDTwDfKOq20TkhyJyRaAbJBf/wWNfcflLWG+bQB9+cdO9JT1PRKJEJAb/h2akiMQE9tIRkUEi0iRwOx3/ge0vizz9FaAV0FdV886XT/xigKjA/ZjANjg9PyIw3w24A/OL+xZowoAV+hpIVf8EPI5/VMwJ/KNiagHXqurpvet/4x/qtw34DH+/+OnnrwX+jP/bwX6gHf5RMWX1V6C/iBwVkb8FHrsH/x7zYfwHTb+6kNdWgsH4+9y/xT+K5iEAVf0C+B3wH/zHJ5oBAwPPScT/reUo/u6dw/j7xMG/19w60G31QRBzgn9b5+E/aP5q4PbpD6bWwFcikot/e6/Hv91OF/5fAO2BfUW+PQwKzG8SuN8k0FZ6oO3TB2vzAu2d9tvAY4/iP9idF3jMhCFRtStM1XSBPfIn8Y/R3uF0HmNMcFmhNwCIyGDAo6pZpS5sjAkrVuiNMaaasz56Y4wJdSKNEXkRka8ROYmIUvLpLc7hyFH05ORkzcjIqNJ15ubmEhcXV6XrDBbL7gzL7pxwzl8Z2TNbtuSZLVtYUlhYcL1/CHRZf8Hu58R5FzIzM7WqzZo1q8rXGSyW3RmW3TnhnL9Ssnu9qqoKLFYYqqAKGWrnujHGGOcVen0cP+XB66vA8VBXxUq1/QDCGGOCLL/Qy/RVe3ll9mY2HsghwiUU+pSWDeK5t3sz+rRrSHSEu8ryWKE3xpggyvN4ueKPX+Lx+sgt8ALg8fr35tfvz+G3U1fzxEdref3uTlyWVqdKMlnXjTHGBMmKndlsOZhLdp7nTJE/W26Bl+w8DwNfXcCKndlVkssKvTHGBEF+oZchry3EV8bfJuV5/MvnFxb/gRBMFS70gZMdLRSRFYEz9j0RjGDGGBNOpq/ai8frK33BIjxeHzNW7St9wQoKxh59PnCNql6G/4RKvUSkcxDaNcaYsPHK7M0ldteUJLfAyyuzN1VSou9UuNAHRnnmBO5GBiY7r4Ixpsbw+pSNB3KKnecrLCB//5YSn7vhQE7Zhl5OmcKdUBfIDDzSG5H+iHQr7alBOdeNiLiBJfgv5fayqv66mGWG4b+wMikpKZlZWVV77qycnBzi4ytygSPnWHZnWHbnhFt+nypr955AVUmJhf15oN5CfHnH8eYcASAq9ZwrXQIgIrRumIBLzn/dnu49epQ0aw6q3c/33KCe1CxwFZ+pwAOqurqk5Tp27KiLFy8O2nrLYvbs2XTv3r1K1xkslt0Zlt054Zbf61Oa/2Y6Pp/yyCU5PPbcS+RtW47m+y/v4K6dQoNbf0dU/YxznisCm/7YB7er9Au0icgSVe1Y3nxBHXWjqtnALKBXMNs1xphQlptzguj1n7H3n/dRePwg0Wlt0Xz/pZIv+vkrxDRqRcHeDcU+t2WD+DIV+YoIxqib+qevxykiscB1+K/kY4wx1dqaNWsYOuxeGjZuwv4Ny0i6/l4ik5uQPed1QGl414tEJqcR1bAF+Xs3nvP8uCg3w7sX36UTTMH4ZWxD4PVAP70LmKyqHwehXWOMCTkej4cPPviA5/7yV75dv4HoNtdRd/BfiUhIBqDgwEbUc4rUIS8Q1eBiAKJSW5K7ZvY5bUW6XfRul1rpmStc6FV1JdAhCFmMMSZk7dmzh/GvTGD8319Faqfiat2Lel1/hbi/K6O7J/wcfvcrUu94nugiB1+jUi7Gc3gHWuhBIiIBiI108/rdnarknDd2rhtjjCmBqjJ37lzGjvsbX37xBXGtfkRM3+IPqu755wgKj+0nMqkR0Y0SvjfPFRlDRJ2GFBzcSlJ6KyLdrio9140VemOMOcuJEyd44403eP6FFzmam4+77Q3Uv2cirujiLyiy941f4jm0gwa3PY1ExRa7THTDltTN3ckfb/4pvdul2tkrjTHGKZ988im39O9PXm4Osc07UatDF1y1EvEc3oUrNgFXbCKu6FqI+Mey7Hv7UQr2bqTBgCeIzWgPFBbbblRqcy6LOcxNHRpV3YsJsEJvjDFFtGnTmuf/9ByHDh1i74FD7Nu/hwOHVnL48GGOHjnCiewj5OflEROfyMnjRwGo3eQSdMvXHN+5Em/TGzmxYinu0x8KsYm4YxOIatCUBd/8w5HXZIXeGGOKSEtLY8SIEeddxuPxkJKaykngjjvu4NZbb+XQoUMcPnyY5FoRXJ98jP0HNnFoh//D4djRo+Qez6ZBw4tQVaSUX8EGmxV6Y4wpp+YtWnD0yBFGjBjByy+//L15s2fP5p1/v37Oc3w+Hx6Pp8qLPNj56I0xplxatmzJju3bueuuu84p8ufjcrmIjo6uxGTnWbcjazXGmDDUrl07Nm7cyMCBA3nttdecjlNmVuiNMaYMOnbsyOrVq7npppt45513nI5TLlbojTGmFF26dGXJkiXccEMvpk6d6nSccrNCb4wx59Gz57V89dV8unXrziefzHA6zgWxQm+MMSX48Y9vZObML+ncuTOzZ89yOs4Fs0JvjDHFGDBgANOn/5f2HTrw9ddfOx2nQqzQG2PMWQYPHsyUKVNo1aoVy5YudTpOhVmhN8aYIoYNG8abb75Js2bNWLt2rdNxgsIKvTHGBIwcOZKJEyfSOC2NTZs2OR0naKzQG2MMMGrUKF588UVSUlLZuWOH03GCKhjXjE0TkVkislZE1ojIg8EIZowxVWX06NGMHTuWunWT2Ldvr9Nxgi4YJzUrBB5R1aUikgAsEZHPVbV6dG4ZY6q1MWPG8OSTT5KQkMiRI4edjlMpKrxHr6p7VXVp4PYJYB1Q9WfWN8aYcnrhhRd4/PHHiYmN5fjxY07HqTSiqsFrTCQDmAu0VdXjZ80bBgwDSElJyczKygraessiJyeH+Pj4Kl1nsFh2Z1h251RF/oMHD7Jjxw5EhMsvvzxo7VZm9h49eixR1Y7lfqKqBmUC4oElwC2lLZuZmalVbdasWVW+zmCx7M6w7M6p7PyvvfaaAup2u4PedmVmBxbrBdTnoIy6EZFI4D/AW6r6fjDaNMaYypCVlcXdd9+NiFBYWPz1XaubYIy6EeCfwDpV/UvFIxljTOWYOnUqt99+O0CNKfIQnHH0XYDBwDUisjww9QlCu8YYEzQzZszglltuAcDr9eJy1ZyfEVV4eKWqzgOq/iKIxhhTRjNnzqRPH//+Z00r8mC/jDXGVHPz58+nZ8+eQM0s8mCF3hhTjS1atIiuXbsC4PF4amSRByv0xphqauXKlXTq1AmA/Px8IiKCcSKA8GSF3hhT7axbt47LLrsMgLy8PKKiohxO5Cwr9MaYamXz5s20bt0agNzcXGJiYhxO5Dwr9MaYamPHjh00b94cgGPHjlGrVi2HE4UGK/TGmGphz549pKenA3D06FESExMdThQ6rNAbY8LegQMHaNTIf9LcgwcPUqdOHWcDhRgr9MaYsHb06FFSUlIA2Lt3L8nJyQ4nCj1W6I0xYev48eMkJSUBsH37dlJTUx1OFJqs0BtjwtLJkyepXbs2AJs2baJJkyYOJwpdVuiNMWHn1KlTxMXFAf4x882aNXM4UWizQm+MCSsFBQXExsYCsGLFCi655BKHE4U+K/TGmLBRWFhIdHQ04D+PzaWXXupwovBghd4YExZ8Ph+RkZEAzJs3j44dy3/p1JrKCr0xJuT5fD7cbjcAX375JV26dHE4UXixQm+MCWlFi/z06dO55pprHE4UfoJ1cfDXROSAiKwORnvGGAPfL/JTp06ld+/eDicKT8Hao58E9ApSW8YYA3DmHPJZWVncdNNNzoYJY0Ep9Ko6FzgSjLaMMQbAHRGBqjJp0iRuu+02p+OENVHV4DQkkgF8rKptS5g/DBgGkJKSkpmVlRWU9ZZVTk4O8fHxVbrOYLHszrDsztm/fz+7du0iPT097M5dU5nbvkePHktUtdzDjars2lqq+irwKkDHjh21e/fuVbVqAGbPnk1VrzNYLLszLLszYmNr8fTTT+F2u+nfv7/TccotFLd9zb2IojEm5CQkJHLqVB6NGjVi4MCBTsepNmx4pTEmJCQl1SMn5wSjR4+2s1AGWbCGV74DfA38QER2icjPg9GuMaZmSElJ4ejRI4waNYo//OEPTsepdoLSdaOqtwejHWNMzdM4LY0DBw4wcuRInnvuOafjVEvWdWOMcUyzZs3ZvWsX99xzD3/961+djlNtWaE3xjiiVevWbNmymcGDB/Pqq686Hadas0JvjKly7Tt04Nt16xgwYABvvPGG03GqPSv0xpgqdUXnzqxYvpwf//hGJk+e7HScGsEKvTGmynTv3p2F33xDz57X8vHH05yOU2NYoTfGVIkbbujFnDlz6NKlK1988bnTcWoUK/TGmEp3000389lnn9KxY0fmzfuf03FqHCv0xphKdfvtt/Phhx/Qrl07Fi1a5HScGskKvTGm0tx1111kZWXRsmVLVq5c6XScGssKvTGmUgwfPpxJkyaRnpHB+vXrnY5To1mhN8YE3cMPP8yECRNoeNFFbNu61ek4NZ4VemNMUD3++OOMGzeO5ORk9uze7XQcgxV6Y0wQPfnkk4wZM4batetw8OBBp+OYACv0xpigGDt2LKNHj6ZWXBzZ2UedjmOKsEJvjKmwl156iVGjRhEdHUNuTo7TccxZrNAbYyrkH//4Bw888AARERGcOpXndBxTDCv0xpgL9uabb3LPPffgcrnweDxOxzElCNalBHuJyHoR2SQijwajTWNMaJsyZQqDBw8GsCIf4ipc6EXEDbwM9AZaA7eLSOuKtmuMCV0fffQRAwYMAMDr9eJyWedAKAvGv04nYJOqblHVAiAL6BeEdo0xIejTTz+lXz//f3Er8uFBVLViDYj0B3qp6tDA/cHAFap6/1nLDQOGAaSkpGRmZWVVaL3llZOTQ3x8fJWuM1gsuzMse/Htnj6dQWZmZtDbL7oe2/bn6tGjxxJV7VjuJ6pqhSagP/CPIvcHAy+d7zmZmZla1WbNmlXl6wwWy+4My/59X331lQIKqMfjCXr7Rdm2Lx6wWC+gTgfjO9duIK3I/caBx4wx1cTSpUu56qqrAP+B14iICIcTmfIIRqFfBLQQkYtFJAoYCHwUhHaNMSFg9erVZ7pp8vPzrciHoQr/i6lqoYjcD3wKuIHXVHVNhZMZYxy3ceNG2rVrB0BeXh5RUVEOJzIXIigfzao6HZgejLaMMaFh69attGzZEoATJ04QExPjcCJzoWxclDHmHLt27aJp06YAHDt2LGxHwBg/K/TGmO/Zt28faWn+8RWHDx8mMTHR4USmoqzQG2POOHLkCA0bNgRg//79JCUlOZzIBIMVemMMANnZ2dSrVw+AnTt30qBBA4cTmWCxQm+MIScnh7p16wKwbds2Gjdu7HAiE0xW6I2p4U6ePElCQgIAGzZsID093eFEJtis0BtTgxUUFBAXFwf4fxjVokULhxOZymCF3pgaqqCggOjoaACWLVtGmzZtHE5kKosVemNqoMLCwjNFfsGCBbRv397ZQKZSWaE3pobx+XxERkYCMGfOHK644gqHE5nKZoXemBrE5/PhdrsB+Oyzz7j66qsdTmSqghV6Y2qIokV+2rRpXHfddQ4nMlXFCr0xNUDRIv/ee+9x4403OpzIVCUr9MbUABGBPvl///vf9O/f3+E0pqpZoTemmouMjER9PiZOnMgdd9zhdBzjACv0xlRj0TExFBYW8uKLLzJ06FCn4xiHWKE3ppqKi4+nID+f559/nvvvv9/pOMZBFSr0IjJARNaIiE9EOgYrlDGmYpYvX87J3FyeeuopHnnkEafjGIdVdI9+NXALMDcIWYwxQZCcXB+v18vjjz/Ob3/7W6fjmBBQoWvGquo6ABEJThpjTIU0vOgiDh8+REpKiu3JmzOsj96YaiI9I4N9e/cyYsQIO5+8+R5R1fMvIPIFkFrMrN+o6oeBZWYD/6eqi8/TzjBgGEBKSkpmVlbWhWa+IDk5OWF7gWPL7oxwyr569Wry8/NJTk4mPT09rLIXJ5zzV2b2Hj16LFHV8h8PVdUKT8BsoGNZl8/MzNSqNmvWrCpfZ7BYdmeES/a2bdsqoLfffvuZx8Ile0nCOX9lZgcW6wXUaOu6MSaMdezYkdWrV3PzzTfz9ttvOx3HhKiKDq+8WUR2AVcC/xWRT4MTyxhTmi5durJkyRJuuKEX77//vtNxTAir6KibqcDUIGUxxpRRz57X8tVX8+nevTuffDLD6TgmxFnXjTFh5sc/vpGZM7+kc+crmTVrltNxTBiwQm9MGBkwYADTp/+X9u3b8/XXXzkdx4QJK/TGhInBgwczZcoUWrVuzbJly5yOY8KIFXpjwsCwYcN48803adqsGWvXrHE6jgkzVuiNCXEjR45k4sSJpKU1YfOmTU7HMWHICr0xIWzUqFG8+OKLpKSksmPHdqfjmDBlhd6YEDV69GjGjh1L3aQk9u3b63QcE8as0BsTgsaMGcOTTz5JfEICRw4fdjqOCXNW6I0JMS+88AKPP/44sbGxnDh+3Ok4phqwQm9MCHnllVf45S9/SWRUFCdPnnQ6jqkmrNAbEyL+9a9/MWLECNxuNwX5+U7HMdWIFXpjQkBWVhZ33303IkJhYaHTcUw1Y4XeGIdNnTqV22+/HcCKvKkUVuiNcdCMGTO45ZZbAPB6vbhc9l/SBJ+9q4xxyMyZM+nTpw9gRd5ULntnGeOA+fPn07NnT8CKvKl89u4ypootWrSIrl27AuDxeKzIm0pn7zBjqtDKlSvp1KkTAPn5+UREVOgib8aUSUWvGTtWRL4VkZUiMlVE6gQplzHVzrp167jssssAyMvLIyoqyuFEpqao6B7950BbVb0U2AA8VvFIxlQ/mzdvpnXr1gDk5uYSExPjcCJTk1So0KvqZ6p6euDvAqBxxSMZU73s2LGD5s2bA3Ds2DFq1arlcCJT04iqBqchkWnAu6r6ZgnzhwHDAFJSUjKzsrKCst6yysnJIT4+vkrXGSyW3RnByO7xeFi5ciUA7du3x+12ByNaqcJ5u0N456/M7D169Fiiqh3L/URVPe8EfAGsLmbqV2SZ3wBTCXxwlDZlZmZqVZs1a1aVrzNYLLszKpp9//79CiigBw8eDE6oMgrn7a4a3vkrMzuwWMtQY8+eSj3kr6rXnm++iNwJ3Aj0DAQxpsbLzs4mJSUFgL1795KcnOxwIlOTVWhsl4j0AkYB3VTVzqlqDHD8+HHq1q0LwPbt20lNTXU4kanpKjrq5iUgAfhcRJaLyIQgZDImbJ08eZLatWsDsGnTJpo0aeJwImMquEevqs2DFcSYcHfq1Cni4uIA+Pbbb2nWrJnDiYzxs1/GGhMEBQUFxMbGArBixQp+8IMfOJzImO9YoTemggoLC4mOjgb857G59NJLHU5kzPfVrEI/ZQrceiukp0NsLPzgB/DYY3DihNPJTJjy+XxERkYCMG/ePDp2LP8QZ2MqW80q9M8/D243PPMMfPIJDB8Or7wC110HPp/T6UyY8fl8Z34A9eWXX9KlSxeHExlTvJp16rxp06B+/e/ud+sGSUkwZAjMng3XXONYNBNeihb56dOnc429d0wIC6s9+kKvj+OnPHh9F/i7rKJF/rQf/tD/d/fuCw9mapSiRX7q1Kn07t3b4UTGnF/I79HnF3qZvmovr8zezMYDOUS4hEKf0rJBPPd2b0afdg2JjqjA+UPmzPH/bdUqOIFNtXf6HPJZWVncdNNNzoYxpgxCutAv35nNna8txOP1kVvgBcDj9e/Nr9+fw2+nruaJj9by+t2duCytTvlXsHs3/P73cO21YAfRTBm4IyJQVSZNmsRtt93mdBxjyiRku25W7Mzm9lcXkJ3nOVPkz5Zb4CU7z8PAVxewYmd2+VaQkwP9+kFEBPzrXxUPbKq9qKgofF4vEyZMYMiQIU7HMabMQrLQ5xd6GfLaQvI8xRf4s+V5/MvnF5ZtefLyoG9f2LIFPv0UGttp9M35xcbWwuPxMG7cOH7xi184HceYcgnJQj991V483vINd/R4fcxYta8MC3qgf39YvBimT4d27S4wpakpEhISOXUqjzFjxvDQQw85HceYcgvJQv/K7M0ldteUJLfAyyuzN51/IZ8PBg2CmTPhgw+gc+cLD2lqhKSkeuTknGD06NE8+uijTscx5oKE3MFYr0/ZeCCn2Hk+TwE7/3ILEfXSqH3VQGKatCMiPunM/A0HcvD6FLdLim/8vvvgvffgN7+BuDhYsOC7eY0bWxeO+Z4VK1Zw9OgRRo0axR/+8Aen4xhzwUKu0OcWFBLhkjOja77H5cIVn0Th4Z0cnjYWxH8/Nq0t0Y1aEdekNcfz8qkbV8KFl2fM8P/94x/9U1GjR4P9ZzYBjdPS+OVDDzFy5Eiee+45p+MYUyEhV+jjoiIoLOEHUS53BI2GTmDf24/ijqvLqa1L8J04RO7a2XhP5XB8yTQuynqM9pd35Loe3fhR1y507tyZhIQEfwPbtlXdCzFhq1mz5uzetYvk5GQeeeQRp+MYU2Eh10fvdgktGpR8YV1XdC0a9B+N5/AOkvv+HykDnwHg1JbFFB7ZRUzmrWyq35WXPl/Lz4b/iuQGqTRv1Y7t23fw9ttvs337duyKh6YkrVq3ZsuWzQwePJj09HSn4xgTFCFX6AGGd29GXFTJv3aNSKhHg/6jOfLlRHC5SP/1xzQa/i/cCclkz32dQ+8/Tf6xA9S66Q+k3v8WuR3v5Gi+8tCzE2h9WSb1UhrSp98tjBs3jkWLFuHxeKrw1ZlQ1b5DB75dt44BAwbwxhtvOB3HmKCp6DVjnwL6AT7gAHCnqu6paKg+7RryxEdrgZJH3kTVzyC57684+MGzpP5sDJH10mg8YhJaWMCh/47j5No57Fg7h4jaKaQM+hPuuNrU6vNrYlUpzN7Hwt3rWPj2F/j+Mp6Th/bQ5rIOXNv9R1z9o65ceeWVJCUllbhuU/1c0bkzK5Yv58c/vpHJkyc7HceYoKroHv1YVb1UVdsDHwO/r3gkiI5w8/rdnYiNPP85bGIz2lO3+50ceO8PeHOPAiARUdTv92uajJpG3Z73UHhsP7vHD6Fg3yZO7VyNiBBZtyHxba8hvucIEge9QINfvMauJtczcd427nx4NA0bp9GkWUvuuPNuJk2aZHv81Vz37t1Z+M03XHvttXz88TSn4xgTdBW9ZuzxInfjgKB1fl+WVoesYZ0Zcta5boqKi3JTp1Nv+jWN4P1pz5Bw61O4Iv0jbkSExI79SOzYj1M7VwOw/23/OOi619xDQsefIOIfhilRsUTUSUUL8/G43CRG1WLP+oW88+83WPD11/Tv3//MxSVM9XLDDb2YM2cOXbt25fPPP3c6jjGVQip6YFJE/gj8P+AY0ENVD5aw3DBgGEBKSkpmVlZWmdpX4Fieh4Mn8jnl8SIiqCoxkW7qJ0RTOzYSAbZs3cqJkwW46qQW205KLOzLKcRzZBd4C0+nwh0Vhc/jwR0RQWxsDLVq1aJWbCyxsbHExMSc+TBwUk5ODvHxJR+gDmWhnH3z5s1kZ2dTKy6OVpdccs78UM5emnDODuGdvzKz9+jRY4mqlv8MjKp63gn4AlhdzNTvrOUeA54orT1VJTMzUy9Eodenx/IKtNDrO2defn6+du7aTet16qfpv/5YGz+YpSk/e1aTrh+h9X54o77w4ssaE5+oiXWT9fLOXTStSRPF/zmiIqILFiy4oExVYdasWU5HuGChmn3gwIEKaLt27UpcJlSzl0U4Z1cN7/yVmR1YrGWosWdPpXbdqOq1ZfzMeAuYDowu30dN2bldQmJM8V0oUVFRzJj2AZd3upL9E4bg8+TTrOUP6ND+Mjr2uZam6Wls37yRBg0anHmOz+dj4MCBvPfee3QOnA7hz3/+Mw8//HBlvQQTAu666y6ysrJo2bIlK1eudDqOMZWuQgdjRaRFkbv9gG8rFqdi6tSpw4olC1mzfDEnc46zetli/v2vf/Lggw+SkJDwvSIP4HK5mDx5MqrKK6+8AsAjjzyCiHDjjTfis+vIVjvDhw9n0qRJpGdksH79eqfjGFMlKjrq5lkRWS0iK4HrgQeDkKlCEhISyMjIKHff+r333ouqsmzZMiIjI/nvf/+L2+3moosu4siRI5WU1lSlhx9+mAkTJtDwoovYtnWr03GMqTIVKvSqequqtlX/EMu+qhr2F15t3749BQUFnDhxgmbNmrF3717q1auH2+1m/vz5TsczF+jxxx9n3LhxJCcns8euD2xqmJD8ZWwoiI+PZ9OmTagqP/vZz/D5fHTt2hUR4dlnn3U6nimHJ598kjFjxlC7dh0OHix2UJgx1ZoV+jJ46623UFUmTpwIwGOPPYaI0KtXL+vHD3Fjx45l9OjRxMXFkZ191Ok4xjjCCn05DB06FFVl1apVREVF8+mnn+J2u0lJTeXQoUNOxzNneemllxg1ahTR0THk5BR/jQNjagIr9Begbdu25OefIjc3l5YtW3Jg/37q16+Py+Vi7ty5TsczwD/+8Q8eeOABIiIiOXUqz+k4xjjKCn0F1KpVi/Xr16OqDBkyBFWlW7duiAhPP/200/FqrDfffJN77rkHl8uFx1PgdBxjHGeFPkgmTZqEqjJp0iQAfve73yEi9OzZ0/rxq9CUKVMYPHgwgJ2MzpgAK/RBdnrPfu3atcTExDBz5kzcbjf169fnwIEDTser1j766CMGDBgAgNfrxeWyt7cxYIW+0rRq1Yq8vDzy8vJo3bo1hw4dIiUlBRHhyy+/dDpetfPpp5/Sr18/wIq8MWez/w2VLCYmhjVr1qCqDB06FIBrr70WEeH3vw/K6ftrvLlz59KrVy/AirwxxbH/EVVo4sSJqCpvvvkmAE899RQiQrdu3awf/wJ9/fXXdOvWDfD3yVuRN+Zc9r/CAYMGDUJV2bBhA7G1ajF37lzcbjdJ9eqxd+9ep+OFjaVLl3LVVVcB/iIfEVGh6+gYU21ZoXdQixYtOJmbS35+Pu3atePokSNcdNFFiAiffPKJ0/FC2urVq8nMzAQgPz/firwx52GFPgRERUWxcuVKVJXhw4cD0Lt3b0SExx57zOF0oWfjxo20a9cOgLy8PKKiohxOZExos0IfYsaPH4+q8u677wLw7LPPsmTJEq666ioKCwtLeXb1t3XrVlq2bAnAiRMniImJcTiRMaHPCn2I+ulPf4qqsnnzZlwuF19//TWRkZHUqVOHXbt2OR3PEbt27aJp06YAHDt2LGyvKWpMVbNCH+KaNm1Khw4dyM/Pp0OHDhw7doy0tDREhI8//tjpeFVm3759pKWlAXD48GESExMdTmRM+LBCHyaioqJYunQpqsrIkSMB6Nu3LyLCr371K4fTVa4jR47QsGFDAPbv309SUpLDiYwJL0Ep9CLyiIioiCQHoz1zfn/9619RVd5//30Ann/+eUSETp06Vbt+/OzsbOrVqwfAzp07z7nurzGmdBUu9CKShv96sTsqHseUx80334yqsm3bNhISEli0aBGRkZEkJiayfft2p+NVWE5ODnXr1gVg27ZtNG7c2OFExoSnYOzRjwNGARqEtswFSE9P5/jx43g8Hjp16sSJEyfOXCB96tSpTse7ICdPniQhIQGADRs2kJ6e7nAiY8KXqF54fRaRfsA1qvqgiGwDOqpqsZdaEpFhwDCAlJSUzKysrAte74XIyckJ21EaF5J9165d7N+//8z9Bg0anDmYWZUuJLuqsnTpUgDatGnj2BDKmvaeCSXhnL8ys/fo0WOJqnYs9xNV9bwT8AWwupipH/ANUDuw3DYgubT2VJXMzEytarNmzarydQZLRbJ/+OGHiv/blgLaoUMHzc/PD164UpQ3e35+/pmsy5Ytq5RMZVVT3zOhIJzzV2Z2YLGWocaePZXadaOq16pq27MnYAtwMbAisDffGFgqIqnl/rQxleYnP/kJqsrOnTupXacOy5YtIzo6mvj4eLZs2eJ0vO8pLCwkOjoagAULFtC+fXtnAxlTTVxwH72qrlLVBqqaoaoZwC7gclXdF7R0JmgaN25M9tGjeDwerrrqKnJzc2nWrBkiwuTJk52Oh8/nIzIyEoA5c+ZwxRVXOJzImOrDxtHXMBEREcyfPx9V5dFHHwXgtttuQ0S47777Sn1+QUEBq1at4p133uHXjz7GmGf/dLqL73vGT/g7KY2asGbtWq7q3pOb+t/GvSMe4IknnuDll18mKyuLL774gmXLlrF9+3bcbjcAn332GVdfffX3G9u5E/r3h9q1ITERbrkFdtggL2PKKmin/Avs1ZswMmbMGMaMGcOMGTPo06cP48ePZ/z48bRr146FCxeye/duVq9ezcpVq/hmyXJWr17Nnh3biKuXSmRyOoW10/BunUKht5Df/ebx77V9MucE+XWb4nLF8evN+xgZ24YT60+hK9YR4cnFVZCDnsrBe/IY2Qe+OzXzbYMGU7tOEvWS61G/XjKNk+owZtqHEB3NqjvvJCExkbbvvktMjx6wciXExVX1ZjMm7Ni5XWs4VaV9+/Z89tlnzJ49m2eeeYZVq1YRGxtLRK3a1Em/BG+dNFxJ6UR17UbDeo2RiO/OFll46fX8adyvaZqRwaBBPzvzeKdOnYgYP4m0Deu5cddapnW6lS8uK747JiHw1+c5hS/vBPl5J9iZd5zteSdo/8086h7N5oq23Vkx8TUK83K47couZC1cAH//Ozz8cGVuHmOqBSv0NcyiRYtYunQpi5ctZ8nylWz8dh1er4+4hk3ROmkk3XAfkcnpRNZrjDu29PPJRCTUI6Hfb/nFfQ+Qltb4TLdLhw4dOL5nC80WLUCBGzZ8xRctzt/v7oqMwRUZA4n1zzz24yUfsTghiW93riQzM5NfP/wgffv2hZ494cMPrdAbUwZW6GuYYcNHsHzJYhAh5uLLSeh5H7FNMxF35AW3GVU/g7heD9P3pltY+PV8Mpo154tNx4hIrE/GskUI0HPzQlAFkVLbU1VObV+Bd/UMWuz5lvWtWrP4P1O45JJLvluoTRt4770LzmxMTWKFvoZZumgh27dvZ/78+cyc8z9mz53Krul/IbFxC7wNWhLR8BKiG7Uq0958UbEZ7fFdeQddr7meBj8bC7GJtKmdSsQx/4+2YgoLaH54J5uSm5TYhu9UDrmrv8S75lPqJcbxfw+NpMGDC0npeyMULfIASUlw9Gi5X78xNZEV+hpGRMjIyCAjI4NBgwYBcPz4cb755hvm/m8en8+aw4oZ44ipk0xE6iVoSkuiG7UiIqkxUsreeFy768jO3s+mt35Hyu3PcCuC+LwAuHw+emxeVGyhLziwlYKVM8hdP4/rr7+e/5v8Jl27dvWv78GRwd8IxtQwVugNiYmJXHfddVx33XU8BXi9XlatWsW8efP4Yvb/+Oq/73M0J4f4Jq3xJLcgulErolJb4IqMPqet2l0HUXhsP4emPc8tR3YTETibZozXw43fzmPiFbcCoF4PJ9d/hW/NJ5BzkPuH38vwj/9JaupZv7erW7f4PfcjR/zzjDGlskJvzuF2u2nfvj3t27fn/vvvB2D37t189dVXzJrzP2bNfZct//mWKe4I+p48UWI7+Wd9A7jk4Fa2PXdj8Qs/8QdYtRL+85/vP96mDaxZc+7ya9dC69bleVnG1Fj2gylTJo0aNWLAgAGMf+lvrFu5lKNHDtHg7+PZl9qQPJe72OdEn/VDqmhvCefKj4uDDh3g2WfPnfeTn8CCBVD0dA3btsH8+f55xphSWaE3F6RWrVpccccdpO7aSfRzz5IXEUVhGUbUFFUoLvIiovE98QQsXgwtWpy70D33QEYG9OvnH0750Uf+22lp8ItfBOfFGFPNWaE3FeN2k3P/g/Qd+hLr61/MyWL67YtzMjKab+tn0Hfoi+TcNxJcJbwV4+Jg5kxo2RIGD4ZBg+Dii/2PhelpbI2patZHbyosLiqCzXUuou+QcQxfMIUHvsoixuspcflT7khe7vxTxl85AFwu4qJKeRs2aXJu370xpsxsj95UmNsltGgQj8/lZkP9dDyl/PjK445kff0MVFy0bBCP21W+Lh9jTPlYoTdBMbx7M+Ki3Nyw4SviCvLOu2xcQZ5/uSg3w7s3r6KExtRcVuhNUPRp15BIl9Bz0yJcRS4f7BNX4EDtd281F0rPzQuJdAm929l1aoypbFboTVBER7jJ6ppItLfgzGMnI6M51CSde275Hd/Wz/jegdqYwgKyflSb6Ijih2YaY4LHCr0JmkuWzSNGwBsYNvnnrnfw7lN/Zt7FHfjJkHH8pesg8iKi8YqLGJd/eWNM5bNCb4Jn8mRchR7kskuZN+Vz5vUdDC4XkW5B3W7+13cI86Z8jlzaDpfHAyFwCUNjagIbXmmCJzUVxo7F9dBDXOdycR0we/ZsltzahbioiO9G19y4BF54AWbPdjCsMTWHFXoTPNOmFftwYsxZwy3dbnjkEf9kjKl0UtyFnSt9pSIHge1VvNpk4FAVrzNYLLszLLtzwjl/ZWZPV9X6pS/2fY4UeieIyGJV7eh0jgth2Z1h2Z0TzvlDMbsdjDXGmGrOCr0xxlRzNanQv+p0gAqw7M6w7M4J5/whl73G9NEbY0xNVZP26I0xpkayQm+MMdVctS30IpIkIp+LyMbA37olLOcVkeWB6aOqznlWll4isl5ENonIo8XMjxaRdwPzvxGRDAdiFqsM2e8UkYNFtvVQJ3KeTUReE5EDIrK6hPkiIn8LvK6VInJ5VWc8nzLk7y4ix4ps999XdcbiiEiaiMwSkbUiskZEHixmmZDd9mXMHzrbXlWr5QT8CXg0cPtR4LkSlstxOmsghxvYDDQFooAVQOuzlhkBTAjcHgi863TucmS/E3jJ6azFZL8auBxYXcL8PsAMQIDOwDdOZy5n/u7Ax07nLCZXQ+DywO0EYEMx75mQ3fZlzB8y277a7tED/YDXA7dfB25yLkqZdAI2qeoWVS0AsvC/hqKKvqYpQE+Rcl6Ru3KUJXtIUtW5wJHzLNIPeEP9FgB1RKRh1aQrXRnyhyRV3auqSwO3TwDrgEZnLRay276M+UNGdS70Kaq6N3B7H5BSwnIxIrJYRBaIyE1VE61YjYCdRe7v4tw3zpllVLUQOAbUq5J051eW7AC3Br6CTxGRtKqJVmFlfW2h7EoRWSEiM0SkjdNhzhboguwAfHPWrLDY9ufJDyGy7cP6pGYi8gVQ3CWKflP0jqqqiJQ0jjRdVXeLSFNgpoisUtXNwc5qmAa8o6r5IvIL/N9MrnE4U02wFP97PEdE+gAfAC2cjfQdEYkH/gM8pKrHnc5TXqXkD5ltH9Z79Kp6raq2LWb6ENh/+mte4O+BEtrYHfi7BZiN/5PZCbuBonu5jQOPFbuMiEQAtYHDVZLu/ErNrqqHVTU/cPcfQGYVZauosvy7hCxVPa6qOYHb04FIEUl2OBYAIhKJv0i+parvF7NISG/70vKH0rYP60Jfio+AIYHbQ4APz15AROqKSHTgdjLQBVhbZQm/bxHQQkQuFpEo/Adbzx4FVPQ19QdmauCoj8NKzX5W3+pP8PdphoOPgP8XGAHSGThWpEsw5IlI6unjOCLSCf//ecd3DgKZ/gmsU9W/lLBYyG77suQPpW0f1l03pXgWmCwiP8d/SuSfAohIR+BeVR0KtAL+LiI+/P8Iz6qqI4VeVQtF5H7gU/yjWF5T1TUi8iSwWFU/wv/G+reIbMJ/AG6gE1nPVsbsI0XkJ0Ah/ux3Oha4CBF5B//oiGQR2QWMBiIBVHUCMB3/6I9NwEngLmeSFq8M+fsDw0WkEMgDBobIzkEXYDCwSkSWBx57HGgCYbHty5I/ZLa9nQLBGGOquercdWOMMQYr9MYYU+1ZoTfGmGrOCr0xxlRzVuiNMaaas0JvjDHVnBV6Y4yp5v4/NpxIcD13dLwAAAAASUVORK5CYII=\n", "text/plain": [ "