1834
Ind. Eng. Chem. Res. 1989,28, 1834-1845
Abbreviations
DMC = dynamic matrix controller MIMC = multivariable internal model control SISO = single-input single-output SLP = sectional linear programming Literature Cited Bamberger, W.; Isermann, R. Adaptive On-line Steady-State Optimization of Slow Dynamic Processes. Automatica 1978, 14, 223-230. Bhattacharya, A,; Joseph, B. On-line Optimization of Chemical Processes. Presented at the First American Control Conference, Arlington, VA, 1982; paper MP5. Garcia, C. E.; Morari, M. Optimal Operation of Integrated Processing Systems. Part I: Open-loop On-line Optimizing Control. AIChE J. 1981,27,960-968. Garcia, C. E.; Morari, M. Optimal Operation of Integrated Processing Systems. Part 11: Closed-loop On-line Optimizing Control. AIChE J. 1984, 30, 226-234. Garcia, C. E.; Morari, M. Internal Model Control 2. Design Procedures for Multivariable Systems. Ind. Eng. Chem. Process Des. Deu. 1985a, 24, 472-484. Garcia, C. E.; Morari, M. Internal Model Control 3. Multivariable Control Law Computation and Tuning Guidelines. Ind. Eng. Chem. Process Des. Dev. 1985b,24, 484-494.
Griffith, R. E.; Stewart, R. A. A Nonlinear Programming Technique for the Optimization of Continuous Processing Systems. Manage. S C ~1961, . 7, 379-392. Gustavsson, I.; Ljung, L.; Soderstrom, T. Identification of Processes in Closed Loop-Identifiabilityand Accuracy Aspects. Automatica 1977, 13, 59-75. Lee, K. S.; Lee, W. On-line Optimizing Control of a Nonadiabatic Fixed Bed Reactor. AZChE J. 1985, 31, 667-675. Ljung, L.; Soderstrom, T. Theory and Practice of Recursive Estimation; MIT Press: Cambridge, MA, 1983. McFarlane, R. C.; Bacon, D. W. Adaptive Optimizing Control of Multivariable Constrained Chemical Processes. Part 2: Application Studies. Znd. Eng. Chem. Res. 1989, following paper in this issue. Morari, M.; Arkun, Y.; Stephanopolous, G. Studies in the Synthesis of Control Structures for Chemical Processes. Part I: Formulation of the Problem. Decomposition and the Classification of the Control Tasks. AIChE J . 1980, 26, 220-232. Prett, D. M.; Gillette, R. D. Optimization and Constrained Multivariable Control of a Catalytic Cracking Unit. Presented at the AIChE National Meeting, Houston, TX, 1979; paper WP5-C. Received for review October 5, 1988 Revised manuscript received July 11, 1989 Accepted August 15, 1989
Adaptive Optimizing Control of Multivariable Constrained Chemical Processes. 2. Application Studies Randall C. McFarlanet and David W. Bacon* Department of Chemical Engineering, Queen’s University, Kingston, Ontario, Canada K7L 3N6
The adaptive optimizing control strategy described in an accompanying paper (part 1)was successfully applied to simulations of two highly nonlinear, multivariable interacting processes. These simulated applications demonstrated the feasibility of using on-line identification of multivariable dynamic input/output behavior as a basis for employing sectional linear programming for on-line optimization. In addition, dynamic information from the identified input/output model was used for on-line updating of a control model for a multivariable internal model controller. Adaptive on-line identification of the system model ensured that the effects of model error on control and optimization performance were minimized.
As described in an accompanying paper (McFarlane and Bacon, 1989), chemical processes that are candidates for optimizing control often involve complex reaction and separation systems. For such systems, development of an optimizing control system based on a mechanistic process model may not be economically justified due to the expense of developing the model. In process control, the lack of adequate mechanistic process models has been overcome by research that has focused on the use of empirical input/output models. The availability of the resulting control designs has led to a dramatic increase in the number of industrial applications. A similar focus on systematic development of empirical methods for on-line optimization and optimizing control might be expected to produce similar results. The primary drawback of empirical steady-state approaches to on-line optimization is slowness of optimizing speed. Optimization methods that utilize direct search approaches (e.g., those based on simplex pattern searches) require that the process reach steady state at each operating condition in the search pattern, resulting in slow optimizing speeds. Similar criticisms apply to approaches
* Author to whom correspondence should
be addressed. Current affiliation: Amoco Corp., Amoco Research Center, Naperville, IL 60566. 0888-5885f 8912628-l834$0l.50/0
based on response surface methodology (e.g., Biles and Swain, 1980; Box and Draper, 1987). While these methods could be applied on-line, they would be useful only in situations where the movement in the location of the constrained optimum was sufficiently slow to allow it to be tracked by a relatively slow optimizer. However, the need for on-line optimization is more often motivated by a rapidly shifting optimum requiring faster optimizing speeds. The optimizing control strategy developed in our accompanying paper (McFarlane and Bacon, 1989) utilizes on-line identification of dynamic empirical process models to develop fast optimizing speeds. The steady-state gain of a dynamic input/output model provides information for optimization without actually requiring that the process ever reach steady state during the search. A related approach has been reported by Cutler and Hawkins (1987) in an application to a hydrotreater reactor. In that study, steady-state predictions from a step response model of the process provided information for on-line application of sectional linear programming. The step response parameters were identified off-line using test process data. No on-line updating of the process model was performed. The optimizing control strategy developed in the present study has a number of advantages over the approach described by Cutler and Hawkins. The use of a parsimonious
0 1989 American Chemical Society
Ind. Eng. Chem. Res., Vol. 28, No. 12, 1989 1835
FA Ft3
--i-]~qL~ :; xC
The object of steady-state optimization for the system is to maximize the yield of desired product P: production rate of P yield of P = total feed rate of A and B
X
=- xPFR
X
FA XG
+ FB
Figure 1. CSTR system.
where
input/output model, rather than a general step or impulse response model, leads to more efficient statistical identification of the process dynamics, which can be performed conveniently on-line. Furthermore, adaptive on-line identification and periodic updating of the control and optimization models can be expected to lead to more efficient optimizing control, since the effects of model error will be reduced. Application of the optimizing control strategy to two highly nonlinear, multivariable, and interacting processes is described in the following sections. The first system is a CSTR supporting a multiple reaction, and the second system is a fluid catalytic cracker. Descriptions of Simulated Systems CSTR System. The CSTR model is part of a larger test system introduced by Williams and Otto (1960). The Williams and Otto model has been used frequently in the literature for simulation testing of control and optimization systems (e.g., Di Bella and Stevens, 1965; Jung et al., 1971; Adelmann and Stevens, 1972; Davison and Chadha, 1972; Luus and Jaakola, 1973; Morari and Stephanopolous,1980; Rangaiah, 1985). In the present study, the CSTR system was simulated with different combinations of equality and inequality constraints on inputs and outputs in order to demonstrate the capabilities of the proposed optimizing control strategy. The CSTR system, shown in Figure 1, supports the following multiple reaction: kl
A+B-C ka
C+B-P+E kS
(4)
= xp
(1)
(2)
P+C-G (3) The desired product is P, while G, C, and E are byproducts subject to specified quality and environmental constraints, discussed below. Reactants A and B enter as pure components in separate streams with flow rates F A and FB, respectively. For the purposes of simulation, the feed rate F A is manipulatable while FB is considered to be a disturbance variable (it can be type 1or 2 (Morari et al., 1980) depending on its assigned properties). Briefly, type 1 disturbances are nonstationary and create the need for steady-state on-line optimization, whereas type 2 disturbances have short-lived effects and are irrelevant to long-term optimization. A control system is available for regulation of the reaction temperature with setpoint, TR.It is assumed that the closed-loop dynamics of this temperature loop are sufficiently fast that the thermal dynamics of the system are at pseudo steady state with respect to the concentration dynamics. With this assumption, the thermal dynamics are not modeled. The equations describing the kinetic behavior of the above reactions and the dynamic mass balances for the CSTR represent a coupled set of nonlinear algebraic and ordinary differential equations. A description of these equations has been provided by Williams and Otto (1960).
FR =FA
+FB
subject to specified constraints, described below. Measured outputs are the mass fractions of the products, xp, xE, xc, and xG, and the inputs are the feed rate, F A , and the reaction temperature setpoint, TR. For all simulation runs on the CSTR system, the ~ a m pling interval for identification and the control interval were 0.02 h. An independent normally distributed error was added to each output with a standard deviation of 0.10 for xp and xE and 0.02 for xc and xG. To ensure system identifiability, an independent pseudorandom binary sequence (PRBS) signal with the following properties was added to each input: uprbs = (3.0, 3.0)T N,,= 4 (0.08 h)
N,,,= 2’ - 1 Fluid Catalytic Cracking System. Fluid catalytic cracking (FCC) has been recognized as a candidate for optimizing control (e.g., Davis et al., 1974; Webb et al., 1978; Prett and Gillette, 1979). It is a complex and highly nonlinear process that represents a challenging test of the multivariable adaptive features of the proposed optimizing control system. The FCC simulation model employed here is based on the mechanistic process model developed by Kurihara (1967). The Kurihara model described an FCC unit in which only partial conversion of CO to COPoccurs in the regenerator bed. In this mode of operation, afterburning of CO to C02 above the regenerator bed will lead to unacceptably high temperatures in the cyclones unless the level of excess O2 in the stream leaving the bed is constrained to a very low level. A summary of the model, which consists of a set of highly coupled nonlinear ordinary differential equations and nonlinear algebraic equations, is provided in the Appendix section. A complete derivation of the model was described by Kurihara (1967). A simplified schematic representation of the FCC system is shown in Figure 2. The feed is a gas oil stream which undergoes a cracking reaction for molecular weight reduction to more valuable products, used for blending into gasoline. Gas oil feed is passed through a preheater, contacted with hot regenerated catalyst in the riser, and then fed to the reactor where the catalyzed cracking reaction occurs. The reaction is endothermic, with the heat of reaction supplied by the hot regenerated catalyst. Carbon is deposited on the catalyst surface as a result of the cracking reaction so that continuous regeneration of the catalyst is required. Spent catalyst is passed to the regenerator where it is contacted with atmospheric air. Oxygen reacts with deposited carbon to produce mostly carbon monoxide and carbon dioxide in an exothermic reaction. As mentioned above, only partial combustion of CO to C 0 2 occurs in the regenerator bed. A profit function for the FCC system was described by Kurihara: p = 60& (D,Yppp + Yglpgl + Ycopco - pt3 (5)
1836 Ind. Eng. Chem. Res., Vol. 28, No. 12, 1989 REGENEPATOR
CRACKING
28.0 Product
Gas
20.0
hg Level Controller
lo3
Iz,o
1 =l\E
11.0
r -4.0
5.0
E
2.5
h$ x lo3 I = XG -2.5 Feed
s
-5.0
f
Fc f 01f
Atmospheric Air
12.0
Air Blower
8.0
Figure 2. Fluid catalytic cracking system
A model for estimating the yields Y,,, Ygl, and Yc0as functions of the overall feed conversion was also given by Kurihara. Various process constraints limit the throughput and conversion capacity of FCC units: upper limits on the reactor and regenerator temperatures, upper limit on the catalyst circulation rate, maximum rate of air supply to the regenerator (combustion and lift air blower capacity), and upper limit on the oxygen content in the flue gas (to prevent afterburning of CO to COz in the cyclones). In the current simulation, the reactor pressure is not modeled, and therefore, the capacity constraint on the wet gas compressor, downstream of the FCC main fractionator, is not considered. System constraints are discussed further below. The principal type 1 disturbance considered by the Kurihara model is the coke formation factor, FCf,of the feed, which varies with the quality of the feedstock. The effect of F,, on the positions and shapes of constraint boundaries was illustrated by Arkun and Stephanopolous (1980). At different levels of Fcb various combinations of constraints on the inputs and outputs become active, resulting in different locations of the constrained optimum of feed conversion. Another significant type 1disturbance is ambient temperature. A decrease in air density with increased ambient temperature reduces the capacity of the air blowers, thereby reducing the throughput capacity of the FCC unit by limiting the amount of coke on spent catalyst that can be burned off in the regenerator. For all simulation runs on the FCC system, an independent normally distributed error was added to each output with standard deviations of 1.5 K for TI, and TI, and 0.002 mol % for 0,.
Simulation Results for the CSTR System Recursive Identification and MIMC Controller Performance. As described in our accompanying paper (McFarlane and Bacon, 1989), an empirical multivariable dynamic input/output model can be written in the following general form: A ( q ) y ( t )= B ( q ) u ( t )+ c
where, for the CSTR system,
th
h: x lo3 u.o
+ e(t)
(6)
1 =
F
q0
-11.0 12.0
8.0
h i x lo3
4.0
1=rp -4.0
P
10
20
30
110
50
k
Figure 3. Comparison of real and, identified values of Markov parameters for the CSTR system. hy, i = x p , xc, zG,xE, j = FA, k = l, 2, ..., 50; (-) real; (0, 0, A) identified after 50, 100, and 150 time periods of recursive estimation, respectively.
An off-line analysis based on a sum of squares criterion indicated that no significant improvement in fit was obtained with values greater than n = 2 and m = 2, and therefore, these values were selected for the structural parameters of eq 6. A recursive least-squares (RLS) estimation was performed under closed-loop conditions (controlled outputs, nc, xG; inputs, FA, TR), with the RLS algorithm initialized by using Oi(0)= 0
i = 1-4
Pi(0) = 1051 (rows of the input/output model) X(t) =
1.0
t = 1, 2,
...
where Oi, Pi, and X ( t ) are as defined in eq 13, 18, and 19, respectively, of the accompanying paper (McFarlane and Bacon, 1989). Independent PRBS signals were added to the inputs with properties as given previously. The estimated values of the impulse response (Markov) parameters (see eq 20 in McFarlane and Bacon (1989)) obtained after specified periods of RLS estimation are compared to their corresponding real values in Figure 3 (Markov parameters relating outputs to FA) and Figure 4 (Markov parameters relating outputs to TR).In these figures, the real Markov parameter values, connected by a continuous line to facilitate the comparison, were obtained from step tests on the deterministic simulation model. Identified Markov parameters were calculated after
Ind. Eng. Chem. Res., Vol. 28, No. 12, 1989 1837
h$ x lo3
70.0
2.00
50.0
1.50
30,0
XG
t
1.00
I = XE 10.0
0.50
-10.0 7.0
5.0
h{ x lo3
3.0
I = XG
0 2.00
6
1.50
xc
1.0 -1.0
12.0
1.00
0.50 I
0
r
78.00
11.0
h$ x lo3
I
71.00 L
-9.0
TR
I = xc -12.0
611.00 51.00
r -20.0
q Xlo3
50.00
32.0
36.00
22.0
30.00
12.0
FA
24.00
I =xp
2.0
18.00
-8.0
12.00 0
10
30
20
90
50
k Figure 4. Comparison of real and,identified values of Markov parameters for the CSTR system. &', i = xp, x ~ x , ~ x , ~ j' , = TR,k = 1,2, ..., 50; (-) real; ( O , O , A) identified after 50,100, and 150 time periods of recursive estimation, respectively.
50, 100, and 150 time periods of recursive estimation. Even after only 50 periods of estimation, it can be seen that reasonable estimates of the Markov parameters are available. The four-output/two-input dynamic model with n = m = 2 contains 20 parameters to be estimated on-line compared to 400 Markov parameters (with N = 50). Direct estimation of the Markov parameters would clearly require considerably more data and computation than estimation of the inputfoutput model considered here. Although an off-line analysis was performed in this study for structural identification of the model, many industrial systems can be adequately represented by a model with assumed values of 2 or 3 for m and n. The performance of an MIMC controller based on estimated Markov parameters was then studied. A two-input (FAand TR)two-output (xc and x ~ MIMC ) controller was obtained by using the following set of tuning parameters and the Markov parameter estimates obtained after 150 time periods of recursive estimation: p h = 10; M I= M 2 = 3; GI = 501,1= 1-3; Gi = 401,1= 4* Gi = 201,1= 7-10; Bl = 101,1= 1-3; Bi = 51,1= 4-6; Bi = 21,1= 7-10, where p h is the optimization horizon, MI and M2 are input suppression parameters for inputs 1and 2, and Bl and GIare input and output weighting matrices. p h , Mi, G I ,and Bi are the primary tuning parameters in the formulation of a predictive optimization problem, the solution of which gives the MIMC control law. Stated briefly, the optimization minimizes the sum of squares of deviations between predicted and desired future output values by adjusting the inputs. The inputs are allowed to
0
1.0
2.0
3.0
U.0
time ( h ) Figure 5. Servo performance of MIMC for the CSTR system based on identified Markov parameters with ai = 0.3,i = 1,2. Step change in xzt from 0.66 to 1.16 at t = 1.0 h.
vary over Mi, i = 1,2, ...,p, time steps into the future. The optimization is performed over Phtime steps into the future. Bl and GIare time-varying matrices that allow different weightings to be placed on inputs and outputs in the optimization problem. The effects of these parameters on the performance and robustness of the controller are described by Garcia and Morari (1985). For practical implementation, a filter, F(q),is placed in the feedback path of the controller to improve the robustness of the controller to model error. Usually a diagonal exponential filter is used:
Increasing the values of ai improves the robustness, but more sluggish controller performance results. Some difficulty was experienced in finding a satisfactory set of weighting matrices Gi and Bi for the CSTR system, finally resulting in those listed above. Insufficient weighting on the inputs resulted in oscillatory and nearly unstable closed-loop behavior, exacerbated by smaller control intervals, suggesting the presence of a zero near the unit circle. In the end, considerable weighting on the inputs was used along with a filter F(q) with constants ai = 0.3, i = 1,2, resulting in moderately sluggish controller performance. The servo performance of this controller for a step change in xgt at t = 1 h is shown in Figure 5. The observed performance was very similar to that of an MIMC
1838 Ind. Eng. Chem. Res., Vol. 28, No. 12, 1989
.I,.;
3.00
3.00
2.50
2 . 50
2.00
2.00
1.50
1.50
1.00
1.00
U.50
4.50
q.00
4.05
3.50
3. 50
3.00
3.CO
2.50
2.50
72.00
72.35
t
y,.
65.00
rp
CA.
t
65,DC 58.0C
56.00 51.00
E l , 3c
4u.00
4U.0C
70.00
7ii.0:
6U.00
64.0:
56.00
5 8 . 0(
52.00
52.0[
t
4 6 . 0(
116.00 0
1.0
3.0
2.0
time
u. 0
0
I O
2. 0
controller based on the real Markov parameters of the system. Figure 6 shows the servo performance (for a step change in xgt) of the MIMC controller at the operating point (XG, x c ) = (3.2, 2.01, with the controller derived from Markov parameter values identified a t the operating point (XG, xc) = (0.81, 1.0). Model error results in extreme oscillatory behavior of the closed-loop system. Increasing the filter constants to cyi = 0.95, i = 1, 2, removed the oscillations, as shown in Figure 7, but at the cost of extremely sluggish servo performance and deteriorated decoupling. In optimizing control of a nonlinear system, a fixed parameter controller will experience increasing effects of model error as the optimization steps are taken further away from the operating region within which the control model was originally identified. With large enough filter constants, F(q) is capable of stabilizing a closed-loop system when significant model/system mismatch occurs (Garcia and Morari, 1985) but the necessary trade-off in performance may be unacceptable. The approach employed here is to update the control law on-line as optimization steps are taken based on adaptive identification of a system model. This explicit approach to adaptive control can be computationally burdensome when the control law is updated at each time step. However, in an optimizing control application, there is no need for such frequent updating of the control law; in this study, the control law was updated only once every Nedtime steps, recognizing that optimization steps are the major source of model error. Performance of the Optimizing Control Strategy.
113
rime ( h
ti)
Figure 6. Servo performance for the CSTR system at ( X G , x c ) = (3.2, 2.0) with MIMC based on true Markov parameter values obtained a t ( x G , x c ) = (0.81, 1.0) with ai = 0.3, i = 1, 2. Step change in x$ from 2.0 to 2.5 at t = 1.0 h.
? i
Figure 7. Servo performance for CSTR system with ai = 0.95, i = 1, 2.
Test runs on the CSTR system were performed with the system model (eq 6) and the MIMC controller as defined in the preceding section. For all runs, the first optimization step was taken after 150 discrete time steps (3 h) of RLS estimation. Thereafter, optimization steps were taken after each adaptation period, Ned = 75 (1.5 h). At the beginning of an adaptation period, X ( t ) was reinitialized using A, = 0.98, h(0) = 0.98, and = 1.0 (see eq 19 in McFarlane and Bacon (1989)). A t the end of each adaptation period, the MIMC control law and the local steady-state optimization model were updated based on current estimates of the system model parameters. Independent PRBS signals, described previously, were added to the system inputs to ensure closed-loop identifiability of the system model. Run 1. For this run, two controlled inequality con~ , one uncontrolled strained outputs, yIC= (xc x ~ ) and inequality constrained output, y : = xE, were defined, resulting in the following optimizing control problem:
Xu)
max xp I t
Lt
IC ,IC
subject to
xFt I X e m
XE
5
xEt I ~ p l "
FA
5 FTax
X
p
TR I TW" where xF" = 1.4, FZ" = 78.0, x E m = 9.5, TE- = 95.0, and xg" = 28.0.
Ind. Eng. Chem. Res., Vol. 28, No. 12, 1989 1839 30.00
12.00 10.00
.rp
T,
E4
25.00
. r y
.
U
.YE
8.00 6.00
15.00
9.00
10.00
94.00
8.00
78.00
6.00
'G
62.00
4.00
46.00
2.00
30.00
0
74.00
2.00
58.00
I . 50
xc
42.00
1.00
0.50
26.00
0
10.00 4.0
0
8.0
12.0
16.0
20.C
The local linear programming problem (see eq 33 in McFarlane and Bacon (1989)) was defined as follows:
subject to rxP
1"
F f -< Fkw -< F > l V
nplv
T P 5 TLw 5
.pet,lv
x~t,l,lv< p t J v
< X"t,u,h - c - c -< xset,lv C < X"t,u,lv G
G
x'v
< nu - c
< xu
- G
< E
E-
F" IF: Ti I T ;
where x t = 1.2, FX = 74.0, x z = 8.0, T k = 94.0, and x $ = 26.0. The region of validity of the local steady-state model (see eq 26 in McFarlane and Bacon (1989)) was based on the amplitudes of the PRBS signals added to the system inputs: AF, = 3.0, AT, = 3.0. Changes to setpoints defining optimization steps were restricted in size (see eq 41 and 42 in McFarlane and Bacon (1989)) with Axrt = 0.5 and A x f = 1.0. The initial feasible point for run 1was (xEt, xgt) = (0.42, 2.0). The dynamic history of the run is shown in Figure 8. Optimization occurred in unconstrained space from t = 3 to t = 9 h. In this region, the optimization policy was to increase xFt while maintaining xFt approximately constant. The requested setpoint changes were implemented by the controller by increasing both inputs FA and TR. In determining the setpoint changes, the optimizer ensured that predicted steady-state levels of input variables resulting from a calculated optimization step were physically
L 0
time ( h ) Figure 8. Dynamic history of input/output variables for run 1, CSTR system.
,gt,1,"
r
u.0
8.0
12.0
16.0
20.0
time'(h)
realizable by the controller (i.e., would not violate constraints). At t = 9 h, the buffer constraint for XE (the uncontrolled inequality constrained output) became active. Feasible search directions were generated by increasing xZt and decreasing xEt. These setpoint changes were implemented by the controller by increasing FA and decreasing TR. As the optimization proceeded along the buffer constraint boundary for xE, improvements in the objective xp were not as pronounced as they had been during the optimization in unconstrained space. A t t = 15 h, the buffer constraint for xc became active and no further improvements in the objective were possible, as a local constrained optimum had been located. The optimizer then correctly maintained constant set points for x c and xG. At t = 18 h, the PRBS signals were removed. Dynamic histories for selected parameter estimates for row 4 of eq 6 (corresponding to output xE) are shown in Figure 9. These estimates clearly show adaptation as the run proceeded. Although not shown here, parameter estimates in other rows of the system model demonstrated similar behavior. From Figure 9, it is apparent that the parameter estimates had not yet converged when the first optimization step was taken at t = 3 h. Nevertheless, these estimates still adequately predicted the steady-state input/output behavior of the system, and a successful optimization step resulted. Not all of the parameter adaptation activity illustrated in Figure 9 can be attributed to system nonlinearity. The structural form of eq 6 is expected to result in some overparameterizationfor many systems, resulting in a high correlation among some of the parameters. Estimates of these parameters might vary over significant ranges without appreciably decreasing the sum of squares function. Run 2. The constraints for this run were identical with those for run 1 except for the boundaries on TR: Tgu = 86.0; T i = 81.0. The initial feasible point for run 2 was (xFt, xg;) = (0.3, 1.3). The dynamic input/output history of this run is
1840 Ind. Eng. Chem. Res., Vol. 28, No. 12, 1989
xFt at a greater rate than during the period 5 h It I 10 h. This optimization policy was implemented by the controller by reducing T R away from its constraint boundary and further increasing FA. At t = 16 h, constraints on both X E and xc were active, defining a constrained optimum, and the optimizer held the operating point fixed. Run 3. The constraint structure for this run was identical with that for run 1 except that the inequality constraint on x c was replaced with a single equality constraint xFt = x T = 1.2. Therefore, for this run, there was one variable adjusted for optimization, the setpoint of the single controlled inequality constrained output, y:" = xG, one controlled equality constrained output, y i = xc, and one uncontrolled inequality constrained output, y: = x E . The initial feasible point for run 3 was xEt = 1.37 (and f7jt = 1.2). The dynamic input/output history for this run is shown in Figure 11. Due to the equality constraint on xc, xFt was held constant at xFt = x y for the entire run. Referring to Figure 11, the optimization policy was to increase xEt beginning at the first optimization step at t = 3 h. This policy was followed until the upper buffer constraint on xE became active at t = 10 h. This operating point defined a local constrained optimum, and no further optimization steps were taken. During this run the controller successfully maintained xc at its fixed setpoint with zero offset and with little disruption from the setpoint changes made to xc. In klruns, contriller model updating was successful; no deterioration of controller performance was observed as optimization steps were taken across the feasible operating region,
xEt and increasing "1 a44
-1.00
1
2.00 I
2'43
I '
-1.00
-2.00
2
1.20
242 -2.00
II
.2
"1 a41
-1.00
!
I
'
0
5.0
10.0
15.0
1
20.0
time ( h ) Figure 9. Estimates of parameters from row 4 of the CSTR system model.
shown in Figure 10. At t = 5 h, the upper buffer constraint boundary for TRwas approached. With this constraint active, the optimization Policy was to increase the setpoints for x c and X G . This resulted in a steady-state optimization path along TW, and only occasional transient violations of TF" occurred. The levels of F A and TR plotted in Figure 10 are those calculated by the controller. The actual levels of F A and TRapplied to the system were not allowed to exceed FT" and TFaX. At t = 10 h, the buffer constraint for xE became active. Feasible search directions were generated by decreasing
for the FCC System Performance of the MIMC Controller. The system model described by eq 6 was defined for the FCC system as follows:
Y ( t )=
TR-
94 00
x u
800
TR
6200
xG
46 00
: 2M) 0
3000
FnU
2w 150
r G
78 00
L
-
1
x -
t
C
X U
C 0 0
40
80
120
160
200
time ( h ) Figure 10. Dynamic history of input/output variables for run 2, CSTR system
0
40
80
120
time ( h l
160
200
Ind. Eng. Chem. Res., Vol. 28, No. 12, 1989 1841
~p
12.00
30.00
10.00
25.00
XE
8.00
20.00
6.00
15.00
4.00
10.00
94.00
8.00
78.00
6.00
U
TR
62.00
4.00
46.00
2.00
30.00
0
74.00 58.00
FA
XG
t
tX
2.00
A
1.50
42.00
Xc
26.00
1.00 0.50 0
10.00 0
4.0
8.0
12.0
0
16.0
4.0
12.0
8.0
16.0
time ( h )
time ( h )
Figure 11. Dynamic history of input/output variables for run 3, CSTR system.
for n = m = 2. The system model was identified adaptively by using RLS estimation with a sampling interval equal to 0.1 h and under closed-loop conditions (controlled outputs, yIC( t)= [T,,(t) T (t)lT; inputs, u ( t ) as defined above). An independent FRBS signal with the following properties was added to each input: uprbs= (2.0, 2.0)T; Ncl = 2 (0.2 h); N ~ = ,27 - 1. At the end of 100 time steps of recursive estimation, Markov parameters were calculated. The estimated Markov parameter values were in excellent agreement with the real values as shown in Figure 12. The identified Markov parameters were used as the basis for derivation of a two-input/twc+output (with inputs and outputs as given previously) MIMC controller with the following tuning parameters: Ph= 10; MI = M , = 3; Gi = 201,l = 1-5; GI = 101,1= 6-10; BI = 41, I = 1,2, ..., 10. A filter F(q) with ai = 0.3, i = 1,2, was employed in the performance evaluations of the controller. The servo performance of the controller for a step change in Tt: from 779.4 to 780.9 K at t = 2 h is shown in Figure 13. The controller forced T , to its new setpoint level smoothly with little disruption to Trg. The regulatory performance of the controller for load disturbances (step changes) in feed temperature Tu was also tested and found to be satisfactory. Performance of the Optimizing Control Strategy. All optimizing control runs on the FCC system were performed with a two-input/ two-output MIMC and system model as described in the previous section. The system model was identified adaptively using a variable forgetting factor A( t ) reinitialized at the beginning of each adaptation period using A, = 0.98, X(0) = 0.98, and Xu) = 1.0 (see eq 19 in McFarlane and Bacon (1989)). For all runs, the first optimization step was taken after 100 time steps of RLS estimation and thereafter with an adaptation period Nad= 25 (2.5 h). An independent PRBS signal, with properties as described in the previous section, was added to each system input to ensure closed-loop identifiability of the system model.
The optimizing control problem was defied as follows: McFarlane max P T=*T:
where P is defined in eq 5 subject to
T $ ITE-
R,, IRRm
T $ IT T
Rai IRZ-
0, Io
r
where Tg-' = 788, REm = 304, TZm = 960, R Y = 28.40, and 0;- = 0.25. The local LP problem (see eq 33 in McFarlane and Bacon (1989)) was defined as m a T~,T=:
subject to
PIv
1842 Ind. Eng. Chem. Res., Vol. 28, No. 12, 1989 0. I 0
h
-0.1
i=T 'g
-0.2 -0.3 0.6
h klJ
i = T,,
0.4
0.2
0 -0.2 0
5
10
15
20
25
30
20
25
30
k
(4 0.6 0.4
h,'i i=T
0.2 rg
0 -0.2
0.3 0.2
hkY i = Tra
o,l 0 -0.1
0
5
10
15
mained active, resulting in optimization steps in which both T:Zt and Tiit were increased. At approximately t = 36 h, the buffer constraint boundary for Trpwas encountered, representing a constrained optimum. No further optimization steps were taken, and at t = 40 h, the PRBS signals were removed. The constrained optimum had been located in approximately26 h, resulting in an improvement in the objective function of $1.6/min. Run 2. The optimizing control problem was modified for this run by reducing the saturation limit on the input Rai, with the resulting inequality constraint on Raigiven by Rai IRE"" = 27.0. The buffer boundary for Rai was defined as R:i = 26.5. The initial feasible point for run 2 was TEt = 773.9 and Tiit = 939.3. Referring to Figure 15, the first buffer constraint boundary to be encountered was that for Of at approximatelyt = 24 h. The search direction was modified in order to follow the buffer boundary of 0, until the buffer boundary of Raiwas encountered at approximately t = 32 h. Transient excursions of Ri exceeded this buffer boundary between t = 24 h and t = 32 h. During this period, optimization steps continued to be taken, which required the controller to increase the steady-state level of R ~ This . was the correct action because the predicted steady-state levels of Rd remained below the level defining its buffer boundary. Excessive saturation of Rai would indicate that the buffer boundary should be decreased. The levels of Ri plotted in Figure 15 are those calculated by the controller. Actual levels of Raiapplied to the system were not allowed to exceed the saturation level, RZ".
k
(b)
Figure 12. Comparison of real,,and identified Markov parameter values for the FCC system. (a) h f , i = T,(t), T9(t)?J = ul(t) = R,(t), k = 1, 2, ..., 30; (-) real; (0)identified. (b) h f , z = Tra(t),T,,(t), j = u z ( t ) = 10Rai(t),k = 1, 2, ..., 30; (-) real; (0) identified.
where T &= 786, TYg= 958, R:c = 300, R:i = 28.0, and O g = 0.20. The region of validity of the local steady-state model was defined by the amplitudes of the PRBS signals: AR,, = 2.0, ARai = 0.2. Changes in setpoints defining an optimization step were restricted in size (see eq 41 and 42 in = 1.5 and McFarlane and Bacon (1989)) with = 1.5. Run 1. The initial feasible point for run 1 was T:Zt = 775 and T:it = 943. The dynamic history for this run is shown in Figure 14. The first optimization step was taken at t = 10 h. Between t = 10 h and t = 17.4 h, the optimizer called for increases in Tiit while T;it was held approximately constant. At t = 17.5 h, the search direction was modified since the operating point had approached the buffer constraint boundary for 0,. This constraint re-
Conclusions The empirical optimizing control strategy developed in the accompanying paper (McFarlaneand Bacon, 1989) was successfully applied to two highly nonlinear and multivariable constrained systems. The simulations demonstrated the utility of closed-loop on-line identification of empirical dynamic input/output models as a basis for on-line application of sectional linear programming. The fast optimizing speed was achieved with this dynamic identification approach. Model accuracy was maintained by adaptive on-line identification of the process model and periodic updating of the control and optimization models. The proposed optimizing control strategy is distinctive in that dynamic information for control and steady-state information for optimization are obtained from a single dynamic process model. The simulations demonstrated the feasibility of this approach to optimizing control. Protection against transient constraint violation due to disturbances d 2 is provided in the proposed optimizing control strategy by defining buffer constraint boundaries at specified distances from the actual inequality constraint
21 00
,
957.00
2688
I
956.00
2 b 7635 2h
- 1
Trg
-
-
955.00
-
941.00
-
953.00
-
781.00
Tra
780.00 77900
.
-
77800 0
10
20
30
40
50
60
time ( h )
Figure 13. Servo performance of MIMC for the FCC system based on identified Markov parameters. Step change in Tzt from 779.4 to 780.9 K a t t = 2 h.
Ind. Eng. Chem. Res., Vol. 28, No. 12, 1989 1843
P
53.00
0.40
52.00
0.30
Ofg
51.00
0
49.00
max
29.00
970.00
28.00
962.00
Trg
Rai 27.00
954.00 946.00
26.00 25.00
938.00
310.00
792.00
r
R
IF
r
r m
786.00
M0.m
Rlc
0.20
0.10
50.00
T,,
290.00
780.00
280.00
714.00
nom
768.00
0
6.0
120
18.0
24.0
30.0
36.0
42.0
I 0
48.0
6.0
120
18.0
time(h)
24.0
30.0
36.0
42.0
48.0
time ( h )
Figure 14. Dynamic history of input/output variables for run 1, FCC system. 53.00
r
0.40
0.30
5200
p
ofg
51.00 50.00 49.00
966.00 958.00
Trg
26.00
950.00 942.00
25.00
I
934.00 792.00
310.00
m.OO
Rrc
fs
0
27.00
24.00
fg O U
0.10
t
28.00
Rai
0-
0.m
786.00
T,,
290.00
280.00
780.W 774.00
270.00
768.00
0
6.0
120
18.0
24.0
m.0
36.0
420
0
48.0
6.0
120
18.0
24.0
30.0
36.0
420
48.0
time(h)
time(h)
Figure 15. Dynamic history of input/output variables for run 2, FCC system.
boundaries. The use of a constrained multivariable predictive controller, rather than the unconstrained MIMC controller used in the present study, would provide further protection against constraint violations and allow buffer boundaries to be placed closer to the actual constraint boundaries. Constrained controllers of this type have been described, for example, by Garcia and Morshedi (19861, Brosilow and Zhao (1988), and McFarlane and Reineman (1990). Many industrial systems can be adequately represented by the general dynamic system model expressed by eq 6. If the noise sequence le ( t ) )is not highly correlated, nonstochastic controllers such as Internal Model Control and Dynamic Matrix Control (Cutler and Ramaker, 1980) can be expected to provide adequate control performance. But eq 6 would not provide an adequate representation for systems affected by disturbances which are highly autocorrelated. For such systems, the controller design should include identification of a disturbance model in addition to a model describing the deterministic system behavior. One such design is the class of linear quadratic controllers based on transfer function models and Wiener-Hopf design principles (e.g., Harris and MacGregor, 1987). An adaptive version of this class of multivariable controllers which should be useful in optimizing control applications has been described by McFarlane (1986).
Acknowledgment Financial support from the Natural Sciences and Engineering Research Council of Canada is gratefully acknowledged.
Nomenclature a = parameter in the system model A(q) = matrix polynomial in the system model B = time varying weighting matrix B(q) = matrix polynomial in the system model c = constant vector in the system model C1 = empirical constant in eq A-13, value = 2.42 kg of oxygen/kg of coke C2 = empirical constant in eq A-15, value = 6.34 X lo4 s/(kg of air-kg of catalysbkPA) C3 = empirical constant in eq A-16, value = 1.16 X kg of oxygen/(s-kPa-kgof catalyst) Ccat= concentration of catalytic carbon on catalyst C,, = concentration of carbon on regenerated catalyst C, = concentration of carbon on spent catalyst Ctf = feed conversion, volume fraction Dtf = feed density, value = 892 kg/m3 e = random error vector in the system model FCf= coke formation factor of the feed, value = 0 (kg of carbon/s)/ (m3/s) Fi= flow rate of component i, kg/s
1844 Ind. Eng. Chem. Res., Vol. 28, No. 12, 1989
F(q) = filter transfer function in MIMC G = time varying weighting matrix h = impulse response parameter Hra= cracking reactor catalyst holdup, value = 1.85 X lo4kg HrB= regenerator catalyst holdup, value = 1.53 x lo5 kg I = identity matrix k = reaction rate constant in the CSTR model, s-l kcc = constant in rate expression A-10, value = 2.66 X loT4 kPa-'.s-' kcr = constant in rate expression A-6, value = 9.32 X m3/(kPa.kgs) Kcc = reaction rate constant in eq A-8, k P a - W K,, = reaction rate constant in eq A-5, m3/(kPa.kgs) K d = diffusion rate constant for oxygen in eq A-14, kPa-'.s-' KO,= reaction rate constant for oxygen in eq A-14, kg of oxygen/ (s.kPa.kg of catalyst) m = degree of B(q) M, = input suppression parameter for input i in MIMC rz = degree of A(q) N = number of parameters for each input/output pair in the impulse response model Nad = number of discrete time steps defining an adaptation period Ncl = minimum number of discrete time steps between changes in level of PRBS signal Nper= PRBS sequence length Of, = concentration of oxygen in the flue gas, mol % p = number of inputs P = objective function P = variance-covariance matrix of parameter estimates Pco = cycle oil value, value = $129/m3 Pgl= gasoline value, value = $173/m3 PBs= gas value, value = $O.l48/kg Ph = optimization horizon in MIMC P,, = cracking reactor pressure, value = 211.7 kPa Prg= regenerator pressure, value = 254.4 kPa Ptf= feed value, value = $118/m3 R = ideal gas constant, 8.314 J/(mol.K) Rad = rate of additive carbon deposition, kg/s Re = air flow rate, kg/s Rcb = rate of carbon burning, kg/s R, = rate of catalytic carbon deposition, kg/s Rcf = total carbon forming rate, kg/s R, = gas oil cracking rate, kg/s R,, = catalyst recirculation rate, kg/s R,, = spent catalyst flow rate, kg/s Ru = feed rate, value = 0.046 m3/s Sa = specific heat of air, value = 1130 J/(kgK) S, = specific heat of catalyst, value = 1047 J/(kg.K) Sf = specific heat of feed, value = 3140 J/(kgK) t = time, h T = temperature, K Te = air temperature, value = 394 K TI, = cracking reactor temperature, K T,, = regenerator temperature, K Ttf = feed temperature, value = 492.8 K u = vector of input variables uprbs = vector of amplitudes of PRBS signals V = reactor volume, value = 2.63 m3 x , = mass fraction of component i, kg of i/kg of mixture (scaled: xc X IO, T E X 100, X G X 100, xp X 100) y = vector of output variables y* = vector of controlled equality constrained outputs y 1 = vector of uncontrolled inequality constrained outputs yIC= vector of controlled inequality constrained outputs Y,, = cycle oil yield, volume fraction Ygl = gasoline yield, volume fraction Y,, = gas yield, weight fraction
AEm = activation energy of catalytic carbon formation, value = 4.18 X lo4 J/mol
AE,,= activation energy of catalytic cracking, value = 6.28 x lo4 J/mol
mor= activation energy of oxygen reaction, value = 1.47 X
lo5 J/mol AHcr= heat of cracking, value = 4.65 X lo5 J/kg AHh = heat of vaporization of feed, value = 1.74 X 105J/kg AH,,= heat of coke burning, value = 3.02 X lo7 J/kg 0 = vector of parameters X = forgetting factor in recursive estimation Subscripts ss = steady-state value R = reactor
Superscripts 1 = lower buffer constraint boundary lv = locally valid model max = absolute constraint boundary prbs = pseudorandom binary sequence set = setpoint spec = specification value T = transpose u = upper buffer constraint boundary
Abbreviations FCC = fluid catalytic cracker MIMC = multivariable internal model control PRBS = pseudorandom binary sequence RLS = recursive least squares SISO = single-input/single-output
Appendix Summary of the Fluid Catalytic Cracking Model. The following is a summary of the FCC model equations derived by Kurihara (1967). The simulations in the current study were run at the operating conditions described by Schuldt and Smith (1971) for a 0.046 m3/s (25000 bbl/day) unit. Reactor Model. energy balance
(A-1) total carbon balance
catalytic carbon balance
Rate Expressions. gas oil cracking rate Roc
= Dtf Rtf Ctf
Ctf -1 - c,
Greek Symbols a
= parameter in the filter transfer function
(A-4)
where
total carbon forming rate
- Kcr P r a H r a Rtf
(-4-5)
Ind. Eng. Chem. Res., Vol. 28, No. 12, 1989 1845
(A-7) where
Regenerator Model. energy balance
Brosilow, C. B.; Zhao, G. Q. A Linear Programming Approach to Constrained Multivariable Process Control. Control Dynam. Syst.: Adv. Theory Appl. 1988,27, 141-180. Cutler, C. R.; Hawkins, R. B. Constrained Multivariable Control of a Hydrocracker Reactor. Presented at the American Control Conference, Minneapolis, MN, 1987; paper TA7. Cutler, C. R.; Ramaker, B. L. Dynamic Matrix Control-A Computer Control Algorithm. Presented at the 86th AIChE National Meeting, Houston, TX, 1980; paper 51b. Davis, T. A,; Griffin, D. E.; Webb, P. U. Cat Cracker Optimization and Control. Chem. Eng. Prog. 1974, 70,53-57. Davison, E. J.; Chadha, K. J. On the Control of a Large Chemical Plant by Using Modal Analysis. Automatica 1972, 8, 263-273. Di Bella, C. W.; Stevens, W. F. Process Optimization by Non-Linear Programming. Znd. Eng. Chem. Process Des. Dev. 1965,4,16-20. Garcia, C. E.; Morari, M. Internal Model Control. 3. Multivariable Control Law Computation and Tuning Guidelines. Ind. Eng. Chem. Process Des. Dev. 1985,24,484-494. Garcia, C. E.; Morshedi, A. M. Quadratic Programming Solution of Dynamic Matrix Control (QDMC). Chem. Eng. Commun. 1986, 46, 73-87.
Harris, T. J.; MacGregor, J. F. Design of Multivariable Linear Quadratic Controllers Using Transfer Functions. AZChE J . 1987,
carbon balance
(A-12) carbon burning rate (A-13)
(A-14)
(A-15)
33,1481-1506.
Jung, B. S.; Mirosh, W.; Ray, W. H. Large Scale Process Optimization Techniques Applied to Chemical and Petroleum Processes. Can. J . Chem. Eng. 1971,49,844-852. Kurihara, H. Optimal Control of Fluid Catalytic Cracking Processes. Sc.D. Dissertation, M.I.T., Cambridge, MA, 1967. Luus, R.; Jaakola, T. H. A Direct Approach to Optimization of a Complex System. AZChE J . 1973,19,645-648. McFarlane, R. C. Empirical Strategies for Optimizing Control of Chemical Engineering Systems. Ph.D. Dissertation, Queen’s University at Kingston, Ontario, 1986. McFarlane, R. C.; Bacon, D. W. Adaptive Optimizing Control of Multivariable Constrained Chemical Processes. 1. Theoretical Development. Znd. Eng. Chem. Res. 1989, preceding paper in this issue. McFarlane, R. C.; Reineman, R. C. MultivariableOptimizing Control of a Model IV Fluid Catalytic Cracker. Paper to be presented at the 1990 AIChE Spring National Meeting, Orlando, FL, March 18-22.
(A-16)
Hold-Up Regulation. It was assumed that a SISO controller was available for regulation of reactor holdup, H,,, using spent catalyst flow rate, R,, as the manipulated input and that this regulation resulted in negligible variation in Hra. Since regenerator holdup, Hrg,is related to Hraby dH,,/dt = -dHw/dt, it could also be assumed that Hrgwas regulated with negligible variation. Literature Cited Adelmann, A.; Stevens, W. F. Process Optimization by the Complex Method. AIChE J . 1972, 18, 20-24. Arkun, Y.; Stephanopolous, G. Studies in the Synthesis of Control Structures for Chemical Processes. Part IV: Design of Steady State Optimizing Control Structures for Chemical Process Units. AZChE J . 1980,26, 975-991. Biles, W. E.; Swain, J. J. Optimization and Industrial Experimentation; Wiley-Interscience: New York, 1980. Box, G. E. P.; Draper, N. R. Empirical Model-Building and Response Surfaces; Wiley: New York, 1987.
Morari, M.; Stephanopolous, G. Studies in the Synthesis of Control Structures for Chemical Processes. Part 11: Structural Aspects and the Synthesis of Alternate Feasible Control Schemes. AZChE J . 1980,26, 232-246. Morari, M.; Arkun, Y.; Stephanopolous, G. Studies in the Synthesis of Control Structures for Chemical Processes. Part I: Formulation of the Problem. Decomposition and the Classificationof the Control Tasks. AIChE J . 1980,26, 220-232. Prett, D. M.; Gillette, R. D. Optimization and Constrained Multivariable Control of a Catalytic Cracking Unit. Presented at AIChE 86th National Meeting, Houston, TX, 1979; paper WP5-C. Rangaiah, G. P. Studies in Constrained Optimization of Chemical Process Problems. Comput. Chem. Eng. 1985,9,395-404. Schuldt, S. B.; Smith, F. B. An Application of Quadratic Performance Synthesis Techniques to a Fluid Cat Cracker. Proc. Joint Autom. Control Conf. 1971, 270-276. Webb, P. U.; Lutter, B. E.; Hair, R. L. Dynamic Optimization of Fluid Cat Crackers. Chem. Eng. Prog. 1978, 74,72-75. Williams, T. J.; Otto, R. E. A Generalized Chemical Processing Model for the Investigation of Computer Control. AZEE Trans. 1960, 79, 458-473.
Received f o r review October 5 , 1988 Revised manuscript received July 11, 1989 Accepted August 15, 1989
