Chemical Reactor

We consider a single first-order, irreversible chemical reaction in an isothermal CSTR

\[A \to B.\]

The material balances and the system data are provided in [DAR11] and is given by the nonlinear model

\begin{equation*} \begin{split} c_{A}(k+1)&=c_{A}(k)+h\left(\frac{Q(k)}{V}(c_f^{A}-c_{A}(k))-k_r{c_{A}(k)}\right)\\ c_{B}(k+1)&=c_{B}(k)+h\left(\frac{Q(k)}{V}(c_f^{B}-c_{B}(k))+k_r{c_{B}(k)}\right), \end{split} \end{equation*}

in which \(c_A\geq 0\) and \(c_B\geq 0\) are the molar concentrations of \(A\) and \(B\) respectively, and \(Q\leq 20\) (L/min) is the flow through the reactor. The constants and their meanings are given in table below.

Reactor constants
		value	unit
feed concentration of \(A\)	\(c_f^{A}\)	1	mol/L
feed concentration of \(B\)	\(c_f^{B}\)	0	mol/L
volume of the reactor	\(V_R\)	10	L
rate constant	\(k_r\)	1.2	L/(mol min)
equilibrium	\((c_e^{A},c_e^B,Q_e)\)	\((\frac 1 2, \frac 1 2, 12)\)
initial value	\((c_0^{A},c_0^B)\)	\((0.4, 0.2)\)

To initialize the system dynamics a function that implements \(f(x,u)\), where \(x = (c_{A},c_{B})^T\) and \(u=Q\) has to be defined.

V = 10.
cf_A = 1.
cf_B = 0.
k_r = 1.2

def f(x,u):
    y = nmpyc.array(2)
    y[0] = x[0] + 0.5*((u[0]/V) *(cf_A - x[0]) - k_r*x[0])
    y[1] = x[1] + 0.5*((u[0]/V) *(cf_B - x[1]) + k_r*x[1])
    return y

After that, the nMPyC system object can be set by calling

system = nmpyc.system(f, 2, 1, system_type='discrete')

In the next step, the objective is defined by using the stage cost given by

\begin{align*} \ell (c_{A}(k),c_{B}(k),Q(k))&=\frac 1 2\vert c_A(k)-\frac 1 2\vert^2+\frac 1 2 \vert c_B(k)-\frac 1 2\vert^2+\frac 1 2 \vert Q(k) -12 \vert^2\\ \end{align*}

Since we do not need terminal cost, we can initialize the objective directly using the following implementation.

def l(x,u):
    return 0.5 * (x[0]-0.5)**2 + 0.5 * (x[1]-0.5)**2 + 0.5 * (u[0]-12)**2

objective = nmpyc.objective(l)

In terms of the constraints we assume that

\[\begin{split}0 &\leq x_1(k) & < \infty & \quad & \text{for } k=0,\ldots,N \\ 0 &\leq x_2(k) & < \infty & \quad & \text{for } k=0,\ldots,N \\ 0 &\leq u(k) & \leq 20 & \quad & \text{for } k=0,\ldots,N-1.\end{split}\]

This can be realized in the code as follows:

constraints = nmpyc.constraints()
lbx = nmpyc.zeros(2)
ubu = nmpyc.ones(1)*20
lbu = nmpyc.zeros(1)
constraints.add_bound('lower','state', lbx)
constraints.add_bound('lower','control', lbu)
constraints.add_bound('upper','control', ubu)

Moreover, we consider the equilibrium \((c_e^{A},c_e^B,Q_e)\) as th terminal condition for our optimal control problem, which is implemented as

xeq = nmpyc.array([0.5,0.5])
def he(x):
    return x - xeq
constraints.add_constr('terminal_eq', he)

After all components of the optimal control problem have been implemented, we can now combine them into a model and start the MPC loop. For this Purpose, we define

\[x(0) = (0.4,0.2)^T\]

and set \(N=15\), \(K=100\).

model = nmpyc.model(objective,system,constraints)
x0 = nmpyc.array([0.4,0.2])
res = model.mpc(x0,15,100)

Following the simulation we can visualize the results by calling

res.plot()

which generates the plot bellow.