lp_intro.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "019bf8b2",
   "metadata": {},
   "source": [
    "# Linear Programming"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c6d0bc44",
   "metadata": {},
   "source": [
    "## Overview\n",
    "\n",
    "**Linear programming** problems either maximize or minimize\n",
    "a linear objective function subject to a set of  linear equality and/or inequality constraints.\n",
    "\n",
    "Linear programs come in pairs:\n",
    "\n",
    "- an original  **primal** problem, and  \n",
    "- an associated **dual** problem.  \n",
    "\n",
    "\n",
    "If a primal problem involves **maximization**, the dual problem involves **minimization**.\n",
    "\n",
    "If a primal problem involves  **minimization**, the dual problem involves **maximization**.\n",
    "\n",
    "We provide a standard form of a linear program and methods to transform other forms of linear programming problems  into a standard form.\n",
    "\n",
    "We tell how to solve a linear programming problem using [SciPy](https://scipy.org/).\n",
    "\n",
    "We describe the important concept of complementary slackness and how it relates to the dual problem.\n",
    "\n",
    "Let’s start with some standard imports."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d15c8775",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from scipy.optimize import linprog\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib.patches import Polygon"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6ea11fa0",
   "metadata": {},
   "source": [
    "## Objective Function and Constraints\n",
    "\n",
    "We want to minimize a **cost function** $ c'x = \\sum_{i=1}^n c_i x_i $ over  feasible values of $ x = (x_1,x_2,\\dots,x_n)' $.\n",
    "\n",
    "Here\n",
    "\n",
    "- $ c = (c_1,c_2,\\dots,c_n)' $ is  a **unit cost vector**,  and  \n",
    "- $ x = (x_1,x_2,\\dots,x_n)' $ is a vector of **decision variables**  \n",
    "\n",
    "\n",
    "Decision variables are  restricted to satisfy a set of linear equality and/or inequality constraints.\n",
    "\n",
    "We describe the constraints with the following  collections of  $ n $-dimensional vectors $ a_i $ and scalars $ b_i $  and associated sets indexing the equality and inequality constraints:\n",
    "\n",
    "- $ a_i $ for $ i \\in M_i $, where $ M_1,M_2,M_3 $ are each  sets of indexes  \n",
    "\n",
    "\n",
    "and a collection of  scalers\n",
    "\n",
    "- $ b_i $ for $ i \\in N_i $, where $ N_1,N_2,N_3 $  are each sets of indexes.  \n",
    "\n",
    "\n",
    "A linear programming can be stated as [[Bertsimas, 1997](https://python.quantecon.org/zreferences.html#id43)]:\n",
    "\n",
    "\n",
    "<a id='equation-linprog'></a>\n",
    "$$\n",
    "\\begin{aligned}\n",
    "\\min_{x} \\ & c' x \\\\\n",
    "\\mbox{subject to } \\ & a_i' x \\ge b_i, & i \\in M_1 \\\\\n",
    "& a_i' x \\le b_i, & i \\in M_2  \\\\\n",
    "& a_i' x = b_i, & i \\in M_3  \\\\\n",
    "& x_j \\ge 0, & j \\in N_1 \\\\\n",
    "& x_j \\le 0, & j \\in N_2 \\\\\n",
    "& x_j\\ \\text{unrestricted}, & j \\in N_3 \\\\\n",
    "\\end{aligned} \\tag{20.1}\n",
    "$$\n",
    "\n",
    "A vector $ x $ that satisfies all of the constraints is called a **feasible solution**.\n",
    "\n",
    "A collection of all feasible solutions is called  a **feasible set**.\n",
    "\n",
    "A feasible solution $ x $ that minimizes the cost function  is called an **optimal solution**.\n",
    "\n",
    "The corresponding value of cost function $ c'x $ is called the  **optimal value**.\n",
    "\n",
    "If the feasible set is empty, we say that solving  the  linear programming problem is **infeasible**.\n",
    "\n",
    "If, for any $ K \\in \\mathbb R $, there exists a feasible solution $ x $ such that $ c'x < K $, we say that the problem is **unbounded** or equivalently that the optimal value is $ -\\infty $."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3e61dd68",
   "metadata": {},
   "source": [
    "## Example 1: Production Problem\n",
    "\n",
    "This example was created by [[Bertsimas, 1997](https://python.quantecon.org/zreferences.html#id43)]\n",
    "\n",
    "Suppose that a factory can produce two goods called Product $ 1 $ and Product $ 2 $.\n",
    "\n",
    "To produce each product requires both material and labor.\n",
    "\n",
    "Selling each product generates revenue.\n",
    "\n",
    "Required per unit material and labor  inputs and  revenues  are shown in table below:\n",
    "\n",
    "||Product 1|Product 2|\n",
    "|:-------------------------------:|:-------------------------------:|:-------------------------------:|\n",
    "|Material|2|5|\n",
    "|Labor|4|2|\n",
    "|Revenue|3|4|\n",
    "30 units of material and 20 units of labor available.\n",
    "\n",
    "A firm’s problem is to construct a  production plan that uses its  30 units of materials and 20 unites of labor\n",
    "to maximize its revenue.\n",
    "\n",
    "Let $ x_i $ denote the quantity of Product $ i $ that the firm produces.\n",
    "\n",
    "This problem can be formulated as:\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "\\max_{x_1,x_2} \\ & z = 3 x_1 + 4 x_2 \\\\\n",
    "\\mbox{subject to } \\ & 2 x_1 + 5 x_2 \\le 30 \\\\\n",
    "& 4 x_1 + 2 x_2 \\le 20 \\\\\n",
    "& x_1, x_2 \\ge 0 \\\\\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "The following graph illustrates the firm’s constraints and iso-revenue lines."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "74ddc1c5",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=(8, 6))\n",
    "ax.grid()\n",
    "\n",
    "# Draw constraint lines\n",
    "ax.hlines(0, -1, 17.5)\n",
    "ax.vlines(0, -1, 12)\n",
    "ax.plot(np.linspace(-1, 17.5, 100), 6-0.4*np.linspace(-1, 17.5, 100), color=\"c\")\n",
    "ax.plot(np.linspace(-1, 5.5, 100), 10-2*np.linspace(-1, 5.5, 100), color=\"c\")\n",
    "ax.text(1.5, 8, \"$2x_1 + 5x_2 \\leq 30$\", size=12)\n",
    "ax.text(10, 2.5, \"$4x_1 + 2x_2 \\leq 20$\", size=12)\n",
    "ax.text(-2, 2, \"$x_2 \\geq 0$\", size=12)\n",
    "ax.text(2.5, -0.7, \"$x_1 \\geq 0$\", size=12)\n",
    "\n",
    "# Draw the feasible region\n",
    "feasible_set = Polygon(np.array([[0, 0],\n",
    "                                 [0, 6],\n",
    "                                 [2.5, 5],\n",
    "                                 [5, 0]]),\n",
    "                       color=\"cyan\")\n",
    "ax.add_patch(feasible_set)\n",
    "\n",
    "# Draw the objective function\n",
    "ax.plot(np.linspace(-1, 5.5, 100), 3.875-0.75*np.linspace(-1, 5.5, 100), color=\"orange\")\n",
    "ax.plot(np.linspace(-1, 5.5, 100), 5.375-0.75*np.linspace(-1, 5.5, 100), color=\"orange\")\n",
    "ax.plot(np.linspace(-1, 5.5, 100), 6.875-0.75*np.linspace(-1, 5.5, 100), color=\"orange\")\n",
    "ax.arrow(-1.6, 5, 0, 2, width = 0.05, head_width=0.2, head_length=0.5, color=\"orange\")\n",
    "ax.text(5.7, 1, \"$z = 3x_1 + 4x_2$\", size=12)\n",
    "\n",
    "# Draw the optimal solution\n",
    "ax.plot(2.5, 5, \"*\", color=\"black\")\n",
    "ax.text(2.7, 5.2, \"Optimal Solution\", size=12)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3711fdf2",
   "metadata": {},
   "source": [
    "The blue region is the feasible set within which all constraints are satisfied.\n",
    "\n",
    "Parallel orange lines are iso-revenue lines.\n",
    "\n",
    "The firm’s objective is to find the  parallel orange lines to the upper boundary of the feasible set.\n",
    "\n",
    "The intersection of the feasible set and the highest orange line delineates the optimal set.\n",
    "\n",
    "In this example, the optimal set is the point $ (2.5, 5) $."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5ae6cc82",
   "metadata": {},
   "source": [
    "## Example 2: Investment Problem\n",
    "\n",
    "We now consider a problem posed and solved by  [[Hu, 2018](https://python.quantecon.org/zreferences.html#id44)].\n",
    "\n",
    "A mutual fund has $ \\\\$ 100,000 $ to be invested over a three year horizon.\n",
    "\n",
    "Three investment options are available:\n",
    "\n",
    "1. **Annuity:**  the fund can  pay a same amount of new capital at the beginning of each of three years and receive a payoff of 130% of **total capital** invested  at the end of the third year. Once the mutual fund decides to invest in this annuity, it has to keep investing in all subsequent  years in the three year horizon.  \n",
    "1. **Bank account:** the fund can deposit any amount  into a bank at the beginning of each year and receive its capital plus 6% interest at the end of that year. In addition, the mutual fund is permitted to borrow no more than \\$20,000 at the beginning of each year and is asked to pay back the amount borrowed plus 6% interest at the end of the year. The mutual fund can choose whether to deposit or borrow at the beginning of each year.  \n",
    "1. **Corporate bond:** At the beginning of the second year, a  corporate bond becomes available.\n",
    "  The fund can buy an amount\n",
    "  that is no more than $ \\\\$ $50,000 of this bond at the beginning of the second year and  at the end of the third year receive a payout of 130% of the amount invested in the bond.  \n",
    "\n",
    "\n",
    "The mutual fund’s objective is to maximize total payout that it owns at the end of the third year.\n",
    "\n",
    "We can formulate this  as a linear programming problem.\n",
    "\n",
    "Let  $ x_1 $ be the amount of put in the annuity, $ x_2, x_3, x_4 $ be  bank deposit balances at the beginning of the three years,  and $ x_5 $ be the amount invested  in the corporate bond.\n",
    "\n",
    "When $ x_2, x_3, x_4 $ are negative, it means that  the mutual fund has borrowed from  bank.\n",
    "\n",
    "The table below shows the mutual fund’s decision variables together with the timing protocol described above:\n",
    "\n",
    "||Year 1|Year 2|Year 3|\n",
    "|:-----------------------:|:-----------------------:|:-----------------------:|:-----------------------:|\n",
    "|Annuity|$ x_1 $|$ x_1 $|$ x_1 $|\n",
    "|Bank account|$ x_2 $|$ x_3 $|$ x_4 $|\n",
    "|Corporate bond|0|$ x_5 $|0|\n",
    "The  mutual fund’s decision making proceeds according to the following timing protocol:\n",
    "\n",
    "1. At the beginning of the first year, the mutual fund decides how much to invest in the annuity and\n",
    "  how much to deposit in the bank. This decision is subject to the constraint:  \n",
    "  $$\n",
    "  x_1 + x_2 = 100,000\n",
    "  $$\n",
    "1. At the beginning of the second year, the mutual fund has a bank balance  of $ 1.06 x_2 $.\n",
    "  It must keep $ x_1 $ in the annuity. It can choose to put $ x_5 $ into the corporate bond,\n",
    "  and put $ x_3 $ in the bank. These decisions are restricted by  \n",
    "  $$\n",
    "  x_1 + x_5 = 1.06 x_2 - x_3\n",
    "  $$\n",
    "1. At the beginning of the third year, the mutual fund has a bank account balance equal\n",
    "  to $ 1.06 x_3 $. It must again invest  $ x_1 $ in the annuity,\n",
    "  leaving it with  a bank account balance equal to $ x_4 $. This situation is summarized by the restriction:  \n",
    "  $$\n",
    "  x_1 = 1.06 x_3 - x_4\n",
    "  $$\n",
    "\n",
    "\n",
    "The mutual fund’s objective function, i.e., its wealth at the end of the third year is:\n",
    "\n",
    "$$\n",
    "1.30 \\cdot 3x_1 + 1.06 x_4 + 1.30 x_5\n",
    "$$\n",
    "\n",
    "Thus, the mutual fund confronts the linear program:\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "\\max_{x} \\ & 1.30 \\cdot 3x_1 + 1.06 x_4 + 1.30 x_5 \\\\\n",
    "\\mbox{subject to } \\ & x_1 + x_2 = 100,000\\\\\n",
    " & x_1 - 1.06 x_2 + x_3 + x_5 = 0\\\\\n",
    " & x_1 - 1.06 x_3 + x_4 = 0\\\\\n",
    " & x_2 \\ge -20,000\\\\\n",
    " & x_3 \\ge -20,000\\\\\n",
    " & x_4 \\ge -20,000\\\\\n",
    " & x_5 \\le 50,000\\\\\n",
    " & x_j \\ge 0, \\quad j = 1,5\\\\\n",
    " & x_j \\ \\text{unrestricted}, \\quad j = 2,3,4\\\\\n",
    "\\end{aligned}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "775ff7f9",
   "metadata": {},
   "source": [
    "## Standard Form\n",
    "\n",
    "For purposes of\n",
    "\n",
    "- unifying linear programs that are initially stated in superficially different forms, and  \n",
    "- having a form that is convenient to put into black-box software packages,  \n",
    "\n",
    "\n",
    "it is useful to devote some effort to describe a **standard form**.\n",
    "\n",
    "Our standard form  is:\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "\\min_{x} \\ & c_1 x_1 + c_2 x_2 + \\dots + c_n x_n  \\\\\n",
    "\\mbox{subject to } \\ & a_{11} x_1 + a_{12} x_2 + \\dots + a_{1n} x_n = b_1 \\\\\n",
    " & a_{21} x_1 + a_{22} x_2 + \\dots + a_{2n} x_n = b_2 \\\\\n",
    " & \\quad \\vdots \\\\\n",
    " & a_{m1} x_1 + a_{m2} x_2 + \\dots + a_{mn} x_n = b_m \\\\\n",
    " & x_1, x_2, \\dots, x_n \\ge 0 \\\\\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "Let\n",
    "\n",
    "$$\n",
    "A = \\begin{bmatrix}\n",
    "a_{11} & a_{12} & \\dots & a_{1n} \\\\\n",
    "a_{21} & a_{22} & \\dots & a_{2n} \\\\\n",
    "  &   & \\vdots &   \\\\\n",
    "a_{m1} & a_{m2} & \\dots & a_{mn} \\\\\n",
    "\\end{bmatrix}, \\quad\n",
    "b = \\begin{bmatrix} b_1 \\\\ b_2 \\\\ \\vdots \\\\ b_m \\\\ \\end{bmatrix}, \\quad\n",
    "c = \\begin{bmatrix} c_1 \\\\ c_2 \\\\ \\vdots \\\\ c_n \\\\ \\end{bmatrix}, \\quad\n",
    "x = \\begin{bmatrix} x_1 \\\\ x_2 \\\\ \\vdots \\\\ x_n \\\\ \\end{bmatrix}. \\quad\n",
    "$$\n",
    "\n",
    "The standard form LP problem can be expressed concisely as:\n",
    "\n",
    "\n",
    "<a id='equation-lpproblem'></a>\n",
    "$$\n",
    "\\begin{aligned}\n",
    "\\min_{x} \\ & c'x \\\\\n",
    "\\mbox{subject to } \\ & Ax = b\\\\\n",
    " & x >= 0\\\\\n",
    "\\end{aligned} \\tag{20.2}\n",
    "$$\n",
    "\n",
    "Here, $ Ax = b $ means that  the $ i $-th entry of $ Ax $  equals the $ i $-th entry of $ b $ for every $ i $.\n",
    "\n",
    "Similarly, $ x >= 0 $ means that  $ x_j $ is greater than $ 0 $ for every $ j $."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "52060198",
   "metadata": {},
   "source": [
    "### Useful Transformations\n",
    "\n",
    "It is useful to know how to transform a problem that initially is not stated in the standard form into one that is.\n",
    "\n",
    "By deploying the following steps, any linear programming problem can be transformed into an  equivalent  standard form linear programming problem.\n",
    "\n",
    "1. **Objective Function:** If a problem is originally a constrained **maximization** problem, we can construct a new objective function that  is the additive inverse of the original objective function. The transformed problem is then a **minimization** problem.  \n",
    "1. **Decision Variables:** Given a variable $ x_j $ satisfying $ x_j \\le 0 $, we can introduce a new variable $ x_j' = - x_j $ and subsitute it into original problem. Given a free variable $ x_i $ with no restriction on its sign, we can introduce two new variables $ x_j^+ $ and $ x_j^- $ satisfying $ x_j^+, x_j^- \\ge 0 $ and replace $ x_j $ by $ x_j^+ - x_j^- $.  \n",
    "1. **Inequality constraints:** Given an inequality constraint $ \\sum_{j=1}^n a_{ij}x_j \\le 0 $, we can introduce a new variable $ s_i $, called a **slack variable** that satisfies $ s_i \\ge 0 $ and replace the original constraint by $ \\sum_{j=1}^n a_{ij}x_j + s_i = 0 $.  \n",
    "\n",
    "\n",
    "Let’s apply the above steps to the two examples described above."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "392570fa",
   "metadata": {},
   "source": [
    "### Example 1: Production Problem\n",
    "\n",
    "The original problem is:\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "\\max_{x_1,x_2} \\ & 3 x_1 + 4 x_2 \\\\\n",
    "\\mbox{subject to } \\ & 2 x_1 + 5 x_2 \\le 30 \\\\\n",
    "& 4 x_1 + 2 x_2 \\le 20 \\\\\n",
    "& x_1, x_2 \\ge 0 \\\\\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "This problem is equivalent to the following problem with a standard form:\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "\\min_{x_1,x_2} \\ & -(3 x_1 + 4 x_2) \\\\\n",
    "\\mbox{subject to } \\ & 2 x_1 + 5 x_2 + s_1 = 30 \\\\\n",
    "& 4 x_1 + 2 x_2 + s_2 = 20 \\\\\n",
    "& x_1, x_2, s_1, s_2 \\ge 0 \\\\\n",
    "\\end{aligned}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3f1576f2",
   "metadata": {},
   "source": [
    "### Example 2: Investment Problem\n",
    "\n",
    "The original problem is:\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "\\max_{x} \\ & 1.30 \\cdot 3x_1 + 1.06 x_4 + 1.30 x_5 \\\\\n",
    "\\mbox{subject to } \\ & x_1 + x_2 = 100,000\\\\\n",
    " & x_1 - 1.06 x_2 + x_3 + x_5 = 0\\\\\n",
    " & x_1 - 1.06 x_3 + x_4 = 0\\\\\n",
    " & x_2 \\ge -20,000\\\\\n",
    " & x_3 \\ge -20,000\\\\\n",
    " & x_4 \\ge -20,000\\\\\n",
    " & x_5 \\le 50,000\\\\\n",
    " & x_j \\ge 0, \\quad j = 1,5\\\\\n",
    " & x_j \\ \\text{unrestricted}, \\quad j = 2,3,4\\\\\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "This problem is equivalent to the following problem with a standard form:\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "\\min_{x} \\ & -(1.30 \\cdot 3x_1 + 1.06 x_4^+ - 1.06 x_4^- + 1.30 x_5) \\\\\n",
    "\\mbox{subject to } \\ & x_1 + x_2^+ - x_2^- = 100,000\\\\\n",
    " & x_1 - 1.06 (x_2^+ - x_2^-) + x_3^+ - x_3^- + x_5 = 0\\\\\n",
    " & x_1 - 1.06 (x_3^+ - x_3^-) + x_4^+ - x_4^- = 0\\\\\n",
    " & x_2^- - x_2^+ + s_1 = 20,000\\\\\n",
    " & x_3^- - x_3^+ + s_2 = 20,000\\\\\n",
    " & x_4^- - x_4^+ + s_3 = 20,000\\\\\n",
    " & x_5 + s_4 = 50,000\\\\\n",
    " & x_j \\ge 0, \\quad j = 1,5\\\\\n",
    " & x_j^+, x_j^- \\ge 0, \\quad j = 2,3,4\\\\\n",
    " & s_j \\ge 0, \\quad j = 1,2,3,4\\\\\n",
    "\\end{aligned}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "db5ed619",
   "metadata": {},
   "source": [
    "## Computations\n",
    "\n",
    "The package *scipy.optimize* provides a function ***linprog*** to solve linear programming problems with a form below:\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "\\min_{x} \\ & c' x  \\\\\n",
    "\\mbox{subject to } \\ & A_{ub}x \\le b_{ub} \\\\\n",
    " & A_{eq}x = b_{eq} \\\\\n",
    " & l \\le x \\le u \\\\\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    ">**Note**\n",
    ">\n",
    ">By default $ l = 0 $ and $ u = \\text{None} $ unless explicitly specified with the argument ‘bounds’.\n",
    "\n",
    "Let’s apply this great Python tool to solve our two example problems."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "949731bb",
   "metadata": {},
   "source": [
    "### Example 1: Production Problem\n",
    "\n",
    "The problem is:\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "\\max_{x_1,x_2} \\ & 3 x_1 + 4 x_2 \\\\\n",
    "\\mbox{subject to } \\ & 2 x_1 + 5 x_2 \\le 30 \\\\\n",
    "& 4 x_1 + 2 x_2 \\le 20 \\\\\n",
    "& x_1, x_2 \\ge 0 \\\\\n",
    "\\end{aligned}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "31326d65",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "# Construct parameters\n",
    "c_ex1 = np.array([3, 4])\n",
    "\n",
    "# Inequality constraints\n",
    "A_ex1 = np.array([[2, 5],\n",
    "                  [4, 2]])\n",
    "b_ex1 = np.array([30,20])\n",
    "\n",
    "# Solve the problem\n",
    "# we put a negative sign on the objective as linprog does minimization\n",
    "res_ex1 = linprog(-c_ex1, A_ub=A_ex1, b_ub=b_ex1)\n",
    "\n",
    "res_ex1"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6ea229a7",
   "metadata": {},
   "source": [
    "The optimal plan tells the  factory to produce 2.5 units of Product 1 and 5 units of  Product 2; that  generates a maximizing value of  revenue of 27.5.\n",
    "\n",
    "We are using the *linprog* function as a **black box**.\n",
    "\n",
    "Inside it, Python first  transforms the problem into  standard form.\n",
    "\n",
    "To do that, for each inequality constraint it generates one slack variable.\n",
    "\n",
    "Here the vector of slack variables is a two-dimensional NumPy array that  equals $ b_{ub} - A_{ub}x $.\n",
    "\n",
    "See the [official documentation](https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.linprog.html#scipy.optimize.linprog) for more details.\n",
    "\n",
    ">**Note**\n",
    ">\n",
    ">This problem is to maximize the objective, so that we need to put a minus sign in front of parameter vector c."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0a785aee",
   "metadata": {},
   "source": [
    "### Example 2: Investment Problem\n",
    "\n",
    "The problem is:\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "\\max_{x} \\ & 1.30 \\cdot 3x_1 + 1.06 x_4 + 1.30 x_5 \\\\\n",
    "\\mbox{subject to } \\ & x_1 + x_2 = 100,000\\\\\n",
    " & x_1 - 1.06 x_2 + x_3 + x_5 = 0\\\\\n",
    " & x_1 - 1.06 x_3 + x_4 = 0\\\\\n",
    " & x_2 \\ge -20,000\\\\\n",
    " & x_3 \\ge -20,000\\\\\n",
    " & x_4 \\ge -20,000\\\\\n",
    " & x_5 \\le 50,000\\\\\n",
    " & x_j \\ge 0, \\quad j = 1,5\\\\\n",
    " & x_j \\ \\text{unrestricted}, \\quad j = 2,3,4\\\\\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "Let’s solve this problem using *linprog*."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "028d0f50",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "# Construct parameters\n",
    "rate = 1.06\n",
    "\n",
    "# Objective function parameters\n",
    "c_ex2 = np.array([1.30*3, 0, 0, 1.06, 1.30])\n",
    "\n",
    "# Inequality constraints\n",
    "A_ex2 = np.array([[1,  1,  0,  0,  0],\n",
    "                  [1, -rate, 1, 0, 1],\n",
    "                  [1, 0, -rate, 1, 0]])\n",
    "b_ex2 = np.array([100000, 0, 0])\n",
    "\n",
    "# Bounds on decision variables\n",
    "bounds_ex2 = [(  0,    None),\n",
    "              (-20000, None),\n",
    "              (-20000, None),\n",
    "              (-20000, None),\n",
    "              (  0,   50000)]\n",
    "\n",
    "# Solve the problem\n",
    "res_ex2 = linprog(-c_ex2, A_eq=A_ex2, b_eq=b_ex2,\n",
    "                  bounds=bounds_ex2)\n",
    "\n",
    "res_ex2"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "be5834c7",
   "metadata": {},
   "source": [
    "Python tells us that  the best investment strategy is:\n",
    "\n",
    "1. At the beginning of the first year, the mutual fund should buy $ \\\\$24,927.75 $ of the annuity. Its bank account balance should be $ \\\\$75,072.25 $.  \n",
    "1. At the beginning of the second year, the mutual fund should buy $ \\\\$50,000 $ of the corporate bond and keep invest in the annuity. Its bank account balance should be $ \\\\$ 4,648.83 $.  \n",
    "1. At the beginning of the third year, the mutual fund should borrow $ \\\\$20,000 $ from the bank and invest in the annuity.  \n",
    "1. At the end of the third year, the mutual fund will get payouts from the annuity and corporate bond and repay its loan from the bank. At the end  it will own $ \\\\$ $141018.24, so that it’s total net  rate of return over the three periods is 41.02%.  "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "98300002",
   "metadata": {},
   "source": [
    "## Duality\n",
    "\n",
    "Associated with a  linear programming of form [(20.1)](#equation-linprog) with $ m $ constraints and $ n $ decision variables,\n",
    "there is an **dual** linear programming problem that takes the form (please see [[Bertsimas, 1997](https://python.quantecon.org/zreferences.html#id43)])\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "\\max_{p} \\ & b' p \\\\\n",
    "\\mbox{subject to } \\ & p_i \\ge 0, & i \\in M_1 \\\\\n",
    "& p_i \\le 0, & i \\in M_2  \\\\\n",
    "& p_i\\ \\text{unrestricted}, & i \\in M_3  \\\\\n",
    "& A_j' p \\le c_j, & j \\in N_1 \\\\\n",
    "& A_j' p \\ge c_j, & j \\in N_2 \\\\\n",
    "& A_j' p = c_j, & j \\in N_3 \\\\\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "Where $ A_j $ is $ j $-th column of the $ m $ by $ n $ matrix $ A $.\n",
    "\n",
    ">**Note**\n",
    ">\n",
    ">In what  follows, we shall  use $ a_i' $ to denote the $ i $-th row of $ A $ and $ A_j $ to denote the $ j $-th column of $ A $.\n",
    "\n",
    "$$\n",
    "A = \\begin{bmatrix}\n",
    "a_1' \\\\\n",
    "a_2' \\\\\n",
    "\\ \\\\\n",
    "a_m' \\\\\n",
    "\\end{bmatrix}.\n",
    "$$\n",
    "\n",
    "To construct the  dual of linear programming problem [(20.1)](#equation-linprog), we proceed as follows:\n",
    "\n",
    "1. For every constraint $ a_i' x \\ge(\\le or =) b_i $, $ j = 1,2,...,m $, in the primal problem, we construct a corresponding dual variable $ p_i $. $ p_i $ is restricted to be positive if $ a_i' x \\ge b_i $ or negative if $ a_i' x \\le b_i $ or unrestricted if $ a_i' x = b_i $. We construct the $ m $-dimensional vector   $ p $ with entries $ p_i $.  \n",
    "1. For every variable $ x_j $, $ j = 1,2,...,n $, we construct a corresponding dual constraint $ A_j' p \\ge(\\le or =) c_j $. The constraint is $ A_j' p \\ge c_j $ if $ x_j \\le 0 $, $ A_j' p \\le c_j $ if $ x_j \\ge 0 $ or $ A_j' p = c_j $ if $ x_j $ is unrestricted.  \n",
    "1. The dual problem is to **maximize**  objective function $ b'p $.  \n",
    "\n",
    "\n",
    "For a **maximization**  problem, we can first transform it to an equivalent minimization problem and then follow the above steps above to construct the  dual **minimization** problem.\n",
    "\n",
    "We can easily verify that  **the dual of a dual problem is the  primal problem**.\n",
    "\n",
    "The following table summarizes relationships between objects in primal and dual problems.\n",
    "\n",
    "|Objective: Min|Objective: Max|\n",
    "|:------------------------------------------------:|:------------------------------------------------:|\n",
    "|m constraints|m variables|\n",
    "|constraint $ \\ge $|variable $ \\ge $ 0|\n",
    "|constraint $ \\le $|variable $ \\le $ 0|\n",
    "|constraint $ = $|variable free|\n",
    "|n variables|n constraints|\n",
    "|variable $ \\ge $ 0|constraint $ \\le $|\n",
    "|variable $ \\le $ 0|constraint $ \\ge $|\n",
    "|variable free|constraint $ = $|\n",
    "As an example, the dual problem of the standard form [(20.2)](#equation-lpproblem) is:\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "\\max_{p} \\ & b'p \\\\\n",
    "\\mbox{subject to } \\ & A'p \\le c\\\\\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "As another example, consider a linear programming problem with form:\n",
    "\n",
    "\n",
    "<a id='equation-linprog2'></a>\n",
    "$$\n",
    "\\begin{aligned}\n",
    "\\max_{x} \\ & c'x \\\\\n",
    "\\mbox{subject to } \\ & A x \\le b\\\\\n",
    "& x \\ge 0\\\\\n",
    "\\end{aligned} \\tag{20.3}\n",
    "$$\n",
    "\n",
    "Its dual problem is:\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "\\min_{p} \\ & b'p \\\\\n",
    "\\mbox{subject to } \\ & A' p \\ge c\\\\\n",
    "& p \\ge 0\\\\\n",
    "\\end{aligned}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "76519ff0",
   "metadata": {},
   "source": [
    "## Duality Theorems\n",
    "\n",
    "Primal and dual problems are linked by powerful **duality theorems** that have  **weak** and **strong** forms.\n",
    "\n",
    "The duality theorems provide the foundations of enlightening economic interpretations of linear programming problems.\n",
    "\n",
    "**Weak duality:** For linear programming problem [(20.1)](#equation-linprog), if $ x $ and $ p $ are feasible solutions to the primal and the dual problems, respectively, then\n",
    "\n",
    "$$\n",
    "b'p \\le c'x\n",
    "$$\n",
    "\n",
    "**Strong duality:** For linear programming problem [(20.1)](#equation-linprog), if the primal problem has an optimal solution $ x $, then the dual problem also has an optimal solution. Denote an optimal solution of the dual problem as $ p $. Then\n",
    "\n",
    "$$\n",
    "b'p = c'x\n",
    "$$\n",
    "\n",
    "According to strong duality, we can find the optimal value for the primal problem by solving the dual problem.\n",
    "\n",
    "But the dual problem tells us even more as we shall see next."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6d8b778f",
   "metadata": {},
   "source": [
    "### Complementary Slackness\n",
    "\n",
    "Let $ x $ and $ p $ be feasible solutions to the primal problem [(20.1)](#equation-linprog) and its dual problem, respectively.\n",
    "\n",
    "Then $ x $ and $ p $ are also  optimal solutions of the primal and dual problems if and only if:\n",
    "\n",
    "$$\n",
    "p_i (a_i' x - b_i) = 0, \\quad \\forall i, \\\\\n",
    "x_j (A_j' p - c_j) = 0, \\quad \\forall j.\n",
    "$$\n",
    "\n",
    "This means that $ p_i = 0 $ if $ a_i' x - b_i \\neq 0 $ and $ x_j = 0 $ if $ A_j' p - c_j \\neq 0 $.\n",
    "\n",
    "These are the celebrated **complementary slackness** conditions.\n",
    "\n",
    "Let’s interpret them."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dd3b5874",
   "metadata": {},
   "source": [
    "### Interpretations\n",
    "\n",
    "Let’s take a version of problem [(20.3)](#equation-linprog2) as a production problem and consider its associated dual problem.\n",
    "\n",
    "A factory produce $ n $ products with $ m $ types of resources.\n",
    "\n",
    "Where  $ i=1,2,\\dots,m $ and $ j=1,2,\\dots,n $, let\n",
    "\n",
    "- $ x_j $ denote quantities of product $ j $ to be produced  \n",
    "- $ a_{ij} $ denote required amount of resource $ i $ to make one unit of product $ j $,  \n",
    "- $ b_i $ denotes the  avaliable amount of resource $ i $  \n",
    "- $ c_j $ denotes the revenue generated by producing one unit of product $ j $.  \n",
    "\n",
    "\n",
    "**Dual variables:** By  strong duality, we have\n",
    "\n",
    "$$\n",
    "c_1 x_1 + c_2 x_2 + \\dots + c_n x_n = b_1 p_1 + b_2 p_2 + \\dots + b_m p_m.\n",
    "$$\n",
    "\n",
    "Evidently, a one unit change of $ b_i $ results in $ p_i $ units  change of revenue.\n",
    "\n",
    "Thus, a dual variable can be interpreted as the **value** of one unit of resource $ i $.\n",
    "\n",
    "This is why it is often called the  **shadow price** of resource $ i $.\n",
    "\n",
    "For feasible but not optimal primal and dual solutions $ x $ and $ p $, by weak duality, we have\n",
    "\n",
    "$$\n",
    "c_1 x_1 + c_2 x_2 + \\dots + c_n x_n < b_1 p_1 + b_2 p_2 + \\dots + b_m p_m.\n",
    "$$\n",
    "\n",
    ">**Note**\n",
    ">\n",
    ">Here, the expression is opposite to the statement above since primal problem is a minimization problem.\n",
    "\n",
    "When a strict inequality holds, the solution is not optimal because  it doesn’t fully utilize all valuable resources.\n",
    "\n",
    "Evidently,\n",
    "\n",
    "- if a shadow price $ p_i $ is larger than the market price for Resource $ i $, the factory should buy more Resource $ i $ and expand its scale to generate more revenue;  \n",
    "- if a shadow price $ p_i $ is less than the market price for Resource $ i $, the factory should sell its  Resource $ i $.  \n",
    "\n",
    "\n",
    "**Complementary slackness:** If there exists $ i $ such that $ a_i' x - b_i < 0 $ for some $ i $, then $ p_i = 0 $ by complementary slackness. $ a_i' x - b_i < 0 $ means that  to achieve its optimal production, the factory doesn’t require as much  Resource $ i $ as it has.\n",
    "It is reasonable that the shadow price of Resource $ i $ is  0:  some of its resource $ i $ is redundant.\n",
    "\n",
    "If there exists $ j $ such that $ A_j' p - c_j > 0 $, then $ x_j = 0 $ by complementary slackness. $ A_j' p - c_j > 0 $ means that  the value of all resources used when producing one unit of product $ j $ is greater than its cost.\n",
    "\n",
    "This means that  producing another product that  can more efficiently utilize these resources is a better choice than producing product $ j $\n",
    "\n",
    "Since producing product $ j $ is not optimal, $ x_j $ should equal $ 0 $."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f8ed8aac",
   "metadata": {},
   "source": [
    "### Example 1: Production Problem\n",
    "\n",
    "This problem is one specific instance of the problem [(20.3)](#equation-linprog2), whose economic meaning is interpreted above.\n",
    "\n",
    "Its dual problem is:\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "\\min_{x_1,x_2} \\ & 30 p_1 + 20 p_2 \\\\\n",
    "\\mbox{subject to } \\ & 2 p_1 + 4 p_2 \\ge 3 \\\\\n",
    "& 5 p_1 + 2 p_2 \\ge 4 \\\\\n",
    "& p_1, p_2 \\ge 0 \\\\\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "We solve this dual problem by using the  function *linprog*.\n",
    "\n",
    "Since parameters used here are defined before when solving the primal problem, we won’t define them here."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f647aafc",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "# Solve the dual problem\n",
    "res_ex1_dual = linprog(b_ex1, A_ub=-A_ex1.T, b_ub=-c_ex1)\n",
    "\n",
    "res_ex1_dual"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "315c01a8",
   "metadata": {},
   "source": [
    "The optimal value for the dual problem equals  27.5.\n",
    "\n",
    "This equals the optimal value of  the primal problem, an illustration of  strong duality.\n",
    "\n",
    "Shadow prices for materials and labor are 0.625 and 0.4375, respectively."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f21bd4e9",
   "metadata": {},
   "source": [
    "### Example 2: Investment Problem\n",
    "\n",
    "The dual problem is:\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "\\min_{p} \\ & 100,000 p_1 - 20,000 p_4 - 20,000 p_5 - 20,000 p_6 + 50,000 p_7 \\\\\n",
    "\\mbox{subject to } \\ & p_1 + p_2 + p_3 \\ge 1.30 \\cdot 3 \\\\\n",
    "& p_1 - 1.06 p_2 + p_4 = 0 \\\\\n",
    "& p_2 - 1.06 p_3 + p_5 = 0 \\\\\n",
    "& p_3 + p_6 = 1.06 \\\\\n",
    "& p_2 + p_7 \\ge 1.30 \\\\\n",
    "& p_i \\ \\text{unrestricted}, \\quad i = 1,2,3 \\\\\n",
    "& p_i \\le 0, \\quad i = 4,5,6 \\\\\n",
    "& p_7 \\ge 0 \\\\\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "We solve this dual problem by using the function *linprog*."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5109b2e2",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "# Objective function parameters\n",
    "c_ex2_dual = np.array([100000, 0, 0, -20000, -20000, -20000, 50000])\n",
    "\n",
    "# Equality constraints\n",
    "A_eq_ex2_dual = np.array([[1, -1.06,     0,  1,  0,  0,  0],\n",
    "                          [0,     1, -1.06,  0,  1,  0,  0],\n",
    "                          [0,     0,     1,  0,  0,  1,  0]])\n",
    "b_eq_ex2_dual = np.array([0, 0, 1.06])\n",
    "\n",
    "# Inequality constraints\n",
    "A_ub_ex2_dual = - np.array([[1, 1, 1, 0, 0, 0, 0],\n",
    "                            [0, 1, 0, 0, 0, 0, 1]])\n",
    "b_ub_ex2_dual = - np.array([1.30*3, 1.30])\n",
    "\n",
    "# Bounds on decision variables\n",
    "bounds_ex2_dual = [(None, None),\n",
    "                   (None, None),\n",
    "                   (None, None),\n",
    "                   (None,    0),\n",
    "                   (None,    0),\n",
    "                   (None,    0),\n",
    "                   (   0, None)]\n",
    "\n",
    "# Solve the dual problem\n",
    "res_ex2_dual = linprog(c_ex2_dual, A_eq=A_eq_ex2_dual, b_eq=b_eq_ex2_dual,\n",
    "                       A_ub=A_ub_ex2_dual, b_ub=b_ub_ex2_dual, bounds=bounds_ex2_dual)\n",
    "\n",
    "res_ex2_dual"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c51164b7",
   "metadata": {},
   "source": [
    "The optimal value for the dual problem is 141018.24, which equals the value of the primal problem.\n",
    "\n",
    "Now, let’s interpret the dual variables.\n",
    "\n",
    "By strong duality and also our  numerical results, we have that  optimal value is:\n",
    "\n",
    "$$\n",
    "100,000 p_1 - 20,000 p_4 - 20,000 p_5 - 20,000 p_6 + 50,000 p_7.\n",
    "$$\n",
    "\n",
    "We know if $ b_i $ changes one dollor, then the optimal payoff in the end of the third year will change $ p_i $ dollars.\n",
    "\n",
    "For $ i = 1 $, this means if the initial capital changes by one dollar, then the optimal payoff in the end of the third year will change $ p_1 $ dollars.\n",
    "\n",
    "Thus, $ p_1 $ is the potential value of one more unit of initial capital, or the shadow price for initial capital.\n",
    "\n",
    "We can also interpret $ p_1 $ as the prospective  value in the end of the third year coming from having  one more  dollar to invest at the beginning of the first year.\n",
    "\n",
    "If the mutual fund can raise money at a cost lower than $ p_1 - 1 $, then it should raise more money to increase its revenue.\n",
    "\n",
    "But if it bears a cost of funds higher than $ p_1 - 1 $, the mutual fund shouldn’t do that.\n",
    "\n",
    "For $ i = 4, 5, 6 $, this means that if the amount of capital that the fund is permitted to borrow from the bank changes by one dollar, the optimal pay out  at the end of the third year will change $ p_i $ dollars.\n",
    "\n",
    "Thus, for $ i = 4, 5, 6 $, $ |p_i| $ indicates the value of one dollar that the mutual fund can borrow from the bank at the beginning of the $ i-3 $-th year.\n",
    "\n",
    "$ |p_i| $ is the shadow price for the loan amount. (We use absolute value here since $ p_i \\le 0 $.)\n",
    "\n",
    "If the interest rate is lower than $ |p_i| $, then the mutual fund should borrow to increase its optimal payoff;  if the interest rate is higher, it is better to not do this.\n",
    "\n",
    "For $ i = 7 $, this means that if the amount of the corporate bond the mutual fund can buy changes one dollar, then the optimal payoff will change $ p_7 $ dollars at the end of the third year.  Again, $ p_7 $ is the shadow price for the amount of the corporate bond the mutual fund can buy.\n",
    "\n",
    "As for numerical results\n",
    "\n",
    "1. $ p_1 = 1.38 $, which means one dollar of initial capital is worth $ \\\\$ 1.38 $ at the end of the third year.  \n",
    "1. $ p_4 = p_5 = 0 $, which means the loan amounts at the beginning of the first and second year are worth nothing. Recall that the optimal solution to the primal problem, $ x_2, x_3 > 0 $, which means at the beginning of the first and second year, the mutual fund has a postive bank account and borrows no capital from the bank. Thus, it is reasonable that the loan amounts at the beginning of the first and second year are valueless. This is what the  complementary slackness conditions mean in this setting.  \n",
    "1. $ p_6 = -0.16 $, which means one dollar of the loan amount at the beginning of the third year is worth $ \\\\$ 0.16 $. Since $ |p_6| $ is higher than the interest rate 6%, the mutual fund should borrow as much as possible at the beginning of the third year. Recall that  the optimal solution to the primal problem is  $ x_4 = -20,000 $ which means the mutual fund borrows money from the bank as much as it can.  \n",
    "1. $ p_7 = 0.0015 $, which means one dollar of the amount of the corporate bond that the mutual fund can buy is worth $ \\\\$ 0.0015 $.  "
   ]
  }
 ],
 "metadata": {
  "date": 1710733973.1921318,
  "filename": "lp_intro.md",
  "kernelspec": {
   "display_name": "Python",
   "language": "python3",
   "name": "python3"
  },
  "title": "Linear Programming"
 },
 "nbformat": 4,
 "nbformat_minor": 5
}