{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "LRtLKOk7Z1XH"
   },
   "source": [
    "# Introduction to Disjunctive Programming"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "-XCVZJk2__cT"
   },
   "source": [
    "## Installations and imports"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "id": "v_-eaAb3Z0np"
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib as mpl\n",
    "import pandas as pd\n",
    "\n",
    "import shutil\n",
    "import sys\n",
    "import os.path\n",
    "\n",
    "if not shutil.which(\"pyomo\"):\n",
    "    !pip install -q pyomo\n",
    "    assert(shutil.which(\"pyomo\"))\n",
    "\n",
    "if not (shutil.which(\"cbc\") or os.path.isfile(\"cbc\")):\n",
    "    if \"google.colab\" in sys.modules:\n",
    "        !apt-get install -y -qq coinor-cbc\n",
    "    else:\n",
    "        try:\n",
    "            !conda install -c conda-forge coincbc \n",
    "        except:\n",
    "            pass\n",
    "\n",
    "assert(shutil.which(\"cbc\") or os.path.isfile(\"cbc\"))\n",
    "from pyomo.environ import *\n",
    "from pyomo.gdp import *\n",
    "import pandas as pd"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "0AFDvlnda6br"
   },
   "source": [
    "## Problem statement"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "eApCyKRVAPy0"
   },
   "source": [
    "### Component data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 142
    },
    "executionInfo": {
     "elapsed": 648,
     "status": "ok",
     "timestamp": 1603203431565,
     "user": {
      "displayName": "Jeffrey Kantor",
      "photoUrl": "https://lh3.googleusercontent.com/a-/AOh14Gg_n8V7bVINy02QRuRgOoMo11Ri7NKU3OUKdC1bkQ=s64",
      "userId": "09038942003589296665"
     },
     "user_tz": 240
    },
    "id": "YQXecWx9ZwRp",
    "outputId": "6ea35497-ac06-480c-d0bd-71611bc78f7f"
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>cost</th>\n",
       "      <th>Vit A</th>\n",
       "      <th>Vit B</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>A</th>\n",
       "      <td>2.0</td>\n",
       "      <td>0.5</td>\n",
       "      <td>0.2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>B</th>\n",
       "      <td>2.0</td>\n",
       "      <td>0.4</td>\n",
       "      <td>0.1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>C</th>\n",
       "      <td>5.0</td>\n",
       "      <td>0.3</td>\n",
       "      <td>0.3</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   cost  Vit A  Vit B\n",
       "A   2.0    0.5    0.2\n",
       "B   2.0    0.4    0.1\n",
       "C   5.0    0.3    0.3"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# load data as dictionary of components\n",
    "# component data consists of cost and composition \n",
    "comp_data = {\n",
    "    \"A\": {\"cost\": 2.0, \"Vit A\": 0.5, \"Vit B\": 0.2},\n",
    "    \"B\": {\"cost\": 2.0, \"Vit A\": 0.4, \"Vit B\": 0.1},\n",
    "    \"C\": {\"cost\": 5.0, \"Vit A\": 0.3, \"Vit B\": 0.3},\n",
    "}\n",
    "\n",
    "# use pandas to create a nice display\n",
    "pd.DataFrame.from_dict(comp_data, orient='index')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "6uQWK9n4a933"
   },
   "source": [
    "### Product Composition Requirements\n",
    "\n",
    "Find the lowest cost blend\n",
    "\n",
    "* Vit A: less than 0.4\n",
    "* Vit B: greater than 0.2\n",
    "\n",
    "Your code should be able to accept alternative specification for data and product requirements.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 111
    },
    "executionInfo": {
     "elapsed": 293,
     "status": "ok",
     "timestamp": 1603203434898,
     "user": {
      "displayName": "Jeffrey Kantor",
      "photoUrl": "https://lh3.googleusercontent.com/a-/AOh14Gg_n8V7bVINy02QRuRgOoMo11Ri7NKU3OUKdC1bkQ=s64",
      "userId": "09038942003589296665"
     },
     "user_tz": 240
    },
    "id": "kaiRnfEQdJm1",
    "outputId": "28b73323-e120-4e70-e976-8884b0e195b1"
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>lb</th>\n",
       "      <th>ub</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>Vit A</th>\n",
       "      <td>0.0</td>\n",
       "      <td>0.4</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Vit B</th>\n",
       "      <td>0.2</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "        lb   ub\n",
       "Vit A  0.0  0.4\n",
       "Vit B  0.2  1.0"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "prod_req = {\n",
    "    \"Vit A\": {\"lb\": 0.0, \"ub\": 0.4},\n",
    "    \"Vit B\": {\"lb\": 0.2, \"ub\": 1.0},\n",
    "}\n",
    "pd.DataFrame.from_dict(prod_req, orient='index')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "w0MKxgb0cKtx"
   },
   "source": [
    "### Component Compatibility\n",
    "\n",
    "For this application, we consider an additional type of constraint specifying the incompatability of certain blends of components. For example, suppose we have a constraint:\n",
    "\n",
    "* A and B cannot be mixed together in the final product\n",
    "\n",
    "The constraint is specified by creating a list of incompatabile pairs."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "id": "YA1C87la47CS"
   },
   "outputs": [],
   "source": [
    "excl_pairs = [(\"A\", \"B\")]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "O1yHzjPYUAaD"
   },
   "source": [
    "## Version 0: Neglecting the compatibility requirments"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 68
    },
    "executionInfo": {
     "elapsed": 288,
     "status": "ok",
     "timestamp": 1603203567142,
     "user": {
      "displayName": "Jeffrey Kantor",
      "photoUrl": "https://lh3.googleusercontent.com/a-/AOh14Gg_n8V7bVINy02QRuRgOoMo11Ri7NKU3OUKdC1bkQ=s64",
      "userId": "09038942003589296665"
     },
     "user_tz": 240
    },
    "id": "Enqmwye9TzVB",
    "outputId": "1dd776c2-2b1b-4009-a24c-e7e0dc25e842"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "A = 0.33333333\n",
      "B = 0.33333333\n",
      "C = 0.33333333\n"
     ]
    }
   ],
   "source": [
    "m = ConcreteModel()\n",
    "\n",
    "# define sets that will be used to index decision variables and constraints\n",
    "# remember to use initialize keyword\n",
    "m.comp = Set(initialize=comp_data.keys())\n",
    "m.req = Set(initialize=prod_req.keys())\n",
    "\n",
    "# decision variables\n",
    "m.x = Var(m.comp, domain=NonNegativeReals)\n",
    "\n",
    "# objective function\n",
    "m.cost = Objective(expr=sum(m.x[c]*comp_data[c][\"cost\"] for c in m.comp), sense=minimize)\n",
    "\n",
    "# structural constraints\n",
    "m.massfraction = Constraint(expr=sum(m.x[c] for c in m.comp)==1)\n",
    "\n",
    "# composition constraints\n",
    "m.lb = Constraint(m.req, rule=lambda m, r: sum(m.x[c]*comp_data[c][r] for c in m.comp) >= prod_req[r][\"lb\"])\n",
    "m.ub = Constraint(m.req, rule=lambda m, r: sum(m.x[c]*comp_data[c][r] for c in m.comp) <= prod_req[r][\"ub\"])\n",
    "\n",
    "solver = SolverFactory('cbc')\n",
    "solver.solve(m)\n",
    "\n",
    "for c in m.comp:\n",
    "    print(f\"{c} = {m.x[c]()}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "byX-oOwvc5cg"
   },
   "source": [
    "## Version 1: Including compatibility requirements with Big-M\n",
    "\n",
    "The challenge of this problem are the disjunctive constraints associated with the component incompatability data. Here we associated a boolean variable for each pair, then use the boolean variable to determine which member of the pair to keep in the blend."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 68
    },
    "executionInfo": {
     "elapsed": 498,
     "status": "ok",
     "timestamp": 1603203949801,
     "user": {
      "displayName": "Jeffrey Kantor",
      "photoUrl": "https://lh3.googleusercontent.com/a-/AOh14Gg_n8V7bVINy02QRuRgOoMo11Ri7NKU3OUKdC1bkQ=s64",
      "userId": "09038942003589296665"
     },
     "user_tz": 240
    },
    "id": "Pb1QoE0dZ4LO",
    "outputId": "b9183569-eaec-4544-c98e-52130a6a34a7"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "A = 0.0\n",
      "B = 0.5\n",
      "C = 0.5\n"
     ]
    }
   ],
   "source": [
    "m = ConcreteModel()\n",
    "\n",
    "# define sets that will be used to index decision variables and constraints\n",
    "# remember to use initialize keyword\n",
    "m.comp = Set(initialize=comp_data.keys())\n",
    "m.req = Set(initialize=prod_req.keys())\n",
    "\n",
    "# define a set to that includes the excluded pairs\n",
    "m.pairs = Set(initialize=excl_pairs)\n",
    "\n",
    "# decision variables\n",
    "m.x = Var(m.comp, domain=NonNegativeReals)\n",
    "\n",
    "# for each excluded pair, create a boolean variable. The value of the boolean\n",
    "# variable will determine which member of the pair is allowed in the product\n",
    "m.y = Var(m.pairs, domain=Boolean)\n",
    "\n",
    "# objective function\n",
    "m.cost = Objective(expr=sum(m.x[c]*comp_data[c][\"cost\"] for c in m.comp), sense=minimize)\n",
    "\n",
    "# structural constraints\n",
    "m.massfraction = Constraint(expr=sum(m.x[c] for c in m.comp)==1)\n",
    "\n",
    "# composition constraints\n",
    "m.lb = Constraint(m.req, rule=lambda m, r: sum(m.x[c]*comp_data[c][r] for c in m.comp) >= prod_req[r][\"lb\"])\n",
    "m.ub = Constraint(m.req, rule=lambda m, r: sum(m.x[c]*comp_data[c][r] for c in m.comp) <= prod_req[r][\"ub\"])\n",
    "\n",
    "# component incompatability constraints\n",
    "M = 100\n",
    "m.disj = ConstraintList()\n",
    "for pair in m.pairs:\n",
    "    a, b = pair\n",
    "    m.disj.add(m.x[a] <= M*m.y[pair])\n",
    "    m.disj.add(m.x[b] <= M*(1-m.y[pair]))\n",
    "\n",
    "solver = SolverFactory('cbc')\n",
    "solver.solve(m)\n",
    "\n",
    "for c in m.comp:\n",
    "    print(f\"{c} = {m.x[c]()}\")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "lXTRureoeN8-"
   },
   "source": [
    "## Version 2. Disjunctive Constraints"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 68
    },
    "executionInfo": {
     "elapsed": 335,
     "status": "ok",
     "timestamp": 1603204548854,
     "user": {
      "displayName": "Jeffrey Kantor",
      "photoUrl": "https://lh3.googleusercontent.com/a-/AOh14Gg_n8V7bVINy02QRuRgOoMo11Ri7NKU3OUKdC1bkQ=s64",
      "userId": "09038942003589296665"
     },
     "user_tz": 240
    },
    "id": "QqrlorbPw3Yw",
    "outputId": "c8ed97b4-6d5c-415e-f7e5-fe3d4f391627"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "A = 0.0\n",
      "B = 0.5\n",
      "C = 0.5\n"
     ]
    }
   ],
   "source": [
    "m = ConcreteModel()\n",
    "\n",
    "# define sets that will be used to index decision variables and constraints\n",
    "# remember to use initialize keyword\n",
    "m.comp = Set(initialize=comp_data.keys())\n",
    "m.req = Set(initialize=prod_req.keys())\n",
    "\n",
    "# define a set to that includes the excluded pairs\n",
    "m.pairs = Set(initialize=excl_pairs)\n",
    "\n",
    "# decision variables\n",
    "m.x = Var(m.comp, domain=NonNegativeReals, bounds=(0, 1))\n",
    "\n",
    "# objective function\n",
    "m.cost = Objective(expr=sum(m.x[c]*comp_data[c][\"cost\"] for c in m.comp), sense=minimize)\n",
    "\n",
    "# structural constraints\n",
    "m.massfraction = Constraint(expr=sum(m.x[c] for c in m.comp)==1)\n",
    "\n",
    "# composition constraints\n",
    "m.lb = Constraint(m.req, rule=lambda m, r: sum(m.x[c]*comp_data[c][r] for c in m.comp) >= prod_req[r][\"lb\"])\n",
    "m.ub = Constraint(m.req, rule=lambda m, r: sum(m.x[c]*comp_data[c][r] for c in m.comp) <= prod_req[r][\"ub\"])\n",
    "\n",
    "# component incompatability constraints\n",
    "m.disj = Disjunction(m.pairs, rule=lambda m, a, b: [m.x[a] == 0, m.x[b] == 0])\n",
    "\n",
    "# apply transformations\n",
    "TransformationFactory('gdp.hull').apply_to(m)\n",
    "\n",
    "# solve\n",
    "solver = SolverFactory('cbc')\n",
    "solver.solve(m)\n",
    "\n",
    "for c in m.comp:\n",
    "    print(f\"{c} = {m.x[c]()}\")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "Sbpc_IPaC_5X"
   },
   "source": [
    "## Analysis "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "id": "keeCKOCDDnIV"
   },
   "outputs": [],
   "source": [
    "comp_data = {\n",
    "    \"A\": {\"cost\": 2.0, \"Vit A\": 0.5, \"Vit B\": 0.2},\n",
    "    \"B\": {\"cost\": 2.0, \"Vit A\": 0.4, \"Vit B\": 0.1},\n",
    "    \"C\": {\"cost\": 4.0, \"Vit A\": 0.3, \"Vit B\": 0.3},\n",
    "}\n",
    "\n",
    "prod_req = {\n",
    "    \"Vit A\": {\"lb\": 0.0, \"ub\": 0.4},\n",
    "    \"Vit B\": {\"lb\": 0.2, \"ub\": 1.0},\n",
    "}\n",
    "\n",
    "excl_pairs = [(\"A\", \"B\")]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "tBg0H7vFEPXK"
   },
   "source": [
    "$$\n",
    "\\begin{align*}\n",
    "x_A + x_B + x_C & = 1 \\\\\n",
    "0.5 x_A + 0.4 x_B + 0.3 x_C & \\leq 0.4 \\\\\n",
    "0.2 x_A + 0.1 x_B + 0.3 x_C & \\geq 0.2 \\\\\n",
    "\\end{align*}\n",
    "$$\n",
    "\n",
    "Solving for x_C\n",
    "\n",
    "$$\n",
    "\\begin{align*}\n",
    "x_C & = 1 - x_A - x_B\n",
    "\\end{align*}\n",
    "$$\n",
    "\n",
    "Substitution\n",
    "\n",
    "$$\n",
    "\\begin{align*}\n",
    "0.2 x_A + 0.1 x_B & \\leq 0.1 \\\\\n",
    "-0.1 x_A - 0.2 x_B & \\geq -0.1 \\\\\n",
    "\\end{align*}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "id": "iiyVCd3Pzik9"
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<pyomo.common.factory.Factory at 0x7fb04e500100>"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "TransformationFactory"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "colab": {
   "collapsed_sections": [],
   "name": "Introduction-to-Disjunctive-Programming.ipynb",
   "provenance": [],
   "toc_visible": true
  },
  "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.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}