{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Generation/Consumption Analysis by Linear Programming\n",
    "\n",
    "by Jeffrey Kantor  (jeff at nd.edu). The latest version of this notebook is available at [https://github.com/jckantor/CBE20255](https://github.com/jckantor/CBE20255). \n",
    "\n",
    "### Summary\n",
    "\n",
    "Generation/Consumption analysis combines technologically feasible chemical reactions into processes for the economic production of chemical goods. This notebook demonstrates generation/consumption analysis using linear programming techniques in Python. The example uses the `PuLP`, a linear programming toolkit for linear programming."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Problem Description\n",
    "\n",
    "### Solar Thermochemical Production of Ammonia\n",
    "\n",
    "Ammonia has a wide range of commercial uses, including an essential role as a fertilizer for the global agriculture. The principle means of production is the Haber/Bosch process where hydrogen and nitrogren react at high temperature and pressure over an iron catalyst, and the required hydrogen is produced by the steam reforming of methane from natural gas. Globally, this process consumes 5% of all natural gas produced, and 2% of all energy resources.\n",
    "\n",
    "Recently, several research groups have proposed alternative routes to the production of ammonia using solar radiation as the major source of process energy. Among them, the [two step production of ammonia via Al_2O_3/AlN\\\\)]() offers the advantage of atmospheric operation production of a synthesis gas by-product. The process consists of the highly endothermic carbothermic reduction of $Al_2O_3$.\n",
    "\n",
    "$$Al_2O_3 + 3\\,CH_4 + N_2 \\longrightarrow 2\\,AlN + 6\\,H_2 + 3\\,CO$$\n",
    "\n",
    "followed by the mildly exothermic hydrolysis of $AlN$\n",
    "\n",
    "$$2\\,AlN + 3\\,H_2O \\longrightarrow Al_2O_3 + 2\\,NH_3$$\n",
    "\n",
    "The by-product of these two reactions is a mixture of carbon monoxide and hydrogen which is generally referred to as synthesis gas.\n",
    "\n",
    "### Production of Liquid Transportation Fuel from Synthesis Gas\n",
    "\n",
    "A variety of liquid fuel products can be made from synthesis gas. Here we'll consider the production dimethyl ether $CH_3OCH_3$ (DME) which can be used an ultraclean alternative to diesel or as a blending component for gasoline. DME can be produced in a single stage over a bi-functional catalyst via the reactions\n",
    "\n",
    "$$CO_2 + 3\\,H_2 \\longrightarrow CH_3OH + H_2O$$\n",
    "$$CO + H_2O \\longrightarrow CO_2 + H_2$$\n",
    "$$2\\,CH_3OH \\longrightarrow CH_3OCH_3 + H_2O$$\n",
    "\n",
    "### Process Objectives\n",
    "\n",
    "Considering just these five reactions, is it possible to produce ammonia, and a mixture of DME and methanol without any by-product carbon dioxide? What is the net process stoichiometry?\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Solution by Linear Programming"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "| Species  | MW |$R_1$|$R_2$|$R_3$|$R_4$|$R_5$| Net      |\n",
    "| :-----: | :-: | :-: | :-: | :-: | :-: | :-: | :-:      |\n",
    "| $Al_2O_3$   | | -1  |  1  |     |     |     | 0        |\n",
    "| $AlN$       | |  2  | -2  |     |     |     | 0        |\n",
    "| $CH_4$      | | -3  |     |     |     |     | $\\leq$ 0 |\n",
    "| $N_2$       | | -1  |     |     |     |     |          |\n",
    "| $H_2$       | |  6  |     | -3  |  1  |     | $\\geq$ 0 |\n",
    "| $CO$        | |  3  |     |     | -1  |     | 0        |\n",
    "| $H_2O$      | |     | -3  |  1  | -1  |  1  |          |\n",
    "| $NH_3$      | |     |  2  |     |     |     | $\\geq$ 0 |\n",
    "| $CO_2$      | |     |     | -1  |  1  |     | 0        |\n",
    "| $CH_3OH$    | |     |     |  1  |     | -2  | $\\geq$ 0 |\n",
    "| $CH_3OCH_3$ | |     |     |     |     |  1  | 1        |   "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Create a dictionary for the stoichiometric coefficients in each reaction."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 205,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'R1': {'Al2O3': -1, 'AlN': 2, 'CH4': -3, 'CO': 3, 'H2': 6, 'N2': -1},\n",
       " 'R2': {'Al2O3': 1, 'AlN': -2, 'H2O': -3, 'NH3': 2},\n",
       " 'R3': {'CH3OH': 1, 'CO2': -1, 'H2': -3, 'H2O': 1},\n",
       " 'R4': {'CO': -1, 'CO2': 1, 'H2': 1, 'H2O': -1},\n",
       " 'R5': {'CH3OCH3': 1, 'CH3OH': -2, 'H2O': 1}}"
      ]
     },
     "execution_count": 205,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "nu = dict()\n",
    "\n",
    "nu['R1'] = {'Al2O3':-1, 'CH4':-3, 'N2':-1, 'AlN':2, 'H2':6, 'CO':3}\n",
    "nu['R2'] = {'AlN':-2, 'H2O':-3, 'Al2O3':1, 'NH3':2}\n",
    "nu['R3'] = {'CO2':-1, 'H2':-3, 'CH3OH':1, 'H2O':1}\n",
    "nu['R4'] = {'CO':-1, 'H2O':-1, 'CO2':1, 'H2':1}\n",
    "nu['R5'] = {'CH3OH':-2, 'CH3OCH3':1, 'H2O':1}\n",
    "nu"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "From the stoichiometric dictionary we can create a set of reactions and species."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 206,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "rxns = set(nu.keys())\n",
    "species = set([s for r in rxns for s in nu[r].keys()])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 207,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "R1: N2 + 3 CH4 + Al2O3 --> 3 CO + 6 H2 + 2 AlN\n",
      "R2: 3 H2O + 2 AlN --> Al2O3 + 2 NH3\n",
      "R3: CO2 + 3 H2 --> H2O + CH3OH\n",
      "R4: H2O + CO --> CO2 + H2\n",
      "R5: 2 CH3OH --> H2O + CH3OCH3\n"
     ]
    }
   ],
   "source": [
    "def printRxn(nu):\n",
    "    c2s = lambda n: '' if n==1 else str(round(n,3))+' '\n",
    "    lhs = ' + '.join([\"{0:s}{1:s}\".format(c2s(-nu[s]),s) for s in nu.keys() if nu[s] < 0])\n",
    "    rhs = ' + '.join([\"{0:s}{1:s}\".format(c2s(nu[s]),s) for s in nu.keys() if nu[s] > 0])\n",
    "    print(' --> '.join([lhs,rhs]))\n",
    "    \n",
    "for r in sorted(rxns):\n",
    "    print(r,end=': ')\n",
    "    printRxn(nu[r])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 227,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "           R1   R2   R3   R4   R5\n",
      "Al2O3      -1    1               \n",
      "AlN         2   -2               \n",
      "CH3OCH3                         1\n",
      "CH3OH                 1        -2\n",
      "CH4        -3                    \n",
      "CO          3             -1     \n",
      "CO2                  -1    1     \n",
      "H2          6        -3    1     \n",
      "H2O             -3    1   -1    1\n",
      "N2         -1                    \n",
      "NH3              2               \n"
     ]
    }
   ],
   "source": [
    "print('{0:8s}'.format(''),end='')\n",
    "print(''.join(['{0:>5s}'.format(r) for r in sorted(rxns)]))\n",
    "for s in sorted(species):\n",
    "    print('{0:7s} '.format(s),end='')\n",
    "    for r in sorted(rxns):\n",
    "        ps = '{0:5d}'.format(nu[r][s]) if s in nu[r].keys() else '     '\n",
    "        print(ps,end='')\n",
    "    print()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 234,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Requirement already satisfied: cvxpy in /Users/jeff/anaconda/lib/python3.5/site-packages\r\n",
      "Requirement already satisfied: ecos>=2 in /Users/jeff/anaconda/lib/python3.5/site-packages (from cvxpy)\r\n",
      "Requirement already satisfied: scs>=1.1.3 in /Users/jeff/anaconda/lib/python3.5/site-packages (from cvxpy)\r\n",
      "Requirement already satisfied: multiprocess in /Users/jeff/anaconda/lib/python3.5/site-packages (from cvxpy)\r\n",
      "Requirement already satisfied: fastcache in /Users/jeff/anaconda/lib/python3.5/site-packages (from cvxpy)\r\n",
      "Requirement already satisfied: six in /Users/jeff/anaconda/lib/python3.5/site-packages (from cvxpy)\r\n",
      "Requirement already satisfied: toolz in /Users/jeff/anaconda/lib/python3.5/site-packages (from cvxpy)\r\n",
      "Requirement already satisfied: numpy>=1.9 in /Users/jeff/anaconda/lib/python3.5/site-packages (from cvxpy)\r\n",
      "Requirement already satisfied: scipy>=0.15 in /Users/jeff/anaconda/lib/python3.5/site-packages (from cvxpy)\r\n",
      "Requirement already satisfied: CVXcanon>=0.0.22 in /Users/jeff/anaconda/lib/python3.5/site-packages (from cvxpy)\r\n",
      "Requirement already satisfied: dill>=0.2.5 in /Users/jeff/anaconda/lib/python3.5/site-packages (from multiprocess->cvxpy)\r\n"
     ]
    }
   ],
   "source": [
    "!pip install cvxpy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 316,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import cvxpy as cvx\n",
    "\n",
    "# extents of reaction\n",
    "X = {r:cvx.Variable() for r in rxns}\n",
    "\n",
    "# net stoichiometric coefficients\n",
    "V = dict()\n",
    "for s in species:\n",
    "    V[s] = sum([nu[r][s]*X[r] for r in rxns if s in nu[r].keys()])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 323,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "obj = cvx.Minimize(V['CO2'])\n",
    "\n",
    "constraints = [\n",
    "    V['Al2O3'] == 0,\n",
    "    V['AlN'] == 0,\n",
    "    V['CH4'] <= 0,\n",
    "    V['H2'] >= 0,\n",
    "    V['CO'] == 0,\n",
    "    V['CH3OH'] == 0,\n",
    "    V['CH3OCH3'] == 3,\n",
    "    V['NH3'] >= 0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 324,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'optimal'"
      ]
     },
     "execution_count": 324,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "prob = cvx.Problem(obj,constraints)\n",
    "prob.solve()\n",
    "prob.status"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 325,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'R1': 2.0, 'R2': 2.0, 'R3': 6.0, 'R4': 6.0, 'R5': 3.0}"
      ]
     },
     "execution_count": 325,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "{r:round(X[r].value,4) for r in rxns}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 326,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'Al2O3': -0.0,\n",
       " 'AlN': 0.0,\n",
       " 'CH3OCH3': 3.0,\n",
       " 'CH3OH': -0.0,\n",
       " 'CH4': -6.0,\n",
       " 'CO': 0.0,\n",
       " 'CO2': -0.0,\n",
       " 'H2': -0.0,\n",
       " 'H2O': -3.0,\n",
       " 'N2': -2.0,\n",
       " 'NH3': 4.0}"
      ]
     },
     "execution_count": 326,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "{s:round(V[s].value,4) for s in species}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 327,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[True, True, True, True, True, True, True, True]"
      ]
     },
     "execution_count": 327,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "[c.value for c in constraints]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [default]",
   "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.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}