convex/linear/product_mix.ipynb

{
 "cells": [
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Product mix"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Overview\n",
    "\n",
    "In this problem, consider a set of resources $I$ used to produce a set of products $J$. Each product demands predefined ammounts of each resource. Profit margins per unit of each product are known.\n",
    "\n",
    "Let us denote $c_{j}$ the margin associated to each product $j \\in J$ and $b_{i}$ each reactant $i \\in I$ availability. The ammount of $i$ used per unit of $j$ is given by $f_{i, j}$. Our objective is to maximize operational profit margin."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\begin{align}\n",
    "    \\text{max} \\quad & \\sum_{j \\in J} c_j x_j \\\\\n",
    "    \\text{s.t.} \\quad & \\sum_{j \\in J} f_{i, j} x_{j} \\leq b_{i} & \\forall \\; i \\in I \\\\\n",
    "    & x_{j} \\geq 0 & \\forall \\; j \\in J \\\\\n",
    "\\end{align}\n",
    "$$"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example\n",
    "\n",
    "In this example, we consider six chemical species:\n",
    "- Reactants {A, B, C}\n",
    "- Products {D, E, F}.\n",
    "\n",
    "Table 1: Reactants\n",
    "| Reactant | Availability (lb/day) |\n",
    "| --- | --- |\n",
    "| A | 8000 |\n",
    "| B | 3000 |\n",
    "| C | 2500 |\n",
    "\n",
    " \n",
    "Table 2: Products\n",
    "| Product | Unitary margin ($/lb) |\n",
    "| --- | --- |\n",
    "| D | 2.15 |\n",
    "| E | 1.34 |\n",
    "| F | 1.72 |\n",
    "\n",
    "Table 3: Proportions\n",
    "| Reactant \\ Product | D | E | F\n",
    "| --- | --- | --- | --- |\n",
    "| A | 7/10 | 1/3 | 1/2 |\n",
    "| B | 1/5 | 2/3 | 1/6 |\n",
    "| C | 1/10 | 0 | 1/3 |\n",
    "\n",
    "\n",
    "Our goal is to maximize daily profit.\n",
    "\n",
    "Let us denote $c_{j}$ the margin associated to each product $j \\in J$ and $b_{i}$ each reactant $i \\in I$ availability. The ammount of $i$ used per unit of $j$ is given by $f_{i, j}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from scipy.optimize import linprog\n",
    "import pyomo.environ as pyo"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Scipy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Pseudo costs\n",
    "margins = np.array([2.15, 1.34, 1.72])\n",
    "c = - margins\n",
    "\n",
    "# A matrix\n",
    "A = np.array([\n",
    "    [7/10, 1/3, 1/2],\n",
    "    [1/5, 2/3, 1/6],\n",
    "    [1/10, 0.0, 1/3]\n",
    "])\n",
    "\n",
    "b = np.array([8000.0, 3000.0, 2500.0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[7105.26315789 1026.31578947 5368.42105263]\n"
     ]
    }
   ],
   "source": [
    "sol = linprog(c, A_ub=A, b_ub=b)\n",
    "print(sol.x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Pyomo"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Initialize model\n",
    "model = pyo.ConcreteModel()\n",
    "\n",
    "# Sets of reactants and products\n",
    "model.I = pyo.Set(initialize=[\"A\", \"B\", \"C\"])\n",
    "model.J = pyo.Set(initialize=[\"D\", \"E\", \"F\"])\n",
    "\n",
    "# Availability\n",
    "model.b = pyo.Param(model.I, initialize=dict(zip(model.I, b)))\n",
    "\n",
    "# Margins (now as a maximization objective)\n",
    "model.c = pyo.Param(model.J, initialize=dict(zip(model.J, margins)))\n",
    "\n",
    "# Proportions\n",
    "proportions = {}\n",
    "for k, i in enumerate(model.I):\n",
    "    for l, j in enumerate(model.J):\n",
    "        proportions[i, j] = A[k, l]\n",
    "\n",
    "model.f = pyo.Param(model.I, model.J, initialize=proportions)\n",
    "\n",
    "# Decision variables\n",
    "model.x = pyo.Var(model.J, within=pyo.NonNegativeReals, name=\"x\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Availability constraint\n",
    "def availability_rule_cstr(model, i):\n",
    "    return sum(model.f[i, j] * model.x[j] for j in model.J) <= model.b[i]\n",
    "\n",
    "model.cstr_available = pyo.Constraint(model.I, rule=availability_rule_cstr, name=\"cstr_available\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Objective function\n",
    "def obj_func(model):\n",
    "    return sum(model.x[j] * model.c[j] for j in model.J)\n",
    "\n",
    "model.obj = pyo.Objective(rule=obj_func, sense=pyo.maximize)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "cbc = pyo.SolverFactory(\"cbc\")\n",
    "sol = cbc.solve(model, tee=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "x : Size=3, Index=J\n",
      "    Key : Lower : Value     : Upper : Fixed : Stale : Domain\n",
      "      D :     0 : 7105.2632 :  None : False : False : NonNegativeReals\n",
      "      E :     0 : 1026.3158 :  None : False : False : NonNegativeReals\n",
      "      F :     0 : 5368.4211 :  None : False : False : NonNegativeReals\n"
     ]
    }
   ],
   "source": [
    "model.x.display()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "cstr_available : Size=3\n",
      "    Key : Lower : Body               : Upper\n",
      "      A :  None :  8000.000056666666 : 8000.0\n",
      "      B :  None : 3000.0000233333335 : 3000.0\n",
      "      C :  None :         2500.00002 : 2500.0\n"
     ]
    }
   ],
   "source": [
    "model.cstr_available.display()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "obj : Size=1, Index=None, Active=True\n",
      "    Key  : Active : Value\n",
      "    None :   True : 25885.263344\n"
     ]
    }
   ],
   "source": [
    "model.obj.display()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3.10.6 ('venv': venv)",
   "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"
  },
  "orig_nbformat": 4,
  "vscode": {
   "interpreter": {
    "hash": "74cc81cb9fe1c1e193d3e0c529c6f5e8b7e2394a8ce57de487f07b5929c4a8c2"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}