{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "(user_guide_optimization_lp)="
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Linear Programming (LP)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "\ud83d\udcd6 {ref}`Optimization API Reference <api_optimization>`"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "The library offers a solver for a wide range of linear optimization problems, including continuous, integer and mixed-integer problems. The {py:class}`~sigmaepsilon.math.optimize.lp.LinearProgrammingProblem` class uses the [linprog module from SciPy](https://docs.scipy.org/doc/scipy/tutorial/optimize.html#linear-programming-linprog) (which eventually calls the [HIGHS](https://highs.dev/) solver) for the solution, combined with the power of SymPy for being able to pose problems directly in symbolic form. The class handles all the preprocessing required to bring the problem into the form required by SciPy."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First, import some stuff. Note that the notebook requires matplotlib to be installed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "from sigmaepsilon.math.function import Function, Relation\n",
    "from sigmaepsilon.math.optimize import LinearProgrammingProblem as LPP\n",
    "import sympy as sy"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "editable": true,
    "raw_mimetype": "",
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "One of the great features of the optimization module is that it handles symbolic functions pretty well. Problems can be defined using `SymPy` expressions, or simple strings. To solve a problem, call the `solve` method, which returns an instance of [scipy.optimize.OptimizeResult](https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.OptimizeResult.html#scipy.optimize.OptimizeResult)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "x1, x2, x3, x4 = variables = sy.symbols(['x1', 'x2', 'x3', 'x4'])\n",
    "obj = Function(3*x1 + 9*x3 + x2 + x4, variables=variables)\n",
    "eq1 = Relation(x1 + 2*x3 + x4 - 4, variables=variables)\n",
    "eq2 = Relation(x2 + x3 - x4 - 2, variables=variables)\n",
    "bounds = [(0, None), (0, None), (0, None), (0, None)]\n",
    "lpp = LPP(obj, [eq1, eq2], variables=variables, bounds=bounds)\n",
    "solution = lpp.solve()\n",
    "solution.x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "solution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "type(solution)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If the bounds for the variables are all the same, you can afford to input a single bound:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x1, x2, x3, x4 = variables = sy.symbols(['x1', 'x2', 'x3', 'x4'])\n",
    "obj = Function(3*x1 + 9*x3 + x2 + x4, variables=variables)\n",
    "eq1 = Relation(x1 + 2*x3 + x4 - 4, variables=variables)\n",
    "eq2 = Relation(x2 + x3 - x4 - 2, variables=variables)\n",
    "bounds = (0, None)\n",
    "lpp = LPP(obj, [eq1, eq2], variables=variables, bounds=bounds)\n",
    "lpp.solve().x"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Furthermore, if all the bounds are in the form $\\geq$, you can omit the arguments for the parameter `bounds`, since it is the default value."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "variables = sy.symbols(['x1', 'x2', 'x3', 'x4'])\n",
    "obj = Function(\"3*x1 + 9*x3 + x2 + x4\", variables=variables)\n",
    "eq1 = Relation(\"x1 + 2*x3 + x4 = 4\", variables=variables)\n",
    "eq2 = Relation(\"x2 + x3 - x4 = 2\", variables=variables)\n",
    "lpp = LPP(obj, [eq1, eq2], variables=variables)\n",
    "lpp.solve().x"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2d example with `matplotlib`"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Consider the following problem:\n\n$$\n\\begin{aligned}\n    & \\text{minimize}& \\quad  x_1 + x_2 &\\\\\n    & \\text{subject to} & & \\\\\n    & & x_1 &\\geq 1, \\\\\n    & & x_2 &\\geq 1, \\\\\n    & & x_1 + x_2 &\\leq 4, \\\\\n    & & x_i &\\geq 0, \\qquad i=1, 2.\n\\end{aligned}\n$$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Setting up the problem is as follows."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "variables = sy.symbols(['x1', 'x2'])\n",
    "f = Function(\"x1 + x2\", variables=variables)\n",
    "ieq1 = Relation(\"x1 >= 1\", variables=variables)\n",
    "ieq2 = Relation(\"x2 >= 1\", variables=variables)\n",
    "ieq3 = Relation(\"x1 + x2 <= 4\", variables=variables)\n",
    "bounds = [(0, None), (0, None)]  # x1 >= 0, x2 >= 0\n",
    "lpp = LPP(f, [ieq1, ieq2, ieq3], variables=variables, bounds=bounds)\n",
    "x = lpp.solve().x\n",
    "print(\"The optimal solution is x1 = {0}, x2 = {1}.\".format(x[0], x[1]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now plotting with `Matplotlib`, with the feasible side of inequalities visualized by hatching, the ticks being on the infeasible side."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "nbsphinx-thumbnail"
    ]
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib import patheffects\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(6, 6))\n",
    "\n",
    "nx = 100\n",
    "ny = 100\n",
    "xvec = np.linspace(0.0, 4.0, nx)\n",
    "yvec = np.linspace(0.0, 4.0, ny)\n",
    "x1, x2 = np.meshgrid(xvec, yvec)\n",
    "\n",
    "obj = x1 + x2\n",
    "g1 = x1 - 1\n",
    "g2 = x2 - 1\n",
    "g3 = x1 + x2 - 4\n",
    "\n",
    "cntr = ax.contour(x1, x2, obj, colors='black')\n",
    "ax.clabel(cntr, fmt=\"%2.1f\", use_clabeltext=True)\n",
    "\n",
    "cg1 = ax.contour(x1, x2, g1, [0], colors='sandybrown')\n",
    "cg1.set(path_effects=[patheffects.withTickedStroke(angle=-135)])\n",
    "\n",
    "cg2 = ax.contour(x1, x2, g2, [0], colors='orangered')\n",
    "cg2.set(path_effects=[patheffects.withTickedStroke(angle=-60)])\n",
    "\n",
    "cg3 = ax.contour(x1, x2, g3, [0], colors='mediumblue')\n",
    "cg3.set(path_effects=[patheffects.withTickedStroke()])\n",
    "\n",
    "ax.plot(x[0], x[1], 'ms', markersize=15)\n",
    "\n",
    "ax.set_xlim(0, 4)\n",
    "ax.set_ylim(0, 4)\n",
    "\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Linear Mixed-Integer Programs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If some of the design variables are constrained to integers, you can use the `integrality` parameter. As an example, consider the following problem:\n\n$$\n\\begin{aligned}\n    & \\text{minimize}& \\quad  3 x_2 + 2 x_3 &\\\\\n    & \\text{subject to} & & \\\\\n    & & 2 x_1 + 2 x_2 - 4 x_3 &= 5, \\\\\n    & & x_i &\\geq 0, \\qquad i=1, 2, 3, \\\\\n    & & x_i &\\in \\mathbb{Z}, \\qquad i=1, 3.\n\\end{aligned}\n$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "variables = x1, x2, x3 = sy.symbols([\"x1\", \"x2\", \"x3\"])\n",
    "obj = Function(3 * x2 + 2 * x3, variables=variables)\n",
    "eq = Relation(2 * x1 + 2 * x2 - 4 * x3 - 5, op=\"=\", variables=variables)\n",
    "bounds = (0, None)\n",
    "integrality = [1, 0, 1]\n",
    "lpp = LPP(obj, [eq], variables=variables, bounds=bounds, integrality=integrality)\n",
    "solution = lpp.solve()\n",
    "solution.x"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Another way to express the same integrality constraints is using assumptions when creating the SymPy variables of the problem."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x1, x3 = sy.symbols([\"x1\", \"x3\"], integer=True)\n",
    "x2 = sy.symbols(\"x2\")\n",
    "variables = x1, x2, x3\n",
    "obj = Function(3 * x2 + 2 * x3, variables=variables)\n",
    "eq = Relation(2 * x1 + 2 * x2 - 4 * x3 - 5, op=\"=\", variables=variables)\n",
    "bounds = (0, None)\n",
    "lpp = LPP(obj, [eq], variables=variables, bounds=bounds)\n",
    "solution = lpp.solve()\n",
    "solution.x"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The solver in SciPy is capable of understanding more complex integrality constraints, colsult their API Reference for more details [here](https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.linprog.html#scipy.optimize.linprog). If you want to impose constraints like these, you must use the `integrality` parameter."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "sigmaepsilon-math-NMVzm02N-py3.12",
   "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.12.7"
  },
  "widgets": {
   "application/vnd.jupyter.widget-state+json": {
    "state": {},
    "version_major": 2,
    "version_minor": 0
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}