Temperature Optimization for Parallel Reaction Selectivity

Background

In a batch reactor, one feed species can form both a desired product and an unwanted by-product. Temperature changes both reaction rates, but it need not change them equally: the channel with the larger activation energy becomes relatively more important as temperature rises. Heating can therefore improve conversion while degrading selectivity.

This example chooses a smooth heat-hold-cool trajectory for the irreversible parallel network

\[ A\longrightarrow B,\qquad A\longrightarrow C, \]

where \(B\) is desired. The batch must reach 92% conversion in 90 minutes and return to its charging temperature with zero temperature ramp.

Kinetics and dynamic states

For channel \(i\in\{B,C\}\), the first-order Arrhenius rate constant is

\[ k_i(T)=A_i\exp\left(-\frac{E_i}{R_gT}\right). \]

Let \(c_A,c_B,c_C\) be concentrations. Temperature is represented by a scaled state \(\vartheta\in[0,1]\),

\[ T=T_{\min}+(T_{\max}-T_{\min})\vartheta, \]

and its physical ramp \(r_T=\dot T\) is also a state. The control is ramp acceleration \(a_T=\dot r_T\). The full dynamics are

\[ \begin{aligned} \dot c_A&=-[k_B(T)+k_C(T)]c_A,\\ \dot c_B&=k_B(T)c_A,\\ \dot c_C&=k_C(T)c_A,\\ \dot\vartheta&=\frac{r_T}{T_{\max}-T_{\min}},\\ \dot r_T&=a_T. \end{aligned} \]

The kinetic data are

\[ \begin{array}{lll} A_B=2.0\times10^8\ \mathrm{min^{-1}},& E_B=60\,000\ \mathrm{J/mol},\\ A_C=5.0\times10^{10}\ \mathrm{min^{-1}},& E_C=75\,000\ \mathrm{J/mol},& R_g=8.314\ \mathrm{J/(mol\,K)}. \end{array} \]

Because \(E_C>E_B\), high temperature increases the unwanted-to-desired rate ratio \(k_C/k_B\).

Objective and constraints

The initial and terminal conditions are

\[ \begin{aligned} (c_A,c_B,c_C,\vartheta,r_T)(0) &=(1,0,0,0,0),\\ c_A(90)&=0.08\ \mathrm{mol/L},\\ \vartheta(90)&=0,\qquad r_T(90)=0. \end{aligned} \]

Final \(c_B\) and \(c_C\) are free optimization outcomes. The path bounds are

\[ \begin{aligned} 0&\le c_A,c_B,c_C\le1\ \mathrm{mol/L},\\ 298.15&\le T(t)\le343.15\ \mathrm K,\\ |r_T(t)|&\le2.0\ \mathrm{K/min},\\ |a_T(t)|&\le0.35\ \mathrm{K/min^2}. \end{aligned} \]

The objective is

\[ \min J=c_C(t_f)+w_r\int_0^{t_f} \left[ \left(\frac{r_T}{2.0}\right)^2 +0.15\left(\frac{a_T}{0.35}\right)^2 \right]\,\mathrm dt, \qquad w_r=2.0\times10^{-4}\ \mathrm{\frac{mol}{L\,min}}. \]

At exactly fixed conversion, conservation gives \(c_B(t_f)+c_C(t_f)=0.92\ \mathrm{mol/L}\), so minimizing \(c_C\) alone would maximize desired yield and the ratio \(c_B/c_C\). The implemented objective is not strictly equivalent to maximum selectivity, because it adds a finite temperature-ramp and ramp-acceleration regularizer. It deliberately trades a small amount of product performance for a smoother thermal program.

Variables and units

Symbol	Meaning	Unit
\(t\)	Batch time	min
\(c_A,c_B,c_C\)	Reactant, desired-product, and by-product concentrations	mol/L
\(T,\vartheta\)	Physical and scaled reactor temperature	K, dimensionless
\(r_T\)	Temperature ramp	K/min
\(a_T\)	Temperature-ramp acceleration	K/min²
\(k_B,k_C\)	First-order rate constants	1/min
\(A_B,A_C\)	Arrhenius pre-exponential factors	1/min
\(E_B,E_C\)	Activation energies	J/mol
\(J\)	By-product plus smoothness objective	mol/L

Modeling choices and limits

Temperature ramp is a state and ramp acceleration is the control, so both the thermal trajectory and its first derivative are continuous and bounded. A \(200\times2\) Lobatto mesh gives piecewise-linear acceleration samples and resolves the smooth Arrhenius kinetics.

This is an ideal, spatially uniform constant-volume batch model. It assumes irreversible first-order reactions and omits energy balance, heat-transfer area, jacket or heater dynamics, reagent depletion other than \(A\), secondary reactions, volume change, uncertainty, and safety constraints such as heat release. The optimized \(T(t)\) should be read as a set-point trajectory for an ideal thermal actuator, not a complete reactor-control design.

Run the example

python -m examples.parallel_reaction_selectivity

Save the figure without opening a window:

python -m examples.parallel_reaction_selectivity --save parallel-reaction-selectivity.png --no-show

Independent numerical validation

The script reconstructs concentrations, temperature, ramp, and acceleration at 6,001 times. It independently evaluates the exact reactant relation

\[ c_A(t)=c_A(0)\exp\left[ -\int_0^t(k_B(T(s))+k_C(T(s)))\,\mathrm ds \right], \]

checks \(c_A+c_B+c_C=c_A(0)\), integrates both product formation rates, checks the temperature-ramp balance, and evaluates every endpoint and path bound.

Quantity	Verified value
Objective \(J\)	\(0.38924799\ \mathrm{mol/L}\)
Desired-product yield \(c_B(t_f)/c_A(0)\)	\(0.533074\)
Desired/unwanted ratio \(c_B(t_f)/c_C(t_f)\)	\(1.377718\)
Conversion	\(0.920000\)
Peak temperature	\(320.047293\ \mathrm K\)
Maximum analytical reactant error	\(1.819\times10^{-5}\ \mathrm{mol/L}\)
Maximum species-conservation error	\(2.220\times10^{-16}\ \mathrm{mol/L}\)
Maximum dense product-rate balance error	\(3.319\times10^{-6}\ \mathrm{mol/L}\)
Dense temperature-balance error	\(0\) at printed precision
Maximum dense path-bound violation	\(0\) at printed precision

Source code

See the complete runnable example: examples/parallel_reaction_selectivity.py.