Copyright © IFAC Identification and System Parameter Estimation 1982 , Washington D.C., USA 1982
VALIDATION OF PREDICTION ERROR IDENTIFICATION RESULTS BY MUL TIVARIABLE CONTROL IMPLEMENTATION G. A. Van Zee* and O. H. Bosgra** * Wageningen University of Agn'culture, Department of Physics and Meteorology, Duivendaal1, 6701 AP Wageningen, The Netherlands * *Delft University of Technology, Department of Mechanical Engineen'ng, Mekelweg 2, 2628 CD Delft, The Netherlands Abstract. This paper discusses the validation of models obtained by prediction error identification of an 8-input, 8-output pilot scale heat exchanging process. Approximating models, obtained from closed-loop experiments, have been used to design multivariable feedback control laws for the process. The in-line digital implementations of these control laws on the actual process have been evaluated by comparing the actual closed-loop behaviour with model predictions. The results have been used as a measure for model accuracy and demonstrate certain deficiencies of the usual onestep ahead prediction error identification criterion. A modified criterion giving more weight to low frequency model errors is shown to yield models with significantly improved actual closed-loop behaviour. Keywords. control.
Model validation, Prediction error identification, Multivariable
1. INTRODUCTION The accuracy of system identification methods has been studied by statistical analysis of the properties of the obtained estimates (Gustavsson, 1972; Mehra, 1974; Soderstrom, Ljung, Gustavsson, 1975) and has been evaluated in case studies (Gustavsson, 1969; Eklund, Gustavsson, 1973; Olsson, 1973). However, there appears to be little knowledge about the relation between the accuracy which is pursued in system identification methods and the suitability of the resulting models for the design of control systems. In this paper we will evaluate this relationship on the basis of identificat ion and multivariable control experiments on a complex pilot scale heat exchanging process operating under inline di g ital control. The process has some properties which are characteristic for industrial processes, such as the occurrence of spatially distributed heat transfer phenomena, the interaction between the process variables, and the presence of stochastic disturbances. Due to its complexity, the process can only be tractably described by approximate models. The paper discusses the validation of models which have been obtained by prediction error identification, a method which has received much attention in the recent literature because of its expected suitability to compute approximate models (Ljung, 1978; Caines, 1978; Anderson, Moore, Hawkes, 1978; Ljung, Van Overbeek, 1978). By designing linear quadratic optimal control laws for the process using identified process' models, and by the implementation of these control laws on the actual process, we are able to compare the actual closed-loop behaviour with model pre643
dictions. The results will serve as an ultimate means of model validation leading to conclusions with regard to the identification method used. The paper is organized as follows. In section 2 prediction error identification is discussed and applied to the pilot scale process which is described in the appendix. In section 3 the resulting models are validated by the design and implementation of multivariable control systems. The validation results suggest the use of a modified prediction error identification criterion which is described and evaluated in section 4. 2. PREDICTION ERROR IDENTIFICATION The principle of prediction error identification is to minimize a criterion of the form N-1 J =
L
eT(k)Qe(k)
(2.1)
k=O
T
.
where N is the number of samples, Q=Q >0 is a weighting matrix, and e(k) is the prediction error (2.2)
e(k) = y(k) - y(k)
In (2.2) y(k) is the true system output and y(k) is a one-step ahead prediction for y{k) which follows by the assumed model. If the process model can be written as x(k+1) = Ax(k) + Bu(k) + Kw(k) , x(O) y(k) Cx(k) + w(k)
x °(2.3)
644
G. A. Van Zee and O. H. Bosgra
with G(z)=C(zI-A)-l K+I strictly minimum phase, and with w(k) a wide sense stationary zero mean white noise with covariance W>O, then a stationary linear minimum variance estimate y(k) follows by (Sage, Melsa, 1971) x(k+1) = (A-KC)x(k) + Bu(k) + Ky(k) , x(O) = Xo y(k) = Cx(k) (2.4) The properties of prediction error estimation have been extensively studied (Ljung, 1976; Caines, 1976); Goodwin, Payne, 1977). In Soderstrom, Ljung, Gustavsson (1976) the identifiability of the parameter vector defining the true system, is discussed for the case of closed-loop experimental conditions. In the associated consistency analyses for the parameter estimates it is assumed that the considered set P of parameter vectors contains the true parameter vector. It also has been shown that if this assumption is not valid, then there still exists a converging sequence of estimates, yielding an approximate model which is at a minimal distance from the exact model in the sense of the identification criterion being used (Ljung, 1978 ; Caines, 1978; Anderson, Moore, Hawkes, 1978). A problem which occurs in the application of prediction error identification is formed by the computational costs of numerical parameter optimization. In Van Zee (1981) a method is derived to compute the gradient g of the criterion J with respect to the parameters {A,B,K,C,xO} in (2.3), using a dynamic programming formulation. In this approach the computational costs to compute g are the same as the costs to compute J, and hence are gene rally significantly lower than the costs of the commonly used sensitivity functions approach (Gupta, Mehra, 1974). If in addition a variable metric method is used to obtain an approximate Hessian, then the calculation of the sensitivity functions is completely circumvented, and a computationally optimal overall procedure ~s obtained. However, even if an efficient numerical procedure is available, then for the type of processes which we consider in this paper the computational costs of parameter optimization are still intractably high compared to the costs of explicit identification methods. A drawback of explicit methods is that these ge nerally are less accurate and hence a useful compromise is to use an explicit method to obtain initial estimates for a parameter optimization method (Caines, Rissanen, 1974; Tse, Weinert, 1975; Akaike, 1976; Van Zee, 1981). In our applications we have computed initial process models in state-space form by approximate stochastic realization of the joint input-output signals, which involves the consecutive steps of spectral factorization, decomposition into open-loop models, and realization (Van Zee, 1981). The joint input-output signals of the pilot scale process have been measured under the closed-loop experimental conditions which are described in the appendix. The model order has been chosen in the realization step by an analysis of singu-
lar values. We have used a canonical parametrization which is based on a nice selec~ion (Denham, 1974) of the rows of the observability matrix of the initial model. A procedure to determine the computationally most favourable nice selection is described in Van Zee (1981). For the considered model, with the dimensions as given in the appendix, and with order n=10, the resulting number of parameters is 250, including the initial state xo. The number of samples used for the parameter optimization is 1024, where each sample consists of an input and output vector. For these dimensions one complete iteration of the optimization procedure, involving an average of 2 evaluations of J and g, appears to cost approximately 36 seconds of CPU time (IBM 370/158, PL / I). Figure 2 .1 shows the path of J as a function of the iteration number i, where i=O represents the initial model. The parameter optimization has been stopped before convergence to a minimum point has been obtained, because our purpose is mainly to judge the improvement of the model, by comparing the suitability of the initial model MO and the final model M for the design of multivariable control systems. 3. VALIDATION BY MULTIVARIABLE CONTROL IMPLEMENTATION In the practice of process control, it is often impossible to maintain a desired operating point by purely proportional feedback control, due to the adverse effects of constant disturbances which act upon the system. This problem can be solved by the inclusion of integrating action in the control system, an approach which also solves the problem to track a piecewise constant setpoint. For the design of multivariabl e control systems we take the linear quadratic optimal estimation and control approach, because this approach makes a direct and expli ci t use of the deterministic and stochastic part of the models. Assuming that the state x(k) of a process model is available, we consider the linear optimal regulator problem to determi ne Lp, LI in u
c
(k) = [-L
P
I
-L
I
0.1)
1 [X(k)] a (k)
such that the criterion
L
J = {uT(k)Q u (k) + a k=O c u c 0.2)
~s
r
minimized, subject to the model equations
X(k+1)] La (k+1)
=
[A 0J[X(k)] C I a (k) s
B ]u (k), rX(O)J La( O) [ OC c
+
0.3)
Validation of Prediction Error Identification Results
Ol [x(k)l IJ a(k)j
645
pose to obtain a better correlation of the criterion, with the accuracy of the control performance predictions by the models.
If the system (3.3) is controllable and observable, then there exists an asymptotically stable control law (3.1) which can be obtained via the stabilizing solution of the associated algebraic Riccati equation (Kwakernaak, Sivan, 1972). The index c in ucand Bc is used to distinguish the control input Uc from the input u=[uE : ua)T of the identified models. The weighting matrices, which we restrict to be diagonal with positive entries, are the design parameters of the control system. The state x(k) of the model can be estimated from the observed variables u(k), y(k) by the minimum variance estimator (2.4), which is directly obtained from the model (2.3). The index s in Ys and Cs is used to define the part of the output which should track a piecewise constant setpoint s(k), in spite of the occurring disturbances. Here we discuss only the simplest of the studied control problem~, being the problem with T
y =[F , F , L , Ls ' Ls) , s P s P 1 2 using the symbols as detined ~n the appendix. As a reference case for this control problem, we use the performance of the experimentally tuned single-loop PI controllers which are described in the appendix. Due to the limited space we will only show the setpoint step responses in the figures, as e.g. in figure 3.1 for the single-loop case. In the multivariable control system design we have chosen Qu' Qy, Qa such that the simulated closed-loop behaviour improves the singleloop results. Figure 3.2 shows the simulation results for the initial model MO. However, both the model MO and M appeared to be insufficiently accurate, i.e. the implemented control systems appeared to be unstable. To obtain control systems which are more robust for the present model errors, the desired quality of control had to be reduced. Stable control system implementations could only be obtained by decreasing Qo , i.e. by reducin g the integrating action. This leads to the conjecture that the instability of the implemented control systems is mainly due to low frequency model errors. The predicted and the actual performances of the control systems as obtained for the models MO and M by reducing Qo to .01 Cl , are shown in figures 3.3 - 3.6. A rem~rkable conclusion from these results is that the accuracy of the control performance prediction by the final model M of. the parameter optimization is worse than in the case of the initial model MO. This is apparently in conflict with the accuracy measure on which the prediction error identification method is based. A possible explanation of this conflict is that the minimization of onestep ahead prediction errors improves mainly the model accuracy in the high frequency range, whereas for the solution of the stated control problem the model accuracy .in the low frequency range is more important. Therefore we propose in the following a modified prediction error identification criterion, with the pur-
4. A MODIFIED PREDICTION ERROR IDENTIFICATION CRITERION Consider the modified criterion function
J
N-1 = [eT(k)Qe(k) + ET(k)Q E(k») E k=O E
I
(4. 1 )
with E(k) defined on the basis of e(k) in (2.2) by dk+1) = dk) + e(k),
dO) = 0
(4.2)
The role of E(k) is to enhance the contribution of the low frequencies in the residuals e(k) to the criterion function. The possibility to use a criterion of the proposed form has also been mentioned by Astrom (1980). It is shown in Van Zee (1981) that with some minor modifications, the criterion (4.1) can be minimized by the same parameter optimization procedure as the criterion (2.1), at only slightly higher computational costs. Figure 4.1 shows the path of J E as a function of the iteration number i, using the same initial model as before. The difference between the final models M and ME resulting from the minimization of J and J E respectively, appears from a comparison of the residuals. In approximate modelling the residuals are useful in order to investi gate the frequency region in which the largest model errors occur. Figure 4.2 shows the diagonal of a smoothed estimate of the spectral matrix of the residuals for both of the final models. The conclusion from these results is that the minimization of J yields smaller model errors in the low fre~uencies, whereas the minimization of J yields smaller errors in the high frequencies. For the model ME the originally stated control problem ha s been satisfactorily solved, which is demonstated by figures 4.3 and 4.4 showing the predicted and actual performance of the control system. Similar results have been obtained for other control problems, e.g. the control of y =[F , L , To ' To Tp)' s s P 1 2' 2 and for other datasets (Van Zee, 1981). Hence we conclude that the modified criterion function J correlates better than the original criterion J with the accuracy of the control performance predictions by the models. In principle there are some other possible solutions to handle the low frequency model errors in prediction error identification. One could for example enhance the low frequency model errors by replacing the one-step ahead prediction errors in the criterion function (2.1) by multi-step ahead prediction errors. Another possible approach is to leave the prediction error identification method unmodified, but to design a proportional plus integral observer by the same procedure as in
646
G. A. Van Zee and O. H. Bosgra
the design of proportional plus integral control. Further research in these directions is necessary to obtain a deeper understanding of the role of low frequency model errors. CONCLUSIONS The accuracy of prediction error identification methods has been investigated by comparing the predicted and actual control performance of a multivariable control system designed on the basis of the identification results. This is the most direct approach to model validation if the model is to be used in control system design. It appears that the accuracy of the initial model, obtained by stochastic realization, deteriorates when its parameters are optimized with respect to the usual one-step ahead prediction error criterion. This demonstrates that the considered criterion may be less suitable in the case of approximate modelling. We have shown that this problem is caused by the low frequency model errors which are important when considering the suitability of the models for the design of control systems incorporating integral feedback, but which are not reduced by the considered prediction error approach. Therefore we have proposed a modified prediction error identification criterion which has been demonstrated to reduce the low frequency model errors efficiently. It is only when using this procedure that satisfactory applications of multivariable optimal control have been obtained. With respect to the size of the considered parameter optimization problem (250 parameters in an 8-input, 8-output model) yielding computing times in the orde.r of hours, we conclude that much larger systems cannot be handled unless they are divided into subsystems which are modelled separately. APPENDIX Figure A.l shows a diagram of the pilot scale process which consists of 3 water flows which are thermally coupled by heat exchangers with variable levels. The energy is transmitted from the primary to the tertiary flow via the closed secundary circuit. Table A.l lists all the relevant variables with their values in the operating point for the experiments. The input u and output y are defined as variations around the operating point. u consists of the control input Uc and the measured disturbances ud' The sampling interval has been set to 4 [sl on the basis of preliminary response measurements. To prevent the hydraulic variables from drifting away from the operating point, the identification experiments have been performed under closed-loop conditions, using experimentally tuned single-loop PI controllers. The control loops have been defined b y feedback of the first 5 elements in y to the corresponding elements in u c ' For the test signals, which have been simultaneously added to all process inputs, we have used
computer-generated sequences which approximate normally distributed indepent white noise signals, having standard deviations tuned such that these yield standard deviations of the process outputs between 2 and 5 percent of their full range. These experimental conditions are sufficient to have the concerning identifibility conditions with respect to the experiment satisfied. REFERENCES Akaike, H. (1976), Canonical correlation analysis of time series and the use of an information criterion, In: R.K. Mehra, D.G. Lainiotis (Eds.), System identification: advances and case studies, Academic Press, New York, 27-96. Anderson, B.D.O., J.B. Moore, R.M. Hawkes (1978), Model approximations via prediction error identification, Automatica, 14, 615-622. Astro~ K.J. (1980), Maximum likelihood and prediction error methods, Automatica, 16, 551-574. CaineS: P.E., J. Rissanen (1974), Maximum likelihood estimation of parameters in multivariable Gaussian stochastic processes, IEEE Trans. Inform. Theory, 20, 102-104. Caines, P.E. (1976), Identification methods for stationary stochastic processes, IEEE Trans. Autom. Control, 21,500-505. Caines, P.E. (1978), Stationary linear and nonlinear system identification and predictor set completeness, IEEE Trans. Autom. Control, 23, 583-594. Denham, M.J. (1974), Canonical forms for the identification of linear multivariable systems, IEEE Trans. Autom. Control, ~, 646-656. Eklund, K., I. Gustavsson (1973), Identification of drum boiler dynamics, 3rd IFAC symposium on Identification and system parameter estimation, The Hague, NorthHolland Publ. Comp., Amsterdam, 87-108. Goodwin, G.C., R.L. Payne (1977), Dynamic system identification, Academic Press, New York. Gupta, N.K., R.K. Mehra (1974), Computational aspects of maximum likelihood estimation and reduction in sensitivity function calculations, IEEE Trans. Autom. Control, 19, 774-783. Gustavsson, I. (1969), Maximum likelihood identification on the Xgista reactor and comparison with results of spectral analysis, Report 6903, Lund Institute of Technology, Lund. Gustavsson, I. (1972), Comparison of different methods for identification of industrial processes, Automatica, 8, 127-142. Kwakernaak, H., R. Sivan (1972), Linear optimal control systems, Wiley & Sons, New York. Ljung, L. (1976), On the consistency of prediction error identification methods, In: R.K. Mehra, D.G. Lainiotis (Eds.), System identification: advances and case studies, Ac. Press, New York, 121-164.
647
Validation of Prediction Error Identification Results Ljung, L. (1978), Convergence analysis of parametric identification methods, IEEE Trans. Autom. Control, 23, 770-78-3-.-Ljung, L., A.J.M. Van Overbeek (1978), Validation of approximate models obtained from prediction error identification, Proc. 7th Triennial World Congress IFAC Helsinki 1978, Pergamon Press, Oxford, 1899-1906. Mehra, R.K . (1974), Optimal input signals for parameter estimation in dynamic systems - survey and new results, IEEE Trans. Autom. Control, 19, 753-768. Olsson, G. (1973), Modelling and identification of nuclear power reactor dynamics from multivariable experiments, 3rd IFAC Symposium on Identification and system parameter estimation, The Hague, NorthHolland Publ. Comp., Amsterdam, 375-384. Sage, A.P., J.L. Melsa (1971), Estimation theory with applications to communications and control, McGraw-Hill, New York. Soderstrom, T., L. Ljung, I. Gustavsson (1975), On the accuracy problem in identification, Proc. of the 6th IFAC world congress, Boston, paper 18.1. Soderstrom, T., L. Ljung, I. Gustavsson (1976), Identifiability conditions for linear multivariable systems operating under feedback, IEEE Trans. Autom. Control, 21, 837-840. Tse, ~, H.L. Weinert (1975), Structure determination and parameter identification for multivariable stochastic linear systems, IEEE Trans. Autom. Control, 20, 603-613. -Van Zee, G.A.,(1981), System identification for multivariable control, Doctoral dissertation, Delft University of Technology , The Netherlands.
., -'-'
'., 1.
,. ,! _
--
-_.
,.' -----I.
.- ----.:-:::..:--
Fig. 3.1
Closed-loop step respons e s of single-loop control
,
'"
,,:1_ _ _-
" :1_ --_,~---,~ ----
...t,j ____
-Fig. 3.2
Simulated closed-loop step responses of multivariable control for model MO
"
,, :1
< --
-
,. ,
I.
,,",
-----------
Fig. 3.3 120
Simulated c losed-l oop s tep resp onses of multivariabl e control wit h reduced Q for model MO a
40 "
O~----------L-----------L---------~
o
25
50
,,: I~ ---
75
" :1-----Fig. 2.1
J as a function of the iteration number i
-----
,:1 - - - ... ~ ----
Fig. 3.4
Actual closed-loop step responses of multivariable control with reduced Q for model
a
Mo
648
C. A. Van Zee and O. H. Bosgra , "
"
: ~ ----
" :1____ ,..J
~
~----
_ __
--
.... j - - -
Fig. 3.5
Fig. 3.6
Simulated closed-loop step responses of multivariable control with reduced Q for model M a
Actual closed-loop step responses of multivariable control with reduced Q for model M a
\
.,
~....
" .~~J;J.).;~""" ,'~
l02L-________~________~ o 25 50
Fig. 4.1
--
'0,,---
..
'.' I - -- /
~ ;\~"":"""'''''''' ' ' '
r")_':~'IA."'''''~
/~~;~,
,. :\~.J.·._ .".•.i''''/
._ ~_~J""
~
'" -~J' ~
~
.
u~
J I, ' ..,., ..Jv..,""\ J-..._
,: r~., .....
1'i1" 1
_-.-
"'l
.----', j..•,Aw
. - .!!!
Fig. 4.2
as a function of the iteration n6mber i
J
~
t
. ,Jo. ., !."" VI
I
"
Normalized amplitude of diagonal of sample spectral matrix of residuals a: model M, b: model ME
.''.
.... -- - - ~
,•.
....' .,1 •
Fig. 4.3
Simulated c los ed-loop step responses of multivariable contr o l for model M
E
Fig. 4.4 Actual closed-loop step responses of multivariable control for model M E
".",;
.'
-- - - - - - - - - - --, " C .
f
'.,
T" T,
..'
t "
..
T.,
r:;
",
",
~!
L
s
T
I
IFig. A.I
T.,
- - -----
I
__ -.J
T
:2
·.,1
Diagram of the pilot scale process
Table A.I
Defini ti on of process variables