APPLICATION OF THE SENSITIVITY ALGORITHM IN THE OPTIMIZATION OF CHEMICAL REACTORS

A. I. Ruban

(Presented by the scientific seminar on automatic control of the Department of PMAKhP)

Problems of analytically calculating optimal conditions for the course of processes in chemical reactors can be posed and solved only when an adequate mathematical description of the reactors is available. Therefore we assume that such a description exists. For example, in the steady-state regime, ideal-displacement reactors are described by a system of ordinary differential equations (ODEs) with respect to the concentrations of the reagents. In this type of reactor, the principal variable with respect to which optimization is carried out is the temperature profile [1—3].

For definiteness, in what follows we shall consider only the indicated type of reactor, although the main results presented are set forth in general form and may be applied to a broad class of chemical-engineering objects (heat exchangers, rectification and absorption columns, etc.).

1. Statement of the Problem

We assume that the process under investigation is described by a system of ODEs

\[ \frac{d\bar{x}}{dt}=\bar{f}(\bar{x}(t), \bar{u}(t)), \tag{1,1} \]

where \(\bar{x}(t)\) is an \(n\)-dimensional column vector of phase coordinates;
\(\bar{u}(t)\) is an \(m\)-dimensional column vector of control actions belonging to the region \(U\). For simplicity, suppose that the beginning and the end of the variation of the variable \(t\) are fixed \((t \in [0, t_k])\), and that the initial values of the phase coordinates are given,

\[ \bar{x}(0)=\bar{x}_0. \tag{1,2} \]

It is necessary to choose the vector of control actions in such a way that the first component of the vector \(\bar{x}\) at the point \(t=t_k\) assumes the minimum value

\[ x_1(t_k)=\min_{\bar{u}(t)\in U}. \tag{1,3} \]

If, in addition to the initial conditions (1,2), final values are also specified for certain components of the vector \(\bar{x}\), then in the corresponding adjoint equations (1,5) below, boundary conditions will not be prescribed for the corresponding \(\bar{\lambda}\).

In accordance with the maximum principle [1–3], the optimal control actions \(\bar u(t)\) are found from the criterion

\[ H=\bar\lambda^{T}(t)\bar f(\bar x(t),\bar u(t))= \max_{\bar u(t)\in U}. \tag{1.4} \]

Here \(T\) denotes matrix transposition, and the variable \(\bar\lambda(t)\), having the dimension of the vector \(x(t)\), satisfies the adjoint system

\[ \frac{d\bar\lambda}{dt} =-\left(\frac{\partial \bar f}{\partial \bar x}\right)^{T}\bar\lambda =-\left(\frac{\partial H}{\partial \bar x}\right)^{T},\qquad \bar\lambda(t_k)= \begin{pmatrix} -1\\ 0\\ \vdots\\ 0 \end{pmatrix}, \tag{1.5} \]

which is solved jointly with the DE (1.1) under conditions (1.2).

If \(\bar u\) lies inside the domain \(U\), then from the necessary condition for the maximum of \(H\)

\[ \frac{dH}{d\bar u}=0 \tag{1.6} \]

one can sometimes obtain an explicit dependence of the optimal controls on the variables \(\bar x\) and \(\bar\lambda\)

\[ \bar u^{*}(t)=\bar u^{*}(\bar x,\bar\lambda). \tag{1.7} \]

Taking (1.7) into account and solving the boundary-value problem (1.1), (1.2), (1.5), we obtain the desired optimal control actions \(\bar u^{*}(t)\).

Thus, through the use of the maximum principle, the original optimization problem has been reduced to solving a system of nonlinear DEs of order \(2n\) with conditions specified at the beginning \((t=0)\) and at the end of the trajectory \((t=t_k)\), i.e., to the so-called nonlinear two-point boundary-value problem. Obtaining reliable numerical methods for solving this problem is one of the principal problems in control theory. Among all existing methods, the most universal and sufficiently effective is the quasilinearization method [4], successfully developed abroad by R. Bellman and his collaborators. In the present paper, for solving the two-point boundary-value problem another method [5] is used, the linearization (sensitivity) method. It is computationally simpler than the quasilinearization method and, in a number of cases, provides a faster rate of convergence.

2. Method for solving the boundary-value problem

We write the system of DEs (1.1), (1.2), (1.5) in compact form

\[ \frac{d\bar z(t)}{dt}=\bar g(\bar z), \tag{2.1} \]

where \(\bar z(t)\) is an \(N\)-dimensional column vector of variables representing \(\bar x(t)\) and \(\bar\lambda(t)\), with linear boundary conditions

\[ (\bar z(0),\bar a_i)=b_i,\quad i=1,2,\ldots,r, \]

\[ (\bar z(t_k),\bar a_i)=b_i,\quad i=r+1,\ldots,N. \tag{2.2} \]

Here \((\bar z,\bar a)\) is the scalar product of vectors, and \(\bar a_i\), \(b_i\) are known constants. It is required to obtain the solution of system (2.1), \(\bar z(t)\).

If, for the DE (2.1), the values of the vector \(\bar z\) were specified either at the initial point \(t=0\) or at the final point \(t=t_k\), then this would be

Cauchy problem, and it could be solved by any standard numerical method on a digital computer. The use of the linearization method precisely makes it possible to pass from conditions (2.2) to initial conditions for \(\bar z(t)\), i.e., it makes it possible to reduce the boundary-value problem (2.1), (2.2) to a Cauchy problem.

Since equations (2.1) are nonlinear, the solution procedure will naturally be iterative and reduces to the following. At the \(l+1\)-st iteration we jointly solve the system of differential equations with known initial conditions

\[ \frac{d\bar z^{\,l}(t)}{dt}=\bar g(\bar z^{\,l}(t)),\qquad \bar z^{\,l}(0)=\bar z_0^{\,l}, \]

\[ \frac{d\bar Z^{\,l}(t)}{dt}=J(\bar z^{\,l})\bar Z^{\,l}(t),\qquad \bar Z^{\,l}(0)=I, \tag{2.3} \]

where

\[ J(\bar z^{\,l})=\left(\frac{\partial g_i}{\partial z_j}\right) \text{ is the Jacobi matrix,} \tag{2.4} \]

\(I\) is the identity matrix,

\(\bar Z^{\,l}(t)\) is the matrix of sensitivity functions of size \(N\times N\),

and we compute the \((l+1)\)-st approximation of the initial conditions \(\bar z^{\,l+1}(0)\) from the system of linear inhomogeneous algebraic equations

\[ (\bar z^{\,l+1}(0),\bar a_i)=b_i,\qquad i=1,\ldots,r, \]

\[ (\bar Z^{\,l}(t_k)\bar z^{\,l+1}(0),\bar a_i) = b_i-(\bar z^{\,l}(t_k)-\bar Z^{\,l}(t_k)\bar z^{\,l}(0),\bar a_i), \qquad i=r+1,\ldots,N. \tag{2.5} \]

The success of applying the described iterative procedure depends to a considerable extent on the success of choosing the initial approximation of the initial conditions for system (2.1), \(\bar z^{\,0}(0)\). It can be shown that if \(\bar z^{\,0}(0)\) is chosen in a neighborhood of the unknown true initial conditions \(\bar z(0)\), then the sequence \(\bar z^{\,0}(0),\bar z^{\,1}(0),\ldots,\bar z^{\,l}(0),\ldots\), obtained from equations (2.3), (2.5), converges to \(\bar z(0)\) and has quadratic convergence

\[ \|\bar z^{\,l+1}(0)-\bar z(0)\| \le \delta \|\bar z^{\,l}(0)-\bar z(0)\|^2, \tag{2.6} \]

where \(\|\bar z^{\,l}(0)-\bar z(0)\|\) denotes the norm of the deviation of the vector \(\bar z^{\,l}(0)\) from the vector \(\bar z(0)\).

3. Example. Consider an ideal displacement chemical reactor in which a reaction of the type

\[ \begin{array}{c} & P\\[-2mm] & \nearrow\\[-1mm] A & \\[-1mm] & \searrow\\[-2mm] & D \end{array} \tag{3.1} \]

takes place, where \(P\) is the useful reaction product and \(D\) is the by-product. We shall denote by \(x_1(t)\) and \(x_2(t)\), respectively, the concentrations of the products \(P\) and \(A\). Here \(t\) is the length of the vector, which varies from 0 to \(t_k\). We also assume that the concentration of product \(A\) is specified at the beginning and at the end of the reactor, i.e., \(x_2(0)=x_2^0\), \(x_2(t_k)=x_2^k\), and the concentration of product \(P\) at the beginning of the reactor, i.e., \(x_1(0)=x_1^0\).

It is necessary to choose the temperature profile of the reactor so that the yield of the useful product is maximal:

\[ x_1(t_k)=\max_{T(t)}. \tag{3.2} \]

The mathematical description of the reactor under consideration is the system of DEs

\[ \frac{d x_1}{dt}=\kappa_1 x_2,\quad x_1(0)=x_1^0, \tag{3.3} \]

\[ \frac{d x_2}{dt}=-(\kappa_1+\kappa_2)x_2,\quad x_2(0)=x_2^0,\quad x_2(t_k)=x_2^k, \]

where the rate constants \(\kappa_1,\ \kappa_2\) are related to the reactor temperature at each point by the Arrhenius law

\[ \kappa_1=\kappa_{10}\exp\left\{-\frac{E_1}{RT(t)}\right\},\quad \kappa_2=\kappa_{20}\exp\left\{-\frac{E_2}{RT(t)}\right\}. \tag{3.4} \]

In the constants \(\kappa_{10},\ \kappa_{20}\) we have included the volumetric flow rate of the reactants and the cross section of the reactor.

In accordance with the maximum principle, let us write the function \(H\) (1.4)

\[ H=\lambda_1\kappa_1 x_2-\lambda_2(\kappa_1+\kappa_2)x_2 \tag{3.5} \]

and construct the system of DEs conjugate to the original one (3.3):

\[ \frac{d\lambda_1}{dt}=-\frac{\partial H}{\partial x_1}=0,\quad \lambda_1(t_k)=1, \]

\[ \frac{d\lambda_2}{dt}=-\frac{\partial H}{\partial x_2}=-\lambda_1\kappa_1+\lambda_2(\kappa_1+\kappa_2). \tag{3.6} \]

In the latter system (3.6) there is no boundary condition at the point \(t=t_k\) for \(\lambda_2\), since for the corresponding coordinate \(x_2\) the condition at the point \(t=t_k\) has already been specified in (3.3). \(\lambda_1(t_k)=1\), since the maximum of the product \(x_1(t_k)\) is sought (in contrast to (1.3)). From the first equation of system (3.6) we obtain that

\[ \lambda_1(t)=1. \tag{3.7} \]

Taking this result into account,

\[ H=[\kappa_1(1-\lambda_2)-\kappa_2\lambda_2]x_2. \tag{3.8} \]

We note that (3.8) contains only the coordinates \(x_2(t)\) and \(\lambda_2(t)\), and \(x_2(t)\) does not depend on \(x_1(t)\) (3.3). Therefore, for simplicity the subscript 2 may be omitted and the system of equations considered as

\[ \frac{dx}{dt}=-(\kappa_1+\kappa_2)x,\quad x(0)=x^0,\quad x(t_k)=x^k, \]

\[ \frac{d\lambda}{dt}=-\kappa_1+\lambda(\kappa_1+\kappa_2), \tag{3.9} \]

associated with the common temperature profile.

For the time being we shall assume that no restrictions are imposed on the reactor temperature. Then one can use the necessary condition for an extremum of \(H\)

\[ \frac{dH}{dT}=\frac{1}{RT^2}\,[E_1\kappa_1(1-\lambda)-E_2\kappa_2\lambda]x=0. \tag{3.10} \]

Hence we find the explicit dependence of the optimal temperature \(T^*\) on \(x(t)\) and \(\lambda(t)\)

\[ T^*(t)=\frac{E_2-E_1}{R\ln\frac{\kappa_{20}E_2}{\kappa_{10}E_1}\frac{\lambda(t)}{1-\lambda(t)}}. \tag{3.11} \]

Substituting this dependence into (3.9), we obtain a boundary-value problem of the type (2.1), (2.2)

\[ \frac{dx}{dt} = -\kappa_1(t)\left(\frac{E_1(1-\lambda)}{E_2\lambda}+1\right)x, \quad x(0)=x^0, \tag{3.12} \]

\[ \frac{d\lambda}{dt} = \kappa_1(t)\left(-\frac{E_1(1-\lambda)^2}{E_2\lambda}+\lambda\right), \quad x(t_k)=x^k, \]

where

\[ \kappa_1(t) = \kappa_{10} \left( \frac{\kappa_{20}E_2\lambda(t)} {\kappa_{10}E_1(1-\lambda(t))} \right)^{\frac{E_1}{E_1-E_2}} . \tag{3.13} \]

We solve the boundary-value problem by the linearization (sensitivity) method (see § 2). For this purpose, at the \(l+1\)-st iteration, by any numerical method we find the solution of the Cauchy problem

\[ \frac{dx}{dt} = -\kappa_1(t)\left(\frac{E_1(1-\lambda)}{E_2\lambda}+1\right)x, \quad x(0)=x^0; \]

\[ \frac{d\lambda}{dt} = -\kappa_1(t)\left(\frac{E_1(1-\lambda)^2}{E_2\lambda}-\lambda\right), \quad \lambda(0)=\lambda^l(0), \tag{3.14} \]

\[ \frac{da}{dt} = -\kappa_1(t)\left[ \frac{E_1 b x_2}{(E_1-E_2)\lambda^2(1-\lambda)} + \left(\frac{E_1(1-\lambda)}{E_2\lambda}+1\right)a \right], \quad a(0)=0, \]

\[ \frac{db}{dt} = \kappa_1(t)\left[ \frac{E_1(-1+2\lambda)}{(E_1-E_2)\lambda^2(1-\lambda)} + \left(\frac{E_1(1-\lambda)}{E_2\lambda}+1\right) \right]b, \quad b(0)=1, \]

and the next approximation for the missing initial conditions \(\lambda(0)\) in (3.12) is computed from the equation

\[ \lambda^{l+1}(0) = \frac{x^k-x(t_k)+a(t_k)\lambda^l(0)}{a(t_k)} . \tag{3.15} \]

Let us also consider a simpler method for solving the above optimization problem. Using the necessary conditions for an extremum of the function \(H\) (3.8), we obtain

\[ \lambda = \frac{E_1\kappa_1}{E_1\kappa_1+E_2\kappa_2}. \tag{3.16} \]

Substituting this equation into \(H\), we obtain the maximum value of the function \(H=H^*\),

\[ H^* = x(E_2-E_1)\frac{\kappa_1\kappa_2}{E_1\kappa_1+E_2\kappa_2}. \tag{3.17} \]

Since it is known [3] that \(H\geq 0\), the optimal temperature profile can be found only if

\[ E_2 \geq E_1. \tag{3.18} \]

Moreover, it is known that on the optimal trajectory the function does not depend on the variable \(t\), that is,

\[ \frac{dH^*}{dt}=0. \]

From this condition, taking (3.9) into account, it is not difficult to obtain the DE for the optimal temperature.

\[ \frac{dT}{dt}=\frac{RT^2}{E_1E_2}(E_1\kappa_1+E_2\kappa_2). \tag{3.19} \]

In this equation the initial value of the temperature \(T(0)\) is unknown; and if DE (3.19) is considered together with the first equation of system (3.9), then we again obtain a boundary-value problem, whose solution at the \((l+1)\)-st iteration reduces to the joint solution of the system

\[ \frac{dx}{dt}=-(\kappa_1+\kappa_2)x,\quad x(0)=x^0, \]

\[ \frac{dT}{dt}=\frac{RT^2}{E_1E_2}(E_1\kappa_1+E_2\kappa_2);\quad T(0)=T^l(0), \tag{3.20} \]

\[ \frac{da}{dt}=-(\kappa_1+\kappa_2)a-\frac{1}{RT^2}(\kappa_1E_1+\kappa_2E_2)bx,\quad a(0)=0, \]

\[ \frac{db}{dt}=\left[2RT(E_1\kappa_1+E_2\kappa_2)+(E_1^2\kappa_1+E_2^2\kappa_2)\right]\frac{b}{E_1E_2},\quad b(0)=1 \]

with subsequent computation of the \((l+1)\)-st approximation for \(T(0)\)

\[ T^{l+1}(0)=T^l(0)+\frac{x^*-x(t_k)}{a(t_k)}. \tag{3.21} \]

The latter algorithm was tested on a digital computer and showed a high rate of convergence.

REFERENCES

1. Yu. M. Volin, G. M. Ostrovskii, M. G. Slin'ko. Application of the maximum principle for determining the optimal regime of exothermic processes. Kinetics and Catalysis, vol. 4, no. 5, Moscow, 1963, (760–767).

2. G. M. Ostrovskii, Yu. M. Volin. On one method for calculating optimal systems. — Izv. AN SSSR. Technical Cybernetics, No. 2, 1968, (174–185).

3. A. I. Boyarinov, V. V. Kafarov. Optimization Methods in Chemical Technology. “Khimiya,” Moscow, 1969.

4. R. Bellman, R. Kalaba. Quasilinearization and Nonlinear Boundary-Value Problems. “Mir,” Moscow, 1968.

5. A. I. Ruban. Some questions in the mathematical description of dynamic objects (dissertation). Tomsk State University, 1969.