Identification and Uncertainty Modeling for Multivariable Processes

Copyright@ IFAC Advanced Control of Chemical Processes Pisa, Italy, 2000 '

IDENTIFICATION AND UNCERTAINTY MODELING FOR MUL TIVARIABLE PROCESSES

Department of Mechanical Engineering, Bogazi~i University, Istanbul, Turkey E-mail: [email protected]

Abstract: The use of orthononnal functions and consistency relationships to identify models of multivariable processes and their uncertainty descriptions is discussed in this paper. The nonn-bounded uncertainty can be obtained from the covariance matrix of the parameters. Using a-priori consistency relationships, unstructured uncertainty may be converted to a structured one. The ~ analysis shows that this uncertainty description is much less conservative. Copyright © 2000 IFAC Keywords: Identification, Uncertainty mode ling, Consistency relationships

1. INTRODUCTION Identifying multivariable models of process systems and the uncertainty bounds for these models, is a crucial step in design of multi input multi output controllers. In many cases, the models of multi variable processes need to be more accurate than their SISO counterparts because of the possible problem of ill conditioning (Skogestad et.al. 1988), which has no SISO counterpart. Unlike the SISO systems, small errors in the gains of multivariable systems may cause instability with the designed controllers, especially for poorly conditioned systems. Element by element identification of the transfer function matrix using SISO identification techniques can give models which have wrong directionality, not suitable for controller design. Therefore, it is advisable to perfonn identification using MIMO identification methods (Andersen et.a!.1989). In addition, the design of the input signals and/or the use of iterative closed loop identification to account for the directionality of the process is might be a remedy (Koung and MacGregor, 1994), (Gaikwad and Rivera, 1996) (Cooley and Lee, 1999). Another approach is to try to identify both the process model and its inverse accurately through repeated experiments (Li and Lee, 1996), which improves the errors due to ill conditioning. The concept of 'pseudo-singular values' (Featherstone and Braatz, 1998) has also been used to detennine validity of the identified ill-conditioned models. However, the parameterization is a difficult issue with MIMO transfer function matrix models, and a

749

priori knowledge of the system structure is needed to parameterize such models (Ljung, 1987). In addition, most of the methods for estimation of the transfer function parameters being iterative, are not guaranteed to converge to the global minimum of the cost function. Recently, a new breed of algorithms called subspace methods, which identify the state space models of linear systems from input output data have been developed (De Moor and Overschree, 1995). Subspace methods have the advantage of not requiring the parameterization of the system matrices and are not iterative. However, it is not easy to obtain a measure of model quality (i.e. parameter covariance matrix) with these methods. The orthogonal functions, recently introduced to system identification literature by Wahlberg (Wahlberg, 1991), (Wahlberg, 1994) have some nice properties and do not suffer from the disadvantages discussed above. Especially, the Laguerre functions are very useful for identifying process systems with overdamped responses. Accurate a-priori knowledge of the dominant time constant of the process (which is not too hard to obtain for most process systems) improves the quality of the estimates significantly. The use of the powerful tools of robust control theory with identified models necessitates an estimate of the quality of the model, in terms of nonn bounds and uncertainty structures. The estimation of these from finite, noisy data alone is a very difficult problem. If hard model error bounds are sought for, the 'worst case identification' approach (Miikila et.a!' 1995) can be used, however the conservativeness of the model error bounds might be a problem with this approach. Stochastic embedding approach developed by Ninnes

and co-workers (Ninnes, Goodwin, 1995) generates probabilistic (soft) model error bounds, but to date, no straightforward multivariable extension of this approach has been published (Toffner-Clausen, 1996). In addition, it is not easy to translate the classical confidence intervals of parameters to norm bounds, used in robust control. Motivated by the recent developments on probabilistic robustness analysis (Stengel and Ray, 1991), (Tempo et.a1.1996), the pragmatic approach of using Monte-Carlo simulations to construct the model uncertainty bounds, will be taken in this paper. With this approach, it is possible to make statements such as 'we can say with 99% confidence that the model uncertainty bound developed, contains 99% of the plant perturbations'. Such accuracy, it is believed, is adequate for most applications. To improve the model quality, the available a-priori information about the system might be combined with the identification process. This can be done, as described in (E~kinat, 1995), by expressing the available information about the system as constraint equations and then solving a constrained minimization problem to obtain the system parameters. Of special interest for process control is the use of consistency relations (Skogestad, 1991) as a-priori information. In addition to improving the parameter estimates, a-priori information allows us to impose a structure on the uncertainties, thus making them more realistic and tighter.

(5)

and:

f/>(/) = [L 1 (q)u(t)

Il> = [f/>(l)

.. .

...

LM (q)u(t)y

(6)

f/>(N)Y

Assuming that the input u(t) is persistently exciting, the parameters e may be estimated through least squares as: 8 = R-1Il>Ty (7) R=Il>TIl>

(8)

The transfer function matrix estimate, G(q,8) will be a linear function of 8:

W(q) = [L 1 (q)

...

(9)

LM (q)y

The covariance matrix of the estimated parameters is given by:

P = Cov(l}) = Cov(vec(8» = ~2 ®R~2

=diag(a

2 j

)

1

(10)

i = l...ny

® denotes the Kronecker product. a/ can be estimated as the sum of the prediction errors for output y;(t). It is well known that (Ljung, 1987), the quantity,

For identification, we assume that the data generating system can be accurately modeled by a finite expansion of orthonormal basis functions:

i

is distributed with np=nyxnuxM degrees of freedom. This allows us to define confidence ellipsoids for scalar transfer functions, if the bias error is negligible. This does not easily carry over to multi variable case. In addition, the norm bounds used in robust control cannot be easily estimated from (11), singular values being complicated functions of matrix elements. One approach to estimate uncertainty bounds is Monte-Carlo simulation: Generate a large set of random vectors c5 of size np, with covariance of pI, replace by + c5, and estimate the additive model uncertainty bound Wa as:

M

L BkLk (q)u(/) + e(/), 1= l..N

(4)

Y = [y(l)

G(q,8)=8 T W(q)

2. IDENTIFICATION OF MODELS AND UNCERT AINTIES

Y(/) =

where:

(1)

k=l

where y and u are the output and input vectors of size ny and nu; Bk (nyxnu) is the matrix of parameters to be estimated, and e(t) are the disturbances of zero mean and variance of a,2 each (i=l ... ny). For

e e

Laguerre models, the fixed basis functions are given by: (2)

(f denotes the maximum singular value. Equation (12) gives a soft bound on the uncertainty, but this is unavoidable with classical identification methods. From (Tempo, et.al. 1996), a bound for the number of Monte-Carlo samples NMC of c5, which guarantee:

a is the dominant time constant and T. is the sampling time. With some manipulation, (Ninnes, et.aI.1995) the above equation (1) can be put into the form : Y=1l>8+E (3)

750

D=V-L

; B=F+L-V.

(19)

Substituting (17) to (18) and taking derivatives with respect to L and V, we obtain:

is given by: (14)

For example, for £,=£}=0.01, we obtain N MC =458, independent of the number of parameters np. This approach requires low computational effort with good accuracy.

3. CONSISTENCY RELATIONSHIPS The consistency relationships, first used by (Haggblom and Wailer, 1988) later elaborated by (Skogestad, 1991) provide means for establishing relationships between steady state gains of multi variable processes. Very simple material and energy balances are used in order to obtain this information. As an example, consider a counterflow heat exchanger. Suppose the outputs of the system are the exit temperatures of the hot and cold fluids ~ nI ~ and the inputs are the mass flow rates of the two streams ;,. and m,. If no heat is lost to the environment, an energy balance on the two streams give:

.

a7;.2 + . am a7;.2 + a me h

mh Cph - -

(21)

These relations may be incorporated into the identification process. With some manipulation, (21) can be written as:

e

ph

is the where covariance of if is:

h

m e Cpc - -.-

f3

(22)

unconstrained

solution. The

Cov(if) = (p-I - p-I AT (AP- 1AT) -I AP-')

me

e.

The partial derivative terms represent the steady state gains of the process at the operating point. As a second example, consider the binary distillation process. The following material balances always hold for a distillation column:

F=D+B

f3

Then, the above linear equations can be enforced in identification, using constrained least squares estimation. The constrained estimates obeying (22) are given by:

a7;,2 = C (7;.1 - 7;.2) am (16) a7;,2 =Cpc(7;,1 - 7;,2) a

m e Cpc - -.-

G(I)T a =

A() =

Differentiating this equation with respect to the mass flow rates gives the following two equations: mh C ph - -.-

Again, we have two linear relationships between the steady state gains of the process. With the equations above, two relationships between the four steady state gains of the process are established. Although none of the gains is known, the above relations may be used to improve the quality of the estimates. The linear relationships may be expressed (for the discrete model) as:

(24)

where P is the co variance of A more efficient solution for if, using the null space of A is also possible. See (E~kinat, 1995) for details.

4. MODELING OF UNCERTAINTY

(17)

The consistency relationships may also be used to determine the structure of uncertainty. Assuming additive unstructured uncertainty ~au, the consistency relationships for the uncertain system (for continuous time) may be expressed by:

where F is the feed rate, D and B are the distillate and bottoms flow rates, Xo and XB are the top and bottom compositions and ZF is the feed composition. Suppose the conventional L- V (reflux-boilup) structure is selected for the control of the compositions Xo and XB. With the equimolal overflow assumption, the external flows D and B are related to the internal flows Land Vby:

75 1

If G obeys (21), then, the following equation is valid for ~a :

In this case, some of the elements of !!."., will be dependent and can be solved in terms of the remaining terms. Using equation (26), unstructured uncertainty can be converted to a structured one ~as· For example, for distillation process, if Wa is diagonal with identical elements, ~au may be expressed as:

~ Hence, Wa ~ar

=[-DIE 1

DU

° °]

-BID][llJ 1

II

= Cll

2

ar

(27)

the steady state case, it can be extended to dynamic case, with some assumptions. The fundamental balances such as (15) and (18) are always valid, whan an accumulation term is added to them. Therefore, when the accumulation inside the control volume is teken into account, it should be expected that the result obtained in (27) should carry over to all frequencies . Linearizing expressions (15) and (18), using harmonic inputs u(t) = Uo e iWl , and assuming outputs to be harmonic with the same frequency with u, the following expression can be obtained:

/3

(28)

'tj

Consider the simple model of the ill-conditioned distillation column thoroughly analyzed in (Skogestad and Postlethwaite, 1996): G(s)

1 75s

[0.878 - 0.864]

+ 1 1.082 -1.096

(30)

This is a model ofa distillation column withxD=0.99, XB=O.OJ, D=B=0.5, L=2.706, V=3.206. The consistency relations require the column model to obey:

G(iro) T(75iro + 1)[°.5] =[ 0.98 ] 0.5 -0.98

(31)

This system was simulated by two independent PRBS inputs of magnitude 1 and discrete white noise of variance 0.01 was added to both outputs. Sampling time is taken as unity. 500 simulations of 3000 data points each, with different noise realizations were done. The results presented are the averages of the simulations. Laguerre models with 4 parameters for each transfer function element (total 16 parameters) was estimated. Laguerre time constant was taken as a=0.98, which corresponds to 't = 50 in continuous time. Despite the error in the time constant, the models compare very well with the actual process. Since the commercial software available for evaluating robustness uses continuous time models, a conversion of discrete time models to continuous time was necessary. Conversion between discrete and continuous models is done using Tustin approximation, as done in (Toffner-Clausen, 1996). Conversion to continuous time results in eighth order state-space models.

Therefore, the analysis in equations (25)-(27) applies for all frequencies. a(im) in equation (28) usually has the form a(iw) = ion i + 1, where

where S is the sensitivity function and Wp is the performance weight. ,1 is obtained by augmenting ,1ar with the full block ,1p. Robust stability of the closed loop system may be checked by finding JlRS of M J J = -WaKS, with respect to the uncertainty ~as and robust performance by JlRP of M with respect to ,1 = diag(,1ar, ,1pJ. In a sense, models which have JlRP
the uncertainty may be formed as where Wa = WaC . Although this result is for

G(iro)T a(iro) =

(29)

is the

dominant time constant for the i'th loop. 'tj can be roughly estimated either from an identification experiment or using analytical expressions as done in (Skogestad and Morari, 1987). Using expression (28), the unstructured uncertainty may be converted to structured one. In general, for a nXll system, it is possible to eliminate n elements of the uncertainty description, using consistency relations. This is most useful in 2x2 systems, where the uncertainty can be converted to a diagonal one. The structured uncertainty gives a much tighter uncertainty description, by enforcing the identified model to obey some fundamental balances such as (15) and (18). As is well known, at higher frequencies the model uncertainty for any plant becomes totally unstructured, but this effect can be modeled by another uncertainty block.

A comparison of two different models obtained using the same data is made. The first model (G J) is estimated using the simulation data only; in the second one (G2Y , the consistency relations are used along with the data. The estimated steady state gains are as follows :

The complex structured singular value Jl (Packard and Doyle, 1993) can be used to evaluate the robustness of the system with a designed controller K. Using additive uncertainty, the augmented plant M in the M-,1 structure (Skogestad and Postlethwaite, 1996) becomes:

752

G == [0.8754 I

- 0.865] -1.093

1.08

G

== [ 0.878

2

- 0.864 ] 1.0806 -1 .0936

The estimated gains are very close to the actual ones, so it can be concluded that the bias error is negligible. Both of the models have the correct input directions. (Actually, some of the 500 simulations give models with wrong directions, but the results above are their averages.) The additive uncertainty bounds estimated using 500 Monte-Carlo simulations and equation (12), is given in Figure 1. It can be seen that the use of consistency relations somewhat improves the uncertainty bound. Additive uncertainty weight Wo is obtained by fitting a third order model to the curves in Figure 1. An inverse based controller of the form : K(s)

= 0.1(7Ss + 1) [6(0)J1

(32)

s is designed and is applied on the identified models. It is well known that, inverse based controllers are not robust to input uncertainty (Skogestad et.al. 1988), but our aim is to evaluate robustness to model uncertainty. Inverse based controllers provide a simple and effective mean for this aim. The performance weight Wp(s), compatible with the designed controller, is selected as:

w (s) == O.5(lOs + 1) p

lOs

(33)

because of the sensItIvIty to possible input uncertainty. In identification, the main difficulty is the possibility of identifying the wrong direction associated with the weak inputs (i.e., internal flows in a distillation column). This problem is usually circumvented by the design of inputs so that the weak input direction is properly excited (Koung and MacGregor, 1994). In this paper, another means for improving the models and their uncertainty descriptions is presented. Some fundamental mass and energy balances exist in process systems. The consistency relations derived from these balances are always satisfied for analytically derived models. For identified models, they are not automatically satisfied, unless they are enforced on the identification process. In this paper, these relationships are combined with the data, using constrained least squares estimation. More importantly, the uncertainty descriptions for these systems become much tighter when the consistency relations are imposed. It was mentioned previously that, the variation in gains of many multi variable processes is dependent (Skogestad, 1991), (Jacobsen and Skogestad, 1994). With this paper, it is shown that this is true for the dynamic case as well. This fact allows us to define better uncertainty descriptions to be used in controller design. Uncertainties are tighter because they are made to obey physical conservation laws. It is believed that this approach will be useful not only for identification, but for modeling for robust control in general.

REFERENCES

The robustness analysis is performed on the augmented plant M. using the Jl-toolbox. For the first model, the model uncertainty is taken as the full block Liau • For the second, it is structured with a diagonal Lias , with the uncertainty weight taken as Wo = Wo C , as defined in equation (27). Results are given in Figures 2 and 3. It can be seen that there is a drastic difference between the plots in Figures 2 and 3. In Figure 2, representing the model and the uncertainty description where the consistency relations are not used, JlRP is around 12, and the dominating effect is the model uncertainty, as can be seen from the robust stability plot. The model seems to be unacceptable in this case. Whereas, if the consistency relations are used, a more realistic picture of the uncertainty can be obtained, as given in Figure 3. In this case, JlRP is less than 1, and the dominating effect is not the model uncertainty, but rather the performance weight. This result agrees with the conclusions reached by (Skogestad, et.al., 1988) about the effect of model uncertainty.

Andersen H.W., Kummel M., Jorgensen S. (1989), Dynamics and Identification of a Binary Distillation Column. Chem. Eng. Sci. Vol. 44, pp. 2571-2581. Cooley, B., Lee J.H. (1998) Integrated Identification and Robust Control, Journal of Process Control, Vol. 8, pp.431-440. Cooley, B., Lee J.H. (1999) Control Relevant Experiment Design for Multivariable Systems Described by Expansions in Orthonormal Bases, Manuscript Submitted to Automatica De Moor B., Van Overschree P. (1995), Numerical Algorithms for Subspace State Space Identification. In "Trends in Control, a European Perspective" A. Isidori (Ed.) pp. 385-422. E~kinat E. (1995) System Identification Using Constrained Estimation. Proceedings of 3 'rd European Control Conference, Rome pp. 856861. Featherstone A., Braatz R.D. (1998) Integrated Robust Identification and Control of Large Scale Processes Ind.Eng.Chem.Res. Vo1.37, pp.97-106. Gaikwad S.V., Rivera D.E. (1996) Control Relevant Input Signal Design for Multivariable System Identification: Application to High Purity

6. CONCLUSIONS Ill-conditioned processes are usually difficult to identify and control. Control is difficult mainly

753

Distillation, Proc. IFAC World Congress, San Francisco, pp. 349-354. Haggblom K.E., Wailer K.V. (1988) Transformations and Consistency Relationships of Distillation Control Structures. AIChE Journal, Vol 34, pp. 1634-1648 Jacobsen E. , Skogestad S. (1994) Inconsistencies in Dynamic Models for Ill-Conditioned Plants: Application to low order models of Distillation Columns. Ind.Eng.Chem.Res. Vo1.33, pp.631640. Koung c., MacGregor 1. (1994) Identification for Robust Multivariable Control: the Design of Experiments. Automatica, Vol. 30, pp. 15411554. Li W., Lee J.H. (1996) Frequency Domain Closed Loop Identification of Multivariable Systems for Feedback Control, AIChE Journal Vol. 42. pp. 2813-2827. Ljung L. (1987) System Identification Prentice Hall. Miikilii P., Partington J., Gustafsson T. (1995) Worst Case Control Relevant Identification, Automatica, 31 , pp. 1799-1819. Ninnes B., Goodwin G. (1995) Estimation of Model Quality. Automatica, 31, pp.I771-1797. Ninnes B., Gomez J.C, Well er S. (1995) MIMO System Identification Using Orthonormal Basis Functions. Proc. CDC'95. Packard A., Doyle J. (1993) The Complex Structured Singular Value. Automatica, Vol.29 pp.71-109 Skogestad S., Morari M. (1987) The Dominant Time Constant for Distillation Columns. Comp. Chem. Engng., Vol.ll, pp.607-617 Skogestad S., Morari M., Doyle J. (1988) Robust Control of III Conditioned Plants: High Purity Distillation. IEEE Trans. Automatic Control. AC-33 , pp. 1092-1105 Skogestad S. (1991) Consistency of Steady-State Models Using Insight about Extensive Variables. Ind. Eng. Chem. Res., Vol. 30, pp. 654-661. Skogestad S., Postlethwaite (1996) Multivariable Feedback Control, Wiley. Stenge! R., Ray L. (1991) Stochastic Robustness of Linear Time Invariant Control Systems, IEEE Trans. Automatic Control V. 36 pp. 82-87. Tempo R., Bai E., Dabenne F. (1996) Probabilistic Robustness Analysis: Explicit Bounds for the Minimum Number of Samples. Proceedings of the 35 th CDC Kobe, Japan pp. 3424-3428. Toffner-Clausen S. (1996) System Identification and Robust Control. Springer Verlag - Series in Advances in Industrial Control. Wahlberg B. (1991) System Identification Using Laguerre Models. IEEE Trans. Automatic. Control. Vol. 36, pp. 551-562. Wahlberg B. (1994) Laguerre and Kautz Models. Proc. IFAC System Identification, Denmark, pp. 965-976.

754

0 . 14r---~---~--~------'

01 2 0.1 'W.

).08 0.06 0 .04

0.02

Figure 1. Model uncertainty bounds for G j (dashed) and G2 (solid).

10 8 6 4

HP

._ --

_._- _

.._

'-- '-

'.--: :-- - - " :.

oL-----~-----~---~~--~ 1~ 1~ 1~ 1if 1~

w

Figure 2. ~ plots for the model identified without consistency relations. Solid line: J1.RP, Dotted line: J1.RS, Dash-dotted line: J1.NP·

0.8 ~

RP

\ . ./' - - .-=---

0.6 HP

./

0.4

02

0 10"

RS

10.2

10·

10'

Figure 3. ~ plots for the model identified with the consistency relations. Solid line: J1.RP, Dotted line: J1.RS, Dash-dotted line: J1.NP·