Ind. Eng. Chem. Res. 1989,28, 1812-1822
1812
analysis of mass-transfer characteristics will be discussed in a later publication. Acknowledgment The authors thank Kazuo Toyomoto of the Industrial Membrane Division of Asahi Chemical Industry Co., Ltd., for his help in providing the starting hollow fiber. Registry No. Co, 7440-48-4. Literature Cited Ang, C. H.; Garnett, J. L.; Levot, R.; Long, M. A. Polyfunctional Monomers as Additives for Enhancing the Radiation Copolymerization of Styrene with Polyethylene, Polypropylene, and PVC. J. Appl. Polym. Sci. 1982,27, 4893. Egawa, H.; Nonaka, T.; Maeda, H. Studies on Selective Adsorption Resins. XXI. Preparation and Properties of Macroreticular Chelating Ion Exchange Resins Containing Phosphoric Acid Groups. J. Appl. Polym. Sci. 1985,30, 3239. Idol, W. K.; Anderson, J. L. Effects of Adsorbed Polyelectrolytes on Convective Flow and Diffusion in Porous Membrane. J. Men.
Sci. 1986, 28, 269. Kalal, J.; Svec. F.; Marousek, V. Reactions of Epoxide Groups of Glycidyl Methacrylate Copolymers. J. Polym. Sci. 1974,47, 155. Okamoto, J. Radiation Synthesis of Functional Polymer. Radiat. Phys. Chem. 1987,29,469. Omichi, H.; Katakai, A.; Sugo, T.; Okamoto, J. A New Type of Amidoxime-Group-Containing Adsorbent for the Recovery of Ura-
nium from Seawater. Sep. Sci. Technol. 1985,20, 163. Saito, K.; Yamada, S.; Furusaki, S.; Sugo, T.; Okamoto, J. Characteristics of Uranium Adsorption by Amidoxime Membrane Synthesized by Radiation-Induced Graft Polymerization. J. Mem. Sci. 1987a, 34, 307. Saito, K.; Hori, T.; Furusaki, S.; Sugo, T.; Okamoto, J. Porous Amidoxime-Group-ContainingMembrane for the Recovery of Uranium from Seawater. Znd. Eng. Chem. Res. 1987b, 26, 1977. Saito, K.; Uezu, K.; Hori, T.; Furusaki, S.; Sugo, T.; Okamoto, J. Recovery of Uranium from Seawater Using Amidoxime Hollow Fibers. AZChE J. 1988, 34, 411. Saito, K.; Kaga, T.; Yamagishi, H.; Furusaki, S.; Sugo, T.; Okamoto, J. Phosphorylated Hollow Fibers Synthesized by Radiation Grafting and Cross-Linking. J. Mem. Sci. 1989,43, 131. Skiens, W. E. Hollow-Fiber Membranes. In Encyclopedia of Polymer Science and Technology 15; Mark, H. F., Gaylord, N. G., Bikales, N. M., Eds.; Interscience Publishers: New York, 1971. Sugo, T.; Saito, K. Development of Hollow Fiber Adsorbent for the Recovery of Uranium from Seawater. Maku 1988, 13, 272. Xu, Z.-L.; Wang, G.-H.; Wang, H.-I.; Cian, G.; Ni, M.-H. Radiation Induced Grafting of Styrene and Its Mixture with Divinylbenzene onto F46 Film. Radiat. Phys. Chem. 1983,22,939. Yamagishi, H.; Saito, K.; Furusaki, S.; Sugo, T.; Okamoto, J. Effect of Vapor- and Liquid-Phase Radiation Grafting on Water Permeability of Porous Hollow-Fiber Membrane. Nippon Kagaku Kaishi 1988, 212. Received for review November 21, 1988 Revised manuscript received August 7 , 1989 Accepted August 30, 1989
PROCESS ENGINEERING AND DESIGN Multistep Nonlinear Predictive Controller David D. Brengel a n d W a r r e n D. Seider* Department of Chemical Engineering, University of Pennsylvania, Philadelphia, Pennsylvania 19104
A new model predictive formulation to control MIMO nonlinear processes is presented. It directly extends the model predictive concepts that have been exploited for the control of linear systems (e.g., IMC and Q/DMC). This new algorithm uses a multistep predictor, unlike the algorithms by Economou et al. and Li and Biegler, which were limited to a single predictive step. This paper derives the multistep predictor through a linearization of the ordinary differential equations at several instants within a sampling interval, leading to recursive, algebraic equations that relate the predicted outputs to future and past values of the manipulated inputs. This multistep capability is shown to result in superior performance over single-step methods. Furthermore, this multistep algorithm does not require iterative calculations, as required by Economou et al. and is more efficient for each predictive step. The algorithm also easily handles constraints involving the state variables and manipulated inputs. The algorithm provides robust servo-control for the fermentation reactor modeled by Agrawal et al., and simulation examples illustrate improved controller performance as compared with the algorithm of Li and Biegler. In addition, the algorithm stabilizes unstable steady states, providing excellent servo-control and disturbance rejection for these processes. With process control receiving increased attention in the chemical process industries over the past decade, control algorithms and design techniques have become increasingly sophisticated. This is the result of several factors, including the availability of powerful, interactive workstations at low cost and more realistic process models. These workstations, which are often networked to mainframe
* Author to whom correspondence should be addressed.
computers that perform the steady-state optimization of increasingly complex process designs, execute improved control algorithms that permit the operation of processes closer to the intersection of their design constraints, in the vicinity of their economic optima. Modern control algorithms often utilize a dynamic model of the process in the control loop. Normally, the model is a straight-forward extension of the steady-state model prepared for design optimization. Its inclusion allows for
0888-588518912628-l812$01.50/00 1989 American Chemical Society
Ind. Eng. Chem. Res., Vol. 28, No. 12, 1989 1813 the more “intelligent” selection of the values of the manipulated inputs, for both servo-control and regulator control, than can be achieved with traditional analogue (e.g., PID) controllers, which are initially tuned by using a linearized process model but in operation have no ”knowledge” of the process dynamics. This advantage is greatest when operating a process in the vicinity of design and operating constraints. These controllers, often referred to as model predictive controllers (MPCs), predict the process outputs several sampling intervals in the future, as a function of the optimized values of the manipulated input for the next sampling interval and those proposed for future sampling intervals. Thus, they can detect an impending constraint violation and adjust the manipulated input to maintain the process at or drive it to the desired set point along a n optimal trajectory. As a consequence of its structure, a MPC is a feedforward controller for known process changes and a feedback controller for those unknown. Thus, MPCs can reject measured disturbances more rapidly than conventional controllers, by anticipating their impact on the process, and achieve set point changes efficiently, through their ability to predict an optimal sequence of manipulated input values to be implemented. The feedback element of a MPC compensates for the effects of unmeasured disturbances on the process outputs and deviations between the model outputs and those measured (Le., process/model mismatch). The costa and effort required to implement an advanced control algorithm, such as MPC, include the development of models to describe the process dynamics, the dedication of one or more computers, and the tuning of more parameters relative to analogue controllers. Thus, the implementation of a MPC is normally justified only for processes that cannot be adequately controlled by using simpler algorithms. These processes typically include the production of chemicals in high purities, chemical reactors with multiple steady-state and periodic attractors, extraction processes with narrow two- and three-phase regions (e.g., supercritical extraction), distillation towers with temperature and concentration fronts sensitively coupled to the reflux ratio (e.g., azeotropic distillation towers), and processes required to operate near many design and operating constraints. In the sections that follow, this paper briefly reviews the status of model predictive control, introduces a new algorithm for the control of nonlinear systems, and illustrates its performance for the control of a bioreactor and a catalytic CSTR.
Model Predictive Control The first modern control structures applied to chemical processes related the desired output to the manipulated input (for single-input-single-output (SISO) systems), through a linear step-response model, and to the unmeasured disturbances. These included Model Algorithmic Control (MAC) by Richalet et al. (1978), Dynamic Matrix Control (DMC) by Cutler and Ramaker (1979), and the inferential control methodology of Brosilow and Tong (1978). The principal advantage of these structures lies in their use of the linear step-response model to predict the process output as a function of past and proposed values of the manipulated input. An objective function, involving the differences between the predicted process output and the desired set point over several future sampling intervals, was minimized, subject to equality constraints that linearly relate the output and the manipulated input.
krP
k
k+l
k+2
k+3
k+M-1
k+M
k+P-1
k+P
Sampling Instant
Figure 1. Output horizon P and input horizon M for model predictive control. k is the index of the current sampling instant. x and u are the output and input deviations from the initial steady state.
These results were unified and extended in the framework of Internal Model Control (IMC) by Garcia and Morari (1982). In the absence of modeling errors (no process/model mismatch), they showed how the IMC control structure results in feedforward servo-control and feedback regulator control. Furthermore, for an open-loop (uncontrolled), stable process, stability of the controller guarantees stability of the closed-loop system. Also presented were several theorems which give the designer useful qualitative methods of tuning the controller. Most chemical processes are, however, nonlinear. While these control algorithms may be sufficiently robust for mildly nonlinear processes (or in the proximity of a steady state about which a linearization applies), they can be expected to yield poor (or even unstable) performance for highly nonlinear processes. Consider the more general MPC formulation, where the multiple-input-multiple-output (MIMO) process is modeled with nonlinear ordinary differential equations (ODES; eq 2) and the time horizons are illustrated in Figure 1. The nonlinear model is used to predict the vector of n outputs, 3, at each of P future sampling instants as a function of the vector of n manipulated inputs, g, at M future sampling instants, where k is the index of the current sampling instant. Note that the optimization algorithm adjusts the M vectors of manipulated inputs to minimize the control objective, more specifically, the Euclidean norm of the weighted deviations of the process outputs from the desired set points and changes in the manipulated inputs. It is minimized subject to constraints that arise in the process model, as linear and nonlinear inequalities involving the process variables (often “design” constraints) and as upper and lower bounds on the manipulated inputs. Mathematically, the following NLP can be formulated:
subject to
1814 Ind. Eng. Chem. Res., Vol. 28, No. 12, 1989
where E is a vector of n weights on the state variables, C = diag [ c ] , d is a vector of m weights on the manipulated inputs, and I) = diag [ d]. rj and @j are weighting coefficients that vary at each sampling instant over the predictive horizon. g and are vectors of equality and inequality constraints on the process operation. During the next sampling interval, only 4 is implemented. However, the values of &+I, ...,gk+M-land s +..., ~s+p , are useful in projecting the performance of the controlled system and providing early warnings when the system is approaching process design and operating constraints. Note that the NLP can be further generalized for processes modeled by differential-algebraicequations (DAEs). For these models, the algebraic constraints can be added to eq 3 and the first term of the objective function can be augmented to include the “algebraic” variables. This IMC framework has been the basis for many control algorithms, including Industrial Model-Predictive Control by Garcia and Prett (1986) and an algorithm by Economou et al. (1986) using operator theory. Other extensions for nonlinear systems include a Newton-based method by Li and Biegler (1988), Universal Dynamic Matrix Control (UDMC) by Morshedi (1986), collocation-based methods by Renfro et al. (1987) and Patwardhan et al. (1988),and a nonlinear, inferential algorithm by Parrish and Brosilow (1988). There are, however, several deficiencies inherent in the approximations to the dynamic constraints (2) used by the above algorithms. Among them are the following: (i) The Industrial Model-Predictive Control algorithm (Garcia and Prett, 1986) uses a nonlinear model that is periodically linearized to generate the step-response model used to optimize the values of the manipulated input, but the frequency of the relinearizations is not clearly stated. Furthermore, it estimates unmeasured disturbances by the law of superposition, which does not apply for nonlinear systems. Because the algorithm does not adequately account for the nonlinear behavior, its performance is suboptimal. (ii) The nonlinear, predictive controller proposed by Economou et al. (1986) computes an optimal vector of manipulated inputs to be implemented over the next sampling interval to produce the desired process outputs. This is a single-step method ( P = M = 1). Economou et al. solve the control problem without process constraints (see eq 3 and 4) by implementing a Newton-type convergence method, while Li and Biegler use an Armijo line search for the problem with constraints. The performance of these algorithms can be expected to be poor since a linear, multistep, predictive controller always gives a sluggish and suboptimal process response when limited to a single vector of manipulated inputs ( M = 1). The restriction to a single predictive step ( P = 1) can cause the algorithm to experience difficulties in avoiding impending process constraints, the most attractive feature of a MPC. (iii) Collocation-based methods (Renfro et al., 1987; Patwardhan et al., 1988) require the selection of trial functions to convert a system of ODEs to a system of algebraic NLEs. Solution of these NLEs is nontrivial, as is the effect of their convergence tolerances on the overall optimization algorithm. The problem is compounded by the need to limit the number of iterations to reduce the computational burden per predictive step. Furthermore, for highly nonlinear processes, the selection of the trial functions significantly impacts the performance of the algorithm. (iv) Morshedi extends the concepts of (Q)DMC to nonlinear processes by applying the mean value theorem over
each sampling interval. The solution of the resulting linear differential equations provides interpolating functions for each interval. For MIMO systems, this requires the calculation of eigenvectors and eigenvalues. Additionally, the UDMC algorithm also updates its predictive model, when process/model mismatch and unmeasured disturbances are encountered, by the law of superposition.
A New Control Algorithm The focus of our research is to develop a general algorithm that can efficiently control nonlinear processes near and within operating regimes characterized by hysteresis and periodic or even chaotic operation. To meet this goal, a new algorithm has been developed which constructs an easily calculable, explicit expression that relates the process outputs to the manipulated inputs. The NLP above is simplified by replacing the ODEs (2) with ..., %+p as a function of the maan expression for nipulated inputs at past and future sampling instants. With the dynamic constraints replaced, the NLP can easily be solved, using programs such as OPT (Lang and Biegler, 1987) and MINOS (Murtagh and Saunders, 1983), for the optimal values of the manipulated inputs and the associated process trajectory, given a set of weighting coefficients and tuning parameters. This result assumes that there is no process/model mismatch and that all of the disturbances are measurable. (For a nonlinear process, because the principle of superposition does not hold, an unmeasured disturbance is analogous to process/model mismatch.) In practice, however, the nonlinear model only approximates the process performance, sensor noise exists, and unmeasured disturbances normally enter the process. To solve the nonlinear control problem, the control algorithm needs to compensate for the process/model mismatch and sensor noise. Since the model is expected to deviate from the process, it is not important to obtain very accurate values of Q + ~ ,..., s +Instead, ~ it should be more efficient to obtain approximate solutions that display the proper dynamic trends while requiring few computations. To achieve this, the following approximation method has been developed. The reader should note how this method differs from those reviewed above (construction of the model inverse and the use of collocation) and numerical integration, which are implemented so as to achieve uniform accuracy across the sampling interval and, consequently, are computationally expensive. Let and a-lbe the deviation variables from the steady state at the current sampling ingtant, just prior to a step change in the inputs, and expand the right-hand side of eq 2 about % and L I - ~ in a first-order Taylor series:
Defining another set of deviation variables gives
Equation 7 can be expressed as
-A X = fJO,Ol + Jzg - + J,U_ Taking the Laplace Transform and simplifying gives
(8)
Ind. Eng. Chem. Res., Vol. 28, No. 12, 1989 1815 Then, for a E-step change in the input vector at the current sampling instant, 1
-u(sJ= -u S-
(10)
and
By use of a general expansion, the inverse can be expressed analytically (see the derivation in the Appendix), and
Sampllng Instant
k+J-l
k+J
Subinstant
2
1
1seg-l
Figure 2. Relinearization instants within a sampling interval.
where is a vector of n eigenvalues of 4-’&.Substituting, eq 12 becomes
whereas for the second and subsequent subintervals eq 17 becomes (L(k+])2
=
kk+j)l
+
( z k + ] )3
=
(xk+, 12
+ &k+]
(xk+J ) f + l
Equation 13 expresses the linearized dynamic constraints in the continuous Laplace domain. To obtain the output vectors at future sampling instants as a function of the past and future manipulated inputs, eq 13 is inverted. The resulting expression in the time domain, evaluated at t = T, is
Replacing the deviation variables gives
Equation 15 can be generalized to apply for projections of the output vector j sampling instants into the future: n
zk+j
= zk+j-l +
[B ~ ~ i e x i T l [ ~ u ~ ~-k !k+j-2) + j - l - + i-1-
+ fk+j-ll j = 1, ..., P (16)
which can be written as zk+j
= %+j-l + i k + j - l ( ! k + j - l
- Xk+j-2)
+*k+j-l
1 = 1, . . . I
p (17)
where
(19)
Equation 15 can be derived using the state transition matrix. However, this approach requires calculation of the eigenvectors of b-’Jz,a step not required in this algorithm. Normally, the sampling time is on the order of the smallest time constant of the process model. However, variations in the smallest time constant may require the linearizations to be performed more often. An alternative to reducing the sampling time, to obtain more accuracy and better control, is to implement the local linearization at equally spaced times within the sampling intervals, as illustrated schematically in Figure 2. Note that for the first subinterval eq 17 can be written as @k+j ) 1
= Zk+j-l + 9-k + j - l ( G k + j - 1 - !k+j-2)
+ $“k+j-1
(20)
?k+j
=
=
(Zk+J )I’
@k+])l
+
)Z
(ik+j
@k+j ) ~ s e g - l + &k+l
11’
)~eag-l
(21)
Over these subintervals, the manipulated inputs do not i’ = 1, ..., iseg-1, are calvary; &+J-l, &+J-l, and culated with a subinterval time, T/iseg. It is noted that the computational requirements increase linearly with the number of linearizations, and hence, the number of linearization instants per sampling interval (iseg) is chosen to balance the competing objectives of efficiency and accuracy. The use of repeated linearizationsof the nonlinear ODES differs significantly from the aforementioned methods of controlling nonlinear processes in which the model is linearized once and the linear model is used over the entire predictive horizon. (The latter is analogous to the linearized constraints for a single iteration of the SQP algorithm.) For highly nonlinear processes, this assumption can be expected to result in poor, or even unstable, control. Even when the aforementioned control algorithms perform satisfactorily, better responses are expected with the control algorithm introduced in this work. The problem of modifying the model to more accurately represent the process dynamics is difficult to address theoretically. When no measurement noise exists and the dynamic equations have the correct terms, the parameters can be updated to eliminate the process/model mismatch. With measurement noise or an incorrect structure for the model, it may be prudent to restrict the number of parameters that are updated, the rate at which they can change, or some combination of the two, to obtain a more effective updating procedure. Such strategies are analogous to those employed for parameter estimation and adaptive control. They have performed satisfactorily in industrial applications and have been applied successfully to update the models in the examples that follow.
Performance of the Nonlinear Control Algorithm The new control algorithm was tested by using both fermentation and catalytic reaction systems in CSTRs. In this section, parametric results are presented, in the absence of process/model mismatch, for various weights on the state variables ( c ) , predictive horizons ( P , M ) , and relinearization instants per sampling interval (iseg). Then, the results are revised to show the effects of process/model mismatch due to errors in the model parameters. Finally,
1816 Ind. Eng. Chem. Res., Vol. 28, No. 12, 1989 lz
3.s
-
3.0
-
2.5
-
2.0
-
1.5
-
0- 1.0
-
0.5
-
/
X, s
Figure 3. Schematic diagram of the bioreactor.
the results are presented using algorithms to compensate for the occurrences of mismatch. Example 1. Bioreactor Process. The growth rate of most cells in fermentation processes is inhibited by large substrate concentrations. The effect of this inhibition on steady-state and periodic solutions was studied by Agrawal et al. (1982), who utilized the one-hump growth model:
-9v
. . . ,. .. .. . . . .. . .. .... .. .. . , , . . ... .. .. .. .
0.0 -
p{S/ =
kSe-S/K (22) where 1.1 is the specific growth rate (maximum at S = K ) and S is the substrate concentration. The specific substrate consumption rate, u, was related to the specific growth rate through the yield coefficient, Y , in the form
For the CSTR in Figure 3, the cell and substrate mass balances are F = --x I.1(S]X dt V
+
g
=
$s* - S ) - u ( S ) X
(25)
-0.5
0.0
I
I
I
I
I
I
I
I
I
1
0.2
0.4
0.6
0.8
7.0
1.2
1.4
1.8
7.8
2.0
of periodic oscillations. A designer might select Da = 1.30 and examine the controllability of the process in response to changes in the set point (CTp) and unmeasured disturbances. A t Da = 1.30, the reactor is marginally stable with eigenvalues of -0.077 965 f 2.28741'. To explore the controllability about such a steady state, it is convenient to introduce deviation variables, X l ( t )= Cl(t)- Cls, X 2 ( t ) = Cz(t)- C2s, and u ( t )= F / F s - 1. Then, eq 26 and 27 become
Dedimensionalizing, these mass balances become
where C1 and C2 are the dimensionless cell mass [X/ (sFy(sF))] and substrate conversion [(SF - s)/sF], p= a/(bSF),y = K/SF, r = t ( F / V ) ,and Da is the ratio of the specific growth rate to the residence time [g{SF)/(F/ V)]. To identify the regimes of hysteresis and periodic operation, Agrawal et al. varied the parameters Da, 8, and y. Figure 4,in which the nonlinear steady-state and periodic branches of the solution diagram with y = 0.48 and 0 = 0.02 are shown, is illustrative of the results they obtained. Note that the stable steady state that yields the highest cell mass occurs at the Hopf bifurcation point, where Da = 1.206 25. Beyond that, the steady-state attractors are unstable and the periodic attractors are stable with large periods, corresponding to the experimental observations of Borzani et al. (1977) for the aerobic fermentation of Saccharomyces cerevisiae on a sugar cane-molasses medium. It is noteworthy that this semiempirical model does not oscillate when the yield coefficient is constant (b = 0, p = m). Agrawal et al. (1982) expressed concern about the range of substrate concentrations over which eq 2 applies. If justified, for design and control optimizations, a more complete model that accounts for the rate of oxygen consumption might be prepared. Such a model, to our knowledge, has not yet appeared in the literature. In the design stage, it seems desirable to operate at the highest yields of cell mass but not too close to the region
1
Das(l - X 2 - CzS)e(Cs+X2)/r 1+ p -
I
R
L ' t J
xz - c,s (ClS + Xl) (29)
where I S is t ( F s / V ) ,Das = p(SFJ/(FS/V), and u is the dimensionless manipulated flow rate, expressed as a deviation from the steady state (u = 0). In one simulation, it was decided to explore the use of MPC to respond to a set point change that increases the dimensionless cell mass by 0.05. Figure 5a shows that, for the uncontrolled process, this corresponds to a periodic attractor and two unstable steady-state attractors at u = 0.38456, with eigenvalues 0.21226 f 2.0835i, and u = 0.606 16, with eigenvalues 0.367 79 and -1.2426. For each of the unstable steady-state attractors, the substrate conversion is shown in Figure 5b. To implement MPC, a sampling time of T = 0.5 time unit was selected, based upon typical transient responses in the vicinity of the initial steady state. Figure 6 shows the responses with P = M = 4,iseg = 3, and various values of E. Here, the response with c2 = 0.25 was judged to be the best. Figure 7 shows the responses with c = (1.0,0.25), iseg = 3 and P ( = M ) = 1, ..., 4. Note that the cases with P = 1and 2 are suboptimal and the results with P = 3 and 4 show no improvement when larger predictive horizons
Ind. Eng. Chem. Res., Vol. 28, No. 12,1989 1817 6.0
2.5
2.0
5.0
1.5 4.0
1 .o
-5
3.0
h
h
0.5
9
v
2.0
v
x-
0.0
X1
.o
-0.5
0.0
-1.0
-1 3
I 3 -0.8
I
1
I
I
I
-0.6
-0.4
-0.2
0.0
0.2
I
I
0.4
0.8
I 0.8
1
-1.0
1.0
I
I
I
I
I
I
I
0.5
1.0
1.5
2.0
2.5
3.0
U
o.2
I 3.5
I
I
I
4.0
4.5
5.0
7s
1
-0.2
-0.4
2 -0.6
...' ..'
...'
,,.."
-0.8
-1'8 -2.0
-1.0 -1.5
-1.0
-0.5
0.0
0.5
1.0
1.5
2.0
2.5
x, (x lo+') Figure 5. Steady-state ((-1 stable, (.-.) unstable) and periodic (0) solutions for the Agrawal bioreactor with Das = 1.30, y = 0.48, and 0 = 0.02. u = 0.0 corresponds to steady-state operation. (a, top) Cell mass as a function of the rate of flow. (b,bottom) Substrate conversion as a function of cell mass.
(P 1 3) are used. Furthermore, other results show that iseg 2 2 is necessary to obtain good performance. With P = 3 and iseg = 2,35.7 CPU seconds were required using a Micro VAX I1 computer to traverse 10 sampling intervals. Note that this computation time includes approximately 4 CPU seconds for simulation of the process using LSODE and about 2 CPU seconds for initialization; MINOS was used to solve the NLP for all of the cases in this manuscript. Although Chang and Chen (1984) demonstrated that simple PID controllers can be tuned to stabilize these steady states, their responses were much more sluggish and oscillatory and exhibited a greater sensitivity to the tuning parameters than the MPC algorithm. Figure 8 illustrates the impact of errors (12%)in the model parameter y when the MPC is not augmented with a mismatch compensation algorithm. Here, ymdelis fixed at 0.48 and the process is assumed to operate with an
1' -1.0
I 0.0
I 1.0
I 2.0
I 3.0
I 4.0
I 5.0
I 6.0
x, (x in+*) Figure 6. Transient response to step change in ClSp with T = 0.5, P = M = 4, iseg = 3, and cl = 1. (a, top) Cell mass. (b, bottom) Substrate conversion-cell mass phase plane.
effective yprocess ranging from 0.47 to 0.49. Note the s i g nificant offset in the resulting steady state when process/model mismatch occurs. Similar results are obtained for errors on the order of 20% in the model parameter 0. Unfortunately, the offset cannot be eliminated because the process is undercontrolled; that is, the two outputs (XI and X , ) are adjusted by only one input (u). When ypf ymodel, with ymdeland Das estimated to eliminate the process/model mismatch, the offset is eliminated, as illustrated in Figure 9. This is because both eq 28 and 29 in the process model at steady state can be satisfied by the same u when XI = XIsp and X z = XzSp,as and Das vary, with Pmodelfixed. When fiprocess # Bmodelhphowever, different values of u are found when XI = X I and X , = X,Sp are substituted in the two equations, and consequently, the offset cannot be eliminated. To eliminate the offset, a second feed stream, of different concentration, can be added to provide a second manipulated input. Example 2. Catalytic CSTR. In this example, the controllability of a catalytic CSTR that converts A + B
1818 Ind. Eng. Chem. Res., Vol. 28, No. 12, 1989
0.5
0.0
1.5
2.0
2.5
3.0
3.5
4.0
4.5
5.0
0.0
I
I
I
I
I
I
I
I
I
1
0.5
1.0
1.5
2.0
2.5
3.0
3.5
4.0
4.5
5.0
-
-0.2
-
-0.4
-
-0.6
-
-0.8
-
-1.0
-
-1.2
-
-i,
1.0
h
.-(
v
X-
-1.4
I I j
-1.0
0.0
1.0
2.0
3.0
4.0
5.0
6.0
-1.4
-7.0
x, (x
10'2)
Figure 7. Transient response to step change in Clsp with T = 0.5, iseg = 3, and c = (1.0,0.25). (a, top) Cell mass. (b, bottom) Substrate conversion-cell mass phase plane.
-
P is examined. A is in excess, and the rate of consumption of B is
The reaction is exothermic, but an excessive rate of cooling is assumed to allow isothermal operation of the reactor. For this design, Figure 10 illustrates a configuration proposed to control the product concentration, c b , and the liquid level, h, by adjusting the flow rate, wl,of the concentrated feed ( c b l ) and the flow rate, wz,of the dilute feed (cbz). This configuration, which was presented by Li and Biegler (1988), has the dynamic model -dh=
dt
w1 + w 2 - 0.2h0.6
(31)
0.0
1.0
2.0
3.0
4.0
5.0
6.0
7.0
x, (x 10'2) Figure 8. Transient response to step change in Clsp with and = 0.48, T = 0.5, P = M = without process/model mismatch. -ymdel 4, iseg = 3, and = (1.0,0.25). (a, top) Cell mass. (b, bottom) Substrate conversion-cell mass phase plane.
Figure 11illustrates the three-dimensional, steady-state surface of C b as a function of the inlet flow rates with k1 = k2 = 1, c b l = 24.9, and Cbz = 0.10. As shown, three steady-state attracton (two stable and one unstable) exist when w1 = w2 = 1. Li and Biegler chose to demonstrate, using their single-step method, the controllability of the reactor in reaching the unstable steady state, hs = 100.0 and Cbs = 2.7927, from an initial condition of low volume and concentration (h = 40.0 and C b = 0.100). Note that at the set point the model is unstable with respect to the concentration only, since it has two eigenvalues, -0.0100 and 0.012 860. Figure 1 2 shows the response of the reactor, under the new control algorithm, to a simulated start-up from the initial condition. The tuning parameters were set at P ( = M ) = 1, ..., 4, c = (1,400), and keg = 3. The sampling time, T = 1.0, was chosen to facilitate a direct comparison with the results of Li and Biegler. In addition, bounds on
Ind. Eng. Chem. Res., Vol. 28,No. 12, 1989 1819 0.0
25
5.0
........................
.........
...............................
4.0
3.0
0 h
0
93
2.0
0
X1 .o
Figure 1.
Cb
response surface.
1 .C
0.0
-1.0
I
I
I
I
I
I
I
1
1
1
0.5
1.0
1.5
2.0
2.5
3.0
3.5
4.0
4.5
5.0
0.0
-0.2
-0.4
-0.0 h
i
;-0.0
-3.0
v
I 0.0
x"
1
1
1
1
1
1
1
1
1.0
2.0
3.0
4.0
5.0
0.0
7.0
8.0
1
t
-1.0 1 .o -1
.a 0.5
-1.4
I I
0.0
I 1.0
x,
(x
T
1
I 2.0
3.0
4.0
1 >
0.0
0.0
10'2)
Figure 9. Transient response to step change in ClSp with pro= 0.47) showing the effect "/model mismatch ( y d = 0.48, yof mismatch compensation. (a, top) Cell mass. (b, bottom) Substrate conversion-cell mass phase plane.
-0.5
2 -
u
-1.0
U" -1.5
-2.0
-2.5
-+ W
Figure 10. Schematic diagram of the catalytic CSTR.
the inlet flow rates of 0.0 and 10.0 were implemented. Note that the new control algorithm responds optimally for P 1 2. A comparison of these results with those reported by Li and Biegler, and reproduced in Figure 13 (curve B), shows the marked improvement of the new algorithm in reducing both the overshoot and settling time. It should be noted that the deviations of the process outputs from
-3.0
.
1
1
1
1
1
1
-00.
-50.
-40.
-30.
-20.
-10.
h
1 0.
10,
- hsp
Figure 12. Transient response to reactor start-up with T = 1.0, iseg = 3, and E = (1,400). Inlet flow rates are bounded between 0 and 10. (a, top) Concentration. (b, bottom) Concentration-liquid level phase plane.
their steady-state values are plotted in Figure 12,while the process outputs alone are plotted in Figure 13. Li and Biegler required approximately 30 CPU seconds on a Micro VAX I1 computer to control and simulate the re-
1820 Ind. Eng. Chem. Res., Vol. 28, No. 12, 1989 9
2.0
8
1 .5
7
t
6
.o
0.5
5
0.0
u" '
-0.5
G
3
-1.0
2 -1.5 1
-2.0
0 20
-2.5
LO
' t
-3.0
I
I
I
I
1
I
I
I
1
1.0
2.0
3.0
4.0
5.0
6.0
7.0
8.0
t
IZ0 110
1
Figure 14. Transient response to reactor start-up with and without process/model mismatch. Cbl,mdel = 24.9, T = 1.0, P = M = 4, iseg = 3, and g = (1,400).
100
I I / d
I
50
1/ ::::/
'f -0.5
40
I
0
I
20
I
i LO
t Figure 13. Transient response to reactor start-up (Li and Biegler, 1988). The upper bounds on the inlet flow rates are 5 (A), 10 (B), and (C). The lower bound is 0. (a, top) Concentration. (b, bottom) Liquid level.
-
U '
-1.0
-2.5
-3.0
actor over 50 sampling intervals. Further examination of the cases with P = 1 and 2 and iseg = 1, 2, and 3, show far better performance with CPU times that vary from 7.5 to 45 CPU seconds over 8 sampling intervals. Figure 14 illustrates the impact of 10% and 20% errors in cbl when MPC is not augmented with a mismatch is fixed at 24.9 and compensation algorithm. Here, Cblsnodel the concentration of inlet flow 1, Cbl, -, varies from 19.9 to 29.9. Note that the offset is smafi but that the effects on the transient are significant. Figure 15 illustrates the improved performance of the control algorithm when it is augmented by a parameter estimation algorithm. Finally, Figure 16 shows the ability of the control algorithm to reject a 20% disturbance in Cbl. Note that the mismatch compensation algorithm completely eliminates the offset.
Conclusions The following are concluded. (1)The multistep MPC algorithm derived in this paper provides excellent servo-control and disturbance rejection in attraction to unstable steady states, as demonstrated
4/
I
0.0
I
I
I
I
I
I
I
I
1.0
2.0
3.0
4.0
5.0
6.0
7.0
8.0
t
Figure 15. Transient response to reactor start-up with process/ = 27.4) showing the effect model mismatch (Cbl@&, = 24.9, Cb1,pof mismatch compensation.
for the fermentation and catalytic reactors. (2) The multistep algorithm performs more effectively than its single-step counterpart and is easily tuned to give excellent performance. (3) As reviewed in this paper, many of the existing methods use a single predictive step ( M = P = 1). These methods are comparable to deadbeat control for linear systems and should be far less effective and efficient than the new algorithm. This is demonstrated in a comparison with the results of Li and Biegler (1988) for the catalytic reactor. (4)Although the multistep method is less accurate in approximating the dynamic model across an entire sampling interval, it requires no iterations and appears to be sufficiently accurate. The results for the fermentation and catalytic reactors show excellent performance and high efficiencies.
Ind. Ecg. Chem. Res., Vol. 28, No. 12, 1989 1821 e’o
Acknowledgment
1
Partial funding was provided by the Design Theory and Methodology Program of the NSF under Grant DMC8613484 and is gratefully acknowledged.
5.0
4.0
5.0
Nomenclature b = see eq 23 = n x n coefficient matrices c = vector of n weights on the state variables C = diag [SI C1 = dimensionless cell mass, X/(SFY(SF)) Cz = dimensionless substrate conversion, (SF- S)/SF C b = concentration of species B d = vector of m weights on the manipulated inputs Q = diag [dl Da = Damkohler number, P ( S F J / ( F / V ) Das = Damkohler number at steady state, ~{SF)/(P~/V) f = see eq 2 k = feed flow rate = vector of equality constraints on process operation $; = liquid level h = vector of inequality constraints on process operation iseg = number of linearizations per sampling interval J = Jacobian matrix k = growth model coefficient k = vector of rate constants in eq 30 and 32 K = maximum of one-hump growth model M = manipulated input horizon Pi = coefficient matrices in eq A l , i = 0, ..., n - 1 P = predicted output horizon r = intrinsic rate of reaction s = parameter in Laplace transform S = substrate concentration t = time T = sampling interval u = manipulated input; F / ( F s - 1) u = deviation variable, u-lk-l V = volume of fermenter w = vector of inlet flow rates x = output = deviation variable, -xk X = cell mass X1 = dimensionless cell mass deviation from steady state, C1 a,
4, €
-1.0
-2.0
-3.0
-4.0 0.0
I
I
1
I
I
I
I
I
1.0
2.0
3.0
4.0
5.0
6.0
7.0
0.0
t
0.5
7
0.4
0.3
0.2
3
B
6-
8
6B
0.1
0.0
-0.1
u
-0.2
- O -0.4 -0.5
I I
I
0.0
1.0
I 2.0.
I
I
I
I
I
1
3.0
4.0
5.0
6.0
7.0
0.0
t
- CIS
X z = dimensionless substrate conversion deviation from steady state, Cz - C,S Y = yield coefficient Greek Symbols sum of all products of n eigenvalues taken i at a time, eq A5 @ = weighting coefficient, a/(bSF) y = weighting coefficient, K/SF A = vector of n eigenvalues 4 = coefficient matrix p = specific growth rate u = specific substrate consumption rate T = t(F/V) @ = see eq 18 $ = see eq 19 Subscripts b = chemical species B F = feed k = index of current sampling instant S = steady state 0 = initial value ai =
-4.0
I -1.0
I -0.5
1 0.0
I
I
0.5
I 1.5
1.0
-
h h”
(X
1 2.0
I 2.5
f 3.0
10”)
Figure 16. Rejection of a 20% disturbance in Cbl with T = 1.0, P = M = 4, iseg = 3, and c = (1,400). Inlet flow rates are bounded between 0 and 10. (a, top) Concentration. (b, middle) Manipulated inputs with mismatch compensation. (c, bottom) Concentrationliquid level phase plane.
Superscripts H = upper bound
1822 Ind. Eng. Chem. Res., Vol. 28, No. 12, 1989
L = lower bound SP = set point i=l
Appendix To invert SA - J x , postulate that the inverse exists and can be expressed in the expansion Ifl
- Ec
Gl
Note that eq A&A10 apply when Xi, i = 1, ...,n, are unique. Modifications are necessary when the eigenvalues are not unique. Literature Cited
det
[sl- 4-lJx] - - = fi(s - hi)
(A3)
i=1
where is a vector of n eigenvalues of A-'Jx. Furthermore, det
[sl_ _4-'Jr] - = s" + Pial + ... + sa,-1 + a ,
(A41
e-'&.eigen-
where aiis the sum of all of the products of values taken i at a time:
Then combining (A2) and (A4) and rearranging gives
i= -PI - = ( a -l l + cl-lJx)iy - - -
1, ..., n - 1
(A71
Then, to obtain s-l(s4 - JX)-', first the eigenvalues of are computed (eq A3), g is determined (eq A5), and Po, ..., Pn-'are computed. Combining eq A1 and A3 gives
where
Agrawal, P.; Lee, C.; Lim, H. C.; Ramkrishna, D. Theoretical Investigations of Dynamic Behavior of Isothermal Continuous Stirred Tank Biological Reactors. Chem. Eng. Sci. 1982, 37(3), 453. Borzani, W.; Gregori, R. E.; Vairo, M. L. R. Some Observations on Oscillatory Changes in the Growth Rate of Saccharomyces cerevisiae in Aerobic Continuous Undisturbed Culture. Biotechnob Bioeng. 1977, 19, 1363. Brosilow, C. B.; Tong, M. Inferential Control of Processes: Part 11. The Structure and Dynamics of Inferential Control Systems. AIChE J. 1978, 24(3), 492. Chang, H.-C.; Chen, L.-H. Bifurcation Characteristics of Nonlinear Systems Under Conventional PID Control. Chem. Eng. Sci. 1984, 39(7/8), 1127. Cutler. C. R.: Ramaker. B. L. Dynamic Matrix Control-A Computer Control Algorithm, Presented at the 86th AIChE National Meeting, Houston, TX; American Institute of Chemical Engineering: New York, 1979. Economou, C. G.;Morari, M.; Palsson, B. 0. Internal Model Control. 5. Extension to Nonlinear Systems. Ind. Eng. Chem. Process Des. Deu. 1986,25,403. Garcia, C. E.; Morari, M. Internal Model Control. 1. A Unifying Review and Some New Results. Ind. Eng. Chem. Process Des. Deu. 1982, 21, 308. Garcia, C. E.; Prett, D. M. Advances in Industrial Model-Predictive Control. Third International Conference on Chemical Process Control; Morari, M., McAvoy, T. J., Eds.; CACHE, Elsevier: New York, 1986. Lang, Y.-D.; Biegler, L. T. A Unified Algorithm for Flowsheet Optimization. Comp. Chem. Eng. 1987, 11(2), 143. Li, W. C.; Biegler, L. T. Process Control Strategies for Constrained Nonlinear Systems. Ind. Eng. Chem. Res. 1988, 27, 1421. Morshedi, A. M. Universal Dynamic Matrix Control. Third International Conference on Chemical Process Control;Morari, M., McAvoy, T. J., Eds.; CACHE, Elsevier: New York, 1986. Murtagh, B. A.; Saunders, M. A. MINOS 5.0 User's Guide;Systems Opt. Lab, Stanford University: Palo Alto, CA, 1983. Parrish, J. R.; Brosilow, C. B. Nonlinear Inferential Control. AIChE J. 1988, 34, 633. Patwardhan, A. A.; Rawlings, J. B.; Edgar, T. F. Model Predictive Control of Nonlinear Processes in the Presence of Constraints, Presented a t the 80th AIChE Annual Meeting, Washington, DC; American Institute of Chemical Engineering: New York, 1988. Renfro, J. G.; Morshedi, A. M.; Asbjornsen, 0. A. Simultaneous Optimization and Solution of Systems Described by Differentiallhlgebraic Equations. Comput. Chem. Eng. 1987, 11, 503. Richalet, J. A,; Rault, A.; Testud, J. D.; Papon, J. Model Predictive Heuristic Control: Applications to Industrial Processes. Automatica 1978, 14, 413. Received for review March 27, 1989 Accepted August 11, 1989
