Chemical Process Controller Design Using Genetic Programming Dominic Searson
Mark Willis
Gary Montague
Advanced Control Group Chemical Engineering Dept. University of Newcastle Newcastle upon Tyne NE1 7RU, UK
[email protected]
Advanced Control Group Chemical Engineering Dept. University of Newcastle Newcastle upon Tyne NE1 7RU, UK
[email protected]
Advanced Control Group Chemical Engineering Dept. University of Newcastle Newcastle upon Tyne NE1 7RU, UK
[email protected]
ABSTRACT In this article, we show that genetic programming (GP) can be used to design discretetime dynamic controllers that offer similar performance to standard Proportional Integral Derivative (PID) controllers for a specific class of control objectives. Two processes are used as examples: an AutoRegressive eXogeneous (ARX) system and a simulated nonlinear Continuous Stirred Tank Reactor (CSTR). Additionally, some of the generalisation properties of the evolved controllers are examined and discussed.
1. Introduction Most common chemical process control systems can be classed as either reference control systems or regulatory control systems (Isermann, 1989). In simple reference control systems the objective is to get the process output to follow a reference signal (or setpoint) as closely as possible, preferably without using excessive control effort. In the case of regulatory control systems the objective is usually to hold the process output at a desired value in the presence of disturbances. Often the same controller can be used to accomplish both tasks for a given process. Given the above objectives there are, broadly speaking, two ways to design an appropriate control system. The most frequently used method is to preselect a controller structure and then to tune the parameters of this controller so that the desired closedloop response is obtained. This is referred to as a parameter optimised control system, the most well known example of which is probably the Proportional Integral Derivative (PID) controller. The other approach is the use of structure optimal control systems, where both the structure and parameters of the controller are adapted to those of the process model (Isermann, 1989). In practice, however, the use of the latter method is severely restricted because exact dynamic term cancellation is required in order
to produce the optimal controller structure. This is usually not possible for various reasons, e.g. of the lack of an appropriate process model, nonlinearities and physical constraints on the process variables. Genetic programming offers an interesting new perspective on structure optimised controller design because of its potential ability to match controller algorithms to the process model by implicitly utilising the information about constraints and nonlinearities that the closedloop simulation provides. GP is also an attractive proposition for control system design in that the performance criterion used to perform the optimisation does not have to be continuous or differentiable. This means that, potentially, the designer can introduce more ‘real world’ performance measures to be used directly by the GP algorithm. Recently there have been a number of articles describing the use of GP in designing control algorithms of various types. Howley (1996) used GP to derive near optimal control laws for simulated specific spacecraft attitude maneuvers but found that they did not always generalise well to arbitrary maneuvers. Howley (1997) also used GP to design algorithms for a two link payload manipulator. Dracopoulos (1997) used GP to derive an attitude control law to detumble a spinning satellite, he also used a Lyapunov function to show that the derived law was stable, even though stability was not explicitly considered in the search criteria. Alba et al. (1996) used GP based fuzzy logic controllers to tackle the cartcentering control problem described and solved by Koza (1992). This paper is structured as follows: Section 2 briefly reviews a common reference control problem and gives examples of which performance criteria may be used in attempting to solve it. The operation of the widely used PID controller is also described as it is used for comparison against the controllers designed using GP. Sections 3 and 4 are used to describe the characteristics of the two case study processes: a constrained second order ARX system and a nonlinear CSTR. Section 5 details the formulation of the control problem and performance criteria in relation to the GP. Section 6 details the implementation and configuration of the GP. Section 7 presents some results and Section 8 concludes with a brief discussion and future work.
2. Dynamic Reference Control Problem In this article, GP is used to design recursive discretetime reference control systems. In both case studies the objective of the GP derived algorithm is to ensure that the process output follows the setpoint in a desirable manner through step changes of differing magnitudes across a predefined operating region. Figure 1 shows the type of control configuration used in both process case studies. disturbances Setpoint

Actuator constraints
Controller
Process
+
Process output
can be derived from Equation 2 by the use of backward differencing.
u(k ) = u(k − 1) + q0e(k ) + q1e(k − 1) + q2e(k − 2) Where:
Td q0 = K 1 + Ts q2 = K
(4) q 1 = − K 1 + 2 Td − Ts Ts Ti
Td Ts
(5) (6)
A typical PID closed loop response to a series of step changes in setpoint is shown in Figure 2. The response shown is for illustrative purposes and is not considered to be welltuned. 70
Figure 1 Control configuration used in GP design
65
t
60 55 Output
A controller that is used in many real chemical process control problems of the type shown in Figure 1 is the PID controller, described in its original continuous time form by Equation 1.
1 de ( t ) + u ss u ( t ) = K e (t ) + ∫ e (τ ) dτ + Td T dt i 0
(3)
50 45 40 35 30
(1)
Where u(t) is the controller output at time t, e(t) is the setpoint error at time t and uss is the controller output at some reference steady state. The PID controller is parameterised by the proportional gain K, the integral time Ti and the derivative time Td. These are adjusted (tuned) in order to give the desired closed loop response. The Ke(t) term produces a controller output proportional to the size of the error but use of proportional control alone produces tracking offset when applied to most processes i.e. those without natural integrating properties. To counter this, the integral term in Equation 1 is used to automatically ‘reset’ the offset. Additionally, the derivative term may be used to provide anticipatory control by introducing a control effort proportional to the rate of change of setpoint error. This term becomes dominant when the setpoint changes suddenly, e.g. in a step change, and so provides a faster initial response. However, many designers choose to set the derivative term to zero as it has an amplifying effect on noisy measurements. With the increasing use of digital computers for chemical process control the discrete form of the algorithm is now prevalent. For small sample times, Ts, the discretised version of the PID algorithm may be obtained by applying rectangular integration to Equation 1 to yield Equation 2. Ts k −1 Td (2) u ( k ) = K e ( k )+ ∑ e (i ) + (e ( k ) − e( k − 1) ) + usss i s T T i =0 The control signal u(k) is now calculated every sample time and, in the case of a continuous process, a zeroorder hold element is usually used to hold the control signal steady between sampling times. The recursive form of the algorithm, Equation 3, is often preferred because it is computationally simpler and because no special provision need be made for smooth switching from manual mode (open loop) to automatic mode (closed loop). Equation 3
0
50
100
150
200 T im e s t e p ( k )
250
300
350
400
Figure 2 Example of PID reference control. It can be seen that, although the response rises quickly to the setpoint, there is some degree of overshoot and oscillation. This example is to illustrate that when tuning a PID controller one must tradeoff the various transient response properties, e.g. settling time, rise time, overshoot, and degree of oscillation. The desired properties of the response depend largely on the preferences of the designer and the problem context so one cannot easily specify a generic best response that holds over all processes. For example, the designer of a reboiler control system for an industrial distillation column would also want to consider the cost of the control effort (i.e. steam usage) in deciding the best parameter values. In cases such as these, aggregate cost functions containing terms pertaining to actuator effort, and a performance measure such as IAE (Integral of the Absolute Error) or ITAE (Integral of the Time weighted Absolute error) can be utilised, and the parameters can then be optimised using numerical procedures. However, any performance criterion that reflects the designer’s a priori qualitative description of a ‘good’ closedloop response may be used.
3. ARX Process Description The ARX process used is of the linear form shown in Equation 7. (7) Ay (k ) = Bu (k − 1) + Cε (k ) Where y(k) and u(k1) are the process output at time k and the controlled input at time k1 respectively. A, B and C −1 are polynomials in the backshift operator z and ε (k ) is a stationary white noise sequence. The polynomial −1 coefficients used (in order of ascending powers of z ) were A= {1 –0.99 0.05}, B = {0.1 0.01} For the purposes of
this study C was set to zero to yield an entirely deterministic system. In practice, the response of the output sequence y(k) to u(k) can most easily be calculated by rearranging Equation 7 and, after substitution of the appropriate coefficients, proceeding in the recursive manner described by Equation 8. y(k ) = 0.99y(k − 1) − 0.05y(k − 2) + 0.1u(k − 1) + 0.01u(k − 2) (8)
Additionally, when this simulating this process in the configuration shown in Figure 1 the following constraints were imposed on the system: •
Setpoint changes: constrained between limits of 60 units (high) and 20 units (low). Process output also limited between these constraints.
•
Process Input u(k): constrained between hard limits of 90 units (high) and 0 units (low) and a ratelimit of a maximum of magnitude 4 units change per discrete time step.
− EA
Vr
dC A = − K A e RTr C AVr + QF (C AF − C A ) dt
Vr
dCB = KAe RTr CAVr − KBe RTr CBVr + QF (CBF − CB ) dt
−EA
Feed
− EA
Heat exchanger Product
Figure 3 Schematic of CSTR and heat exchange system. The following state variables are defined: Tr is the temperature of the reactor product stream, TRC is the temperature of the reactor product recycle stream and TCO is the temperature of the coolant stream on leaving the heat exchanger. CA and CB are the concentrations of Reactant A and Reactant B in the reactor product stream. Equations 9 and 10 are material balances on reactants A and B respectively. Equation 11 describes the thermal energy balance over the reactor. Equation 12 describes the energy balance for the reactor product recycle stream and Equation 13 describes the energy balance for the coolant passing through the heat exchange unit.
(11)
Additionally, the following constraints were imposed on the system: •
Reactor temperature, Tr : constrained to follow setpoints generated between limits of 373 K and 325 K.
•
Cooling water flowrate, QCO: constrained between hard limits of 0 litres/sec and 10 litres/sec. A ratelimit of a maximum cooling water magnitude change of 5 litres/sec per 5second interval (i.e. every discrete time step) was also imposed.
The physical parameter values for this system are shown in Table 1 and were fixed for all the simulations carried out. Table 1 Model Parameters for CSTR system Symbol
Description
Unit 3
Reactor volume Rate coefficient AÆB
m
Rate coefficient BÆC
mol/sec
Activation energy AÆB Activation energy BÆC
J
R
Gas constant
QF
Reactor feed flowrate Heat exchange recycle flowrate Conc. of reactant A in feed stream Conc. of reactant B in feed stream Process fluid density Process fluid heat capacity CSTR heat transfer coefficient Heat exchange area Heat of reaction AÆB Heat of reaction BÆC Heat exchanger process
J/mol K m3/ sec
Vr KA
EB
Cooling Water Q CO
− EB
dT RC = ρ c pQ RC (Tr − T RC ) (12) dt + U H A H (TCO − T RC ) dT ρ c pVCL CO = ρ c pQCO (TCIN − TCO ) + U H AH (TRC − TCO) (13) dt
EA
Tr
(10)
ρ c pV RC
KB
Reactor
−EB
dTr =∆H A K A e RTr CAVr + ∆H B K B e RTr dt + ρ c p Q F (TF − Tr ) + ρ c p Q RC (TRC − Tr )
ρ c p Vr
4. CSTR Process Description The simulated continuous stirred tank reactor system used in the second case study is shown in Figure 3. A feed stream that has a high concentration of reactant A enters the reactor. Then, within the reactor, the following exothermic reaction takes place: A → B → C. Reactant C is regarded as the undesired byproduct of the reaction and at high temperatures most of reactant B is converted to C. The reactor is cooled by means of an external pumparound heat exchanger. The reactor product stream temperature, Tr, and hence the concentration of the various products within it, is controlled by manipulating the flowrate of cooling water, QCO, to the heat exchanger.
(9)
QRC CAF CBF
ρ Cp UH AH ∆HA ∆HB VRC
mol/sec
Value 3 7x101 1
9x101 1
J
90000 10000 0 8.314
m3/ sec
0.015 0.025
mol/m3
2500
mol/m3
0
Kg/m3
1000 4140 3000
J/kg K W/m2K m2 J/mol J/mol m3
100 80000 40000 0.2
fluid volume 293 Cooling water entry K temperature TF K 320 Reactor feed temperature VCL 0.2 Heat exchanger cooling m3 fluid volume The CSTR system equations (McKay et al., 1997) were simulated using a RungeKutta 4th order method with a fixed step size of 5 seconds, both the integration method and step length were determined in a trial and error manner. TCIN
5. GP Problem Formulation The GP algorithm, for both case studies, is set up to construct a discrete recursive feedback control law of the following form: u(k )= u(k1)+ [output of individual tree] For a population size of M the output of the ith {i=1,2,3…M} controller at time k is equal to the output at timestep k1 plus some correction term applied by the ith GP individual. Hence, the controller can be regarded to have a state that is updated by the GP individual (based on current information about the process supplied by an appropriately designed terminal set). The controller output is only limited by the physical constraints peculiar to each process and need not take an on/off only value. The recursive formulation of the controller was chosen because it simplifies computer simulation of the closedloop systems and has several other properties that make it a more practical implementation than the nonrecursive form.
5.1 Fitness Function The fitness of an individual, Fi, is constructed as follows:
Fi = ∑ j =1 P
∑
n k =0
k e ( k ) + r∆ u 2 ( k ) ∆S j
(14)
This is calculated by simulating the closed loop behaviour of the ith controller over P independent and randomly generated setpoint changes: ∆ S j {j = 1,2,3…P } Each closedloop step change response was evaluated from 0 to n discrete timesteps. The value of n was determined before the GP run by a trial and error procedure to give sufficient time for a reasonable controller to settle to its steady state value. The performance criterion, Fi , contains an ITAE term that is calculated with respect to some ‘desirable’ reference trajectory (in both cases this was set to equal the setpoint trajectory delayed by one discrete time step). The ITAE weighting factor for the sum of absolute errors was set to be simply the value of the current timestep k. The performance criterion also contains a weighted (with a constant r) penalty term for excessive control effort. In practice this was determined, again, in a trial and error manner before the GP run.
The inner summation term is then scaled according to the size of the step change ∆ S j in order to ensure that contributions from ‘small’ setpoint changes are not dominated by those from ‘large’ changes.
5.2 Randomisation of Training Cases In each case study the setpoint trajectories were generated randomly from generation to generation, thus giving rise to a nonstationary fitness function. The rationale behind this was the attempt to produce controllers that generalise well over the operating region of the process. Some preliminary runs using the GP with the same number of fixed setpoint trajectories throughout the entire run almost always led to very brittle solutions, i.e. controllers that performed badly on step changes they were not explicitly trained upon. This effect appeared to be more deleterious in the case of the CSTR than the ARX system. A possible explanation for this is that the CSTR’s inherent nonlinearities means that the limited subset of training cases utilised in the GP run are less likely to accurately encapsulate the dynamic behaviour across the entire operating region of the process. This leads to poor generalisation properties when the GP is trained on small fixed training set of step changes.
6. GP Configuration and Implementation Aspects Both the GP and the process simulations were implemented within the Matlab programming environment and runs were performed on a desktop PC. The GP algorithm is of a standard nature, using the operators of direct reproduction, mutation and crossover as defined by Koza (1992). Table 2 summarises the configuration and algorithm control parameters used in both case studies. They were decided upon on the basis of available computing time and the results of some preliminary runs. Selection using linear ranking was chosen to avoid problems with fitness scaling, which would be difficult to implement due to the nonstationary fitness measure. The probabilities of the reproductive operators were based on preliminary runs as well as reflecting the fact that contributions in the literature, e.g. Luke and Spector (1997), suggest that the proportions used work well for a number of problems. The functional set is fairly standard except for the ‘backshift operator’ designated as the Z operator in this article and the variable nonlinearity sigmoid, which is designated NLTF (NonLinear Transfer Function). The backshift operator has arity 1 and outputs the value of the input subtree evaluated at the previous timestep. This operator, when nested, has the potential to select terminals and subtrees from several time steps into the past and so any controllers created need not rely on just the past 2 values of the process output and the setpoint error as stated in Table 2. The variable sigmoid NLTF (Sharman et al. (1995) has arity 2 and implements the function described by Equation 15:
1 − e ( − xy ) 1 + e ( − xy )
(15)
Where the values of x and y at any given time step determine the nonlinearity of the sigmoid, which is similar in properties to the basis functions used in many Artificial Neural Networks. The terminal set consists of those terms most likely to be correlated with controller performance i.e. past and present values of errors and process output. Past values of controller output were not included because a number of preliminary runs using them failed to yield any controllers that could produce an acceptable response over a validation sequence of step changes. Table 2 GP run configuration Population size M 200 Termination criterion 100 generations elapsed Fitness cases: Fitness measure Number of discrete steps n Selection method P (Crossover) P (Direct reproduction) P (Mutation) Elite individuals copied to next generation. Max program size Terminal set
Functional set GP runs per case study Wrapper
48 stepchanges generated randomly each generation. The performance measure Fi as defined in Section 5.1 ARX:100 CSTR: 140 Linear ranking. 0.8 0.1 0.1 10 % Approx. 100 nodes. Current and 2 past values of setpoint error e(k). Current and 2 past values of process output y(k). Random constants in range: [100Æ100]. +  % * EXP PLOG PSQRT IFLTE Z NLTF (see below) 25 Case specific actuator constraints.
7. Results After the GP runs were performed each best of run individual (BRI) was tested against a predefined fixed validation sequence of setpoint step changes. This enabled controllers from different runs to be compared in a qualitative manner (i.e. the shape of the responses.) In order to provide a general yardstick for the behaviour of the controllers the validation sequence responses were also compared with those generated by recursiveform PID controllers. These were manually tuned to give a fast response with a low overshoot, and a Simplex optimisation was performed to finetune them further. The manually
obtained PID parameters were used as starting values and Fi (over the validation sequence) as the objective function. (The process examples considered in this paper can be controlled in a near optimal manner, given the actuator constraints, by PID controllers, so we would not really expect any GP derived controllers to do substantially better on these simple test systems.) In addition to the validation step sequence tests it was decided to observe the behaviour of the evolved controllers over a range of conditions that were not present during the training phase, i.e. to see how well they generalised. Two different tests were carried out on each BRI: (1) Setpoint tracking with respect to a differently shaped reference signal (in this case a ramp signal.) and measurement noise on the output. (2) The rejection of a series of load disturbances applied at the process output.
7.1 ARX System Results Of the 25 runs performed, several produced BRIs that could be simplified to linear recursive form PI or PID controllers, but most produced controllers that were not of the PID type and often included the nonlinear functionals such as NLTF. Additionally, when tested over the validation step change sequence, nearly half the runs produced BRIs that were deemed unacceptable due to anomalous behavior; e.g. they had large offsets or very oscillatory responses. Figure 4 shows a typical closedloop response of the ARX process using a nonlinear BRI and the nominal model. The response of the PID controller is also shown for comparison. Qualitatively, it can be seen from the responses that the evolved controller offers similar performance to that of the PID controller. In the case of a ramped reference signal and noise most of the evolved nonlinear controllers generalised well and gave PID levels of performance (see Figure 5.) 65 E v o lv e d C o n tro lle r P ID C o n tro lle r S e tp o in t
60 55 50
Output Value
NLTF ( x, y ) =
45 40 35 30 25 20 15
0
200
400
600
800
T im e s te p
Figure 4 ARX: nominal process response.
1000
1200
4 4 .8
375
3 0 4 5 E 0 P S 1 2 4 3 4 v .0 ID e t8 C p 2 4 6 o vo l n e i td ro e lC ro n tro e l r
370
4 4 .6
365
Reactor Temp. K
Output Value
4 4 .4 4 4 .2 44 4 3 .8 E v o lv e d C o n t r o lle r P I D C o n t r o lle r S e t p o in t
4 3 .6
355
E v o lv e d C o n tro lle r P I D C o n tro lle r S e tp o in t
350 345 340 335
4 3 .4 4 3 .2
360
330
0
100
200
300
400
500
T im e s te p
0
100
200
300
400
T im e s te p
500
600
700
800
700
800
Figure 7 CSTR: nominal process response
Figure 5 ARX: ramp reference input with noise.
10
90
9 E v o lv e d C o n tr o lle r P ID C o n tro lle r S e tp o in t
80
8 7
Cooling Water Flow (litres/sec)
Output value
70
60
50
40
30
20
0
50
100
150
200 T im e ste p
250
300
350
400
Figure 6 ARX: Load disturbances on output. In the disturbance rejection tests, however, it became apparent that the evolved controllers’ robustness did not extend to classes of control objectives not explicitly described during the training phase. All of the evolved ARX controllers (except the linear ones) became unstable or performed badly at rejecting output load disturbances (see Figure 6.)
7.2 CSTR System Results The proportion of runs that generated acceptable controllers was similar to that obtained for the ARX system and again a proportion of the evolved controllers could be reduced to the form of linear controllers. Figures 7 and 8 show the response of the nominal CSTR system with a nonlinear evolved controller (with a PID controller for comparison) and the corresponding cooling water flowrates. It can be seen that although the controllers have very similar performance, the evolved controller uses less control effort. On generalisation testing, the BRIs exhibited no real difference to the ARX case, and behaved in a very similar manner on the ramp tracking problems. The nonlinear controllers evolved for the CSTR also tended to go unstable, or offer a very poor response, when presented with the regulatory problem of load disturbances on the output, although due to space limitations, this is not illustrated.
6 5 4 3 2 1 0
E v o lv e d C o n tro lle r P ID C o n tro lle r 0
100
200
300
400
500
600
T im e s te p
Figure 8 CSTR: manipulated flowrate
8. Conclusions It has been shown that the GP framework is capable of producing dynamic recursive controllers for two example simulated processes for a particular class of control problem. It was found that the control algorithms seemed to be able to generalise across the process operating regions without any difficulty when a random fitness case based scheme was used during training. However, as expected, these control algorithms are brittle in the face of some situations they have not encountered e.g. disturbance rejection problems. A major problem is that of controller stability, even for the nominal case there does not appear to be any simple way of incorporating a stability criterion into the GP procedure, or even establishing stability of the controllers a posteori. These issues will have to be addressed if GP is ever going to be used for real controller applications. The fact that GP can offer similar performance to the PID controllers in these cases is encouraging, however. Future work will concentrate on the robustness issues presented by GP derived control laws (e.g. can controller robustness be improved by a more intelligent scheme for presenting the training cases?) Additionally, the technique will be applied to more challenging control problems, i.e. those in which traditionally designed controllers fail to operate effectively.
Bibliography Alba, E., Cotta, C. and Troya, J.M. (1996). ‘Typeconstrained Genetic Programming for RuleBase Definition in Fuzzy Logic Controllers’, GP ‘96Proceedings of the 1st Annual Conference, July 2831, pp. 255260. Dracopoulos, D.C. (1997). ‘Evolutionary Control of a Satellite’, GP ’97 Proceedings of the 2nd Annual Conference, July 1316, pp. 7781. Howley, B. (1996). ‘Genetic Programming of NearMinimumTime Spacecraft Attitude Manoeuvres’, GP ‘96 Proceedings of the 1st Annual Conference, July 2831, pp. 98106. Howley, B. (1997). ‘Genetic Programming and Parametric Sensitivity: a Case Study in Dynamic Control of a TwoLink Manipulator’, GP ’97 Proceedings of the 2nd Annual Conference, July 1316, pp. 180185. Isermann, R. (1989).‘Digital Control Systems’, Vol. 1, 2nd Edition, SpringerVerlag. Koza, J. (1992). ‘Genetic programming: On the programming of computers by means of natural selection’, The MIT press. Luke, S. and Spector, L. (1997). ’A Comparison of Crossover and Mutation in Genetic Programming.’ GP ‘97 Proceedings of the 2nd Annual Conference, July 1316, pp. 240248. McKay, B., Willis, M. and Barton, G. (1997). ‘Steady state modelling of chemical process systems using genetic programming.’ Computers & Chemical Eng., Vol. 21, No.9, pp 981986. Sharman, K., EsparciaAlcazar, A. I. and Li, Y. (1995).‘Evolving Signal Processing Algorithms using Genetic Programming’, IEE conference on GAs in Engineering Systems, Conf. Number 414, September1214.