1658
I n d . Eng. Chem. Res. 1988,27, 1658-1664
Ku, H. M.; Karimi, I. A. "Scheduling in Multistage Serial Batch Processes with Finite Intermediate Storage, Part 1: MILP Formulation" Presented a t the AIChE Annual Meeting, - Miami, 1986a, Paper 72e. Ku. H. M.: Karimi. I. A. "Scheduling in Multistage Serial Batch Processes with Finite Intermediate storage, Part 12: Approximate Algorithms" Presented at the AIChE Annual Meeting, Miami, 198613, Paper 72e. Ku, H. M.; Rajagopalan, D.; Karimi, I. A. "Scheduling in Batch Processes" Chem. Eng. Prog. 1987,83(8), 35. Kuriyan, K.; Reklaitis, G . V. "Approximate Scheduling Algorithms for Network Flowshops" Znd. Chem. Eng., Symp. Ser. 1985,92, 79. Parakrama, R. "Improving Batch Chemical Processes" Chem. Eng. 1985,Sept, 24.
Rajagopalan, D.; Karimi, I. A. "Scheduling in Serial Mixed-Storage Multiproduct Processes with Transfer and Set-up Times" Proceedings of the 1987 Foundation of Computer-Aided Operations Conference, Park City, UT, July 5-10, 1987a. Rajagopalan, D.; Karimi, I. A. "Completion Times in Serial MixedStorage Multiproduct Processes with Transfer and Set-up Times". Comp. Chem. Eng. 1987b,in press. Rajagopalan, D.; Karimi, I. A. "Scheduling in Serial Mixed-Storage Multiproduct Processes with Transfer and Set-up Times" Technical Report 8703,1987~; Center for Manufacturing Engineering, Northwestern University, Evanston, IL.
Received for review September 14, 1987 Revised manuscript received March 30, 1988 Accepted April 22, 1988
An Adaptive Estimation Algorithm for Inferential Control Mohammad T. Guilandoust, A. Julian Morris, and Ming T. Tham* Department of Chemical and Process Engineering, University of Newcastle upon Tyne, Newcastle upon Tyne NE1 7RU, U.K.
An estimator relationship for inferring infrequently measured process outputs, from other more rapidly sampled secondary outputs, is derived. The structure of the estimator is not application dependent, and its parameters can be continuously estimated and updated. As a result, slow variations in plant or disturbance characteristics can be tracked. Closed-loop control schemes using the estimated output data are shown to exhibit superior performance over those using the slowly measured output data. In contrast to some other inferential control strategies suggested in the literature, the method proposed requires minimal effort for estimator design and secondary output selection. In digital control systems, the infrequent measurement of some process outputs, determined by sampling limitations, prevents the early detection of load disturbances. This results in large deviations from the desired output and consequently long disturbance recovery times. Often, these adverse effects cannot be acceptably overcome by the use of existing advanced control algorithms and can lead to unsatisfactory control system performance. Examples of industrial situations where this can occur are in product composition control of distillation columns (Patke et al., 1982) and chemical reactors (Wright et al., 1977). In these cases, the sampling delay is a direct result of the long cycle time of on-line composition analyzers. Because of potential problems due to infrequent sampling, the control of product quality of many industrial multicomponent columns is commonly achieved by maintaining an a priori chosen tray temperature near to its set-point value. However, this type of single-temperature feedback control is not always effective since maintaining a constant tray temperature does not necessarily result in constant product composition (Patke et al., 1982). The problem of controlling infrequently measured process outputs has long been studied, and publications in this area date back to the early 1970s. There have been essentially two approaches to the problem. One is to design special controllers for the infrequently sampled outputs. For example, Soderstrom (1980) formulated a number of minimum variance controllers enabling the manipulated (control) input to be changed between the sampling intervals of the primary process output. These control algorithms were only developed for first-order plant models, and no comparative results were presented. Parrish and Brosilow (1985) also proposed a controller design method based upon the philosophy of reconstructing the effects of disturbance inputs. The controller parameters are de-
termined on-line by heuristic tuning rules. Their simulation results indicated superior performance of their control strategy over that achieved by using conventional PID control. A second approach to deal with the problem of controlling infrequently sampled process outputs is to use the information provided by other more easily measurable variables. For example, this information can be used to provide an estimate of the controlled output. The estimated values of the output can then be used for overall control of the plant. Control schemes based on the feedback of estimated outputs are often termed "inferential control schemes". An ideal situation arises when the plant states are completely observable from the secondary outputs. Under such circumstances, Kalman filtering techniques can be employed to estimate plant states using the secondary output measurements. Estimates of the controlled output can then be computed by using its relationship with the states. Control of the plant is achieved by feedback of either the state estimates or the output estimates to appropriate controllers. Published literature on the above methods is extensive. The papers by Morari and Stephanopoulos (1980a,b) and Morari and Fung (1982) are amongst the more recent contributions describing the application issues of these techniques. However, the use of Kalman filters is confined to situations where the plant is completely observable from the secondary outputs. For most plants, such a set of secondary outputs can be difficult to determine or, in some cases, may not even exist. Brosilow and co-workers (Cwiklinski and Brosilow, 1977; Joseph and Brosilow, 1978) have suggested an estimator design technique using an input-output representation of the plant. The design is approached by obtaining a least-squares-based static estimator, which can be used to infer the controlled output from secondary measurements
0888-5885/88/2627-1658$01.50/0 0 1988 American Chemical Society
Ind. Eng. Chem. Res., Vol. 27, No. 9, 1988 1659 at steady state. The estimator is then additionally made applicable to the transient period by incorporating heuristically derived lead-lag elements into its structure. Methods for minimizing the steady-state estimation error by appropriate choice of the secondary measurements were also proposed. However, a set of secondary measurements for which this error is satisfactorily small may not exist. In fact, the evaluation work carried out by Patke et al. (1982) indicated that an inferential control scheme based on the output estimation technique proposed by Brosilow and co-workers can result in significant offsets. Moreover, to apply the technique, it is necessary to know the gains and approximate time constants of the controlled output and all secondary outputs, for all plant disturbances and manipulated inputs. The control strategies described so far rely either on infrequent measurements of the controlled output or on the use of secondary measurements. It is, however, possible to use both types of measurements. One technique is to set up a parallel cascade control strategy (Luyben, 1973). Here measurements of the controlled output are fed to a controller whose output acts as the set point to a secondary output controller. Although no offset problems arise with this control configuration, the transient performance can at times be poor (Patke et al., 1982). A second technique for using both measurements of the controlled output and secondary output is within an adaptive inferential control framework. In such a scheme, the infrequent measurements of the controlled output are used only during parameter estimation while plant control is achieved using a secondary output at its faster sampling rate. To our knowledge, very little investigative work has been carried out in this area, although D'Hulster and van Cauwenberghe (1981) have suggested an algorithm using a similar approach. Their algorithm is, however, restricted to the assumption of a first-order plant model, an equal number of disturbance inputs and secondary measurements, and the use of a dead-beat control law. In the present paper, an adaptive algorithm for estimating plant-controlled output from a secondary measurement is developed. The algorithm is derived from an unknown general input-input representation of the plant. Attention has been confined to the SISO case, with extensions to MIMO situations being the subject of future work. A second estimation technique has also been derived from a state space representation of the plant. This approach is described elsewhere (Guilandoust et al., 1987) and has been briefly compared with the input-output method to be described here (Guilandoust and Morris, 1985).
Plant Description Let the plant dynamics be represented by the following discrete time relationships:
measurable at each time step. ml and m2 are the time delays in the responses of the controlled and secondary outputs, respectively, to changes in the manipulated input, u ( t ) . w ( t ) E R' is the vector of stationary random and unmeasurable load disturbances. z-l is the backward shift operator. Gl(z-') and G2(z-') and all elements of row vectors L l(z-l) and L 2(z-1) are polynomial ratios with the order of the numerator less than or equal to that of the denominator. It is assumed that a change in any one of the components of w ( t )affects both y o ( t ) and uo(t); otherwise, u o ( t ) will not be suitable for estimation of y o ( t ) . Consequently, an entry of L l(z-l) cannot be zero, unless the same entry of L2(z-')is also zero and vice versa.
Adaptive Estimation Algorithm In this section, a suitable structure for an output estimator is obtained in a systematic manner. Based on this estimator structure, an adaptive algorithm for estimating the controlled output at the rate the secondary output is sampled is derived. Equations 1and 2 can be written in the following form: y ( t + d ) = Gl(z-l)u(t - ml)
+ [ ( t )+ 4 t + d )
+ {(t)
u ( t ) = Gz(z-')u(t - m2)
(3)
where
+
~ ( t d ) A zdv1(t) - v2(t)
(44
5 ( t ) = L,(z-')w(t)+ v 2 ( t )
(4b)
{(t)= Lz(z-l)w(t)+ v d t )
(4c)
and The physically nonexisting sequence ~ ( tmay ) be noted to be zero mean and white. This is because zdvl(t)is a zero mean, white sequence and so is its algebraic sum with v&). Since v 2 ( t ) is a zero mean, white sequence and w ( t ) is a row vector of stationary signals, eq 4b and 4c may according to the spectral factorization theorem (Astrom, 1970) be represented as
[ ( t )= H,(z-')o(t) { ( t )= H2(z-l)&)
(5)
where H1(z-') and If2(%-')are stable and proper polynomial ratios and the sequence o ( t ) is zero mean and white. Notice that both [ ( t )and { ( t )are generated by the same white noise sequence w ( t ) . If the driving noise for one of the signals was different from the other, then the two signals would clearly be uncorrelated. This cannot be the case, because variations in load disturbances are assumed to affect both controlled and secondary outputs (see eq 3 and 4). Moreover, neither Hl(z-') nor H,(z-') can be zero without the other being zero too. Substituting from eq 5 into eq 3 and eliminating o(t)results in the relationship ~ (+td ) = [ G ~ ( z - ~ ) z - ~-~ G + ~" ( Z - ~ ) H ~ ( Z - ' ) / H ~X( Z - ' ) ] u ( t - mz) + [H1(z-')/Hz(z-')]u(t)+ 4t + d ) (6) for ml I m2. In the case where ml Im2, eq 6 can be rewritten as
where y o ( t ) and u o ( t ) are the controlled output and secondary output at time t. y ( t ) and u ( t ) are their corresponding observations contaminated with measurement noises vl(t) and vz(t), respectively, assumed to be zero mean white sequences. d is the analyzer delay associated with y ( t ) and is assumed to be an integer multiple of the discretization time step. Therefore, values of y ( t )only become available every d time steps. u ( t ) is assumed to be
~ (+ td ) = [ G ~ ( z - '-) Z-"~+"~G,(Z-')H~(Z-~)/H~(Z-') X u ( t - ml) + [Hl(z-')/H2(z-')]u(t)+ 4t + d ) (7) Equations 6 and 7 relate measurements of the controlled output, y ( t ) ,to those of the secondary output, u ( t ) , and the values of the manipulated input, u ( t ) . If G1(z-l), Gz(z-l), H1(z-l),and H2(z-') were all known, u ( t ) and u ( t ) could be used in eq 6 or 7 to compute
1660 Ind. Eng. Chem. Res., Vol. 27, No. 9, 1988
jqt + d ) = y ( t + d ) - t ( t + d ) (8) which is an estimate of the controlled output, yo(t). The variance of the deviations of 9(t + d ) from the controlled output is, by relationship 4a, limited to the sum of the variances of vl(t) and u2(t). The variance of the latter noise can, by appropriate choice of the secondary variable, be kept to a minimum. Thus, given a sufficiently accurate plant model, eq 6 or 7 may be used for recursive computation of delay free estimates of the controlled output at the rate u ( t ) is sampled. The accuracy of the estimates will be close to that of the measurements of the controlled output, y ( t ) . In many instances, due to the presence of a large number of disturbance inputs and complicated plant dynamics, the polynomial ratios G1(z-l), G2(z-l), and H2(z-') can only be approximately obtained using off-line system identification. In this paper, however, these polynomial ratios are not assumed to be known a priori. Instead, it is proposed that eq 6 or 7 be approximated on-line as described below. Based on the definitions of G1(z-l), G2(z-l),H1(z-l),and HZ(z-'), it would seem reasonable to choose an approximation of the following structure for eq 6 or 7 to serve as an estimator for computing the controlled output: B(z-') C(z-1) y ( t + d ) = -u ( t - m ) + -u ( t ) + 4t + d ) (9) A (2-l)
A(2-l)
where
m = min (m1,m2) A(z-') = 1 + ulz-' B(2-I)
= blz-l
C(z-1) = co
+ u 2 f 2+ ... + u g n + b 2 f 2 + ... + b , , ~ - ~
(10) The value of the integer n is chosen to achieve the required accuracy of approximation. It should be pointed out that eq 9 is not of the ARMAX type, although its structure may at first sight suggest otherwise. This is because u ( t ) is not an unknown random signal as in the case of an ARMAX relationship. Instead it is the measurements of the secondary output which, like y(t), respond to both the control input, u ( t ) ,and the unknown random disturbances. The parameters of eq 9, in its present form, cannot be straightforwardly estimated by using recursive techniques because y ( t ) is only measurable every d time steps. This problem may, however, be resolved. By use of the relationship 8, eq 9 can be written in the form (11)
Consider initially the case where n = 2. In this case, A@') in eq 11 may be factorized as A(2-l) = 1 ~ 1 z - l u2Y2 = (1 + 71.Y1)(1 72Z-l) (12)
+
+
+
For generality, assume that A(2-l) has complex roots; i.e., T~ and T~ are a complex conjugate pair. Now let both sides of eq 11 be multiplied by polynomials
Pl(Z-1)
= [l -
T12-1
+ ... + (-Tl)d-lZ-d+l 1
(13)
and
Pz(z-') = [l - 722-1
+ ... +
1
(14) It may readily be noted that the resulting relationship may be written in the following form: [1 ( - ~ l ) ~ z - ~+] [( 1- ~ 2 ) ~ 2 - ~ ] 54-( td ) = Pl(z-')P2(z-1)[(blz-1 + b2z-2)u(t- m ) + (co + clz-l + ~ ~ z - ~ ) u ((15) t)] (-Tz)d-1Z-d+1
n(z-d)= 1
+
+ (Y2dZ-2d + ( - - T ~ )and ~ azd=
(16)
with red coefficients a d = ( - T ~ ) ~ ( 7 ' ~ ~ ) ~ . Since both A(2-l) and a(z4) are polynomials with real is also coefficients, it follows that the product Pl(z-1)P2(z-1) a polynomial of order 2(d - 1)with real coefficients. As a result, eq 15 can be rewritten as (1
+
+ (uzdz-2d)9(t + d ) = (p1Z-l + pzz-' + ... + - m) + (yo + ylz-l + y g 2 + ... +
@2dZ-2d)U(t
y2dz-2d)u(t) ( 17)
The dynamics of eq 11 do not change as a result of the above multiplication, and eq 11 and 17 are equivalent. However, it may be observed that, in contrast to eq 11, the parameters of eq 17 can be recursively estimated every d time steps when a measurement of y ( t ) is available. The parameter estimation steps will be further discussed at a later stage. The above reformulation of eq 11 for values of n other than 2 is now readily deduced. All that is required is to note that the polynomial A(2-l) in eq 11can always be factorized as the product of first-order polynomials. For each first-order polynomial, both sides of eq 11should be multiplied by a polynomial similar to those in definitions 13 and 14. I t may thus be observed that, for a given value of n, eq 11 can be reformulated as
(1 + (plz-l
+ c1z-1 + c2z-2 + ... + c,z-n
A(z-l)g(t+ d ) = B(z-')u(t - m) + C(z-')u(t)
As the dth power of the complex conjugate pair, T~ and r2 is another complex conjugate pair, the product of the square bracketed terms on the left-hand side of eq 15 is a polynomial with real coefficients. It is, therefore, of the form
+ + andz-"d)9(t + d ) = + p2z-2 + ... + PndZ-"d)u(t - m) + (yo+ ylz-' + yzz-' + + yndz-nd)u(t)(18) ..e
Since eq 18 and 11 are equivalent, both constitute an estimator for recursive computation of the controlled output. However, eq 18 has the desired form because its parameters, unlike those of eq 11, may be recursively estimated. In order to estimate the parameters of eq 18, it is rewritten as y ( t ) = -a&(t - d ) - ... - nn&(t - nd) + &u(t - m - d - 1) + ... + &&[t - m - (n + l)d] + + y n d U [ t - (n + 1)d] + t(t) (19) you(t - d ) + This can then be expressed in a more compact form y ( t ) = eTV(t- d ) + e ( t )
(20)
where
eT = [ a d , ..*, a n d , p1,
**a,
Pnd, 70,
***,
cp(t - d ) = [-jl(t - d ) , ..., u ( t - m - d -l),
Yndl
..., u ( t - d ) , ...I
and e ( t ) is the equation error. Now let a value of y ( t ) be available at time t. This new value can @eused to update the current estimate of the parameters, 8. To do this, we first update the elements of the data vector, rp(t - d), using saved values oj9, u , and u. After updating the estimate of parameters 9, we can also compute the a posteriori value of g ( t ) = OTv(t- d ) for updating cp(t - d ) at time t + d. During the period t to t + d, the estimated parameters can be used in the estimator relationship 18 for the recursive computation of jl(t + d ) to g ( t + 2 4 . These are estimates of the actual values of the controlled output, yo(t),during this period (see eq 2 ) . Equation 20 has n(2d + 1) + 1 unknown parameters. Our experience indicates that a value of n between 1 and 3 is normally sufficient for good estimation. d may typically lie in the range 2-6. A first-order estimator ( n = 1)
Ind. Eng. Chem. Res., Vol. 27, No. 9, 1988 1661 COOLING WATER
Disturbance Inputs (w)
I
Parameter
9
out u t Esgmator output Ea timator
Figure 1. Adaptive inferential feedback control scheme with constant parameter controller.
may therefore have between 6 and 14 unknown parameters and a second-order estimator between 11 and 27. These figures may at first sight give the impression that the tuning in period for larger values of d will be inconveniently long. This is not, however, the case. The reason is that, although the number of unknown parameters increases with d, so does the number of data (intermediate values of both u and u ) supplied to the estimation algorithm at each sample time of the controlled output (see eq 20). The results of simulation studies, some of which are presented later, have confirmed that the initial tuning in period increases only slightly for larger values of d (analyzer delays of up to six sample intervals have been looked at). Estimates of the controlled output can often be obtained with sufficient accuracy by using a first-order estimator. In this case, the parameter a d of the first-order estimator can easily be seen to be related to ul as ad
=
(-al)d
As may be noted from eq 6, lull is normally less than unity (except for the special case where H,(z-') has a dominant positive zero). Therefore, for larger values of d , it follows that a d
