Obviously, the above procedure can be applied to many chemical engineering problems. This approach is especially suited for on-line updating of parameters of a computer-controlled process. Another area of application is the estimation of chemical reaction rate constants from raw kinetic data with certain assumed or known mechanism of reaction. Since the estimation procedure is based on the quasilinearization technique, this approach is particularly suited for systems of differential equations with mixed boundary conditions. T h e above procedure can be generalized in various ways. I t can be used to estimate parameters or coefficients whose values are not constants but are functions of the independent variable, t. Instead of Equation 10, the following weighted criterion can be used
Nomenclature = integration constant = experimental value of x a t position t s = mean mass axial dispersion coefficient = Jacobi matrix, defined by Equation 13 = specific chemical reaction rate = total length of reactor = total number of experimental values of x =
= = =
= = = = = = =
dimensionless Peclet group, or unknown constant parameter defined by Equation 10 reaction rate group dimensionless reactor length, or independent variable final value o f t integration step size flow velocity of reaction mixture concentration of reactant A, or dependent variable concentration of reactant A before entering reactor experimental or measured value of x dependent variable
SUBSCRIPTS h = homogeneous solution n = nth iteration, assumed known 1 = (n 1)st iteration, assumed unknown n p = particular solution
+
+
literature Cited
-
and [xs bs]’ A[x, - b,] represents the quadratic form associated with the matrix A . The superscript T indicates the transpose of the vector within the bracket. T h e expansion of this quadratic form leads to a weighted sum of squares of the elements x - b, with the weighting determined by the elements of A . If A is a unit matrix, this quadratic form reduces to Equation 10. This generalized criterion has been used by Kalman (Kalman and Bucy, 1961).
Bellman, R., Kagiwada, H., Kalaba, R., “Quasilinearization, System Identification, and Prediction,” RAND Corp., Santa Monica, Calif., RM-3812-PR (August 1963). Bellman, R., Kalaba, R., “Quasilinearization and Nonlinear Boundary-Value Problems,” American Elsevier, New York, 1965. Danckwerts, P. V., Chem. Eng. Sct. 2, l(1953). Kalman, R . E., Bucy, R. S., J . Basic Eng. 83, 95 (1961). Lee, E. S., A.Z.Ch.E. J . 10, 309 (1964). Lee, E. S., Chem. Eng. Scz. 21, 183 (1966). Lee, E. S., “Quasilinearization and Invariant Imbedding,” Academic Press (in preparation), 1968. LVehner, J. F., FVilhelm, R. H., Chem. Eng.Scz. 6 , 89 (1956). RECEIVED for review December 5, 1966 ACCEPTEDAugust 8, 1967
LINEAR P R O G R A M M I N G 5-Transform Design of Digital
Controllersf o r Regulator Tvstems CARL I . H U B E R ’ A N D RICHARD I . KERMODEZ Department of Chemical Engineering, Carnegie-Adellon University, Pittsburgh, Pa. 152 13
Digital controllers for linear regulator control systems can readily be designed b y linear programming coupled with the Z-transform method. This method permits designs based on performance criteria and classes of deterministic input functions which are not suitable for conventional design techniques. The method is described in detail and illustrated by an example. N THE past 10 years a considerable amount of work has been l d o n e on the analysis and design of digital or sampled-data control systems. Articles by Koppel (1966a, 1966b), Min and Williams (1961), Mosler et al. (1966), and Williams and Min (1959) indicate that significant work is being done by chemical engineers in the design of sampled-data control sys-
1 Present address, E. I . du Pont de Nemours 8r Co., Wilmington, Del. 2 Present address, University of Kentucky, Lexington, Ky.
158
l&EC FUNDAMENTALS
tems. Most of the design methods for digital systems utilize the Z-transform method. This method can be used to reduce the linear differential equations describing a sampled-data system to linear algebraic equations. Kuo (1963), Ragazzini and Franklin (1958), and Tou (1959) have shown how to apply familiar design tools such as the Kyquist plot, Bode diagram, and root locus plot to the design of digital systems. Recent papers by Fegley (1964, 1965), Porcelli and Fegley (1964), Propoi (1963), Torng (1964), Whalen (1962), and Zadeh (1962) discuss various aspects of the linear programming
synthesis of optimal digital controllers for linear servo systems. Chemical engineers fr'equently encounter regulator control systems, and this paper describes the linear programmingZ-transform synthesis of optimal digital controllers for such systems.
fies some performance criteria for a specified disturbance or load function. T h e transfer function of the required digital controller is then obtained from Equation 5 as:
Typical Digitally Compensated Regulator System
Figure 1 is the block diagram of a typical digital control system, where all samplers are assumed to operate in synchronism with period T , gl(s) and g2(s) are the transfer functions of the plant, D ( z ) is the transfer function of the digital controller that is to be designed, and h(s) is the transfer function of a zero-order hold device. T h e Z-transform equations relating the system output and inputs are derived from Figure 1 ; they are Equations 1 and 2.
T h e form of Equation 1 is very similar to the usual closedloop transfer functions for continuous systems, and indicates that a n explicit transfer function for servo operation, c ( z ) / r ( z ) , can be obtained. An explicit transfer function for regulator operation, c ( z ) / L ( z ) , cannot be obtained as L(z)gl(z) = [Lgl](z). This situation is typical of digital or sampled-data systems.
Linear programming coupled with the 2-transform method offers a considerable increase in design flexibility over conventional design tools. This method permits the design of digital controllers for performance criteria and classes of input functions not suitable for conventional design methods. T h e designer can put explilcit constraints on the sampled system errors and in effect ''lailor" the system response. Physical constraints on the manipulated variable can readily be included in the design. Constraints can also be placed on the form of the digital controller to be synthesized. T h e method will be illustrated by its application to the design of the digital controller for the regulator system of Figure 1. For a constant set point Equations 1 and 2 in terms of the system error are as follows:
e(.) = - c ( z )
=
-m ( Z ) g ( z )
(3)
-
In any control system design the over-all response and the response of each element of the system must be physically realizable. This implies that no element of the system responds before it is disturbed. T h e following general forms of g(z), D ( z ) ,and K ( z ) satisfy this criterion:
The term 2-" in Equation 7 indicates a time delay of a sampling periods in the forward path of the system. Basic System Constraints
T h e 2-transform of a time function, f ( t ) , is defined by
Linear Programming Z-Transform Design Method
e(.) = - - c ( z ) = -d(z) K ( z )
Physical Realizability
(4)
where
F(z)
A
=
"
f(iT)zpi
(10)
i=O
Using Equation 10, Equations 3 and 4 can be written as follows:
+ . . + @(nT)z-, -($(OT) + $(lT)z-' + . . . + @ [ ( n- k) T]~-c"-')] (A1z-I + A~L-' + . . . + A ~ z - ~-] {e(OT)+ e(1T)z-1 + . . . + e ( n T ) z - n ) (11) +
@(OT) c$(lT)z-1
=
,
and
+
+ . . + $(IT)z-l = -(m(OT) + m(1T)z-1 + . . . + m [ ( I - a)T]z-('-")} (g(aT)z-rn + g [ ( a + l)T]z-("+l) + . . . + g ( l T ) z - J ] (e(OT)+ e(1T)z-1 + . . . + e ( I T ) z - ' ) (12)
@(OT) @(lT)2--1
,
where Equation 11 is written for n samples and Equation 12 for I samples. T h e pulse transfer function is given by
K(z)
= 1
+ AIZ-' + Azz-' + . . . + Akt-'
(1 3)
1
1
+ D(z)~hg1gzI(z) ( 5 )
K ( z ) will be referred to as the over-all pulse transfer function, while @ ( z ) is called the load function. I n general terms, the design procedure is to synthesize K ( z ) so that the system is stable and physically realizable, and satis-
Since constants A i are as yet unspecified, it can be assumed that K ( z ) is of finite length and order k . T h e length of K ( z ) determines the complexity of the controller transfer function D ( z ) , as shown by Equation 6. Expanding Equations 11 and 12 and equating coefficients of like powers of 2-l yield the following set of linear algebraic equations:
L I t).L(s )
1 Figure 1.
4%= -&A - e,
(14)
&
(15 )
and c($ c (s)
Typical digitally compensated regulator system
- (32- El where I$,, d ~E,, , E l , A,and -l? are column vectors: =
VOL. 7
NO. 1
FEBRUARY 1968
159
&
= [+(Ot), + ( I T ) ,
. . ., + ( i T ) ] - l ,
i = n or 1
81 = [e(OT),e ( l T ) , . . ., e ( j T ) ] - I , j
=
n or 1
(16)
Inclusion of Equation 1 9 in the design formulation ensures that the resulting digital controller, D ( t ) , will not cause the system to saturate. Additional constraint equations may be written to limit all or some of the sampled errors, e(iT),as follows:
where E , is the column vector of errors to be constrained, and 2,is the column vector of constraints. Thus, inclusion of Equation 2 0 in the design formulation allows the designer to tailor the system error response to some desired form.
Stability and Zero Steady-State Error Constraints
T h e necessary constraints to guarantee a stable, reliable digital controller and over-all system response are well known for servo operation. Stability and realizability constraints for regulator operation can be obtained in a similar manner and are stated below. A. K ( z ) must contain as its zeros all the unstable poles of the plant function g ( z ) = [hglgz](t). B. 1 - K ( z ) must contain as its zeros all the finite unstable zeros of the plant function, g ( z ) .
. . . @(OT)
0 0 0
d.T)
Neither A nor B implies that the system is being made stable by canceling out unstable poles or zeros. These two constraints ensure that the choice of K ( z ) be such that the system response is stable, and the controller transfer function contains no unstable poles or zeroes, thereby ensuring a stable and realizable controller. I t may be necessary to ensure that the sampled system error, e(iT), be zero a t steady state. T h e final value theorem for Z-transforms and Equation 3 indicate that for zero steadystate error the following must be true: If the load function, @ ( z ) ,contains ajth-order pole at z = 1, then K ( z ) must contain a zero a t z = 1 of o r d e r j or greater or
K(z)
(1
=
- 2-1)i
(21)
F(2)
where F ( z ) is a polynomial in 2-l. Thus the formulation of the stability, realizability, and zero steady-state error conditions reduces to the problem of constraining K ( z ) , or 1 - K ( z ) , to contain certain factors of the form (1 - u z - l ) . For each distinct factor, (1 - uiz-l), to be contained in K ( z ) the following equation can be written :
K ( z ) = 1 f A1z-I f A ~ z - ' f (1
Equations 1 4 and 15 represent the basic system constraint equations relating the variable vectors A, 8, and A?. The can readily be obtained from elements of the vectors, +(iT), 4 ( z ) once the disturbance, L, is specified. If 4 ( z ) is expanded into a power series in z-l, then, from Equation 10, the coefficient of 2 - j is $ ( j T ) . T h e elements of matrices 6 and 6 can be generated in the same manner. Control Effort and Error Constraints
In any real control system it is always possible that the control effort will saturate. T o avoid the nonlinear complication and to obtain a realistic design, inequality constraints are placed on the manipulated variable, m. Assuming symmetric bounds, these constraints are: 160
ILEC FUNDAMENTALS
-
UJ-')[~
+
b,Z-'
f biz-'
,
=
.. f
+ ...
f bk-iZ-('-')]
(22)
For the steady-state error condition 61 is 1. For the stability condition u i is the value of the pole of g ( z ) to be canceled. Expanding Equation 22 and equating the coefficients of like powers of z-l yields k equations in the A's and b's. Elimination of the b's gives a single linear algebraic equation in the A's, of the form : Alut'-'
+ A ~ u , ' - ~+
...
+ Ak-1~1 f
At =
-u.' I
(23)
This procedure is repeated for each additional distinct factor, yielding a set of constraints on the coefficients of K ( z ) . T h e same procedure is used to constrain 1 - K ( z ) to contain the necessary factors. If a total of j distinct factors are to be included in K ( z ) or 1 - K ( L ) ,the resultingj constraints can be written as
fi
where 8 is a constant column vector of length j and is a constant k X j matrix. For example, for the zero steady-state error condition ut is 1 and Equation 23 becomes:
... +A,=
Al+Az+
-1
(25)
and Equation 24 becomes:
-1
= [ l , 1, 1,
. . ., l ] [ A i ,
A l , . . ., A k ] . - '
(26)
Performance Criteria and Objective Function
Fegley (1964), Fegley and Hsu (1965), and Porcelli and Fegley (1964) have givm several useful performance criteria as particular applications of the general criterion of minimization of the integral of the weighted absolute value of a variablethat is, Minimize W ; JV =
1
w(t)
1 ~ ' ( tI )dt
(27)
where v ( t ) is the variable of interest and w(t) is the weighting function. I n discrete systems the integral is replaced by a summation and Equation 27 becomes
TV = 5 w i i T )
~ u ( t ~ = ) l
i;j~l
(28)
a=O
where & and p a r e row vectors. If w(iT) is replaced by iT and u(iT) is replaced by e(iT), then Equation 28 becoines the discrete version of the integral of the time-weighted absolute error (ITAE). If w(iT) is replaced by iT and u(iT) by m ( i T ) , then Equation 28 becomes a time-weighted control effort criterion. Clearly the weighting sequence need not be time, iT,but could be a cost sequence associated with product quality for the former case or the cost of control effort for the latter. I n some cases it may be desirable to formulate 1he objective function in terms of both system error and control effort, as follons:
,.
and .I?+*- contain the T h e vectors p+s-, A+,-, variables of the system; all other vectors and matrices are known if the disturbance is specified. Equation 31 is a more convenient form of the objective function and follows from a fundamental theorem of linear programming. Conversion of the inequalities to equalities and application of a linear programming algorithm to Equations 31 and 32 completes the design method. A typical linear programming solution of these equations yields the sampled error vecttrs, &+ and &-, the sampled control effort vectors, &If+ and -11-, and the coefficient vectors, A^+ and A-, which minimize W. T h e transfer function of the required digital controller is obtained from Equation 6. load Function
T h e necessary requirements on the input load function, q ( z ) , for the above method are: T h e sampled load function must be deterministic and mathematically representable as a power series in 2-1. T h e sampling frequency of the digital controller used must be sufficiently high so that the Z-transform or power series representation of the load function is a good representation of the actual continuous function. This condition will usually be satisfied if the sampling frequency is a t least twice the highest frequency of the continuous function. If the load change L ( t ) is known, $ ( z ) can usually be obtained from tables of Z-transform. T h e power series form of $(z) is obtained from the Z-transform by dividing the denominator into the numerator. Where a tabular L ( t ) is available, the coefficients of d(z) must be obtained by either approximating the data with a function, or numerically solving the differential equation corresponding to $(s) = L ( s ) g(s). This approximation procedure leads to a linear prograniming design which is also approximate. This is not the case for servo systems, where exact design is possible even when the set point change is known only in tabular form. Numerical Example
linear Programming Formulation
Equations 14, 15, 19, 20, 24, and 28 describe the control system and objective function, and provide constraints on the error response, control (effort, and form of K ( z ) . These equations are in a form suitable for linear programming. I n all linear programming alizorithms, the variables must be nonnegative. Thus, all variables are replaced by the difference between two nonnegative variables as presented by Hadley (1962), as follows:
Q
+
Q+
- Q-
(30)
where Q+ and Q - are nonnegative variables and Q is the variable unrestricted in sign. Using the above procedure the objective function and system constraints can be compactly summarized in vector-matrix form as follows:
T h e system shown in Figure 2 is a master-slave regulator system for an idealized model of a continuous stirred-tank chemical reactor. The master, or outer loop, operates on the measured concentration error, e l , and drives the set point of the slave loop. The slave loop operates on the measured reactor temperature error, e 2 , and changes the manipulated variable, coolant flow rate, q c . D ( z ) is to be designed to regulate the concentration of component A,%, in the face of a disturbance in the feed concentration of A,xi. T h e slave loop controller, g,, is a continuous proportional controller, and was independently tuned to provide a reasonably fast response to set point changes. T h e mathematical model is the same as that used by Kermode and Stevens (1961) and includes the following assumptions :
-
The reactant mass is perfectly mixed and of constant volume. T h e reaction is first-order, exothermic, and irreversible, A products, and the reaction rate constant is given by the Arrhenius equation. xi
Figure 2.
xi
Master-slave regulator control system VOL. 7
NO.
1
FEBRUARY 1 9 6 8
161
T h e energy balance is based on an arithmetic mean temperature gradient between the coolant and reactant mass. All flow rates except coolant rate are constant. T h e transfer functions for this system are x
($1
qc(S)
.
0.0340
gl(s) = __ = (3
+ (s + b ) ' a)
X(S)
g/l(s) = __ = 4s)
0.1800(~- 0.3829)
(f
+ + 6) a)(s
the coolant rate is sufficiently attenuated so that explicit constraints are not needed. The steady-state sampled error must be zero. For this 2-1) in its denominator. example + ( z ) contains the factor (1 The sum of the time weighted absolute error is to be minimized. The disturbance is a 5% step increase in the inlet composition, x i , a t zero time.
-
T h e system constraints and objective function for the design are of the same form as Equations 31 and 32.
Table I. Pulse Transfer Function and Corresponding Digital Controller CASE1 k = 7, n = 18, T = 0.1 min.
K(z)
where a = -0.2111; b = 0.1596; gl,g2, g21,g12,and gi are reactor transfer functions, g , is the transfer function of the slave controller, and h is the transfer function of the hold. The basic constraints on the system of Figure 2 are the same as Equations 3 and 4 where
1
- 0.11172-' - 0.2241~-' - 0.23842-3 0.2537~ ~0 ~. 2 6 0 5 ~ -~ 0 . 0 5 5 4 ~ -+ ~ 0.14332-7 -4939.97 - 780.662-' 4- 3638.622-2 -
+
+
+
2 0 5 . 0 3 ~ - ~ 2 0 7 . 1 5 ~ - ~ 9293.472-6 927.442-6 1 3 9 6 1 , 1 6 0 5508.862-8 D(z) = 1 0.84172-' - O.3301~-*- 0.4520~-' 0 . 4 8 0 9 ~ -~ 0 . 5 0 2 2 ~ -~ 0.30362-6 0.09052-7 0.13652-8 ITAE = 0.00690
+
+
-
+
+
CASE2 k = 5 , n = 18, T = 0 . 2 m i n . K ( z ) = 1 - 0.21 852-l - O.4440~-' - 0.48442-3 0.110i2-4 0.2570~-6 -2547.25 - 7 8 9 . 9 8 ~ - ' $ 1 3 6 2 . 2 7 ~ - ~ + 4 5 6 9 . 9 1 2 - ~ + 9 7 9 . 0 4 ~ --~ 6120.952-5 2242.480 D(r) = 1 0,68952-1 - 0.64232-2 - 0.88752-8 0.23332-6 0.5499t-'+ 0.15702-6 ITAE = 0.0081221
+
+ +
+
CASE3 k = 5, n = 18, T = 0 . 4 min. K ( z ) = 1 - 0.42392-' - 0.8112t-*
For this system there are two different load functions, + ( z ) and + f i ( z ) , rather than one as in the previous example. This situation causes no additional problems, since both functions are known once the disturbance in x i is specified. T h e design of D(z) is based on the following specifications: T h e magnitudes of the first 1 control efforts, q,, are constrained, so that qc 0.18. I t is assumed that after I samples,
I I
