L10-LP-NLP.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "distinguished-badge",
   "metadata": {},
   "source": [
    "# Lecture 10\n",
    "_____________\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fancy-qualification",
   "metadata": {},
   "source": [
    "Import libraries"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "grateful-deficit",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np \n",
    "import matplotlib.pyplot as plt\n",
    "plt.rcParams['figure.dpi'] = 600\n",
    "from scipy import optimize"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "enabling-familiar",
   "metadata": {},
   "source": [
    "## _NLP examples_\n",
    "\n",
    "Non-linear programming examples for a simple constrained problem"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "level-agency",
   "metadata": {},
   "source": [
    "2-user water allocation example\n",
    "\n",
    "_First approach:_ `scipy.optimize.minimize`"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "foster-peeing",
   "metadata": {},
   "source": [
    "__Objective function__ (x is decision variable vector)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "dried-proxy",
   "metadata": {},
   "outputs": [],
   "source": [
    "def f(x):\n",
    "    return -1*(30*x[0] - 4*x[0]**2 + 10*x[1] - 2*x[1]**2) # quadratic benefit function (maximize)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "western-integrity",
   "metadata": {},
   "source": [
    "If constraint function returns > zero, constraint satisfied. otherwise not satisfied"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "partial-celebrity",
   "metadata": {},
   "outputs": [],
   "source": [
    "def f_constraint(x):\n",
    "    return -(x[0] + x[1] - 5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "textile-ground",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "     fun: -66.66666666666687\n",
      "     jac: array([-3.33333302, -3.33333397])\n",
      " message: 'Optimization terminated successfully'\n",
      "    nfev: 7\n",
      "     nit: 2\n",
      "    njev: 2\n",
      "  status: 0\n",
      " success: True\n",
      "       x: array([3.33333333, 1.66666667])\n"
     ]
    }
   ],
   "source": [
    "constraint = {'type': 'ineq', 'fun': f_constraint}\n",
    "\n",
    "sol = optimize.minimize(f, x0 = [0,5], constraints=[constraint])\n",
    "print(sol) # does not return duals"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "subjective-conditioning",
   "metadata": {},
   "source": [
    "Decision variables `sol.x`, objective function `sol.fun`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "intellectual-italic",
   "metadata": {},
   "outputs": [],
   "source": [
    "import cvxpy as cvx"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "id": "recovered-carol",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Objective = 66.666667\n",
      "X1 = 3.333333\n",
      "X2 = 1.666667\n",
      "Dual (x1 + x2 <= 5.0) = 3.333333\n",
      "Dual (0.0 <= x1) = 0.000000\n",
      "Dual (0.0 <= x2) = 0.000000\n"
     ]
    }
   ],
   "source": [
    "x1 = cvx.Variable(name='x1')\n",
    "x2 = cvx.Variable(name='x2')\n",
    "Q = 5 # units of water\n",
    "\n",
    "obj = cvx.Maximize(30*x1 - 4*x1**2 + 10*x2 - 2*x2**2)\n",
    "\n",
    "constraints = [x1 + x2 <= Q, x1 >= 0, x2 >= 0] # magic?\n",
    "\n",
    "prob = cvx.Problem(obj, constraints)\n",
    "prob.solve()\n",
    "\n",
    "print('Objective = %f' % obj.value)\n",
    "print('X1 = %f' % x1.value)\n",
    "print('X2 = %f' % x2.value)\n",
    "\n",
    "for c in constraints:\n",
    "  print('Dual (%s) = %f' % (c, c.dual_value))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "terminal-vessel",
   "metadata": {},
   "source": [
    "------------"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "prescribed-edwards",
   "metadata": {},
   "source": [
    "## _LP farm example_\n",
    "\n",
    "Linear programming examples for crop choice irrigation example problem"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "demanding-infection",
   "metadata": {},
   "source": [
    "_First approach:_ `scipy.optimize.linprog`\n",
    "\n",
    "Need matrix form: minimize c^T * x, subject to Ax <= b"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "id": "crazy-insertion",
   "metadata": {},
   "outputs": [],
   "source": [
    "c = [-5, -3] # negative to maximize (b/c the default is minimize)\n",
    "A = [[10,5], [1,1.5], [2,2], [-1,0], [0,-1]]\n",
    "b = [20, 3, 4.5, 0, 0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "id": "previous-rochester",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Scipy Output:\n",
      "     con: array([], dtype=float64)\n",
      "     fun: -10.249999999652971\n",
      " message: 'Optimization terminated successfully.'\n",
      "     nit: 4\n",
      "   slack: array([ 8.90686636e-10,  5.00000000e-01, -1.84900983e-11,  1.75000000e+00,\n",
      "        5.00000000e-01])\n",
      "  status: 0\n",
      " success: True\n",
      "       x: array([1.75, 0.5 ])\n"
     ]
    }
   ],
   "source": [
    "sol = optimize.linprog(c, A, b)\n",
    "print('Scipy Output:')\n",
    "print(sol)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "handmade-department",
   "metadata": {},
   "source": [
    "The x's are the decision variables"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fifteen-helen",
   "metadata": {},
   "source": [
    "_Second approach:_ `cvxpy`\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "id": "extensive-tobago",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "10.249999998189105"
      ]
     },
     "execution_count": 38,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "xc = cvx.Variable(name='xc')\n",
    "xb = cvx.Variable(name='xb')\n",
    "pc = 5\n",
    "pb = 3\n",
    "\n",
    "obj = cvx.Maximize(pc*xc + pb*xb)\n",
    "\n",
    "constraints = [10*xc + 5*xb <= 20, # cvxpy recognizes the logical comparison signs\n",
    "               xc + 1.5*xb <= 3,\n",
    "               2*xc + 2*xb <= 4.5,\n",
    "               xc >= 0,\n",
    "               xb >= 0]\n",
    "\n",
    "prob = cvx.Problem(obj, constraints)\n",
    "prob.solve()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "id": "prescription-reconstruction",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "cvxpy Output:\n",
      "Objective = 10.250000\n",
      "xc = 1.750000\n",
      "xb = 0.500000\n",
      "Dual (10.0 @ xc + 5.0 @ xb <= 20.0) = 0.400000\n",
      "Dual (xc + 1.5 @ xb <= 3.0) = 0.000000\n",
      "Dual (2.0 @ xc + 2.0 @ xb <= 4.5) = 0.500000\n",
      "Dual (0.0 <= xc) = 0.000000\n",
      "Dual (0.0 <= xb) = 0.000000\n"
     ]
    }
   ],
   "source": [
    "print('\\ncvxpy Output:')\n",
    "print('Objective = %f' % obj.value)\n",
    "print('xc = %f' % xc.value)\n",
    "print('xb = %f' % xb.value)\n",
    "\n",
    "for c in constraints:\n",
    "    print('Dual (%s) = %f' % (c, c.dual_value))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "contemporary-restriction",
   "metadata": {},
   "source": [
    "Should match the answer above"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "deluxe-gossip",
   "metadata": {},
   "source": [
    "The units of the dual values come from the objective function and constraint value ($ / constraint value)\n",
    "\n",
    "Constraints with nonzero dual values are binding"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "hungarian-techno",
   "metadata": {},
   "source": [
    "-----------"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d4c1d696-b2a9-4429-95ec-53e7278e2d50",
   "metadata": {},
   "source": [
    "#### Interactive widget"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "id": "919656c6-66af-459e-a5e3-dddc442bfb7a",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "def lpplot(pc=5, pb=3, W=20, L=3, F=4.5):\n",
    "    A = np.array([[10,1,2,-1,0], [5,1.5,2,0,-1 ]]).T\n",
    "    b = np.array([ W, L, F, 0, 0])\n",
    "    c = np.array([ -pc, -pb ])\n",
    "    sol = optimize.linprog(c,A,b)\n",
    "\n",
    "    x = np.arange(0,5.01,0.01)\n",
    "    X1,X2 = np.meshgrid(x,x)\n",
    "    Z = -c[0]*X1 + -c[1]*X2\n",
    "\n",
    "    # Contour lines in the background\n",
    "    plt.contour(X1,X2,Z,50,cmap=plt.cm.cool)\n",
    "    plt.colorbar()\n",
    "\n",
    "    # Plot constraint lines\n",
    "    plt.plot(x, (W-10*x)/5, color='k')\n",
    "    plt.plot(x, (L-x)/1.5, color='k')\n",
    "    plt.plot(x, (F-2*x)/2, color='k')\n",
    "\n",
    "    # And the optimal point ...\n",
    "    plt.scatter(sol['x'][0], sol['x'][1], s=40, color='r', zorder=5)\n",
    "    plt.xlim([0,3])\n",
    "    plt.ylim([0,3])\n",
    "    plt.xlabel('$x_c$')\n",
    "    plt.ylabel('$x_b$')\n",
    "    plt.show()\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "id": "0bb14cd0-9735-4920-a17e-e998d2671ea5",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "8abbf14006b74b07ae4f1ae2ba1a5ffc",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(FloatSlider(value=5.0, description='pc', max=12.0), FloatSlider(value=3.0, description='…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "<function __main__.lpplot(pc=5, pb=3, W=20, L=3, F=4.5)>"
      ]
     },
     "execution_count": 41,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from ipywidgets import interact\n",
    "\n",
    "interact(lpplot, pc=(0.,12.), pb=(0.,12.), W=(15.,25.), L=(1.,5.), F=(3.,6.))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c822fcfc-6459-4043-a9f3-4c2fe19cc1b3",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}