On-line dynamic optimization and control strategy for improving the performance of batch reactors

Chemical Engineering and Processing 44 (2005) 101–114

On-line dynamic optimization and control strategy for improving the performance of batch reactors A. Arpornwichanop a , P. Kittisupakorn a , I.M. Mujtaba b,∗ b

a Department of Chemical Engineering, Chulalongkorn University, Bangkok 10330, Thailand School of Engineering, Design & Technology (EDT 3), University of Bradford, Bradford, West Yorkshire BD7 1DP, UK

Received 23 December 2002; received in revised form 25 August 2003; accepted 28 April 2004 Available online 2 July 2004

Abstract Since batch reactors are generally applied to produce a wide variety of specialty products, there is a great deal of interest to enhance batch operation to achieve high quality and purity product while minimizing the conversion of undesired by-product. The use of process optimization in the control of batch reactors presents a useful tool for operating batch reactors efficiently and optimally. In this work, we develop an approach, based on an on-line dynamic optimization strategy, to modify optimal temperature set point profile for batch reactors. Two different optimization problems concerning batch operation: maximization of product concentration and minimization of batch time, are formulated and solved using a sequential optimization approach. An Extended Kalman Filter (EKF) is incorporated into the proposed approach in order to update current states from their delayed measurement and to estimate unmeasurable state variables. A nonlinear model-based controller: generic model control algorithm (GMC) is applied to drive the temperature of the batch reactor to follow the desired profile. A batch reactor with complex exothermic reaction scheme is used to demonstrate the effectiveness of the proposed approach. The simulation results indicate that with the proposed strategy, large improvement in batch reactor performance, in term of the amount of a desired product and batch operation time, can be achieved compared to the method where the optimal temperature set point is pre-determined. This paper gathers and integrates efficient well-known methodologies such as EKF, GMC and on-line optimization together, resulting in an applicable, reliable and robust control technique for batch reactors. © 2004 Elsevier B.V. All rights reserved.

1. Introduction In many chemical industries, there is an increasing trend to place a consideration on the production of high value products (e.g. polymers, pharmaceuticals, and specialty chemicals) in batch processes. As an important main unit in such processes, a batch reactor is generally involved in manufacturing of these products. The use of batch reactors offers many advantages. Firstly, a batch reactor is quite flexible, it can adapt to small volume production of various products, which are greatly submitted to the rapid changes in market conditions and the advent of new technology. Secondly, a batch reactor provides the natural way to scale-up processes from laboratory experiment where complex chemical synthesis is studied, to industrial manufacturing. Finally, it is ∗ Corresponding author. Tel.: +44 1274 233645; fax: +44 1274 235700. E-mail address: [email protected] (I.M. Mujtaba).

0255-2701/$ – see front matter © 2004 Elsevier B.V. All rights reserved. doi:10.1016/j.cep.2004.04.010

also especially suitable to carry out reactions where materials involved are dangerous and difficult to handle [1,2]. As batch reactors are used to produce a wide variety of expensive products, it is known that this process involves several competing reactions which may cause undesired product or waste. As a result, there is a great deal of interest to enhance batch operation to achieve high quality and purity products while minimizing the conversion of undesired by-products. Recently, the use of process optimization in the control of batch reactors has received much attention in the literature. This provides a useful tool for operating batch reactors efficiently and optimally. For this purpose, it is desirable to optimize the process conditions during its operation in order to meet the desired product and safety specifications. A control system is an essential part to ensure that the desired operating conditions can be maintained as close as possible during the course of a batch run. However, achieving such a proposed method for an optimal operation of batch reactors is quite difficult, and still

102

A. Arpornwichanop et al. / Chemical Engineering and Processing 44 (2005) 101–114

provides challenging and interesting problems. This is mainly due to the inherent complexity of the batch reactors which can be characterized by (i) highly nonlinear behavior resulting from the dependence of reaction rates on concentrations and temperature, (ii) time-varying system; the process variables (e.g. concentration, temperature) and parameters change with time, (iii) no steady state operating condition, (iv) imperfect model; complex kinetic reactions occurring within batch reactors are rarely well understood that leads to an inaccuracy in developing the system model, and (v) lack of measurement information; the product qualities or key properties to be controlled (e.g. molecular weight) cannot be measured until the end of batch run or even if they can be measured sometime (e.g. concentration), there is a significant time delay. Only a few of physical quantities such as temperature and pressure are available for direct on-line measurement. That makes direct control of product properties difficult [2]. Although, in recent years there have been significant advances in developing new sensors for measuring these product properties, they have rarely been used in industrial processes due to high operating cost and expensive investment on the measurement devices. Thus, the usual practice is to control other variables that can be measured rapidly in order to obtain desired product properties instead. In general, optimal batch reactor operation can be carried out by two-step approach; firstly, determining an optimal set point profile of key operating process variables such as temperature [3]and secondly, tracking the desired profile by a control system [4]. The optimal profile can usually be determined off-line by solving an optimal control problem. However, as mentioned above, because of the complexity of chemical reaction schemes, modeling errors are always present and in addition, process disturbances can occur during the process operation. Due to the existence of this error and disturbance, the final product may significantly differ from the desired value, even though the pre-specified optimal profile is tracked perfectly [5]. To realize this fact, it is necessary to recalculate the optimal profile as an on-line optimization strategy whenever new feedback information is available. This strategy could compensate the modeling error leading to process operation improvement. In order to perform the on-line optimization strategy, the knowledge of current state variables and/or parameters in the process models is required. Due to the fact that some of these variables cannot be known exactly or sometime can be measured with time delay, it is essential to include an on-line estimator to estimate these process variables using available process measurements as well. The sequence of an estimation and optimization procedure is known as an estimation-optimization task [6]. As in several estimation techniques, an Extended Kalman Filter (EKF) has become increasingly popular because it is relatively easy to implement. It has been found that the EKF can be applied to a number of chemical process applications with great success. Once the estimate of unknown process variables is deter-

mined and then the model is updated, the optimization is performed on-line to generate a new optimal input profile. With the modified optimal profile, a designed controller is used to control the system to follow this profile until the new one is available. Apart from specifying the optimal set point profile, a control system used to implement such a profile is another important issue to be considered. This is because the deviation from the desired profile may cause an off-spec product. However, since it is well known that the control of batch reactors is difficult due to the inherently nonlinear behavior, the use of a linear control technique may give a poor performance. For this reason, many advanced control techniques have been developed and applied to the control of batch reactor. These include, for example, nonlinear feedforward-feedback control [7], generic model control [8], adaptive control [1], globally linearizing control [9], dynamic matrix control [10], linear model predictive control [11], or inverse model control [12]. Reviews on the progress in control techniques that have been applied to the control of batch reactors as well as their importance and performance are given by [13]. Among these advanced control methods, a generic model control (GMC) technique is one of the most studied control algorithm. This is because nonlinear process models can be interpreted straightaway in the GMC control algorithm so that they do not need to be linearized. Furthermore, its implementation is relatively easy when compared to other model-based control methods; consequently, the application of this control technique appears in many chemical processes. In this work, an approach based on an on-line dynamic optimization strategy to modify optimal temperature set point profile is developed for improving batch operation performance. To demonstrate the effectiveness of the developed approach, the batch reactor studied by [8], where two parallel exothermic reactions occur, is chosen here as a case study. For solving the on-line optimization problem, it needs the knowledge of the current states of the system. Although most physical quantities (e.g. temperature, flowrate) can be measured frequently and are available as on line measurement, some other properties (e.g. concentration) are measured infrequently with time delay. To overcome this difficulty, the Extended Kalman Filter (EKF) is incorporated into the proposed strategy in order to estimate the concentrations from their delayed measurements. The resulting optimal control problem is solved by the sequential solution and optimization method which is often referred to as the control vector parameterization (CPV) method. Once the optimal temperature profile is modified, a controller based on generic model control algorithm (GMC) is applied to control the batch reactor temperature following the desired profile. In the GMC formulation, the EKF is also used to estimate the heat released from the chemical reactions. It should be noted that the proposed strategy would be one of several strategy studied to promote the applicability of on-line dynamic optimization with set point tracking for improving a batch process; the aim of this work is to integrate efficient methods

A. Arpornwichanop et al. / Chemical Engineering and Processing 44 (2005) 101–114

including the EKF and GMC and on-line optimization leading to an applicable, reliable and successful strategy for an actual implementation.

with


k11

k1 = exp

k12 − Tr + 273.15




k22 k2 = exp k21 − Tr + 273.15

2. Batch reactor A reactor system considered by [8] which consists of a batch reactor and jacket cooling system is chosen here as a case study. The typical diagram of this system is shown in Fig. 1. It is assumed that two parallel highly exothermic reactions are carried out in the reactor k1

A + B− →C

103





W = MWA MA + MWB MB + MWC MC + MWD MD M r = MA + M B + M C + M D Cpr =

(CpA MA + CpB MB + CpC MC + CpD MD ) Mr

Qr = −H1 (k1 MA MB ) − H2 (k2 MA MC )

k2

→D A + C−

Qj = UA(Tj − Tr )

where A and B are raw materials, C and D are desirable product and undesirable by-product, respectively. The rate constant k1 and k2 are temperature dependence according to the Arrhenius relation.

A=

2.1. Reactor model The batch reactor is modeled by the following equations: Material balances in the reactor dMA = −k1 MA MB − k2 MA MC (1) dt dMB = −k1 MA MB (2) dt dMC = +k1 MA MB − k2 MA MC (3) dt dMD = +k2 MA MC (4) dt Energy balances around the reactor Qr + Q j dTr = (5) dt Mr Cpr Fj ρj Cpj (Tjsp − Tj ) − Qj dTj = dt Vj ρj Cpj

(6)

TRC

TJC

Tr

Tj

Cooling fluid Heating fluid Heat exchanger

Batch reactor

Fig. 1. Batch reactor system.

2W ρr

where Mi is the amount of mole of component “i”, Tr is the reactor temperature, Tj is the jacket temperature, and Tjsp is a set point value of the jacket temperature control system. The meaning of other variables and parameters are explained in the nomenclature. The dynamic behavior of the reactor can be simulated by solving Eqs. (1)–(6). The differential-algebraic solver DASSL [14] is used to give the solution of these equations. The initial conditions for MA , MB , MC , MD used in all simulation studies are 12, 12, 0, and 0 kmol, respectively. The initial values of both reactor and jacket temperature are set to 20 ◦ C. Other process parameter values used in the reactor models are listed in Table 1. In this work, it was assumed that all temperature can be measured frequently without delay. The sampling time of Table 1 Process parameter values MWA = 30 kg/kmol MWB = 100 kg/kmol MWC = 130 kg/kmol MWD = 160 kg/kmol k11 = 20.9057 k12 = 10000 k21 = 38.9057 k22 = 17000 r = 0.5 m Fj = 0.348 m3 /min U = 40.842 kJ/(min m2 ◦ C) CpA = 75.31 kJ/(kmol ◦ C) CpB = 167.36 kJ/(kmol ◦ C) CpC = 217.57 kJ/(kmol ◦ C) CpD = 334.73 kJ/(kmol ◦ C) H1 = −41840 kJ/kmol H2 = −25105 kJ/kmol ρ = 1000 kg/m3 ρj = 1000 kg/m3 Cpj = 1.8828 kJ/(kg ◦ C) Vj = 0.6912 m3

104

A. Arpornwichanop et al. / Chemical Engineering and Processing 44 (2005) 101–114

0.2 min is used for temperature measurements. It is noted that measurement noise is also included in the temperature sensors having Gaussian noise with zero mean and 1 ◦ C standard deviation. In addition, it was assumed that the concentration (amount of mol) of reactants in the reactor is measured infrequently and has a sampling time and measurement delay of 10 min.

3. On-line dynamic optimization strategy The aim of a dynamic optimization problem is to determine a control profile minimizing (or maximizing) a given objective function subject to process constraints. With the optimal control policy, the controlled system is driven from the initial state to a final desired state in an optimal way. However, in the presence of modeling error, the pre-specified control profile may lose its optimal character [2]. For this reason, an on-line optimization strategy through a receding horizon scheme is employed in this work to compensate such an error. The basic concept is to compute the optimal control profile based on current feedback information. However, only the initial value of the optimal trajectory is sent to the system as a set point for the controller. After new information of states is available from either measurement or estimation, the optimization is repeated again to generate updated optimal set point profile at the next time interval. The method proposed for improving the batch operation can be divided into two phases: on-line modification of the reactor temperature trajectory and on-line tracking of the desired temperature trajectory. The first phase involves determining an optimal temperature set point profile by solving the on-line dynamic optimization problem and will be described in this section. The other phase involves designing a nonlinear model-based controller to track the obtained temperature set point and will be presented in the next section. Since both the on-line dynamic optimization and the model-based control strategy rely on process models, the knowledge of current states and/or model parameters is required. However, in most industrial processes, state variables are not all measurable and some parameters are not known exactly. As a consequence, there is a need for estimating these states and parameters. In this work, two Extended Kalman Filters (EKF) are implemented. The first one is applied to predict the reactant concentration, which will be used for on-line dynamic optimization, from its delayed measurement. The other one is applied to estimate the unknown heat of reaction, which will be used for model-based controller, from the frequently available measurements of temperature. 3.1. Problem statements

the more efficient approach to obtain a maximum yield in a minimum time or minimum cost and also to reach the specific final conditions of the products in terms of quality and quantity. In general, an objective function in the optimization problem can be chosen, depending on the nature of the problem. Here, two practical optimization problems related to batch operation: maximization of product concentration in a fixed batch time and minimization of batch operation time given amount of desired product, are considered to determine an optimal reactor temperature profile. The first problem formulation is applied to a situation where we need to increase the amount of desired product while batch operation time is fixed. This is due to the limitation of complete production line in a sequential processing. However, in some circumstances, we need to reduce the duration of batch run to allow the operation of more runs per day. This requirement leads to the minimum time optimization problem. These problems can be described in details as follows. 3.1.1. Maximum product concentration problem (P1) In this type of problem, the objective is to compute the optimal temperature policy maximizing the amount of a desired product concentration for a given fixed batch time subject to bounds on the reactor temperature. The problem can be written mathematically as max J = X(tf ) T(t)

subject to x˙ = f(x(t), T, p, t) x(t0 ) = x(0)

process model

initial conditions

TL ≤ T ≤ TU tf = tf∗ where X is the amount of the desired product at a given final batch time, x is the vector of state variables, x˙ is the derivative of x with respect to time (t), T is the reactor temperature, p are process parameters, tf∗ is the fixed batch time, and TL and TU are lower and upper bounds of the reactor temperature. 3.1.2. Minimum batch time problem (P2) The purpose of this optimization problem is to determine the optimal temperature profiles to achieve the desired final product concentration in minimum batch time, thus the performance index is the final time whereas the desired production concentration is defined as a terminal constraints. The formulation of the minimum batch time problem can be shown as min J = tf T(t)

As a batch reactor is utilized for the production of a wide variety of high value products, an optimization of batch operating conditions, e.g. temperature, operating time, etc. is

subject to x˙ = f(x(t), T, p, t)

process model

A. Arpornwichanop et al. / Chemical Engineering and Processing 44 (2005) 101–114

x(t0 ) = x(0)

initial conditions

TL ≤ T ≤ TU X(tf ) = X∗ where X∗ is the desired product concentration at the end of batch run and tf is final batch time. 3.2. Solution of dynamic optimization problems Computational techniques for the solution of a dynamic optimization problem as formulated above have been an active area of research. There are a number of different techniques that have been proposed in the literature to solve this class of problems. In general, they are mainly classified into three classes. The first one is based on a classical variation method. This approach is also known as an indirect method as it focuses on obtaining the solution of the necessary conditions rather than solving the optimization directly. Solution of these conditions often results in a two-point boundary value problem (TPBVP), which is accepted that it is difficult to solve [15]. Although several numerical techniques have been developed to address the solution of TPBVP, e.g. control vector iteration (CVI) and single/multiple shooting method, these methods are generally based on an iterative integration of the state and adjoint equations and are usually inefficient [16]. Another difficulty relies on the fact that it requires an analytical differentiation to derive the necessary conditions. The second class of solutions is based on dynamic programming. Unlike the variation method, this approach applies the principle of optimality to formulate an optimization problem, leading to the development of the Hamilton–Jacobi–Bellman partial equations that determine the solution of the optimal control problem. However, this approach is quite limited to a simple control problem [17] because of a difficulty in obtaining the solution of the optimality equations. [18] extended the idea of the optimality principle to develop an alternative technique, named as iterative dynamic programming (IDP). Although the implementation of the IDP for solving many optimal control problems can be found in the literature, it is known that the IDP algorithm would be slower than most other gradient-based algorithms [19]. The last one is based on discretization techniques, received major attention and considered as an efficient solution method. The concept of this approach is to transform the original optimal control problem into a finite dimensional optimization problem, typically a nonlinear programming problem (NLP). Then, the optimal control solution is given by applying a standard NLP solver to directly solve the optimization problem. For this reason, the method is known as a direct method. The transformation of the problem can be made by using discretization technique on either only control variables (partial discretization) or both state and control variables (complete discretization). Based on this con-

105

sideration, this approach can be divided into two categories: sequential and simultaneous strategy. In the sequential strategy, a control (manipulated) variable profile is discretized over a time interval. The discretized control profile can be represented as a piecewise constant, a piecewise linear, or a piecewise polynomial function. The parameters in such functions and the length of time subinterval become decision variables in optimization problem. This strategy is also referred to a control vector parameterization (CVP). One advantage in the sequential approach is that only the parameters that are used to discretize the control variable profile are considered as the decision variables. The optimization formulated by this approach is a small scale NLP that makes it attractive to apply for solving the optimal control with large dimensional systems that are modeled by a large number of differential equations. In addition, this approach can take the advantage of available IVP solvers. However, the limitation of the sequential method is a difficulty to handle a constraint on state variables (path constraint). This is because the state variables are not directly included in NLP. In contrast to the sequential solution method, the simultaneous strategy solves the dynamic process model and the optimization problem at one step. This avoids solving the model equations at each iteration in the optimization algorithm as in the sequential approach. In this approach, the dynamic process model constraints in the optimal control problem are transformed to a set of algebraic equations which is treated as equality constraints in NLP problem [20]. To apply the simultaneous strategy, both state and control variable profiles are discretized by approximating functions and treated as the decision variables in optimization algorithms. The main advantage of the simultaneous approach is a capability in handing constraints on the state variables. This is because these constraints can be dealt with by including them directly in the NLP as additional constraints. However, due to the discretization on both state and control variables, this leads the simultaneous approach to a large scale optimization problem consisting of a large set of algebraic constraints and decision variables and needs a special solution strategy. As the ability of the sequential approach to handle large systems without the need to solve excessively large optimization problem, this approach is utilized in this study to solve the optimal control problem. The formulation of the optimal control problem as a nonlinear programming is described below. 3.2.1. Sequential approach Consider a dynamic process model as in the form of an implicit function f(t, x˙ (t), x(t), u(t), p) = 0

[t0 , tF ]

(7)

where x(t) is the set of all state variables, x˙ denotes the derivative of x(t) with respect to time, u(t) is a vector of con-

106

A. Arpornwichanop et al. / Chemical Engineering and Processing 44 (2005) 101–114

trol variables, and p is a vector of time invariant parameters. The time interval of interest is [t0 , tF ] and the function f is assumed to be continuously differentiable with respect to all its arguments. For given initial condition x(t0 ), the optimal control of a dynamic process represented by Eq. (7) can be determined by computing u(t) maximizing (or minimizing) an objective function J of the form J = F(˙x(tF ); x(tF ); u(tF ); tF )

(8)

subject to bounds on u(t) and terminal constraints. To transform such a dynamic optimization problem into a nonlinear programming problem via the sequential approach, the control u(t) is approximated by a finite dimensional representation. The time interval [t0 , tF ] is divided into a finite number of subintervals (P). In each subinterval, the control u(t) is represented by a set of basis functions involving a finite number of parameters u(t) = φj (t, zj ),

t ∈ [tj−1 , tj ]

with j = 1, 2, . . . , J (9)

where tJ = tF . The control profile is defined by the parameters zj and switching times tj . In this study, the piecewise constant control is assumed and used because the form of the solution is ideally suited for implementation on a digital computer. Thus the set of decision variables for the nonlinear program can be written as y = {z1 , z2 , . . . , zJ , t1 , t2 , . . . , tJ }

3.3. State and parameter estimation The implementation of the on-line optimization strategy requires the knowledge of current states and/or parameters in nonlinear process models in order to modify a new optimal profile defined as a set point for a controller. It is known that some measurements i.e. concentration are available at low sampling rate with significant time delay. To overcome this difficulty, state and parameter estimation is incorporated into the proposed on-line optimization algorithm. In this work, an Extended Kalman Filter (EKF), an extension of the Kalman Filter, is designed to reconstruct the current state variables from the delayed state measurement. The advantage of the EKF is that it requires information only from the previous sampling time and allows prior knowledge of a system via process models to be used for the estimation. The algorithm of the EKF can be seen in Appendix A.

(10)

The resulting nonlinear programming problem of the form min(or max ) J(y)

(11)

y

subject to equality constraints (process model), inequality constraints, bounds on control, etc. is solved using an effective decomposition algorithm named the successive reduced quadratic programming developed by [21]. Computational procedure of the sequential approach for the optimization problem (Eq. (11)) is illustrated in Fig. 2. With the initial guess of the decision variables (y), an integrator based on Set initial conditions Guess initial control parameters

Model Solver Evaluate Objective function Constraints Gradient

no

NLP

Gear’s type method is used for solving the process model which provides the value of the objective function and constraints. Gradient information of the objective function and constraints with respect to the decision variables is evaluated in an efficient way using adjoint variable approach. Then, a NLP solver determines a new set of control parameters and sends it back to the model solver. This procedure is repeated until the optimal value is found satisfying a specified accuracy.

3.3.1. Application to the batch reactor Since the concentrations of reactants (MA , MB , MC and MD ) in batch reactor are assumed to be measurable with a delay of one sampling time; that is, at time k of the reactor, only information at time k − 1 is available. Thus, the EKF is applied to estimate the value of reactant concentration at current time k from their delayed measurements at sampling time k − 1. However, since it is expected to exhibit uncertainty in reaction rate constant (i.e. k21 and k22 ) in real plant, the EKF is also used to estimate these uncertain parameters. The following equations, therefore, are appended for parameter estimation dk21 =0 dt

(12)

dk22 =0 dt

(13)

Eqs. (1)–(4), (12) and (13) correspond to Eq. (A.1) in the EKF algorithm. Based on the estimate of the current information, the dynamic optimization problem is resolved to generate a new optimal temperature trajectory.

Check tolerance yes

Optimal control parameters

Fig. 2. Optimal control solution via the sequential approach.

4. Generic model control (GMC) This section presents the design of a controller to control the reactor temperature following a desired temperature tra-

A. Arpornwichanop et al. / Chemical Engineering and Processing 44 (2005) 101–114

107

jectory. It is accepted that the use of linear control techniques in highly nonlinear and time variant chemical processes (e.g. batch reactors) is quite limited to their performances and may give a poor control response. Therefore, in this work, a nonlinear control technique based on generic model control (GMC) algorithm is utilized. This model-based control methodology has been received much interest during the last decade and a number of applications of GMC to the control of batch processes have been reported in the literature (e.g. [4,8,22,23], etc.). However, most of these works focus on using the GMC to track the pre-determined optimal profile (off-line calculation) of a reactor temperature. No effort has been made to apply the GMC to implement an on-line optimal set point profile. Therefore, in this work, the performance of the GMC controller with optimal temperature set points determined by on-line optimization strategy is also evaluated and compared to that of the GMC controller with pre-determined set point.

is negligible compared to the heat transferred in the reactor, the energy balance equation becomes

4.1. Control algorithm

Replacing these equations in Eq. (17), we have  Wr Cpr K1 (Trsp − Tr ) Tj = Tr + U r Ar   t Qr + K2 (Trsp − Tr ) dt − Ur A r 0

Let us consider a process based on the following model equations: dx = f(x, p, t) + g(x, t)u dt

(14)

y = h(x)

(15)

where x is a vector of state variables, y is a vector of output variables, u is a vector of input variables, p is vector of process parameters, and f, g, and h are generally nonlinear function vectors. The general form of GMC algorithm can be written as  tf dy = K1 (ysp − y) + K2 (ysp − y) dt (16) dt 0 The GMC control response can be designed via the tuning parameters K1 and K2 based on the tuning curve given by [24]. It should be noted that the GMC approach is a special case of the global input output linearizing control technique in which a transformed control action is chosen properly with the external PI controller. The use of Eq. (16) forces y toward its set point, ysp , with zero offset. If Eq. (15) is differentiated, and the Eq. (16) is substituted into the resulting equation, the GMC control law is  [K1 (ysp − y) + K2 (ysp − y) dt − (dh/dx)f(x, d, t)] u= (dh/dx)g(x, t) (17) 4.2. Application of GMC controller to the batch reactor To implement the GMC, an energy balance around the reactor is required; it gives the relation between the reactor temperature (controlled variable) and the jacket temperature (manipulated variable). Based on the assumption that the amount of the heat accumulated in the walls of the reactor

Qr + Ur Ar (Tj − Tr ) dTr = dt Wr Cpr

(18)

where Ur is the heat transfer coefficient, Ar is the heat transfer area, Wr is the mass of the reactor contents, Cpr is the mass heat capacity of the reactor content, and Qr is the heat released by the reactions. Rearranging the Eq. (18) as in the form of GMC algorithm, the following functions, f, g, and h can be defined f(x, p, t) = g(x, t) =

Qr − Ur Ar Tr Wr Cpr

Ur Ar Wr Cpr

h(x) = Tr

(19) (20) (21)

(22)

In order to make the GMC control law available for an on-line control implementation, the integral term in Eq. (22) is approximated by numerical integration. This leads to the discrete-time form of the GMC algorithm as given in the following equation.
Wr Cpr Tj (k) = Tr (k) + K1 (Trsp − Tr (k)) U r Ar  k  Qr (k) (23) + K2 (Trsp − Tr (k)) t − Ur A r 0

where t is the sampling time chosen to be equal to the frequency of the temperature measurement. It is sufficiently small to approximate the continuous-time form without requiring a large number of sampled values. However, Eq. (23) gives the actual jacket temperature (Tj (k)) required at the next sampling time to control the reactor temperature (Tr (k)) at the desired trajectory (Trsp ). As in usual practice, the reactor temperature control is cascaded with the jacket temperature control (heating and cooling system); the output of reactor temperature controller (master loop) is the set point of the jacket temperature controller (slave loop), as demonstrated in Fig. 3. In addition, as the model of the heat exchanger system for heating and cooling is not included in Eq. (18), if Tj (k) is applied directly as the set point for jacket temperature control system without considering on its dynamic, the resulting control response would be sluggish. To accommodate such an effect, it is reasonable to assume that the dynamic of the jacket control system

108

A. Arpornwichanop et al. / Chemical Engineering and Processing 44 (2005) 101–114 Initial setpoint

Tj

Dynamic compensation

GMC

Tjsp

Tr(k) Tj(k)

REACTOR

Mi (K-1) Qr(k)

EKF with simplified reactor model

EKF

Updated Mi (K) Trsp (K) New reactor temperature setpoint

On-line dynamic optimization problem

Fig. 3. The proposed strategy for on-line update and control of reactor temperature profile.

can be approximated by a first order model with time constant (τ j ) [25]. Consequently, the Tjsp (k) can be computed by   Tj (k) − Tj (k − 1) Tjsp (k) = Tj (k − 1) + τj (24) t

states as shown below is used instead. Qr + Ur Ar (Tj − Tr ) dTr = dt Wr Cpr

With this jacket temperature set point, a jacket temperature controller (setting as a PI controller) through a heat exchanger system opens or closes control valve reflecting to the flowrate of heating and cooling fluid. However, in reality the ability of the heat exchanger in adjusting the jacket temperature is always limited, thus, in this work the jacket temperature is bounded between 0 ◦ C and 120 ◦ C. The tuning parameters of GMC are given in Table 2. As shown in Eq. (23), the heat released (Qr (k)), which cannot be measured, is needed in the GMC algorithm. Here, the EKF algorithm, as used in on-line optimization strategy, coupled with the simplified reactor model, given by [26], is also applied to estimate the heat released (Qr (k)). The reason of using the simplified model, not the exact model of the plant, is because if the exact model were used, too many uncertain/unknown parameters as well as too many unmeasurable states would be involved. That may lead to poor performance of the EKF. Hence, the simplified model with less uncertain/unknown parameters and unmeasurable

dN (27) = −bNTr dt dQr dTr dN =N + Tr (28) dt dt dt db =0 (29) dt where N = −bMr (H), b is a pseudo reaction rate constant, Mr is the total reactant concentration, and H is heat of reaction. It should be noted that the variable N representing two unknown parameters (Mr and H) can be estimated instead of these parameters so that the number of state equations for estimation decrease and an error of estimation corresponding to the uncertainty of each parameters can be reduced. Eqs. (25)–(29) correspond to Eq. (A.1) in the EKF algorithm. Once the reactor and jacket temperature measurement are available, the EKF with simplified model estimates the heat release of reactions (Qr (k)). Table 3 summarizes the initial conditions and tuning parameters of the EKF used in this simulation work.

Table 2 Parameters in GMC algorithm Wr = 1560 kg Cpr = 1.8828 kJ/(kg ◦ C) Ur = 40.842 kJ/(min m2 ◦ C) Ar = 6.24 m2 τ j = 2 min K1 = 2.4 K2 = 10−4

Fj ρj Cpj (Tjsp − Tj ) − Ur Ar (Tj − Tr ) dTj = Vj ρj Cpj dt

(25) (26)

5. Simulation results 5.1. Maximum conversion problem (P1) All simulation results given here are based on the optimization problem P1. The objective in the problem formulation is to find the optimal reactor temperature profile, such

A. Arpornwichanop et al. / Chemical Engineering and Processing 44 (2005) 101–114 Table 3 Parameters and initial conditions in EKF

Estimate Qr Tr (0) = 20 ◦ C Tj (0) = 20 ◦ C N(0) = 1.8462 Qr (0) = 0 b(0) = 1.8386 × 10−6

98

P = diag[100 10 100 100 100 2000] Q = diag[100 1 100 1 500 5000] R = diag[10 10 10 10]

P = diag[1 1 100 20 10] Q = diag[10 10 2500 100 100] R = diag[10 10]

96 5

Temperature (C)

Estimate Mi , k21 and k22 MA (0) = 12 kmol MB (0) = 12 kmol MC (0) = 0 kmol MD (0) = 0 kmol k21 (0) = 38.9057 k22 (0) = 17000

94 4 92

that the amount of mole of product C is maximized in a fixed batch time with respect to a constraint on the temperature due to the reason of safe operation. In this case study, the specified final batch time (tf ) of 200 min is used and the reactor temperature is bounded according to 20 ≤ T (◦ C) ≤ 120. 5.1.1. Temperature set point profile determined off-line with perfect tracking First, simulation studies have investigated the case where the optimal temperature profile is determined by off-line computation and perfect tracking of such a profile is assumed. This results in the maximum amount of product C (maximum conversion) that can be achieved at the end of batch run and is served as a reference to be compared with results obtained from the proposed strategy. The optimal control problems have been solved using time interval with equal length varied from one to 40 intervals to discretize the profile. The switching time is fixed and the length of each interval is specified by dividing the batch operation time (tf ) by a number of time intervals (P). Thus, the problem is to seek an optimal temperature value (decision variables) in each subinterval. Simulation results with different time interval (P) are reported in Table 4. Optimal control policy in reactor temperature for each case is shown in Fig. 4. As shown in Table 4, when one time interval (P = 1) is used, the amount of product C obtained at the final time (tf = 200 min) is 7.0171 kmol and the optimal temperature (isothermal operation) set point is 88.01 ◦ C whereas using P = 20, the amount of product C achieved is 7.0379 kmol. It was found from Table 4 Summary of the results: off-line optimization and perfect tracking Time interval

Product C (kmol)

By-product D (kmol)

CPU time∗ (s)

1 5 10 20 40

7.0171 7.0281 7.0339 7.0379 7.0402

1.3464 1.3605 1.3594 1.3585 1.3579

0.1591 1.2600 1.9778 2.8066 5.0519

Run on Pentium III/850 computer.

3

90

1

2

88

86

∗

109

0

20

40

60

80

100

120

140

160

180

200

Time (min) Fig. 4. Optimal temperature profile: 1 interval (1), 5 intervals (2), 10 intervals (3), 20 intervals (4), 40 intervals (5)

these results that the amount of the desired product C increases as a number of interval increases. This is due to that as the number of intervals enlarges (more degrees of freedom in the optimization), the approximated optimal profile with piecewise constant policy is closer to the actual optimal profile. 5.1.2. Temperature set point profile determined on-line with GMC controller Next, the proposed strategy using an on-line dynamic optimization to update optimal temperature set point profile is implemented. Rather than assuming the reactor temperature trajectory can be tracked perfectly as in previous studies, the GMC controller is applied here to drive the system to follow the desired trajectory. Considering the time elapsed in the determination of the optimal control problem, with P = 20, the temperature set point profile is updated every 10 min. This means that, to apply this strategy on-line, the computational time for the updated temperature set point profile must be less than 10 min. Here, the computational time based on Pentium III/850 mHz is approximately 3 s. As a result, this strategy is applicable for on-line implementation. Regarding to the GMC control performance, it was found that the GMC controller can drive the system from current set point value to a new one and maintain it at this set point. Therefore, the GMC controller can be used for tracking the profile obtained from the proposed strategy. Simulation studies carried out to compare the amount of desired product C obtained from on-line dynamic optimization strategy with that from off-line strategy, are cases where the perfect model (all parameters correctly specified) is used (nominal case), and where plant/model mismatch is introduced by changing parameters in actual plant i.e. pre-exponential rate constant (k0 ) decreased by 50% and activation energy (Ea) increased by 20% from their nominal values, as shown in Table 5.

110

A. Arpornwichanop et al. / Chemical Engineering and Processing 44 (2005) 101–114

Table 5 Comparisons of the results obtained from off-line and on-line optimization strategy with GMC controller (problem P1) Case studies

Product C (kmol) Off-line

On-line

1. Nominal case All parameters specified correctly

6.9478

2. Plants/model mismatch case −50% k0 of reaction 2 in plant model (k21 = 38.2125)

7.6751

7.8851

+20% Ea of reaction 2 in plant model (k22 = 20400) −50% k0 and +20% Ea of reaction 2 in plant model

8.5825 8.5827

10.2131 10.2137

6.9584

Note: k0 = exp(k21 ) and Ea = (k22 )(R) where k0 is pre-exponential rate constant, Ea = activation energy, and R is ideal gas constant.

In the nominal case, the product obtained from the off-line strategy (C = 6.9478) is close to that obtained by the on-line strategy (C = 6.9584). Fig. 5(a) and (b) shows the response of GMC controller to track the reactor temperature trajectory that is pre-specified by off-line calculation, and the compar-

Fig. 5. (a) Control response (off-line, nominal case): Trsp (dash), Tr (solid), Tj (dot). (b) Heat released (off-line, nominal case): actual (solid), estimated (dot).

ison of the actual and estimated heat released by reactions, respectively. It can be seen that the EKF provides an excellent estimation of the heat released. With this heat released, the GMC controller gives very good temperature control. Similarly, in the case that optimal temperature is modified via the on-line optimization strategy based on the current information of Mi , k21 and k22 obtained from delayed measurement of Mi , the GMC controller is able to successfully track the reactor temperature (Fig. 6(b)). The performance of the EKF to predict the amount of MA , MB , MC , and MD at current time from their measurement with time delay of 10 min is illustrated in Fig. 6(a). Also, the EKF gives good estimations of k21 and k22 as shown in Fig. 6(c). It is interesting to note that since the initial reactor temperature starts at 20 ◦ C which is below the optimal temperature set point, the GMC controller attempts to drive the reactor temperature from this condition to the specified set point. However, due to the dynamics of the reactor temperature, the reactor temperature differs from the optimal set point during the first 20 min, causing that the amount of the product C obtained from the on-line strategy with the GMC controller (C = 6.9584) is less than that obtained from the off-line strategy with perfect tracking (see Table 4 for P = 20, C = 7.0379) For the mismatch in k21 , the value of the desired product C = 7.8851 can be achieved at the end of batch for the on-line optimization strategy which is higher than that obtained from off-line strategy where the mismatch is not noticed (C = 7.6751). Similar results can be observed under the case of plant model mismatch in k22 as shown in Table 5. These results indicate clearly that the performance of batch reactor operation is improved via the proposed strategy. Due to similarity in their control responses, only the result for change in k21 is shown in Fig. 7. Finally, with a change in both k21 (−50% k0 ) and k22 (+20% Ea) in plant model, the results using the on-line optimization strategy show that the GMC controller is able to accommodate this change very well as can be seen in Fig. 8(a). Fig. 8(b) presents the performance of the EKF for estimation of k21 and k22 . Since the EKF estimates these parameters close to the true values, the mismatch is eliminated. That leads to high product C obtained at the final batch time (C = 10.2137) compared to the value of C = 8.5827 obtained from the off-line optimization strategy.

A. Arpornwichanop et al. / Chemical Engineering and Processing 44 (2005) 101–114

111

Fig. 6. (a) Product profile (on-line, nominal case): actual (solid), estimated (䊐). (b) Control response (on-line, nominal case): Trsp (dash), Tr (solid), Tj (dot). (c) Estimate of k21 and k22 (on-line, nominal case): actual (solid), estimated (䊐).

5.2. Minimum time problem (P2)

Fig. 7. Control response (on-line, mismatch in k21 ): Trsp (dash), Tr (solid), Tj (dot).

The results presented here correspond to the case where the objective is to minimize the batch time of operation subject to the terminal constraints on the desired amount of mole of product C (MC (tf ) = 6.00 kmol). The reactor temperature constraints are same as in problem P1. Several simulations have been carried out under process parameter uncertainties e.g. in pre-exponential rate constant (k0 ) and activation energy (Ea). In all case studies we considered 10 time intervals when reactor temperature and switching time are optimized while minimizing the final batch operation time. Results, reported in the value of minimum batch time to obtain the desired product C and the amount of the desired product C at the end of batch operation, from on-line dynamic optimization strategy are also compared with those from the off-line strategy. With an 20% increase of parameter k21 in plant model, it can be seen from Table 6 that the final batch time needed to achieve the desired product C from the proposed on-line

112

A. Arpornwichanop et al. / Chemical Engineering and Processing 44 (2005) 101–114

Fig. 9. Control response (off-line, mismatch in k21 , problem P2): Trsp (dash), Tr (solid), Tj (dot).

Fig. 8. (a) Control response (on-line, mismatch in k21 and k22 ): Trsp (dash), Tr (solid), Tj (dot). (b) Estimate of k21 and k22 (on-line, mismatch in k21 and k22 ): actual (solid), estimated (䊐).

modification of temperature set point profile (tf = 49.4 min) is shorter compared to the result with the off-line strategy (tf = 64.8 min). This is because the EKF can acknowledge this parameter uncertainty, so that the temperature set point profile is updated corresponding to a modified parameter value close to the actual value. Fig. 9 shows, for this mismatch, the control response of GMC controller to deliver the reactor temperature from initial condition to the temperature set point determined off-line. For the on-line temperature set

Fig. 10. Control response (on-line, mismatch in k21 , problem P2): Trsp (dash), Tr (solid), Tj (dot).

point modification, as expected, GMC controller can also control the reactor corresponding to the changes in temperature set point as can be seen in Fig. 10 The results for the remaining case studies are summarized in Table 6. An important aspect obtained from these results is that in all cases, the minimum batch time to obtain the desired product concentration in the on-line set point

Table 6 Comparisons of the results obtained from off-line and on-line optimization strategy with GMC controller (problem P2) Case studies

(1) −50% (1) +20% (1) −50%

k0 of reaction 2 in plant model (k21 = 38.2125) Ea of reaction 2 in plant model (k22 = 20400) k0 and + 20% Ea of reaction 2 in plant model

Final time (min)/product C (kmol) Off-line

On-line

64.8/6.0049

49.4/6.0089

56.4/6.0147 56.7/6.0144

44.6/6.0104 45.2/6.0121

A. Arpornwichanop et al. / Chemical Engineering and Processing 44 (2005) 101–114

modification strategy decreases. This points out the effectiveness of the proposed method to improve the operation of batch reactor.

6. Conclusions In this work, the method using an on-line dynamic optimization and control strategy to enhance batch reactor operation has been proposed. The on-line dynamic optimization by the idea of receding horizon scheme is performed. The approach is based on the updated current information of states of the system which are estimated from the delayed measurement of the amount of reactants in the reactor using an Extended Kalman Filter (EKF) technique, at specified time interval to provide a new updated optimal reactor temperature set point profile. Two types of optimization formulation related to batch operation (maximum concentration and minimum time problem) were considered in the proposed on-line set point modification strategy. The obtained reactor temperature set point was implemented using a generic model control (GMC). The EKF was also incorporated into the GMC algorithm in order to estimate the heat released by reactions using the direct measurement of reactor and jacket temperature. A batch reactor with highly exothermic reactions was used as a simulation case study to demonstrate the effectiveness of the proposed approach. Simulation studies have been carried out in both nominal case and plant/model mismatch case and the results show that the performance of the batch reactor in terms of the amount of a desired product and batch operation time can be improved significantly by the proposed strategy. In addition, they also clearly indicate the capability of the GMC controller to control the reactor temperature along the specified trajectory and that of the EKF to estimate the states and parameters of the system.

Acknowledgements The financial support to A. Arpornwichanop through local graduate scholarship from the national science and technology development agency (NSTDA) is gratefully acknowledged. A. Arpornwichanop also acknowledges the support of the Department of Chemical Engineering, University of Bradford, UK during his stay as a visiting researcher at the University in 2000.

Appendix A. The EKF algorithm Since in this work the process models used are in time continuous form and measurements are in discrete form, the EKF in discrete/continuous formulation [27] is used. The basic algorithm of the EKF can be summarized as follows:

113

For nonlinear systems, the process model can be described by differential equations x˙ = F(x(t), u(t), t) + ζ(t)

(A.1)

y = G(x(t)) + η(t)

(A.2)

where F is a vector of system function, G is a vector of measurement function, ζ is a zero mean Gaussian process noise with covariance Q, and η is a zero mean Gaussian measurement noise with covariance R. The equations for the EKF are given by a set of correction and prediction equations as shown in the following: Correction phase: Correct the prior estimates of state at k − 1 and update the weighting matrix Kk−1 = Pk−1/k−2 Ck−1 T [Ck−1 Pk−1/k−2 Ck−1 T + R]−1 (A.3) xˆ k−1/k−1 = xˆ k−1/k−2 + Kk−1 [yk−1 Ck−1 xˆ k−1/k−2 ] (A.4) Pˆ k−1/k−1 = [I − Kk−1 Ck−1 ]Pˆ k−1/k−2 [I − Kk−1 Ck−1 ]T + Kk−1 RKTk−1

(A.5)

Prediction phase: Integrate the nonlinear state and covariance equations from time k − 1 to k in order to acquire estimate xˆ k/k−1 and Pˆ k/k−1 x˙ˆ = F(x, u)

(A.6)

ˆ Tk−1 + Q P˙ˆ = Ak−1 Pˆ + PA

(A.7)

where xˆ k/k−1 denotes the estimate of state x at t = k from information at t = k − 1, K is Kalman gain matrix, P is covariance matrix, and matrices

   ∂F  ∂G  Ak−1 = and C = ∂x xˆ k−1|k−1 ∂x xˆ k−1|k−2 are the Jacobians of function F and G with respect to the state vector, respectively.

Appendix B. Nomenclature ki ki1 ki2 t t u x y A Cp Cpi

rate constant for reaction i (kmol−1 /s) rate constant 1 for reaction i rate constant 2 for reaction i time (min) sampling time (min) input variables state variables output variables heat transfer area (m2 ) mass heat capacity (kJ/(kg ◦ C)) molar heat capacity of component i (kJ/(kmol ◦ C))

114

F Hi K1 , K2 Mi MWi P Q Qr R Ri T U V W ρ

A. Arpornwichanop et al. / Chemical Engineering and Processing 44 (2005) 101–114

flowrate (m3 /min) heat of reaction of reaction i (kJ/kmol) GMC controller parameters number of moles of component i (kmol) molecular weight of component i (kg/kmol) covariance matrix covariance matrix of process noise heat released from reactions (kJ/min) covariance matrix of measurement noise rate of reaction i (kmol/min) reactor temperature (◦ C) heat transfer coefficient (kJ/(min m2 ◦ C)) reactor volume (m3 ) Reactor content (kg) density (kg/m3 )

Subscripts f filter j jacket r reactor sp set point References [1] G.E. Rotstein, D.R. Lewin, Control of an unstable batch chemical reactor, Comp. Chem. Eng. 16 (1) (1992) 27–49. [2] D. Bonvin, Optimal operation of batch reactors – a personal view, J. Proc. Cont. 8 (5/6) (1998) 355–368. [3] N. Aziz, I.M. Mujtaba, Optimal operation policies in batch reactors, Chem. Eng. J. 85 (2002) 313–325. [4] N. Aziz, M.A. Hussain, I.M. Mujtaba, Performance of different types of controllers in tracking optimal temperature profiles in batch reactors, Comp. Chem. Eng. 24 (2000) 1069–1075. [5] C. Loeblein, J.D. Perkins, B. Srinivasan, D. Bonvin, Performance analysis of on-line batch optimization systems, Comp. Chem. Eng. 21 (Suppl) (1997) s867–s872. [6] D. Ruppen, D. Bonvin, D.W.T. Rippin, Implementation of adaptive optimal control operation for a semi-batch reaction system, Comp. Chem. Eng. 22 (1998) 185–189. [7] C. Kravaris, R.A. Wright, J.F. Carrier, Nonlinear controller for trajectory tracking in batch processes, Comp. Chem. Eng. 13 (1/2) (1989) 73–82. [8] B.J. Cott, S. Macchietto, Temperature control of exothermic batch reactors using generic model control, Ind. Eng. Chem. Res. 28 (1989) 1177–1184.

[9] Z.H. Liu, S. Macchietto, Model based control of a multipurpose batch reactor – an experimental study, Comp. Chem. Eng. 19(Suppl) (1995) s477–s482. [10] S. Yuce, A. Hasaltun, S. Erdogan, M. Alpbaz, Temperature control of a batch polymerization reactor, Trans. IChemE 77 (Part A) (1999) 413–420. [11] P. Kittisupakorn, M.A. Hussain, A. Arpornwichanop, Temperature control of an exothermic batch reactor by model predictive exothermic batch reactor by model predictive control, in: Proceedings of the Third Asian Control Conference, Shanghai, China, 4–7 July 2000, pp. 2506–2511. [12] N. Aziz, M.A. Hussain, I.M. Mujtaba, Optimal control of batch reactor: comparison of neural network based GMC and inverse model control approach, in: Proceedings of the Sixth World Congress of Chemical Engineering, Melbourne, Australia, 23–27 September 2001. [13] R. Berber, Control of batch reactors: a review, Trans. IChemE 74 (Part A) (1996) 3–20. [14] L.R. Petzold, A description of DASSL: a differential/algebraic system solver, SAND82-8637, Sandia National Laboratories, 1982. [15] W.H. Ray, Advanced Process Control, McGraw Hill, New York, 1981. [16] J.G. Renfro, A.M. Morshedi, O.A. Asbjornsen, Simultaneous optimization and solution of systems described by differential/algebraic equations, Comp. Chem. Eng. 11 (5) (1987) 503–517. [17] V. Nevistic, Ph.D. thesis, Constrained control of nonlinear systems, Swiss Federal Institute of Technology, 1997. [18] R. Luus, Optimal control by dynamic programming using systematic reduction in grid size, Int. J. Control 51 (1990) 995–1013. [19] S.A. Dadebo, K.B. Mcauley, Dynamic optimization of constrained chemical engineering problems using dynamic programming, Comp. Chem. Eng. 19 (5) (1995) 513–525. [20] J.E. Cuthrell, L.T. Biegler, On the optimization of differentialalgebraic process systems, AIChE J. 33 (8) (1987) 1257– 1270. [21] C.L. Chen, Ph.D. thesis, A class of successive quadratic programming methods for flowsheet optimization, University of London, 1988. [22] L.S. Kershenbaum, P. Kittisupakorn, The use of a partially simulated exothermic (PARSEX) reactor for experimental testing of control algorithms, Trans. IChemE. 72 (Part A) (1994) 55–63. [23] J.X. Shen, M.S. Chiu, Q.G. Wang, A comparative study of modelbased control techniques for batch crystallization process, J. Chem. Eng. Jpn. 32 (4) (1999) 456–464. [24] P.L. Lee, G.R. Sullivan, Generic model control (GMC), Comp. Chem. Eng. 12 (6) (1988) 573–580. [25] B.G. Liptak, Controlling and optimization chemical reactors, Chem. Eng. Magazine (May 1986) 69–81. [26] P. Kittisupakorn, Ph.D. thesis, The use of nonlinear model predictive control techniques for the control of a reactor with exothermic reactions, University of London, 1995. [27] A. Gelb, Applied Optimal Estimation, The MIT Press, Cambridge, MA, 1974.