{
"cells": [
{
"cell_type": "markdown",
"id": "804a9656",
"metadata": {},
"source": [
"\n",
"<a id='lp-intro'></a>"
]
},
{
"cell_type": "markdown",
"id": "90cec011",
"metadata": {},
"source": [
"# Linear Programming\n",
"\n",
"In this lecture, we will need the following library. Install [ortools](https://developers.google.com/optimization) using `pip`."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "ffca306f",
"metadata": {
"hide-output": false
},
"outputs": [],
"source": [
"!pip install ortools"
]
},
{
"cell_type": "markdown",
"id": "168afdce",
"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/) and [Google OR-Tools](https://developers.google.com/optimization).\n",
"\n",
"In another lecture, we will employ the linear programming method to solve the\n",
"[optimal transport problem](https://tools-techniques.quantecon.org/opt_transport.html).\n",
"\n",
"Let’s start with some standard imports."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "48957850",
"metadata": {
"hide-output": false
},
"outputs": [],
"source": [
"import numpy as np\n",
"from ortools.linear_solver import pywraplp\n",
"from scipy.optimize import linprog\n",
"import matplotlib.pyplot as plt\n",
"from matplotlib.patches import Polygon"
]
},
{
"cell_type": "markdown",
"id": "054aa4b3",
"metadata": {},
"source": [
"Let’s start with some examples of linear programming problem."
]
},
{
"cell_type": "markdown",
"id": "4725d17d",
"metadata": {},
"source": [
"## Example 1: production problem\n",
"\n",
"This example was created by [[Bertsimas, 1997](https://intro.quantecon.org/zreferences.html#id63)]\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 units of labor to maximize its revenue.\n",
"\n",
"Let $ x_i $ denote the quantity of Product $ i $ that the firm produces and $ z $ denote the total revenue.\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.\n",
"\n",
"Iso-revenue lines show all the combinations of materials and labor that produce the same revenue."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "108fea28",
"metadata": {
"hide-output": false
},
"outputs": [],
"source": [
"fig, ax = plt.subplots()\n",
"# Draw constraint lines\n",
"ax.set_xlim(0,15)\n",
"ax.set_ylim(0,10)\n",
"x1 = np.linspace(0, 15)\n",
"ax.plot(x1, 6-0.4*x1, label=\"$2x_1 + 5x_2=30$\")\n",
"ax.plot(x1, 10-2*x1, label=\"$4x_1 + 2x_2=20$\")\n",
"\n",
"\n",
"# Draw the feasible region\n",
"feasible_set = Polygon(np.array([[0, 0],[0, 6],[2.5, 5],[5, 0]]), alpha=0.1)\n",
"ax.add_patch(feasible_set)\n",
"\n",
"# Draw the objective function\n",
"ax.plot(x1, 3.875-0.75*x1, label=\"iso-revenue lines\",color='k',linewidth=0.75)\n",
"ax.plot(x1, 5.375-0.75*x1, color='k',linewidth=0.75)\n",
"ax.plot(x1, 6.875-0.75*x1, color='k',linewidth=0.75)\n",
"\n",
"# Draw the optimal solution\n",
"ax.plot(2.5, 5, \".\", label=\"optimal solution\")\n",
"ax.set_xlabel(\"$x_1$\")\n",
"ax.set_ylabel(\"$x_2$\")\n",
"ax.legend()\n",
"\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"id": "5d1e301c",
"metadata": {},
"source": [
"The blue region is the feasible set within which all constraints are satisfied.\n",
"\n",
"Parallel black lines are iso-revenue lines.\n",
"\n",
"The firm’s objective is to find the parallel black lines to the upper boundary of the feasible set.\n",
"\n",
"The intersection of the feasible set and the highest black line delineates the optimal set.\n",
"\n",
"In this example, the optimal set is the point $ (2.5, 5) $."
]
},
{
"cell_type": "markdown",
"id": "94ebe540",
"metadata": {},
"source": [
"### Computation: using OR-Tools\n",
"\n",
"Let’s try to solve the same problem using the package `ortools.linear_solver`.\n",
"\n",
"The following cell instantiates a solver and creates two variables specifying the range of values that they can have."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "d7c5aa79",
"metadata": {
"hide-output": false
},
"outputs": [],
"source": [
"# Instantiate a GLOP(Google Linear Optimization Package) solver\n",
"solver = pywraplp.Solver.CreateSolver('GLOP')"
]
},
{
"cell_type": "markdown",
"id": "c0166c9f",
"metadata": {},
"source": [
"Let’s create two variables $ x_1 $ and $ x_2 $ such that they can only have nonnegative values."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "7c9a02bc",
"metadata": {
"hide-output": false
},
"outputs": [],
"source": [
"# Create the two variables and let them take on any non-negative value.\n",
"x1 = solver.NumVar(0, solver.infinity(), 'x1')\n",
"x2 = solver.NumVar(0, solver.infinity(), 'x2')"
]
},
{
"cell_type": "markdown",
"id": "8b78dbe8",
"metadata": {},
"source": [
"Add the constraints to the problem."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "c3d0b797",
"metadata": {
"hide-output": false
},
"outputs": [],
"source": [
"# Constraint 1: 2x_1 + 5x_2 <= 30.0\n",
"solver.Add(2 * x1 + 5 * x2 <= 30.0)\n",
"\n",
"# Constraint 2: 4x_1 + 2x_2 <= 20.0\n",
"solver.Add(4 * x1 + 2 * x2 <= 20.0)"
]
},
{
"cell_type": "markdown",
"id": "0af01e93",
"metadata": {},
"source": [
"Let’s specify the objective function. We use `solver.Maximize` method in the case when we want to maximize the objective function and in the case of minimization we can use `solver.Minimize`."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "ef1bdc12",
"metadata": {
"hide-output": false
},
"outputs": [],
"source": [
"# Objective function: 3x_1 + 4x_2\n",
"solver.Maximize(3 * x1 + 4 * x2)"
]
},
{
"cell_type": "markdown",
"id": "7390e422",
"metadata": {},
"source": [
"Once we solve the problem, we can check whether the solver was successful in solving the problem using its status. If it’s successful, then the status will be equal to `pywraplp.Solver.OPTIMAL`."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "143a3772",
"metadata": {
"hide-output": false
},
"outputs": [],
"source": [
"# Solve the system.\n",
"status = solver.Solve()\n",
"\n",
"if status == pywraplp.Solver.OPTIMAL:\n",
" print('Objective value =', solver.Objective().Value())\n",
" print(f'(x1, x2): ({x1.solution_value():.2}, {x2.solution_value():.2})')\n",
"else:\n",
" print('The problem does not have an optimal solution.')"
]
},
{
"cell_type": "markdown",
"id": "fae978a5",
"metadata": {},
"source": [
"## Example 2: investment problem\n",
"\n",
"We now consider a problem posed and solved by [[Hu, 2018](https://intro.quantecon.org/zreferences.html#id64)].\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": "0984a106",
"metadata": {},
"source": [
"### Computation: using OR-Tools\n",
"\n",
"Let’s try to solve the above problem using the package `ortools.linear_solver`.\n",
"\n",
"The following cell instantiates a solver and creates two variables specifying the range of values that they can have."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "40011430",
"metadata": {
"hide-output": false
},
"outputs": [],
"source": [
"# Instantiate a GLOP(Google Linear Optimization Package) solver\n",
"solver = pywraplp.Solver.CreateSolver('GLOP')"
]
},
{
"cell_type": "markdown",
"id": "3bcb3c2e",
"metadata": {},
"source": [
"Let’s create five variables $ x_1, x_2, x_3, x_4, $ and $ x_5 $ such that they can only have the values defined in the above constraints."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "ef3980f1",
"metadata": {
"hide-output": false
},
"outputs": [],
"source": [
"# Create the variables using the ranges available from constraints\n",
"x1 = solver.NumVar(0, solver.infinity(), 'x1')\n",
"x2 = solver.NumVar(-20_000, solver.infinity(), 'x2')\n",
"x3 = solver.NumVar(-20_000, solver.infinity(), 'x3')\n",
"x4 = solver.NumVar(-20_000, solver.infinity(), 'x4')\n",
"x5 = solver.NumVar(0, 50_000, 'x5')"
]
},
{
"cell_type": "markdown",
"id": "ab1466a9",
"metadata": {},
"source": [
"Add the constraints to the problem."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "5f2258da",
"metadata": {
"hide-output": false
},
"outputs": [],
"source": [
"# Constraint 1: x_1 + x_2 = 100,000\n",
"solver.Add(x1 + x2 == 100_000.0)\n",
"\n",
"# Constraint 2: x_1 - 1.06 * x_2 + x_3 + x_5 = 0\n",
"solver.Add(x1 - 1.06 * x2 + x3 + x5 == 0.0)\n",
"\n",
"# Constraint 3: x_1 - 1.06 * x_3 + x_4 = 0\n",
"solver.Add(x1 - 1.06 * x3 + x4 == 0.0)"
]
},
{
"cell_type": "markdown",
"id": "6c1f1733",
"metadata": {},
"source": [
"Let’s specify the objective function."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "286429f8",
"metadata": {
"hide-output": false
},
"outputs": [],
"source": [
"# Objective function: 1.30 * 3 * x_1 + 1.06 * x_4 + 1.30 * x_5\n",
"solver.Maximize(1.30 * 3 * x1 + 1.06 * x4 + 1.30 * x5)"
]
},
{
"cell_type": "markdown",
"id": "29567693",
"metadata": {},
"source": [
"Let’s solve the problem and check the status using `pywraplp.Solver.OPTIMAL`."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "69d1da0a",
"metadata": {
"hide-output": false
},
"outputs": [],
"source": [
"# Solve the system.\n",
"status = solver.Solve()\n",
"\n",
"if status == pywraplp.Solver.OPTIMAL:\n",
" print('Objective value =', solver.Objective().Value())\n",
" x1_sol = round(x1.solution_value(), 3)\n",
" x2_sol = round(x2.solution_value(), 3)\n",
" x3_sol = round(x1.solution_value(), 3)\n",
" x4_sol = round(x2.solution_value(), 3)\n",
" x5_sol = round(x1.solution_value(), 3)\n",
" print(f'(x1, x2, x3, x4, x5): ({x1_sol}, {x2_sol}, {x3_sol}, {x4_sol}, {x5_sol})')\n",
"else:\n",
" print('The problem does not have an optimal solution.')"
]
},
{
"cell_type": "markdown",
"id": "94199388",
"metadata": {},
"source": [
"OR-Tools 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.755 $ of the annuity. Its bank account balance should be $ \\\\$75,072.245 $. \n",
"1. At the beginning of the second year, the mutual fund should buy $ \\\\$24,927.755 $ of the corporate bond and keep invest in the annuity. Its bank balance should be $ \\\\$24,927.755 $. \n",
"1. At the beginning of the third year, the bank balance should be $ \\\\$75,072.245 $. \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 $ \\\\$141,018.24 $, so that it’s total net rate of return over the three periods is $ 41.02\\% $. "
]
},
{
"cell_type": "markdown",
"id": "17b7bcb7",
"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 linear programming 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 \\geq 0\\\\\n",
"\\end{aligned} \\tag{37.1}\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 \\geq 0 $ means that $ x_j $ is greater than equal to $ 0 $ for every $ j $."
]
},
{
"cell_type": "markdown",
"id": "cd6394bb",
"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 substitute 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": "21e485e2",
"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": "7f4782e1",
"metadata": {},
"source": [
"### Computation: using SciPy\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",
"$ A_{eq}, b_{eq} $ denote the equality constraint matrix and vector, and $ A_{ub}, b_{ub} $ denote the inequality constraint matrix and vector.\n",
"\n",
">**Note**\n",
">\n",
">By default $ l = 0 $ and $ u = \\text{None} $ unless explicitly specified with the argument `bounds`.\n",
"\n",
"Let’s now try to solve the Problem 1 using SciPy."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "8da8a785",
"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])"
]
},
{
"cell_type": "markdown",
"id": "363842f9",
"metadata": {},
"source": [
"Once we solve the problem, we can check whether the solver was successful in solving the problem using the boolean attribute `success`. If it’s successful, then the `success` attribute is set to `True`."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "9d58796c",
"metadata": {
"hide-output": false
},
"outputs": [],
"source": [
"# 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",
"if res_ex1.success:\n",
" # We use negative sign to get the optimal value (maximized value)\n",
" print('Optimal Value:', -res_ex1.fun)\n",
" print(f'(x1, x2): {res_ex1.x[0], res_ex1.x[1]}')\n",
"else:\n",
" print('The problem does not have an optimal solution.')"
]
},
{
"cell_type": "markdown",
"id": "05d4f5a1",
"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": "ebe47379",
"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": "code",
"execution_count": null,
"id": "f8460c20",
"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([100_000, 0, 0])\n",
"\n",
"# Bounds on decision variables\n",
"bounds_ex2 = [( 0, None),\n",
" (-20_000, None),\n",
" (-20_000, None),\n",
" (-20_000, None),\n",
" ( 0, 50_000)]"
]
},
{
"cell_type": "markdown",
"id": "0d2ffa57",
"metadata": {},
"source": [
"Let’s solve the problem and check the status using `success` attribute."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "beef5839",
"metadata": {
"hide-output": false
},
"outputs": [],
"source": [
"# Solve the problem\n",
"res_ex2 = linprog(-c_ex2, A_eq=A_ex2, b_eq=b_ex2,\n",
" bounds=bounds_ex2)\n",
"\n",
"if res_ex2.success:\n",
" # We use negative sign to get the optimal value (maximized value)\n",
" print('Optimal Value:', -res_ex2.fun)\n",
" x1_sol = round(res_ex2.x[0], 3)\n",
" x2_sol = round(res_ex2.x[1], 3)\n",
" x3_sol = round(res_ex2.x[2], 3)\n",
" x4_sol = round(res_ex2.x[3], 3)\n",
" x5_sol = round(res_ex2.x[4], 3)\n",
" print(f'(x1, x2, x3, x4, x5): {x1_sol, x2_sol, x3_sol, x4_sol, x5_sol}')\n",
"else:\n",
" print('The problem does not have an optimal solution.')"
]
},
{
"cell_type": "markdown",
"id": "25ec1e85",
"metadata": {},
"source": [
"SciPy 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 $ \\\\$141,018.24 $, so that it’s total net rate of return over the three periods is $ 41.02\\% $. \n",
"\n",
"\n",
">**Note**\n",
">\n",
">You might notice the difference in the values of optimal solution using OR-Tools and SciPy but the optimal value is the same. It is because there can be many optimal solutions for the same problem."
]
},
{
"cell_type": "markdown",
"id": "cc3cb65d",
"metadata": {},
"source": [
"## Exercises"
]
},
{
"cell_type": "markdown",
"id": "666f08bf",
"metadata": {},
"source": [
"## Exercise 37.1\n",
"\n",
"Implement a new extended solution for the Problem 1 where in the factory owner decides that number of units of Product 1 should not be less than the number of units of Product 2."
]
},
{
"cell_type": "markdown",
"id": "91d0e13e",
"metadata": {},
"source": [
"## Solution to[ Exercise 37.1](https://intro.quantecon.org/#lp_intro_ex1)\n",
"\n",
"So we can reformulate the problem 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 \\ge x_2 \\\\\n",
"& x_1, x_2 \\ge 0 \\\\\n",
"\\end{aligned}\n",
"$$"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "91b5e317",
"metadata": {
"hide-output": false
},
"outputs": [],
"source": [
"# Instantiate a GLOP(Google Linear Optimization Package) solver\n",
"solver = pywraplp.Solver.CreateSolver('GLOP')\n",
"\n",
"# Create the two variables and let them take on any non-negative value.\n",
"x1 = solver.NumVar(0, solver.infinity(), 'x1')\n",
"x2 = solver.NumVar(0, solver.infinity(), 'x2')"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "205b19b1",
"metadata": {
"hide-output": false
},
"outputs": [],
"source": [
"# Constraint 1: 2x_1 + 5x_2 <= 30.0\n",
"solver.Add(2 * x1 + 5 * x2 <= 30.0)\n",
"\n",
"# Constraint 2: 4x_1 + 2x_2 <= 20.0\n",
"solver.Add(4 * x1 + 2 * x2 <= 20.0)\n",
"\n",
"# Constraint 3: x_1 >= x_2\n",
"solver.Add(x1 >= x2)"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "a6811470",
"metadata": {
"hide-output": false
},
"outputs": [],
"source": [
"# Objective function: 3x_1 + 4x_2\n",
"solver.Maximize(3 * x1 + 4 * x2)"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "df85e591",
"metadata": {
"hide-output": false
},
"outputs": [],
"source": [
"# Solve the system.\n",
"status = solver.Solve()\n",
"\n",
"if status == pywraplp.Solver.OPTIMAL:\n",
" print('Objective value =', solver.Objective().Value())\n",
" x1_sol = round(x1.solution_value(), 2)\n",
" x2_sol = round(x2.solution_value(), 2)\n",
" print(f'(x1, x2): ({x1_sol}, {x2_sol})')\n",
"else:\n",
" print('The problem does not have an optimal solution.')"
]
},
{
"cell_type": "markdown",
"id": "7c331ed8",
"metadata": {},
"source": [
"## Exercise 37.2\n",
"\n",
"A carpenter manufactures $ 2 $ products - $ A $ and $ B $.\n",
"\n",
"Product $ A $ generates a profit of $ 23 $ and product $ B $ generates a profit of $ 10 $.\n",
"\n",
"It takes $ 2 $ hours for the carpenter to produce $ A $ and $ 0.8 $ hours to produce $ B $.\n",
"\n",
"Moreover, he can’t spend more than $ 25 $ hours per week and the total number of units of $ A $ and $ B $ should not be greater than $ 20 $.\n",
"\n",
"Find the number of units of $ A $ and product $ B $ that he should manufacture in order to maximise his profit."
]
},
{
"cell_type": "markdown",
"id": "b4583e37",
"metadata": {},
"source": [
"## Solution to[ Exercise 37.2](https://intro.quantecon.org/#lp_intro_ex2)\n",
"\n",
"Let us assume the carpenter produces $ x $ units of $ A $ and $ y $ units of $ B $.\n",
"\n",
"So we can formulate the problem as:\n",
"\n",
"$$\n",
"\\begin{aligned}\n",
"\\max_{x,y} \\ & z = 23 x + 10 y \\\\\n",
"\\mbox{subject to } \\ & x + y \\le 20 \\\\\n",
"& 2 x + 0.8 y \\le 25 \\\\\n",
"\\end{aligned}\n",
"$$"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "920fab64",
"metadata": {
"hide-output": false
},
"outputs": [],
"source": [
"# Instantiate a GLOP(Google Linear Optimization Package) solver\n",
"solver = pywraplp.Solver.CreateSolver('GLOP')"
]
},
{
"cell_type": "markdown",
"id": "f6738c3e",
"metadata": {},
"source": [
"Let’s create two variables $ x_1 $ and $ x_2 $ such that they can only have nonnegative values."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "4867c832",
"metadata": {
"hide-output": false
},
"outputs": [],
"source": [
"# Create the two variables and let them take on any non-negative value.\n",
"x = solver.NumVar(0, solver.infinity(), 'x')\n",
"y = solver.NumVar(0, solver.infinity(), 'y')"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "8c684c27",
"metadata": {
"hide-output": false
},
"outputs": [],
"source": [
"# Constraint 1: x + y <= 20.0\n",
"solver.Add(x + y <= 20.0)\n",
"\n",
"# Constraint 2: 2x + 0.8y <= 25.0\n",
"solver.Add(2 * x + 0.8 * y <= 25.0)"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "b6a14b90",
"metadata": {
"hide-output": false
},
"outputs": [],
"source": [
"# Objective function: 23x + 10y\n",
"solver.Maximize(23 * x + 10 * y)"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "ff060f7d",
"metadata": {
"hide-output": false
},
"outputs": [],
"source": [
"# Solve the system.\n",
"status = solver.Solve()\n",
"\n",
"if status == pywraplp.Solver.OPTIMAL:\n",
" print('Maximum Profit =', solver.Objective().Value())\n",
" x_sol = round(x.solution_value(), 3)\n",
" y_sol = round(y.solution_value(), 3)\n",
" print(f'(x, y): ({x_sol}, {y_sol})')\n",
"else:\n",
" print('The problem does not have an optimal solution.')"
]
}
],
"metadata": {
"date": 1750037600.741068,
"filename": "lp_intro.md",
"kernelspec": {
"display_name": "Python",
"language": "python3",
"name": "python3"
},
"title": "Linear Programming"
},
"nbformat": 4,
"nbformat_minor": 5
}