484
Ind. Eng. Chem. Process Des. Dev. l985?24, 484-494
that the CL do not cross the real axis to the left of -l/(l - a ) with a* 5 a < 1. Literature Cited Alevlsakis, 0.; Seborg, D. E. Int. J. Control 1973. 3 , 541. Bristol, E. H. “A New Process Interaction Concept: Pinned Zeros”; Foxboro Company: Boston, MA, 1980. Brosllow, C. B. Paper presented at the Proceedings of the Joint Automatic Control Conference. Denver, CO, 1979. Calller, F. M.;C. A. Desoer “Muithrariable Feedback Systems”; Springer-Verlag: New York, 1982. Cutler, C. R.; Ramaker, B. L. Paper presented at the AIChE 86th National Meeting, Houston, April 1979, and the Proceedings of the Joint Automatic Control Conference, San Francisco, 1980. Desoer, C. A.; Chen, M. J. IEEE Trans. A m m . Conhol1981, AC-26, 408. Desoer, C. A.; Gustafson, C. L. Memorandum No. UCBlERL M82/80, l982a, Electronics Research Laboratory, Unhwsity of Callfornia, Berkeley, CA. Desoer, C. A.; Gustafson, C. L. Memorandum No. UCBlERL M82-90, 1982b, Electronics Research Laboratory, University of California, Berkeley, CA. Doyle, J. C.; Stein, G. IEEE Trens. Autom. Control 1981, AC-26, 4. Economou, C. G.; Morarl, M. AIChE J., submitted for publtcatlon. Foss, A. S. AIChEJ. 1973, 19, 209. Foss, A. S.;Edmunds, J. M.; Kowarltakls, B. Ind. Eng. Chem. Fundam., 1980, 19, 109. Frank, P. M. “Entwwf von Regelkreisen mlt vorgeschrlebenem Verhalten”; G. Braun Veriag: Karlsruhe, 1974. Garcia, C. E. Ph.D. Thesis, Unlv. of Wisconsin-Madison, 1982. Garcia, C. E.; Morari, M. Ind. Eng. Chem. Process Des. Dev. 1982, 21, 308. Holt, B. R.; Morari, M. Chem Eng . Sci., in press. Jafarey, A.; McAvoy, T. J. Ind. Eng. Chem. Process Des. Dev. 1978, 17, 485. Kestenbaum, A,; Shinnar, R.; Thau. F. E. Ind. Eng. Chem. Process Des. D e v . 1978, 15, 2.
Kwakemaak, H.; Slvan, A. “Linear Optimal Control Systems”; Wilay-Interscience: New York, 1972. Lee, W.; Weekman, V. W., Jr. AIChE J. 1976, 22, 27. MacFarlane, A. G. J.; Kowarltakis. B. Int. J . Control 1977, 25, 837. Mehra, R. K.; Rouhani, R.; Eterno, J.; Rlchalet, J.; Rautt, A. I n “Chemical Process Control 2”; Edgar, T. F., Seborg, D. E., Eds.; United Engineering Trustees Inc.: New York, 1982, p 287. Morarl, M. I n ”Chemical Process Control 2”; Edgar, T. F., Seborg, 0.E., Eds.; United Engineering Trustees, Inc.: New York, 1982, p 467. Morarl, M. Chem. Eng. Sci. 1983, 38, 1881. Mora!i, M.; Grlmm, W.; Oglesby, M. J.; Prosser, I. D. Chem . Eng . Sci.. submMed for publication. Newton, G. C.; GouM, L. A.; Kaiser, J. F. “Analytical Design of Linear Feedback Controls”; Wlley: New York, 1957. Ogunnake, 8. A.; Ray, W. H. A I C M J. 1979. 25, 1043. Ogunnake, B. A. Ph.D. Thesis, Unlv. of Wisconsin-Madison, 1981. Palmor, 2. J.; Shinnar, R. Ind. Eng. Chem. Process Des. Dev. 1979, 18, 8. Palmor, 2. J.; Powers, D. V. Paper presented at the Proceedings of the Joint Automatlc Control Conference, Charlottesvllle, VA, 1981. Postiethwaite, I.; MacFarlane, A. 0.J. “A Complex Variable Approach to the Analysis of Linear Multivariable Feedback Systems”; Springer-Verlag: New York, 1979. Rlchalet, J. A.; Rautt, A.; Testud, J. D.; Papon, J. Automatica 1978, 14, 413. Rosenbrock, H. H. “Computer-Aided Control System Design”; Academic Press: New York, 1974. Smith, 0. J. M. Chem. Eng. Progr., 1957, No. 5 , 217. Soliman, M. A.; Ray, W. H. Int. J. Control 1972, 15, 609. Wood, R. K.; Berry, M. W. Chem. Eng. Sci. 1973, 28, 1707. Zames, G. IEEE Trans. Autom. Control 1981, AC-26.
.
Received for review April 25, 1983 Revised manuscript received May 23, 1984 Accepted July 23, 1984
Internal Model Control. 3. Multivariable Control Law Computation and Tuning Guidelines Carlos E. Garcia’ and Manfred Morarl’ Chemical Engineering Department, University of Wisconsln, Madison, Wisconsin 53706
A model predlctive control law formulation is suggested for the computation and implementation of the Internal Model Controller. The advantages of this approach are that the controller parameters can be calculated explicitly without inversion or factorizatbn of polynomial matrices and that the IMC perfect controller is obtained as a special case. Furthermore, this formulation can be extended in a straightforward manner to take general linear process constraints into account when the control Inputs are computed. I n W e I y plausible tuning parameters are introduced and theorems are derived to elucidate their effect on Stability and performance. Simulation examples show the available flexibility in imposing deslrable dynamics on the response.
1. Introduction In part 2 of this series, a new multivariable control technique based on the Internal Model Control (IMC) structure was introduced. The fundamental stability and performance relationships were derived, and basic design guidelines were given. As a first step of the design procedure, the controller is assumed ta be the realizable and stable part of the system inverse. This controller yields ”optimal” system response characteristics. In the second step, a filter is added to provide robustness, i.e., stability
* Address all correspondence to this author a t the Chemical Engineering Department, 206-41, California Institute of Technology, Pasadena, CA 91125. Current affiliation: Shell Development Co., Houston, TX 77001. 0196-4305/85/1124-0484$01.50/0
and acceptable performance deterioration in the presence of model/ plant mismatch. A series of simulation examples in part 2 demonstrated that, in general, the results obtained by this technique are very good and that the structure makes on-line tuning very simple: However, it has two drawbacks. For one, to obtain the realizable and stable part of the system inverse by factoring, a matrix whose elements are rational functions of polynomials is awkward and full of numerical difficulties. Furthermore even when this step is successful, zeros close to the unit circle (UC) can give rise to undesirable response characteristics (“rippling”). To alleviate these two potential problems, we discuss here-similarly as for SISO systems in part 1-a model predictive control law formulation: The system inputs are calculated in an open-loop fashion such that the system response predicted on the basis of a process model has desirable character0 1985 Amerlcan Chemical Society
Ind. Eng. Chem. Process Des. Dev., Vol. 24, No. 2 , 1985 485
istics. This procedure yields in general an approximate inverse and as a special case the exact inverse of the invertible part of the process model. We will discuss how the choice of tuning parameters affects the approximation and through it the form and speed of response.
I
FILTER
CONTROLLER
d
PLANT
+ G
2. Model Formulation As in part 2, we employ the discrete-time linear formulation r
x(k
r
+ 1) = k ( k ) + i=l C j=1 C Tijm(k-
7ij)
+ w(k)
which describes the behavior of a linear continuous system at the discrete times kT, k = 0, 1, ...,where T i s the sampling interval, x E R” is the system state vector, y E R‘ is the system outputs, m E Rr is the system inputs, and w E R” is a vector of disturbances. Taking z transforms, model 1 becomes y ( z ) = C(z1 - !B)-l(CC Tijz-Wn(z) i=l j=1
+ w ( z ) )= G(z)m(z) + 4%) (2)
where d(z) is the new disturbance vector and G(z) is the system transfer matrix. The system is assumed to be open-loop stable (i.e., the roots of det (XI - a) = 0 lie inside the unit circle). Thus, one can equivalently represent G in “causal sequence” or “convolution” form by expanding (zI in an infinite series
G(z) = CC@l-l 1=1
Tij~-Ti~-l =
i s 1 j=1
Figure 1. Internal model control structure.
In the following, a short review of the IMC design procedure is given. For more details, the reader is referred to the companion paper (part 2).
~ ( k=) Cx(k)
r
I N T E R N A L MODEL
(1)
x(0) = 0
r
d
HI + H g l + H g 2 + ... (3)
where 8 1 E Rrxrare the impulse response matrices. Note that G(z) can also be put in the form
where each hi,(%)is an impulse response polynomial sequence whose coefficients compose the matrices 8,;
clearly for stable systems lim hij, = 0 V i , j I--
and therefore, each element of G(z) in (4) can be approximated by a finite number of terms Nijto any desired degree of accuracy. Here Nijdenotes the truncation order of hij(z) and N = maxij Nip Nonparametric models, e.g., impulse response models, have advantages and disadvantages. Our main reason for using them here is that the control law computation becomes extremely simple. The “optimal” control law can be calculated explicitly while the standard parametric-state space models lead to nonlinear Riccati equations. The storage requirements for nonparametric models are larger, but with rapidly increasing chip capacity, this should not be of serious concern. Nonparametric models are not parsimonious and should not be used for process identification and parameter estimation unless the data are essentially noise-free.
3. IMC in a Nutshell IMC has the feedback structure shown in Figure 1-where 4%)= set point for output y(z), G(z) = plant, G(z) = internal model, G,(z) = IMC controller, and F(z) = filter. The closed-loop expression relating inputs and outputs is Y(Z)
= G(z)[I + G,(z)F(z)(G(z) - G(z))I-’G,(z)F(z) X ( s ( z ) - d k ) ) + 4%) (6)
which for the case G = G becomes ~ ( 2= ) G(z)G,(z)F(z)(s(z) - d(z)) + 4%) (7) Thus, in the absence of modeling errors and disturbances, IMC behaves like the feedforward controller G,(z)F(z). The following properties can be derived from (6) and (7). Property 1. Dual Stability. Assume G = G and G, G,, and F stable. Then IMC is closed-loop stabje. Property 2. Perfect Control. Assume G, = G-l; F = I; IMC is closed-loop stable. Then y ( z ) = sk), V d(z). Property 3. Zero Offset. Assume G,(l) = G W ; F(1) = I. IMC is closed-loop stable; the disturbance and the set point are asymptotically constant (lim(t a) s ( t ) = so). Then lim(t a) y ( t ) = sg. It is important to emphasize that property 1 holds only in the absence of modeling errors (G = G). If the model does not describe the-plant precisely, the system can be unstable even when G, G, G,, and F are stable. In part 2, we explained how the filter F can be used to stabilize a system in the presence of model/plant mismatch. In this paper we concentrate on the design of G, and the issue of dynamic performance. Therefore, an exact model is assumed throughout part 3. In terms of performance, G, = G-l would be most desirable because it yields perfect set point following. However, this controller can usually not be implemented because of the following reasons. (i) Inversion of time delays in G leads to prediction and thus to a controller which is cot causal. (ii) Transmission zeros of G outside the UC become unstable controller poles. (iii) If transmission zeros of G lie inside but close to the UC, severe rippling behavior of the inputs and, consequently, strong oscillations of the outputs between samples are likely to occur. In order to circumvent these problems, it is possible to factor the transfer matrix into two parts as was suggested in the companion paper and to use as a controller G, the inverse of the “well behaved” part only. The complexity of the matrix excludes an analytical factorization and inversion in all but trivial cases. Numerical techniques for this purpose are error-prone. Thus, a simple procedure is needed which can provide the exact inverse or a suitable approximation.
-
-
486
Ind. Eng. Chem. Process Des. Dev., Vol. 24, No. 2, 1985
4. Approximate Inverses (1) General Discussion of Model Predictive Methods. As was suggested in part 1 for SISO systems, it is possible to generate approximate inverses with almost any desired characteristics by solving a simple optimization problem. The process model can be employed to predict the outputs resulting from a series of inputs, or, alternatively, desired outputs can be prescribed and the inputs could be calculated such that the predicted outputs follow the prescribed outputs in some “optimal” manner. If one requires the predicted values to agree with the prescribed ones exactly, the system inputs resulting from the solution of this matching problem will be the same as would be obtained by an inversion of the process model. If one requires the predicted values only to be close to the desired ones in the least-squares sense, the solution of the optimization problem will provide an approximate inverse of the process model. The characteristics of the approximate inverse can be affected by the choice of weighting matrices in the least-squares objective function. The IMC design procedure proposed in part 2 suggests to factor the process model as G(z) = G+i(z)G+dz)G-(z) (8)
where G+,(z) is selected to make G-lG+l(z) realizable and G+zto make G-lG+z stable. We saw that when the controller is selected as G, = G--l, an optimal time delay compensator is introduced automatically into the feedback path. Our objective is to propose a formulation of the model predictive problem such that this desirable feature is preserved and that-at least in the absence of zeros outside the UC (G+z= 1)-the exact inverse controller G, = G--l can be obtained as a special case. This is straightforward for SISO systems but, as we will see, quite difficult for multivariable systems. Let us pose the following problem to be solved at time
significantly more storage than necessary, however. As we will explain below, it is more efficient to start the optimization period T intervals later (because the first T outputs cannot be affected anyway) and to eliminate the first 7 impulse response matrices which are equal to zero. Case 111. Different Time Delays in the Different Transfer Function Elements. The closed-loopresponse of IMC using the inverse controller is ~ ( 2= ) G+(z)(~(z) - d(z)) + d(z)
As discussed in part 2, G+(z)is a matrix of delays selected to make G,(z) = G-(z)-’ = G(z)-lG+(z) realizable. In general G+(z)is not unique. There is always a unique G+(z),however, which gives rise to a decoupled system response G+*(z) = diag (z-(‘lt+1)z-(‘zt+1)...%-(‘~++1)) (10) In case 11, the inputs had no effect on the outputs before the very time step when all of them reached the set point. In case 111, the inputs can generally cause changes in the outputs before these are able to settle. Let us illustrate this with a simple example:
Both outputs are affected after two time steps but because G+*(z) = diag ( z - ~ , z - ~ ) the system can settle only after five time steps. On the other hand, for the system
R
min
5 ~ l l s+( ~1)
m,(K) ...,m,(h+P) l=1 i=l,r
+
-y ( ~
+ Ilm(6 + 1 - 1)llBl%,21
~l~)llpl~r~
G + * ( z )= diag ( z ” , 5’)
(9)
where llxll$ = x T Q x , P = controller horizon (P 1 l),s ( k ) = desired trajectory (set points), rl,BI= weighting matrices, y ( k + 1IR) = model predicted system output, and m(k) = system input. Depending on the type of time delay structure of the system, we can distinguish the following three cases. (To simplify the discussion, we will assume that there are no zeros outside the UC.) Case I. No Time Delays. Obviously in this case, there exist inputs which make the predicted system outputs agree with the prescribed set points at all future time steps. The actions of the inverse controller are found from (9) by selecting Bl = 0 and rl # 0. The minimum value of the objective function is identically equal to zero. Case 11. Equal Time Delays ( 7 ) in All Transfer Function Elements. The input actions have no effect on the outputs for 7 time steps. All the outputs can be made equal to the set points after that. Thus, the actions of the inverse controller are found by setting rl # 0 and Bl = 0 for 1 = 1,P - 7 and Bl # 0 otherwise. (The last 7 inputs m(k + P - T ) , ..., m(k P - 1)have no effect on the objective function, and a unique solution is not possible unless the corresponding weights are made nonzero.) In cases I and 11, MIMO systems behave essentially like SISO systems and no special difficulties occur. In the suggested form, the optimization problem can require
+
and therefore the outputs can settle as soon as they are affected by one of the inputs, i.e., after two time steps. Based on these observations, it is useful to introduce the be the earliest output i can following definitions: Let qmin be affected by one of the inputs, i.e., =min
7imin
I
7ij
(13)
where 7 i jare the delays in the individual transfer function elements. We call the time delays “balanced”if the system can settle as soon as the outputs are affected by the inputs, Le., if 7imin
for i = 1, r
=
(14)
Otherwise the delays are “unbalanced” and a measure of imbalance is provided by 70
= max ( T ~ +- 7imin)
(15)
1
Obviously from the definition 70
30
(16)
Thus, 7 0 measures the difference between the time an output is affected first and the time it can settle, or, alternatively, T~ is the time interval during which some of the inputs have to be prevented from moving in order to yield a decoupled response. Holt and Morari (1983) showed that the time delays are balanced ( 7 0 = 0) if and only if the rows and columns of the transfer matrix can
Ind. Eng. Chem. Process Des. Dev., Vol. 24, No. 2, 1985 487
be rearranged such that the smallest time delay in each row is on the diagonal. It is, therefore, reasonable that similar to case 11, when T~ = 0, the inverse controller can be found from the problem 9 with the weighting matrices r1 f 0, Bl,ii= 0 for 1 = 1, P - T~~~~
Bl,ii# 0 otherwise Again the weights on the inputs which have no effect on the outpub have to be nonzero to allow a unique solution. As was mentioned above, when 70 # 0, the inputs are held back artificially to yield a decoupled response. Therefore, it is intuitively clear that a decoupled response cannot be obtained by solving a least-squares problem, and, vice versa, for systems with T~ # 0, decoupling is not optimal in the sense of least squares. We will substantiate this claim below. These findings on the optimality of decoupling are new and were only possible thanks to the transparency offered by the IMC approach. (2) Formulation of the Model Predictive Problem. We concluded in the last section that if the system contains time delays, the least-squares problem 9 is not formulated conveniently for two reasons. Because of the time delays, the outputs during the first few time intervals cannot be affected at all by the present and future inputs; to include them in the objective function complicates the problem unnecessarily. Furthermore, it is desirable to obtain the "optimal" time delay predictor solution (G, = G-*-l) as a special case. When the formulation 9 is used, this is generally difficult and-as we will demonstrate laterimpossible when T~ > 0. In the following, one possible formulation is presented which removes the drawbacks of (9). Let us introduce the factorization
G ( z ) = G+*(z)G-*(z) (17) where G+*(z)was defined in (10). Next we define the time-shifted vector of future outputs y+(z): G+*-'(z)y(z)= G-*(z)m(z) (18) G-*(z) can be found from ( 3 ) G-*(z) = Hlzro + H ~ z ' o -+~ ...+ H,,z + Hro+l+ ... + H,,+Nz-~+' (19) (18) and (19) imply the following relationship between future outputs and past and future inputs. For the nuy+(z)
lz(k"+' "] y,(k t
y + ( k + 1)=
7,+
t
Y,(k + 7; + 1) H T o + , m ( k+) H70+2m(k 1 ) +
-
according to which + T ~ ++ 1, ya not
T ~ +
k
before time k + T ~ ++ 1, etc. Thus, we conclude that because of the definition of G+*, the matrices H1, H2, ..., H, are singular and H,o+l is H2, ..., H,, express the freedom which nonsingular. H1, exists in the system to affect specific outputs earlier than others. In order to obtain a decoupled response, the inputs multiplied by H1, ..., H,, have to be set to zero. If the decoupling constraint is relaxed and a least-squares problem is solved, these inputs will generally not turn out to be zero-an indication that decoupling is only optimal when T~ = 0. With these definitions, we can repose problem 9. At the present time R, the control sequence m(k),m(R + l),..., m(R + P T~ - 1)is to be found which makes the model prediction y + ( k + 1)follow a future trajectory s+(k + 1) over a horizon P + T~ This is expressed as the solution of
+
mi(R+M,-l) i=l,r
IlmG
+ 1 - 1)112B,%,)
(22)
subject to
y+(klk) = yM+(k)+ d+(klk) = N+ro
Hp(k j=l
+ T~ - j ) + d+(klk)
... = mi(E + P + T o - 1) i = 1, r where 11x11~~ = xTQx, P = controller horizon (P2 l),rland mi(k + Mi - 1) =
Bl = time-varying weighting matrices, Mi = input suppression parameter [the number of discrete intervals into the future that input mi is allowed to vary (Mi I P+T ~ &Mi = M 5 r(P + T ~ ) ) ]yM+(k) , = future values of vector Y&), [yM+(z)= G+*-'y&)] and y+(klR) = output prediction based on the impulse response model yM+(k)and a prediction of the disturbances d+(klk). In the absence of future information on d, the best prediction of the disturbance vector is d + ( kI h) =
i
d(5)
k
>h
d(k)
k
0) and when a decoupled system response is sought. Difficulties should arise in the solution of (26) under these conditions because we know
Ind. Eng. Chem. Process Des. Dev., Vol. 24, No. 2, 1985
from our previous discussion that the decoupled solution does not minimize the sum of the squares of the errors. For example, assume that Mi = P T~ (i.e., Tu = I) and rl = I, then rhTM = A. A has to be singular when T~ > 0 because H1, ..., H, are singular. Consider the following example:
489
from UM(k),the control law is given by
+
G(z) = (1- 0.5.Z') 1
[""-' ]1:] %-z
=E,.i2'-i
m ( h- 4 ) . . .
m(k - 3 ) -
(33)
(36)
Taking z transforms and assuming an infinite truncation order, we have
An optimal factorization was found in part 2 for this system which yields the inverse control law 35. Note that 52 and A are not unique. We can choose yielding G-*(z) = G,*(z)-'G(z) m
-
(34)
HiZT0-i+1
i= 1
[: 11% 1] [::: 1
=
+
[112
0
+
t,,]z-l+
where we have made the first row of FATM identically equal to zero R = [OiI,] Then we obtain the control law
which shows how the Hjmatrices follow from Hjand that To = 1. Using the rational transfer matrix 33, we can compute the "perfect" controller easily as
However, had the transfer matrix been given in impulse response form 34 with large truncation order, computation of the inverse of this polynomial matrix would have been extremely tedious. Let us assume that Nij >> 4. Employing the discretetime domain formulation 22 with Bl = 0, rl = I, P = 1,and Mj= 2, we have that
Note that since T~ = 1, HIis singular and FATMis not of full rank. Performing a Gaussian elimination, we fiid that
m(k - 2) -
m ( k - 1)0
0
(36) and (37) are two possible realizations of the same inverse control law 35. In (36), (37), and all other realizations of (35), the input m(h) depends on both present and past disturbances. This is always the case when T~ > 0. Formulation 9 does not include past set points and disturbances, and therefore it is impossible to obtain from it the decoupled perfect controller G, = G-*-l regardless of the choice of optimization parameters. (4) Control Law Computation and Implementation. The solution of (22) provides a sequence of inputs to be implemented in the future. However, in order to compensate for disturbances, only the present inputs m(R)are implemented and the problem 22 is "solved" again at subsequent intervals. In practice, the iterative solution is not necessary, but a constant recursive relation for the inputs of the following form can be obtained
m(R) +
N+ro-l 1=1
P+ro
Dlm(k - I ) = C Elt(k + 1 - TO) = 1=1
70
C E , [ s + ( k+ - T O ) - d(k
1=1
P
and thus K = 3. When 32 is used and only m(R)is selected
+ 1 - TO)]+
C El+,[s+(k + 1) - d(k)l (38)
1=1
490
Ind. Eng. Chem. Process Des. Dev., Vol. 24, No. 2, 1985
where Dj and E, E Rrxr. Note that the control law is not only a function of past inputs and present and predicted errors determined from set point and disturbance projections. It also depends on past errors when T~ > 0. The control algorithm 38 can be easily implemented on any process computer. The procedure for computing Dj and E,is straightforward. For example, in case A-l exists, they are obtained by suitably partitioning the matrices in (30) as follows: [DID~...DN+,~-~] = bTA-'TMTATrTr\k
(39)
[E1E2.. .EP+ro]= bTA-lTMTA *rTr
(40)
bT = [I,iO]
When the Gaussian elimination algorithm suggested above is required, similar relations follow from eq 32. A computer program has been developed which carries out all these operations in a user friendly manner. Expression (38) in z transforms, the following controller transfer function is obtained m(z) = DC(z)-'(PJp(~) + N~(z)ls(z)- [Npk) + N~(l)ld(z)l(41) with the corresponding block diagram given in Figure 2. 5. Stability Theorems Property 1 indicates that in the absence of modeling errors and for a stable plant, the stability of the controller is sufficient for the overall stability of IMC. In terms of the formulation above, the roots of det [D,(z)] = 0 (42) are required to lie inside the UC for stability. In the following, several theorems will be stated which show the effects of the tuning parameters on the roots of (42). These results will then be used as guidelines for controller design and tuning. In general, the tuning parameters can be assumed to influence the dynamics of the closed-loop system in a similar way as their counterparts in the SISO algorithm. Let us first show how the exact inverse G-*(z) can be obtained from (22) or equivalently (26) in a very simple manner. Theorem 1. Let rr= I, Br = 0 , l = 1,P + T ~ M, , =P T ~ and , j = 1, r. Then the solution of (22) yields the perfect control law G, = G-*(z)-'. Proof. See Appendix section. The next result shows that for any system, a stable controller can be obtained by making the matrices B1sufficiently large. Theorem 2. Assume det Br # 0 , l = 1, P + 70. Then there exists a finite p* > 0 such that if llTMrB%TMll> /3*, the IMC controller is stable. Proof. See Garcia (1982). Another way to obtain a stable control law is to increase the horizon. Theorem 3. Assume = I,, Br= 0,1= 1,P + 7,,. Then for a sufficiently small M and a sufficiently large P > N chosen subject to the condition that rank [I'ATM] = M, the controller 38 is stable. Proof. See Garcia (1982).
+
6. Tuning Procedures According to theorem 1,the tuning parameters can be selected so as to obtain the perfect controller from (22). This control law should be tried first if it is stable. Should
the perfect controller be unstable or exhibit extreme rippling, other parameter settings and consequently different approximations to G-*(z)-l have to be tried to bring the
Figure 2. Block diagram of multivariable IMC with predictive cptroller. Offset compensator: GOFF(Z)= [G(l)GC(l)]-' = IG( 1)DC(-)-'[Np(1) + N ~ ( 1 )I-'. l
roots of (42) inside the unit circle or closer to the origin. As in part 1,we should distinguish between two cases. Case A. G-*(z) has transmission zeros outside the UC and consequently G-*(z)-l is unstable. Case B. All the transmission zeros of G-*(z) lie inside the UC. Since computation of the characteristic polynomial is quite tedious, a computer routine is imperative. A convenient procedure is to triangularize D,(z) in order to compute its determinant. An algorithm to perform the necessary manipulations is given by Kucera (1979). Once obtained, the largest root of the resulting polynomial in z is found by any suitable procedure. (1) Sampling Time (T).Case A. A continuous system with RHP zeros produces a G-*(z) with roots outside the unit circle for T small. Intuitively and from the SISO results (Garcia and Morari, 19821, it is to be expected that the zeros outside the UC disappear as T is increased. Case B. Stability of IMC is not affected by T. However, a small T tends to produce strong variations of the manipulated variables. On the other hand, high sampling frequencies are preferable in order for IMC to be able to handle frequent disturbance changes. The optimal value of T should be selected according to the particular application. (2) Input Suppression Parameters (Mi).Case A. As shown by theorem 3, Mi's small relative to P produce a stable controller. Case B. Extreme excursions of the manipulated variable are reduced by using small Mi's. (3) Input Penalty Matrices (Br). Case A. By selecting Brnonsingular and sufficiently large, stability of IMC is ensured for any system, as shown in theorem 2. Case B. A sluggish controller response is produced by increasing B1.In addition, D,(l)-'[Np(l)
+ NF(1)] #
G-*(l)-'
which leads to offset. This is corrected by the multivariable offset compensator shown in Figure 2. (4) Output Penalty Parameter (I'J. Individual output responses can be influenced directly by suitably selecting the elements in the matrices rl. For example, a large relative weight on a particular output decouples its response from the rest of the outputs. This will be demonstrated later in an example. As for SISO systems, time-varying r/s can improve the response. Unfortunately, unlike the SISO case, for multivariable systems, guidelines for the "optimal" selection are not available to date. (5) Optimization Horizon ( P ) . Case A. As shown in theorem 3, any system can be stabilized by increasing P sufficiently. Case B. As P is increased, the controller is generally made more cautious. 7. Simulation Examples The IMC formulation introduced in this paper allows
No. 2, 1985 491
Ind. Eng. Chem. Process Des. Dev., Vol. 24,
the designer to obtain as the controller either the exact model inverse or an approximate inverse whose properties depend on the particular choice of tuning parameters. The following examples demonstrate the effect of these parameters on response and stability of the controlled system. As we shall see, the designer has much freedom in shaping the system dynamics according to the needs of the particular application. (1) Shaping of Dynamics: Wood/Berry Column. In part 2, the transfer matrix describing the response of a methanol/water column (Wood and Berry, 1973) was given. When a sampling time of 1 min is used, the discrete-time transfer matrix has the form
1
G ( z )=
The diagonal time-delay factorization matrix 10 for this system is
Because T~ = 0, this factorization is optimal in terms of the sum of the squares of the errors or in terms of the response time. The perfect controller for this system is found by setting P = Mj = 1, B1= 0, and rl = I. In Figure1 3 we show the response to a step change in the overhead composition set point for the perfect controller with an exponential filter F=diag(-)
1 - “i
e
19
2e
3e
e
ie
28
30
50
40
58
HII4UTES
Figure 3. Wood/Berry column. Change in overhead composition set point (-) P = Mj = 1,rl = I, B, = 0, ai = 0.8;(---) P = Mi = 1, rl = I, B = diag (1.25,3.0);(-.-.-) P = 10, Mj = 2, r i = I , B = 0.
-
-
1
0
5
0
ie
15
20
5 1s
1 28
i=l,2
1 - aiz-1 (ai = 0.8) added. This perfect controller produces an exponential approach to the set point and a perfectly decoupled response, since G+* and F are both diagonal. Instead of the filter input weights, B = diag (1.25,3.0) can be used to obtain a smooth approach of the overhead composition to its set point. However, since the input moves are restrained, decoupling is lost. As an alternative to penalizing the inputs, the number of input moves can be made smaller than the horizon ( P = 10, M j = 2, Bi = 0, and r i = I) which has a similar effect. Figure 3 shows that the overhead reaches the set point fast, but due to the input restriction, closed-loop interactions arise. In this case ( T ~= 0), the perfect controller can also be found without factoring the delays by solving (9), with P = 4 = M j , rj = I, and B = 10-21. Any P < 4 does not produce the desired response. Note that now an 8 X 8 matrix needs to be inverted to compute the control law as compared to a 2 X 2 matrix previously. Also, the proper B needs to be selected by trial and error. Though it is clearly possible to shape the response by using the tuning parameters as we can see from this example and though it is done in practice (cf. Dynamic Matrix Control), we do not recommend this procedure unless it is absolutely necessary, for example, when severe rippling occurs or when the inverse is unstable. We feel that in all other cases, the filter time constant is a much simpler and much more direct tool to shape the response and to dampen severe input fluctuations. (2) Effect of Time Delay Factorization. This example will demonstrate that IMC allows the designer to make conscious decisions on the trade-offs between speed of response and decoupling which are involved in time delay compensation.
-
e 0
.
5
5
10
TIME
Figure 4. Effect of factorization: (-) perfect controller with delay factorization; (- - -) without factorization with P = 5, Mj = 5, I’ = I, B = W31;(- - -) without factorization with P = Mj = 5, r = diag U O , ~ ) ,B = 10-31.
If we artifically alter the Wood/Berry model time delays z- h
G ( z )=
[
(2)
1
2-5 hi2( z )
~ - ~ h , , ( z~)- ~ h ~ , ( z )
the time delay factorization yields G+* =
:[
.-PI
For this case ( r 0 = 2), G+* is not “optimal” and consequently there is freedom in reducing the output delays. When P = 1, Mj = P + T~ = 3, ri = I, and Bi = 0 are selected, the perfect controller is obtained after performing the Gaussian elimination procedure described previously. The responses to a change in the set point of output y1are shown in Figure 4. Note the expected delay of four time units and the absence of interactions. The decoupling constraint can be relaxed, and a control law which meets a minimum square error criterion can be found either from (9) or (22). For the simulations shown here, the delays were not factored and (9) was used. Be-
492
Ind. Eng. Chem. Process Des. Dev., Vol. 24, No. 2, 1985 2 ,
I
I
I
I
-, , ,,
,
I
,
,
.
SECOMS
SECONDS
Figure 5. Hydrogen/oxygen fixed-bed reactor with sample time of T = 20 a: (-) exact inverse; ( - - - ) P = 20, M j = 3, r f = I, Bf = 0; (- - -) state deadbeat.
Figure 6. Hydrogen/oxygen fixed-bed reactor with sampling time of T = io 8: (-1 P = 1 = M ~r,= I, B = 0.021; ( - - - ) P = 20, M~= 4; r = I, B = 0; (- - -1 state deadbeat.
cause the smallest delay in the system is 1, the inherent sampling delay, (22) does not offer any storage advantages. When T M = I, problem 9 does not have a unique solution (det A = 0) unless det Bl # 0; 1 = 1, P. Figure 4 shows the optimal responses for P = 5, M, = 5, B = 10-21,rI= I, and rl = diag (10,l). In both cases, a faster approach to the set point of y1 is obtained than with the decoupling constraint. However, the strong interactions which are observed in y2 can be undesirable for some practical applications. (3) Fixed-Bed Reactor. A hydrogen oxidation fixedbed reactor system model (Foss et al., 1980) with seventh order dynamics was described in part 2. For a sampling time of T = 20 s, the perfect controller for this delay-free system is stable. The response of the perfect IMC controller to a step in the outlet temperature (TOUT) set point is given in Figure 5. Note that both TOUT and COUT (outlet concentration) are at 1.0 and 0, respectively, at every sampling interval. However, significant oscillations occur between the sampling times. The reason is that the inverse controller forces only the outputs to zero at the sampling times but not the states. The oscillations can be dampled by restricting the input moves by choosing P = 20, MI = 3, l' = I, and B = 0 (Figure 5). Also, one can select the following: P = 7 = M,, I'l = 0, 1 = 1,3; rl = I, 1 = 4,7; Bl = 10-21,1 = 1,4; Bl = 1031,1 = 5, 7. Similarly to the SISO case discussed in part 1, this parameter choice produces a state deadbeat controller (Reid et al., 1979) which makes the system settle in at most seven steps. The response is also shown in Figure 5. Note that the faster response is accomplished a t the cost of larger excursions in TOUT and stronger interactions. For a sampling time of T = 10, the reactor exhibits nonminimum-phase characteristics and therefore the perfect controller is unstable. The system can be stabilized by either increasing B or decreasing Mirelative to P. The responses to a step in the set point of COUT are shown for both choices in Figure 6. Note that the M, reduction does a better job here in terms of eliminating interactions. We can also try the state deadbeat controller for this sampling time. Again, stronger interactions and larger overshoots are observed. Contrary to what we observed for SISO systems in part 1, the state deadbeat controller does not necessarily produce responses which are acceptable in practice. (4) Effect of Zeros Outside the UC. As was discussed
in part 2, zeros outside the UC can be factored and included in G+2
-
(3%) = G+i(z)G+2(z)G-(z) so that G--'(z) is stable. In this paper we have only discussed the factorization of time delays, and, therefore, the inverse of C+,(z)-'G(z) could be unstable. As explained above, one can stabilii such systems by the proper choice of tuning parameters, e.g., penalizing the inputs. However, it is clear that since ~ ( 2 = )
G(z)Gc(z)ys= G+(z)ys
no stable Gc can cancel the zeros of G and thus they have to show up in the closed-loop response. The following system which has one zero outside the UC a t E = 1 + (0.3)'12 was considered in part 2:
r
1
Z -
G(z)=
1
0.6 0.5 -~
0.4 0.5
1.2(-) Z -
0.5
z - 0.5 0.6
Z -
0.4
It was found that if decoupling is imposed (i.e., by choosing G+ diagonal), then this zero must appear in both output responses. However, if one-way coupling is allowed, the inverse response behavior can be limited to only one of the outputs. When the least-squares formulation (9) or (22) is used and the tuning parameters are selected appropriately, the response can be shaped as desired. Obviously it is impossible to obtain a perfectly decoupled response in this manner because decoupling would introduce another zero outside the UC, and this is not optimal in terms of least squares.
But through appropriate tuning, the effect of the zero of G can be shifted to different outputs. The system above has very fast response, and therefore P = 20 is for all practical purposes "iafinit"horizon length. I t follows from optimal control theory (Kwakernaak and Sivan, 1972) that if we select P = 20 = Mj, this system can be stabilized with any input weight; in fact B = 10-41
Ind. Eng. Chem. Process Des. Dev., Vol. 24, No. 2, 1985 493 I
I
l
l
-
1
1
1
1
50
0
-2
I
I
I
l
l
I
I
I
l
l
1
I
I
I
I
l
l
1
280
150
100
I
I
I
,
1
1
1
I
l
l
coupled and decoupled controllers with equal ease. In part 3, we concentrated on the design of the controller which is, in general, a stable approximate inverse of the process model. The impulse response model which we employed appears to make the inversion particularly simple and allows to affect the characteristics of the approximation transparently through the choice of tuning parameters. Furthermore, it yields directly a recursive expression for the controller which is readily implemented on any process computer. An interactive computer program has been prepared through which the controller can be easily calculated by starting from the impulse response model and the tuning parameters. This design program is coupled with a simulation program where the effect of the tuning parameters on the response can be observed.
1
Acknowledgment Financial support from the National Science Foundation, the Department of Energy, and Shell Development Co. is gratefully acknowledged.
produces a stable controller. In order to demonstrate the shifting of the zeros of GI different output weights will be selected. For rr = diag (1.0, lO.O), deviations in y 2 are severely penalized. Thus, the resulting closed-loop response should have a “triangular” structure. G+=
3
This is shown in Figure 7. For a step in the set point of y l , interactions are minimal, while for a step in the set point of yz extreme interactions are observed. Also severe inverse response in y 1 is observed. Therefore, by penalizing y a , we have selectively shifted the less desirable dynamic characteristics to y l . The opposite occurs if rl = diag (10.0, 1.0) as shown in Figure 7, where now y 2 exhibits strong interactions and inverse response. 8. Conclusions IMC consists of three parts: (1)internal model to predict the effect of the manipulated variables on the output; (2) filter to achieve a desired degree of robustness; and (3) control algorithm to compute values of the manipulated variables based on present and past errors and set-point trajectories. This structure has important theoretical and practical advantages: (1)The internal model gets rid of the issue of closed-loop stability. Even for nonlinear controllers (input saturation), the stability is guaranteed. (2) A model-based prediction of the output is attractive in the presence of output constraints. Any constraint violations can be anticipated, and corrective action can be taken. (3) The filter allows simple on-line tuning of multivariable controllers by the operating personnel. (4) Through the fiiter, the robustness of the system can be affected directly and not in a roundabout fashion through the choice of weighting matrices-as is the case for the linear quadratic Gaussian problem. (5) The structure and parameters of the “optimal” controller are known a priori: it is the inverse of the invertible part of the system model. This “target prediction” makes it simple to find suitable approximations to the inverse for practical implementation. (6) IMC automatically includes an optimal multivariable time delay compensator. (7) IMC allows one to obtain
Appendix Proof of Theorem 1. The perfect controller
m(z) = G-*-’(z)t(z)
(AI)
yields a response Y(Z)
G+l*(z)s(z)+ [I - ( 3 + 1 * ( ~ ) ] d(z)
This response is decoupled and the sequences of inputs which produce it are unique. These inputs must satisfy e(z) = G-*(z)m(z) or ~ ( k I) = Hlm(k 70) H,,+lm(k)
+
+
+ ... + H,,m(k + 1) + + ... + H,,+pm(k + P - 1) (A2)
For the tuning parameter selection rr= I, BI = 0; 1 = 1,P r0,problem 22 is equivalent to solving the system of equations
+
EP(i
t
1)=
I iI 1 € ( h+
l€(h
1)
= y M + ( k + 1)=
P)
AU(h) t $V(k-
I) (A3)
(cf. eq 34) which is just (A2) written for time k = k - r0 + P - 1. Without loss of generality, let us assume d = 0. Now consider a step change in set point from 0 to s. Then (A3) becomes
to k = lz
As explained previously, for 7o > 0, A is not of full rank and (A41 has generally no solution. Because HI, ...,H, are singular and H,+l is nonsingular, it can be easily shown
Ind. Eng. Chem. Process Des. Dev. 1985, 24, 494-498
494
that the.last PP rows of A are linearly independent. Thus, there exists an infinite number of solutions to (A4) which satisfy the last r.P equations. By carrying out a Gaussian elimination on the top and permuting the elements of U(L) if necessary, we obtain -
0
I 1
I
A HTo+i
I
. . . H,
U ( Z )+
Hi
-HTo+p . . . . . . . . . .
Q V ( Z - 1) (A5)
where the nonzero last rows of the ney matrix h’ are linearly independent. All solutions U ( k )which satisfy the nontrivial bottom half Of the eq A5 produce the decoupled response y i ( k ) = Si k 2
~ i +
+1
But this is the response of the perfect controller, and since
it is unique, it must be true that all these solutions are equal to the perfect control law.
Literature Cited Foss, A. S.; Edmunds. J. M.; Kouvaritakis, B. Mol. Eng. Chem. Fundam. 1980, 19, 109. Gantmacher, F. R. “The Theory of Matrices”; Chelsea: New York, 1959; voi. 1. Garcia, C. E. Ph.D. Thesis, University of Wisconsin-Madison, 1982. Garcla, C. E.; Morari, M. Ind. Eng. Chem. Process Des. D e v . 1982, 2 1 , 308. Hoit, 6.; Morari, M. Chem. Eng. Sci., in press. Kucera, V. “Dlscrete Linear Control”; Wiley: New York, 1979. Kwakernaak. H.: Sivan. R. “Linear ODtimai Control Systems”; Wiley-Interscience: New York, 1972. Penrose, R. proc, Cambridge phi/os, sot. 1955, 51, 406, ReM, J. G.; Mehra, R. K.; Kirkwood, E. Proc. IEEE Conf. Dec. Control 1979, 307. Wood, R. K,; Berry, M, W. Chem, Eng. sei. 1973, 28,1707,
Received for review April 25,1983 Revised manuscript received May 23, 1984 Accepted July 23, 1984
COMMUNICATIONS Improved Algorithm for Obtaining Liquid-Liquid EquMibrium Parameters by Use of a Penalty Term An algorithm is proposed to improve optimization on liquid-liquid equilibrium data to obtain liquid solution model parameters. The method increases the efflciency of the Marquardt technique coupled with a penalty function. The algorithm was tested with the UNIQUAC and NRTL liquid solution models.
Introduction Due to the increased interest in extraction and multiple-phase distillation, numerous papers concerning liquid-liquid equilibria have appeared in the recent literature. Several molecular thermodynamic models, for example NRTL (Renon and Prausnitz, 1968), UNIQUAC (Abrams and Prausnitz, 1975), and UNIFAC (Fredenslund, Jones, and Prausnitz, 1975), can be used for LLE calculations. The predictive ability of these semiempirical representations of the liquid state is strongly dependent on the method used to determine the model parameters from experimental data. Simonetty et al. (1982) have shown that regressing LLE tie line data alone yields the best model parameter set for type I liquid-liquid systems. LLE correlations with these models are quite nonlinear because of the exponential interaction energy term. Regression on experimental tie line data with a standard numerical technique will yield parameters that give a good fit to the experimental data included in the regression set. However, in subsequent calculations unstable liquid phases may be obtained. Ma,ny techniques to avoid this have been developed. Serrensen (1980) adds a penalty term to the objective fuhction to prevent the model parameters from assuming
large numerical values. This work proposes an algorithm to automatically adjust the penalty function to improve optimization of the model parameters.
Criterion of Equilibrium Three conditions must be met for an LLE system (Sarensen et al., 1979). First, the material balance must be preserved. Second, the chemical potentials of each compound in the two liquid phases must be equal. Finally, the Gibbs energy must be at a minimum at the system temperature and pressure. The last restriction is a necessary and sufficient condition, whereas the second restriction is only necessary but not sufficient. The first and second restrictions are used due to their simplicity and are normally written as
EX;’= 1 i
(1)
Since these are necessary but not sufficient conditions, computations based on eq 1 through 3 can lead to false
0 1985 American Chemical Society Q196-43Q5/85/1124-Q4945Q1.5Q/Q
