1.2 LP Problem - DMC - 3x3 and 4x4 (Generalized)-Copy1.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "460e40f9",
   "metadata": {},
   "source": [
    "# Linear Programming: DMC Formulation (3x3)\n",
    "- Author: Siang Lim, Shams Elnawawi\n",
    "- Last Updated: June 7th 2022\n",
    "- Created March 2nd 2022\n",
    "\n",
    "## References\n",
    "- Morshedi, A. M., Cutler, C. R., & Skrovanek, T. A. (1985). Optimal solution of dynamic matrix control with linear programing techniques (LDMC). In 1985 American Control Conference (pp. 199-208). IEEE.\n",
    "- Sorensen, R. C., & Cutler, C. R. (1998). LP integrates economics into dynamic matrix control. Hydrocarbon Processing, 77(9), 57-65.\n",
    "- Ranade, S. M., & Torres, E. (2009). From dynamic mysterious control to dynamic manageable control. Hydrocarbon Processing, 88(3), 77-81.\n",
    "- Godoy, J. L., Ferramosca, A., & González, A. H. (2017). Economic performance assessment and monitoring in LP-DMC type controller applications. Journal of Process Control, 57, 26-37."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c4f71c0d",
   "metadata": {},
   "source": [
    "## Before starting\n",
    "\n",
    "If this is your first time using this notebook, you may need to install a couple of packages to be able to run the code. If you end up with a ```ModuleNotFoundError```, go to the below cell (with the ```pip install``` lines) and uncomment both lines, then run the cell. \n",
    "\n",
    "You will have to restart the kernel after running installation commands, under the \"Kernel\" tab in the toolbar."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "6ecd4b52",
   "metadata": {},
   "outputs": [],
   "source": [
    "# pip install pulp\n",
    "# pip install matplotlib-label-lines"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "7e804c41",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Just importing libraries and tweaking the plot settings\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib widget\n",
    "\n",
    "from matplotlib.ticker import StrMethodFormatter\n",
    "import matplotlib.gridspec as gridspec\n",
    "from labellines import labelLine, labelLines\n",
    "\n",
    "from ipywidgets.widgets.interaction import interact\n",
    "import ipywidgets.widgets as widgets\n",
    "from ipywidgets import Layout\n",
    "\n",
    "# Import PuLP modeler functions\n",
    "from pulp import *\n",
    "\n",
    "fsize = 8\n",
    "tsize = 12\n",
    "tdir = 'in'\n",
    "major = 5.0\n",
    "minor = 3.0\n",
    "lwidth = 0.8\n",
    "lhandle = 2.0\n",
    "plt.style.use('default')\n",
    "plt.rcParams['font.size'] = fsize\n",
    "plt.rcParams['legend.fontsize'] = tsize\n",
    "plt.rcParams['xtick.direction'] = tdir\n",
    "plt.rcParams['ytick.direction'] = tdir\n",
    "plt.rcParams['xtick.major.size'] = major\n",
    "plt.rcParams['xtick.minor.size'] = minor\n",
    "plt.rcParams['ytick.major.size'] = 5.0\n",
    "plt.rcParams['ytick.minor.size'] = 3.0\n",
    "plt.rcParams['axes.linewidth'] = lwidth\n",
    "plt.rcParams['legend.handlelength'] = lhandle"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "368d8fa3",
   "metadata": {},
   "source": [
    "## MV-CV equations\n",
    "\n",
    "$$\n",
    "G =\n",
    "    \\begin{bmatrix}\n",
    "        -0.200 & -0.072 & 0.0774 \\\\\n",
    "        0.125 & -0.954 &  0.0063 \\\\\n",
    "        0.025 & 0.101 & −0.0143\n",
    "    \\end{bmatrix}\n",
    "$$\n",
    "\n",
    "Using the gain matrix, the CV relationship can be written in terms of its MVs, starting with:\n",
    "\n",
    "$$\n",
    "\\Delta \\text{CV}_{1} = G_{11} \\Delta \\text{MV}_{1} + G_{12} \\Delta \\text{MV}_{2} + G_{13} \\Delta \\text{MV}_{3} \\\\ \n",
    "\\Delta \\text{CV}_{2} = G_{21} \\Delta \\text{MV}_{1} + G_{22} \\Delta \\text{MV}_{2} + G_{23} \\Delta \\text{MV}_{3} \\\\\n",
    "\\Delta \\text{CV}_{3} = G_{31} \\Delta \\text{MV}_{1} + G_{32} \\Delta \\text{MV}_{2} + G_{33} \\Delta \\text{MV}_{3}\n",
    "$$\n",
    "\n",
    "We can impose upper and lower limits on the MVs:\n",
    "\n",
    "$$\n",
    "\\text{MV}_{1, \\text{Lo}} \\leq \\text{MV}_{1} \\leq \\text{MV}_{1, \\text{Hi}}\\\\ \n",
    "\\text{MV}_{2, \\text{Lo}} \\leq \\text{MV}_{2} \\leq \\text{MV}_{2, \\text{Hi}}\\\\\n",
    "\\text{MV}_{3, \\text{Lo}} \\leq \\text{MV}_{3} \\leq \\text{MV}_{3, \\text{Hi}}\\\\\n",
    "$$\n",
    "\n",
    "As well as the CVs:\n",
    "\n",
    "$$\n",
    "\\text{CV}_{1, \\text{Lo}} \\leq \\text{CV}_{1} \\leq \\text{CV}_{1, \\text{Hi}}\\\\ \n",
    "\\text{CV}_{2, \\text{Lo}} \\leq \\text{CV}_{2} \\leq \\text{CV}_{2, \\text{Hi}}\\\\\n",
    "\\text{CV}_{3, \\text{Lo}} \\leq \\text{CV}_{3} \\leq \\text{CV}_{3, \\text{Hi}}\\\\\n",
    "$$\n",
    "\n",
    "Since the CVs are related to the MVs by the gain matrix, we can substitute the equations to get CV limits in terms of MV movements:\n",
    "\n",
    "$$\n",
    "G_{11} \\Delta \\text{MV}_{1} + G_{12} \\Delta \\text{MV}_{2} + G_{13} \\Delta \\text{MV}_{3} \\leq \\Delta \\text{CV}_{1, \\text{Hi}}\\\\\n",
    "G_{11} \\Delta \\text{MV}_{1} + G_{12} \\Delta \\text{MV}_{2} + G_{13} \\Delta \\text{MV}_{3} \\geq \\Delta \\text{CV}_{1, \\text{Lo}}\\\\\n",
    "G_{21} \\Delta \\text{MV}_{1} + G_{22} \\Delta \\text{MV}_{2} + G_{23} \\Delta \\text{MV}_{3} \\leq \\Delta \\text{CV}_{2, \\text{Hi}}\\\\\n",
    "G_{21} \\Delta \\text{MV}_{1} + G_{22} \\Delta \\text{MV}_{2} + G_{23} \\Delta \\text{MV}_{3} \\geq \\Delta \\text{CV}_{2, \\text{Lo}}\\\\\n",
    "G_{31} \\Delta \\text{MV}_{1} + G_{32} \\Delta \\text{MV}_{2} + G_{33} \\Delta \\text{MV}_{3} \\leq \\Delta \\text{CV}_{3, \\text{Hi}}\\\\\n",
    "G_{31} \\Delta \\text{MV}_{1} + G_{22} \\Delta \\text{MV}_{2} + G_{33} \\Delta \\text{MV}_{3} \\geq \\Delta \\text{CV}_{3, \\text{Lo}}\\\\\n",
    "$$\n",
    "\n",
    "Let's see what this looks like:"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "90f372d2",
   "metadata": {},
   "source": [
    "# Gains and CV Limits"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "6873757d",
   "metadata": {},
   "outputs": [],
   "source": [
    "G11 = -0.200\n",
    "G12 = -0.072\n",
    "G13 = 0.0774\n",
    "G14 = 0.0574\n",
    "\n",
    "G21 = 0.125\n",
    "G22 = -0.954\n",
    "G23 = 0.0063\n",
    "G24 = 0.0374\n",
    "\n",
    "G31 = 0.025\n",
    "G32 = 0.101\n",
    "G33 = -0.0143\n",
    "G34 = -0.1143\n",
    "\n",
    "G41 = 0.120\n",
    "G42 = 0.150\n",
    "G43 = -0.1143\n",
    "G44 = -0.0943\n",
    "\n",
    "CV1Lo = -6\n",
    "CV1Hi = 6\n",
    "\n",
    "CV2Lo = -10\n",
    "CV2Hi = 10.5\n",
    "\n",
    "CV3Lo = -3\n",
    "CV3Hi = 3.5\n",
    "\n",
    "CV4Lo = -4.5\n",
    "CV4Hi = 4\n",
    "\n",
    "cost_MV1 = -1\n",
    "cost_MV2 = -1\n",
    "cost_MV3 = -1\n",
    "cost_MV4 = 1\n",
    "\n",
    "limits = 10\n",
    "plot_limits = 20\n",
    "\n",
    "MV1Lo = -limits\n",
    "MV1Hi = limits\n",
    "MV2Lo = -limits\n",
    "MV2Hi = limits\n",
    "MV3Lo = -limits\n",
    "MV3Hi = 15\n",
    "MV4Lo = -limits\n",
    "MV4Hi = -limits"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "562bb6e4",
   "metadata": {},
   "outputs": [],
   "source": [
    "G           = [(G11, G12, G13, G14), \n",
    "               (G21, G22, G23, G24), \n",
    "               (G31, G32, G33, G34), \n",
    "               (G41, G42, G43, G44)]\n",
    "\n",
    "CV_values   = [(CV1Lo, CV1Hi), \n",
    "               (CV2Lo, CV2Hi), \n",
    "               (CV3Lo, CV3Hi)] \n",
    "#                (CV4Lo, CV4Hi)]\n",
    "CV_init_vals= [(-2.5,2.0), \n",
    "               (-4,4.5), \n",
    "               (-0.9,0.2)] \n",
    "#                (-4,4.5)]\n",
    "MV_costs    = [cost_MV1, \n",
    "               cost_MV2, \n",
    "               cost_MV3]\n",
    "#                cost_MV4]\n",
    "MV_values   = [(MV1Lo, MV1Hi), \n",
    "               (MV2Lo, MV2Hi), \n",
    "               (MV3Lo, MV3Hi)] \n",
    "#                (MV4Lo, MV4Hi)]\n",
    "\n",
    "# (x,y,c) triplets of MVs, c for constant, i.e. \n",
    "# (MV1, MV2, MV3) for the first plot, \n",
    "# (MV1, MV3, MV2) for the second, \n",
    "# (MV2, MV3, MV1) for the third\n",
    "# plot_MV_indices = [(0,1), (0,2), (1,2)]\n",
    "\n",
    "nCVs = len(CV_values)\n",
    "nMVs = len(MV_values)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "a7fdf889-850d-40db-93d0-a93bd835d92e",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "{'CV1 Limits': (-2.5, 2.0), 'CV2 Limits': (-4.0, 4.5), 'CV3 Limits': (-0.9, 0.2), 'MV1 Cost': -1.0, 'MV2 Cost': -1.0, 'MV3 Cost': -1.0, 'MV1 Limits': (-10.0, 10.0), 'MV2 Limits': (-10.0, 10.0), 'MV3 Limits': (-10.0, 15.0)}\n"
     ]
    }
   ],
   "source": [
    "CV_widgets   = []\n",
    "MV_widgets   = []\n",
    "cost_widgets = []\n",
    "stepsize = 0.2\n",
    "\n",
    "# Make CV sliders\n",
    "for i in range(nCVs):\n",
    "    widget = widgets.FloatRangeSlider(\n",
    "                value=CV_init_vals[i], min=CV_values[i][0], max=CV_values[i][1], step=stepsize,\n",
    "                description=f'CV{i+1} Limits', continuous_update=False)\n",
    "    CV_widgets.append(widget)\n",
    "\n",
    "# Make MV sliders\n",
    "for i in range(nMVs):\n",
    "    widget = widgets.FloatRangeSlider(\n",
    "                value=MV_values[i], min=MV_values[i][0], max=MV_values[i][1], step=stepsize,\n",
    "                description=f'MV{i+1} Limits', continuous_update=False)\n",
    "    MV_widgets.append(widget)\n",
    "    \n",
    "# Make cost sliders\n",
    "for i in range(nMVs):\n",
    "    widget = widgets.FloatSlider(\n",
    "                value=MV_costs[i], min=-3, max=3, step=stepsize/2, # finer step for costs\n",
    "                description=f'MV{i+1} Cost', continuous_update=True)\n",
    "    cost_widgets.append(widget)\n",
    "    \n",
    "sliders = CV_widgets + cost_widgets + MV_widgets\n",
    "values = [slider.value for slider in sliders]\n",
    "values_dict = {}\n",
    "for slider in sliders:\n",
    "    values_dict[slider.description] = slider.value\n",
    "print(values_dict)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "cb507e02-fab8-41d7-a338-89f60cbe41f2",
   "metadata": {},
   "outputs": [],
   "source": [
    "# display ui\n",
    "ncols = 3\n",
    "\n",
    "nrows = nCVs // ncols + 1\n",
    "CV_row = []\n",
    "for i in range(nrows):\n",
    "    CV_row.append(widgets.VBox(CV_widgets[(ncols)*i:(ncols)*i+ncols]))\n",
    "    \n",
    "nrows = nMVs // ncols + 1\n",
    "MV_row = []\n",
    "cost_row = []\n",
    "for i in range(nrows):\n",
    "    MV_row.append(widgets.VBox(MV_widgets[(ncols)*i:(ncols)*i+ncols]))\n",
    "    cost_row.append(widgets.VBox(cost_widgets[(ncols)*i:(ncols)*i+ncols]))\n",
    "\n",
    "ui = widgets.VBox([widgets.VBox(CV_row),\n",
    "              widgets.VBox(MV_row),\n",
    "              widgets.VBox(cost_row)])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "337d52ee",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solve the LP\n",
    "def run_lp(values_dict):\n",
    "    prob = LpProblem(\"DMC_problem\",LpMinimize)\n",
    "    \n",
    "    # How many MVs\n",
    "    MVs = []\n",
    "    for i in range(nMVs):\n",
    "        MVs.append(LpVariable(f\"MV{i+1}\",-limits))\n",
    "    \n",
    "    # the objective function\n",
    "    obj = 0\n",
    "    for indx, MV in enumerate(MVs):\n",
    "        obj += values_dict[f\"MV{indx+1} Cost\"]*MV\n",
    "        \n",
    "    prob += obj, \"Cost function of MVs\"\n",
    "    \n",
    "    # constraint formulation in terms of MV1 and MV2\n",
    "    CV_contraint_lo = []\n",
    "    CV_contraint_hi = []\n",
    "    \n",
    "    for i in range(nCVs):\n",
    "        c = 0\n",
    "        for indx, MV in enumerate(MVs):\n",
    "            c += G[i][indx]*MV\n",
    "        prob += c <= values_dict[f'CV{i+1} Limits'][1], f'CV{i+1} High Limit'\n",
    "        prob += c >= values_dict[f'CV{i+1} Limits'][0], f'CV{i+1} Low Limit'\n",
    "    \n",
    "#     prob += G[i][0]*MV1 + G[i][1]*MV2 <= values_dict[f'CV{i+1} Limits'][1], f'CV{i+1} High Limit'\n",
    "#     prob += G[i][0]*MV1 + G[i][1]*MV2 >= values_dict[f'CV{i+1} Limits'][0], f'CV{i+1} Low Limit'\n",
    "    \n",
    "    for indx, MV in enumerate(MVs):\n",
    "        prob += MV <= values_dict[f'MV{indx+1} Limits'][1], f'MV{indx+1} High Limit'\n",
    "        prob += MV >= values_dict[f'MV{indx+1} Limits'][0], f'MV{indx+1} Low Limit'        \n",
    "    \n",
    "    if (prob.solve(PULP_CBC_CMD(msg=0)) == 1):\n",
    "#         print([v.varValue for v in prob.variables()])\n",
    "#         print([v for v in prob.variables()])\n",
    "        return [v.varValue for v in prob.variables()], value(prob.objective)\n",
    "    else:\n",
    "        print(\"NOT SOLVED - Infeasibility!\")\n",
    "        return np.zeros(nMVs), 0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "f0e50a42",
   "metadata": {},
   "outputs": [],
   "source": [
    "d = np.linspace(-plot_limits, plot_limits, 1000)\n",
    "x,y = np.meshgrid(d,d)\n",
    "\n",
    "# Recalculate shaded regions\n",
    "constraints = []\n",
    "for i in range(nCVs):\n",
    "    c_hi = G[i][0]*x + G[i][1]*y <= values_dict[f'CV{i+1} Limits'][1]\n",
    "    c_lo = G[i][0]*x + G[i][1]*y >= values_dict[f'CV{i+1} Limits'][0]\n",
    "    constraints.append(c_lo)\n",
    "    constraints.append(c_hi)\n",
    "\n",
    "# the MVs\n",
    "for indx, var in enumerate([x,y]):\n",
    "    c_hi = var <= values_dict[f'MV{indx+1} Limits'][1]\n",
    "    c_lo = var >= values_dict[f'MV{indx+1} Limits'][0]\n",
    "    constraints.append(c_lo)\n",
    "    constraints.append(c_hi)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "3b8daec8",
   "metadata": {},
   "outputs": [],
   "source": [
    "def handle_slider_change(change):\n",
    "    ## grab slider vals\n",
    "    values_dict = {}\n",
    "    for slider in sliders:\n",
    "        values_dict[slider.description] = slider.value\n",
    "    \n",
    "    # Find the LP soln\n",
    "    soln, V = run_lp(values_dict)\n",
    "    \n",
    "    for indx, key in enumerate(axs_dict):\n",
    "        ax = axs_dict[key]\n",
    "        \n",
    "        # remove previous line labels (TODO could be more efficient, modify labelLine library to return handles)\n",
    "        # for txt in ax.texts:\n",
    "        #     txt.remove()\n",
    "        \n",
    "        # figure out which MVs we are plotting, and which ones are constant\n",
    "        x_MV = key[1] # column is x-axis\n",
    "        y_MV = key[0]+1 # row is y-axis\n",
    "        c_MV = [a for a in range(nMVs) if a != x_MV and a != y_MV]        \n",
    "\n",
    "        # Recalculate shaded regions \n",
    "        constraints = []\n",
    "        \n",
    "        # the obj func\n",
    "        y_obj = (V - sum([values_dict[f'MV{cmv+1} Cost']*soln[cmv] for cmv in c_MV]) - values_dict[f'MV{x_MV+1} Cost']*d)/values_dict[f'MV{y_MV+1} Cost']\n",
    "                \n",
    "        # Update CV constraint lines\n",
    "        for i in range(nCVs):\n",
    "            cv_lim_lo = values_dict[f'CV{i+1} Limits'][0]\n",
    "            cv_lim_hi = values_dict[f'CV{i+1} Limits'][1]\n",
    "\n",
    "            # for CV constraint lines\n",
    "            y_lo = (cv_lim_lo - G[i][x_MV]*d - sum([G[i][cmv]*soln[cmv] for cmv in c_MV]))/G[i][y_MV]\n",
    "            y_hi = (cv_lim_hi - G[i][x_MV]*d - sum([G[i][cmv]*soln[cmv] for cmv in c_MV]))/G[i][y_MV]\n",
    "            lines_handler_dict[key]['CV_lines_lo'][i].set_data(d, y_lo)\n",
    "            lines_handler_dict[key]['CV_lines_hi'][i].set_data(d, y_hi)\n",
    "\n",
    "            # for shading\n",
    "            c_hi = G[i][x_MV]*x + G[i][y_MV]*y + sum([G[i][cmv]*soln[cmv] for cmv in c_MV]) <= cv_lim_hi\n",
    "            c_lo = G[i][x_MV]*x + G[i][y_MV]*y + sum([G[i][cmv]*soln[cmv] for cmv in c_MV]) >= cv_lim_lo\n",
    "            constraints.append(c_lo)\n",
    "            constraints.append(c_hi)\n",
    "            \n",
    "            # search for valid xvals where the yvals are within plots limits\n",
    "            # x_valid_lo = d[(y_lo < plot_limits) & (y_lo > -plot_limits)]\n",
    "            # x_valid_hi = d[(y_hi < plot_limits) & (y_hi > -plot_limits)]\n",
    "            # offset_lo = max(0,min(10, len(x_valid_lo)-3))\n",
    "            # offset_hi = max(0,min(10, len(x_valid_hi)-3))\n",
    "            # if len(x_valid_lo) > 0:\n",
    "            #     txt_cv_lo = labelLine(line_lo, x_valid_lo[offset_lo], fontsize=5, zorder=2.5)\n",
    "            # if len(x_valid_hi) > 0:\n",
    "            #     txt_cv_hi = labelLine(line_hi, x_valid_hi[-offset_hi], fontsize=5, zorder=2.5)            \n",
    "\n",
    "        # Update MV constraint lines\n",
    "        lines_handler_dict[key]['v_lo'].set_data([values_dict[f'MV{x_MV+1} Limits'][0],values_dict[f'MV{x_MV+1} Limits'][0]], [-limits, limits])\n",
    "        lines_handler_dict[key]['v_hi'].set_data([values_dict[f'MV{x_MV+1} Limits'][1],values_dict[f'MV{x_MV+1} Limits'][1]], [-limits, limits])\n",
    "        lines_handler_dict[key]['h_lo'].set_data([-limits, limits], [values_dict[f'MV{y_MV+1} Limits'][0],values_dict[f'MV{y_MV+1} Limits'][0]])\n",
    "        lines_handler_dict[key]['h_hi'].set_data([-limits, limits], [values_dict[f'MV{y_MV+1} Limits'][1],values_dict[f'MV{y_MV+1} Limits'][1]])    \n",
    "\n",
    "        # the 4 MV limits for this ax\n",
    "        constraints.append(x >= values_dict[f'MV{x_MV+1} Limits'][0])\n",
    "        constraints.append(x <= values_dict[f'MV{x_MV+1} Limits'][1])\n",
    "        constraints.append(y >= values_dict[f'MV{y_MV+1} Limits'][0])\n",
    "        constraints.append(y <= values_dict[f'MV{y_MV+1} Limits'][1]) \n",
    "\n",
    "        # Shade the right regions\n",
    "        lines_handler_dict[key]['im'].set_data((np.logical_and.reduce(constraints)).astype(float))\n",
    "        \n",
    "        # the quiver/vector field\n",
    "        z_obj = (values_dict[f'MV{x_MV+1} Cost'] * xv) + (values_dict[f'MV{y_MV+1} Cost'] * yv) + sum([values_dict[f'MV{cmv+1} Cost']*soln[cmv] for cmv in c_MV])\n",
    "        lines_handler_dict[key]['quiver'].set_UVC(-values_dict[f'MV{x_MV+1} Cost'],-values_dict[f'MV{y_MV+1} Cost'], z_obj)\n",
    "\n",
    "        # the soln\n",
    "        lines_handler_dict[key]['soln_marker'].set_data(soln[x_MV], soln[y_MV]);\n",
    "        lines_handler_dict[key]['soln_text'].set_position((soln[x_MV], soln[y_MV]));\n",
    "        lines_handler_dict[key]['soln_text'].set_text(\"({:.1f}, {:.1f})\".format(soln[x_MV], soln[y_MV]))\n",
    "        lines_handler_dict[key]['soln_func'].set_data(d, y_obj);\n",
    "    \n",
    "    fig.canvas.draw()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "8236d3a0-005d-440d-9482-15ec484d788a",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "9fb33089b9d34fc9b14df25515cf3d64",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(Output(), VBox(children=(HTMLMath(value='<p>Use the controls to interact with this linear progr…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Initialize plots\n",
    "d = np.linspace(-plot_limits, plot_limits, 100)\n",
    "x,y = np.meshgrid(d,d)\n",
    "\n",
    "gs = gridspec.GridSpec(nMVs-1, nMVs-1)\n",
    "gs.update(wspace=0.05, hspace=0.05)\n",
    "gs_indices = [(row, col) for row in range(nMVs-2,-1,-1) for col in range(row+1)]\n",
    "axs_dict = {}\n",
    "\n",
    "# mesh for vector field\n",
    "dvec = np.linspace(-plot_limits + (0.1*plot_limits), plot_limits - (0.1*plot_limits), 12)\n",
    "xv, yv = np.meshgrid(dvec, dvec)\n",
    "\n",
    "# init soln\n",
    "soln, V = run_lp(values_dict)\n",
    "\n",
    "# initialize line handlers\n",
    "CV_lines_lo = []\n",
    "CV_lines_hi = []\n",
    "MV_lines_lo = []\n",
    "MV_lines_hi = []\n",
    "lines_handler_dict = {}\n",
    "\n",
    "colors = ['r', 'b', 'y', 'g'] # TODO generalize this using a cmap instead of defining manual colors??!\n",
    "\n",
    "# plot as widget\n",
    "output = widgets.Output()\n",
    "with output:\n",
    "    fig = plt.figure(figsize=(6,6), dpi=100, facecolor='white')\n",
    "    plt.show()\n",
    "    \n",
    "for r,c in gs_indices:\n",
    "    # build the shared axes correctly and handle the tick labels\n",
    "    if (r == nMVs-2 and c == 0): # no shared axis for the ax, bottom left corner\n",
    "        ax = plt.subplot(gs[r,c])\n",
    "    elif(c != 0): # all non-first columns share a y-axis with stuff to the left of it\n",
    "        ax = plt.subplot(gs[r,c], sharey=axs_dict[(r,0)])\n",
    "        plt.setp(ax.get_yticklabels(), visible=False)\n",
    "    elif(r != nMVs-2): # all non-first rows share a x-axis with stuff below it\n",
    "        ax = plt.subplot(gs[r,c], sharex=axs_dict[(nMVs-2,c)])\n",
    "        plt.setp(ax.get_xticklabels(), visible=False)\n",
    "\n",
    "    # add the labels\n",
    "    if(r == nMVs-2):\n",
    "        ax.set_xlabel(f'$\\Delta MV_{c+1}$')\n",
    "    if(c == 0):\n",
    "        ax.set_ylabel(f'$\\Delta MV_{r+2}$')\n",
    "\n",
    "    axs_dict[(r,c)] = ax\n",
    "\n",
    "for indx, key in enumerate(axs_dict):\n",
    "    ax = axs_dict[key]\n",
    "    lines_handler_dict[key] = {}\n",
    "\n",
    "    # figure out which MVs we are plotting, and which ones are constant\n",
    "    x_MV = key[1] # column is x-axis\n",
    "    y_MV = key[0]+1 # row is y-axis\n",
    "    c_MV = [a for a in range(nMVs) if a != x_MV and a != y_MV]\n",
    "\n",
    "    # plot the MV limits\n",
    "    v_lo = ax.axvline(x=values_dict[f'MV{x_MV+1} Limits'][0], color='gray', lw=1, label=f'MV{x_MV+1} Lo')\n",
    "    v_hi = ax.axvline(x=values_dict[f'MV{x_MV+1} Limits'][1], color='gray', lw=1, label=f'MV{x_MV+1} Hi')\n",
    "    h_lo = ax.axhline(y=values_dict[f'MV{y_MV+1} Limits'][0], color='gray', lw=1, label=f'MV{y_MV+1} Lo')\n",
    "    h_hi = ax.axhline(y=values_dict[f'MV{y_MV+1} Limits'][1], color='gray', lw=1, label=f'MV{y_MV+1} Hi')\n",
    "\n",
    "    # labelLines([v_lo, v_hi], fontsize=4)\n",
    "    # labelLine(h_hi, fontsize=4)\n",
    "        \n",
    "    # store the line handlers per ax in a dict\n",
    "    lines_handler_dict[key]['v_lo'] = v_lo\n",
    "    lines_handler_dict[key]['v_hi'] = v_hi\n",
    "    lines_handler_dict[key]['h_lo'] = h_lo\n",
    "    lines_handler_dict[key]['h_hi'] = h_hi\n",
    "\n",
    "    # Recalculate shaded regions \n",
    "    constraints = []\n",
    "\n",
    "    # calculate CV lines\n",
    "    lines_handler_dict[key]['CV_lines_lo'] = []\n",
    "    lines_handler_dict[key]['CV_lines_hi'] = []    \n",
    "    for i in range(nCVs):\n",
    "        cv_lim_lo = values_dict[f'CV{i+1} Limits'][0]\n",
    "        cv_lim_hi = values_dict[f'CV{i+1} Limits'][1]\n",
    "\n",
    "        # for CV constraint lines\n",
    "        y_lo = (cv_lim_lo - G[i][x_MV]*d - sum([G[i][cmv]*soln[cmv] for cmv in c_MV]))/G[i][y_MV]\n",
    "        y_hi = (cv_lim_hi - G[i][x_MV]*d - sum([G[i][cmv]*soln[cmv] for cmv in c_MV]))/G[i][y_MV]\n",
    "        line_lo, = ax.plot(d, y_lo, f'--{colors[i]}', label=f'$CV_{i+1}$ Lo');\n",
    "        line_hi, = ax.plot(d, y_hi, f'-{colors[i]}', label=f'$CV_{i+1}$ Hi');\n",
    "        lines_handler_dict[key]['CV_lines_lo'].append(line_lo)\n",
    "        lines_handler_dict[key]['CV_lines_hi'].append(line_hi)\n",
    "        \n",
    "        # # search for valid xvals where the yvals are within plots limits\n",
    "        # x_valid_lo = d[(y_lo < plot_limits) & (y_lo > -plot_limits)]\n",
    "        # x_valid_hi = d[(y_hi < plot_limits) & (y_hi > -plot_limits)]\n",
    "        # offset_lo = max(0,min(10, len(x_valid_lo)-3))\n",
    "        # offset_hi = max(0,min(10, len(x_valid_hi)-3))\n",
    "        # if len(x_valid_lo) > 0:\n",
    "        #     labelLine(line_lo, x_valid_lo[offset_lo], fontsize=5, zorder=2.5)\n",
    "        # if len(x_valid_hi) > 0:\n",
    "        #     labelLine(line_hi, x_valid_hi[-offset_hi], fontsize=5, zorder=2.5)\n",
    "        \n",
    "        # for shading\n",
    "        c_hi = G[i][x_MV]*x + G[i][y_MV]*y + sum([G[i][cmv]*soln[cmv] for cmv in c_MV]) <= cv_lim_hi\n",
    "        c_lo = G[i][x_MV]*x + G[i][y_MV]*y + sum([G[i][cmv]*soln[cmv] for cmv in c_MV]) >= cv_lim_lo\n",
    "        constraints.append(c_lo)\n",
    "        constraints.append(c_hi)\n",
    "            \n",
    "\n",
    "    # the 4 MV limits for this ax\n",
    "    constraints.append(x >= values_dict[f'MV{x_MV+1} Limits'][0])\n",
    "    constraints.append(x <= values_dict[f'MV{x_MV+1} Limits'][1])\n",
    "    constraints.append(y >= values_dict[f'MV{y_MV+1} Limits'][0])\n",
    "    constraints.append(y <= values_dict[f'MV{y_MV+1} Limits'][1])    \n",
    "\n",
    "    # the obj function        \n",
    "    y_obj = (V - sum([values_dict[f'MV{cmv+1} Cost']*soln[cmv] for cmv in c_MV]) - values_dict[f'MV{x_MV+1} Cost']*d)/values_dict[f'MV{y_MV+1} Cost']\n",
    "    soln_func, = ax.plot(d, y_obj, '--k')\n",
    "    soln_marker, = ax.plot(soln[x_MV], soln[y_MV], 'ok', zorder=10)\n",
    "    soln_text = ax.text(soln[x_MV], soln[y_MV], '({:.2f},{:.2f})'.format(soln[x_MV], soln[y_MV]))\n",
    "\n",
    "    lines_handler_dict[key]['soln_func'] = soln_func\n",
    "    lines_handler_dict[key]['soln_marker'] = soln_marker\n",
    "    lines_handler_dict[key]['soln_text'] = soln_text\n",
    "    \n",
    "    z_obj = (values_dict[f'MV{x_MV+1} Cost'] * xv) + (values_dict[f'MV{y_MV+1} Cost'] * yv) + sum([values_dict[f'MV{cmv+1} Cost']*soln[cmv] for cmv in c_MV])\n",
    "    lines_handler_dict[key]['quiver'] = ax.quiver(xv, yv,-values_dict[f'MV{x_MV+1} Cost'],-values_dict[f'MV{y_MV+1} Cost'], z_obj, cmap='gray', headwidth=4, width=0.003, scale=40, alpha=0.5)\n",
    "\n",
    "    ax.plot(0,0,'kx');\n",
    "    ax.set_aspect('equal')\n",
    "    im = ax.imshow(np.logical_and.reduce(constraints).astype(int), extent=(x.min(),x.max(),y.min(),y.max()), origin=\"lower\", cmap=\"Blues\", alpha=0.1)\n",
    "    lines_handler_dict[key]['im'] = im\n",
    "        \n",
    "####################################################\n",
    "# END OF INIT FUNC\n",
    "####################################################\n",
    "\n",
    "# register slides\n",
    "for widget in sliders:\n",
    "    widget.observe(handle_slider_change, names='value')\n",
    "\n",
    "widgets.HBox([output, \n",
    "              widgets.VBox([widgets.HTMLMath(\n",
    "                            value=r\"<p>Use the controls to interact with this linear program</p><br><br>\"), \n",
    "                            ui])])\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c752d9b1-788f-4c3d-b3ae-4f98877f6d83",
   "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.8.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}