Identification of Distributed Parameter Systems Using Polynomial Identification Techniques

Copyright © IFAC Identification and System Parameter Estimation. Budapest. Hungary 1991

DISTRIBUTED PARAMETER SYSTEM

IDENTIFICATION OF DISTRIBUTED PARAMETER SYSTEMS USING POLYNOMIAL IDENTIFICA TION TECHNIQUES P. J. Gawthrop* and M. T. Nihtilii** *Department o[ Mechanical Engineering, The University, Glasgow G 12 8QQ, U K **HeLlinki University o[Technology, Otakaari 5A, SF 02150 Espoo, Finland

Abstract Certilin distributed parameter systems may be approximated by finite dimensional systems which are, how~ Vt'r polynomial in unknown parameters. Algorithms are given for estimated unknown parameters of such ,ys tems

Keywords

Shinskey [5) makes the distinction between interacting and non-interacting lags in the context of process control. R

R

R

v~

Introduction

Typical industrial processes have transfer functions arising from the effect of many sub-processes. Such trillls[er [unctions are necessarily of high order. For the purposes of modelling, identification and control, it has been customary to use reduced-order models comprising rational transfer functions of low order. Low order approximate modelling clearly leads to a 31llallllumber of parameters to be adjusted in the modelling and identification phases, and this is an important motivation for its use. However, the converse is not. true: a small number of parameters does not imply low order modelling. As discussed in the sequel, N non-interacting lags can be conveniently modelled by the one parameter transcendental transfer function

(1) and N interacting lags by the one-parameter transcendental transfer function

e- V7i

Industrially Relevant Transcclldental Transfer Functions

~ T T

Parilllle ter estimation; recurstve estimation; distributed-parameter systems; least-squares 2st.illlat.ion; polynomial-in-the-parameter systems.

1

2

Figure 1: Interacting RC circuit

Figllfe 2: Non-illteractillg RC circuit The electrical analogy he uses is the sequence of N RC circuits in Figures 1 and 2 respectively. In the former case, the individual RC circuits do not interact; in the latter case they do interact: the current in each capacitor depends on the voltage across both adjacent resistors.

(2)

Although such models are parsimonious in the parameters, they are not linear in the parameters and so are not amenable to linear least-squares algorithms. The pllI'pose of this paper is to present parameter identification algorithms appropriate to these transcendental transfer functions. The estimation of time delay (e-,T) systems has been considered by a number of authors [1) [2) [3). This paper builds on the polynomial approach of Nihtilii, [4) used in [3) as a basis for the estimation of time delay systems, to provide a general algorithm for parameter estimation for single-input single-output. Systems containing one or other of transcendental transfer functions described above. The outline of the paper is as follows: Section 2 considers some typical industrial processes which give rise t.o transcendental transfer functions. Section 3 provides a general setting for the parameter estimation problem and section 4 gives some specific algorithms. Section 5 contains implementation details including variable step size integration routines and the implementat.ion of time-varying transcendental transfer functions. Section 6 contains some simulations and section 7 concludes the paper.

1145

Figure 3: Interacting tanks The hydraulic analogy that Shinskey uses is a sequence of N tanks containing liquid; the non-interactive version appears in Figure 4, the interactive version appears in 3. Again, tlte net. outflow frolll each tank depends on the relative levels in both adjacent tanks in the latter case. In each case, the Nth order differential equation can hc ilpproXilltatC'd Ity it 1'1l1'linl differential equation and hence by a transcendental transfer function. The two cases will he treat.ed separately. In each case, the electrical analogy will be considered but the same equations resnlt from hydraulic systems.

2.1

Non-interacting lags

Letting Vi be the volt.age on the ith capacitor then Cdvj 1 ---;It "R(Vi-l - v;) (3)

=

·.J--L

3

----~':, '

The algorithm presellted in carlicr papers [4] [2] [3] refers to the general error model

Li

e(t) ::.:. "

De fining the continuous space coordinate x(O wit.h increment ~ 1 J' = ~; ~ = N

J(O(t), t) =

2

( 4)

(5)

RC

(6)

and T = RC, the time constant of one lag, Setting v = V(x)e't, the partial differential equation (5) becomes th e following ordinary differential equation

+ dV dx

= 0

(7)

(8)

The transfer function relating V(l) to V(O) is thus:

= e-,T, V(O)

(9)

Tlllls, as N becomes large, the transfer function of N non-interacting lags can be approximated by a pure time-delay with value

T\ = NT where T is the time constant of each lag ,

2.2

(10)

Interacting Lags

With the same notation as before: CdVi 1

-;u-

=

n(Vi - 1 -

Vi

+ Vi+1

-

Vi)

(11)

Taking limits as before , the following second-order partial differential equation (PDE) is obtained ov

02v

T2Ft + ox 2 = 0

(12)

where

T2 = RC ~2

where f3 is the non -lIcgat.ive scalar fOT'ge lling factor:

f3'2,0

(19)

and So is the posit.ivc defillit.e IlIat.rix initial cost weight-

(20)

So> 0

ALGORITHM d -

-

-B(t) = -JI(O(t) , I) dt

-I

> 1, Ji (O( t), I)

and for i equation

III(t)

(21 )

is gi ven hy the differential

where

H(t) •

= ~ oig2(y(t), ~t(t), 0(1)) 2

DOi The only non-zero init.i a l co ndit.i o ns a re:

(23)

= Ba

(24)

h(O, 0) = Su

(25)

0(0)

As discussed in t.hat pap er [3] this algorithm involves 2N differential equat.ions alld th e ith equation involves i dimensiollal matrices where N is o rd er of the highest order non-zero derivative of g(y(t), u(l), !J(t)) with respect to !J,

An cxtcnded CITOl' lll<)(ld In this paper , th e [lolynollli a lme t.hod is extended to include single-input single-output syst.ems desc ribed by a rational transfer functi o n ~ in series with a (possibly transcendental) t.ransfer functi o n of t.h e special form

(13)

= V(x)e't

yes)

we have

d2V

+ dx 2 = 0

(14)

Strictly speaking, the second order differential cquation (14) has two solutions, But from causality considerations, one of these solutions may be discarded to gIve

\'(1)

(18)

0

D(sT) :

= N 2T

Once again, setting V

s1 2V

(17)

dJi(O(t), t) , dO(t) dt +f3J.(O(l), t) = Hi(t)+J i + 1 (!J(t), t)di(22

and the solution is: V = e-,T,x

V(I)

00 ]

is minimised by the following rec ursive

= IS: = N RC = NT

s7 1V

~e -!1 t[O(t) - oaf 5 0 [O( t) -

mg:

Du

Ft + ox = 0

where

TI

( 16)

" • +-I1t e- !1(t - T) g-(Y(T) , !leT) , O{l))dT

< x < 1)

alld taking the limit as N -+ 00, we obtain the firstorder partial differential equation (PDE) 1\

= g(y(t), lI(t) , O(t))

where e(t) is the scalar e rror , yet) , !let) are the scalar system input and out.put. and B(l) is a vector containing the unknown paralllct.e rs to b e es t.imat.ed, Furthermore, to obtain a finit. e climensional optimal algorithm, g must be a polynomial in B, This general model has been applied t.o the part.icular case o f t.im e-delay identification, [3] As shown in [4] [2] [3] t.h e least-squarc cost function

Figure 4: Non-interacting tanks

• . iJu

The General Least-squares AIgoritluTI

= e-,;;-r,V(O)

( 15)

1146

= ~~:~ D(sT)ft(s)

(26)

Note that this impli es zero system initial conditions, The key idea in applying th e polynomial method to th e new transcen
In other words, for the purposes of taking derivatives, D(sT) is approximated by an Nth degree polynomial in T; however, when D(sT) itself is to be used, an exact representation is used. As discussed in [3]. the basic model (26) must be transformed into a form more suitable for application of the polynomial method . This is achieved by filtering equation (26) by ~ to give:

_ A(s) _ B(s) _ c(s) = C(s) yes) - C(s) v(s) \V

= D(sT)u(s)

D(sT)

= e-·

go(t)

(29)

[}'D(sT) = (-s)iD(sT) [}T'

In the case where VesT) = e-..J7'f

[}D(sT) = -~(s/T)1 D(sT) (31 )

III prillciple, the method is applicable to other systems of this form. It. is convenient to reparameterise the error model as:

(32)

where

(43)

Distributed Lag

(30)

+ Or Xb

(41 )

Time delay In the case where V(st) = e-'1'

and in the case of the distributed lag

c(s) = OT X = O~ Xa

= g(y(t), tI(t), B(t»

(40)

where L - I dellotes the illverse Laplace t.ransform.

The term e( s) corresponds to the measurement and modelling errors not included in equation (26). Both models presented in section 2 can be rewritten in this form. In the case of the time delay : T

'1 !.

f3ij = J.'1('1 - J')1. and the scalars 9i(t) are given by

(28)

here

,'(s)

where f3ij is thejth coefficient of the ith degree binomial polynomial:

aT

2

(44)

Higher derivatives are more complicated and can more conveniently obtained hy a computer algebra package such as Reduce. [G] For exalllple, the first 3 derivatives can be obtained lIsillg the following Reduce interaction: REDUCE 3.3, 15-Jan-88 ...

1: off nat; (33)

of =

[b o; bl ; .. . bn ]

(34)

Xa = V(S)y(S)

(35)

D := l/E**SQRT(S*T)

.\'b = V(s)D(sT)u(s)

(36)

3:

and \'(s) is the state-variable filter transfer function:

f / C( s)

\ (oS) = [s"; s" - 1; .. . ; 1

Particular Algorithms

As discussed in the previous paper [3] the exact algorit.hm is not implementable as an infinite number of equations are required; however, approximate algorit.hms were shown to have good performance. With this in mind, selected algorithms devised in [3] for the special case of time delay identification (D(sT) = e-· T ) are extended to the more general case. In each case, the algorithm is stated for arbitrary D(sT) and then specialised to the two particular cases of D(sT). More details of the assumptions and approximations are given in [3]. The algorithms are numbered to correspond to

[3].

4.1

Algorithm 2 (Identification of T only)

= !ti(t)

(38)

and t.he scalars h,(t) are given by:

1 [}'gO(t)2 l' h,(t) = 2 [}T' = 2 Lf3ij9j9,-j

(2+SQRT(T)*SQRT(S)*E**SQRT(S*T» 4: df(D,T,2); (SQRT(S) + SQRT(T)*S)/ (4*SQRT(T)*E**SQRT(S*T)*T) 5: df(D,T,3); - (S*(3+SQRT(T)*SQRT(S) + S+T + 3»/ (8+SQRT(T)*SQRT(S)*E**SQRT(S*T)*T**2) The determination of 9i(t) from 9i(S) is considered in the next section.

Derivation As 0 = T is a scalar, the formula of Leibnitz [7] may be used to give hi(t) in terms of 96(t). Otherwise the algorithm corresponds to the general algorithm of section 3.

4.2

Algorithm The algorithm is a scalar version on the general algorithm of section 3 where:

lIi(t)

df (D, T, 1) ;

- S/

(37)

The order of the filter polynomial C(s) depends not only 011 the order of the rational part of the system, but also on the order N of the approximation in (26).

4

2: D := exp(- \sqrt(s*T»;

Algorithm 4 (Identification of T and A(s))

Algorithm Defining

B=

(~)

( 45)

HI(s) is partitiolled as: (39)

H (s) = ~ ae(~)2 = ( I

;=0

1147

2

[}O

hl_(s) ) Xa90(S)

(46)

where the scalar hl(S) is as defined in Algorithm 2. Di ITe rcnt.iating once again with respect to 0:

Lernma Assuming that the inverse Laplace transform fT(t) of the causal t.rilnsfer fllnction F(sT) exists then:

( 47)

h(t) For i :::: 2, the only component of Hi(t) that is not zero is that corresponding to T and is given by:

_ ( h;(s) lli(s) - ()X g,-1 S a

( 48)

ell

.

+j3Ji ((}(t), t)

Proof By definition:

h(t) =

.

dT

= Hi(t)+J i +1((}(t), t)Tt(49)

(56)

where f(t) is the inverse Laplace transform of F(s).

Thus, in this case, equation (22) becomes:

dJi(O(t), t)

t

1

= 1'f(1')

100 e't F(sl')ds

Let s' = sT and t' = 1

f oo

f T () t = l' io

The det.ermination of Hi(t) from H;(s) is considered in Ill<' 11.'xt section.

,f ,' t'

c

(57)

t.h e n ,

F(s )d,;

,

= ~f(t') = ~J(t /'/') Th e formulae for H;(t) follow from repeated differentiation of e(t)2 with respect to O(t). As in algorit.hm 2, the fo rmula of Leibnitz is used to obtain an expression for hi(l).

5

Implementation

ElTor derivatives illll'ir11lcnl(llion of the algorithm requires t.h e nu IJlcrical evaluat.ion of the derivatives of I.he error wit.h r,'sl"'ct. 1.0 I.he parameter vector e. To illustril.i.e the problellls involved, and their solution, consider identification of T when:

= 1; B(s) = 1; D(sT) =

..1(s)

e-V?i'

(50)

ogo(s)

!J 1 S

=----ar-

(51 )

_ -1

s 1 -V?i'C(s) ..;sre u(s)

(52)

- 2

D(sT) = c-·,T Following the sallle

cb!

F(sT)

= _l_e-V?i'

(53)

..;sr Alt.hough F(sT) represents a linear filter there are t.wo implementation problems. 1. The filter is irrational due to both the exponential

The first problem is solved by implementing

= F(sT)u(s)

F( sT) = e-'l'

(62)

thus and

= h(t) * u(t) = l~-oo h(t -

= 8(t -

1)

(63)

so 1 t h(t) = 1'6(1' - I) = b(l - 'f')

(61)

In this particlllar case, t.llell, the convollltion is trivial

h(t)

* 11(/)

= u(l - 'J')

(65)

This procedure generalises to the following Derivative HIgol'i HIIlI 1. Symbolically derive !Ji(S) by taking the itl! deriva-

tive of 90(S) with res pect to T 2. Rewrite t.his in tile forll1:

(66)

3. Take the illverse Fomier transform of each Fi(s) to give fi(t) (i.e. for '1 '=1). 4. Obtain fi1,(i) = t/;(I/T)

(54)

5 . Obtain each !/i(l) by creating u'(l) by convolution of fiT(t) and u(t) and passing 11'(1) through the linear SVF with transfer funct.ion V(s).

T)U(T)dT (55)

Numerical aspects The algorithms have been implemented in Matlab [91 code. A particular IlUlllerical feat.ur e of interest is the choice of integrat.ion algorit.hlllto int.egrat.e the differen .. tial equations describing t.he algorit.hl11 'Ve have fOllnd that it is advantageolls t.o llse a variable step size integration algorithm to t.uk e account of the initial rapid , and final slow, changes to parameter values. This is illustrated in t.he simulation

in the time domain by

l£'(t)

as before ,

(61 )

9i(S) = Fi(sT)u(s)

and square root term. 2. The filter is a function of the estimated parameter T.

ii'(s)

(60) ar~lIl1wllt.S

(-s r -s ·/' . gi(S) = C(s) e u(s)

f(t)

This can be regarded as two filters in series operating 011 the signal u(s). The first filter is which can be easily implemented as a state variable filter. [8] [3] The second filter is

(59)

In other words, the function f(t) (which is independent of T) can be ge nerat.ed once and for all at the commencement of th e es t.imation. Th e function fT(t) may then be created by interpolation at each time instant corresponding (.0 t.h e current est.imat.e of T. In what sense is this a generalisation of our previous results pert.aining t.o the pure delay given in the following equation?

In t.h e s domain:

.. ()

(58)

The impulse response h(t) can be readily found from F(sT) by the inverse (fast) Fourier transform. The second problem is now that fT would have to be reevaluated at each time step as the estimate of T varies - a time consuming process. However, this process can be accomplished without using the inverse Fourier transform by utilising the following simple lemma.

1148

6

SiInulation

TI1t' r,dluwil1!,; sil1ll1l .. t.ion is designed t,o indicate the im-

portance or a variable step size algorithlll for this sort of algorit.hm. These simulations programmed ill Mat"d) , illld implemented on a SUN 3/60 workstation using I'RO-l\IATLAB. The integ ration algorithm is the modified 'odw15'. The system considered IS

2.5

(67)

y(s) = e-'ii(s)

and the algorithm 2 is used with N=4. In each case, the exponent (J was 0.2. The polynomial C(s) = (1+0.2s)5.

1.5

J

°:1~~

~r- ·r -~-~ -~~ -~~~-

II 05 0 -

()

~

""-~/

10

J

/

~tl'P

25

20

15

Integration

r

inJex

Figure 7: Integration step size: SO

= 1e-7

·0.5 I ·

095 16

IX

20

0 .9

Time

0.H5 . 0.8

....

Figure 5: Input and output data

0.75 0.7 0 .65 0.6 0 .55

O'9·~t

0.5

2

12

10

14

16

Time

09

OHS

Figure 8: Estimated delay: SO = 1e-9

0.8 ·· ....

0.75 -

3.5

0.55 O . 5~--:--~---:---:,::,----:,::----:-:--...,.".----:,::-----: 2 10 12 14 16 lR 20

2.5

Time

Figure 6: Estimated delay: SO = 1e-7

1.5

=

Fignres 6 and 7 correspond to So 10- 7 , Figures 8 and 9 correspond to So 10- 9 . Figure 5 shows the raw data used; Figures 6 and 8 show the evolution of the delay estimate l' with respect to t.illle. Figures 7 and 9 show t.he integrat.ion algorithm st.ep size. Note that the convergence is much Illore rapid in the latter case (So 10- 9 ) . But this rapid convergence reqllircs it very small step size (0.001) illitially, With it fixed step size algorithm, 20/0.004 = 5000 iterat.ions

=

0.5 -

IllIegrati on stt"r indc:'l.

=

1149

Figure 9: Integration step size: SO = 1e-9

18

20

lI'ollld be needed; the variable step size algorithm reqllire's only 18 steps.

7

Conclusion

The method of Nihtilii and Gawthrop [2] [3] for identification of time-delay systems e- 3T has been generali s (~ d to include identification of systems containing distributed lags e-F?i'. In particular the use of convollltion and symbolic manipulation to generate error derivatives has been introduced. The industrial relevallce of this has been noted. ,t.;illlIJiat.ioll results obt.ained IIsing MatLab s ll ()w t.he i Ilq )()J't.ance of variable step-size integration algori thllls.

References [1]

J . E. Marshall, Identification Strategies for Timedelay Systems, in 11 . .Unbehauen, editor, Methods and Applications in Adaptive Control, pages 141ISO , Springer, 1980.

[~]

1' ..1. Gawthrop and M.T. Nihtila, Identification o f time-delays using a polynomial identification Illdhod, Systems and Control Letters, 5:267 - 271, 1!)1'5.

[:J] P ..!. Gawthrop, M.T. Nihtila, and A. BesharatiHad, Recursive parameter estimation of continuoustime systems with unknown time delay, C- TA 1", 15(3),1989.

[4] I\[T. Nihtila, Optimal finite dimensional Recursi ve Identification in a Polynomial Output Mapping Class, Systems and Contro l Letters , 3:341-348, 1983.

P',]

F .e;. Shinskey, Process Control Systems, lIill , ]979 .

[li]

C;.

Mce~raw­

H.ayna, Reduce: Software for Algebraic tation, Springer, 1987.

COlll]lU-

[7] W.L. Ferrar, Differential Calculus, Oxford Universi ty Press, 1960.

[8] P.C . Young , Parameter Estimation for Continuoustill1e Models H)81

[a]

c.

A Survey, A utomatica, 17(1) :23- 39 ,

Moler, J. Little, and S. Bangert, user 's Guide, Mathworks Inc., 1987.

MATLAB

1150