I n d . Eng. C h e m . R e s . 1987,26, 577-581 Denbigh, K. G.; Turner, J. C. R. Chemical Reactor Theory, 3rd ed.; Cambridge University Press: London, 1984. Fournier, C. D.; Groves, F. R. Chem. Eng. 1970,2, 21-25. Levenspiel, 0. Chemical Reaction Engineering, 2nd ed.; Wiley: New York, 1972. Millman, M. C.; Katz, S. Ind. Eng. Chem. Process Des. Dev. 1967, 6, 447.
577
Rosenbrock, H. H. Comput. J . 1960, 3, 175-184. Scott, S. K . Chem. Eng. Sci. 1983, 38, 1701-1708. Wnek, W. J.; Trvlarides, L. L.; Haas, W. R. AZChE J . 1974, 20, 707-712.
Receiued for review March 25, 1986 Accepted September 17, 1986
Internal Model Control Using the Linear Quadratic Regulator Theory Mario A. Rotea and Jacinto L. Marchetti* Institute of Technological Development for the Chemical Industry-INTEC Nacional del Litoral), 3000 Santa Fe, Argentina
(CONZCET, Universidad
The Linear Quadratic Regulator theory is applied to the design of Internal Model Control for single-input, single-output systems. The proposed technique provides an easy and sound alternative to obtain the best stable inverse of the plant model that limits the achievable control quality and yields a stable open-loop controller which is optimal when the model fits the plant exactly. Constraints in the control variable can be handled by a single tuning parameter that is adapted online to yield maximum performance to set-point changes. This study demonstrates that when open-loop controllers are used, like IMC, in the presence of constraints, a better performance is obtained if the controller is adjusted to account for the constraints instead of ignoring them. 1. Introduction Internal Model Control, IMC, was developed by Garcia and Morari (l982,1985a,b)following initial suggestions by Brosilow (1979) and Palmor and Shinnar (1981) on the use of Smith predictors. The motivation behind the IMC technique is the need of an effective control structure for chemical processing units. The basic IMC structure for open-loop stable plants, G,(z), is shown in Figure 1. The design of G,(z), the IMC controller, is equivalent to the design of an open-loop controller subject to the following conditions: (1) G,(z) must be realizable; (2) G,(z) must be stable;_(3) 1 - G,(l)G,(l) = 0; (4) 1 - G,(z)G,(z) f 0, where G,(z) is the transfer function model. The first three conditions have been extensively discussed by Garcia and Morari (1982,1985a). However, the last condition has never been presented before and implies structural realizability of the IMC loop. The structural realizability can be analyzed on the IMC loop indicated in Figure 1. Two relationships can be written from the block diagram 4 2 ) = G,(z)[r(z)- &)I (1) d ( z ) = y ( z ) - d,(z)u(z) (2)
Equations 1 and 2 can be rearranged as a linear system having two equations and two unknowns (3) G,(z,u(z)
+ &z)
= y(z)
(4)
Since loop realizability is required, the following con, dition must be satisfied: given a pair of inputs y ( ~ )r(z), there _must be a unique solution for the pair of variables u(z), d(z). More clearly, this requires
or 1 - G,(z)d,(z) f 0
(5)
In this work, the Linear Quadratic Regulator (LQR) theory is used to find the IMC controller G,(z). The development is made by applying the LQR theory to a nominal “internal model” of the plant, where the state space coordinate system is chosen so that the state vector is the input to the system,, i.e., the manipulated variable. It will be shown that this particular realization of the internal model allows the design of an optimal open-loop controller in a simple manner. This open-loop controller is then fitted into the IMC structure, and the conditions mentioned above are verified. Finally, the proposed design technique is adapted to accept constraints on the control variable while maintaining optimality. 2. Optimal IMC Design 2.1. Plant Representation. A finite weighting sequence model is used to describe the plant dynamic
The advantages of using the above representation for modeling stable chemical plants have been exposed by Garcia and Morari (1982,198513)and Richalet et al. (1978). The analysis including time delay in the plant model is discussed later in section 2.3. The following development is for undelayed plants only, i.e., the case where hl # 0. In order to apply the LQR theory for the derivation of the IMC controller, G,(z), the following state space model is used to realize the transfer function in eq 6: x(k + 1) = Ax@) + bu(k) (7) Y ( k ) = cTx(k)
(8)
where x ( k ) E RN,since X y k ) = [u(k - N) u(k - N
* Author t o whom
all correspondence should be addressed.
+ 1) ... u(k - l ) ]
and A E R N x Nhas the canonical form
0888-5885/87/2626-0571~0~.5~/0 0 1987 American Chemical Society
(9)
[; 1
;]
578 Ind. Eng. Chem. Res., Vol. 26, No. 3, 1987 0
A =
...
0 1
0
0
............................................ 0 ... 0 1 0 0 ... 0 0
, IMC loop ------------(10)
The control vector b is given by
bT = [0 0
......... 0
11
E RN
(11)
while the dynamic information of the plant is in the observation vector, cT =
[ h N hN-1 h N - 2
... h,]
(12)
Note in eq 9 that the state variables are the last N inputs to the system. This particular realization allows that any state feedback control law designed on the basis of these variables to result in an open-loop controller. Note also that the system (A,b,c)defined in eq 10-12 is controllable and observable. 2.2. Linear Quadratic Regulator. The open-loop control law to be used in the IMC structure is obtained by solving a standard quadratic optimization problem. Like in other controller designs,j t is first assumed that there is no modeling error; i.e., G,(z) = G,(z). For simplicity, it is also assumed that the unmeasurable disturbance, d, is zero. Hence, the basic control objective is to force the model output, 9, to follow the set-point variable, P, while preventing excessive movements in the control variable. Then, by use of the realization defined in eq 7-12 to describe the model, the problem to be solved can be written as
2
min { [ g ( k + 1)- FI2 u(k),k a o k = O
+ p k=O 5 [ u ( k )- GI2)
(13)
subject to
x(k
+ 1) = Ax@) + bu(k)
(14)
g ( k ) = cTx(k)
(15)
G;W
(16)
where ii =
Figure 1. Basic IMC structure.
where the state vector, x ( k ) , has been already defined in eq 9. It is appropriate to confirm now that the design technique suggested here yields a controller that satisfies the conditions listed at the beginning of section 1. After taking z transforms in eq 9 and 20, the following transfer function is obtained: G,(z) = G;l(l)[l
where -yT(z) = [ z - z-N+l ~ ... z-']. Inspection of eq 21 leads to the following conclusions: (i) The controller defined by eq 21 is realizable since the 1term in the denominator confirms this property (Luyben, 19731. Moreover, eq 21 allows structural realizability since 1 - G,(z)G,(t) f 0 is apparent. (ii) The controller defined by _eq 21 satisfies the zerooffset conditions since G,(l) = G,,-l(l). (iii) The controller defined by eq 21 is stable for all p > 0. This property is the main theoretical support behind the technique proposed to design G,(z). Since the system (A,b, c )defined by eq 7-12 is observable and controllable, the optimal control law in eq 17 is asymptotically stable for all p > 0 (Kwakernaak and Sivan, 1972). 2.3. Time Delay Systems. Many real systems have time delays. Therefore, let us consider the following process model:
This problem can be reduced to the LQR problem for infinite time (Kwakernaak and Sivan, 1972). Then, after some rearrangements the solution is given by u ( k ) = -fTx(k) + [l + fT-y(l)]&p-l(l)r(k) (17)
where ~ ~ (=1[l) 1 ... 11, r(k) = FH(k), and f calculated from
fT
=
[p
C C ~(19)
Since the system (A,b,c) defined in eq 10-12 is controllable and observable, the solution P always exists for
> 0.
Next, to account for modeling errors as well as unmeasurable disturbances, the open-loop control law defined in eq 17 must be included into the IMC structure. Note, from Figure 1,,that to accomplish this task the estimated disturbance, d ( k ) ,must be added to the set-point variable. Therefore, the IMC control law proposed in this work can be written as u ( k ) = -fTx(k)
i=l
hl # 0
(22)
where r is a positive integer defined by
rT, d time delay < (7 + 1)T,
(23)
The performance function, for this case, can be written (18)
The P matrix in eq 18 is a positive semidefinite matrix obtained by solving the algebraic Riccati equation
p
N
G,(z) = z-'Chiz-'
E RN is
+ bTPb]-'bTPA
P = ATPA - ATPb[p + bTPb]-'bTPA +
+ fTy(l)][l + f T - y ( ~ ) ] - ~(21)
+ [l + fT~(l)](?p-~(l)[r(k) - d(k)J
(20)
as m
k =O
[g(k +
m
+ 1)- TI2 + pC [ ~ ( k-) i i I 2 k=O
(24)
Note that the change of variable f ( k ) = 3(r + k) gives the optimization problem m
m
min { [9(k + 1) - TI2 + pC [u(k)- iiI2J (25) k=O u ( k ) ,k 3 0 k=O subject to N
Y ( k ) = Ch,u(k - i) i=l
(26)
Equations 25 and 26 define a problem completely equivalent to the undelayed systems. This fact allows us to conclude that G&) can be designed assuming 7 = 0.
Ind. Eng. Chem. Res., Vol. 26, No. 3, 1987 579 u2 I
U,
U,
Figure 3. Sketch of control actions for set-point changes. Figure 2. Basic IMC structure including constraints on the controller output.
3. Constraints on t h e Control Variable Constraints on the control variable are found very frequently in process control problems. They arise due to actual physical limitations such as the limits of a control valve actuator or due to operating or safety margins desired by the practitioner. Garcia and Morari (1985a) have already indicated that constraints on the control variable do not affect the closed-loop stability when (i) the control structure has the configuration shpwn in Figure 2 and (ii) there is no mismatch between G,(z) and G (2). It can also be noted that if (i) and (ii) are satisfied, tke controller G,(z) is not sensitive to the constraints. Then, even in the case where the stability is guaranteed, the risk of an unpredictable control quality is present when the first two or three control actions are chopped by the constraints. This is not recommendable for most chemical processes. This section proposes a method based on adaptive controller gains computed online to satisfy input restrictions. The following discussion assumes that the controller gains are updated any time a set-point change is made, but the application can be readily extended to the case where the estimated disturbance, d ( k ) , jumps significantly. Hence, in the following development it is assumed that d ( k ) = 0; i.e., there is no mismatch between the plant and model, and the external disturbance, d, is zero. It can be shown from previous studies (Kwakernaak and Sivan, 1972; Priel and Shaked, 1983) that when the weighting factor, p , goes to zero, the controller of eq 21 tends to the r e e a b l e part of the inverse of the minimum phase image of G,(z), i.e., the perfect controller defined by Garcia and Morari (1982). Typically, this controller produces excessively energetic control movements and, most frequently, violates realistic restrictions. The opposite limit for the weighting factor, Le., p going to infinity, produces a step change in the controller output when a set-point change is made. Observe, from eq 18 and 19, that if p goes to infinity the gain vector f tends to zero. Hence, the control law takes the form
u(k) = G;l(l)r(k)
(27)
i.e., the controller output jumps from the old stationary value to the new one, showing neither overshoot nor oscillations. Therefore, for a given set-point change, an intermediate p value can be found in order to obtain an optimal sequence u(k),k 3 0, that yields the fastest plant response without violating the constraints on the control variable. From eq 20 it can be shown that the first control action, after a set-point change, AT, satisfies Au(0) = [ l
+ fTr(l)]Aii
(28)
where Aii is the difference between steady-state control values; Le., Aii = a,,, - iiold.
It has been observed in many simulation exercises that when the control law in eq 20 is used, the first control action is the one that exceeds the new stationary value most, i.e., (29) ( ~ ( 0-) a,,l 2 lu(k) - anewl k 3 0 Note that, as long as the condition in eq 29 is satisfied, the constraints will not be violated if the first control action, Au(O), is appropriately prescribed. For instance, if the control variable follows the pattern indicated in Figure 3, Au(0) must be chosen such that lAu(O)I =
(AG( + min (ua - anew, u,,, - u l Jif u1 < anew< u2
iIAal if anew1
u2 or
a,,,
5 u1
(30)
where u1 and u2 are the limiting control values. Therefore, the problem is to find the p value and the optimal gain vector, f, that yield the desired Au(0). Shaked (1979) and Priel and Shaked (1983) have presented a useful relationship R(z-')(p b p b ) R ( z )= p + Gp(z-')Gp(z) (31)
+
where P is the solution of the Riccati equation and R ( z ) is the return difference transfer function given by R ( z ) = 1 + P ( z I - A)-'b (32) Substituting eq 10 and 11 into eq 32 gives R ( z ) = 1 fTr(z) (33)
+
and from eq 28, 31, and 33, the following result can be obtained: P p
+ Gp2(1)
Au(0)
+ b T b =[-E]
(34)
The importance of eq 34 becomes apparent when it is solved simultaneouslywith eq 19 for a given value of Au(0) obtained from eq 30. The result is a gain vector f that produces the fastest optimal plant response without violation of input restrictions. Notice that the simultaneous solution of eq 19 and 34 does not lead to a proper solution, P 3 0, p > 0, for all arbitrary values of Au(0). This is equivalent to indicating that 1 + fTr(l), in eq 28, is restricted to the results obtained for positive values of p only. Therefore, a prescribed Au(0) value yields an optimal solution only if (35)
where = lim [I + fTr(1)l
(36)
P-O+
The condition in eq 35 is not an obstacle. Actually, it only indicates that there is a maximum admissible value for the prescribed Au(O), which can be calculated offline from eq 36 for each particular system.
580 Ind. Eng. Chem. Res., Vol.-26, No. 3, 1987
Note that if eventually eq 29 is not satisfied, the whole sequence of control actions that follows a set-point change can be checked through the open-loop control law in eq 20, previous to executing the change and with little additional effort. The logic to prescribe lAu(O)(that suits a particular application best can be determined by using additional information obtained by simple inspection of the system behavior. For instance, for the second set-point change in Figure 3, it would be possible to prescribe (Au(O)(= (Aiil + anew- u1 instead of (Aii(O)( = IAiil + u2- i,, given by eq 30, if the next control action satisfies lAu(l)l C (Au(0)I u2 - sold. In other words, the user may propose alternative methods to specify Au(0). The routine to compute the gain vector f remains unchanged.
+
0
40
120
80
160
200
240
160
200
240
TIME
-10
"
~
0
40
80
I 120
TIME
4. Computational Algorithm
Today microcomputers are capable to execute the above control technique and to adjust the controller gains right after a set-point change is commanded. Summarizing,the essentials of the proposed algorithm are included in the following steps: (i) Compute the change of the stationary value of the control variable, Aii = (?C1(1)AF (37) (ii) Compute lAu(O)I according to eq 30 or a user's proposition. (iii) If (Au(O)I = IAiil, then f = 0 . (iv) If lAu(O)I 2 ~ $ ~ J A i i l then , f = f+,,. (v) If IAiil C lAu(O)I C C$mm(Aiil,compute the gain vector by using the following recursive method: 1. Take k = 0 and P(0) = c.c*. 2. Compute recurrently p = [d,2(1) - 92brp(k)b][92- 11-1 (38) where C$ = Au(O)/Aii and P(k + 1) = ATP(k)A - ATP(k)b[p + bTP(k)b]-'bTp(k)A
+C C ~ (39)
until
+
IIP(k + 1) - P(k)ll < E 1). 3. Compute the gain vector by using
Then P = P(k eq 18. The simplicity of the matrix A and the vector b defined in eq 10 and 11, respectively, leads to important simplifications in the computation of p , P, and f when eq 39 and 18 are written down entry by entry. The resulting routine does not need matrix inversions, and the memory requirement is small. The convergence of the above algorithm is guaranteed. Recall from section 2.1 that the system ( A , b , c )defined in eq 10-12 is complete; hence, the difference Riccati equation, eq 39, converges as long as p > 0 (Kwakernaak and Sivan, 1972). Since 4 satisfies eq 35, only positive p values are generated by the algorithm. Furthermore, large amplitude disturbances might trigger the adjusting rouiine, without operator intervention, when the signal Ir(k) - d(k)l overpasses a selected minimum value at any sampling instant. However, additional steps must be included in the above routine in order to avoid too frequent changes of the controller gains. 5. Simulation Results
In order to demonstrate the characteristics and the quality of the proposed technique, several simulation examples were run on a digital computer. A number of common conditions were established for simplicity; i.e., the
I0
">
I 1-1
I
Table I. Simulation Results for Three Different Plants IAE ISE ITAE (a) Overdumped Plant, GJs) = e-208/(10s+ 1)2;See Figure 4 automatic tuning 68 68 7970 81 75 9813 const p = 0.02 (b) Underdumped Plant, G,(s) = e-205/(100s2+ 12s + 1); See Figure 5 automatic tuning 63 64 7219 77 73 9371 const p = 0.02 (c) Nonminimum-Phase Plant, G,(s) = [(I - 5s)e-20s]l/(10s+ 1)2; See Figure 6 automatic tuning 78 81 9182 87 85 10603 const p = 0.02
model order ( N = lo), the sampling time ( T , = 10 units), the constraints on the control variable (-1, l ) , and the precision required to compute the optimal gains (E = lo+) are the same in all the examples. The same test is applied to three different plants that cover most of the dynamic behaviors typically found in process control. The simulations start from a steady-state condition with the setpoint value r(0) = 0.0; then, two set-point changes, to r(10) = 0.5 and to r(120) = -0.9, are commanded. Table I and Figures 4-6 compare the performance obtained when automatic tuning is used to satisfy input constraints with that obtained for p values that produce saturation on the control variable. The parameter NI gives the number of iterations of the algorithm of section 4, step v, to obtain convergence within precision t. Inspection of these results indicates the benefit of adjusting the optimal gains to avoid chopped control actions. This study demonstrates that when open-loop controllers are used, like IMC, in the presence of constraints, a better performance is obtained if the controller is adjusted to account for the constraints instead of ignoring them. 6. Conclusions A new method is proposed to design an optimal openloop controller G,(z) that fits into the IMC structure developed by Garcia and Morari (1982). The most important features of the design technique presented here are the following:
Ind. Eng. Chem. Res., Vol. 26, No. 3, 1987 581 - aulomalic
I)
--_p = 0 0 2
10
0
80
40
120
160
variable to the constraints in an optimal fashion. A simple technique to accomplish this task through online automatic tuning is presented, and the benefits are demonstrated by using digital simulations of typical plant dynamics in process control.
luning
200
Acknowledgment
240
TIME 10
,
I I/ )
This research was supported by Universidad Nacional del Litoral and Consejo Nacional de Investigaciones Cientificas y Tecnicas (CONICET), Argentina.
u2
p=007 NI
15
p=396 NI.6
:
Nomenclature (A,b,c) = state space realization d = disturbance variable f = optimal gain vector G(s) = transfer function in s domain G ( t ) = transfer function in z domain H ( k ) = step function hi = weighting sequence coefficient k = discrete time N = model order NI = number of iterations of the algorithm in section 4, step
-1 0
O
O
200
160
240 L
€T;i
40
I 1-1.
-10
I
V
P = solution of the algebraic Riccati equation P(k) = solution of the difference Riccati equation R(z) = return difference transfer function in z domain r = reference variable T, = sampling time u = manipulated variable x = state vector y = controlled variable
-oulomallc tunlng
1)
---p=oo2
10
l
b
0
i
1
40
8
4
80
8
A
120
~
