Constrained nonlinear multivariable control of a fluid catalytic cracking process

Journal of Process Control 10 (2000) 539±555

www.elsevier.com/locate/jprocont

Constrained nonlinear multivariable control of a ¯uid catalytic cracking process R.M. Ansari*, M.O. Tade School of Chemical Engineering, Curtin University of Technology, GPO Box U1987, Perth 6845, Western Australia

Abstract In this paper, a nonlinear constrained optimization strategy is proposed and applied to the reactor-regenerator section of a ¯uid catalytic cracking (FCC) unit. A nonlinear dynamic model of the ¯uid catalytic cracking process was used for the dynamic analysis of the plant and nonlinear multivariable control system. The model realistically simulates the riser-reactor and the one stage regenerator by assembling the mass and energy balances on the system of reactions. The model results were tested in a real-time application and the results were used to provide the initial values for the nonlinear control system design. A dynamic parameter update algorithm was used to reduce the e€ect of large modelling errors by regularly updating the model parameters. The constrained nonlinear optimization algorithm and strategies were tested in real-time on the ¯uid catalytic cracking reactor-regenerator. The results compared favourably to those from a linear multivariable controller. # 2000 Elsevier Science Ltd. All rights reserved. Keywords: Nonlinear process control; Dynamic model; Parameter update; Optimization

1. Background Several dynamic models have been presented in the literature, for the reactor-regenerator system (Kurihara [1], Zheng [2] and Moro and Odloak [3]). Although the engineering principles of all FCC models are the same, the dynamic e€ects vary depending on the geometric con®guration of the system. The scope of the dynamic modelling in this paper is not to describe the complex kinetics of the cracking reactions or the intricate hydrodynamics of the ¯uidized regenerator beds. The scope is to include in the model enough details to capture the control relevant dynamics, without sacri®cing important aspects such as the description of the nonlinearities and the interactions. The FCC unit is a typical example of a constrained multivariable process, where various predictive controllers have been applied, commercially, by Caldwell and Dearwater [4] with reported good results. In the majority of the practical cases, the most pro®table operating point lies on the interception of several FCC constraints. The predictive controllers are capable of including constraints in the formulation of the control law. In the dynamic matrix control * Corresponding author at present address TSU 215, Riyadh Re®nery, PO Box 50, Riyadh 11383, Saudi Arabia. Tel.: +966-014984000/1719; fax: +96-01-4984000/2211. E-mail address: [email protected] (R.M. Ansari).

(DMC) approach presented by Cutler and Ramaker [5] and the model algorithm control by Rouhani and Mehra [6], constraints can be added to the predicted error equations. The controller calculates the manipulated variables that minimize the errors in a least square sense. Only quality constraints can be handled with this strategy. In the quadratic DMC approach, the control problem is formulated as a quadratic-programming problem. The cost function is the square of the distance from the predicted to the reference trajectories. The constraints on the controlled and manipulated variables can be explicitly included and be described rigorously as well. Therefore, both linear programming by Chang and Seborg [7] and quadratic programming by Little and Edgar [8] have been applied in model predictive control for processes with linear constraints. Internal Model Control (IMC) by Garcia and Morari [9±11] as a control framework, is also a linear model-based control. Linear programming developed by Brosilow et al. [12] and quadratic programming by Ricker [13] were also applied for constraint treatment in a similar manner to that applied for the model predictive control algorithms. Economou et al. [14] extended IMC to nonlinear lumped parameter systems by an operator approach. The key issue was the inversion of the nonlinear process model in the control law, and a Newton-type method was adopted. Li and Biegler [15] have extended this operator approach to deal with linear input and state

0959-1524/00/$ - see front matter # 2000 Elsevier Science Ltd. All rights reserved. PII: S0959-1524(99)00059-1

540

R.M. Ansari, M.O. Tade / Journal of Process Control 10 (2000) 539±555

variable constraints by a successive quadratic programming (SQP) strategy. In a di€erent approach, generic model control (GMC), a control framework for both linear and nonlinear systems, was developed in the time domain by Lee and Sullivan [16]. The control law employs a nonlinear process model directly within the controller Lee [17]. An integral feedback term was also included, so that the closed-loop response exhibits zero o€-set. For each of the output variables, two controller performance parameters specify the shape of the closed-loop system response. These parameters are selected a priori by considering the open-loop characteristics such as the process time constant and process deadtime, together with the sampling time interval. GMC has been extended to compensate for process/model mismatch by Lee et al. [18] and deadtime compensation by Lee et al. [19], but the strategies for constraint handling within the GMC framework have not yet been explored extensively, and have not been applied to complex processes in petroleum re®neries. Brown et al. [20] have applied these control strategies to two nonlinear simulated systems. The ®rst system was a forced circulation singlestage evaporator and the second system was a stirred tank reactor from the paper of Li and Biegler [15]. The authors did not consider process-model mismatch in the control strategy. Their work did not address the nonlinear controller implementation issues such as real-time application and on-line tuning of the controllers. The authors did not also address the complexities of handling the multiple constraints and process interaction, which are the characteristics of the FCC process addressed in this paper. 2. Introduction The (FCC) unit converts gas oils into a range of hydrocarbon products of which gasoline is the most valuable. The amount of low market-value feedstocks available for catalytic cracking is considerable in any re®nery. The ability of a typical FCC unit to produce gasoline from low market value feedstocks gives the FCC unit a major role in the overall economic performance of the re®nery. Therefore it is a prime candidate for any type of advanced control applications. The challenge stems from the characteristics and from the fact that the controlled variables must be kept not so much at target but against constraints, where the process often exhibit high nonlinearity. The economic objectives on the manipulated variables, such as the need to maximize either charge rate or cracking severity, or both, complete with the operating constraints further complicate the control problem (Prett and Gillette [21]). This paper provides details of the application of the constraint nonlinear multivariable control and solves a

nonlinear optimization problem on the reactor-regenerator section of a ¯uid catalytic cracking unit. The main contribution of this work is to directly incorporate the nonlinear FCC process dynamic model within the nonlinear multivariable control algorithm and to develop a nonlinear constrained optimization strategy and its application to a highly nonlinear FCC process with multiple constraints and process interaction. In addition, dynamic process models were identi®ed from open-loop step test data and a linear model predictive control algorithm was developed in MATLAB using the features from established model-based control algorithm such as DMC and model algorithm control (MAC). In simulation studies, the linear controller results were compared to the nonlinear controller and it has been shown that the nonlinear controller outperforms the linear controller with regards to its ability to control the process with inverse response. In real-time application, the constrained nonlinear optimization problem was constructed using the simpli®ed FCC process model with dynamic parameter update system. The nonlinear controller was tested and compared to the DMC type controller built with the dynamic models obtained from process identi®cation test. The main advantage of the proposed nonlinear approach is that a single-time step control law resulted in a much smallerdimensional nonlinear programme compared to the other previous methods. The other advantage of this approach is the minimization of the time spent to do the plant testing since it uses nonlinear process models. Therefore it needs much less maintenance and re-identi®cation than the traditional linear controller that often have to be re-identi®ed as the process moves into a new region of operation. 3. The FCC process control overview 3.1. Process description An overview of the FCC process is given in Fig. 1. The FCC process converts heavy oils into lighter and more valuable products. The unit consists of a reactorregenerator section, main fractionator and gas processing facilities. The fresh feed is heated in a ®red heater and combined with regenerated catalyst in the reactor riser. The reaction temperature is controlled at 530 C by manipulating the ¯ow of the hot regenerated catalyst to the reactor-riser. The catalyst and the reacted hydrocarbon vapours ¯ow up the riser and are separated in the reactor cyclone. The catalyst circulation is achieved by burning o€ the coke deposit in a ¯uidized bed inside the regenerator. A steam turbine driven air blower supplies the oxygen needed to burn the coke deposit. The spent catalyst is held up in a small ¯uidized bed in the stripping section of the reactor before being returned to

R.M. Ansari, M.O. Tade / Journal of Process Control 10 (2000) 539±555

541

Fig. 1. Overview of the ¯uid catalytic cracking process.

the regenerator. The air ¯ow for combustion was controlled by adjusting the blade angle of the blower. The ¯ue gas from the regenerator is sent to a waste heat boiler before being discharged to the atmosphere. The FCC unit is operated in full burn mode and the ¯ue gas, therefore, contains a small amount of excess oxygen, typically 1%. Vapour products from the reactor are sent to the bottom of the main fractionator where various boiling point fractions are withdrawn such as distillate, light cycle oil (LCO) and heavy cycle oil (HCO), etc. The feed rate, the reaction temperature, and the regenerator air rate are the main operating variables used for control. There are typically two to four feeds that include fresh feeds from the crude and vacuum units, and the recycle feed from the main fractionator. All of the feeds may be used for control, however, the process response is the same for two or more feeds. 3.2. Economic objectives and nonlinear control strategies The economic objective of the FCC unit is often to maximize the feed rate at constant riser outlet temperature. A more sophisticated approach has been applied here by using the dynamic model of the process to calculate economic optimal values of the key operating variables. Optimal values are calculated for the feed rates, the riser temperature, and the feed preheat temperature and are enforced subject to the operating constraints on the process. In addition, dynamic model parameter update algorithm is also used to update some key operating variables.

The nonlinear multivariable control strategies have been applied to the FCC because of its ability to handle the multivariable interactions and constraints, its ease of con®guration and implementation, adaptability and its robustness. Although it can be applied to the complete FCC process, this paper will focus only on the most complex multivariable controls, those for the reactorregenerator section. In general, the FCC reactor-regenerator nonlinear control strategy is divided into two modules: the regenerator loading control module and the FCC severity control module. The regenerator loading control module maintains the oxygen in the ¯ue gas by manipulating the ¯ow rate of air to the regenerator. The primary control variable is the ¯ue gas composition subject to temperature constraints. In the complete combustion operation, the loading control regulates the ¯ue gas oxygen concentration. In the partial combustion operation, the ¯ue gas CO/CO2 ratio is controlled. On-line analysis is used for the ¯ue gas composition control. The control compensates for changes in the feed and the recycles ¯ow rates, the riser outlet temperature and the feed temperature. Fig. 2 shows the variables of the nonlinear control system for the FCC unit with various system constraints. The FCC severity control module operates the reactor-regenerator at the most desirable point while observing the various process constraints. Optimum operation of the FCC units usually occurs at multiple constraints (Martin et al. [22]). The module manipulates the riser outlet temperature, the feed ¯ow rate and the

542

R.M. Ansari, M.O. Tade / Journal of Process Control 10 (2000) 539±555

feed temperature to control the unit within its constraints while trying to satisfy a speci®ed operating objective. Fig. 3 shows the nonlinear multivariable control strategy on the FCC reactor-regenerator system. In this control strategy, the severity and loading control have been combined and solved simultaneously. The combined strategy also works well if the regenerator becomes air limited.

4. The dynamic model of the FCC process 4.1. Model development The FCC unit consists of a cracking reactor where the desired reactions include cracking of the high boiling gas oil fractions into the lighter hydrocarbons. The undesired reactions include carbon formation reactions,

Fig. 2. Schematic of the nonlinear multivariable control problem for the catalytic cracking reactor-regenerator system.

Fig. 3. Nonlinear multivariable control strategies on the FCC reactor-regenerator system.

R.M. Ansari, M.O. Tade / Journal of Process Control 10 (2000) 539±555

and a regenerator, where the carbon removal reactions take place. Detailed discussions on the features of the FCC unit are given in Huq et al. [23] and McFarlane et al. [24]. The nonlinear process model considered here takes the following simpli®ed form of the model given by Denn [25]. The simpli®ed process model in this work includes enough details to capture the control relevant dynamics, without sacri®cing the nonlinearities and interactions of the FCC process. 4.1.1. The riser and the reactor section The coke generated by the cracking reaction is usually known as catalytic coke. The material balance of the catalytic coke in the reaction section (riser-reactor) can be written as: Vra

dCcat ˆ ÿ60 Frc Ccat ‡ 50 Rcf dt

…1†

Vra

dCsc ˆ ÿ60 Frc …Crc ÿ Csc † ‡ 50 Rcf dt

…2†

where Cca ; Csc ; Crc denote the concentrations of catalytic carbon on spent catalyst, the total carbon on spent catalyst, and carbon on the regenerated catalyst. Vra denotes the holdup of the reactor, Frc denotes the ¯ow rate of catalyst from the reactor to the regenerator and Rcf denotes the reaction rate. In modelling the reaction for the riser, it was assumed that there is a parallel gas-solid two-phase plug ¯ow in the rise-reactor. The reaction rate is kinetically limited. The dynamics of the cracking reaction in the riser is negligible when compared to the dominant time constants of the system. Since the residence time of the hydrocarbons and the catalyst in the riser is very small (about 2 s), the material balance is expressed by the static distribution parameter system. The heat balance is handled in accordance with a lumped parameter system. This leads to the following energy balance equation in the riser: Vra

ÿ
ÿ
dTra Sf ˆ60 Frc Trg ÿ Tra ‡ 0:875 Dtf Rtf Tfp ÿ Tra dt Sc ÿ
Hfv Dtf Rtf ‡ 0:875 Sc ÿ
Hcr ‡ 0:50 Roc Sc …3†

where Tra and Trg denote the temperatures in the reactor and the regenerator, Dtf , is the density of the total feed, Sf and Sc denote speci®c heats, Roc denote the reaction rate, Tfp denote the inlet temperature of the feed in the reactor and Rtf denote the total feed rate.
Hfv and
Hcr are the heat of feed vaporization and the heat of reaction, respectively.

543

4.1.2. The regenerator section The equations of the model for the regenerator are: Coke balance in the regenerator Vrg

dCrc ˆ 60 Frc …Csc ÿ Crc † ÿ 50 Rcb dt

…4†

Energy balance in the regenerator Vrg

ÿ
ÿ
dTrg Sa ˆ60 Frc Tra ÿ Trg ‡ 0:5 Rai Tai ÿ Trg dt Sc   ÿ
Hrg ‡ 0:50 Rcb Sc …5†

The analytical expressions for the reaction rates Rcf , Roc and Rcb were obtained from Denn [28].   kcc Vra Pra ÿEcc exp …6† Rcf ˆ Ccat C 0:06 R…Tra ‡ 460:0† rc Roc ˆ

kcr Vra Pra Rtf Dtf Kcr Rtf ‡ Vra Pra Kcr

…7†

Kcr ˆ

  kcr ÿEcr exp Ccat C 0:15 R…Tra ‡ 460:0† rc

…8†

Rcb

ÿ
Rai 21 ÿ Ofg ˆ 200

Ofg ˆ 21 exp

8 > > > > <

…9† ÿVrg Prg Rai (

> > 106 100 ÿEor > >
Crc ‡ exp ÿ : 4:76R2ai kor R Trg ‡ 460:0

)

9 > > > > = > > > > ;

…10† where kcr , kcr and kcr are the pre-exponential kinetic rate constants, Ecc , Ecr , Eor are the activation energies, Ofg is the oxygen in the ¯ue gas, Rtf is the total feed ¯ow and Rai is the air rate. The values of the process parameters and the corresponding initial steady-state values were obtained from Denn [25]. In the model considered here the coke formation rate given by Eq. (6) and the cracking reaction rate given by Eq. (7) are directly a€ected by the gas oil composition. The two important parameters, namely, Crc and Ccat , involved in Eqs. (6) and (8) were updated using the dynamic parameter update system. The system of ®rst order nonlinear ordinary di€erential equations were solved by the method of Runge-Kutta. The RungeKutta numerical integration technique improves the accuracy of a simple integration by means of interpolation that compensates for the nonlinearity of the equation. From a control viewpoint, a process model should represent the dynamic behaviour of the system inside its

544

R.M. Ansari, M.O. Tade / Journal of Process Control 10 (2000) 539±555

operating range. This implies that the dominant time constants and the gains of both plant and model must match closely. Test runs were carried out on the FCC unit and the results predicted by the model were compared with the plant data. The overall results of the FCC process obtained by the model and the plant test runs are given in Table 1. A number of model validations were needed through many sets of real process data to get these results. The table shows that the maximum relative errors in the riser outlet and regenerator bed temperatures between the model prediction and the plant were 0.3% and 0.5% respectively, whereas the average relative error in the ¯ue gas oxygen is 6.8%. There are large excursions in the plant temperatures, with a slightly larger time constant for the regenerator but with very close gains in all cases. In the comparison test a step change of about 5% was introduced to the setpoint of the air¯ow controller, and the air¯ow pattern is shown in Fig. 4. The step change in the simulation model is of the same magnitude and value as the real disturbance. Figs. 5 and 6 show the plant and model responses for the regenerator and riser temperatures with the changes in the air¯ow rate. The model responses are quite satisfactory. This shows that, for control purposes this model can represent the FCC dynamics with an adequate degree of accuracy. In realtime implementation, for some of the controlled variables, the model step responses were used in the nonlinear multivariable controller design, without extensive process identi®cation test. However, process identi®cation test was performed for the development of linear controller design such as DMC. 5. Control algorithm for the FCC reactor-regenerator system + _ Select: A, F, T1, T2, l+ p , lc , lc where A is the combustion air ¯ow, F is the feed ¯ow, T1, T2 are the feed pre-heat and riser outlet temperatures respectively. To minimize:

ÿ ‡ ‡ ÿ ÿ ‡ ‡ J ˆ Wÿ p lp ‡ W lp ‡ W c lc ‡ W c lc ‡
ui…iˆ1;2;3;4†

…11† where W are weighting matrices with elements such that wii 50, wij ˆ 0, i 6ˆ j and i and j denote weighting factors on the ith and jth slack variables. ui are the manipulated variables as de®ned above. In addition, the positive ‡ values of lÿ p or lp represent the di€erence between output response and the GMC reference trajectory. This di€erence is a measurement of the control performance ‡ degradation. Positive values of lÿ c or lc measure the di€erences between the actual rates of changes of the constraints beyond the speci®ed maximum rates of approach towards the constraint bounds. Larger values indicate that the constraints may exceed their bounds in the future. Thus they are the measurement of potential constraint violation. All of these l's are desired to be minimized. The slack variable weighting factors, W, re¯ect the relative importance of the outputs and constraints. Thus, if the constraint control is more important than the quality control, then: ÿ ÿ ÿ W‡ c ; W c >> W p ; W p

…12†

Conversely, if the quality control is more important than the constraint control, then: ÿ ‡ ÿ W‡ p ; W p >> W c ; W c

…13†

This indicates that the functional constraints virtually have no e€ect on the closed-loop responses of the system. Subject to: …t ÿ
ÿ
dy1 ÿ ‡ l‡ ÿ l ˆ k y ÿ y ‡ k y1sp ÿ y dt 11 1sp 21 p1 p1 dt 0 …14† equations for ¯ue gas oxygen concentration

Table 1 Comparison of model prediction with the plant results Variables Riser outlet temperature Regen bed temperature Air ¯ow rate Feed ¯ow rate Oxygen in ¯ue gas

Prediction Plant Prediction Plant Prediction Plant Prediction Plant Prediction Plant

Units

Case 1

Case 2

Case 3

Relative error



515 517 705 708 145 147 100 102 1.2 1.3

520 521 710 712 148 150 105 107 1.5 1.6

530 529 720 725 150 155 110 111 1.8 1.9

0.3%

C C  C  C tons/h tons/h m3/h m3/h % vol. % vol. 

0.5% 2.0% 1.6% 6.8%

R.M. Ansari, M.O. Tade / Journal of Process Control 10 (2000) 539±555

Fig. 4. Air ¯ow pattern on the 5% change in the setpoint of the air ¯ow controller.

545

Fig. 6. Plant and model responses for the reactor riser temperature with changes in the air ¯ow rate.

equation for wet gas compressor suction pressure output dy5 ‡ l‡ c 4K5C …y5U ÿ y5 † dt

…20†

equation for riser outlet temperature controller output AL 4A4AU …input constraint†

…21†


AL 4
A…t ‡
t† Fig. 5. Plant and model responses for the regenerator temperature with changes in the air ¯ow rate.

ÿ

ÿ ‡ ÿ l‡ p 50; lp 50; lc 50; lc 50




dy2 ÿ ‡ l‡ p2 ÿ lp2 ˆk12 y2sp ÿ y dt …t ‡ k22 …y2sp ÿ y†dt

…15†

0

equation for regenerator bed temperature dy2 ÿ lÿ c 4K1C …y2U ÿ y† dt

…16†

dy2 ‡ lÿ c 4K2C …y2 ÿ y2L † dt

…17†

equations for regenerator bed temperature constraints dy3 ‡ l‡ c 4K3C …y3U ÿ y3 † dt

…18†

equation for fuel gas ¯ow dy4 ‡ l‡ c 4K4C …y4U ÿ y4 † dt

ÿ
A…t†4
AU …input movement constraint†

…19†

…22† …23†

The inequality constraints on the controlled variables are the high and low limits placed on the system to represent the plant operating constraints. In all the FCC units operating in the full burn mode, coke combustion is essentially complete and increasing the ¯ow of combustion air to the regenerator can only raise the ¯ue gas oxygen concentration. All other controlled variables are controlled to one-sided zone limits. These limits set bounds on the units throughput and cracking severity. The optimization problem can be solved using a nonlinear constrained optimization algorithm. The form of the optimization problem is well structured since a slack variable is added to each control equation to ensure that a solution to the set of equation does exist. If the control is implemented at a reasonable frequency, the solution of the nonlinear programme (NLP) is very fast since the current control settings and the slack variables provide a good initial estimate of the solution vector. The above optimization problem was applied to the FCC reactor-regenerator section and the solution was obtained by using the Optimization toolbox in

546

R.M. Ansari, M.O. Tade / Journal of Process Control 10 (2000) 539±555

MATLAB. The toolbox contains many commands for the optimization of general linear and nonlinear functions. The optimization toolbox contains CONSTR.M ®le that ®nds the constrained minimum of a function of several variables. The program uses sequential quadratic program (SQP) method in which the search direction is the solution of a quadratic-programming problem. From the nonlinear control design (NCD) toolbox, a nonlinear optimization ®le NLINOPT.M was called to run the optimization algorithm. An alternative approach to solving the nonlinear optimization problem is to use the subroutine SOL/NPSOL described by Gill et al. [26] or to use any equation-based nonlinear optimization software. 6. Updating of the model parameters 6.1. System development In model-based controller, including GMC, there is an element of mismatch between the model and the true process. This process-model mismatch leads to deterioration in control performance. There are two types of model mismatch: structural mismatch occurs when the process and the model are of a di€erent nature (e.g. ®rst order/second order, or linear/nonlinear); parameteric mismatch occurs when the numerical values of the parameters in the model do not correspond to the true values. Lee et al. [18] presented a process-model mismatch compensation algorithm for model-based control. This algorithm compensated for model errors and updated the model parameters at steady state. Signal and Lee [27] have adopted a more practical approach by proposing the generic model adaptive techniques. When the system is minimum phase and in the absence of process model error, the closed loop response by GMC will follow the reference trajectory exactly. However, there are disturbances, which suddenly cause a di€erence between the output and the reference trajectory. The integral term in the GMC control law provides a compensation for the process-model mismatch. The control structure does not need to be changed if the closed loop response by this compensation is satisfactory. However, this method may not be applicable to all the systems because of the stability and the process constraints. It must be noted that the mathematical conditions for accurate tracking and o€set free control are much stronger than the speci®ed reference trajectories. When the mismatch becomes larger, the closed loop response will not be satisfactory. In the process industry, it is often required that the closed-loop response exhibit no overshoot and has a short rise-time. These properties are ensured when a perfect process model is used or the mismatch is small, but these are not guaranteed when the mismatch and any unmeasured

disturbances are large. The basic idea behind the development and the application of a dynamic parameter update system is to cope with these diculties and reduce the e€ect of large modelling errors by regularly updating the model parameters. Consider the process described by dx ˆ f…x; u; d; ; t† dt

…24†

y ˆ g…x; †

…25†

where x is the state vector of dimension m, u is the input vector of dimension n, d is the disturbance vector of dimension l, and y is the output vector of dimension n,  is the parameter vector of dimension k and t is the time. In general f and g are vectors of nonlinear functions. Eq. (24) is the general form of the model of a nonlinear process. Consider the model dx ^ ˆ f…x; u; d; ; t† dt

…26†

where f^ is an approximation to the true model function. Assuming that all states and outputs can be measured, then it follows that @g : y ˆ T f^…x; u; d; ; t† @x

…27†

The model form is given by: ÿ
@g ^ f…x; u; d; ; t† ˆ k1 ysp ÿ y ‡ k2 T @x

…t t0

ÿ


ysp ÿ y dt

…28†

The integral term in Eq. (28) eliminates o€set as discussed in detail by Lee et al. [18]. However, a more desirable approach is to eliminate one of the sources of the o€set, which is process-model mismatch. This would also be useful when process parameters and/or model structures show signi®cant time-dependence. Signal and Lee [27] derived an adaptive algorithm capable of adapting model parameters in a nonlinear model. The algorithm was developed within the GMC framework and reduces the e€ect of large modelling errors by regularly updating the model parameters. The e€ect of a change in the model parameters is related to the change in  given by a model parameter correlation as:
 ˆ 


…29†

where  is the di€erence between the model and the actual process and the matrix  is not necessarily square since the number of parameters which are being adapted can be less than or equal to the number of controlled

R.M. Ansari, M.O. Tade / Journal of Process Control 10 (2000) 539±555

variables. The terms
 and
 in Eq. (29) are de®ned as: ÿ

 ˆ L yt‡T ÿ ypt‡T

…30a†

L…i;i† ˆ ti

…30b†

where t+T is the sampling interval and L is the (kk) diagonal estimation reference trajectory matrix which determines the speed of the parameter estimation, ti is the speed parameter and
 is the di€erence between the predicted and the actual output variables over the sampling interval T. The
 is given by the equation: ÿ
ÿ1
 ˆ T W T W


…31†

We used a weighted least square approach to ®nd the optimal parameter change. The diagonal weighting matrix, W, re¯ects the range of uncertainty in measuring the outputs, y, and is given by:   1 1 …32† W ˆ diag 2 :::: 2 i k where i are the standard deviations of each of the measurements, yi . The expression (TW) in Eq. (31) is non-singular. In order to calculate the new parameters we equate
 from Eq. (30) to
 in Eq. (31) to obtain an expression for the dynamic model parameter update. ÿ
ÿ1  ÿ
 t‡T ˆ t ‡ T W T W L yt ÿ ypt

…33†

From the above, an algorithm for the dynamic parameter update procedure is summarized below and the procedure is given in Fig. 7 for the FCC reactor-regenerator section.

547

6.2. Parameter update algorithm Step 1 Choose the values of the GMC reference trajectory parameters for each controlled variables, calculate the matrices K1, K2, etc. Step 2 Choose ti , the speed parameters for the estimation of each parameter, i to be updated. Calculate the matrix L. Step 3 Estimate the error in the measurement and calculate the weighting matrix W, using Eq. (32). At each sampling time: Step 4 At time t measure the process. Calculate the right hand side of Eq. (28) and calculate the manipulated variable u from Eq. (27). Predict the output at time t+T. Step 5 At time t+T measure the output variable, yt+T, and calculate L (yÿyp) in Eq. (30). Step 6 Calculate , the dependence of the model on the parameters from Eq. (29). Step 7 Calculate the new model parameters from Eq. (33). Step 8 Go to Step 4. 6.3. Application to the FCC process In order to apply the dynamic parameter update program to estimate the Crc and Ccat from Eqs. (6) and (8), the normally distributed, zero mean noise (with standard deviation =0.2 and 0.4, respectively) was added to each of the measured controlled variables. Using the Eq. (32), the weighting matrix W, is given by:   25 0 Wˆ …34† 0 4

Fig. 7. System of dynamic parameter update for the FCC reactor-regenerator.

548

R.M. Ansari, M.O. Tade / Journal of Process Control 10 (2000) 539±555

From Eqs. (6) and (29), the concentration of carbon on the regenerated catalyst Crc , was estimated with the dynamic parameter update system by using a nonsquare  in the parameter estimation. The matrix  is given by:   ÿ0:03  kcc Vra Pra ÿEcc C rc exp …35† ˆ 0 Ccat C 0:03 R … T ‡ 460:0 † ra rc and from Eq. (31), the following expression is obtained: ÿ

T W


ÿ1

Ccat C 0:03 C 0:03 rc rc  ÿEcc kcc Vra Pra exp R…Tra ‡ 460:0†  ÿ0:06  Crc

T W ˆ

…36†

0 In this case, the weighting matrix W, does not a€ect the parameter estimation as the system consists of ®ve controlled variables and only one parameter is estimated in this example. Fig. 8 shows the trajectories of the concentration of carbon on regenerated catalyst. The actual values were plotted against the model's estimated values from the dynamic parameter update system. The sampling time T, was 2 s and L, the speed of parameter estimation was equal to 1. It was observed that while the model updates the parameter Crc , well within the region (0.65±0.70) wt. % of Crc , the estimation is not good below that region. This fact may be related to structural mismatch that can be made negligible if the model parameters are regularly updated. A similar approach can be applied to estimate and update the parameter Ccat . 7. Linear multivariable control The linear model predictive control algorithm applied here utilizes features from established model-based control algorithms such as DMC and MAC. The intent was to use the features of each algorithm as they best apply

to the problems associated with the FCC control. The main characteristics of the control algorithm are summarized here: 1. Discrete step response dynamic models 2. Reference trajectories to de®ne the desired closedloop response 3. Quadratic objective function solved for least squared error 4. Iterative control calculation for constrained input operation These methods are well established and ®eld proven; however, an overall control strategy requires additional programme functions, which are built around the basic MPC algorithm. This is especially true for the FCC control. Additional techniques employed to improve the FCC controls include: signal conditioning, prediction trend correction and control move calculation. The MATLAB toolbox on MPC was used for the linear control application together with the control techniques mentioned above. These techniques were incorporated into the programme and applied to the FCC reactorregenerator section. The results obtained were compared with the constrained nonlinear multivariable control application developed for the FCC process. 7.1. Process identi®cation test The scope of the process identi®cation test was to capture the plant dynamics in the models to develop linear controllers such as DMC. Process identi®cation tests were carried out by making step changes on each manipulated variable (u1, u2, u3 and u4). The step tests for each manipulated variable lasted 2±6 h during which time the variable was repeatedly stepped about its nominal value. In each case, the magnitude of the steps was selected to ensure that a clear response was observed on the process. The transfer function models of the FCC process are given in Table 2. A similar approach was adopted by Grosdidier et al. [28] during process identi®cation tests on a FCC process. 8. Real-time implementation

Fig. 8. Model parameter update for the FCC-estimating Crc, the concentration of carbon on the regenerated catalyst.

There are three important steps to implement the multivariable control applications in real-time. The ®rst step is to follow the exact procedure and sequence of implementation. Fig. 9 illustrates the implementation procedure of the FCC nonlinear model-based multivariable controller. The second step is to con®gure the controller interface to distributed control system and ®nally the on-line controller tuning to meet the speci®ed targets without process deviation and to gain the con®dence of the operators.

R.M. Ansari, M.O. Tade / Journal of Process Control 10 (2000) 539±555

549

Table 2 The FCC process models obtained by process identi®cation test and the FCC process control objectives

Flue gas oxygen conc.

Regen. bed temp.

Combustion air ¯ow (u1)

Feed ¯ow rate (u2)

Feed preheat temp. (u3)

Riser outlet temp. (u4)

Control objective

0:100…1:7s ‡ 1†eÿ2s 18s2 ‡ 7:0s ‡ 1

ÿ0:90 eÿ2s 13s2 ‡ 4:6s ‡ 1

0:030eÿ7s 12s ‡ 1

ÿ0:080…4:8s ‡ 1† 9:0s2 ‡ 3:0s ‡ 1

y1 (setpoint)

0:08 eÿ4s 11s2 ‡ 8:0s ‡ 1

0:45 eÿ4s 23s2 ‡ 8:0s ‡ 1

0:010 eÿ7s 10s ‡ 1

0:80…1:7s ‡ 1†eÿ2s 10s2 ‡ 7:3s ‡ 1

y2 (zone limits)

Fuel gas ¯ow

±

0:18 eÿ11s 40s2 ‡ 8:5s ‡ 1

±

0:37…16s ‡ 1† 50s2 ‡ 20s ‡ 1

y3 (max. zone limit)

WG comp. suc. pres. controller output

±

0:30 eÿ11s 16s2 ‡ 7:0s ‡ 1

±

0:60 3:0s ‡ 1

y4 (max. zone limit)

Riser out temp. controller output

±

0:72 eÿs 2:5s ‡ 1

1:5 2:0s ‡ 1

y5 (max. zone limit)

Fig. 9. Flowchart illustrating the implementation procedure of the FCC nonlinear model based multivariable controller.

8.1. Nonlinear controller tuning The time constants of the temperature controllers residing in a distributed control system (DCS) were estimated from open-loop tests. The time to steady state was found to be about 100 min. The sampling frequency was taken once per minute. The desired trajectories of the process outputs were determined using the tuning rules de®ned by Lee and Sullivan [16]. The control application was ®rst turned on and tested in an openloop mode (i.e. the outputs were suppressed). The open loop testing includes verifying process inputs and checking controller move sizes and directions. Next, the controls were tested in closed-loop mode using tuning

ÿ0:8 eÿ10s 6:0s ‡ 1

constants gathered from the simulation model and step response test. The closed-loop testing was performed within a narrow operating range until con®dence was gained. It was also noted that the form of the optimization problem is well structured since a slack variable was added to each control law equation to ensure that a solution to the set of equations does in fact exist. It is important to note that if the control is implemented at a reasonable frequency, the solution of the NLP is very fast (3 to 4 iterations), since the current control settings and slack variables provide a good initial estimate of the solution vector. Further, when the control system fails or turns itself o€, the basic control system continues to function. In e€ect, the operator does not have to do anything special to put the reactor-regenerator on computer or to take it o€. Two types of tuning parameters were used; the K1C and K2C, specifying the maximum speed of approach towards the constraint bounds, and weighting factors W, re¯ecting the relative importance of the outputs and the constraints. The tuning parameters K1C and K2C can be compared with the move suppression factor, a term used in linear multivariable controller design such as dynamic matrix control (DMC). Changing the move suppression factor on a manipulated variable causes the controller to change the balance between movement of that manipulated variable, and errors in the dependent variables. The role of the weighting factors, W, is similar to the correction factor in a linear multivariable control framework such as DMC. In the regenerator loading controller tuning analysis, the weighting factor for the regenerator temperature (controlled variable) was assigned a value of 1 C. It was observed that 1 C on the regenerator bed temperature was as important as 2 tons/h air ¯ow. If the controller was unable to reach all the setpoints or limits, it would prioritize the errors in the ratio 1 C=2 tons/h.

550

R.M. Ansari, M.O. Tade / Journal of Process Control 10 (2000) 539±555

8.2. Results and discussions The transfer functions of the FCC process models obtained through the process identi®cation tests are given in Table 2, where the process control objectives of each variable is also de®ned. The controller performance was checked for the unconstrained case with minimum manipulated variable move j
MV j=2.75%. The controlled variable response is shown in Fig. 10. This ®gure shows that while the O2 concentration and the riser outlet temperature controller output (RxTout) are least a€ected with the manipulated variable move, the regenerator temperature (RegT) rises signi®cantly.

Fig. 11. shows the response of the controlled variables when process variables are constrained to ‹5%. The O2 and riser outlet temperature controller output are still not a€ected with the manipulated variables move, min. j
MV j=8.50%. The regenerator temperature and fuel gas (FGas) values remain within the constraint limits. This demonstrates the constraint handling capabilities of the multivariable control. Fig. 12 shows the response of the controlled variables when both the process variables and manipulated variables are constrained. The minimum manipulated variable move, min. j
MV j=6.23%. A reduction in the feed resulted in a reduction of the fuel gas, therefore

Fig. 10. Controlled and manipulated variables for the unconstrained case.

Fig. 11. Controlled and manipulated variables for the constrained case (process variables constrained to ‹5).

R.M. Ansari, M.O. Tade / Journal of Process Control 10 (2000) 539±555

lower conversion and less e€ect on the O2 concentration and the riser temperature controller output. This fact may also be observed in the operation of the FCC unit. 8.3. Comparison of nonlinear control with DMC 8.3.1. Simulation results In Figs. 13, 14 and 15, linear control performance was compared with the nonlinear control. The setpoint, the operating limits of the controlled variables and the GMC parameters for the nonlinear controller are given in Table 3. For each output variable in Table 3, the GMC controller parameters K1i and K2i (the ith diagonal elements of K11 and K21 ) from Eqs. (14) to (15) were calculated by the following relationships: K1i ˆ

2i i

…37†

K2i ˆ

1 i 2

…38†

i was estimated from the process identi®cation test (plant test runs) and i was selected from the GMC performance speci®cation curve given by Lee and Sullivan [16]. This pair of parameters has explicit e€ects on the closed-loop response of the output variable. The parameter i , determine the shape of the closed-loop response while i , determine the speed of the response (large i implies slow response). Fig. 13. compares the DMC control with nonlinear control in terms of decoupling capability of the controllers. Since the FCC process is nonlinear and time-variant,

551

decoupling capability of the nonlinear controller provides better performance compared to the linear controller. It can be observed from Fig. 13 that the nonlinear controller displays a good decoupling capability. Consider the response of O2 concentration in Fig. 16. This response was obtained from the transfer function model in Table 2 between the riser outlet temperature (u4) and the ¯ue gas oxygen concentration (y1). The oxygen concentration shows the presence of nonminimum-phase characteristics due to the right-half plane (RHP) zeros and poses dicult control problems in linear systems. For a linear system, the zeroes are the poles of the process inverse and are nonminimum phase if the process inverse is unstable. For a nonlinear system, there is no explicit de®nition for process zeroes and poles. Fig. 14 compares the performance of linear and nonlinear controllers in controlling the process with inverse response. This ®gure shows that tuning the move suppression factor in linear control to avoid inverting RHP zero is dicult. The nonlinear control gives a better performance in this regard. Fig. 15 compares the e€ect of tuning on the performance of linear and nonlinear controllers. This ®gure shows that under linear control small changes in the manipulated variable has a signi®cant impact on the O2 concentration and the regenerator temperature. The control action is much more aggressive and instability can be observed for large changes in the manipulated variable moves. However, the nonlinear control performance under the same conditions is far better and smooth. This also demonstrates that the tuning of the nonlinear controller is much easier than that for the linear controller.

Fig. 12. Constrained controlled and manipulated variables.

552

R.M. Ansari, M.O. Tade / Journal of Process Control 10 (2000) 539±555

Fig. 13. Comparison of the decoupling capability of linear and nonlinear controllers.

Fig. 14. Comparison of linear and nonlinear controllers in controlling the process with inverse response.

R.M. Ansari, M.O. Tade / Journal of Process Control 10 (2000) 539±555

553

Fig. 15. E€ect of tuning parameters on the performance of linear and nonlinear controllers.

Table 3 Setpoint, operating limits and the GMC parameters Output

O2 conc. Reg.T RxTout F. gas

Setpoint/limits

1.0% 714±715 C 80% 15 tons/h

GMC parameters 

 (min)

5 5 5 5

70 70 70 70

8.3.2. Real-time implementation results Figs. 17 and 18 compare the nonlinear controller performance with DMC in real-time application on the FCC process. Fig. 17 shows the regenerator bed temperature control trends. The regenerator bed temperature is required to be controlled between the limits (714±715) C as shown in Table 3. The DMC controller shows that the temperature control sometimes exceeds the desired range, however, the nonlinear controller keeps the temperature well within the range. Both the controllers perform well, however, under nonlinear control, the regenerator bed temperature response is much more stable. This could be because the regenerator bed temperature is strongly a€ected by coke formation on the regenerator bed catalyst. The two parameters de®ned in

Fig. 16. Responses of the controlled variables to a step change in the riser outlet temperature (u4).

the coke formation rate Eqs. (6) and (8) are updated regularly by the parameter update system. This improves the performance of the nonlinear controller over DMC. Fig. 18 shows the oxygen concentration control trends for both controllers in real-time application. The primary control objective is to minimize the oxygen concentration in the ¯ue gas analysis to a setpoint value

554

R.M. Ansari, M.O. Tade / Journal of Process Control 10 (2000) 539±555

Fig. 17. Comparison of nonlinear control performance with DMC: real-time trends of the regenerator bed temperature control.

Fig. 18. Comparison of nonlinear control performance with DMC: real-time trends of the ¯ue gas oxygen concentration control.

of 1% vol. The DMC controller keeps the oxygen concentration between the limits (0.9±1.1)% vol. The response is far better under the nonlinear controller, since the deviation from the setpoint is very small (<0.1%). In Figs. 17 and 18, the prime objective of nonlinear control is to protect the regenerator temperature constraints. By reducing the variation in the regenerator bed temperature in Fig. 17 and keeping it within the constraint limits, nonlinear control helped to increase the life of the catalyst. Better control of oxygen in the ¯ue gas as shown in Fig. 18, also compensates for the changes in the feed rate, the recycle ¯ow rate, the riser outlet temperature and the feed temperature. Typical bene®ts here are 4±6% increases in the fresh feed rate to the unit against active constraints. 9. Conclusions The FCC unit is the most complex and challenging operating process in a modern re®nery. The success of any advanced control application (linear or nonlinear) on this unit depends on the ability to deal with the unit operating constraints and to address its economic

objectives which are de®ned by the economic optimal points for the process operation. A constrained nonlinear optimization strategy for handling the constraints has been developed and applied in real-time to the FCC reactor-regenerator section. The nonlinear control strategy presented in this paper translates operating constraints into setpoint and zone limit objectives for the controlled variables of the process. Because it is more important to satisfy operating constraints than to meet economic objectives, the nonlinear controller moves the manipulated variables in such a way that all the controlled variables remain at their setpoints or within their zone limits. In this manner, the process remains as close as possible to its optimal operating point while at the same time ensures that no operating constraint is violated. A simpli®ed form of nonlinear dynamic model from Denn [25] was used as a reference for the design of the control algorithm. The model results were compared with the plant data and used to improve the dynamic response of the model. The control problem of the FCC reactor-regenerator section was analyzed with an extension of the conventional MPC algorithm using MATLAB control features. The results of the linear MPC were compared with the nonlinear model-based multivariable control. It was shown that the nonlinear multivariable control provided better decoupling ability and control of the process with large deadtime compared to the traditional linear model-based control. It has also been demonstrated that the tuning of the nonlinear controller is much easier than the tuning of the linear controller. The same model was used for optimization and control, minimizing the modelling errors due to process/model mismatch. The main contribution of this work is to combine nonlinear process model with the nonlinear constrained optimization algorithm and to apply it to a highly nonlinear ¯uid catalytic cracking process. The MATLAB toolbox on MPC was used for the linear control application incorporating the control techniques in the program such as control variable compensation and prediction trend correction. For the nonlinear system, constrained nonlinear optimization problem was applied to the FCC reactor-regenerator section and the solution was obtained by using the optimization algorithm described in Section 5. The merit of this technique and the degree of ¯exibility in establishing the proper balance between the violation of the constraint variables and the deterioration of the control performance coupled with the ability to pre-de®ne the response trajectories for both control and constraints have been shown through the optimization formulation. The process model played a vital role in this strategy. It was used to make the outputs follow their trajectories as given in Eqs. (14) and (15) and to avoid the constraint violations of regenerated bed temperature in Eqs. (16) and (17). Therefore the accuracy of

R.M. Ansari, M.O. Tade / Journal of Process Control 10 (2000) 539±555

the model is very important for the success of the strategy. Although, the integral term in the GMC control law provided the compensation for the process-model mismatch, a dynamic model parameter update algorithm was developed for the FCC process to account for any large modelling errors by regularly updating the two important parameters in Eqs. (6) and (8). In real-time application, the constrained nonlinear optimization strategies were tested and compared to the DMC. The main advantage of the nonlinear control approach, besides the improved control performance, is that a single-step control law resulted in a smaller dimensional programme compared to the other methods available and the time spent to do the plant testing to obtain the process models was reduced signi®cantly. The other advantage of this approach is that it minimizes the maintenance and process re-identi®cation e€ort that is often required for linear controllers since the operating conditions change all the time. This is because the same model was used for both optimization and control. The implementation becomes easy if the operators are involved from the modelling and the plant identi®cation stages to the ®nal controller design and tuning. It is also important to ®rst turn the controller on, in an open loop mode, i.e. with the suppressed outputs and then test the controller in closed-loop within a narrow operating range. This procedure is safe and gives more con®dence to the operators while using these applications. Acknowledgements We are grateful to Dr. Weibiao Zhou of University of Western Sydney and Professor Peter Lee of Murdoch University, Australia for useful discussions and information on the GMC. We also acknowledge Ashraf Ghazzawi of Saudi Aramco who helped us to include the DMC features in the MATLAB system. References [1] H. Kurihara, Optimal Control of Fluid Catalytic Cracking Processes. ScD. Thesis, MIT 1967. [2] Y.Y. Zheng, Dynamic modelling and simulation of a catalytic cracking unit, Computers Chem. Engng. 18 (1) (1994) 39±44. [3] L.F. Moro, D. Odloak, Constrained multivariable control of ¯uid catalytic cracking converters, Journal of Process Control 5 (1) (1995) 29±39. [4] J.M. Caldwell, J.G. Dearwater, Model predictive control applied to FCC units. Fourth International Conference on Chemical Process Control. South Padre Island, Texas, 17±22 February, 1991. [5] C.R. Cutler, B.L. Ramaker, Dynamic matrix control Ð A computer control algorithm. Joint Automatic Control Conference Preprints, Paper WP5-13, San Francisco, 1980.

555

[6] R. Rouhani, R.K. Mehra, Model algorithmic control (MAC): Basic theoretical properties, Automatica 18 (4) (1982) 401±414. [7] T.S. Chang, D.E. Seborg, A linear programming approach for multivariable feedback control with inequality constraints, Int. J. Control 37 (3) (1983) 583±597. [8] D.L. Little, T.F. Edgar, Predictive control using constrained optimal control. Proceedings of the American Control Conference, Seattle, 3, pp. 1365±1371, 1986, 18±20 June. [9] C.E. Garcia, M. Morari, Internal model control. 1. A unifying review and some new results, Ind. Eng. Chem. Proc. Des. Dev. 21 (2) (1982) 308±323. [10] C.E. Garcia, M. Morari, Internal model control. 2. Design procedure for multivariable systems, Ind. Erg. Chem. Proc. Des. Dev. 24 (2) (1985a) 472±484. [11] C.E. Garcia, M. Morari, Internal model control. 3. Multivariable control law computation and tuning guidelines, Ind. Eng. Chem. Proc. Des. Dev. 24 (2) (1985b) 484±494. [12] C.B. Brosilow, G.Q. Zhao, K.C. Rao, A linear programming approach to constrained multivariable control. Proceedings of the third American control conference, San Diego, 2, pp. 667± 674, 1984, 6±8 June. [13] N.L. Ricker, Use of quadratic programming for constrained internal model control, Ind. Eng. Chem. Process Des. Dev. 24 (4) (1985) 925±936. [14] C.G. Economou, M. Morari, B.O. Passon, Internal model control Ð extension to nonlinear systems, Ind. Eng. Chem. Process Des. Dev. 25 (2) (1986) 403±411. [15] W.C. Li, L.T. Biegler, Process control strategies for constrained nonlinear systems, Ind. Eng. Chem. Res. 27 (8) (1988) 1421±1433. [16] P.L. Lee, G.R. Sullivan, Generic model control, Computers Chem. Engng. 12 (6) (1988) 573±580. [17] P.L. Lee, Direct use of nonlinear models for process control. Chemical Process Control Ð CPCIV (1991) 973±984. [18] P.L. Lee, G.R. Sullivan, W. Zhou, Process/model mismatch compensation for model-based controllers, Chem. Eng. Comm. 80 (1989) 33±51. [19] P.L. Lee, G.R. Sullivan, W. Zhou, A new multivariable deadtime control algorithm, Chem. Eng. Comm. 91 (1990) 49±63. [20] M.W. Brown, P.L. Lee, G.R. Sullivan, W. Zhou, A constrained nonlinear multivariable control algorithm, Trans IChemE 68 (1990) 464±476. [21] D.W. Prett, R.D. Gillette, Optimization and constrained multivariable control of a catalytic cracking unit. AlChE National Meeting, Houston, Texas, 1979. [22] G.D. Martin, J.D. Mahoney, H.C. Kliesch, Rigorous simulation used to determine FCC computer control strategies. Sixth Annual Fluid Catalytic Symposium, Munich, Germany, 1985, 21±23 May. [23] I. Huq, M. Morari, R.C. Sorensen, Modi®cation to model IV ¯uid catalytic cracking units to improve dynamic performance, AlChE Journal 41 (1995) 1481. [24] R.C. McFarlane, R.C. Reineman, J.F. Bartee, C. Georgakis, Dynamic simulator for a model IV ¯uid catalytic cracking unit, Computers Chem. Engng. 17 (1993) 275±285. [25] M.M. Denn, Process modelling, Longam Inc, New York, 1986. [26] P.E. Gill, W. Murray, M.A. Saunders, M.H. Wright, User's Guide for NPSOL (Version 4.0): A Fortran package for nonlinear programming. Technical Report SOL 86-2, Stanford University, Stanford, California, 1986. [27] P.D. Signal, P.L. Lee, Generic model adaptive control, Chem. Eng. Comm 115 (1992) 35±52. [28] P. Grosdidier, A. Mason, A. Aitolahti, P. Heinonen, V. Vanhamaki, FCC unit reactor-regenerator control, Computers Chem. Engng. 17 (2) (1993) 165±179.