An Optimal Reduction Algorithm for Discrete-Time Systems

Copyright © IFAC 12th Triennial World Congress, Sydney, Australia, 1993

AN OPTIMAL REDUCTION ALGORITHM FOR DISCRETE-TIME SYSTEMS Dingyii Xuet and D.P. Atherton School of Engineering and Applied Sciences. The University of Sussex. Falmer. Brighton BN/9QT, UK

Abstract: In this paper, an optimal model reduction algorithm for linear discrete-time systems is presented. The optimisation criterion is introduced to minimise the weighted mean squared error between the original model and the reduced order one when both of the models are subjected to the same input signal. The accuracy of the algorithm is compared with conventional methods and the advantages of the method are readily seen . Key Words: Model reduction; System order-reduction ; Discrete-time systems; Control system design ; Optimisation

1.

and the reduced model Hr/m(z) is given by

INTRODUCTION

Model reduction techniques play a very important part in the analysis and design of control systems (Chen and Shieh, 1968) . In a recent paper, Xue and Atherton (1991a) introduced an algorithm for finding the optimal reduced order models for SISO linear continuous systems. There has been a lot of work done on finding reduced order models for continuous systems, but much less work has been done on finding the reduced order model of a discretetime system. Shih and Wu (1978) gave a reduction method for discrete-time models using the continuedfraction technique which somet.imes can not, preserve the stability of the original system. Therapos (1984) proposed a method using the discrete stabili ty equation technique to solve the problem . The method proposed by Warwick (1984) employed the concept of 'error polynomial' to match the Markov parameters of the system. The approaches given by Hamidi-Hashemi and Leondes (1988), and Hwang et al (1992) employed a bilinear transformation which preserves the Routh denominator and the coefficients of the numerator were optimised for certain criterion .

(z) _ nlz r

H r/m

-

=

The diagram defining the error e( i) at the sampling instants iT is shown in Fig. 1. The Z-transform of the error function e(i) can be evaluated from

E(z)

(1)

r(~ Sampling period T

Figure 1: Block diagram of reduction error The performance index for the optimal model reduction is defined as 00

J = min[I:w2(i)e2(i,B)] =

CRITERIA FOR OPTIMAL REDUCTION

=

00

min[I:f2(i)]

(2)

;=0

where f( i) w( i)e( i, B), and w( i) is a weighting function . Here we explicitly write the error function e(i,B) as a function of the parameter array B. Assuming that E(z) is the Z-transform of the error signal e( i), then for different weighting functions w( i), the Z-transfer function :F( z) can be evaluated as a function of E( z) .

Assume that the original model H(z) is given by

+ 92 Z n-l + •.• + 9n+l + !2zn-l + ... + fn+l

Zn

zn

Hr/m(z)] R(z)

where R(z) is the Z-transform of the input signal

;=0

= 91

= [H(z) -

r(t) .

Two illustrative examples are given to show the use of the algorithm . The reduction results are compared with the reduced models by other approaches and considerable advantages are shown.

H(z)

+ 112 Z r - l + . . . + nr+l + d 2 Z m - 1 + . . . + d m+l

where for simplicity, we assume that d 1 1. The parameter array B to be optimised is defined as

In the method proposed in this paper, the coefficients of the numerator and denominator are optimised . It also allows the reduction of systems which are not strictly proper, which is sometimes not possible using some of the reduction methods available for the continuous case.

2.

d lZ m

t Pennanent address : Dept. of Automatic Control, NOI·theast University of Technology, Shenyang 110006, P.R.China.

821

It is known (Astrom, 1970) that for an equivalent transfer function .1"(z) B(z)/A(z), where .1"(z) is the Z-transform of I( i), such that .1"( z) is stable (ie., A(z) has all its zeros within a unit circle) the performance index J can be evaluated from

~ 12(i) ~

= _1_ 1 B(z)B(z-l) dz 27rj !r A(z)A(Z-l) z

(3)

fr

1

unbiased

(iT)P

E(z) zPdPE(z)/dz P

o-i

E(o--lz)

low

1 - o-i

E(z) - E(o--I Z )

high

i=O,I, . .. ,k

,61:+1 CLI.:+1 ,I:+I-i

1 N b2 . J=-~~ aN 0

,

L.J ,=0

ai ' 0

(7)

Powell's optimisation algorithm can be used to perform optimal order reduction. The procedure for optimisation is as follows:

Table 1: Evaluation 01.1"( z) Freq. range

o-k+IOk+I,I:+I-i

with o-k+1 =elk+1 ,l:+d £11:+1 ,0, ,61.:+1 = bl:+ l ,k+1/ ak+I,O, aN,i = ai, and bN ,i = bi . The integral can then be evaluated from

where denotes the integral along the unit circle in the Z-plane in the counterclockwise direction. The functions .1"(z) for several commonly used weighting functions are given in Table 1.

.1"(z)

= ak+1,i bk,i = bl:+l,i -

elk,i {

""

w(i)

0,1, . . . , k - 1), can be

where elk,i, and b/;:,i, (i evaluated from

=

l. Select an initial reduced model H~/m(z).

2. 3. 4. 5.

low

Obtain the error function E(z) using Eq(l) Calculate.1"( z) for the chosen w( i) . Evaluate the integral in Eq(3). Use Powell's algorithm to iterate one step for a better estimated model H,~/m (z) .

=

6. Set H~/m(z) H;/m(z), go to item 2 until an optimal reduced model H,~/m(z) is obtained.

A suitable choice of the input function ,·(t) is

r(t)

=

Cl

+ C2e-C31

It is well known that numerical optimisation algorithms such as Powell's method, sometimes fail to converge to a global optimum point if the initial conditions are improperly chosen . Thus, a reasonable initial choice of a reduced model H~/m(z) is sometimes very important. It can ensure the optimisation process converges and/or speed up the calculation process. The selection of a good initial reduced model depends on the experience of the user, but some suggestions are :

(4)

where Cl, C2, and C3 are given constants. The Ztransform of this input signal r(t) is

It is worth mentioning that. the input ,·(t) can be simplified to a step input or impulse input when the coefficients Cl, C2, and C3 are accordingly chosen . It is important to note t.hat when there exists a component of a step function in the input signal r(t), ie., Cl "I 0, in order to ensure there is no steady-state error between the original model and the reduced model, one of the paramet.ers, nr+l, then in the reduced order model Hr/m(z) can be evaluated from ",n+l

",m+1

_ L.Ji=l Ui L.Ji-l

n r+l -

n+l

d

i

The user may transform the original model H (z) into the corresponding continuous one for a sampling interval T, then find the stable reduced order model Gr/m(s) using the available methods for continuous systems (Xue and Atherton 1991b) and then transform G r / m (s) back to the discrete- time case to give Hr/m(z). The reduced order model Hr/m(z) may then be l'egarded as an approximation to the original discrete-time system model. The user can use this model as an initial model to find the optimal model H;/m(z) for the original model.

r

-

L:i=l li

~

L.J ni

(6)

i=1

and the parameter array B is reduced to

4.

3.

In this section , two examples are given to show the uses and applications of the optimal reduction method for a high order linear SISO discrete-time system. In the first example, the reduced models obtained by other reduction methods are also given and compared. In the second one, an optimal second order model is obtained for a high order model. Using the reduced model , a PlO controller is designed which can be applied successfully to the original plant model.

THE OPTIMAL REDUCTION ALGORITHM

Once the error function E(z) is established, from Table 1, the Z-transfer function .1"( z) can be obtained and can be written as .1"(z) B(z)/A(z), where B(z) L:~obizi, and A(z) L:~oaizi . If A(z) has all its zeros within a unit circle, then the system is a stable one and the weighted mean squared error can be evaluated recursively (Astrom, 1970) using the algorithm summarised below:

=

=

=

[Example 1J Consider the fourth-order discretetime model given by vVarwick (1984)

Introduce the polynomials k

Ak(Z)

= Lel/;: ,iZi, i=O

ILLUSTRATIVE EXAMPLES

k

and B/;:(z)

= L\bl:,iZi

H(z)

i=O

822

O.3124 z3-0 .5743z 2 +O .3879z-0.0889 z4-3 .233z3+3. 9S69z 2 -2 .2209z+0 .4723

The four poles of the model are located at 0 .8583 ± jO.1887, 0.7582±jO.1915, and the three zeros are at 0.4893, 0.6745 ± jO .3558, respectively.

9r---~--~--~-----------------. 8

7

The reduced order model using the continued fraction method is

HC(z) =

0.1l36z - 0.2269

Original system response ....... Response of H'(z) I• . Response of H'(z) 16 .. ... Response of H«z) .......- Response of HW(z) ....... Response of H'(z) 1•

z2 - 1.7207z + 0.7693 The reduced order model given by Warwick (1984) is

H

W( ) z

0.3124z - 0 .0298 = -;;--7-:::-::-:-::--~=::: z2 - 1.7369z + 0 .7773 5

The optimal Routh reduced order model, using the method given by Hamidi-Hashemi and Leondes (1988) for step input, is

Hr(z)1

= u

0.7082z - 0.4946 z2 _ 1.7213z + 0.7578

0.1299z + 0.1820 H·(z)lu = 2 z - 1.7431z + 0.7877 6

0.8 ~

§.

0.6 ......

~ 0.4

..s&.

0.2327z + 0 .0431 z2 _ 1.7721z + 0.8115

H(z)

00

49 .7177

HC(z)

0.0595

1.2897

HW(z)

0.1463

1.2266

Hr(z)lu

l.l821

8.6061

H·(z)lu

0 .0304

0.7856

H·(z)16

0.0884

0.4715

Original system resp . ....... Response of H'(z) I. ...... ..... Response of H'(z) 16 . . ... Response of HC(z) _ Response of HW(z) \'::'. ....... Response of H'(z) 1

0.2

.

"',

.0.20~---5':-----:'1'=-0---1'::-5----=-2'=-0------='25 Time (Sec)

Figure 3: Impulse response plots

[Example 2] A digital PID control design application A digital PID control system is shown in Fig. 4 . It can be shown (Kuo 1980) that the Z-transform of a digital PID controller Gc(z) can be written as

Table 2: Performance indices of reduced models Impulse input

35 40 Time (Sec)

o~-----------~~~~--~

The performance indices for step and impulse inputs are given in Table 2 . Note that, in the table, the values for H(z) are the mean squared values of the output signal y(i) for the original model H(z) driven by the corresponding input and the rest of the values are the mean squared errors obtained from Eq(2), with w( i) = 1. Also the optimal reduced order

Step input

30

....~.....: :-: .: :~.

~

and

model

25

Figure 2 : Step response plots

Selecting a weighting function w( i) = 1, the optimal reduced models for step and impulse input signals using the algorithm given in this paper are given by

H.(z)1 =

20

IS

10

Gc(z)

= koz 2 + k1z + k2 2Tz(z - 1)

where ko=K/T2+2Kd-2Kp T, kl=K/T2_2Kp T-4Kd, k2=2Kd, and T is the sampling period of the system .

model H·(z)16 uses an impulse function as the input signal, and Eq(6) has been used to determine nr+l . The step and impulse responses of the reduced order models are shown in Figs . 2, and 3, respectively, together with the responses of the original model. It can be seen from Table 2 and the figures that the step and impulse responses of the optimal reduced order models are very close t.o those of the original system . It is very difficult to define in general the accuracy for model reduction, because different users have different views of what is 'optimal'. It can be seen that the optimal reduced order models from the method given in this paper are better than the other reduced order models, if the criterion is closeness of the matching of the output response .

Figure 4: A digital PID control system For a second order plant model

H

2

() _ nOz +nlz+n2 2/2 z z2 + d1z + d 2

the integral gain of the controller can be found using the method given by Kuo (1980) as J{ /

J{v(l + d 1 + d2 ) = --'-----...:..:.

no

823

+ 111 + 112

(8)

where /{v is a defined ramp-error constant. Assuming that the numerator of Ge(z) can cancel the denominator of H 2/ 2(Z), then two independent linear algebraic equations can be set up for the unknowns /(p, and /{d, namely 2T( l+d~ )Kp +2(2+d 1 )I~d=K:(1-dl:T2 { -2d 2TI\p + 2(1- d 2 )l\d 1\[d 2T-

Table 3: Step response 11

(9)

=

which can easily be solved for

and

/(p

/(d.

Consider the original plant model H(z) given by (Hamidi-Hashimi, and Leondes, 1988), where H

2Z4 + 1.84z3+0 .64z2+0.008z-0 .096 (z) z4-1.2203z3 + 0.2257z 2 +O .1166z+0.0249

Choosing a sampling period of T = 0.15, then the optimal reduced model is

• (z)1 _ 2.1855z 2 +O.9051z+0.4482 H 2/2 6 z2 _ 1.4878z + 0.6060

kT

y( kT)

0 0.15 0.3 0.45 0.6 0.75 0.9 l.05 l.2 l.35 l.5 l.65

0.5264 0.8790 1.020 1.0182 1.0025 0.9987 0.9994 1.0000 l.0001 1.0000 1.0000 l.0000

I

y( kT) 0.5043 0.9174 1.018 l.0027 0.9940 l.0050 1.0042 l.00 16 l.0005 l.0001 0.9995 0.9992

I

e( kT) 0.0221 -0.0384 0.0020 0.0155 0.0085 -0 .0063 -0.0048 -0.0016 -0.0004 -0 .0000 0.0005 0.0008

0.042 -0.044 0.002 0.002 0.008 -0 .006 -0.005 -0.002 -0 .001 -0 .000 0.000 0.001

in this paper and the result.s are compared with the other methods from which it can be seen that the accuracy of the method is better than the other methods. In the second example, a digital PID controller is designed from the parameters of the second order reduced model and applied to the original system. The control result can be considered satisfactory from an engineering viewpoint. The MATLAB based model reduction routine, Xue and Atherton (1991b), has been extended so that the discrete-time models can also be reduced using different reduction methods .

An impulse input signal was used for the modelling since this is a reasonable approximation for the actual error signal c(i) for the input to the plant model in the closed loop . The Nyquist diagrams of the corresponding continuous models of H(z) and H;/2( z) 16, transformed using a zero-order hold, are shown in Fig. 5, from which it can be seen that they are very close.

0

6. REFERENCES -5

-Original system ----- Response of H~ (z) I 6

-10

Astrom, K. J., (1970) . Iniroduction io siochasiic conirol theory, Academic Press. Chen, C. F., and Shieh, L. S., (1968). "A novel approach to linear model simplification" ,/nt. J. Contr., 8, .561-570. Hamidi-Hashemi, H., and Leondes, C. T ., (1988). "Reduction of high-order discrete-time systerns", Int. J. Sysi. Sci ., 19, 1883-1890. Hwang, C., Guo, T.-Y., and Shieh, L.-S., (1992). "Model reduction using new optimal Routh approximation technique", Int . 1. Contr., 55, 989-1007. Kuo, B. C. , (1980). Digital control systems, HRW series in electrical and computer engineering, Holt, Rinehart and \Vinston , Inc . Shih, Y.-P , and \Vu, W.-T. , (1973). "Simplification of z-transfer functions by continued fractions", Int. J. Contr. , 17 , 1089-1094. Therapos, C . P., (1984). "Low order modelling via discrete stability equations", Proc. lEE, D-131, 1811-1813. Warwick, K., (1984). "A new approach to reducedorder modelling" , Proc. lEE, D-131, 74-78. Xue, D., and Athert.on, D. P., (1991a). "An optimal model reduction algorithm for linear systems", Proc . A merican Control Conference, 2, 21282129, Boston , USA. "A Xue, D., and Atherton, D. P., (1991b) . menu-driven model reduction program and its applications" , PrepT'ints of 5th IFAC Symp . on CAD in c07ltT'ol.systems, 316-321 , Swansea, U.K.

-15 -20 -25 -30 -5

0

5

\0

15

20

25

30

Figure 5: Nyquist diagrams comparison Selecting f{v = 12, the controller parameters can be found from Eqs (8), and (9) as f{p = 0.1706, f{[ = 0.4014, and f{d = 0.0463. Applying the PID controller to the original high order plant model, the step response data of the closed loop system is shown in Table. 3, where y(kT), and y(kT) are the outputs of the reduced system and original one, respectively, and e( kT) = y( kT) - y( kT). It can be seen that the output signal y( kT) is very similar for control of the original and reduced order models. 5. CONCLUSION

In this paper, an optimal reduction approach for discrete-time linear systems has been proposed and examined . The fundamental idea of this approach is that the weighted integral squared error between the output of the original system and that of the reduced order one is minimised when both of them are subjected to the same input signal. Two illustrative examples are presented to show the approach given

824