{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# ILP basics\n",
    "\n",
    "_Combinatorial Optimization course, FEE CTU in Prague. Created by [Industrial Informatics Department](http://industrialinformatics.fel.cvut.cz)._\n",
    "\n",
    "During the second lab of the course, we concentrate on the (I)LP modeling basics. The goal is to go through several simple examples and try to implement them using the ILP formalism. \n",
    "In this document, we show the basic syntax on a simple LP problem. \n",
    "Afterward, several NP-hard problems are formally presented, for which the students should fill the ILP formulations.\n",
    "\n",
    "\n",
    "## Introduction: Linear programming\n",
    "\n",
    "In the linear programming, we try to solve the following problem:\n",
    "\n",
    "$$\n",
    "\\min_{x} \\{ c^{\\text{T}} x \\ | \\ Ax \\leq b \\}, \\text{ where } c \\in \\mathbb{R}^{n}, x \\in \\mathbb{R}^{n}, A \\in \\mathbb{R}^{m \\times n}, b \\in \\mathbb{R}^{m}\n",
    "$$\n",
    "\n",
    "i.e., we want to minimize a linear objective function, subject to linear inequality constraints defininy the feasible region (here, convex polytope).\n",
    "\n",
    "### Example\n",
    "\n",
    "Take as an example the following problem:\n",
    "\n",
    "\\begin{align}\n",
    "         \\min_{x,y} {-x + 2y} \\ \\ \\text{s.t.} \\\\\n",
    "        -4x -9y \\leq -18 \\\\\n",
    "        \\frac{3}{2}x - y \\leq \\frac{27}{4} \\\\\n",
    "        \\frac{8}{17}x - y \\geq -2 \n",
    "\\end{align}\n",
    "\n",
    "Rewritten to the matrix form, the same problem would look as follows:\n",
    "\n",
    "\\begin{align}\n",
    "            \\min_{x,y} \\begin{bmatrix}\n",
    "            -1 & 2\n",
    "            \\end{bmatrix} \\cdot \n",
    "            \\begin{bmatrix}\n",
    "            x \\\\ y\n",
    "            \\end{bmatrix} \\ \\text{s.t.} \\\\\n",
    "            \\begin{bmatrix}\n",
    "            -4 & -9 \\\\ \\frac{3}{2} & -1 \\\\ -\\frac{8}{17} & 1 \n",
    "            \\end{bmatrix} \\cdot\n",
    "            \\begin{bmatrix}\n",
    "            x \\\\ y\n",
    "            \\end{bmatrix} \\leq \n",
    "            \\begin{bmatrix}\n",
    "            -18 \\\\ \\frac{27}{4} \\\\ 2\n",
    "            \\end{bmatrix}        \n",
    "\\end{align}\n",
    "\n",
    "For this simple problem, we can visualize the individual constraints, the feasible region and even the optima, which is $x = \\frac{9}{2}, y = 0$.\n",
    "\n",
    "![LP illustration](images/lp.png)\n",
    "\n",
    "Now, let us try to rewrite this problem using the Gurobi syntax. Then, we can call the solver and see if the solution that we have found 'graphically' is, indeed, correct.\n",
    "\n",
    "The model itself will consist of variables (including definition of their domains) and constraints. Here, the implementation is pretty straightforward:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-02-19T14:09:43.100195Z",
     "start_time": "2021-02-19T14:09:43.043823Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Using license file /home/benedond/Apps/gurobi/gurobi910/gurobi.lic\n",
      "Academic license - for non-commercial use only - expires 2021-04-16\n",
      "Gurobi Optimizer version 9.1.0 build v9.1.0rc0 (linux64)\n",
      "Thread count: 2 physical cores, 4 logical processors, using up to 4 threads\n",
      "Optimize a model with 3 rows, 2 columns and 6 nonzeros\n",
      "Model fingerprint: 0xb2d66890\n",
      "Coefficient statistics:\n",
      "  Matrix range     [5e-01, 9e+00]\n",
      "  Objective range  [1e+00, 2e+00]\n",
      "  Bounds range     [0e+00, 0e+00]\n",
      "  RHS range        [2e+00, 2e+01]\n",
      "Presolve removed 3 rows and 2 columns\n",
      "Presolve time: 0.01s\n",
      "Presolve: All rows and columns removed\n",
      "Iteration    Objective       Primal Inf.    Dual Inf.      Time\n",
      "       0   -4.5000000e+00   0.000000e+00   0.000000e+00      0s\n",
      "\n",
      "Solved in 0 iterations and 0.01 seconds\n",
      "Optimal objective -4.500000000e+00\n",
      "x: 4.5\n",
      "y: 0.0\n"
     ]
    }
   ],
   "source": [
    "# Gurobi model for the example problem\n",
    "import gurobipy as g  # import Gurobi\n",
    "\n",
    "# Input data ----------------------\n",
    "# The program above\n",
    "\n",
    "# Model ---------------------------\n",
    "model = g.Model()  # create an empty model\n",
    "# - variables\n",
    "x = model.addVar(vtype=g.GRB.CONTINUOUS, name=\"x\")  # method addVar is used to add variable\n",
    "y = model.addVar(vtype=g.GRB.CONTINUOUS, name=\"y\")\n",
    "# - constraints\n",
    "model.addConstr(-4*x - 9*y <= -18)  # method addConstr is used to add constraint\n",
    "model.addConstr(3/2*x - y <= 27/4)\n",
    "model.addConstr(8/17*x - y >= -2)\n",
    "# - objective\n",
    "model.setObjective(-x + 2*y, g.GRB.MINIMIZE)\n",
    "\n",
    "model.optimize()\n",
    "\n",
    "# Solution ------------------------\n",
    "print(\"x:\", x.X)  # use [variable].X to get the optimal value\n",
    "print(\"y:\", y.X)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**While executing the model, Gurobi solver prints some information:**\n",
    "\n",
    "- `Thread count: 2 physical cores, 4 logical processors, using up to 4 threads`\n",
    "\n",
    "Some parts of the solver can work in parallel - which might be beneficial when solving complex models.\n",
    "\n",
    "- `Optimize a model with 3 rows, 2 columns and 6 nonzeros`\n",
    "\n",
    "Here, we see that the matrix $A$ contains 3 rows (constraints), 2 columns (variables) and 6 non-zero coefficients at total.\n",
    "\n",
    "```\n",
    "Presolve removed 3 rows and 2 columns\n",
    "Presolve time: 0.01s\n",
    "Presolve: All rows and columns removed\n",
    "Iteration    Objective       Primal Inf.    Dual Inf.      Time\n",
    "       0   -4.5000000e+00   0.000000e+00   0.000000e+00      0s\n",
    "\n",
    "Solved in 0 iterations and 0.01 seconds\n",
    "```\n",
    "\n",
    "This parts just says that Gurobi was able to find the solution during the \"Presolve\" phase.\n",
    "\n",
    "- `Optimal objective -4.500000000e+00`\n",
    "\n",
    "The value of the objective function $c^{\\text{T}} x $ is -4.5."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Integer linear programming\n",
    "\n",
    "Now we solve similar problem, but the domains of (some) variables are restricted to be integers.\n",
    "\n",
    "$$\n",
    "\\min_{x} \\{c^{T} x \\ \\vert \\ Ax \\leq b\\}, c \\in \\mathbb{R}^{n}, x \\in \\mathbb{Z}^{n}, A \\in \\mathbb{R}^{m \\times n}, b \\in \\mathbb{R}^{m}\n",
    "$$\n",
    "\n",
    "- Note that solving the problem as LP and rounding the solution **may not lead to optimal solution** (it might even lead to infeasible solution - see the previous example).   \n",
    "- Also note that optimizing over a subset will not give better solution (it may only worsen). Therefore, the optimal objective of LP provides a bound for the optimal objective of ILP (lower bound for the minimization).\n",
    "\n",
    "Examine the following situation and compare the optimal objectives of LP and ILP variants of the same problem:\n",
    "\n",
    "\\begin{align}\n",
    "\\max_{x,y} y \\\\ \\text{ subject to} \\\\\n",
    "-x+y \\leq 1 \\\\\n",
    "2x+3y \\leq 12 \\\\\n",
    "3x+2y \\leq 12 \\\\\n",
    "x \\geq 0 \\\\\n",
    "y \\geq 0\n",
    "\\end{align}\n",
    "\n",
    "![ILP vs LP](images/ILP_LP.png)\n",
    "\n",
    "The area inside of the blue polygon represents the feasible region of the LP formulation, while the red dots represent the feasible solutions of the ILP formulation.\n",
    "\n",
    "**Now, lets enjoy some fun while implementing the models for some well-known NP-hard problems. The goal is to understand the problems and formalize them in Gurobi language.**\n",
    "\n",
    "### Problem 1: Two partition\n",
    "\n",
    "\n",
    "Given multiset $S$ of positive integers, we ask if it can be partitioned to two subsets $S_1$, $S_2$, such that sum of the numbers in $S_1$ equals to the sum of the numbers in $S_2$, and each number is assigned (either to $S_1$ or to $S_2$).\n",
    "\n",
    "**Example 1** Multiset $S = \\{ 1,1,1,2,2,3\\}$ can be partitioned, e.g., to $S_1 = \\{ 2, 3 \\}$,  $S_2 = \\{1,1,1,2\\}$. Note that it is not a unique solution.\n",
    "\n",
    "**Example 2** Multiset $S = \\{2, 4\\}$ cannot be partitioned.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-02-19T14:09:46.562163Z",
     "start_time": "2021-02-19T14:09:46.558236Z"
    }
   },
   "outputs": [],
   "source": [
    "import gurobipy as g\n",
    "\n",
    "# Example input\n",
    "S = [100, 50, 50, 50, 20, 20, 10]\n",
    "\n",
    "# TODO: implement the model and find the solution"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Problem 2: Maximum independent set\n",
    "\n",
    "Given graph $G = (V,E)$, we are trying to find subset $V^{\\star} \\subseteq V$, such that $V^{\\star}$ is as large as possible and for each pair of vertices $u,v \\in V^{\\star}$, $u \\neq v$, it holds that there is no edge $e = \\{u,v\\}$ in the original graph."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-02-19T14:09:47.930754Z",
     "start_time": "2021-02-19T14:09:47.926977Z"
    }
   },
   "outputs": [],
   "source": [
    "import gurobipy as g\n",
    "\n",
    "# Example input\n",
    "n = 5  # number of vertices\n",
    "edges = [(0, 1), (2, 3), (0, 4), (3, 1), (3, 4)]  # list of edges\n",
    "\n",
    "# TODO: implement the model and find the solution"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Problem 3: Minimum vertex cover\n",
    "\n",
    "Given graph $G = (V,E)$, we are looking for subset $V^{\\star} \\subseteq V$, such that $V^{\\star}$ is as small as possible and each edge $e = \\{u,v\\}$ is `covered' by at least one vertex $x \\in V^{\\star}$, i.e. $\\forall e = \\{u,v\\} \\in E \\ \\exists x \\in V^{\\star}$ such that $x = u$ or $x = v$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-02-19T14:09:49.231082Z",
     "start_time": "2021-02-19T14:09:49.223753Z"
    }
   },
   "outputs": [],
   "source": [
    "import gurobipy as g\n",
    "\n",
    "# Example input\n",
    "n = 5  # number of vertices\n",
    "edges = [(0, 1), (2, 3), (0, 4), (3, 1), (3, 4)]  # list of edges\n",
    "\n",
    "# TODO: implement the model and find the solution"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Problem 4: SAT (boolean satisfiability problem)\n",
    "\n",
    "Given formula $\\varphi$ (in CNF), we are looking for truth assignment (True / False) to each variable, such that $\\varphi$ is satisfied (evaluates to True).\n",
    "\n",
    "**Example 1** Formula $\\varphi = (x_1 \\vee x_2) \\wedge (\\neg x_2 \\vee x_3)$ can be satisfied by $x_1 := \\text{True}$, $x_3 := \\text{True}$, $x_2 := \\text{False}$.\n",
    "\n",
    "**Example 2** Formula $\\varphi = (x_1) \\wedge (\\neg x_1)$ cannot be satisfied."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import gurobipy as g\n",
    "\n",
    "# Example input\n",
    "# SAT formula: (x0 or ~x1) and (~x0 or x1 or x2) <- encode directly into the model\n",
    "\n",
    "# TODO: implement the model and find the solution"
   ]
  }
 ],
 "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.8.10"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  },
  "varInspector": {
   "cols": {
    "lenName": 16,
    "lenType": 16,
    "lenVar": 40
   },
   "kernels_config": {
    "python": {
     "delete_cmd_postfix": "",
     "delete_cmd_prefix": "del ",
     "library": "var_list.py",
     "varRefreshCmd": "print(var_dic_list())"
    },
    "r": {
     "delete_cmd_postfix": ") ",
     "delete_cmd_prefix": "rm(",
     "library": "var_list.r",
     "varRefreshCmd": "cat(var_dic_list()) "
    }
   },
   "types_to_exclude": [
    "module",
    "function",
    "builtin_function_or_method",
    "instance",
    "_Feature"
   ],
   "window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}