docs/tutorials/11_using_classical_optimization_solvers_and_models.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Using Classical Optimization Solvers and Models with Qiskit Optimization"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can use classical optimization solvers (CPLEX and Gurobi) with Qiskit Optimization.\n",
    "Docplex and Gurobipy are the Python APIs for CPLEX and Gurobi, respectively.\n",
    "We can load and save an optimization model by Docplex and Gurobipy and can apply CPLEX and Gurobi to `QuadraticProgram`.\n",
    "\n",
    "If you want to use the CPLEX solver, you need to install `pip install 'qiskit-optimization[cplex]'`. Docplex is automatically installed, as a dependent, when you install Qiskit Optimization.\n",
    "\n",
    "If you want to use Gurobi and Gurobipy, you need to install `pip install 'qiskit-optimization[gurobi]'`.\n",
    "\n",
    "**Note**: these solvers, that are installed via pip, are free versions and come with some limitations, such as number of variables. The following links provide further details:\n",
    "\n",
    "- https://pypi.org/project/cplex/\n",
    "- https://pypi.org/project/gurobipy/"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## CplexSolver and GurobiSolver\n",
    "\n",
    "Qiskit Optimization supports the classical solvers of CPLEX and Gurobi as `CplexSolver` and `GurobiSolver`, respectively.\n",
    "We can solve `QuadraticProgram` with `CplexSolver` and `GurobiSolver` as follows."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Problem name: \n",
      "\n",
      "Maximize\n",
      "  x*y\n",
      "\n",
      "Subject to\n",
      "  Linear constraints (1)\n",
      "    x - y <= 0  'c0'\n",
      "\n",
      "  Integer variables (1)\n",
      "    -1 <= y <= 4\n",
      "\n",
      "  Binary variables (1)\n",
      "    x\n",
      "\n"
     ]
    }
   ],
   "source": [
    "from qiskit_optimization.problems import QuadraticProgram\n",
    "\n",
    "# define a problem\n",
    "qp = QuadraticProgram()\n",
    "qp.binary_var(\"x\")\n",
    "qp.integer_var(name=\"y\", lowerbound=-1, upperbound=4)\n",
    "qp.maximize(quadratic={(\"x\", \"y\"): 1})\n",
    "qp.linear_constraint({\"x\": 1, \"y\": -1}, \"<=\", 0)\n",
    "print(qp.prettyprint())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Restricted license - for non-production use only - expires 2025-11-24\n",
      "cplex\n",
      "objective function value: 4.0\n",
      "variable values: x=1.0, y=4.0\n",
      "status: SUCCESS\n",
      "\n",
      "gurobi\n",
      "objective function value: 4.0\n",
      "variable values: x=1.0, y=4.0\n",
      "status: SUCCESS\n"
     ]
    }
   ],
   "source": [
    "from qiskit_optimization.algorithms import CplexOptimizer, GurobiOptimizer\n",
    "\n",
    "cplex_result = CplexOptimizer().solve(qp)\n",
    "gurobi_result = GurobiOptimizer().solve(qp)\n",
    "\n",
    "print(\"cplex\")\n",
    "print(cplex_result.prettyprint())\n",
    "print()\n",
    "print(\"gurobi\")\n",
    "print(gurobi_result.prettyprint())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can set the solver parameter of CPLEX as follows. We can display the solver message of CPLEX by setting `disp=True`.\n",
    "See [Parameters of CPLEX](https://www.ibm.com/docs/en/icos/20.1.0?topic=cplex-parameters) for details of CPLEX parameters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Version identifier: 22.1.1.0 | 2023-02-11 | 22d6266e5\n",
      "CPXPARAM_Read_DataCheck                          1\n",
      "CPXPARAM_Threads                                 1\n",
      "CPXPARAM_TimeLimit                               0.10000000000000001\n",
      "Found incumbent of value 0.000000 after 0.00 sec. (0.00 ticks)\n",
      "Found incumbent of value 4.000000 after 0.00 sec. (0.00 ticks)\n",
      "\n",
      "Root node processing (before b&c):\n",
      "  Real time             =    0.00 sec. (0.00 ticks)\n",
      "Sequential b&c:\n",
      "  Real time             =    0.00 sec. (0.00 ticks)\n",
      "                          ------------\n",
      "Total (root+branch&cut) =    0.00 sec. (0.00 ticks)\n",
      "objective function value: 4.0\n",
      "variable values: x=1.0, y=4.0\n",
      "status: SUCCESS\n"
     ]
    }
   ],
   "source": [
    "result = CplexOptimizer(disp=True, cplex_parameters={\"threads\": 1, \"timelimit\": 0.1}).solve(qp)\n",
    "print(result.prettyprint())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We get the same optimal solution by QAOA as follows."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "objective function value: 4.0\n",
      "variable values: x=1.0, y=4.0\n",
      "status: SUCCESS\n",
      "\n",
      "display the best 5 solution samples\n",
      "SolutionSample(x=array([1., 4.]), fval=4.0, probability=0.024179976472345502, status=<OptimizationResultStatus.SUCCESS: 0>)\n",
      "SolutionSample(x=array([1., 3.]), fval=3.0, probability=0.05819514104929059, status=<OptimizationResultStatus.SUCCESS: 0>)\n",
      "SolutionSample(x=array([1., 2.]), fval=2.0, probability=0.16262615426689528, status=<OptimizationResultStatus.SUCCESS: 0>)\n",
      "SolutionSample(x=array([1., 1.]), fval=1.0, probability=0.1342595428408242, status=<OptimizationResultStatus.SUCCESS: 0>)\n",
      "SolutionSample(x=array([0., 0.]), fval=0.0, probability=0.048555539627953996, status=<OptimizationResultStatus.SUCCESS: 0>)\n"
     ]
    }
   ],
   "source": [
    "from qiskit_optimization.algorithms import MinimumEigenOptimizer\n",
    "\n",
    "from qiskit_algorithms import QAOA\n",
    "from qiskit_algorithms.optimizers import COBYLA\n",
    "from qiskit.primitives import Sampler\n",
    "\n",
    "meo = MinimumEigenOptimizer(QAOA(sampler=Sampler(), optimizer=COBYLA(maxiter=100)))\n",
    "result = meo.solve(qp)\n",
    "print(result.prettyprint())\n",
    "print(\"\\ndisplay the best 5 solution samples\")\n",
    "for sample in result.samples[:5]:\n",
    "    print(sample)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Translators between `QuadraticProgram` and Docplex/Gurobipy\n",
    "\n",
    "Qiskit Optimization can load `QuadraticProgram` from a Docplex model and a Gurobipy model.\n",
    "\n",
    "First, we define an optimization problem by Docplex and Gurobipy."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "// This file has been generated by DOcplex\n",
      "// model name is: docplex\n",
      "// single vars section\n",
      "dvar bool x;\n",
      "dvar int y;\n",
      "\n",
      "maximize\n",
      " [ x*y ];\n",
      " \n",
      "subject to {\n",
      " x <= y;\n",
      "\n",
      "}\n"
     ]
    }
   ],
   "source": [
    "# docplex model\n",
    "from docplex.mp.model import Model\n",
    "\n",
    "docplex_model = Model(\"docplex\")\n",
    "x = docplex_model.binary_var(\"x\")\n",
    "y = docplex_model.integer_var(-1, 4, \"y\")\n",
    "docplex_model.maximize(x * y)\n",
    "docplex_model.add_constraint(x <= y)\n",
    "docplex_model.prettyprint()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\\ Model gurobi\n",
      "\\ LP format - for model browsing. Use MPS format to capture full model detail.\n",
      "Maximize\n",
      " [ 2 x * y ] / 2 \n",
      "Subject To\n",
      " R0: x - y <= 0\n",
      "Bounds\n",
      " -1 <= y <= 4\n",
      "Binaries\n",
      " x\n",
      "Generals\n",
      " y\n",
      "End\n",
      "\n"
     ]
    }
   ],
   "source": [
    "# gurobi model\n",
    "from io import StringIO\n",
    "from tempfile import NamedTemporaryFile\n",
    "from os import remove\n",
    "\n",
    "import gurobipy as gp\n",
    "\n",
    "\n",
    "def gpy_display(mdl):\n",
    "    \"\"\"Convenience function to pretty-print a Gurobipy-Model.\"\"\"\n",
    "    with NamedTemporaryFile(suffix=\".lp\") as tmp_file:\n",
    "        mdl.write(tmp_file.name)\n",
    "\n",
    "        with open(tmp_file.name, \"r\") as f:\n",
    "            print(f.read())\n",
    "\n",
    "\n",
    "gurobipy_model = gp.Model(\"gurobi\")\n",
    "x = gurobipy_model.addVar(vtype=gp.GRB.BINARY, name=\"x\")\n",
    "y = gurobipy_model.addVar(vtype=gp.GRB.INTEGER, lb=-1, ub=4, name=\"y\")\n",
    "gurobipy_model.setObjective(x * y, gp.GRB.MAXIMIZE)\n",
    "gurobipy_model.addConstr(x - y <= 0)\n",
    "gurobipy_model.update()\n",
    "gpy_display(gurobipy_model)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can generate `QuadraticProgram` object from both Docplex and Gurobipy models. We see that the two `QuadraticProgram` objects generated from Docplex and Gurobipy are identical."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "QuadraticProgram obtained from docpblex\n",
      "Problem name: docplex\n",
      "\n",
      "Maximize\n",
      "  x*y\n",
      "\n",
      "Subject to\n",
      "  Linear constraints (1)\n",
      "    x - y <= 0  'c0'\n",
      "\n",
      "  Integer variables (1)\n",
      "    -1 <= y <= 4\n",
      "\n",
      "  Binary variables (1)\n",
      "    x\n",
      "\n",
      "-------------\n",
      "QuadraticProgram obtained from gurobipy\n",
      "Problem name: gurobi\n",
      "\n",
      "Maximize\n",
      "  x*y\n",
      "\n",
      "Subject to\n",
      "  Linear constraints (1)\n",
      "    x - y <= 0  'R0'\n",
      "\n",
      "  Integer variables (1)\n",
      "    -1 <= y <= 4\n",
      "\n",
      "  Binary variables (1)\n",
      "    x\n",
      "\n"
     ]
    }
   ],
   "source": [
    "from qiskit_optimization.translators import from_docplex_mp, from_gurobipy\n",
    "\n",
    "qp = from_docplex_mp(docplex_model)\n",
    "print(\"QuadraticProgram obtained from docpblex\")\n",
    "print(qp.prettyprint())\n",
    "print(\"-------------\")\n",
    "print(\"QuadraticProgram obtained from gurobipy\")\n",
    "qp2 = from_gurobipy(gurobipy_model)\n",
    "print(qp2.prettyprint())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can generate a Docplex model and a Gurobipy model from `QuadraticProgram` too."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "convert docplex to gurobipy via QuadraticProgram\n",
      "\\ Model docplex\n",
      "\\ LP format - for model browsing. Use MPS format to capture full model detail.\n",
      "Maximize\n",
      " [ 2 x * y ] / 2 \n",
      "Subject To\n",
      " c0: x - y <= 0\n",
      "Bounds\n",
      " -1 <= y <= 4\n",
      "Binaries\n",
      " x\n",
      "Generals\n",
      " y\n",
      "End\n",
      "\n",
      "\n",
      "convert gurobipy to docplex via QuadraticProgram\n",
      "\\ This file has been generated by DOcplex\n",
      "\\ ENCODING=ISO-8859-1\n",
      "\\Problem name: gurobi\n",
      "\n",
      "Maximize\n",
      " obj: [ 2 x*y ]/2\n",
      "Subject To\n",
      " R0: x - y <= 0\n",
      "\n",
      "Bounds\n",
      " 0 <= x <= 1\n",
      " -1 <= y <= 4\n",
      "\n",
      "Binaries\n",
      " x\n",
      "\n",
      "Generals\n",
      " y\n",
      "End\n",
      "\n"
     ]
    }
   ],
   "source": [
    "from qiskit_optimization.translators import to_gurobipy, to_docplex_mp\n",
    "\n",
    "gmod = to_gurobipy(from_docplex_mp(docplex_model))\n",
    "print(\"convert docplex to gurobipy via QuadraticProgram\")\n",
    "gpy_display(gmod)\n",
    "\n",
    "dmod = to_docplex_mp(from_gurobipy(gurobipy_model))\n",
    "print(\"\\nconvert gurobipy to docplex via QuadraticProgram\")\n",
    "print(dmod.export_as_lp_string())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Indicator constraints of Docplex\n",
    "\n",
    "`from_docplex_mp` supports indicator constraints, e.g., `u = 0 => x + y <= z` (u: binary variable) when we convert a Docplex model into `QuadraticProgram`. It converts indicator constraints into linear constraints using the big-M formulation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\\ This file has been generated by DOcplex\n",
      "\\ ENCODING=ISO-8859-1\n",
      "\\Problem name: docplex\n",
      "\n",
      "Maximize\n",
      " obj: 3 x + y - z\n",
      "Subject To\n",
      " lc1: x = 1 -> y - z >= 0\n",
      "\n",
      "Bounds\n",
      " 0 <= x <= 1\n",
      " -1 <= y <= 2\n",
      " -1 <= z <= 2\n",
      "\n",
      "Binaries\n",
      " x\n",
      "\n",
      "Generals\n",
      " y z\n",
      "End\n",
      "\n"
     ]
    }
   ],
   "source": [
    "ind_mod = Model(\"docplex\")\n",
    "x = ind_mod.binary_var(\"x\")\n",
    "y = ind_mod.integer_var(-1, 2, \"y\")\n",
    "z = ind_mod.integer_var(-1, 2, \"z\")\n",
    "ind_mod.maximize(3 * x + y - z)\n",
    "ind_mod.add_indicator(x, y >= z, 1)\n",
    "print(ind_mod.export_as_lp_string())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's compare the solutions of the model with an indicator constraint by\n",
    "\n",
    "1. applying CPLEX directly to the Docplex model (without translating it to `QuadraticProgram`. CPLEX solver natively supports the indicator constraints),\n",
    "2. applying QAOA to `QuadraticProgram` obtained by `from_docplex_mp`.\n",
    "\n",
    "We see the solutions are same."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "QAOA\n",
      "objective function value: 6.0\n",
      "variable values: x=1.0, y=2.0, z=-1.0\n",
      "status: SUCCESS\n",
      "-----\n",
      "CPLEX\n",
      "solution for: docplex\n",
      "objective: 6\n",
      "status: OPTIMAL_SOLUTION(2)\n",
      "x=1\n",
      "y=2\n",
      "z=-1\n",
      "\n"
     ]
    }
   ],
   "source": [
    "qp = from_docplex_mp(ind_mod)\n",
    "result = meo.solve(qp)  # apply QAOA to QuadraticProgram\n",
    "print(\"QAOA\")\n",
    "print(result.prettyprint())\n",
    "print(\"-----\\nCPLEX\")\n",
    "print(ind_mod.solve())  # apply CPLEX directly to the Docplex model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<h3>Version Information</h3><table><tr><th>Software</th><th>Version</th></tr><tr><td><code>qiskit</code></td><td>0.45.1</td></tr><tr><td><code>qiskit_optimization</code></td><td>0.7.0</td></tr><tr><td><code>qiskit_algorithms</code></td><td>0.2.1</td></tr><tr><th colspan='2'>System information</th></tr><tr><td>Python version</td><td>3.10.13</td></tr><tr><td>Python compiler</td><td>GCC 13.2.1 20230728 (Red Hat 13.2.1-1)</td></tr><tr><td>Python build</td><td>main, Aug 28 2023 00:00:00</td></tr><tr><td>OS</td><td>Linux</td></tr><tr><td>CPUs</td><td>12</td></tr><tr><td>Memory (Gb)</td><td>31.056411743164062</td></tr><tr><td colspan='2'>Wed Dec 06 10:30:30 2023 CET</td></tr></table>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<div style='width: 100%; background-color:#d5d9e0;padding-left: 10px; padding-bottom: 10px; padding-right: 10px; padding-top: 5px'><h3>This code is a part of Qiskit</h3><p>&copy; Copyright IBM 2017, 2023.</p><p>This code is licensed under the Apache License, Version 2.0. You may<br>obtain a copy of this license in the LICENSE.txt file in the root directory<br> of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.<p>Any modifications or derivative works of this code must retain this<br>copyright notice, and modified files need to carry a notice indicating<br>that they have been altered from the originals.</p></div>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import tutorial_magics\n",
    "\n",
    "%qiskit_version_table\n",
    "%qiskit_copyright"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}