MILP-demo.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. Import necessary module"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.optimize import milp, LinearConstraint, Bounds\n",
    "from scipy.sparse import csc_matrix\n",
    "from qiskit_optimization import QuadraticProgram\n",
    "from qiskit import QiskitError\n",
    "import numpy as np\n",
    "from qiskit_optimization.algorithms import GurobiOptimizer\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "%load_ext autoreload\n",
    "%autoreload 2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. Create a demo problem\n",
    "\n",
    "### Example 1\n",
    "\n",
    "\n",
    "$${\\rm max} \\ x_1+2x_2+2x_3+5x_4 \\\\\n",
    "\n",
    "-x_1+x_2+x_3+x_4 \\leq 1 \\\\ 3x_1+2x_2 - x_3+x_4 \\leq 12 \\\\ x_1+x_2 +x_3+x_4 \\leq 10 \\\\ 2x_1+3x_2 + 3x_3+x_4 \\geq 12 \\\\ \n",
    "\n",
    "x_1, x_2 \\in \\mathcal{Z}\\\\\n",
    "x_4 \\in \\{0, 1\\}\n",
    "$$\n",
    "\n",
    "Let's solve it with the milp solver!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "            fun: -14.0\n",
      "        message: 'Optimization terminated successfully. (HiGHS Status 7: Optimal)'\n",
      " mip_dual_bound: -14.0\n",
      "        mip_gap: 0.0\n",
      " mip_node_count: 1\n",
      "         status: 0\n",
      "        success: True\n",
      "              x: array([4., 1., 3., 1.])\n"
     ]
    }
   ],
   "source": [
    "## constant value in objective\n",
    "\n",
    "c = -np.array([1, 2, 1, 5])  \n",
    "A = np.array([[-1, 1, 1, 1], [3, 2, -1, 1], [1, 1, 1, 1], [2, 3, 3, 1]])\n",
    "b_u = np.array([1, 12, 10, np.inf])\n",
    "b_l = np.array([-np.inf, -np.inf, -np.inf, 12])\n",
    "\n",
    "constraints = LinearConstraint(A, b_l, b_u)\n",
    "integrality = np.array([1, 1, 0, 1])\n",
    "bounds = Bounds([0, 0, 0, 0],[np.inf, np.inf, np.inf, 1])\n",
    "#integrality = np.zeros(4)\n",
    "res = milp(c=c, constraints=constraints, integrality=integrality, bounds=bounds)\n",
    "\n",
    "print(res)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We see the MILP solver tells us the object function is maximized when $x =  [4, 1, 3, 1]$.\n",
    "\n",
    "As the next step, let's create a corresponding `QuadraticProgram` class in `qiskit_optimization`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Problem name: example-1\n",
      "\n",
      "Maximize\n",
      "  x1 + 2*x2 + x3 + 5*x4\n",
      "\n",
      "Subject to\n",
      "  Linear constraints (4)\n",
      "    -x1 + x2 + x3 + x4 <= 1  'c0'\n",
      "    x1 + x2 + x3 + x4 <= 10  'c1'\n",
      "    3*x1 + 2*x2 - x3 + x4 <= 12  'c2'\n",
      "    2*x1 + 3*x2 + 3*x3 + x4 >= 12  'c3'\n",
      "\n",
      "  Integer variables (2)\n",
      "    0 <= x1\n",
      "    0 <= x2\n",
      "\n",
      "  Continuous variables (1)\n",
      "    0 <= x3\n",
      "\n",
      "  Binary variables (1)\n",
      "    x4\n",
      "\n"
     ]
    }
   ],
   "source": [
    "## Consider the possibility of using sparse matrix(to_dict)\n",
    "## 1e20=inf\n",
    "qp = QuadraticProgram('example-1')\n",
    "qp.integer_var(name='x1')\n",
    "qp.integer_var(name='x2')\n",
    "qp.continuous_var(name='x3')\n",
    "qp.binary_var(name='x4')\n",
    "#qp.binary_var_list(2000)\n",
    "qp.maximize(linear={'x1': 1, 'x2': 2, 'x3': 1, 'x4': 5})\n",
    "#qp.minimize(linear={'x2': 1})\n",
    "qp.linear_constraint({'x1': -1, 'x2': 1, 'x3': 1, 'x4': 1},'<=',1)\n",
    "qp.linear_constraint({'x1': 1, 'x2': 1, 'x3': 1, 'x4': 1},'<=',10)\n",
    "qp.linear_constraint({'x1': 3, 'x2': 2, 'x3': -1, 'x4': 1},'<=',12)\n",
    "qp.linear_constraint({'x1': 2, 'x2': 3, 'x3': 3, 'x4': 1},'>=',12)\n",
    "print(qp.prettyprint())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We solve the `qp(QuadraticProgram)` with the Gurobi solver."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Restricted license - for non-production use only - expires 2023-10-25\n",
      "objective function value: 14.0\n",
      "variable values: x1=5.0, x2=0.0, x3=4.0, x4=1.0\n",
      "status: SUCCESS\n"
     ]
    }
   ],
   "source": [
    "gurobi_result = GurobiOptimizer().solve(qp)\n",
    "print(gurobi_result.prettyprint())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 61,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "['__class__',\n",
       " '__delattr__',\n",
       " '__dict__',\n",
       " '__dir__',\n",
       " '__doc__',\n",
       " '__eq__',\n",
       " '__format__',\n",
       " '__ge__',\n",
       " '__getattribute__',\n",
       " '__getitem__',\n",
       " '__gt__',\n",
       " '__hash__',\n",
       " '__init__',\n",
       " '__init_subclass__',\n",
       " '__le__',\n",
       " '__lt__',\n",
       " '__module__',\n",
       " '__ne__',\n",
       " '__new__',\n",
       " '__reduce__',\n",
       " '__reduce_ex__',\n",
       " '__repr__',\n",
       " '__setattr__',\n",
       " '__sizeof__',\n",
       " '__str__',\n",
       " '__subclasshook__',\n",
       " '__weakref__',\n",
       " '_fval',\n",
       " '_raw_results',\n",
       " '_samples',\n",
       " '_status',\n",
       " '_variable_names',\n",
       " '_variables',\n",
       " '_variables_dict',\n",
       " '_x',\n",
       " 'fval',\n",
       " 'get_correlations',\n",
       " 'prettyprint',\n",
       " 'raw_results',\n",
       " 'samples',\n",
       " 'status',\n",
       " 'variable_names',\n",
       " 'variables',\n",
       " 'variables_dict',\n",
       " 'x']"
      ]
     },
     "execution_count": 61,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "dir(gurobi_result)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3. Build a Converter\n",
    "\n",
    "The converter should be able to extract necessary information from the `QuadraticProblem` class and input it to the milp solver."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[]"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "qp.quadratic_constraints"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 90,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "            fun: -14.0\n",
      "        message: 'Optimization terminated successfully. (HiGHS Status 7: Optimal)'\n",
      " mip_dual_bound: -14.0\n",
      "        mip_gap: 0.0\n",
      " mip_node_count: 1\n",
      "         status: 0\n",
      "        success: True\n",
      "              x: array([4., 1., 3., 1.])\n"
     ]
    }
   ],
   "source": [
    "## First, Check if qp is linear\n",
    "## Check objective is linear\n",
    "##from qiskit_optimization import vartype\n",
    "\n",
    "if qp.objective.quadratic.to_dict()!={}:\n",
    "    raise QiskitError('Obejective function is not linear!')\n",
    "\n",
    "## Check constrains are linear\n",
    "if qp.quadratic_constraints:\n",
    "    raise QiskitError('Constraints are not linear!')\n",
    "\n",
    "## Sense: 1 for minimization and -1 for maximization \n",
    "sense = qp.objective.sense.value\n",
    "\n",
    "## cost function for milp solver\n",
    "c_qp = qp.objective.linear.to_array() * sense  \n",
    "\n",
    "## constraints for milp solver\n",
    "for i, constraint in enumerate(qp.linear_constraints):\n",
    "\n",
    "    constraint_array = constraint.linear.to_array()\n",
    "    constraint_sense = constraint.sense.value # 0 for leq, 1 for geq\n",
    "    constraint_value = constraint.rhs\n",
    "    if i==0:\n",
    "        A_qp = [constraint_array]\n",
    "        bl_qp = np.array([constraint_value if constraint_sense==1 else -np.inf])\n",
    "        bu_qp = np.array([constraint_value if constraint_sense==0 else np.inf])\n",
    "    else:\n",
    "        A_qp.append(constraint_array)  ## Not efficient. Use list directly.\n",
    "        bl_qp = np.append(bl_qp, [constraint_value if constraint_sense==1 else -np.inf])\n",
    "        bu_qp = np.append(bu_qp, [constraint_value if constraint_sense==0 else np.inf])\n",
    "\n",
    "constraints = LinearConstraint(A_qp, bl_qp, bu_qp)\n",
    "\n",
    "## integrity for milp solver\n",
    "## Use Enum and VarType to compare\n",
    "integrality = np.array([ 1 if variable.vartype.value==2 or variable.vartype.value==1 else 0 for variable in qp.variables]) ## check on other vartype!\n",
    "m\n",
    "raw_res = milp(c=c_qp, constraints=constraints, integrality=integrality, bounds=bounds)\n",
    "\n",
    "print(raw_res)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "            fun: -14.0\n",
      "        message: 'Optimization terminated successfully. (HiGHS Status 7: Optimal)'\n",
      " mip_dual_bound: -14.0\n",
      "        mip_gap: 0.0\n",
      " mip_node_count: 1\n",
      "         status: 0\n",
      "        success: True\n",
      "              x: array([4., 1., 3., 1.])\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/dranz/Library/Python/3.9/lib/python/site-packages/scipy/sparse/_index.py:103: SparseEfficiencyWarning: Changing the sparsity structure of a csc_matrix is expensive. lil_matrix is more efficient.\n",
      "  self._set_intXint(row, col, x.flat[0])\n"
     ]
    }
   ],
   "source": [
    "## Sparse converter based on to_dict\n",
    "\n",
    "if qp.objective.quadratic.to_dict()!={}:\n",
    "    raise QiskitError('Obejective function is not linear!')\n",
    "\n",
    "## Check constrains are linear\n",
    "if qp.quadratic_constraints:\n",
    "    raise QiskitError('Constraints are not linear!')\n",
    "\n",
    "## Sense: 1 for minimization and -1 for maximization \n",
    "sense = qp.objective.sense.value\n",
    "\n",
    "## cost function for milp solver\n",
    "c_qp = qp.objective.linear.to_array() * sense  \n",
    "\n",
    "\n",
    "A_qp = csc_matrix((len(qp.linear_constraints), len(qp.variables)))\n",
    "bl_qp = np.array([])\n",
    "bu_qp = np.array([])\n",
    "## constraints for milp solver\n",
    "for i, constraint in enumerate(qp.linear_constraints):\n",
    "    constraint_dict = constraint.linear.to_dict()\n",
    "    constraint_sense = constraint.sense.value # 0 for leq, 1 for geq\n",
    "    constraint_value = constraint.rhs\n",
    "    \n",
    "    for var in constraint_dict:\n",
    "        A_qp[i,var]= constraint_dict[var]\n",
    "\n",
    "    bl_qp = np.append(bl_qp, [constraint_value if constraint_sense==1 else -np.inf])\n",
    "    #print(bl_qp)\n",
    "    bu_qp = np.append(bu_qp, [constraint_value if constraint_sense==0 else np.inf])\n",
    "\n",
    "constraints = LinearConstraint(A_qp, bl_qp, bu_qp)\n",
    "\n",
    "\n",
    "## integrity for milp solver\n",
    "## Use Enum and VarType to compare\n",
    "integrality = np.array([ 1 if variable.vartype.value==2 or variable.vartype.value==1 else 0 for variable in qp.variables]) ## check on other vartype!\n",
    "lower_bounds = np.array([var.lowerbound for var in qp.variables])\n",
    "upper_bounds = np.array([var.upperbound for var in qp.variables])\n",
    "bounds = Bounds(lower_bounds,upper_bounds)\n",
    "raw_res = milp(c=c_qp, constraints=constraints, integrality=integrality, bounds=bounds)\n",
    "\n",
    "print(raw_res)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "['fun',\n",
       " 'message',\n",
       " 'mip_dual_bound',\n",
       " 'mip_gap',\n",
       " 'mip_node_count',\n",
       " 'status',\n",
       " 'success',\n",
       " 'x']"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "raw_"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 65,
   "metadata": {},
   "outputs": [],
   "source": [
    "A = csc_matrix((4,4))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 67,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.0"
      ]
     },
     "execution_count": 67,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A[0,0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[-1.,  1.,  1.],\n",
       "       [ 3.,  2., -1.]])"
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A = qp.get_linear_constraint(0).linear.to_array()\n",
    "B = qp.get_linear_constraint(1).linear.to_array()\n",
    "A = np.vstack((A, B))\n",
    "A"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 62,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<ConstraintSense.LE: 0>"
      ]
     },
     "execution_count": 62,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "qp.get_linear_constraint(0).sense"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [],
   "source": [
    "integrality = np.array([ 1 if variable.vartype.value==2 else 0 for variable in qp.variables])"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3.9.6 64-bit",
   "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.9.6"
  },
  "orig_nbformat": 4,
  "vscode": {
   "interpreter": {
    "hash": "31f2aee4e71d21fbe5cf8b01ff0e069b9275f58929596ceb00d14d90e3e16cd6"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}