# Maximizing Concentration of an Intermediate in a Batch Reactor

A desired product $B$ is produced as intermediate in a series reaction

\begin{align}
A \overset{k_A}{\longrightarrow} B \overset{k_B}{\longrightarrow} C
\end{align}

where $A$ is a raw material and $C$ is a undesired b
y-product. The reaction operates at temperature where the rate constants are $k_A = 0.5\ \mbox{min}^{-1}$ and $k_A = 0.1\ \mbox{min}^{-1}$. The raw material is available as a solution with concenration $C_{A,f} = 2.0\ \mbox{moles/liter}$.

A 100 liter tank is available which can be operated as a batch reactor. At what point should the operation be terminated in order to maximize the concentration of $B$?


```
!pip install -q pyomo
!wget -N -q "https://ampl.com/dl/open/ipopt/ipopt-linux64.zip"
!unzip -o -q ipopt-linux64

ipopt_executable = '/content/ipopt'
```

## Mathematical Model

A material balance for an isothermal stirred batch reactor with a volume $V = 40$ liters and an initial concentration $C_{A,f}$ is given by

\begin{align}
V\frac{dC_A}{dt} & = - V k_A C_A \\
V\frac{dC_B}{dt} & = V k_A C_A - V k_B C_B
\end{align}

Eliminating the common factor $V$

\begin{align}
\frac{dC_A}{dt} & = - k_A C_A \\
\frac{dC_B}{dt} & = k_A C_A - k_B C_B
\end{align}

With an initial concentration $C_{A,f}$. A numerical solution to these equations is shown in the following cell.


```
%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
from scipy.integrate import odeint

V = 40     # liters
kA = 0.5   # 1/min
kB = 0.1   # l/min
CAf = 2.0  # moles/liter

def batch(X, t):
    CA, CB = X
    dCA_dt = -kA*CA
    dCB_dt = kA*CA - kB*CB
    return [dCA_dt, dCB_dt]

t = np.linspace(0,30,200)
soln = odeint(batch, [CAf,0], t)
plt.plot(t, soln)
plt.xlabel('time / minutes')
plt.ylabel('concentration / moles per liter')
plt.title('Batch Reactor')
plt.legend(['$C_A$','$C_B$'])
```

## Optimization with `scipy.minimize_scalar`

To find the maximum value, we first write a function to compute $C_B$ for any value of time $t$.


```
def CB(tf):
    soln = odeint(batch, [CAf, 0], [0, tf])
    return soln[-1][1]
```

We gain use `minimize_scalar` to find the value of $t$ that minimizes the negative value of $C_B(t)$.|


```
from scipy.optimize import minimize_scalar
minimize_scalar(lambda t: -CB(t), bracket=[0,50])
```




         fun: -1.3374806339222158
        nfev: 23
         nit: 19
     success: True
           x: 4.023594924340666




```
tmax = minimize_scalar(lambda t: -CB(t), bracket=[0,50]).x

print('Concentration c_B has maximum', CB(tmax), 'moles/liter at time', tmax, 'minutes.')
```

    Concentration c_B has maximum 1.3374806339222158 moles/liter at time 4.023594924340666 minutes.


## Solution Using Pyomo

The variable to be found is the time $t_f$ corresponding to the maximum concentration of $B$. For this purpose we introduce a scaled time

$$\tau = \frac{t}{t_f}$$

so that $\tau=1$ as the desired solution. The problem then reads

\begin{align}
\max_{t_f} C_B(\tau=1)
\end{align}

subject to

\begin{align}
\frac{dC_A}{d\tau} & = - t_f k_A C_A \\
\frac{dC_B}{d\tau} & = t_f(k_A C_A - k_B C_B)
\end{align}

The solution to this problem is implemented as a solution to the following Pyomo model.


```
from pyomo.environ import *
from pyomo.dae import *

V   = 40    # liters
kA  = 0.5   # 1/min
kB  = 0.1   # l/min
cAf = 2.0   # moles/liter

m = ConcreteModel()

m.tau = ContinuousSet(bounds=(0, 1))

m.tf = Var(domain=NonNegativeReals)
m.cA = Var(m.tau, domain=NonNegativeReals)
m.cB = Var(m.tau, domain=NonNegativeReals)

m.dcA = DerivativeVar(m.cA)
m.dcB = DerivativeVar(m.cB)

m.odeA = Constraint(m.tau, 
    rule=lambda m, tau: m.dcA[tau] == m.tf*(-kA*m.cA[tau]) if tau > 0 else Constraint.Skip)
m.odeB = Constraint(m.tau,
    rule=lambda m, tau: m.dcB[tau] == m.tf*(kA*m.cA[tau] - kB*m.cB[tau]) if tau > 0 else Constraint.Skip)

m.ic = ConstraintList()
m.ic.add(m.cA[0]  == cAf)
m.ic.add(m.cB[0]  == 0)

m.obj = Objective(expr=m.cB[1], sense=maximize)

TransformationFactory('dae.collocation').apply_to(m)
SolverFactory('ipopt', executable=ipopt_executable).solve(m)
print('Concentration c_B has maximum', m.cB[1](), 'moles/liter at time', m.tf(), 'minutes.')
```

    Concentration c_B has maximum 1.3374805810221073 moles/liter at time 4.023594178375687 minutes.



```

```

1a9c4af4bb247fdb41b4f1668dc81e82665d064f
ipynb
Jupyter Notebook
notebooks/dae/Maximizing_Concentration_of_an_Intermediate_in_a_Batch_Reactor.ipynb
edgBR/ND-Pyomo-Cookbook
0fb121f7e572e088a83c22c166c159ff55c2efb4
Apache-2.0
notebooks/dae/Maximizing_Concentration_of_an_Intermediate_in_a_Batch_Reactor.ipynb
edgBR/ND-Pyomo-Cookbook
0fb121f7e572e088a83c22c166c159ff55c2efb4
Apache-2.0
notebooks/dae/Maximizing_Concentration_of_an_Intermediate_in_a_Batch_Reactor.ipynb
edgBR/ND-Pyomo-Cookbook
0fb121f7e572e088a83c22c166c159ff55c2efb4
Apache-2.0
2020-09-03T18:53:07.000Z
2020-09-03T18:53:07.000Z
Qwen/Qwen-72B
1. YES
2. YES