- Email: [email protected]
Adaptive—predictive temperature control of semi-batch reactors
Cheudcol Engineering Science, Vol. 48, No. 19, pp. 3373-3382, Printed in Great Britain.
1993
ooo9-2509/93 66.00 + 0.00 Q 1993 Pergamon Press Ltd
ADAPTIVE-PREDICTIVE TEMPERATURE SEMI-BATCH REACTORS G. DEFAYE. Equipe
I.A.G.P.,
N. REGNIER,
URA 503 CNRS, Universite (Received
J. CHABANON,’
CONTROL
L. CARALP
and C. VIDAL
EC& Nationale Sup&ieure de Chimie et Physique Bordeaux I, 33405 Talence Cedex, France
18 Noumber
1992; acceptedfor
publication
OF
de Bordeaux,
30 April 1993)
Abstract-Although reactant addition induces unavoidable temperature changes in semi-batch reactors, it is shown that adaptive-predictive techniques may offer a means of balancing them successfully. To this end, appropriate data filtering must be integrated to the estimation algorithm, so that parameters remain insensitive to the nonrelevant dynamical perturbations generated by semi-batch operations. A pilot reactor is used to illustrate imorovements thus workable while carrying out a moderately exothermic copolymerization reaction.
I. INTRODUCTION
The control of batch reactors is more difficult to achieve than that of continuous processes (Juba and Hamer, 1986) due to the lack of steady state. Indeed, the dynamics of a batch process is essentially timedependent. PID regulators, which are well fitted to maintain a fixed setpoint, often fail to control the batch processes properly. In fact, such a task requires controllers able to face, without adjustment, a wide range of working conditions. Liptak (1986) showed how to improve the performances of PlD regulators, in particular during the initial phase of a batch process. However, the so-called “batch unit” and “dual mode” are nothing much more than a palliative. Better results can be obtained by simply giving up PIDs. Cluett et al. (1985) thus developed a self-tuning algorithm based on a transfer function model which connects the reactor temperature to the flow rate of a circulating thermofluid. Cott and Macchietto (1989) designed another procedure, involving the on-line estimation of the rate at which heat is released by the ongoing reaction. The controller may even include nonlinear models (Lee and Sullivan, 1988). Besides the improvements set up in this way, the prominent outcome of this kind of study is the following conclusion: nothing but adaptive predictive controllers can be expected to provide the tools necessary to overcome the many problems encountered. As a matter of fact, flexibility is the very first property required from batch process controllers, owing to the great diversity of situations they have to cope with. This conclusion holds even more when dealing with semi-batch processes. There, reactant addition results in unavoidable perturbations of the process, whose balancing is a real challenge. For instance, because the added chemicals are not, in general, at the temperature of the reacting medium, the outset and, to a somehow lower extent, the stop of a feed stage always induce temperature variations. Not only tem‘Present address: Loing, France.
VIDEO-GLASS,
77167
Bagneux
sur
3373
perature thus changes, but it does so with a dynamics having little to do with the dynamics of the chemical reaction. As a consequence, any control procedure faces perturbations whose dynamical characteristics are misleading! This is basically the reason why even adaptiveepredictive controllers must be especially designed, if a keen control of semi-batch reactors is intended for good. In the present work, devoted to semi-batch process control, temperature is, as always, the reference variable. In Section 2, the pilot reactor and the copolymerization reaction, that were used to test different control procedures, are briefly presented. Section 3 describes the algorithms settled together so as to build up an adaptive-predictive technique fitted to the specific conditions issuing from semi-batch operations. Section 4 explains how the technique was actually implemented. Finally, experimental results are discussed in Section 5, where a comparison is also drawn with commercially available controllers. 2. EXPERIMENTAL SECTION
2.1. Apparatus A schematic sketch of the experimental device used all over the present study is given in Fig. 1. The core is a 25 I cylindrical stainless-steel reactor. Its temperature is regulated by means of a fluid circulating through a surrounding jacket. Depending on whether the reactor temperature has to be raised or lowered, this fluid is (i) either heated by a set of two resistors whose electric power P, can be varied at will from 0 to 6 kW, (ii) or it is cooled in a tubular cooler whose cooling rate is changed by varying the external water flow F, from 0 to 32 l/min at most. These two setups are represented on the right part of Fig. 1. A vertical stirrer, set along the reactor axis, is continuously rotated at speed rA (usually 300 rpm) in order to keep the medium as homogeneous as possible in both temperature and composition. Finally, should a semibatch recipe be applied, a specific device, located above the reactor (on the right), allows to supply any additional reactant at a rate F, under control. The
3374
G. DEFAYE et al.
condenser
Fig. I. Schematic sketch of the semi-batch reactor
and its surroundings. On-line measured variables are
encircled once; (adjustable) prescribed parameters are encircled twice.
reactions we are interested in are all carried out at atmospheric pressure. However, in the temperature range under consideration, solvent evaporation always takes place to a significant extent. Accordingly, a water-cooled condenser must be used to flow the solvent back to the reactor, as shown in Fig. 1. Apart from the above-mentioned control parameters, five variables are measured on-line. Three of them are related to the circulating thermofluid (Aow rate F,, inlet temperature T, outlet temperature To), whereas the other two characterize the state of the reacting medium itself: temperature T, and dynamic viscosity t]. This last property is determined by means of an industrial sensor (trademark SOFRASER). It measures how the vibration amplitude of a rod immersed in the medium varies while put in motion at constant energy supply.
2.2. Chemical reaction Among the huge number of chemical reactions that could be involved, we chose a copolymerization reaction between vinyl acetate and 2-ethyl-hexyl acrylate. Like many other free radical chain reactions, this one is initiated by benzoyl peroxide. Its commonest solvent is ethyl acetate, whose boiling temperature is 77°C at 1 atm. This reaction is well fitted to our goal, due to the following characteristics. First, it releases heat at a rather high rate. This is useful for testing the
performance of any control procedure which has to respond quickly enough and with an appropriate amplitude to temperature variations. Second, acryiate has a much higher activity than does vinyl acetate: the difference rises to two orders of magnitude. As a consequence, the first reactant must be progressively added to the reacting medium; otherwise, no copolymer is obtained at all, only polyacrylate is formed. This is, of course, the reason why a semibatch recipe is absolutely required to achieve the desired reaction. And, in turn, this is also why the procedure cannot be entirely set free of temperature variations. Finally, the end product is a glue whose adhesiveness is sensitive to what has actually occurred during the course of the process, especially with respect to the temperature profile. Therefore, beside on-line measurements, this provides a good sensor of the overall control. At the same time, it allows one to figure out the improvements arising from various changes that can be applied to the procedure. Nevertheless, this issue which deserves to be addressed in spite of industrial applications is not taken into consideration hereafter. 2.3. Applied procedure Figure 2 schematically presents the current procedure as it stands, together with the expected reactor temperature and medium viscosity. At the beginning,
Adaptive-predictive
temperature
control
of semi-batch
reactor temperature
= 20min
semi-batch
SECTIONS
The background of the analysis is an autoregressive moving-average discrete model, assumed to provide a relevant description of the process under consideration. It is written as follows:
i=,
a&
- i) + g
bju(k
-j
-
profile
time
ml
vinyl acetate, diluted in one half of the solvent, is put in the reactor at room temperature. In our experiments, this means an initial volume of 10 I at about 20°C (“feed 1”). Then the reactor is heated as quickly as possible (Pw = 6 kW) up to a temperature of 76°C that is to say close to the boiling temperature of the solvent. This heating stage takes about 20 min with our device. Once the running temperature is thus reached, the appropriate amount of benzoyl peroxide (i.e. 0.7 1) is suddenly added. A temperature drop results, as indicated, since the initiator is at 20°C. Next, the reactor is fed at a constant flow rate with the second half of the solvent containing the whole amount of 2-ethyl-hexyl acrylate. This stage (“feed 2”) lasts for 2 h at a flow rate of 5 l/h. The second part of the reactant is also sent at room temperature. After completion of the feed stage, one lets the copolymerization continue at 76°C over a variable length of time t,. In fact, one simply waits until the medium viscosity reaches the desired value of 4 P (i.e. 0.4 Pas). Indeed, it is known that the final viscosity at room temperature will thus be near 20 P, the required value for the glue to be used in industrial pregumming machines. As soon as the condition ~(76°C) = 4 P is fulfilled, the medium is quenched, so as to stop immediately the polymerization.
y(k) = z
3375
)
- 2 hours Fig. 2. Standard
3. ALGORITHM
reactors
r) + e(k)
i=1 (1)
procedure.
where k is the discretized time (t = kAt, At is the sampling time interval), y(k) the output variable, u(k) the input (or control) variable, e(k) the model error, r the process time delay, n,, nb the polynomial orders (with respect to input and output) and ai, bj are model parameters. Using the time delay operator q-l, an alternative form of eq. (1) is A(q-‘)y(k)
= B(q-‘)q-‘u(k)
+ e(k)
(2)
where ~(~-1)
= 1-
2 i=l
aiq-i
and
B(q-‘)
= 5
h,q-j.
j=l
It is the calculus of parameters ai and bj at each step k that makes the above model adaptive. In order to be
efficient, the computation algorithm must be able to follow to the nearest point any significant output variation, At the same tine, however, this algorithm also has to remain insensitive to noise. These two requirements are obviously opposite, and a compromise between them has to be settled, which largely depends upon the noise level. We were thus led to couple two estimation procedures with a predictive control method, robust enough to deal with a nonminimum phase, possibly unstable process. We begin with a description of this method. 3.1. Control algorithm It belongs to a category named model-based predictive controls (MBPC), recently reviewed by Garcia et al. (1989) and De Keyser (1990). For the sake of clarity, it is useful to recall briefly the four basic ingredients shared by MBPC methods. (i) A transfer function model is used to predict the output values obtained over a time interval m = 1,
G.
3376
DEFAYE
. . . , H, t r (H, is the prediction horizon), according to y(k + m/k) =
5 aiy(k
+ 5
+ m - i/k)
bjU(k+ffl-j-r)+c(k+m/k).
(3)
j=l
By convention, the standard notation y(k + m/k) means: value of y at time k + m computed from data available at time k. An error term e(k + m/k) must be included in eq. (3), since the model can never account perfectly well for the process. On the other hand, this term is not really computable in advance. So, one has to set it equal to some prescribed value. The error observed at time k is, in general, used to this end, and it is applied all over the time window m = 1, . , H, + r so that e(k + m/k)
= e,(k)
-
= y,,(k) -
5 j=
b+(k
-j
5 a,y,(k i=1
- i)
- r)
1
(subscript p refers to the quantities actually observed from the real process, as opposed to those predicted). (ii) A reference trajectory, labelled s,(k + m), is then introduced. It represents the desired process output, which starts from y(k + r/k) and must end at the setpoint w. Of course, there are many possible trajectories matching these two conditions. A certain knowhow is helpful in choosing a trajectory having the right dynamical properties. Otherwise, the simplest solution consists in taking the first-order trajectory: s,(k + m) = a&k
+ m -
1) + (1 - a,,,)w(k
+ m).
(iii) In order to determine the control variables u that have to be applied, one minimizes the following function: tf,+r J(k)
=
x
[s,(k
+ m) -
y(k + m/k)12
y
(iv) Only the first value u(k) thus determined is actually applied to the process, the whole procedure being iterated at each step. The algorithm that we used (Defaye et al., 1983) involved a particular strategy to evaluate the predicted output y(k + m/k). Details on this calculus, which is slightly different from the standard k-step ahead predictor (Clarke et al., 1987), are given in the appendix. 3.2. Estimation algorithm Identification methods are numerous. A good, fairly recent analysis of these methods is due to Ljung (1987). One of the most popular techniques for updating parameters is presumably the so-called recursive least-squares (RLS) algorithm. It is this kind of method that we used to calculate parameters ai and bj at each step k. However, rather than taking the direct form [eq. (2)]-or, as well, its modified expression according to the CARlMA model (Tuffs and Clarke, 1985)-we found it more appropriate to apply an incremental form of this equation, namely A(q-‘)Ay(k)
= E(q-‘)q-‘Au(k)
+ Ae(k)
(5)
where the symbol A represents 1 - 9-l. In doing so, the origin of the data is shifted at the current (steady) state and troubles issuing from drifts are overcome. As mentioned above, the main problem to be tackled is that of semi-batch-induced dynamical perturbations. Because they develop over their own timescale, their effects can be significantly reduced by a suitable filtering of input and output signals. Furthermore, all measured variables are simultaneously washed out, by serendipity, of their high-frequency noise. The incidence of such a filtering procedure has been discussed by Ljung (1987) and also by Inglis et al. (1991) in a more recent work. In the present study, a second-order filter, whose transfer function is given by
m=Hl
+
et al.
z
m=,
ru(k + m -
1) - u(k + m -
2)-j’.
F(q_‘) (4)
The first sum in the r.h.s. of eq. (4) represents the “distance” (to the square) between the reference trajectory and the predicted output computed at time k according to eq. (3). This distance is integrated over a certain time window which opens at time If,. This parameter can be varied between the two bmits: r < N, < H, + r. The second term in eq. (4) takes into consideration the amplitude variations (also to the square) of the input, up to a time limit H, named control horizon. Necessarily, one has 1 < H, -c H,. Beyond H,, nothing is claimed about the changes that could be imposed on u. Because the function J(k) mixes together two different kinds of “deviation”, it is appropriate to have a means of balancing the two contributions. This is achieved through the control weighting parameter y.
=
0.49 + 0.259-l 1 - 0.49-I
+0.149-z
appeared to be the most appropriate with respect to the various time-scale and delay involved in the experimental procedure. Starting with an input-output vector x(k) and a vector of model parameters 8, x=(k) = [Ay(k -
l), . . . ,Ay(k - n.)
Au@ - I - 1), . . . , Au(k - r - nb)] @ = Cal,.
. aran,, bl,. . . , h.,]
the RLS algorithm writes down in standard notations: B(k) = 8(k -
1) + L(k)&(k)
e(k) = By(k) -
x*(k)B(k
-
1).
The gain vector L(k) = P(k - l)x(k)/[A(k)
+ xT(k)P(k
-
1)x(k)]
Adaptive-predictive temperature control of semi-batch reactors and the covariance matrix P(k) = [P(k -
1) - T@)x=(&)P(k - 1)1/R(k)
both involve a forgetting factor 1. This parameter, which is less than unity, is chosen at will. It allows to balance the data taken into consideration with respect to their “age”. Indeed, it is appropriate to forget progressively the former states of the system at an adjustable rate, characterized by a time constant = equal to i/(1 - A). It was soon realized that a difficulty has to be overcome when applying to the RLS algorithm to filtered data in the form presented above. As a matter of fact, the calculus must be modified in some way if a high enough sensitivity to process variations is to be preserved. Several variants have been proposed over the past years, a critical survey of which has been published by Shah and Cluett (1991). In the present work, we brought into play two approaches aimed at improving the dynamics: (i) According to Fortescue et al. (1981), the parameter I is changed at each step in such a way that a link is established between the amount of forgotten data and the new information supplied by the last measurement. To this end, the new value of 1 at time k + 1 is selected by I(k + 1) = max {l - Cl - xT(k)L(k)]s(k)‘/(rO,
A,,,,,}
where G,, and &,, are two adjustable parameters. Aminrepresents a lower limit that parameter n must not cross (otherwise, the covariance matrix would diverge). On the other hand, the smaller the no the larger is the sensitivity of the algorithm. (ii) As suggested earlier by Irving et al. (1979), one can also give the trace of matrix P a fixed value /I by varying jl properly. In this case, the gain factor is determined without data balancing: L(k) = P(k -
l)x(k)/[l
+ x’(k)P(k
-
1)x(k)].
Then I(k) is calculated so as to satisfy the condition tr P(k) = tr P(k -
1) = p.
Of course, the larger the B the more sensitive is the corresponding estimator. Finally, since it is necessary to have a positivedefinite covariance matrix, the two variants are improved by including a U-D factorization technique (Bierman, 1977) in the computation of P(k).
4. CONTROL
1MPLEMENTING
is the heat release by the chemical reaction which is responsible for the main nonlinearity of the process, due to Arrhenius’s law (exponential dependence of rate constants upon temperature). Viscosity, on the other hand, is rather a “slave” property. Uneasy to manage at will, it has almost no direct influence on the reaction rate, despite its effects on stirring efficiency and, also, on the heat exchange coefficient (Nusselt number). The only relevant control variable one can It
3317
act upon is thus the temperature. The adaptive control algorithms described in the previous section are aimed to cope with this situation and to some other features of semi-batch reactors (e.g. the heat removal capacity changes over the time since the reacting volume doubles in the present case). The goal one must reach is twofold: (i) maintain the medium temperature at a fixed constant value once the setpoint is established and (ii) reduce as much as possible the amplitude of the unavoidable temperature variations resulting from semi-batch operations (especially initiator addition). To this end, one has to vary, properly and in due time, the temperature Ti of the thermofluid. We did it in the following way. The cooling water flow F, is usually set at a fixed mediumrange value (namely 15 l/min), and the thermofluid is warmed up to the desired temperature by the two resistors. The applied electric power P, is thus the control parameter actually managed. However, when the medium temperature T, exceeds the setpoint by too large a difference (say l”C), the resistors are switched off (P, = 0) and the flow F, is raised to its maximum (i.e. 32 l/min). As soon as the temperature T, begins to diminish, the “standard” control conditions (i.e. F, = 15 I/min, P, variable) are restored. In practice, the transfer function model, assumed to represent how the temperature T, depends upon the delivered electric power P,, is the key ingredient of the control. To build up successfully this function during the estimation phase, data must be filtered as indicated above. Moreover, specific care has to be taken in order to eliminate bias in estimated parameters issuing from the presence of the second feed stage. As a matter of fact, the beginning and end of that stage correspond to “extraordinary” events, which will not occur twice. So, one has to cancel the spurious correlation which is established at these times by least-squares computation (Inglis et al., 1991). The time characteristics of our experimental device with respect to T, variations (reactor containing 10 I of solvent, open-loop condition) have the following order of magnitude: time constant _ IO3 s heating/cooling dead time - 180 s. Owing to this, we
used
our algorithms:
-a sampling period of 120 s for parameter estimation (this period may seem rather long; it was chosen essentially according to the incremental model [eq. (S)], which sometimes yields numerical problems when smaller periods are taken); ~ a control performing period also equal to 120 S; -a control data recording period of 20 s for better control watching. 5.
RESULTS
AND DISCUSSION
The main results of this study are summarized in Figs 3 and 4, which represent the variations vs time of three quantities: measured temperature T,, electric
G. DEFAYE
3378
power P, represented as a percentage of 6 kW, cooling water flow F,. Moreover, Fig. 6 displays similar data obtained with commercially available adaptive controllers. In all cases, the applied experimental procedure was that depicted in Section 2.3.. What changes from one experiment to another is uniquely the way in which control is carried out. Consequently, comparison is really meaningful. Figure 3(a) was obtained under “control-free” conditions with respect to T,. In other words, the temperature T, was solely recorded, but not involved in any control algorithm. The heating/cooling setup was used to maintain the temperature 7;: at 76°C independently of what happened in the reactor. Now, during feed 2, the reaction rate was very large and, as a consequence, so was the rate of heat release. A point was T, became higher even reached where the temperature than the boiling temperature of the solvent (i.e. 77°C). Nevertheless, due to solvent evaporation on one hand, and to flow back to the reactor after solvent condensation, the rise in temperature remained slow enough. Even at the end of feed 2, T, did not exceed 79°C. Then it decreased, as did the reaction rate. In Fig. 3(b) a predictive control algorithm was employed. To this end, a transfer function model was determined by open-loop identification with the following experimental conditions: - reactor containing 10 1 of solvent -FF, = 15 l/min -P, = 5 kW + a pseudo-random binary quence superimposed at 600 s time intervals. This led to the first-order T,(k) = O.S9T,(k -
estimation
se-
model
1) + O.O57P,(k
-
2)
then applied once and for all. The controller tuning parameters yielding the best result [i.e. the minimum value of J(k), see eq. (4)] were H, = 5,
H, = 1,
y = 0.06,
a, = 0.74
et al.
corresponding to a time constant of 400 s. In some sense, this experiment provides us with a reference point. As far as the viscosity of the thermofluid is temperature-dependent, as it actually is, the flow rate Ff will vary depending on whether the fluid is warmed or cooled. Accordingly, one must also expect variations in the time delay. This fact has to be taken into account if an adaptive-predictive control procedure is to be applied. To model these unavoidable variations, we found it appropriate to take a third-order polynomial B(q-‘). For the two estimation algorithms described in Section 3.2 (namely Fortescue, and constant trace), we thus referred to the same kind of model given by AT,(k)
= a,AT,(k
1) + blAP,(k
+ bt AP,(k
Fortescue: Constant
&in = 0.8, trace:
-
1)
- 2) + b3 AP,(k
The selected tuning parameters algorithms were, respectively,
of these
- 3). estimation
(r0 = 1
p = 0.04.
Furthermore, the same set of control parameters was used for both algorithms, so as to allow for a fully significant comparison of their performance. As may be seen, the Fortescue estimation [Fig. 4(a)] undoubtedly provided better results than did the constant trace [Fig. 4(b)]. It was obviously more effective, not only in lowering the temperature rise at the beginning of feed 2 but also in damping oscillations around the setpoint. As already mentioned, the input and output data that served to perform identification in RLS algorithms were filtered. The applied second-order filter has a damping factor of 0.7; it washes out perturbations whose time constant is shorter than 540s. Thanks to this, the adaptive control is insensitive to perturbations which should remain unmodelled, in
80
ä
4 78
-
feed2
78
76
20 0 0
50
loo
150
200
250
(a> Fig. 3. Observed
lemperature
300
350
Time (min)
400
66! 0
! 50
103
150
200 (b)
and applied heating/cooling parameters as functions free” condition (fixed Ti); (b) predictive control.
250
300
350
Time (min) of time: (a) “control-
400
0
Adaptive-predictive
temperature
control of semi-batch
reactors
3379
80
78 T@)
76
4
74 -
72
loo
WfJo Pw(%) +* 60
70
- 40
-40
68
II
Fwo/min) *-
66 +
F,#W
20
r0
I-0 50
0
loo
150
200 (a)
250
300
350
400
. . 20
0
50
loo
150
Time (min)
200
250
(b)
300 Tie
350
400
(min)
Fig. 4. Observed temperature and applied heating/cooling parameters as functions of time: (a) Fortescue estimation algorithm; (b) constant trace estimation algorithm.
particular
those
occurring
at the beginning
and at the
end of feed 2. Indeed, the corresponding variations have time constants well below that of thermal events, whereas one has to determine the transfer function linking P, and r,, and nothing else. The effect of data filtering on the model time constant is shown in Fig. 5 in the case of the Fortescue algorithm. As may be seen, a drastic drop in the time constant is observed right after initiator addition and the beginning of feed 2, when data are not filtered. This drop is undoubtedly irrelevant: indeed the time constant is three times smaller than that calculated for a reactor working in open loop. The slight modification resulting from adding initiator is not at all able to account for such a large ratio. To the contrary, no trouble of this sort is met when filtering of data is achieved. This result points out the usefulness and even the need of data filtering at least in the presence of large external perturbations. Trademark adaptive controllers being readily available, it is interesting to look at their control performances when driving a semi-batch reactor. Two such controllers were tested: the Foxboro EXACT and the A.B.B. Novatune. They yielded the results shown in Fig. 6. The Foxboro EXACT (Kraus and Myron, 1984) is basically a discretized version of an analog PID. Owing to this, its sampling period has to be small (0.1 s). To set the initial values of its parameters, we filled the reactor with 10 1 of solvent, as usual, and then we let the controller analyze the transient response of the closed-loop system to setpoint changes. Temperature data filtering was performed by a second-order Butterworth filter, with a 30 s time constant. The Foxboro system is designed so that, during a run, parameter adaptation takes place whenever the difference between the measured temperature (T,) and the setpoint exceeds twice the (allowed) noise amplitude. In the experiment presented in Fig. 6(a), this amplitude was assigned to the lowest adjustable value, namely 0.5”C.
soo.
Fig. 5. Effect of data filtering on the model time constant (Fortescue estimation algorithm).
The nature of the A.B.B. Novatune system (Bengtsson and Egardt, 1984) is much closer to our approach. Indeed, it involves a RLS estimation algorithm (forgetting factor k = 0.98) together with a minimum variance controller. It relies upon the following model: T,(k + H,)
- PLT,(k
+ H, + 1) - (1 - PL)T,(k)
= A(q-‘)AT,(k)
+ B(q-‘)AP,(k)
where PL is the desired pole location in closed loop. After various trials, we selected the following conditions: = 3 (with a 60 s sampling time interval) = 0.86 (so that the model time constant still equal to 400 s) A (q - ‘) = second-order polynomial B(q- ‘) = third-order polynomial.
HC PL
is
It is worth noting that, here too, parameter adaptation is not active at all times, but only when variations
G. DEFAYE er al.
3380 80 78
TR(*c) t
76
76
74 100 72
80 Py&w )
70
60 40
68
20
FwWmin).
0
66 0
50
100
150
200 (a)
250
300
350
0
400
SO
loo
150
UK)
Time (min)
250
300
350
400
Time (min)
@I
Fig. 6. Observed temperature and applied heating/cooling parameters as functions of time: (a) Foxboro EXACT adaptive controller; (b) A.B.B. Novatune adaptive controller.
e control - constant trace
0
50
100
150
200
250
300
350
400
Time @in) Fig. 7. Viscosity vs time plots for the various tested controllers
in either temperature T, or electric power P, exceed an adjustable threshold (respectively, 0.4”C and 2%). Even at a glance at Figs 4 and 6 the improvement brought out by our adaptive controllers is indisputable. Perhaps the most striking feature is the trend towards oscillation exhibited by the process. This trend is clearly emphasized by the last two experiments, and was already visible before. Whatever the criterion taken into consideration, the superiority of the Fortescue algorithm is prominent. It leads to smaller temperature differences with the setpoint and to fewer changes in electric power supply. It also provides a much more efficient damping of oscillations compared to any other algorithm. Although the differences might perhaps seem relatively weak, they are still of consequence, owing to the sensitivity of the end product to temperature vari-
ations having occurred during the synthesis. An overall indication, supporting this conclusion, appears in Fig. 7. There, one can see how the viscosity q has evolved during the course of the reaction. The rapid increase at the end of each curve corresponds simply to the outset of quenching and reflects the final temperature drop. As expected, the most regular increase in viscosity has been obtained with the Fortescue adaptive controller. Curve 3 has no bump at all; its slope always changes in a progressive manner, indicating that the copolymerization proceeds smoothly. It even continues to do so beyond the desired threshold of 4 P. Furthermore, except in the case of the A.B.B. controller where something peculiar has occurred presumably due to the chemicals, the reaction is also more rapid and the targetted viscosity is reached sooner.
:e control of semi-batch Adaptive-predictivetemperatul
6. CONCLUSION this work, we put into shape an adaptive-predictive controller combining the following: In
(i) An estimation algorithm based on RLS techniques; a thorough investigation led us to select two variants of this type, one having a variable forgetting factor (Fortescue) and the other preserving the trace of the covariance matrix (constant trace); (ii) A stage of data filtering aimed at eliminating spurious dynamical effects, in particular those linked to reactants addition(s). In this way, we were able to control successfully the temperature of a reactor where a copolymerization reaction took place despite the perturbations due to initiator addition and to secondary feed. Particularly worth of note is the complete
vanishing
of
temperature
(Fortescue) that commercial failed to achieve.
adaptive
e
FJF. FW HC
HP k
L nn, nb PW P
parameters defined in eq. (1) model error flow rate of thermofluid, l/h flow rate of feed 2, I/h flow rate of cold water, l/min control horizon prediction horizon discretized time gain vector polynomial orders electric power, kW covariance matrix
time delay stirring rate, rpm reference trajectory thermofluid inlet temperature, “C thermofluid outlet temperature, “C reacting medium temperature, “C model input setpoint input-output vector model output
r r.4 b
Ti T, T, u W
x Y
Greek letters
control weighting factor forgetting factor reaction medium viscosity, P vector of model parameters
: I
controllers
Acknowledgements-This research was granted by the CNRS (ARC “Optimal Control and Automation of Batch Reactors”) and by the Regional Council of Aquitaine. The authors would like to acknowledge assistance and fruitful discussions with their colleagues from “Institut du Pin” (Bordeaux University).
ai, bj
3381
oscillations
Seeking for the main source of the observed improvement, we came to the conclusion that data filtering was crucial. Indeed, the cut-off of the dynamical perturbations which result from some sort of “foreign event”, is essential in preserving the right set of estimated parameters. Otherwise, troubles ensue. As is wellknown cut-off is all the more effective when the timescales involved are well separated. Obviously, the process we tested fulfilled this condition. It must be emphasized, however, that the chosen example is in no way peculiar: it is a mere archetype of semi-batch processes. Furthermore, if the time-scale separation be missing, one could still expect this kind of adaptive-predictive controller to be able to do well, thanks to its intrinsic flexibility. As a matter of fact, spurious dynamical perturbations developing on a time-scale comparable to that of the process will no longer induce errors in parameters estimation. Accordingly, a wide range of applications is open, in principle, to this kind of controller, at the expense of an appropriate fitting to the process dealt with.
NOTATION
reactors
REFERENCES Bengtsson, G. and Egardt, B., 1984,Experiences with selftuning control in the process industry. Preprint 9th fFAC World Congress, Budapest, Vol. XI, pp. 132-140. Bierman, G. I., 1977, Factorization Methods for Discrete Sequential Esfimation. Academic Press, New York. Clarke, D. W., Mohtadi, C. and Tuffs, P. S., 1987, Generalized predictive control. Automatica 23, 137-160. Cluett, W. R., Shah, S. L. and Fisher, D. G., 1985, Adaptive control of a batch reactor. Chem. Engng Commun. 38, 67-78. Cott, B. S. and Macchietto, S., 1989, Temperature control of exothermic batch reactors using generic model control. Ind. Engng Cbem. Res. UI, 1177-I 184. Defaye, G., Caraip, L. and Jouve, P., 1983, A simple deterministic predictive control algorithm and its application to an industrial chemical process: a distillation column. Chem. Engng J. 27, 161-166. De Keyser, R. M. C., 1990, Model based predictive control toolbox. Proceedings of CIM Europe Workshop on Computer Integrated Design of ControHed Industrial Systems, Paris, pp. l-56. Fortescue. J. R., Kershenbaum, L. S. and Ydstie, B. E., 1981, Implementation of self-tuning regulators with variable forgetting factor. Automatica 17, 831-835. Garcia, C. E., Prett, D. M. and Morari, M., 1989, Model predictive control: theory and practice-a survey. Automatica 25, 335-348. Inglis, M. P., Cluett, W. R. and Penlidis, A., 1991, Long range predictive control of a polymerization reactor. Can. J. them. Engng 69, 120-129. Irving, E., Barret. J. P., Charcossey, C. and Monville, J. P., 1979, Improving power network stability and unit stress with adaptive generator control. Automatica IS, 31-46. Juba M. R. and Hamer, J. W., 1986, Progress and challenges in batch process control. Proceedings of 3rd Inrernational Conference on Chemical Process Control. Elsevier.New
York.
Kraus, T. W. and Myron, T. J.. 1984, Self-tuning PID controller uses pattern recognition approach. Contr. Engng 6, 106-11 I. Lee, P. L. and Sullivan, G. R., 1988, Generic model control. Comput. them. Engng 12, 573-598. Liptak, B. C., 1986. Controlling and optimizing chemical reactors. Chem. Engng 24, 69-81. Ljung, L., 1987. System Identification: Theory for the User. Prentice-Hall, Englewood ClilTs, NJ. Shah, S. L. and Cluett, W. R., 1991, Recursive least-squares based estimation schemes for self-tuning control. Can J. them. Emma 69.89-96. Tuffs, P.-S.“a&C$larke, D. W., 1985, Self-tuning control of offset: a unified approach. Proc. IEE, Part D 132, pp. _ 100-110.
G. DEFAYE et al.
3382 APPENDIX
We briefly summarize hereafter the way in which we determine the input u(k) that is applied at each step k. The complete derivation of these equations may be found elsewhere (Defaye et al.. 1983). The output values, which are predicted according to eq. (3). gather three contributions which may be written as y(k + m/k)
Ear y,l(k +
m/k)1
+ effect of future inputs
Car YAk +
m/k)1
+ predicted error
[or e(k + m/k)].
= effect of past inputs
Identification
+ i condition
form=I,...,r,
notation.
y=y@+Cu+e with
+ m - ilk)
b&k
+ m - j - r) + e(k)
when summing
over j:
“=a, n=nb-m+r
y,,(k + m/k) = 5 c,u(k + m - i) i= 1
+ H, + r/k)]
II’ = [u(k). . . , u(k + H,, - I)]. C is the lower triangular matrix whose elements Cij are Cii = q+
The input vector u corresponding given by a = (C’C + M)-‘{C’(s,
-j_
to the minimum
of J is
- yri - e) + pu(k - I)}
where
P’=cY.4..~701
c, =b, i-L
. ,nb, Cl= bi + c ca,Ck-,+ b, j-1
fori=n~+l,....H,+r
t ci= x a~c~-~ with /=min(i-
I
sT=[s,(k+r+l),...,s,(k+r+H,)]
the coefficients c, being given by
i=*
. . .y(k + HP + r/k)]
+ r + l/k), . . . ,y,&
e7=[Ce(k+r+l/k),....e(k+H,+r/k)]
(of course, no other contribution of the inputs has to be considered further). The yn are written as follows:
for i = 2,.
H, = H,.
it then becomes
y’ = [y(k + r + I/k),
form=r+I,....r+nb-1,
I,
In vector
yz = b,,(k
j-l
with the twofold
I;
H,=r+
of the y+ yields
I. y,f(k + m/k) = E atyAk 1-1
for i=
adaptiveness comes out from the calculus of the cis at each sampling time k. The purpose is to determine the future input giving the function J(k) [see eq. (4)] its minimum value. Being more specific now, we ascribe to parameters H1 and H, the particular values
1.n.).
In a nonadaptive control algorithm, these coefficients e, are calculated once and for all at the outset. On the other hand,
M=
;?r --y
-Y
0
0
0
2y
--y
0
0
0
--Y
2Y
0
0
. 0
0
0
2Y
-Y
0
0
0
-Y
Y
y being the control weighting parameter which balances the two types of deviations appearing in the criterion J(k) [see =I. (411.
1993
ooo9-2509/93 66.00 + 0.00 Q 1993 Pergamon Press Ltd
ADAPTIVE-PREDICTIVE TEMPERATURE SEMI-BATCH REACTORS G. DEFAYE. Equipe
I.A.G.P.,
N. REGNIER,
URA 503 CNRS, Universite (Received
J. CHABANON,’
CONTROL
L. CARALP
and C. VIDAL
EC& Nationale Sup&ieure de Chimie et Physique Bordeaux I, 33405 Talence Cedex, France
18 Noumber
1992; acceptedfor
publication
OF
de Bordeaux,
30 April 1993)
Abstract-Although reactant addition induces unavoidable temperature changes in semi-batch reactors, it is shown that adaptive-predictive techniques may offer a means of balancing them successfully. To this end, appropriate data filtering must be integrated to the estimation algorithm, so that parameters remain insensitive to the nonrelevant dynamical perturbations generated by semi-batch operations. A pilot reactor is used to illustrate imorovements thus workable while carrying out a moderately exothermic copolymerization reaction.
I. INTRODUCTION
The control of batch reactors is more difficult to achieve than that of continuous processes (Juba and Hamer, 1986) due to the lack of steady state. Indeed, the dynamics of a batch process is essentially timedependent. PID regulators, which are well fitted to maintain a fixed setpoint, often fail to control the batch processes properly. In fact, such a task requires controllers able to face, without adjustment, a wide range of working conditions. Liptak (1986) showed how to improve the performances of PlD regulators, in particular during the initial phase of a batch process. However, the so-called “batch unit” and “dual mode” are nothing much more than a palliative. Better results can be obtained by simply giving up PIDs. Cluett et al. (1985) thus developed a self-tuning algorithm based on a transfer function model which connects the reactor temperature to the flow rate of a circulating thermofluid. Cott and Macchietto (1989) designed another procedure, involving the on-line estimation of the rate at which heat is released by the ongoing reaction. The controller may even include nonlinear models (Lee and Sullivan, 1988). Besides the improvements set up in this way, the prominent outcome of this kind of study is the following conclusion: nothing but adaptive predictive controllers can be expected to provide the tools necessary to overcome the many problems encountered. As a matter of fact, flexibility is the very first property required from batch process controllers, owing to the great diversity of situations they have to cope with. This conclusion holds even more when dealing with semi-batch processes. There, reactant addition results in unavoidable perturbations of the process, whose balancing is a real challenge. For instance, because the added chemicals are not, in general, at the temperature of the reacting medium, the outset and, to a somehow lower extent, the stop of a feed stage always induce temperature variations. Not only tem‘Present address: Loing, France.
VIDEO-GLASS,
77167
Bagneux
sur
3373
perature thus changes, but it does so with a dynamics having little to do with the dynamics of the chemical reaction. As a consequence, any control procedure faces perturbations whose dynamical characteristics are misleading! This is basically the reason why even adaptiveepredictive controllers must be especially designed, if a keen control of semi-batch reactors is intended for good. In the present work, devoted to semi-batch process control, temperature is, as always, the reference variable. In Section 2, the pilot reactor and the copolymerization reaction, that were used to test different control procedures, are briefly presented. Section 3 describes the algorithms settled together so as to build up an adaptive-predictive technique fitted to the specific conditions issuing from semi-batch operations. Section 4 explains how the technique was actually implemented. Finally, experimental results are discussed in Section 5, where a comparison is also drawn with commercially available controllers. 2. EXPERIMENTAL SECTION
2.1. Apparatus A schematic sketch of the experimental device used all over the present study is given in Fig. 1. The core is a 25 I cylindrical stainless-steel reactor. Its temperature is regulated by means of a fluid circulating through a surrounding jacket. Depending on whether the reactor temperature has to be raised or lowered, this fluid is (i) either heated by a set of two resistors whose electric power P, can be varied at will from 0 to 6 kW, (ii) or it is cooled in a tubular cooler whose cooling rate is changed by varying the external water flow F, from 0 to 32 l/min at most. These two setups are represented on the right part of Fig. 1. A vertical stirrer, set along the reactor axis, is continuously rotated at speed rA (usually 300 rpm) in order to keep the medium as homogeneous as possible in both temperature and composition. Finally, should a semibatch recipe be applied, a specific device, located above the reactor (on the right), allows to supply any additional reactant at a rate F, under control. The
3374
G. DEFAYE et al.
condenser
Fig. I. Schematic sketch of the semi-batch reactor
and its surroundings. On-line measured variables are
encircled once; (adjustable) prescribed parameters are encircled twice.
reactions we are interested in are all carried out at atmospheric pressure. However, in the temperature range under consideration, solvent evaporation always takes place to a significant extent. Accordingly, a water-cooled condenser must be used to flow the solvent back to the reactor, as shown in Fig. 1. Apart from the above-mentioned control parameters, five variables are measured on-line. Three of them are related to the circulating thermofluid (Aow rate F,, inlet temperature T, outlet temperature To), whereas the other two characterize the state of the reacting medium itself: temperature T, and dynamic viscosity t]. This last property is determined by means of an industrial sensor (trademark SOFRASER). It measures how the vibration amplitude of a rod immersed in the medium varies while put in motion at constant energy supply.
2.2. Chemical reaction Among the huge number of chemical reactions that could be involved, we chose a copolymerization reaction between vinyl acetate and 2-ethyl-hexyl acrylate. Like many other free radical chain reactions, this one is initiated by benzoyl peroxide. Its commonest solvent is ethyl acetate, whose boiling temperature is 77°C at 1 atm. This reaction is well fitted to our goal, due to the following characteristics. First, it releases heat at a rather high rate. This is useful for testing the
performance of any control procedure which has to respond quickly enough and with an appropriate amplitude to temperature variations. Second, acryiate has a much higher activity than does vinyl acetate: the difference rises to two orders of magnitude. As a consequence, the first reactant must be progressively added to the reacting medium; otherwise, no copolymer is obtained at all, only polyacrylate is formed. This is, of course, the reason why a semibatch recipe is absolutely required to achieve the desired reaction. And, in turn, this is also why the procedure cannot be entirely set free of temperature variations. Finally, the end product is a glue whose adhesiveness is sensitive to what has actually occurred during the course of the process, especially with respect to the temperature profile. Therefore, beside on-line measurements, this provides a good sensor of the overall control. At the same time, it allows one to figure out the improvements arising from various changes that can be applied to the procedure. Nevertheless, this issue which deserves to be addressed in spite of industrial applications is not taken into consideration hereafter. 2.3. Applied procedure Figure 2 schematically presents the current procedure as it stands, together with the expected reactor temperature and medium viscosity. At the beginning,
Adaptive-predictive
temperature
control
of semi-batch
reactor temperature
= 20min
semi-batch
SECTIONS
The background of the analysis is an autoregressive moving-average discrete model, assumed to provide a relevant description of the process under consideration. It is written as follows:
i=,
a&
- i) + g
bju(k
-j
-
profile
time
ml
vinyl acetate, diluted in one half of the solvent, is put in the reactor at room temperature. In our experiments, this means an initial volume of 10 I at about 20°C (“feed 1”). Then the reactor is heated as quickly as possible (Pw = 6 kW) up to a temperature of 76°C that is to say close to the boiling temperature of the solvent. This heating stage takes about 20 min with our device. Once the running temperature is thus reached, the appropriate amount of benzoyl peroxide (i.e. 0.7 1) is suddenly added. A temperature drop results, as indicated, since the initiator is at 20°C. Next, the reactor is fed at a constant flow rate with the second half of the solvent containing the whole amount of 2-ethyl-hexyl acrylate. This stage (“feed 2”) lasts for 2 h at a flow rate of 5 l/h. The second part of the reactant is also sent at room temperature. After completion of the feed stage, one lets the copolymerization continue at 76°C over a variable length of time t,. In fact, one simply waits until the medium viscosity reaches the desired value of 4 P (i.e. 0.4 Pas). Indeed, it is known that the final viscosity at room temperature will thus be near 20 P, the required value for the glue to be used in industrial pregumming machines. As soon as the condition ~(76°C) = 4 P is fulfilled, the medium is quenched, so as to stop immediately the polymerization.
y(k) = z
3375
)
- 2 hours Fig. 2. Standard
3. ALGORITHM
reactors
r) + e(k)
i=1 (1)
procedure.
where k is the discretized time (t = kAt, At is the sampling time interval), y(k) the output variable, u(k) the input (or control) variable, e(k) the model error, r the process time delay, n,, nb the polynomial orders (with respect to input and output) and ai, bj are model parameters. Using the time delay operator q-l, an alternative form of eq. (1) is A(q-‘)y(k)
= B(q-‘)q-‘u(k)
+ e(k)
(2)
where ~(~-1)
= 1-
2 i=l
aiq-i
and
B(q-‘)
= 5
h,q-j.
j=l
It is the calculus of parameters ai and bj at each step k that makes the above model adaptive. In order to be
efficient, the computation algorithm must be able to follow to the nearest point any significant output variation, At the same tine, however, this algorithm also has to remain insensitive to noise. These two requirements are obviously opposite, and a compromise between them has to be settled, which largely depends upon the noise level. We were thus led to couple two estimation procedures with a predictive control method, robust enough to deal with a nonminimum phase, possibly unstable process. We begin with a description of this method. 3.1. Control algorithm It belongs to a category named model-based predictive controls (MBPC), recently reviewed by Garcia et al. (1989) and De Keyser (1990). For the sake of clarity, it is useful to recall briefly the four basic ingredients shared by MBPC methods. (i) A transfer function model is used to predict the output values obtained over a time interval m = 1,
G.
3376
DEFAYE
. . . , H, t r (H, is the prediction horizon), according to y(k + m/k) =
5 aiy(k
+ 5
+ m - i/k)
bjU(k+ffl-j-r)+c(k+m/k).
(3)
j=l
By convention, the standard notation y(k + m/k) means: value of y at time k + m computed from data available at time k. An error term e(k + m/k) must be included in eq. (3), since the model can never account perfectly well for the process. On the other hand, this term is not really computable in advance. So, one has to set it equal to some prescribed value. The error observed at time k is, in general, used to this end, and it is applied all over the time window m = 1, . , H, + r so that e(k + m/k)
= e,(k)
-
= y,,(k) -
5 j=
b+(k
-j
5 a,y,(k i=1
- i)
- r)
1
(subscript p refers to the quantities actually observed from the real process, as opposed to those predicted). (ii) A reference trajectory, labelled s,(k + m), is then introduced. It represents the desired process output, which starts from y(k + r/k) and must end at the setpoint w. Of course, there are many possible trajectories matching these two conditions. A certain knowhow is helpful in choosing a trajectory having the right dynamical properties. Otherwise, the simplest solution consists in taking the first-order trajectory: s,(k + m) = a&k
+ m -
1) + (1 - a,,,)w(k
+ m).
(iii) In order to determine the control variables u that have to be applied, one minimizes the following function: tf,+r J(k)
=
x
[s,(k
+ m) -
y(k + m/k)12
y
(iv) Only the first value u(k) thus determined is actually applied to the process, the whole procedure being iterated at each step. The algorithm that we used (Defaye et al., 1983) involved a particular strategy to evaluate the predicted output y(k + m/k). Details on this calculus, which is slightly different from the standard k-step ahead predictor (Clarke et al., 1987), are given in the appendix. 3.2. Estimation algorithm Identification methods are numerous. A good, fairly recent analysis of these methods is due to Ljung (1987). One of the most popular techniques for updating parameters is presumably the so-called recursive least-squares (RLS) algorithm. It is this kind of method that we used to calculate parameters ai and bj at each step k. However, rather than taking the direct form [eq. (2)]-or, as well, its modified expression according to the CARlMA model (Tuffs and Clarke, 1985)-we found it more appropriate to apply an incremental form of this equation, namely A(q-‘)Ay(k)
= E(q-‘)q-‘Au(k)
+ Ae(k)
(5)
where the symbol A represents 1 - 9-l. In doing so, the origin of the data is shifted at the current (steady) state and troubles issuing from drifts are overcome. As mentioned above, the main problem to be tackled is that of semi-batch-induced dynamical perturbations. Because they develop over their own timescale, their effects can be significantly reduced by a suitable filtering of input and output signals. Furthermore, all measured variables are simultaneously washed out, by serendipity, of their high-frequency noise. The incidence of such a filtering procedure has been discussed by Ljung (1987) and also by Inglis et al. (1991) in a more recent work. In the present study, a second-order filter, whose transfer function is given by
m=Hl
+
et al.
z
m=,
ru(k + m -
1) - u(k + m -
2)-j’.
F(q_‘) (4)
The first sum in the r.h.s. of eq. (4) represents the “distance” (to the square) between the reference trajectory and the predicted output computed at time k according to eq. (3). This distance is integrated over a certain time window which opens at time If,. This parameter can be varied between the two bmits: r < N, < H, + r. The second term in eq. (4) takes into consideration the amplitude variations (also to the square) of the input, up to a time limit H, named control horizon. Necessarily, one has 1 < H, -c H,. Beyond H,, nothing is claimed about the changes that could be imposed on u. Because the function J(k) mixes together two different kinds of “deviation”, it is appropriate to have a means of balancing the two contributions. This is achieved through the control weighting parameter y.
=
0.49 + 0.259-l 1 - 0.49-I
+0.149-z
appeared to be the most appropriate with respect to the various time-scale and delay involved in the experimental procedure. Starting with an input-output vector x(k) and a vector of model parameters 8, x=(k) = [Ay(k -
l), . . . ,Ay(k - n.)
Au@ - I - 1), . . . , Au(k - r - nb)] @ = Cal,.
. aran,, bl,. . . , h.,]
the RLS algorithm writes down in standard notations: B(k) = 8(k -
1) + L(k)&(k)
e(k) = By(k) -
x*(k)B(k
-
1).
The gain vector L(k) = P(k - l)x(k)/[A(k)
+ xT(k)P(k
-
1)x(k)]
Adaptive-predictive temperature control of semi-batch reactors and the covariance matrix P(k) = [P(k -
1) - T@)x=(&)P(k - 1)1/R(k)
both involve a forgetting factor 1. This parameter, which is less than unity, is chosen at will. It allows to balance the data taken into consideration with respect to their “age”. Indeed, it is appropriate to forget progressively the former states of the system at an adjustable rate, characterized by a time constant = equal to i/(1 - A). It was soon realized that a difficulty has to be overcome when applying to the RLS algorithm to filtered data in the form presented above. As a matter of fact, the calculus must be modified in some way if a high enough sensitivity to process variations is to be preserved. Several variants have been proposed over the past years, a critical survey of which has been published by Shah and Cluett (1991). In the present work, we brought into play two approaches aimed at improving the dynamics: (i) According to Fortescue et al. (1981), the parameter I is changed at each step in such a way that a link is established between the amount of forgotten data and the new information supplied by the last measurement. To this end, the new value of 1 at time k + 1 is selected by I(k + 1) = max {l - Cl - xT(k)L(k)]s(k)‘/(rO,
A,,,,,}
where G,, and &,, are two adjustable parameters. Aminrepresents a lower limit that parameter n must not cross (otherwise, the covariance matrix would diverge). On the other hand, the smaller the no the larger is the sensitivity of the algorithm. (ii) As suggested earlier by Irving et al. (1979), one can also give the trace of matrix P a fixed value /I by varying jl properly. In this case, the gain factor is determined without data balancing: L(k) = P(k -
l)x(k)/[l
+ x’(k)P(k
-
1)x(k)].
Then I(k) is calculated so as to satisfy the condition tr P(k) = tr P(k -
1) = p.
Of course, the larger the B the more sensitive is the corresponding estimator. Finally, since it is necessary to have a positivedefinite covariance matrix, the two variants are improved by including a U-D factorization technique (Bierman, 1977) in the computation of P(k).
4. CONTROL
1MPLEMENTING
is the heat release by the chemical reaction which is responsible for the main nonlinearity of the process, due to Arrhenius’s law (exponential dependence of rate constants upon temperature). Viscosity, on the other hand, is rather a “slave” property. Uneasy to manage at will, it has almost no direct influence on the reaction rate, despite its effects on stirring efficiency and, also, on the heat exchange coefficient (Nusselt number). The only relevant control variable one can It
3317
act upon is thus the temperature. The adaptive control algorithms described in the previous section are aimed to cope with this situation and to some other features of semi-batch reactors (e.g. the heat removal capacity changes over the time since the reacting volume doubles in the present case). The goal one must reach is twofold: (i) maintain the medium temperature at a fixed constant value once the setpoint is established and (ii) reduce as much as possible the amplitude of the unavoidable temperature variations resulting from semi-batch operations (especially initiator addition). To this end, one has to vary, properly and in due time, the temperature Ti of the thermofluid. We did it in the following way. The cooling water flow F, is usually set at a fixed mediumrange value (namely 15 l/min), and the thermofluid is warmed up to the desired temperature by the two resistors. The applied electric power P, is thus the control parameter actually managed. However, when the medium temperature T, exceeds the setpoint by too large a difference (say l”C), the resistors are switched off (P, = 0) and the flow F, is raised to its maximum (i.e. 32 l/min). As soon as the temperature T, begins to diminish, the “standard” control conditions (i.e. F, = 15 I/min, P, variable) are restored. In practice, the transfer function model, assumed to represent how the temperature T, depends upon the delivered electric power P,, is the key ingredient of the control. To build up successfully this function during the estimation phase, data must be filtered as indicated above. Moreover, specific care has to be taken in order to eliminate bias in estimated parameters issuing from the presence of the second feed stage. As a matter of fact, the beginning and end of that stage correspond to “extraordinary” events, which will not occur twice. So, one has to cancel the spurious correlation which is established at these times by least-squares computation (Inglis et al., 1991). The time characteristics of our experimental device with respect to T, variations (reactor containing 10 I of solvent, open-loop condition) have the following order of magnitude: time constant _ IO3 s heating/cooling dead time - 180 s. Owing to this, we
used
our algorithms:
-a sampling period of 120 s for parameter estimation (this period may seem rather long; it was chosen essentially according to the incremental model [eq. (S)], which sometimes yields numerical problems when smaller periods are taken); ~ a control performing period also equal to 120 S; -a control data recording period of 20 s for better control watching. 5.
RESULTS
AND DISCUSSION
The main results of this study are summarized in Figs 3 and 4, which represent the variations vs time of three quantities: measured temperature T,, electric
G. DEFAYE
3378
power P, represented as a percentage of 6 kW, cooling water flow F,. Moreover, Fig. 6 displays similar data obtained with commercially available adaptive controllers. In all cases, the applied experimental procedure was that depicted in Section 2.3.. What changes from one experiment to another is uniquely the way in which control is carried out. Consequently, comparison is really meaningful. Figure 3(a) was obtained under “control-free” conditions with respect to T,. In other words, the temperature T, was solely recorded, but not involved in any control algorithm. The heating/cooling setup was used to maintain the temperature 7;: at 76°C independently of what happened in the reactor. Now, during feed 2, the reaction rate was very large and, as a consequence, so was the rate of heat release. A point was T, became higher even reached where the temperature than the boiling temperature of the solvent (i.e. 77°C). Nevertheless, due to solvent evaporation on one hand, and to flow back to the reactor after solvent condensation, the rise in temperature remained slow enough. Even at the end of feed 2, T, did not exceed 79°C. Then it decreased, as did the reaction rate. In Fig. 3(b) a predictive control algorithm was employed. To this end, a transfer function model was determined by open-loop identification with the following experimental conditions: - reactor containing 10 1 of solvent -FF, = 15 l/min -P, = 5 kW + a pseudo-random binary quence superimposed at 600 s time intervals. This led to the first-order T,(k) = O.S9T,(k -
estimation
se-
model
1) + O.O57P,(k
-
2)
then applied once and for all. The controller tuning parameters yielding the best result [i.e. the minimum value of J(k), see eq. (4)] were H, = 5,
H, = 1,
y = 0.06,
a, = 0.74
et al.
corresponding to a time constant of 400 s. In some sense, this experiment provides us with a reference point. As far as the viscosity of the thermofluid is temperature-dependent, as it actually is, the flow rate Ff will vary depending on whether the fluid is warmed or cooled. Accordingly, one must also expect variations in the time delay. This fact has to be taken into account if an adaptive-predictive control procedure is to be applied. To model these unavoidable variations, we found it appropriate to take a third-order polynomial B(q-‘). For the two estimation algorithms described in Section 3.2 (namely Fortescue, and constant trace), we thus referred to the same kind of model given by AT,(k)
= a,AT,(k
1) + blAP,(k
+ bt AP,(k
Fortescue: Constant
&in = 0.8, trace:
-
1)
- 2) + b3 AP,(k
The selected tuning parameters algorithms were, respectively,
of these
- 3). estimation
(r0 = 1
p = 0.04.
Furthermore, the same set of control parameters was used for both algorithms, so as to allow for a fully significant comparison of their performance. As may be seen, the Fortescue estimation [Fig. 4(a)] undoubtedly provided better results than did the constant trace [Fig. 4(b)]. It was obviously more effective, not only in lowering the temperature rise at the beginning of feed 2 but also in damping oscillations around the setpoint. As already mentioned, the input and output data that served to perform identification in RLS algorithms were filtered. The applied second-order filter has a damping factor of 0.7; it washes out perturbations whose time constant is shorter than 540s. Thanks to this, the adaptive control is insensitive to perturbations which should remain unmodelled, in
80
ä
4 78
-
feed2
78
76
20 0 0
50
loo
150
200
250
(a> Fig. 3. Observed
lemperature
300
350
Time (min)
400
66! 0
! 50
103
150
200 (b)
and applied heating/cooling parameters as functions free” condition (fixed Ti); (b) predictive control.
250
300
350
Time (min) of time: (a) “control-
400
0
Adaptive-predictive
temperature
control of semi-batch
reactors
3379
80
78 T@)
76
4
74 -
72
loo
WfJo Pw(%) +* 60
70
- 40
-40
68
II
Fwo/min) *-
66 +
F,#W
20
r0
I-0 50
0
loo
150
200 (a)
250
300
350
400
. . 20
0
50
loo
150
Time (min)
200
250
(b)
300 Tie
350
400
(min)
Fig. 4. Observed temperature and applied heating/cooling parameters as functions of time: (a) Fortescue estimation algorithm; (b) constant trace estimation algorithm.
particular
those
occurring
at the beginning
and at the
end of feed 2. Indeed, the corresponding variations have time constants well below that of thermal events, whereas one has to determine the transfer function linking P, and r,, and nothing else. The effect of data filtering on the model time constant is shown in Fig. 5 in the case of the Fortescue algorithm. As may be seen, a drastic drop in the time constant is observed right after initiator addition and the beginning of feed 2, when data are not filtered. This drop is undoubtedly irrelevant: indeed the time constant is three times smaller than that calculated for a reactor working in open loop. The slight modification resulting from adding initiator is not at all able to account for such a large ratio. To the contrary, no trouble of this sort is met when filtering of data is achieved. This result points out the usefulness and even the need of data filtering at least in the presence of large external perturbations. Trademark adaptive controllers being readily available, it is interesting to look at their control performances when driving a semi-batch reactor. Two such controllers were tested: the Foxboro EXACT and the A.B.B. Novatune. They yielded the results shown in Fig. 6. The Foxboro EXACT (Kraus and Myron, 1984) is basically a discretized version of an analog PID. Owing to this, its sampling period has to be small (0.1 s). To set the initial values of its parameters, we filled the reactor with 10 1 of solvent, as usual, and then we let the controller analyze the transient response of the closed-loop system to setpoint changes. Temperature data filtering was performed by a second-order Butterworth filter, with a 30 s time constant. The Foxboro system is designed so that, during a run, parameter adaptation takes place whenever the difference between the measured temperature (T,) and the setpoint exceeds twice the (allowed) noise amplitude. In the experiment presented in Fig. 6(a), this amplitude was assigned to the lowest adjustable value, namely 0.5”C.
soo.
Fig. 5. Effect of data filtering on the model time constant (Fortescue estimation algorithm).
The nature of the A.B.B. Novatune system (Bengtsson and Egardt, 1984) is much closer to our approach. Indeed, it involves a RLS estimation algorithm (forgetting factor k = 0.98) together with a minimum variance controller. It relies upon the following model: T,(k + H,)
- PLT,(k
+ H, + 1) - (1 - PL)T,(k)
= A(q-‘)AT,(k)
+ B(q-‘)AP,(k)
where PL is the desired pole location in closed loop. After various trials, we selected the following conditions: = 3 (with a 60 s sampling time interval) = 0.86 (so that the model time constant still equal to 400 s) A (q - ‘) = second-order polynomial B(q- ‘) = third-order polynomial.
HC PL
is
It is worth noting that, here too, parameter adaptation is not active at all times, but only when variations
G. DEFAYE er al.
3380 80 78
TR(*c) t
76
76
74 100 72
80 Py&w )
70
60 40
68
20
FwWmin).
0
66 0
50
100
150
200 (a)
250
300
350
0
400
SO
loo
150
UK)
Time (min)
250
300
350
400
Time (min)
@I
Fig. 6. Observed temperature and applied heating/cooling parameters as functions of time: (a) Foxboro EXACT adaptive controller; (b) A.B.B. Novatune adaptive controller.
e control - constant trace
0
50
100
150
200
250
300
350
400
Time @in) Fig. 7. Viscosity vs time plots for the various tested controllers
in either temperature T, or electric power P, exceed an adjustable threshold (respectively, 0.4”C and 2%). Even at a glance at Figs 4 and 6 the improvement brought out by our adaptive controllers is indisputable. Perhaps the most striking feature is the trend towards oscillation exhibited by the process. This trend is clearly emphasized by the last two experiments, and was already visible before. Whatever the criterion taken into consideration, the superiority of the Fortescue algorithm is prominent. It leads to smaller temperature differences with the setpoint and to fewer changes in electric power supply. It also provides a much more efficient damping of oscillations compared to any other algorithm. Although the differences might perhaps seem relatively weak, they are still of consequence, owing to the sensitivity of the end product to temperature vari-
ations having occurred during the synthesis. An overall indication, supporting this conclusion, appears in Fig. 7. There, one can see how the viscosity q has evolved during the course of the reaction. The rapid increase at the end of each curve corresponds simply to the outset of quenching and reflects the final temperature drop. As expected, the most regular increase in viscosity has been obtained with the Fortescue adaptive controller. Curve 3 has no bump at all; its slope always changes in a progressive manner, indicating that the copolymerization proceeds smoothly. It even continues to do so beyond the desired threshold of 4 P. Furthermore, except in the case of the A.B.B. controller where something peculiar has occurred presumably due to the chemicals, the reaction is also more rapid and the targetted viscosity is reached sooner.
:e control of semi-batch Adaptive-predictivetemperatul
6. CONCLUSION this work, we put into shape an adaptive-predictive controller combining the following: In
(i) An estimation algorithm based on RLS techniques; a thorough investigation led us to select two variants of this type, one having a variable forgetting factor (Fortescue) and the other preserving the trace of the covariance matrix (constant trace); (ii) A stage of data filtering aimed at eliminating spurious dynamical effects, in particular those linked to reactants addition(s). In this way, we were able to control successfully the temperature of a reactor where a copolymerization reaction took place despite the perturbations due to initiator addition and to secondary feed. Particularly worth of note is the complete
vanishing
of
temperature
(Fortescue) that commercial failed to achieve.
adaptive
e
FJF. FW HC
HP k
L nn, nb PW P
parameters defined in eq. (1) model error flow rate of thermofluid, l/h flow rate of feed 2, I/h flow rate of cold water, l/min control horizon prediction horizon discretized time gain vector polynomial orders electric power, kW covariance matrix
time delay stirring rate, rpm reference trajectory thermofluid inlet temperature, “C thermofluid outlet temperature, “C reacting medium temperature, “C model input setpoint input-output vector model output
r r.4 b
Ti T, T, u W
x Y
Greek letters
control weighting factor forgetting factor reaction medium viscosity, P vector of model parameters
: I
controllers
Acknowledgements-This research was granted by the CNRS (ARC “Optimal Control and Automation of Batch Reactors”) and by the Regional Council of Aquitaine. The authors would like to acknowledge assistance and fruitful discussions with their colleagues from “Institut du Pin” (Bordeaux University).
ai, bj
3381
oscillations
Seeking for the main source of the observed improvement, we came to the conclusion that data filtering was crucial. Indeed, the cut-off of the dynamical perturbations which result from some sort of “foreign event”, is essential in preserving the right set of estimated parameters. Otherwise, troubles ensue. As is wellknown cut-off is all the more effective when the timescales involved are well separated. Obviously, the process we tested fulfilled this condition. It must be emphasized, however, that the chosen example is in no way peculiar: it is a mere archetype of semi-batch processes. Furthermore, if the time-scale separation be missing, one could still expect this kind of adaptive-predictive controller to be able to do well, thanks to its intrinsic flexibility. As a matter of fact, spurious dynamical perturbations developing on a time-scale comparable to that of the process will no longer induce errors in parameters estimation. Accordingly, a wide range of applications is open, in principle, to this kind of controller, at the expense of an appropriate fitting to the process dealt with.
NOTATION
reactors
REFERENCES Bengtsson, G. and Egardt, B., 1984,Experiences with selftuning control in the process industry. Preprint 9th fFAC World Congress, Budapest, Vol. XI, pp. 132-140. Bierman, G. I., 1977, Factorization Methods for Discrete Sequential Esfimation. Academic Press, New York. Clarke, D. W., Mohtadi, C. and Tuffs, P. S., 1987, Generalized predictive control. Automatica 23, 137-160. Cluett, W. R., Shah, S. L. and Fisher, D. G., 1985, Adaptive control of a batch reactor. Chem. Engng Commun. 38, 67-78. Cott, B. S. and Macchietto, S., 1989, Temperature control of exothermic batch reactors using generic model control. Ind. Engng Cbem. Res. UI, 1177-I 184. Defaye, G., Caraip, L. and Jouve, P., 1983, A simple deterministic predictive control algorithm and its application to an industrial chemical process: a distillation column. Chem. Engng J. 27, 161-166. De Keyser, R. M. C., 1990, Model based predictive control toolbox. Proceedings of CIM Europe Workshop on Computer Integrated Design of ControHed Industrial Systems, Paris, pp. l-56. Fortescue. J. R., Kershenbaum, L. S. and Ydstie, B. E., 1981, Implementation of self-tuning regulators with variable forgetting factor. Automatica 17, 831-835. Garcia, C. E., Prett, D. M. and Morari, M., 1989, Model predictive control: theory and practice-a survey. Automatica 25, 335-348. Inglis, M. P., Cluett, W. R. and Penlidis, A., 1991, Long range predictive control of a polymerization reactor. Can. J. them. Engng 69, 120-129. Irving, E., Barret. J. P., Charcossey, C. and Monville, J. P., 1979, Improving power network stability and unit stress with adaptive generator control. Automatica IS, 31-46. Juba M. R. and Hamer, J. W., 1986, Progress and challenges in batch process control. Proceedings of 3rd Inrernational Conference on Chemical Process Control. Elsevier.New
York.
Kraus, T. W. and Myron, T. J.. 1984, Self-tuning PID controller uses pattern recognition approach. Contr. Engng 6, 106-11 I. Lee, P. L. and Sullivan, G. R., 1988, Generic model control. Comput. them. Engng 12, 573-598. Liptak, B. C., 1986. Controlling and optimizing chemical reactors. Chem. Engng 24, 69-81. Ljung, L., 1987. System Identification: Theory for the User. Prentice-Hall, Englewood ClilTs, NJ. Shah, S. L. and Cluett, W. R., 1991, Recursive least-squares based estimation schemes for self-tuning control. Can J. them. Emma 69.89-96. Tuffs, P.-S.“a&C$larke, D. W., 1985, Self-tuning control of offset: a unified approach. Proc. IEE, Part D 132, pp. _ 100-110.
G. DEFAYE et al.
3382 APPENDIX
We briefly summarize hereafter the way in which we determine the input u(k) that is applied at each step k. The complete derivation of these equations may be found elsewhere (Defaye et al.. 1983). The output values, which are predicted according to eq. (3). gather three contributions which may be written as y(k + m/k)
Ear y,l(k +
m/k)1
+ effect of future inputs
Car YAk +
m/k)1
+ predicted error
[or e(k + m/k)].
= effect of past inputs
Identification
+ i condition
form=I,...,r,
notation.
y=y@+Cu+e with
+ m - ilk)
b&k
+ m - j - r) + e(k)
when summing
over j:
“=a, n=nb-m+r
y,,(k + m/k) = 5 c,u(k + m - i) i= 1
+ H, + r/k)]
II’ = [u(k). . . , u(k + H,, - I)]. C is the lower triangular matrix whose elements Cij are Cii = q+
The input vector u corresponding given by a = (C’C + M)-‘{C’(s,
-j_
to the minimum
of J is
- yri - e) + pu(k - I)}
where
P’=cY.4..~701
c, =b, i-L
. ,nb, Cl= bi + c ca,Ck-,+ b, j-1
fori=n~+l,....H,+r
t ci= x a~c~-~ with /=min(i-
I
sT=[s,(k+r+l),...,s,(k+r+H,)]
the coefficients c, being given by
i=*
. . .y(k + HP + r/k)]
+ r + l/k), . . . ,y,&
e7=[Ce(k+r+l/k),....e(k+H,+r/k)]
(of course, no other contribution of the inputs has to be considered further). The yn are written as follows:
for i = 2,.
H, = H,.
it then becomes
y’ = [y(k + r + I/k),
form=r+I,....r+nb-1,
I,
In vector
yz = b,,(k
j-l
with the twofold
I;
H,=r+
of the y+ yields
I. y,f(k + m/k) = E atyAk 1-1
for i=
adaptiveness comes out from the calculus of the cis at each sampling time k. The purpose is to determine the future input giving the function J(k) [see eq. (4)] its minimum value. Being more specific now, we ascribe to parameters H1 and H, the particular values
1.n.).
In a nonadaptive control algorithm, these coefficients e, are calculated once and for all at the outset. On the other hand,
M=
;?r --y
-Y
0
0
0
2y
--y
0
0
0
--Y
2Y
0
0
. 0
0
0
2Y
-Y
0
0
0
-Y
Y
y being the control weighting parameter which balances the two types of deviations appearing in the criterion J(k) [see =I. (411.