{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "229dff1b-612f-46a6-8ae8-9c7581cb1782",
   "metadata": {},
   "source": [
    "# Production model using disjunctions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bb501ddb-f145-4b79-96f0-69fa1403467c",
   "metadata": {},
   "outputs": [],
   "source": [
    "# install dependencies and select solver\n",
    "%pip install -q amplpy\n",
    "\n",
    "SOLVER = \"highs\"\n",
    "\n",
    "from amplpy import AMPL, ampl_notebook\n",
    "\n",
    "ampl = ampl_notebook(\n",
    "    modules=[\"highs\"],  # modules to install\n",
    "    license_uuid=\"default\",  # license to use\n",
    ")  # instantiate AMPL object and register magics"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cf9db689-23e6-4d6d-8dd3-68149c42301c",
   "metadata": {},
   "source": [
    "## Disjunctions\n",
    "\n",
    "Disjunctions appear in applications where there is choice among discrete alternatives. Given two logical propositions $\\alpha$ and $\\beta$, the \"or\" disjunction is denoted by $\\vee$ and defined by the truth table\n",
    "\n",
    "| $\\alpha$ | $\\beta$ | $\\alpha \\vee \\beta$ |\n",
    "| :-: | :-: | :-: |\n",
    "| False | False | False |\n",
    "| True | False | True |\n",
    "| False | True | True |\n",
    "| True | True | True |\n",
    "\n",
    "The \"exclusive or\" is denoted by $\\veebar$ and defined by the truth table\n",
    "\n",
    "| $\\alpha$ | $\\beta$ | $\\alpha \\veebar \\beta$ |\n",
    "| :-: | :-: | :-: |\n",
    "| False | False | False |\n",
    "| True | False | True |\n",
    "| False | True | True |\n",
    "| True | True | False |\n",
    "\n",
    "This notebook shows how to express disjunctions in AMPL models using the Generalized Disjunctive Programming syntax for a simple production model.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c1e98d83-34eb-45a4-a6e6-5eb733f80fcc",
   "metadata": {
    "tags": []
   },
   "source": [
    "## Multi-product factory optimization\n",
    "\n",
    "A small production facility produces two products, $X$ and $Y$. With current technology $\\alpha$, the facility is subject to the following conditions and constraints:\n",
    "\n",
    "* Product $X$ requires 1 hour of labor A, 2 hours of labor B, and 100€ of raw material. Product $X$ sells for 270€ per unit. The daily demand is limited to 40 units.\n",
    "\n",
    "* Product $Y$ requires 1 hour of labor A, 1 hour of labor B, and 90€ of raw material. Product $Y$ sells for 210€ per unit with unlimited demand. \n",
    "\n",
    "* There are 80 hours per day of labor A available at a cost of 50€/hour.\n",
    "\n",
    "* There are 100 hours per day of labor B available at a cost of 40€/hour.\n",
    "\n",
    "Using the given data we see that the net profit for each unit of $X$ and $Y$ is 40€ and 30€, respectively. The optimal product strategy is the solution to a linear optimization\n",
    "\n",
    "$$\n",
    "\\begin{align*}\n",
    "\\max_{x, y \\geq 0} \\quad & \\text{profit}\\\\\n",
    "\\text{s.t.} \\quad \n",
    "& \\text{profit}  = 40 x + 30 y\\\\\n",
    "& x  \\leq 40 & \\text{(demand)}\\\\\n",
    "& x + y  \\leq 80 & \\text{(labor A)} \\\\\n",
    "& 2 x + y  \\leq 100 & \\text{(labor B)}\n",
    "\\end{align*}\n",
    "$$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "13174feb",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Overwriting multi_product_factory.mod\n"
     ]
    }
   ],
   "source": [
    "%%writefile multi_product_factory.mod\n",
    "\n",
    "# decision variables\n",
    "var profit;\n",
    "var production_x >= 0;\n",
    "var production_y >= 0;\n",
    "\n",
    "# profit objective\n",
    "maximize maximize_profit: profit;\n",
    "    \n",
    "# constraints\n",
    "s.t. profit_expr: profit == 40 * production_x + 30 * production_y;\n",
    "s.t. demand: production_x <= 40;\n",
    "s.t. laborA: production_x + production_y <= 80;\n",
    "s.t. laborB: 2 * production_x + production_y <= 100;"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "0291f0fe-a9f6-4427-900a-70bd0e5ac916",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "HiGHS 1.5.1: \b\b\b\b\b\b\b\b\b\b\b\b\bHiGHS 1.5.1: optimal solution; objective 2600\n",
      "2 simplex iterations\n",
      "0 barrier iterations\n",
      "Profit = 2600.00 €\n",
      "Production X = 20.0\n",
      "Production Y = 60.0\n"
     ]
    }
   ],
   "source": [
    "m = AMPL()\n",
    "m.read(\"multi_product_factory.mod\")\n",
    "m.option[\"solver\"] = SOLVER\n",
    "m.solve()\n",
    "\n",
    "print(f\"Profit = {m.var['profit'].value():.2f} €\")\n",
    "print(f\"Production X = {m.var['production_x'].value()}\")\n",
    "print(f\"Production Y = {m.var['production_y'].value()}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "017223c6-7b4d-42c3-8af3-64e6586594a2",
   "metadata": {},
   "source": [
    "Now suppose a new technology $\\beta$ is available that affects that lowers the cost of product $X$. With the new technology, only 1.5 hours of labor B is required per unit of $X$.\n",
    "\n",
    "The net profit for unit of product $X$ with technology $\\alpha$ is equal to $270 - 100 - 50 - 2 \\cdot 40 = 40$€\n",
    "\n",
    "The net profit for unit of product $X$ with technology $\\beta$ is equal to $270 - 100 - 50 - 1.5 \\cdot 40 = 60$€\n",
    "\n",
    "The decision here is whether to use technology $\\alpha$ or $\\beta$. There are several commonly used techniques for embedding disjunctions into mixed-integer linear optimization problems. The \"big-M\" technique introduces a binary decision variable for every exclusive-or disjunction between two constraints. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3a4f3501-72d0-4908-b95c-1b4c373abc59",
   "metadata": {},
   "source": [
    "## MILO implementation\n",
    "\n",
    "We can formulate this problem as the following mixed-integer linear optimization (MILO) problem:\n",
    "\n",
    "$$\n",
    "\\begin{align*}\n",
    "    \\max_{x, y \\geq 0, z \\in \\mathbb{B}} \\quad & \\text{profit}\\\\\n",
    "    \\text{s.t.} \\quad \n",
    "    & x  \\leq 40 & \\text{(demand)}\\\\\n",
    "    & x + y  \\leq 80 & \\text{(labor A)} \\\\\n",
    "    & \\text{profit} \\leq 40x + 30y + M z \\\\\n",
    "    & \\text{profit} \\leq 60x + 30y + M (1 - z) \\\\\n",
    "    & 2 x + y \\leq 100  + M z \\\\ \n",
    "    & 1.5 x + y \\leq 100 + M (1 - z).\n",
    "\\end{align*}\n",
    "$$\n",
    "\n",
    "where the variable $z \\in \\{ 0, 1\\}$ \"activates\" the constraints related to the old or new technology, respectively, and $M$ is a big enough number. The corresponding AMPL implementation is given by:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "66999b5c",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Overwriting multi_product_factory_milo.mod\n"
     ]
    }
   ],
   "source": [
    "%%writefile multi_product_factory_milo.mod\n",
    "\n",
    "# decision variables\n",
    "var profit;\n",
    "var production_x >= 0;\n",
    "var production_y >= 0;\n",
    "var z binary;\n",
    "param M = 10000;\n",
    "\n",
    "# profit objective\n",
    "maximize maximize_profit: profit;\n",
    "\n",
    "# constraints\n",
    "s.t. profit_constr_1: profit <= 40 * production_x + 30 * production_y + M * z;\n",
    "s.t. profit_constr_2: profit <= 60 * production_x + 30 * production_y + M * (1 - z);\n",
    "s.t. demand: production_x <= 40;\n",
    "s.t. laborA: production_x + production_y <= 80;\n",
    "s.t. laborB_1: 2 * production_x + production_y <= 100 + M * z;\n",
    "s.t. laborB_2: 1.5 * production_x + production_y <= 100 + M * (1 - z);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "5e9a7dc7-d427-4f85-b0de-bb7e6db4af8d",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "HiGHS 1.5.1: \b\b\b\b\b\b\b\b\b\b\b\b\bHiGHS 1.5.1: optimal solution; objective 3600\n",
      "8 simplex iterations\n",
      "1 branching nodes\n",
      "Profit = 3600.00 €\n",
      "Production X = 40.0\n",
      "Production Y = 40.0\n"
     ]
    }
   ],
   "source": [
    "m = AMPL()\n",
    "m.read(\"multi_product_factory_milo.mod\")\n",
    "m.option[\"solver\"] = SOLVER\n",
    "m.solve()\n",
    "\n",
    "print(f\"Profit = {m.var['profit'].value():.2f} €\")\n",
    "print(f\"Production X = {m.var['production_x'].value()}\")\n",
    "print(f\"Production Y = {m.var['production_y'].value()}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "31dfdc49-c671-4cf6-966e-74985de748e1",
   "metadata": {},
   "source": [
    "## Disjunctive programming implementation\n",
    "\n",
    "Alternatively, we can formulate our problem using a disjunction, preserving the logical structure, as follows:\n",
    "\n",
    "$$\n",
    "\\begin{align*}\n",
    "\\max_{x, y \\geq 0} \\quad & \\text{profit}\\\\\n",
    "\\text{s.t.} \\quad \n",
    "& x  \\leq 40 & \\text{(demand)}\\\\\n",
    "& x + y  \\leq 80 & \\text{(labor A)} \\\\\n",
    "& \\begin{bmatrix}\n",
    "    \\text{profit} = 40x + 30y\\\\\n",
    "    2 x + y \\leq 100\n",
    "\\end{bmatrix}\n",
    " \\veebar\n",
    "\\begin{bmatrix}\n",
    "    \\text{profit} = 60x + 30y\\\\\n",
    "    1.5 x + y \\leq 100\n",
    "    \\end{bmatrix}\n",
    "\\end{align*}\n",
    "$$\n",
    "\n",
    "This formulation, if allowed by the software at hand, has the benefit that the software can smartly divide the solution of this problem into sub-possibilities depending on the disjunction. AMPL supports disjunctions, as illustrated in the following implementation:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "1c312cb5",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Overwriting multi_product_factory_milo_disjunctive.mod\n"
     ]
    }
   ],
   "source": [
    "%%writefile multi_product_factory_milo_disjunctive.mod\n",
    "\n",
    "var profit >= -1000, <= 10000;\n",
    "var x >= 0, <= 1000;\n",
    "var y >= 0, <= 1000;\n",
    "\n",
    "maximize maximize_profit: profit;\n",
    "\n",
    "s.t. demand: x <= 40;\n",
    "s.t. laborA: x + y <= 80;\n",
    "s.t. technologies:\n",
    "    ((profit == 40 * x + 30 * y and 2 * x + y <= 100)\n",
    "    or\n",
    "    (profit == 60 * x + 30 * y and 1.5 * x + y <= 100))\n",
    "    and not\n",
    "    ((profit == 40 * x + 30 * y and 2 * x + y <= 100)\n",
    "    and\n",
    "    (profit == 60 * x + 30 * y and 1.5 * x + y <= 100))\n",
    "    ;"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "7e0dcacb-16d3-4031-99b0-e10c7b68050e",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "HiGHS 1.5.1: \b\b\b\b\b\b\b\b\b\b\b\b\bHiGHS 1.5.1: optimal solution; objective 3600\n",
      "12 simplex iterations\n",
      "1 branching nodes\n",
      "Profit = 3600.00 €\n",
      "x = 40.0\n",
      "y = 40.0\n"
     ]
    }
   ],
   "source": [
    "m = AMPL()\n",
    "m.read(\"multi_product_factory_milo_disjunctive.mod\")\n",
    "m.option[\"solver\"] = SOLVER\n",
    "m.solve()\n",
    "\n",
    "print(f\"Profit = {m.var['profit'].value():.2f} €\")\n",
    "print(f\"x = {m.var['x'].value()}\")\n",
    "print(f\"y = {m.var['y'].value()}\")"
   ]
  }
 ],
 "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.10.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}