Identification and model reduction of multivariable continuous systems via a block-pulse functions scheme

Identification and model reduction ofmultivariable continuous systems via a block-pulse functions scheme H o n g Y. Z h a n g a n d L e a n g S. S h i e h

Department of Electrical Engineering, University of Houston, Houston, Texas 77004, USA R o b e r t E. Y a t e s

Guidance and Control Directorate, US Army Missile Command, Redstone Arsenal, Alabama 35809, USA {Received May 1981) This article deals with the identification and model reduction of multiinput and multi-output (MIMO) continuous systems, described by a highdegree vector differential equation or a matrix transfer function, from the discretized input and output data. The approach is based on using the block-pulse functions scheme and the least squares method. A two-shaft gas turbine model is used as an illustration. Key words: multivariable control systems, identification, model reduction Introduction Accurate description of most practical continuous systems often results in a set of high-degree coupled-multivariable differential equations. The linear representation of these systems is a high-degree vector differential equation or a matrix transfer function as follows:

p+l

p

2 AiDi-lx(t)= 2 i=1

BiDi-lu(t)

(la)

i=1

with initial vectors D i - I x ( 0 ) and D i - lu(0), i = 1, 2 . . . . . where x(t) is an n x 1 output vector and u(t) is an n x 1 input vector. A i and B i are n x n constant matrix coefficients and D = d]dt is the differential operator. When all initial vectors are null vectors, the corresponding polynomial matrix description 1 of a MIMO system is:

p+l

p

Z Ai si-lX(s)= ~ Bi si-tU(s) i=1

(lb)

i=1

If the system of interest is a single-input and single-output (SISO) continuous system, or n = 1, equation (1) becomes a pth degree scalar differential equation. Chen and Hsiao, 2 Rao and Sivakumar a and Gopalsami and Deekshatulu a have successfully applied the Walsh functions approach, and blockpulse functions approach for the identification o f the scalar differential equation or the corresponding (scalar) transfer function from the discretized input-output data. Recently, Rao and Sivakumar s have extended their Walsh functions approach a to identify a transfer function matrix but not a matrix transfer function of a MIMO system. In this paper, an extension of the methods of Chen and Hsiao 2 and Rao and Sivakumar a to the multivariable case is carried out by using the combination of the block-pulse functions scheme and the least-squares method. The MIMO systems of interest are formulated in a high-degree vector differential equation or a matrix transfer function. The same approach is then extended to determine the reduceddegree continuous model of the original high-degree continuous system by using the input-output data.

and the matrix transfer function becomes: Model identification

X(S) = [p+l L i~l

Aisi_l] -I [ Lp

Bisi-I ] U(s)

0307-904X/82/050369-05/$03.00 © 1982 Butterworth & Co. (Publishers) Ltd

(lc)

Consider a noise free linear time-invariant system described by a pth-degree vector differential equation in equation

Appl. Math. Modelling, 1982, Vol. 6, October 369

Multivariable continuous systems: H. Y. Zhang et aL (la) with all zero initial vectors. Let Ap+] in equation (la) be In, an n x n identity matrix. When m sets of discretized data (or sampled data) ofx(t) and u(t) are available, x(t) and u(t) can be described, using the block-pulse functions6' 7 as follows: x(t) = [x,(t), x2(t) . . . . . xn(t)]' = Cqb(t)

(2)

u(t) = [u,(t), u2(t) ..... u,(t)]' -~ Dck(t)

(3)

t

The prime in equations (2) and (3) designates transpose. ~(t) contains m block-pulse functions ~i(t) defined as: ¢(t) = [¢l(t), ~bu(t). . . . . ¢,n(t)]'

and

jT

'f

dii= T

ui(t) dt~-½[ui(jT)+ui(]T+ T)] (6b)

(i-Or

i = 1 , 2 . . . . ,n

j=1,2 ..... m

The integration of the set of block-pulse functions ~b(t) can be expressed as: t

f

(4a)

where:

¢(t) dt mHck(t)

(7a)

0

¢i(t)=

1, jT
(4b)

0, otherwise

(4c)

where the m x m operational matrix H is:

T in equation (4) is the sampling period or T ~=t:/m where t.t-is the final time of interest and m is the number of discrete data used. The constant matrices C and D in equations (2) and (3) are: rCll

C12 • • •

Clt~

C = I c21 . . . c2~ . . . ... C2m}

LCnl

(5a)

H= T

The waveforms of the block-pulse functions ¢i(t) for tf= 1 and m = 4 and their integrations are shown in Figure 1. Integrating equation (la)p times from 0 to t and using the relationship in equation (7) yields: p+l

p

i=1

i=1

~" AiCHp+l-i~(t)= ~. BiOHP+l-i~(t)

and d12 . . .

dim 1

... d2m]

D= II2,. . . d22 . . .

(5b)

p

AiCHV+l-i= ~ BiDH p+l-i

dn2 ... dnnO

i=1

The elements of C and D can be determined by using the orthogonal property of the block-pulse functions expressed as follows: l

ci/

T

(9)

i=1

Rearranging equation (9), we have the following matrix equation: C = QZ

iT

I

(8)

Dropping the common vector ¢(t) in both sides of equation (8) gives: p+l

Ld,,1

(7b)

0

cn2 ... Cnmd

dll

½

xi(t) dt~-½[xi(jT)+xi(jT+ 7")1

(6a)

(10a)

where:

Q = [--A~,--A2 .... , - A p , BhB2 ..... Bp]

U-I)T

(10b)

~ ,./'qbi d/"

I

0

I/4

I t I

I

0

~fCzd t

I 0

I/4

1/2

I t I

1/41

I

0

I/4

,

0

1t¢4 I 0

I/2

I

3/4

I

7q

3/4

I

0

0

Figure 1 Block-pulse functions and their integrations

370 Appl. Math. Modelling, 1982, Vol. 6, October

-•

I

I

I

0

-~

I

1

~2

f ~b3dt

0

0

½

1

¢3

fq~4 dt_

0

0

0

½

f~bl dt

~ 1 I T

I

/

it

I/2

I

5/4

I

f~b2dt

=,

Multivariable continuous systems: 14. Y. Zhang et aL Z ' = [(CHP) ', (CHP-i) ', . . . . (CH)', (DHP) ', (DHP-I) ', . . . . (DH)']' (10c)

where:

where: Q E R n x 2 p n , z E R 2pnx'n and C E R nxm. Normally, the number of data m is greater than the total number of unknown parameters 2pn. The weighted leastsquares solution of Q (defined as Q+) that minimized the sum of the squares of residual J, where:

J = ( C - QZ) W ( C - Q Z ) '

(11)

is given by:

Q+ = CWZ'(ZWZ') -l i f r a n k Z = 2pn

(12)

where W is a symmetric m x m positive definite weighting matrix. The existence of the solution in equation (12) heavily depends on the selection of input functions ui(t ) in equation (1). For example, if the functions ul(t) and u2(t) are identical, then from equation (6) we have dli = d2i for ] = 1, 2 . . . . . m. As a result, Z matrix will have two linear dependent rows and rank Z < 2pn. Thus, ZWZ' in equation (12) becomes singular and Q÷ cannot be uniquely determined. The accuracy of the identified parameter matrix Q in equation (10) can be improved upon by selecting binary signals as exciting signals u(t) in equation (1). Since the binary signals are constant in the sampling intervals, the elements dii in equation (6b) are exact values instead of approximate ones. The pseudo-Random Binary Signals (PRBS) are ideal inputs for the proposed identification process. Instead of directly determining the matrix inverse in equation (12) to solve Q, an effective and accurate numerical method called the singular value decomposition method, a is recommended to solve the Q in the normal equation of equation (10a). Note that although discrete data of input u(t) and output x(t) are used to identify the unknown parameters, the obtained parameters are those of a continuous system, not those of a discrete system.

(~= [--A,B]

(14b)

Z' = [(CH)', (DH)']'

(14o)

C, D and H are shown in equations (5) and (7), respectively. The weighted least-squares solution (defined as Q÷) becomes:

(~÷= CWZ'(ZWZ')-1 if rank Z = 2n

(15)

If the static matrix gain (defined as an t7 x n constant matrix K) can be obtained from the available mathematic formulation of a MIMO system or from the experimental data at steady state, then the steady-state values of the unit-step response of the original high-degree system and the reduced-degree model can be exactly matched as follows:

K ~=,4-1B = A-I1BI if rank,4 = rank A 1 = n

(16a)

/7 = .4K

(16b)

or

K in equation (16a) is formulated from equations (13) and (la). Substituting equation (16b) into (14) we have:

c : t-A, AKi

LDH d

:`4 t- c + K D I -

=A~P

(17a)

where:

2 ~= [ - - C + K D ] H

(17b)

The weighted least-squares solution of equation (17) is:

A = CW2'(2WZ) -1 if rankZ = 21l

(18)

and the parameter/7 can be determined from equation (16b). Thus, the reduced-degree model in equation (13) has the same steady-state value of the original high-degree system.

Model reduction A low-degree model of a high-degree MIMO system is often desirable for computer simulation and control system design. In recent years, various model-reduction methods have been proposed and compared 9-1a in both time domain, and frequency domain. In this paper, the proposed identification procedure is extended to determine the reduced-degree model of the original high-degree system by using the discrete data generated from the original high-degree continuous system. Suppose the original system is described by equation (1) with p > 1 and the desirable reduced-degree model is assumed to be a first degree vector differential equation:

2 ( 0 + Ax(t) = Bu(t)

(13a)

2 ( 0 + A22(t) + A1X(t) = Bdt(t)

(19a)

x(O) =.~(0) = 02x,

(19b)

where: [11.8567

-262.21]

A2 = L0.

101.368.1

= r 18.5671 Al

L-0.005214

'

-2622.1

]

136.829.1

[x2(t)J'

,,,:,>:r",<'>l Luz(t)J

(13b)

where A and/7 are unknown n x n constant matrices to be determined by using the input-output data generated from equation (1). Integrating equation (13a) once and using the peculiar property of the block-pulse functions shown in equation (7) yields the following algebraic matrix equation: C = OZ

Consider the following linearized two-shaft gas turbine MIMO system: J2, xa

7]

where the initial vector, x(0) is a null vector. The corresponding matrix transfer function is:

X(s) = (si n + d )-'/Tu (s)

An illustrative e x a m p l e

(14a)

It is desirable to identify a second-degree MIMO model in equation (1 a) and obtain a first-degree reduced-degree model in equation (13a) by using the input-output data generated from equation (19). If the inputs u(t) are chosen as the PRBS, the simulations of equation (19) can be easily conducted by using the equivalent discrete model of the continuous model of equation (19). The equivalent discrete model is determined by using a geometric-series method 14

Appl. Math. Modelling, 1982, Vol. 6, October 371

Multivariable continuous systems: H. Y. Zhang et aL

as follows. The state-space equation of equation (19) is:

)'~(t) = Ay(t) + Bu(t)

(20a)

y(0) = 04x I

(20b)

a

,.O:gS,Jo.ge; o

d

i i I

where:

A=

Reduced-degree model

--A1

~

BI ,

--A2

-1,0--

-2,0

y(t) = Oq,Y2,Ys, Y,]' ~ [xl,x2,£c,,2=]' -3.0

I 25

The eigenvalues of this system are -- 1.341367, -- 1.883743, --10 and --100. The equivalent discrete model is:

y(k + 1) = Oy(k) + Hu(k)

(21a)

y ( 0 ) = 04x1

(21b)

I 50

I 125

I 150

I 175

200

150

175

200

006 0riginol system and Identified model

0.03 oJ x

where:

~-

1

-7/(AT) + 16/---= (AT;

1_2__(At)2]}j

1

X [I4 + ~/.(AT) + 1612

H ~- (G --I4)A-1B

8 (21c)

o

-003

-0.06

25

50

75

2.0302418E-02

3.6594978E-08

9.9903966E-01

7.2521597E-02

9.9918938E 00

1.000266E 00

1.7138524E-05

-4.4977069E-01

B~= L-5.793261E-05

3.9064349E-03

1.8356973E-03]

4.9796707E-11

3.2871977E-03 /

9.5353638E-01

8.5852776E-01 /

3.6004553E-08

6.6582302E-01J

125

Figure 2 PRBS responses of various systems

9.9985380E-01

I_

I00

Number of samples

(21d)

When/" = 4 in equation (21c) and the sampling period T = 0.004 s, the G and H in equation (21) become:

G=

I I00

Number of samples

b

G

I 75

and

(21e)

-7.044263E-03] 1.002301E 00.1

Comparing the parameters in equation (19) with those in (22), we obtain very satisfactory results. In order to measure the goodness of the approximation, a criterion Ji is defined as follows:

J ~ [xiqT)

and

-- Xict(jT)]

2

/=1

7.8748197E-06

2.5298540E-06 ~

5.0899795E-14

7.0187807E-06

H=[3.9064349E-03

1.8369732E-03

L4.9796707E-11

3.2871977E-03

Ji = (210

The simulation results of equation (21), using PRBS signals as the exciting signals ul(k) and u2(k) are shown in Figure 2. The ul(k) and u2(k) are generated from different registers of the PRBS generator. The identified second-degree vector differential equation is:

fc(t) + A*=f~(t) + A+lx(t) = B~u(t)

(22a)

:~(0) = ~ (0) = O=x ]

(22b)

m

i = 1, 2

(23)

E

j=l where xi(t ) is the response of the original system, whereas xia is the response of the identified system. The computed criterion with m = 200 for x,(t) is 0.1342486% and x2(t) is 0.281794%. The PRBS responses of the original system and the identified model are shown in Figure 2. The results are quite satisfactory. The reduced degree model of equation (19) with the static matrix gain in equation (16):

K=(A1)_,BI=[5.4150121E-02 L2.0635025E-06

where:

1.0377272E 00] 7.3481507E-03-1 (24)

can be determined from the proposed method as follows:

A~-

372

-1.438191E-02

1.358917E 02

Appl. Math. Modelling, 1982, Vol. 6, October

Yc*(t) + .4"1x*(t) = B~(u(t)

(25a)

x*(0) = 02x 1

(25b)

Multivariable continuous systems: H. Y. Zhang et al. 12

a

0.9

06 Originol system

0 0.3

/"

0

I 15

I :50

i I I 45 60 75 Number of somples ~10'

I 90

I 105

120

pn x n input matrix, p n ( p n + n) parameters are often required to be identified. However, by using the proposed method, only 2pn parameters are required to be determined. Thus the computational difficulty in the identification can be significantly reduced. Also, the identified parameters are the parameters o f a continuous MIMO system but not those of a discrete MIMO system although discretized input-output data are used. The same method has been employed to construct a reduced-degree model with retention of the static matrix gain o f the original high-degree system. A two-shaft gas turbine model has been used as an example.

012 Acknowledgements

b Reduced-degree model

This work was supported in part by the US Army Research Office, under Grant DAAG 29-80-K-0077 and US Army Missile Research and Development Command, under Contract DAAH01-82-C-A 139.

0.09 oJ x ~.

Q06

d o.o3 References [

00

15

I

30

i

45

I

60

I

75

i

90

I

105

120

Number of somples *10'

Figure 3 Unit-step

responses of original and reduced-degree models

1 2 3 4

where:

5

~ 4 ] = [2.065324E 00

-3.011279E 02]

L4.745973E-03

5.296944E-01J

/~].=[

1.112162E-01 -6.949005E-02

2.580880E-04]

6 7

8.817229E-03 J

The time responses o f x * ( t ) and x * ( t ) having PRBS inputs are shown in Figure 2. The computed criterion with m = 200 for the x * ( t ) is 9.259323% and x * ( t ) is 13.131070%. Also, the unit-step responses o f the original system and the reduced-degree model are compared in Fig~tre 3. The computed criterion with m = 1200 for x * ( t ) is 7.276602% and x * ( t ) is 6.179920%. Again, the results are very satisfactory.

Conclusion A method based on the block-pulse functions has been presented for the identification and model reduction of a MIMO system. For a system described by a first-degree state- space equation with pn x pn system matrix and

8 9 10 11 12 13 14

Kailath, T. 'Linear systems', Prentice-Hall, Englewood Cliffs, 1980 Chen, C. F. and Hsiao, C. H. 'Time-domain synthesis via Walsh functions', Proc. lEE, 1975, 122,565 Rag, G. P. and Sivakumar, L. 'System identification via Walsh functions', Proc. IEE, 1975, 122, 1160 Gopaisami, N. and Deekshatula, B. L. 'Time-domain synthesis via Walsh functions', Proc. lEE, 1976, 123,461 Rag, G. P. and Sivakumar, L. 'Transfer function matrix identification in MIMO systems via Walsh functions', Proc. 1EEE, 1981,69 (4), 465 Chen, C. F. et al. 'Walsh operational matrix for fractional calculus and their application to distributed systems', J. Franklin Inst. 1977, 303, 267 Shieh, L. S. et at 'Solution of state-space equations via blockpulse functions', hzt. J. Control, 1978, 28 (3), 383 Forsythe, G. E. et al. 'Computer methods for mathematical computation', Prentice-Hall, Englewood Cliffs, 1979 Bosely, M. J. and Lees, F. P. 'A survey of simple transfer function derivation from high order state variable models', Automatica, 1972, 8 (6), 765 Power, H. M. 'A review of method reduction',lEE Control Automation Div. Colloquium on: reduced order models and their use in design of dynamic systems, I]1-I/7, 1975 Decoster, M. and Van Cauwenberghe, A. R. 'A comparative study of different reduction methods', J. A, 1976, 17 (2), 68 and (3), 125 Muller, G. S. 'Linear model of 2-shaft turbojet and its properties', Proc. IEE, 1971, 118, 813 Wei, Y. J. and Skieh, L. S. 'Synthesis of optimal block controllers for multivariable control systems and its inverse optimalcontrol problem', Proc. IEE, 1979, 126 (5), 449 Shieh, L. S. et al. 'Discrete-continuous model conversion', Appl. Math. Modelling, 1980, 4,449

A p p l . Math. Modelling, 1982, V o l . 6, O c t o b e r

373