Chemical Reactor Design#
Prepared by: Zhicheng Lu (zlu3@nd.edu, 2024) and Shammah Lilonfe (slilonfe@nd.edu, 2024) at the University of Notre Dame.
Learning Objectives#
Review syntax for Pyomo and Pyomo.DAE
Demonstrate the power of Pyomo in solving and optimizing complex dynamic process systems
Problem Statement: Design of Chemical Reactors for Chemical Transformations#
Chemical reactors are often the most important unit operations in a chemical plant. Reactors come in many forms, however two of the most common idealizations are the continuously stirred tank reactor (CSTR) and the plug flow reactor. The CSTR is often used in modeling studies, and it can be effectively modeled as a lumped parameter system. In this example, we will consider the following reaction scheme known as the Van de Vusse reaction:
A diagram of the system is shown in figure below, where F is the volumetric flowrate. The reactor is assumed to be filled to a constant volume, and the mixture is assumed to have constant density, so the volumetric flowrate into the reactor is equal to the volumetric flowrate out of the reactor. Since the reactor is assumed to be wellmixed, the concentrations in the reactor are equivalent to the concentrations of each component flowing out of the reactor, given by \(C_A\), \(C_B\), \(C_C\), and \(C_D\).
This example (both text and figure) was taken from Hart et al (2017).
Install Packages and Load Modules#
# Importing the libraries needed for the project
import sys
if "google.colab" in sys.modules:
!wget "https://raw.githubusercontent.com/ndcbe/optimization/main/notebooks/helper.py"
import helper
helper.easy_install()
else:
sys.path.insert(0, '../')
import helper
helper.set_plotting_style()
# Importing packages to be used in this project
import pyomo.environ as pyo
import pyomo.dae as dae
import matplotlib.pyplot as plt
Original Problem from the Reference Material#
The goal is to produce product B from a feed containing reactant A. If we design a reactor that is too small, we will obtain insufficient conversion of A to the desired product B. However, given the reaction scheme, if the reactor is too large (e.g., too much reaction is allowed to occur), a significant amount of the desired product B will be further reacted to form the undesired product C. As a result, the goal in this exercise will be to solve for the optimal reactor volume producing the maximum outlet concentration for product B.
The steady-state mole balances for each of the four components are given by,
The known parameters for the system are,
Since the volumetric flowrate F always appears as the numerator over the reactor volume V, it is common to consider this ratio as a single variable, called the space velocity (sv). Our optimization formulation will seek to find the space-velocity that maximizes the outlet concentration of the desired product B.
This original problem was taken from Hart et al (2017).
Optimization Model Formulation#
Mathematically, we want to solve the following optimal reactor volume problem:
Here “s.t.” means “subject to”, i.e., the constraints.
Degrees of Freedom Analysis#
From the optimization model formulated above:
Number of variables = 5
Number of equality constraints = 4
Therefore, the degrees of freedom = 5 - 4 = 1
Pyomo Implementation of the Original Problem#
# function that creates pyomo model for the reaction system
def create_model(k1, k2, k3, caf):
"""Creates reactor model
Argument:
caf: feed concentration of A, gmol/m^3
k1: rate constant for A to B, min^-1
k2: rate constant for B to C, min^-1
k3: rate constant for 2A to D, m^3/(gmol.min)
Returns:
model: pyomo model"""
# Create a Pyomo model
model = pyo.ConcreteModel()
# create the variables
model.sv = pyo.Var(initialize=1.0, within=pyo.PositiveReals) # space-velocity
model.ca = pyo.Var(initialize=5000.0, within=pyo.PositiveReals) # concentration of A
model.cb = pyo.Var(initialize=2000.0, within=pyo.PositiveReals) # concentration of B
model.cc = pyo.Var(initialize=2000.0, within=pyo.PositiveReals) # concentration of C
model.cd = pyo.Var(initialize=1000.0, within=pyo.PositiveReals) # concentration of D
# create the objective
model.obj = pyo.Objective(expr = model.cb, sense=pyo.maximize)
# create the constraints
model.ca_bal = pyo.Constraint(expr = (0 == model.sv * caf - model.sv * model.ca - k1 * model.ca -
2.0 * k3 * model.ca ** 2.0)) # species A mole balance
model.cb_bal = pyo.Constraint(expr=(0 == -model.sv * model.cb + k1 * model.ca - k2 * model.cb)) # species B mole balance
model.cc_bal = pyo.Constraint(expr=(0 == -model.sv * model.cc + k2 * model.cb)) # species C mole balance
model.cd_bal = pyo.Constraint(expr=(0 == -model.sv * model.cd + k3 * model.ca ** 2.0)) # species D mole balance
return model
# define data
k1 = 5.0/6.0 # min^-1
k2 = 5.0/3.0 # min^-1
k3 = 1.0/6000.0 # m^3/(gmol.min)
caf = 10000.0 # gmol/m^3
# reactor model object
m = create_model(k1, k2, k3, caf)
# solve model
solver = pyo.SolverFactory('ipopt')
solver.solve(m, tee=True)
# print the solution of sv and the outlet concentration of B
print(f"The optimum space-velocity is: {round(pyo.value(m.sv), 2)} min^-1")
print(f"The outlet concentration of species B is: {round(pyo.value(m.cb), 2)} gmol/m^3")
# calculate optimum reactor volume for a volumetric flow rate of 1 m^3/min
F = 1
optimum_reactor_volume = F / pyo.value(m.sv)
print(f"The optimum reactor volume is: {round(optimum_reactor_volume, 2)} m^3")
Ipopt 3.13.2:
******************************************************************************
This program contains Ipopt, a library for large-scale nonlinear optimization.
Ipopt is released as open source code under the Eclipse Public License (EPL).
For more information visit http://projects.coin-or.org/Ipopt
This version of Ipopt was compiled from source code available at
https://github.com/IDAES/Ipopt as part of the Institute for the Design of
Advanced Energy Systems Process Systems Engineering Framework (IDAES PSE
Framework) Copyright (c) 2018-2019. See https://github.com/IDAES/idaes-pse.
This version of Ipopt was compiled using HSL, a collection of Fortran codes
for large-scale scientific computation. All technical papers, sales and
publicity material resulting from use of the HSL codes within IPOPT must
contain the following acknowledgement:
HSL, a collection of Fortran codes for large-scale scientific
computation. See http://www.hsl.rl.ac.uk.
******************************************************************************
This is Ipopt version 3.13.2, running with linear solver ma27.
Number of nonzeros in equality constraint Jacobian...: 11
Number of nonzeros in inequality constraint Jacobian.: 0
Number of nonzeros in Lagrangian Hessian.............: 5
Total number of variables............................: 5
variables with only lower bounds: 5
variables with lower and upper bounds: 0
variables with only upper bounds: 0
Total number of equality constraints.................: 4
Total number of inequality constraints...............: 0
inequality constraints with only lower bounds: 0
inequality constraints with lower and upper bounds: 0
inequality constraints with only upper bounds: 0
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
0 -2.0000000e+03 7.50e+03 6.25e-01 -1.0 0.00e+00 - 0.00e+00 0.00e+00 0
1 -1.0801475e+03 5.54e+02 2.46e+00 -1.0 1.39e+03 - 5.54e-01 1.00e+00h 1
2 -1.0763574e+03 8.86e+01 1.87e+02 -1.0 3.66e+02 - 9.31e-01 1.00e+00h 1
3 -1.0727252e+03 1.07e+01 5.26e+00 -1.0 9.35e+01 - 9.71e-01 1.00e+00h 1
4 -1.0726714e+03 4.01e+00 1.22e-01 -1.0 6.03e+01 - 9.68e-01 1.00e+00h 1
5 -1.0724371e+03 3.44e-04 2.93e-05 -2.5 4.19e-01 - 1.00e+00 1.00e+00h 1
6 -1.0724372e+03 1.63e-08 3.93e-09 -3.8 3.85e-03 - 1.00e+00 1.00e+00h 1
7 -1.0724372e+03 5.46e-11 2.70e-11 -5.7 2.24e-04 - 1.00e+00 1.00e+00h 1
8 -1.0724372e+03 4.55e-13 3.56e-14 -8.6 2.78e-06 - 1.00e+00 1.00e+00h 1
Number of Iterations....: 8
(scaled) (unscaled)
Objective...............: -1.0724372001086319e+03 -1.0724372001086319e+03
Dual infeasibility......: 3.5606114692563649e-14 3.5606114692563649e-14
Constraint violation....: 4.5474735088646414e-14 4.5474735088646412e-13
Complementarity.........: 2.5059065225787338e-09 2.5059065225787338e-09
Overall NLP error.......: 2.5059065225787338e-09 2.5059065225787338e-09
Number of objective function evaluations = 9
Number of objective gradient evaluations = 9
Number of equality constraint evaluations = 9
Number of inequality constraint evaluations = 0
Number of equality constraint Jacobian evaluations = 9
Number of inequality constraint Jacobian evaluations = 0
Number of Lagrangian Hessian evaluations = 8
Total CPU secs in IPOPT (w/o function evaluations) = 0.000
Total CPU secs in NLP function evaluations = 0.000
EXIT: Optimal Solution Found.
The optimum space-velocity is: 1.34 min^-1
The outlet concentration of species B is: 1072.44 gmol/m^3
The optimum reactor volume is: 0.74 m^3
Optimal Solution Found for the Original Problem#
The optimal space-velocity is \(\mathrm{1.34\,min^{-1}}\), giving an outlet concentration for B of \(\mathrm{1072.44\,gmol\cdot m^{-3}}\). For a volumetric flow rate of \(\mathrm{1\,m^{3}\cdot min^{-1}}\), the optimum reactor volume is \(\mathrm{0.74\,m^{3}}\).
Extension: Adding a New Species C and Considering Reaction Kinetics#
In this project, species C is the desired product, and we are considering the dynamic state of the system (most often, industrial operations are dynamic). The goal of this work is to determine the volume of the reactor that will maximize the concentration of species C in the effluent after 10 minutes. Furthermore, for economic reasons, the reactor volume should not exceed \(\mathrm{100\,m^3}\).
The mole balance for each of the species in this system under dynamic state is given by:
Where \(C_A\), \(C_B\), \(C_C\) and \(C_D\) are the concentrations of species A, B, C and D, respectively. \(F\) is the feed flow rate, \(C_{Af}\) is the concentration of species A in the feed, \(k_1\), \(k_2\) and \(k_3\) are the rate constants of reaction 1 (species A to B), reaction 2 (species B to C) and reaction 3 (species A to D), respectively.
The data for the system are: $\(\begin{align*} C_{Af}=10^4\,\mathrm{\frac{mol}{m^3}},\quad{k_1=\frac{5}{6}\,\mathrm{min^{-1}}},\quad{k_2=\frac{5}{3}\,\mathrm{min^{-1}}},\quad{k_3=\frac{1}{6000}\,\mathrm{\frac{m^3}{mol\cdot min}}},\quad{F=1\,\mathrm{\frac{m^3}{min}}} \end{align*}\)$
The data was obtained from Hart et al (2017).
Optimization Model Formulation#
Mathematically, we want to solve the following optimal reactor volume problem:
Degrees of Freedom Analysis#
From the optimization model formulated above:
Number of variables = 5
Number of equality constraints = 4
Therefore, the degrees of freedom = 5 - 4 = 1
Pyomo Implementation#
We can use Pyomo.dae
to automatically formulate the collocation equations (i.e., add constraints that numerically integrate the ODE model).
# function that creates pyomo model for the dynamic reaction system
def create_dynamic_model(F, CAF, k1, k2, k3):
"""Creates reactor model
Argument:
F: feed volumetric flow rate in m^3/min
CAF: feed concentration of A, gmol/m^3
k1: rate constant for A to B, min^-1
k2: rate constant for B to C, min^-1
k3: rate constant for 2A to D, m^3/(gmol.min)
Returns:
model: pyomo model"""
# Create a Pyomo model
model = pyo.ConcreteModel()
# Define time as a continuous set
model.t = dae.ContinuousSet(bounds=(0, 10)) # time from 0 to 10 minutes
# Define variables for concentrations and reactor volume
model.CA = pyo.Var(model.t, within=pyo.NonNegativeReals) # concentration of A
model.CB = pyo.Var(model.t, within=pyo.NonNegativeReals) # concentration of B
model.CC = pyo.Var(model.t, within=pyo.NonNegativeReals) # concentration of C
model.CD = pyo.Var(model.t, within=pyo.NonNegativeReals) # concentration of D
model.V = pyo.Var(within=pyo.NonNegativeReals, bounds=(0, 100), initialize=0) # reactor volume
# Define derivatives of concentrations
model.dCA_dt = dae.DerivativeVar(model.CA, wrt=model.t)
model.dCB_dt = dae.DerivativeVar(model.CB, wrt=model.t)
model.dCC_dt = dae.DerivativeVar(model.CC, wrt=model.t)
model.dCD_dt = dae.DerivativeVar(model.CD, wrt=model.t)
# CA differential equation
def ca_rate_ode(m, t):
return m.dCA_dt[t] == (F/m.V)*CAF - (F/m.V)*m.CA[t] - k1*m.CA[t] - 2*k3*m.CA[t]**2
model.CA_ode = pyo.Constraint(model.t, rule=ca_rate_ode)
# CB differential equation
def cb_rate_ode(m, t):
return m.dCB_dt[t] == -(F/m.V)*m.CB[t] + k1*m.CA[t] - k2*m.CB[t]
model.CB_ode = pyo.Constraint(model.t, rule=cb_rate_ode)
# CC differential equation
def cc_rate_ode(m, t):
return m.dCC_dt[t] == -(F/m.V)*m.CC[t] + k2*m.CB[t]
model.CC_ode = pyo.Constraint(model.t, rule=cc_rate_ode)
# CD differential equation
def cd_rate_ode(m, t):
return m.dCD_dt[t] == -(F/m.V)*m.CD[t] + k3*m.CA[t]**2
model.CD_ode = pyo.Constraint(model.t, rule=cd_rate_ode)
# Initial conditions
def init_conditions(m):
yield m.CA[0] == CAF
yield m.CB[0] == 0
yield m.CC[0] == 0
yield m.CD[0] == 0
model.init_conditions = pyo.ConstraintList(rule=init_conditions)
# Objective: Maximize concentration of C at t = 10 minutes
def objective_function(m):
return m.CC[10]
model.objective = pyo.Objective(rule=objective_function, sense=pyo.maximize)
return model
# Parameters from data
conc_A_feed = 10000 # gmol/m^3
k1 = 5/6 # min^-1
k2 = 5/3 # min^-1
k3 = 1/6000 # m^3/gmol/min
feed_flowrate = 1 # m^3/min
# reactor model object
reactor_model = create_dynamic_model(feed_flowrate, conc_A_feed, k1, k2, k3)
# Discretize the model using collocation
discretizer = pyo.TransformationFactory('dae.collocation')
discretizer.apply_to(reactor_model, nfe=20, ncp=10)
# Solve the model
solver = pyo.SolverFactory('ipopt')
results = solver.solve(reactor_model, tee=True)
# print optimal reactor volume
optimal_reactor_volume = pyo.value(reactor_model.V) # m^3
print(f"Optimal reactor volume is: {round(optimal_reactor_volume, 2)} m^3")
# Extract the concentration of C at t = 10 minutes
optimal_concentration_C = pyo.value(reactor_model.CC[10])*10**(-3) # kmol/m^3
print(f"Maximum Concentration of C at optimal reactor volume is: {round(optimal_concentration_C, 2)} kmol/m^3")
# Extract results for plotting
time = [t for t in reactor_model.t] # minutes
CA_values = [(pyo.value(reactor_model.CA[t])*10**(-3)) for t in reactor_model.t] # kmol/m^3
CB_values = [(pyo.value(reactor_model.CB[t])*10**(-3)) for t in reactor_model.t] # kmol/m^3
CC_values = [(pyo.value(reactor_model.CC[t])*10**(-3)) for t in reactor_model.t] # kmol/m^3
CD_values = [(pyo.value(reactor_model.CD[t])*10**(-3)) for t in reactor_model.t] # kmol/m^3
# Plot the results
plt.figure(figsize=(8, 6))
plt.plot(time, CA_values, label='$\mathrm{C_{A}}$')
plt.plot(time, CB_values, label='$\mathrm{C_{B}}$')
plt.plot(time, CC_values, label='$\mathrm{C_{C}}$')
plt.plot(time, CD_values, label='$\mathrm{C_{D}}$')
plt.xlabel('Time (min)', fontsize=16,fontweight='bold')
plt.ylabel('Concentration ($\mathbf{\dfrac{kmol}{m^3}}$)', fontsize=16,fontweight='bold')
plt.xticks(fontsize=14)
plt.yticks(fontsize=14)
plt.title('Concentration Profiles Over Time', fontsize=16,fontweight='bold')
plt.legend()
plt.show()
Ipopt 3.13.2:
******************************************************************************
This program contains Ipopt, a library for large-scale nonlinear optimization.
Ipopt is released as open source code under the Eclipse Public License (EPL).
For more information visit http://projects.coin-or.org/Ipopt
This version of Ipopt was compiled from source code available at
https://github.com/IDAES/Ipopt as part of the Institute for the Design of
Advanced Energy Systems Process Systems Engineering Framework (IDAES PSE
Framework) Copyright (c) 2018-2019. See https://github.com/IDAES/idaes-pse.
This version of Ipopt was compiled using HSL, a collection of Fortran codes
for large-scale scientific computation. All technical papers, sales and
publicity material resulting from use of the HSL codes within IPOPT must
contain the following acknowledgement:
HSL, a collection of Fortran codes for large-scale scientific
computation. See http://www.hsl.rl.ac.uk.
******************************************************************************
This is Ipopt version 3.13.2, running with linear solver ma27.
Number of nonzeros in equality constraint Jacobian...: 12619
Number of nonzeros in inequality constraint Jacobian.: 0
Number of nonzeros in Lagrangian Hessian.............: 1006
Error in an AMPL evaluation. Run with "halt_on_ampl_error yes" to see details.
Error evaluating Jacobian of equality constraints at user provided starting point.
No scaling factors for equality constraints computed!
Total number of variables............................: 1609
variables with only lower bounds: 804
variables with lower and upper bounds: 1
variables with only upper bounds: 0
Total number of equality constraints.................: 1608
Total number of inequality constraints...............: 0
inequality constraints with only lower bounds: 0
inequality constraints with lower and upper bounds: 0
inequality constraints with only upper bounds: 0
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
0 -9.9999900e-03 1.00e+06 1.56e+00 -1.0 0.00e+00 - 0.00e+00 0.00e+00 0
1 -9.9998726e-03 1.00e+06 8.10e+02 -1.0 1.01e+06 - 9.90e-07 7.52e-05h 1
2r-9.9998726e-03 1.00e+06 9.99e+02 6.0 0.00e+00 - 0.00e+00 3.10e-07R 3
3r-3.2550526e+00 8.36e+05 2.18e+09 6.0 1.09e+09 - 1.95e-04 1.95e-04f 1
4 -3.2756934e+00 7.01e+05 6.67e+06 -1.0 8.47e+05 - 5.03e-04 2.84e-01h 1
5 -3.3599969e+00 4.50e+05 4.28e+05 -1.0 2.51e+05 - 6.14e-02 9.85e-01h 1
6 -3.4374817e+00 2.50e+05 2.67e+05 -1.0 4.42e+05 - 2.83e-03 4.99e-01h 1
7 -3.6563863e+00 1.27e+05 2.54e+05 -1.0 2.24e+05 - 7.45e-02 9.90e-01h 1
8 -4.0376969e+00 6.23e+04 2.60e+05 -1.0 5.04e+04 - 4.10e-02 9.90e-01h 1
9 -4.2238169e+00 4.57e+04 2.56e+05 -1.0 1.87e+04 - 1.46e-02 3.26e-01h 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
10 -4.8943726e+00 2.20e+04 2.14e+05 -1.0 1.48e+04 - 1.98e-02 9.98e-01h 1
11 -5.0999886e+00 1.69e+04 2.38e+06 -1.0 4.00e+03 - 1.90e-02 2.92e-01h 1
12 -5.8592444e+00 8.22e+03 8.43e+06 -1.0 3.31e+03 - 3.47e-02 1.00e+00h 1
13 -6.6724932e+00 7.12e+03 8.23e+06 -1.0 1.71e+03 - 3.91e-02 1.54e-01h 1
14 -2.5003886e+01 3.51e+03 8.41e+06 -1.0 1.55e+03 - 8.22e-02 1.00e+00h 1
15 -8.4301553e+01 3.08e+03 8.00e+06 -1.0 1.70e+03 - 4.83e-02 1.40e-01h 1
16 -1.2308262e+02 2.89e+03 7.06e+06 -1.0 1.70e+03 - 1.19e-01 6.44e-02h 1
17 -8.1386603e+02 1.44e+03 5.95e+06 -1.0 1.70e+03 - 1.57e-01 1.00e+00f 1
18 -2.0541698e+03 8.40e+02 3.71e+06 -1.0 1.77e+03 - 3.77e-01 7.14e-01f 1
19 -3.5049209e+03 4.20e+02 1.79e+06 -1.0 1.45e+03 - 5.17e-01 1.00e+00f 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
20 -3.7251819e+03 2.85e+02 2.82e+04 -1.0 4.62e+02 - 9.84e-01 4.76e-01f 1
21 -4.0002085e+03 1.42e+02 2.93e+03 -1.0 2.75e+02 - 8.98e-01 1.00e+00f 1
22 -4.0087092e+03 1.01e+02 2.07e+03 -1.0 7.06e+01 - 1.00e+00 4.15e-01h 1
23 -4.0087453e+03 1.00e+02 2.09e+03 -1.0 1.00e+02 - 7.72e-01 2.93e-03h 1
24 -4.0087456e+03 1.00e+02 5.80e+03 -1.0 1.00e+02 - 4.72e-03 2.92e-05h 1
25 -4.0087598e+03 1.00e+02 2.12e+03 -1.0 3.24e+02 - 4.40e-07 4.37e-05f 1
26 -4.3343427e+03 4.70e+00 1.83e+05 -1.0 3.26e+02 - 4.36e-05 1.00e+00f 1
27 -4.3141499e+03 8.39e-03 9.09e+01 -1.0 2.02e+01 - 1.00e+00 1.00e+00h 1
28 -4.3141284e+03 2.41e-08 1.08e+05 -1.7 2.15e-02 - 9.87e-01 1.00e+00h 1
29 -4.3141354e+03 7.14e-08 2.90e+04 -1.7 7.02e-03 - 7.31e-01 1.00e+00f 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
30 -4.3150872e+03 1.32e-03 2.82e+04 -1.7 9.52e-01 - 2.70e-02 1.00e+00f 1
31 -4.5157522e+03 2.76e+02 2.80e+04 -1.7 2.01e+02 -4.0 7.26e-03 1.00e+00f 1
32 -5.3180445e+03 1.11e+02 1.95e+04 -1.7 8.02e+02 - 3.05e-01 1.00e+00f 1
33 -5.0535790e+03 1.43e+00 1.26e+01 -1.7 2.64e+02 -4.5 1.00e+00 1.00e+00h 1
34 -5.0655151e+03 9.00e+00 2.23e+00 -1.7 6.01e+01 - 1.00e+00 1.00e+00h 1
35 -5.0601188e+03 1.92e-01 7.30e-02 -1.7 1.58e+01 - 1.00e+00 1.00e+00h 1
36 -5.0598090e+03 5.39e-03 4.07e-04 -2.5 1.48e+00 - 1.00e+00 1.00e+00h 1
37 -5.0598054e+03 1.71e-04 1.11e-06 -3.8 2.78e-01 - 1.00e+00 1.00e+00h 1
38 -5.0598053e+03 5.24e-07 2.29e-09 -5.7 1.55e-02 - 1.00e+00 1.00e+00h 1
39 -5.0598053e+03 1.94e-10 1.67e-13 -8.6 1.91e-04 - 1.00e+00 1.00e+00h 1
Number of Iterations....: 39
(scaled) (unscaled)
Objective...............: -5.0598052797570590e+03 -5.0598052797570590e+03
Dual infeasibility......: 1.6745249521991457e-13 1.6745249521991457e-13
Constraint violation....: 1.9372237147763371e-10 1.9372237147763371e-10
Complementarity.........: 2.5063723623922813e-09 2.5063723623922813e-09
Overall NLP error.......: 2.5063723623922813e-09 2.5063723623922813e-09
Number of objective function evaluations = 43
Number of objective gradient evaluations = 40
Number of equality constraint evaluations = 43
Number of inequality constraint evaluations = 0
Number of equality constraint Jacobian evaluations = 41
Number of inequality constraint Jacobian evaluations = 0
Number of Lagrangian Hessian evaluations = 39
Total CPU secs in IPOPT (w/o function evaluations) = 0.351
Total CPU secs in NLP function evaluations = 0.008
EXIT: Optimal Solution Found.
Optimal reactor volume is: 11.35 m^3
Maximum Concentration of C at optimal reactor volume is: 5.06 kmol/m^3
Pyomo Enables Dynamic Chemical Reactor Design#
From the concentration profile above, \(C_{A}\) decreases with time. This makes sense beacause species A is the reactant and it is continuously been consumed to form the products.
Also, it is observed that \(C_{B}\) and \(C_{D}\) increase at the early phase of the reaction, but then start to decrease at a later time. This is because at the early phase of the reaction, only species B and D are produced until a later time when species C begins to form, and since the rate constant of species C (\(k_2\)) is faster than the rate constants of species B and D (\(k_1\) and \(k_3\) respectively), \(C_{C}\) continue to increase whereas \(C_{B}\) and \(C_{D}\) deplete.
The optimization result shows that the optimal reactor volume required to maximize the desired product (species C) is \(\mathrm{11.35\,m^{3}}\) and the corresponding species C concentration is \(\mathrm{5.06\,kmol\cdot m^{-3}}\).
Take Away Messages#
Pyomo is a powerful tool in solving complex dynamic optimization problems
Proper modelling of the problem is key to obtaining quality results
Reference#
Hart, W. E., Laird, C. D., Watson, J. P., Woodruff, D. L., Hackebeil, G. A., Nicholson, B. L., & Siirola, J. D. (2017). Pyomo-optimization modeling in python (Vol. 67, p. 277). Berlin: Springer.