- Email: [email protected]
Optimal model reduction by multipoint Padé approximation
Optimal Model Reduction by Multipoint
by T. N. LUCAS Department of Mathematical and Computer Sciences, Dundee Institute of Technology, Bell Street, Dundee DDI IHG,
U.K.
A novel method of obtaining optimal reduced-order models of linear system transfer functions is presented. It uses the popular multipoint Pad& approximation technique in am iterative way to generate eflciently the optimal models. Central to the method is a new way of calculating Pad& approximants about many points by reducing them to equivalent Taylor series approximants. Optimal reduced-order models for impulse and step inputs are considered, and it is seen how the method may be extended to ramp and other polynomial inputs. Numerical examples are given to demonstrate the method.
ABSTRACT:
I. Introduction Reducing the order of linear time-invariant systems has continued to be a topic of interest in control system design over the past couple of decades or so. Its importance stems from the ever increasing complexity of modern control system models, which tend to be of high-order when linearized about convenient operating points. The use of reduced-order plant models to design an adequate controller is now an accepted practice among control engineers and, consequently, an array of different approaches to the “order reduction” problem has been well documented in the literature, e.g. (l-3). single-output systems the transfer function description For single-input, (frequency-domain) tends to find favour over the state-space description (timedomain) which is, however, more suited to multi-input, multi-output systems than the transfer function. Both approaches to the model reduction problem are equally important and it is found that the two approaches are not entirely divorced from each other, where many examples of links between time- and frequency-domain methods have been cited, e.g. (4, 5). The principal interest of this paper is in the frequency-domain approach, and, in particular, that of multipoint Pad&/continued-fraction approximation methods (68). These are known to be very flexible and computationally efficient methods of reduction and as such have proved to be popular techniques to use. However,
The Frmkhn
hmtutc
0016
0032;93 $j.OO+O
00
79
T. N. their biggest drawback is that they cannot guarantee stability of the reduced-order model if the original system is stable, unless a subset of the system poles are used as expansion points to ensure that these poles are retained in the reduced model (799). To overcome this stability problem, a number of methods have been developed alongside the Pad& methods, ensuring stability of the reduced models. These are based mostly on the Routh table properties utilised by Hutton and Friedland (10) and Chen et al. (11). Such methods have proved to be popular because of their relative ease of computation along with the stability preserving properties. However, it is generally recognised that the reduced-order models produced by these methods tend not to be such good approximations as those of the multipoint Pad& methods when they produce stable models. This is due to twice the number of system parameters being retained in the Pade models compared to those of the “stability preserving” methods, in which the reduced denominator polynomials are derived before calculating the numerators via matching appropriate system parameters (time moments, Markov parameters or frequency points). More recently, interest has developed in finding “optimal” and “sub-optimal” reduced models in the frequency-domain by minimizing the integral-square-error :
I=
sOx
[y(t) -vr(O12dt>
where y(t) and y,(t) are the full and reduced-order models’ impulse (or step) responses, respectively. The sub-optimal transfer function models are usually derived first by finding a reduced denominator polynomial by one of the stability methods, and then choosing the numerator to minimize 1, e.g. (12-14). Such models require a modest amount of computation, having eventually to solve a system of linear equations for the numerator coefficients. Various commendable approaches for finding the highly desirable optimal reduced-order transfer function models have been given by, among others, Luus (15), Vepa (16), and Bryston and Carrier (17). However, the relatively little use made of these methods by control engineers appears to be due to the quite substantial amount of computation involved compared to the other methods already mentioned. In general, these optimal methods require solving sets of nonlinear equations by classical numerical techniques such as the generalized NewtonRaphson method. This paper presents a novel method for obtaining optimal reduced-order models by the application of multipoint Pad& approximation. Tt removes the need for classical numerical methods, and instead uses the efficient PadC technique in an iterative way to converge onto the optimal reduced-order model. Both impulse and step response models may be obtained by the method in a straightforward manner, requiring only the solution of linear sets of equations. The technique is seen to be possible only because of a new way of calculating multipoint PadC approximants by means of multipoint Taylor approximants of polynomials. This result is discussed in Sections IT, III and IV before its use in optimal model reduction is developed in Section V. Numerical examples are given in Section VI to demonstrate the application of the method. 80
Journal of the Franklm lnst~tute Pergamon Press Ltd
Optimal Model Reduction II. Basis of the Method Consider
the nth order transfer G@)
which is to be reduced
=
function,
in usual notation,
b,-,s-‘+.~~_~
~~..
a,s”+
“.
given by
+h,s+h,
(1)
... +a,s+a,
to the kth order model dA_,skm’+ ... +d,s+d, R(s) = s +ek_,s ‘-I+ .:. +e,s+e,
(2)
by Pad& approximation about the 2k points, s = p,, i = 1, 2,. . . ,2k. The problem is to find the 2k parameters d, and e, (v = 0, 1, 2,. . , k- I), such that G(p,) = R(p,), Using Eqs (l), (2) and (3), it is obvious P(S) = Q(S)
at
i = 1,2,. . . ,2k.
(3)
that (3) is equivalent
s = pl,
to
i = 1,2,. . . ,2k,
(4)
where P(S) = (a,?+
... +a,s+ao)(dk_,sk-‘+
... +d,s+d,)
and Q(S) = (b,_,s”-‘+
... +b,,)(s”+ek_,sk-‘+
... +e,s+e”).
This relationship changes the problem to one of finding the coefficients of two (n+k1)th degree polynomials, P(s) and Q(s), such that their values are equal at the 2k expansion points pi. These coefficients will of course be linear combinations ofthed,ande,(r=0,1,2 ,..., k-l). The central idea of the new method is now to find the unique polynomials, of degree (2k - l), which interpolate P(S) and Q(S) at the 2k expansion points s = pi (i = 1, 2, . . . ,2k). By definition, these polynomials must be identical because P(s) and Q(s) have the same values at the expansion points. Hence the 2k coefficients of both polynomials may be equated to give a set of 2k linear equations to solve for the parameters d, and e,. These equations will, in effect, be the multipoint Pade equations of the system. The reduced degree interpolating polynomials can be derived very easily by a Routh-type algorithm (18) and may be thought of as multipoint Taylor approximants of P(S) and Q(S), respectively, about s = pI. Furthermore, the algorithm is well suited to cope with real, complex and multiple expansion points. Lucas (18) shows that dividing a polynomial T(S) of degree n by another, H(s) of degree m (m < n), from highest powers, gives a remainder polynomial of degree (m- 1) which is the Taylor approximation of T(S) about the m roots of H(s). To illustrate this important concept, consider the simple case of reducing the 3rd-order system : 8?+6s+2 G(s) = ~ ~ s3+4s2+5s+2 81
T. N. Lucas to one of 2nd order, given by d,s+d, R(s) = yp s +e,s+erj by multipoint approximation about the four points s = 0, 1, 2, 3. The divisor polynomial in the algorithm is thus s(s-l)(s-2)(s-3) and the polynomials
- s4-66s3+lls2-66s
P(s) and Q(s), as defined by Eq. (4), are given by
P(s) = d,s4 + (4d, + do)s3 + (5d, +4d,)s2 + (2d, + 5d,,)s+ 2d, and Q(s) = 8s4+(8e,
+6)s3+(6e,
The table for the 3rd-degree given as follows (18) :
Taylor
+8e,+2)s2+(2e, approximant
+6e0)s+2eo.
of P(s) about
s = 0, 1, 2, 3 is
d,
(4d, + do)
(5d, +4&J
(2d, + 5do)
2d,,
1
-6
11
-6
0
(lOd, +d,)
(-6d,
+4do)
2do
(8d, +5do)
where the first row coefficients are those of P(s), the second row are those of the divisor (expansion points) polynomial, and the third row are those of the 3rd degree approximant of P(s) calculated from rows 1 and 2 by the usual Routh division algorithm. Similarly, the approximant of Q(s) is given by 8
(8er +6)
1
-6
(8e,+54) Equating
(6e, +8e,+2) 11
(6e,+8eo-86)
like coefficients
(2e,+6eo+48)
in both approximants
(2e, +6eo)
2eo
-6
0 ’
2e,
gives the four linear equations
lOd, +do = 8e, +54 -6d,
+4do = 6e, +8eo-86
8d, +5d,
= 2e, +6eo+48
2d,, = 2eo which solve to give the reduced do = 5.3488,
model parameters
d, = 7.4884,
as :
e. = 5.3488,
e, = 3.2791.
It is easily verified that this model fits G(s) at all four expansion points. Although the above 2nd-order reduced model was obtained in a straightforward manner, it is easy to see that for larger systems requiring more than a single division step for approximating P(s) and Q(s) (or the reduced model is of order greater than 2), then the Taylor approximant tables can become rather cumbersome and therefore unattractive to use. Consequently, it is desirable to produce a more 82
Journal
of the Franklin Institute Pergamon Press Ltd
Optimul Model Reduction compact form of the algorithm the next section.
for computational
purposes.
This is discussed
in
III. Matrix Form In general, the (n+k) coefficients of P(s) and Q(Y), as defined by Eq. (4), may be written in matrix form as the vectors Ad and B[l e] ‘, respectively, where
a, a,-, an-2
0
a, a,,- 1
. an-k+1
’
.
.
.
. . .
A= a0
.
al
.
a0 .
.
.
.
0
4- I ‘6-z
a, . ,
d=
ak-I
d,
ak-2
do
_
. .
a0
0
6-k-l
’
. .
.
.
.
.
.
bn-,
B=
and ho
b, bo
b,
.’ .
.
bk h-1
.
.
.
e =
0 bo
Suppose that the reduced model is obtained which are the roots of the polynomial m(s) = s2k+~2k_iS2k~‘+
by expanding “’ +~is+&,.
about
the 2k points
(5)
Now, from the elementary theory of transformation matrices, it is seen that the third rows in the multipoint Taylor approximant tables of P(s) and Q(s) (after the first division by m(s)) will be given by the entries in the two vectors M,Ad and Vol 330, No. I, pp. 79-93. 1993 Prmted m Great Bntain
83
T. N. Lucas M, B[ 1 e] ‘, respectively,
where
-XZk-
1 0 0 . .
I
-aZk&Z -a2&3
M,
=
-sIo 0 0
1 0 . .
1 .
0 ’ O....Ol () . . .
0
0 .
.
.
0
(n-f- T-1)
1
. .
. .
0 .
o......
<
1 .
.
.
0
I
,
(n+k)
Similarly, the next rows of each approximant table (after the division by m(s) into the third row polynomial) may be represented by the vectors obtained by premultiplying the above matrix products by the transformation matrix M2, which is formed from M, by deleting its last row and column. Hence, in general, as there will be (n-k) division steps, the coefficients of the 2kth degree Taylor approximants of P(s) and Q(s), respectively, are given by the elements of the vectors
where
-x2h-2
1 0
-@2h-3
0
-a2k
I
1 0
Mi = -cl0
0 0
L -
0
0 1
.
.
.
.
.
.
(n+k-i) .
0
O....Ol .
1 . 0
o..... (n+k+l
(6)
-i)
1A >
fori= 1,2 ,..., n-k. Equating these vectors to solve for d and e, gives Fd = [c
H]
(7)
where F = Mplmk. ..M?M,A and [c H] = M,, k.. M2M,B, with c being the first column of this matrix. Equation (7) may be rearranged to give the linear multipoint Padt equations as x4
Optimal Model Reduction
II
[F-H]
t
=c
which can be directly solved by any standard algorithm for linear equations. Notice that complex expansion points, which always occur in conjugate pairs for physical realisation of the reduced model, cause no problems for the entries in the M, matrices as real polynomial coefficients LX,will always result. Notice also that the ordinary Pad6 approximant (about s = 0 only) is obtained by putting CI, = 0, m = 0, 1, 2,. . . ,2k- 1, in the A4, matrices. The whole process therefore consists of simple matrix multiplications and solving 2k linear equations.
IV. Illustrative Consider
Example for Multipoint
the 4th-order
transfer
G(s) =
Method
function
(8) given by
9s’+42s2+31s+10 s4+8S3+21s2+22s+8’
For expansion points such as s = f 2j, 1 + lOj, most existing multipoint continuedfraction/Pad6 methods would not be able to cope without resorting to complex arithmetic, which renders them unattractive to use (the choice of complex expansion points will be seen to be important in obtaining optimal models later). However, the given method quickly produces a second-order model as follows. From the given expansion points, the polynomial m(s) is given by m(s) = (.~~+4)[(s-1)~+100] and the corresponding
rs4-22S3+105s2-8s+404
transformation
matrices
are
2 M,
=
and
1
0
0
0
0
1
0
0
8
0
0
1
0
-404
0
0
0
l_
-105
Mz
with 1
A=
0
9
8
1
21
8
I
2208
Vol. 330, No I. pp. 79-93.1993 Prmted m Great Britam
.
22 218
B=
0
0
42
9
0
31
42
9
10
31
42
0
10
31
0
0
10
85
T. N. Lucas These give the matrix products
in Eq. (7) as
-64 F = M2M,A
10 -84
- 1020
=
-316 1
30
-4040
-396 : -24240 -3156 -6218 - 794 and
-3636 -91482 60
42 31 10 9
=
[c H] = M2MIB
which in turn gives the linear system, Eq. (8), as
with solution: d, = 8.984, d, = 2.510, e, = 3.598, e,, = 3.134. The reduced 2nd-order model is thus 8.984s+2.510 R(s) = ~+-5.59~34 , . which can be verified to match G(s) at the four expansion excellent impulse response approximation).
V. Optimal Reduced-order
points
(and gives an
Models
It will now be shown how the method given in the previous sections can be applied in an iterative way to derive optimal reduced-order models in the sense of minimizing the integral-square-error index : I=
;r [y(t) -.vr(t>l’ s0
dt,
(9)
where pa and yl.(t) are the time responses of the full and reduced models, respectively. Longman (19) has shown that the necessary conditions to minimize Z, for an impulse input, are given by
G(J) = R(s) and at
s = pi,
i-1,2,3
,...,
dG(s) ds
dR(s) =
d,r
k Journal d-the
86
Frankhn Institute Pcrgamon Press Ltd
Optimal Model Reduction where -pi are the poles of the kth order reduced transfer function R(s). It is noticed that the Eqs (10) are essentially the multipoint Pad& equations which ensure that the full and reduced transfer functions are matched about the k points s = pI, along with their first derivatives. To try to solve these nonlinear equations for the parameters e, and d, (r = 0, 1,2,. . , k- 1) in R(s), as given in Eq. (2), is a complex numerical task in general if standard approaches are used (16, 17). However, application of the technique presented in Section III is relatively straightforward, and its adaptation to finding optimal reduced models will now be explained. Suppose that it is required to find the optimal kth order transfer function R(s) such that Eq. (9) is minimized for a unit impulse input. In order for Eq. (10) to be satisfied it is seen from (5) that m(s) = (s--p,)*(s--p2)*.
. .(s--p,)*,
(11)
where any factor of the type (s- /3)” appearing in m(s) means that the function and its first (q- 1) derivatives are matched at s = b (18). An initial estimate must be made for the parameters pZ (i = 1, 2, . . . , k), which is a common problem to all optimal methods, except that most require estimates of all the 2k numerator and denominator coefficients. This might be done in a number of ways, for example, retention of k system poles, or use of a stability preserving method to derive a kth degree denominator, or even using an “educated guess” for the reduced model poles. If such an estimate for the reduced denominator is denoted by D,(s), then Eq. (11) becomes m(s) = Di(-s)
(12)
which is computed very easily once Do(s) is known. Notice that Do(s) is quite likely to contain complex conjugate roots, but this causes no problems with the new method, unlike other methods, because only real coefficients will appear in Do(s) and hence m(s). The new method is now applied to reduce G(s) using m(s) in (12) to give the first iterate R(‘)(s) to the optimal model. The denominator Or(s) of R(‘)(s) is then taken as the next estimate for the optimal denominator, i.e. m(s) is updated to m(s) E Df( -s) and the next iterate of the reduced model R’*‘(s) is calculated by the method. The whole process is repeated until the reduced model parameters converge sufficiently to the required optimal model, using the updated estimates of the reduced denominators D,(s) (i = 0, 1, 2,. . .) for the expansion points polynomial m(s) = D,?( -s) to derive the (i+ 1)th iterate R(‘+ ‘j(s). The rate of convergence of the method will obviously depend on how good the initial estimate D,,(s) is to the optimal denominator polynomial. However, it is noticed that the full- and reduced-order models are matched in value along with their first derivatives at k expansion points for each iteration (in essence, similar to the Newton-Raphson method), hence there is a quadratic convergence for the algorithm in the vicinity of the optimal model. In fact, experience with a range of models has shown the method to be quite robust, even if unstable denominators are generated at some point in the iteration they often stabilize and converge quickly to the optimal model. Of course, convergence cannot be guaranteed for all Vol. 330, No. I. pp. 79-93. 1993 Prmted m Great Britain
87
T. N. Lucus possible types of models, and the question of for which systems the method will not converge remains open to future mathematical investigation, To find the optimal reduced-order model for a unit step input a similar procedure is followed to that for the impulse input, except that only the transient part of the response is considered. For the full nth order system given by Eq. (l), the transfer function of the transient response for a step input is seen to be T(s) =
G(s)
K
s
s
or
T(s)
~ KU,S”
= ~~~~~
‘+(h,,- ,-Ka,,-,)s’7~2+ ,Sn ‘+
u,s”+u,,_
“’
.”
+(h,-Ku,) (13)
+a,s+a,
where K = G(0) = bo/uo. As described above, T(s) can now be reduced by the new method to an optimal transient model, say V(S), and the overall reduced step transfer function will be given by
R(4 s
- V(s)+
J-
(14)
This ensures steady-state matching of the full and reduced models as well as optimality in the sense of least integral-square-error as defined by Eq. (9). It should be noticed that, in general, R(s) given in Eq. (14) will have the numerator and denominator polynomial degrees equal in value and correspond to the approximate step response not starting at zero, which might not be suitable for some applications. These ideas on using the transfer function of the transient response may be extended to finding optimal models for ramp and other polynomial inputs with relative ease. For example, for a unit ramp input, T(s) is given by
T(s) = G$) - $
(first two time moments
- $,
where K, = G’(0)
of system),
=
s
VI. Illustrative
$.
Examples for Optimal Models
Example 1 Consider the 4th-order transfer function
88
K2 = G(0)
giving
RCs) V(s)+ :- +
2
and
system given by Bryson
and Carrier
(17) which has the
Journal of the Franklin Institute Pergamon Press Ltd
Optimal Model Reduction 100s+400 G(s) =~fi9~~+113~*+245s+150’
Suppose it is required input, denoted by
to find the optimal
2nd-order
reduced model for an impulse
The simple Pade approximation about s = 0 only is used to obtain the initial estimate D,(s) for the reduced denominator ; however, any of the standard methods could be used, as suggested previously. The new method is now used to iterate to the optimal model, and values for the parameters and integral-square-errors are obtained as shown in Table I. It is seen that the optimal model to 2 decimal places is obtained after 3 iterations, and 6 iterations are needed for 4 decimal place accuracy. The optimal 2nd-order model is then R(s) =
-0.3223s+7.3018 s* + 3.6104s+2.7601
which has a relative integral-square-error value of 0.154% (this is defined by I/ This optimal model is seen to 1;; y’( t ) d t ), m d’mating an excellent approximation. be identical to that obtained by Bryson and Carrier (17), but involves considerably less computation than their method. TABLE I
i
d,
d,
el
e0
-0.5017 -0.3046 -0.3249 -0.3219 -0.3223 -0.3223 -0.3223 -0.3223
8.2472 7.1799 7.3193 7.2994 7.3022 7.3018 7.3018 7.3018
4.090 1 3.5416 3.6200 3.6090 3.6106 3.6104 3.6104 3.6104
3.0927 2.7187 2.7661 2.7592 2.7602 2.7600 2.7601 2.7601
10)
Ix
7.37 4.22 4.16 4.16 4.16 4.16 4.16 4.16
For an optimal step response reduced model, consider the Laplace transform the transient part of the full system response, as given by Eq. (13), i.e.
of
G(s) 2.6667 T(s) = ____ - _______ s s - 2.6667~~ - 50.6667s’ - 301.3333s-
553.3333
s4+19s3+113s2+245s+150 Vol 330, No. I. pp. 79-93. Prmed m Great Britain
1993
89
T. N. Lucas
Applying
to T(s) for reduction
the new method
V(s) =
to a 2nd-order
model denoted
,f,s+fo s2+g,s+go
and again using the simple Pad& approximation about s = 0 to start gives the values of the parameters and integral-square-errors as shown It is seen that convergence to the optimal model (to 4 decimal places) after 3 iterations. The final reduced step response model, as given by R(s)
2.6667
S
S
by
(2.6437s+
the process in Table II. is achieved Eq. (14), is
11.6379)
~~f4.1489sf3.1558
which has a relative integral-square-error excellent approximation.
value of 0.00033%,
again indicating
an
TABLEII
Example
i
.fi
.fb
.41
,40
0 1 2 3 4
-2.6391 -2.6438 -2.6437 -2.6437 -2.6437
- 11.8440 -11.6314 -11.6381 - 11.6379 - 11.6379
4.2111 4.1468 4.1490 4.1489 4.1489
3.2107 3.1540 3.1559 3.1558 3.1558
Ix
1oj
2.27 1.97 1.97 1.97 1.97
2
In the previous example it is clear that the initial estimates of the optimal impulse and step models, obtained by simple Padt approximation, are very good models and convergence to the optimal models was quick. Now it is demonstrated that even with a bad initial estimate it is possible to obtain convergence to the optimal model quite quickly. Consider the 4th-order transfer function given previously in Section IV, i.e. G(s)
=
9s3+42s2+31s+10 s”Ti~3+~is~.+_~~F-t
s
Suppose it is required to find the optimal 2nd-order reduced model for an impulse input. Pade approximation about s = 0 only for the initial estimate gives an unstable model, however, proceeding with the new method even with this bad initial approximation gives the values in Table III for the reduced models along with their integral-square-errors. It is seen that good reduced-order models are quickly generated after the poor initial estimate and the optimal model (to 3 decimal places) is obtained after only 6 iterations of the new method, giving an excellent relative integral-square-error value of 0.06%. It is also interesting that the poles of the optimal model are at - 1.846 f l.l47j, hence the optimal expansion points are complex conjugates at 1.846i 1.147j (repeated) in the Pade approximant. 90
Journal
of the Franklin Institute Pergamon Press Ltd
Optimal Model Reduction TABLE III
i
d,
do
el
e0
0
-2.922 7.287 8.884 8.909 8.905 8.906 8.906 8.906
-0.464 4.106 4.598 4.341 4.389 4.380 4.382 4.382
- 2.208 3.666 3.689 3.692 3.692 3.692 3.692 3.692
-0.371 3.354 4.106 3.958 3.985 3.980 3.981 3.981
1 2 3 4 5 6 7
Ix
103
unstable 423.2 7.34 7.02 7.01 7.01 7.01 7.01
TABLE IV
Ix 100 0
2 3 4 5 6 7 8 9 10
0.634 0.161 0.177 - 0.472 -0.897 - 1.050 - 1.087 - 1.094 - 1.096 - 1.096 - 1.096
0.144 0.443 1.157 2.152 3.219 3.536 3.633 3.649 3.654 3.654 3.654
0.876 0.205 1.089 1.535 2.574 2.874 2.993 3.007 3.015 3.015 3.015
0.328 0.844 1.521 2.940 4.044 4.557 4.663 4.697 4.701 4.703 4.703
46.51 57.56 22.45 8.25 3.33 2.89 2.87 2.86 2.86 2.86 2.86
Similarly, for a step input using Pad& approximation about s = 0 to start the iteration, the values in Table IV are obtained for 2nd-order models of the transient part. Again, after a poor initial estimate, good reduced models are quickly generated and the optimal model (to 3 decimal places) is obtained after 9 iterations, giving the final model as 1.25 (- 1.096s+3.654) ~R(s)= ~ + ~~ ~_._ _.~. ._~ s2 +3.015s+4.703 s s which yields a very good relutiue integral-square-error value of 4.1%. It is interesting to notice again that the optimal expansion points for this step model are complex conjugate pairs at 1.508 f 1.559j (repeated). VII. Conchions
A novel method for calculating optimal reduced-order models in transfer function form has been given. It is based upon a new multipoint Pad& approximation method, applied iteratively, which allows complex expansion points in the iteration Vol 330, No. I, pp. 79-93. Prmted m Great Britam
1993
91
T. N. Lucas without the need for any complex arithmetic. This allows quick computation of the reduced models when compared to current optimal model techniques. The method as presented produces optimal impulse and step response models, but can be extended to ramp and higher degree polynomial inputs by considering the transfer function of the transient part of the response only. It should be noted that the method can be easily adapted to find suboptimal models in the sense of finding a reduced numerator, given a reduced denominator polynomial, which yields the least integral-square-error in Eq. (9) for the given set of reduced system poles. The method is not iterative in this case as only the first half of the conditions in Eq. (10) will apply, i.e. G(s) = R(s) at the negative values of the poles of R(s). The problem then becomes one of solving Eq. (7) for d only. The question of guaranteed convergence of the method to the optimal model is left open to further mathematical investigation. However, the experience of applying the technique to a number of different types of systems suggests that it is very robust. This is probably due to the quadratic convergence property in the vicinity of the optimal model.- An obvious extension of the method is to discrete-time systems and work is proceeding in this direction.
References
(1) R. Genesio and M. Milanese, “A note on the derivation and use of reduced order models”, IEEE Trans. Autom. Control, pp. 118-I 22, 1976.
(2) A. Bultheel and M. Van Bare& “Pade techniques (3) (4) (5) (6) (7) (8) (9) (10)
(11) (12) (13)
92
for model reduction in linear system theory: a survey”, J. Comput. Appl. Math., Vol. 14, pp. 401438, 1986. M. S. Mahmoud and M. G. Singh, “Large Scale Systems Modelling”, Pergamon Press, Oxford, 198 1. T. N. Lucas and A. M. Davidson, “Frequency-domain reduction of linear systems using Schwarz approximation”, ht. J. Control, Vol. 37, pp. 1167-l 178, 1983. C. P. Therapos, “Balanced minimal realisation of SISO systems”, Electron. Lett., Vol. 19, pp. 424426, 1983. C. P. Therapos and J. E. Diamessis, “Sampling method for linear system reduction”, J. Franklin Inst., Vol. 317, pp. 3599371, 1984. T. N. Lucas, “Continued-fraction expansion about two or more points: a flexible approach to linear system reduction”, J. Franklin Inst., Vol. 321, pp. 49960, 1986. C. Hwang and M. Y. Chen, “A multipoint continued-fraction expansion for linear system reduction”, IEEE Trans. Autom. Control, Vol. AC-3 1, pp. 648-65 1, 1986. Y. Shamash, “Linear system reduction using Pade approximation to allow retention of dominant poles”, ht. J. Control, Vol. 21, pp. 2577272, 1975. M. F. Hutton and B. Friedland, “Routh approximations for reducing order of linear, time-invariant sytems”, IEEE Trans. Autom. Control, Vol. AC-20, pp. 3299337, 1975. of transfer functions by the T. C. Chen, C. Y. Chang and K. W. Han, “Reduction Stability Equation method”, J. Franklin Inst., Vol. 308, pp. 389404, 1979. W. Krajewski, A. Lepschy, G. A. Mian and U. Viaro, “On model reduction by Lzoptimal pole retention”, J. Franklin Inst., Vol. 327, pp. 61-70, 1990. “New reduction technique by step S. S. Lamba, R. Gorez and B. Bandyopadhyay, error minimization for multivariable systems”, ht. J. Syst. Sci., Vol. 19, pp. 9991009, 1988.
Optimal Model Reduction (14) C. P. Therapos, “Low-order modelling via constrained least squares minimisation”, Electron. Lett., Vol. 24, pp. 5499550, 1988. Int. J. Control, Vol. 32, pp. 141~747, (15) R. Luus, “Optimization in model reduction”, 1980. (16) R. Vepa, “On reducing the order of transfer functions”, J. Dyn. Syst. Meas. Control, Vol. 108, pp. 27&272, 1986. (17) A. E. Bryson and A. Carrier, “Second order algorithm for optimal model order reduction”, J. Guid. Control Dyn., Vol. 13, pp. 8877892, 1990. (18) T. N. Lucas, “Finding roots by deflated polynomial approximation”, J. Franklin Inst., Vol. 327, pp. 819-830, 1990. (19) I. M. Longman, “Best rational function approximation for Laplace transform inversion”, SIAM J. Math. Anal., Vol. 5, pp. 574-580, 1974.
Vol 330. No. I. pp. 79-03, Prmtcd in Chat Bnlam
1993
93
by T. N. LUCAS Department of Mathematical and Computer Sciences, Dundee Institute of Technology, Bell Street, Dundee DDI IHG,
U.K.
A novel method of obtaining optimal reduced-order models of linear system transfer functions is presented. It uses the popular multipoint Pad& approximation technique in am iterative way to generate eflciently the optimal models. Central to the method is a new way of calculating Pad& approximants about many points by reducing them to equivalent Taylor series approximants. Optimal reduced-order models for impulse and step inputs are considered, and it is seen how the method may be extended to ramp and other polynomial inputs. Numerical examples are given to demonstrate the method.
ABSTRACT:
I. Introduction Reducing the order of linear time-invariant systems has continued to be a topic of interest in control system design over the past couple of decades or so. Its importance stems from the ever increasing complexity of modern control system models, which tend to be of high-order when linearized about convenient operating points. The use of reduced-order plant models to design an adequate controller is now an accepted practice among control engineers and, consequently, an array of different approaches to the “order reduction” problem has been well documented in the literature, e.g. (l-3). single-output systems the transfer function description For single-input, (frequency-domain) tends to find favour over the state-space description (timedomain) which is, however, more suited to multi-input, multi-output systems than the transfer function. Both approaches to the model reduction problem are equally important and it is found that the two approaches are not entirely divorced from each other, where many examples of links between time- and frequency-domain methods have been cited, e.g. (4, 5). The principal interest of this paper is in the frequency-domain approach, and, in particular, that of multipoint Pad&/continued-fraction approximation methods (68). These are known to be very flexible and computationally efficient methods of reduction and as such have proved to be popular techniques to use. However,
The Frmkhn
hmtutc
0016
0032;93 $j.OO+O
00
79
T. N. their biggest drawback is that they cannot guarantee stability of the reduced-order model if the original system is stable, unless a subset of the system poles are used as expansion points to ensure that these poles are retained in the reduced model (799). To overcome this stability problem, a number of methods have been developed alongside the Pad& methods, ensuring stability of the reduced models. These are based mostly on the Routh table properties utilised by Hutton and Friedland (10) and Chen et al. (11). Such methods have proved to be popular because of their relative ease of computation along with the stability preserving properties. However, it is generally recognised that the reduced-order models produced by these methods tend not to be such good approximations as those of the multipoint Pad& methods when they produce stable models. This is due to twice the number of system parameters being retained in the Pade models compared to those of the “stability preserving” methods, in which the reduced denominator polynomials are derived before calculating the numerators via matching appropriate system parameters (time moments, Markov parameters or frequency points). More recently, interest has developed in finding “optimal” and “sub-optimal” reduced models in the frequency-domain by minimizing the integral-square-error :
I=
sOx
[y(t) -vr(O12dt>
where y(t) and y,(t) are the full and reduced-order models’ impulse (or step) responses, respectively. The sub-optimal transfer function models are usually derived first by finding a reduced denominator polynomial by one of the stability methods, and then choosing the numerator to minimize 1, e.g. (12-14). Such models require a modest amount of computation, having eventually to solve a system of linear equations for the numerator coefficients. Various commendable approaches for finding the highly desirable optimal reduced-order transfer function models have been given by, among others, Luus (15), Vepa (16), and Bryston and Carrier (17). However, the relatively little use made of these methods by control engineers appears to be due to the quite substantial amount of computation involved compared to the other methods already mentioned. In general, these optimal methods require solving sets of nonlinear equations by classical numerical techniques such as the generalized NewtonRaphson method. This paper presents a novel method for obtaining optimal reduced-order models by the application of multipoint Pad& approximation. Tt removes the need for classical numerical methods, and instead uses the efficient PadC technique in an iterative way to converge onto the optimal reduced-order model. Both impulse and step response models may be obtained by the method in a straightforward manner, requiring only the solution of linear sets of equations. The technique is seen to be possible only because of a new way of calculating multipoint PadC approximants by means of multipoint Taylor approximants of polynomials. This result is discussed in Sections IT, III and IV before its use in optimal model reduction is developed in Section V. Numerical examples are given in Section VI to demonstrate the application of the method. 80
Journal of the Franklm lnst~tute Pergamon Press Ltd
Optimal Model Reduction II. Basis of the Method Consider
the nth order transfer G@)
which is to be reduced
=
function,
in usual notation,
b,-,s-‘+.~~_~
~~..
a,s”+
“.
given by
+h,s+h,
(1)
... +a,s+a,
to the kth order model dA_,skm’+ ... +d,s+d, R(s) = s +ek_,s ‘-I+ .:. +e,s+e,
(2)
by Pad& approximation about the 2k points, s = p,, i = 1, 2,. . . ,2k. The problem is to find the 2k parameters d, and e, (v = 0, 1, 2,. . , k- I), such that G(p,) = R(p,), Using Eqs (l), (2) and (3), it is obvious P(S) = Q(S)
at
i = 1,2,. . . ,2k.
(3)
that (3) is equivalent
s = pl,
to
i = 1,2,. . . ,2k,
(4)
where P(S) = (a,?+
... +a,s+ao)(dk_,sk-‘+
... +d,s+d,)
and Q(S) = (b,_,s”-‘+
... +b,,)(s”+ek_,sk-‘+
... +e,s+e”).
This relationship changes the problem to one of finding the coefficients of two (n+k1)th degree polynomials, P(s) and Q(s), such that their values are equal at the 2k expansion points pi. These coefficients will of course be linear combinations ofthed,ande,(r=0,1,2 ,..., k-l). The central idea of the new method is now to find the unique polynomials, of degree (2k - l), which interpolate P(S) and Q(S) at the 2k expansion points s = pi (i = 1, 2, . . . ,2k). By definition, these polynomials must be identical because P(s) and Q(s) have the same values at the expansion points. Hence the 2k coefficients of both polynomials may be equated to give a set of 2k linear equations to solve for the parameters d, and e,. These equations will, in effect, be the multipoint Pade equations of the system. The reduced degree interpolating polynomials can be derived very easily by a Routh-type algorithm (18) and may be thought of as multipoint Taylor approximants of P(S) and Q(S), respectively, about s = pI. Furthermore, the algorithm is well suited to cope with real, complex and multiple expansion points. Lucas (18) shows that dividing a polynomial T(S) of degree n by another, H(s) of degree m (m < n), from highest powers, gives a remainder polynomial of degree (m- 1) which is the Taylor approximation of T(S) about the m roots of H(s). To illustrate this important concept, consider the simple case of reducing the 3rd-order system : 8?+6s+2 G(s) = ~ ~ s3+4s2+5s+2 81
T. N. Lucas to one of 2nd order, given by d,s+d, R(s) = yp s +e,s+erj by multipoint approximation about the four points s = 0, 1, 2, 3. The divisor polynomial in the algorithm is thus s(s-l)(s-2)(s-3) and the polynomials
- s4-66s3+lls2-66s
P(s) and Q(s), as defined by Eq. (4), are given by
P(s) = d,s4 + (4d, + do)s3 + (5d, +4d,)s2 + (2d, + 5d,,)s+ 2d, and Q(s) = 8s4+(8e,
+6)s3+(6e,
The table for the 3rd-degree given as follows (18) :
Taylor
+8e,+2)s2+(2e, approximant
+6e0)s+2eo.
of P(s) about
s = 0, 1, 2, 3 is
d,
(4d, + do)
(5d, +4&J
(2d, + 5do)
2d,,
1
-6
11
-6
0
(lOd, +d,)
(-6d,
+4do)
2do
(8d, +5do)
where the first row coefficients are those of P(s), the second row are those of the divisor (expansion points) polynomial, and the third row are those of the 3rd degree approximant of P(s) calculated from rows 1 and 2 by the usual Routh division algorithm. Similarly, the approximant of Q(s) is given by 8
(8er +6)
1
-6
(8e,+54) Equating
(6e, +8e,+2) 11
(6e,+8eo-86)
like coefficients
(2e,+6eo+48)
in both approximants
(2e, +6eo)
2eo
-6
0 ’
2e,
gives the four linear equations
lOd, +do = 8e, +54 -6d,
+4do = 6e, +8eo-86
8d, +5d,
= 2e, +6eo+48
2d,, = 2eo which solve to give the reduced do = 5.3488,
model parameters
d, = 7.4884,
as :
e. = 5.3488,
e, = 3.2791.
It is easily verified that this model fits G(s) at all four expansion points. Although the above 2nd-order reduced model was obtained in a straightforward manner, it is easy to see that for larger systems requiring more than a single division step for approximating P(s) and Q(s) (or the reduced model is of order greater than 2), then the Taylor approximant tables can become rather cumbersome and therefore unattractive to use. Consequently, it is desirable to produce a more 82
Journal
of the Franklin Institute Pergamon Press Ltd
Optimul Model Reduction compact form of the algorithm the next section.
for computational
purposes.
This is discussed
in
III. Matrix Form In general, the (n+k) coefficients of P(s) and Q(Y), as defined by Eq. (4), may be written in matrix form as the vectors Ad and B[l e] ‘, respectively, where
a, a,-, an-2
0
a, a,,- 1
. an-k+1
’
.
.
.
. . .
A= a0
.
al
.
a0 .
.
.
.
0
4- I ‘6-z
a, . ,
d=
ak-I
d,
ak-2
do
_
. .
a0
0
6-k-l
’
. .
.
.
.
.
.
bn-,
B=
and ho
b, bo
b,
.’ .
.
bk h-1
.
.
.
e =
0 bo
Suppose that the reduced model is obtained which are the roots of the polynomial m(s) = s2k+~2k_iS2k~‘+
by expanding “’ +~is+&,.
about
the 2k points
(5)
Now, from the elementary theory of transformation matrices, it is seen that the third rows in the multipoint Taylor approximant tables of P(s) and Q(s) (after the first division by m(s)) will be given by the entries in the two vectors M,Ad and Vol 330, No. I, pp. 79-93. 1993 Prmted m Great Bntain
83
T. N. Lucas M, B[ 1 e] ‘, respectively,
where
-XZk-
1 0 0 . .
I
-aZk&Z -a2&3
M,
=
-sIo 0 0
1 0 . .
1 .
0 ’ O....Ol () . . .
0
0 .
.
.
0
(n-f- T-1)
1
. .
. .
0 .
o......
<
1 .
.
.
0
I
,
(n+k)
Similarly, the next rows of each approximant table (after the division by m(s) into the third row polynomial) may be represented by the vectors obtained by premultiplying the above matrix products by the transformation matrix M2, which is formed from M, by deleting its last row and column. Hence, in general, as there will be (n-k) division steps, the coefficients of the 2kth degree Taylor approximants of P(s) and Q(s), respectively, are given by the elements of the vectors
where
-x2h-2
1 0
-@2h-3
0
-a2k
I
1 0
Mi = -cl0
0 0
L -
0
0 1
.
.
.
.
.
.
(n+k-i) .
0
O....Ol .
1 . 0
o..... (n+k+l
(6)
-i)
1A >
fori= 1,2 ,..., n-k. Equating these vectors to solve for d and e, gives Fd = [c
H]
(7)
where F = Mplmk. ..M?M,A and [c H] = M,, k.. M2M,B, with c being the first column of this matrix. Equation (7) may be rearranged to give the linear multipoint Padt equations as x4
Optimal Model Reduction
II
[F-H]
t
=c
which can be directly solved by any standard algorithm for linear equations. Notice that complex expansion points, which always occur in conjugate pairs for physical realisation of the reduced model, cause no problems for the entries in the M, matrices as real polynomial coefficients LX,will always result. Notice also that the ordinary Pad6 approximant (about s = 0 only) is obtained by putting CI, = 0, m = 0, 1, 2,. . . ,2k- 1, in the A4, matrices. The whole process therefore consists of simple matrix multiplications and solving 2k linear equations.
IV. Illustrative Consider
Example for Multipoint
the 4th-order
transfer
G(s) =
Method
function
(8) given by
9s’+42s2+31s+10 s4+8S3+21s2+22s+8’
For expansion points such as s = f 2j, 1 + lOj, most existing multipoint continuedfraction/Pad6 methods would not be able to cope without resorting to complex arithmetic, which renders them unattractive to use (the choice of complex expansion points will be seen to be important in obtaining optimal models later). However, the given method quickly produces a second-order model as follows. From the given expansion points, the polynomial m(s) is given by m(s) = (.~~+4)[(s-1)~+100] and the corresponding
rs4-22S3+105s2-8s+404
transformation
matrices
are
2 M,
=
and
1
0
0
0
0
1
0
0
8
0
0
1
0
-404
0
0
0
l_
-105
Mz
with 1
A=
0
9
8
1
21
8
I
2208
Vol. 330, No I. pp. 79-93.1993 Prmted m Great Britam
.
22 218
B=
0
0
42
9
0
31
42
9
10
31
42
0
10
31
0
0
10
85
T. N. Lucas These give the matrix products
in Eq. (7) as
-64 F = M2M,A
10 -84
- 1020
=
-316 1
30
-4040
-396 : -24240 -3156 -6218 - 794 and
-3636 -91482 60
42 31 10 9
=
[c H] = M2MIB
which in turn gives the linear system, Eq. (8), as
with solution: d, = 8.984, d, = 2.510, e, = 3.598, e,, = 3.134. The reduced 2nd-order model is thus 8.984s+2.510 R(s) = ~+-5.59~34 , . which can be verified to match G(s) at the four expansion excellent impulse response approximation).
V. Optimal Reduced-order
points
(and gives an
Models
It will now be shown how the method given in the previous sections can be applied in an iterative way to derive optimal reduced-order models in the sense of minimizing the integral-square-error index : I=
;r [y(t) -.vr(t>l’ s0
dt,
(9)
where pa and yl.(t) are the time responses of the full and reduced models, respectively. Longman (19) has shown that the necessary conditions to minimize Z, for an impulse input, are given by
G(J) = R(s) and at
s = pi,
i-1,2,3
,...,
dG(s) ds
dR(s) =
d,r
k Journal d-the
86
Frankhn Institute Pcrgamon Press Ltd
Optimal Model Reduction where -pi are the poles of the kth order reduced transfer function R(s). It is noticed that the Eqs (10) are essentially the multipoint Pad& equations which ensure that the full and reduced transfer functions are matched about the k points s = pI, along with their first derivatives. To try to solve these nonlinear equations for the parameters e, and d, (r = 0, 1,2,. . , k- 1) in R(s), as given in Eq. (2), is a complex numerical task in general if standard approaches are used (16, 17). However, application of the technique presented in Section III is relatively straightforward, and its adaptation to finding optimal reduced models will now be explained. Suppose that it is required to find the optimal kth order transfer function R(s) such that Eq. (9) is minimized for a unit impulse input. In order for Eq. (10) to be satisfied it is seen from (5) that m(s) = (s--p,)*(s--p2)*.
. .(s--p,)*,
(11)
where any factor of the type (s- /3)” appearing in m(s) means that the function and its first (q- 1) derivatives are matched at s = b (18). An initial estimate must be made for the parameters pZ (i = 1, 2, . . . , k), which is a common problem to all optimal methods, except that most require estimates of all the 2k numerator and denominator coefficients. This might be done in a number of ways, for example, retention of k system poles, or use of a stability preserving method to derive a kth degree denominator, or even using an “educated guess” for the reduced model poles. If such an estimate for the reduced denominator is denoted by D,(s), then Eq. (11) becomes m(s) = Di(-s)
(12)
which is computed very easily once Do(s) is known. Notice that Do(s) is quite likely to contain complex conjugate roots, but this causes no problems with the new method, unlike other methods, because only real coefficients will appear in Do(s) and hence m(s). The new method is now applied to reduce G(s) using m(s) in (12) to give the first iterate R(‘)(s) to the optimal model. The denominator Or(s) of R(‘)(s) is then taken as the next estimate for the optimal denominator, i.e. m(s) is updated to m(s) E Df( -s) and the next iterate of the reduced model R’*‘(s) is calculated by the method. The whole process is repeated until the reduced model parameters converge sufficiently to the required optimal model, using the updated estimates of the reduced denominators D,(s) (i = 0, 1, 2,. . .) for the expansion points polynomial m(s) = D,?( -s) to derive the (i+ 1)th iterate R(‘+ ‘j(s). The rate of convergence of the method will obviously depend on how good the initial estimate D,,(s) is to the optimal denominator polynomial. However, it is noticed that the full- and reduced-order models are matched in value along with their first derivatives at k expansion points for each iteration (in essence, similar to the Newton-Raphson method), hence there is a quadratic convergence for the algorithm in the vicinity of the optimal model. In fact, experience with a range of models has shown the method to be quite robust, even if unstable denominators are generated at some point in the iteration they often stabilize and converge quickly to the optimal model. Of course, convergence cannot be guaranteed for all Vol. 330, No. I. pp. 79-93. 1993 Prmted m Great Britain
87
T. N. Lucus possible types of models, and the question of for which systems the method will not converge remains open to future mathematical investigation, To find the optimal reduced-order model for a unit step input a similar procedure is followed to that for the impulse input, except that only the transient part of the response is considered. For the full nth order system given by Eq. (l), the transfer function of the transient response for a step input is seen to be T(s) =
G(s)
K
s
s
or
T(s)
~ KU,S”
= ~~~~~
‘+(h,,- ,-Ka,,-,)s’7~2+ ,Sn ‘+
u,s”+u,,_
“’
.”
+(h,-Ku,) (13)
+a,s+a,
where K = G(0) = bo/uo. As described above, T(s) can now be reduced by the new method to an optimal transient model, say V(S), and the overall reduced step transfer function will be given by
R(4 s
- V(s)+
J-
(14)
This ensures steady-state matching of the full and reduced models as well as optimality in the sense of least integral-square-error as defined by Eq. (9). It should be noticed that, in general, R(s) given in Eq. (14) will have the numerator and denominator polynomial degrees equal in value and correspond to the approximate step response not starting at zero, which might not be suitable for some applications. These ideas on using the transfer function of the transient response may be extended to finding optimal models for ramp and other polynomial inputs with relative ease. For example, for a unit ramp input, T(s) is given by
T(s) = G$) - $
(first two time moments
- $,
where K, = G’(0)
of system),
=
s
VI. Illustrative
$.
Examples for Optimal Models
Example 1 Consider the 4th-order transfer function
88
K2 = G(0)
giving
RCs) V(s)+ :- +
2
and
system given by Bryson
and Carrier
(17) which has the
Journal of the Franklin Institute Pergamon Press Ltd
Optimal Model Reduction 100s+400 G(s) =~fi9~~+113~*+245s+150’
Suppose it is required input, denoted by
to find the optimal
2nd-order
reduced model for an impulse
The simple Pade approximation about s = 0 only is used to obtain the initial estimate D,(s) for the reduced denominator ; however, any of the standard methods could be used, as suggested previously. The new method is now used to iterate to the optimal model, and values for the parameters and integral-square-errors are obtained as shown in Table I. It is seen that the optimal model to 2 decimal places is obtained after 3 iterations, and 6 iterations are needed for 4 decimal place accuracy. The optimal 2nd-order model is then R(s) =
-0.3223s+7.3018 s* + 3.6104s+2.7601
which has a relative integral-square-error value of 0.154% (this is defined by I/ This optimal model is seen to 1;; y’( t ) d t ), m d’mating an excellent approximation. be identical to that obtained by Bryson and Carrier (17), but involves considerably less computation than their method. TABLE I
i
d,
d,
el
e0
-0.5017 -0.3046 -0.3249 -0.3219 -0.3223 -0.3223 -0.3223 -0.3223
8.2472 7.1799 7.3193 7.2994 7.3022 7.3018 7.3018 7.3018
4.090 1 3.5416 3.6200 3.6090 3.6106 3.6104 3.6104 3.6104
3.0927 2.7187 2.7661 2.7592 2.7602 2.7600 2.7601 2.7601
10)
Ix
7.37 4.22 4.16 4.16 4.16 4.16 4.16 4.16
For an optimal step response reduced model, consider the Laplace transform the transient part of the full system response, as given by Eq. (13), i.e.
of
G(s) 2.6667 T(s) = ____ - _______ s s - 2.6667~~ - 50.6667s’ - 301.3333s-
553.3333
s4+19s3+113s2+245s+150 Vol 330, No. I. pp. 79-93. Prmed m Great Britain
1993
89
T. N. Lucas
Applying
to T(s) for reduction
the new method
V(s) =
to a 2nd-order
model denoted
,f,s+fo s2+g,s+go
and again using the simple Pad& approximation about s = 0 to start gives the values of the parameters and integral-square-errors as shown It is seen that convergence to the optimal model (to 4 decimal places) after 3 iterations. The final reduced step response model, as given by R(s)
2.6667
S
S
by
(2.6437s+
the process in Table II. is achieved Eq. (14), is
11.6379)
~~f4.1489sf3.1558
which has a relative integral-square-error excellent approximation.
value of 0.00033%,
again indicating
an
TABLEII
Example
i
.fi
.fb
.41
,40
0 1 2 3 4
-2.6391 -2.6438 -2.6437 -2.6437 -2.6437
- 11.8440 -11.6314 -11.6381 - 11.6379 - 11.6379
4.2111 4.1468 4.1490 4.1489 4.1489
3.2107 3.1540 3.1559 3.1558 3.1558
Ix
1oj
2.27 1.97 1.97 1.97 1.97
2
In the previous example it is clear that the initial estimates of the optimal impulse and step models, obtained by simple Padt approximation, are very good models and convergence to the optimal models was quick. Now it is demonstrated that even with a bad initial estimate it is possible to obtain convergence to the optimal model quite quickly. Consider the 4th-order transfer function given previously in Section IV, i.e. G(s)
=
9s3+42s2+31s+10 s”Ti~3+~is~.+_~~F-t
s
Suppose it is required to find the optimal 2nd-order reduced model for an impulse input. Pade approximation about s = 0 only for the initial estimate gives an unstable model, however, proceeding with the new method even with this bad initial approximation gives the values in Table III for the reduced models along with their integral-square-errors. It is seen that good reduced-order models are quickly generated after the poor initial estimate and the optimal model (to 3 decimal places) is obtained after only 6 iterations of the new method, giving an excellent relative integral-square-error value of 0.06%. It is also interesting that the poles of the optimal model are at - 1.846 f l.l47j, hence the optimal expansion points are complex conjugates at 1.846i 1.147j (repeated) in the Pade approximant. 90
Journal
of the Franklin Institute Pergamon Press Ltd
Optimal Model Reduction TABLE III
i
d,
do
el
e0
0
-2.922 7.287 8.884 8.909 8.905 8.906 8.906 8.906
-0.464 4.106 4.598 4.341 4.389 4.380 4.382 4.382
- 2.208 3.666 3.689 3.692 3.692 3.692 3.692 3.692
-0.371 3.354 4.106 3.958 3.985 3.980 3.981 3.981
1 2 3 4 5 6 7
Ix
103
unstable 423.2 7.34 7.02 7.01 7.01 7.01 7.01
TABLE IV
Ix 100 0
2 3 4 5 6 7 8 9 10
0.634 0.161 0.177 - 0.472 -0.897 - 1.050 - 1.087 - 1.094 - 1.096 - 1.096 - 1.096
0.144 0.443 1.157 2.152 3.219 3.536 3.633 3.649 3.654 3.654 3.654
0.876 0.205 1.089 1.535 2.574 2.874 2.993 3.007 3.015 3.015 3.015
0.328 0.844 1.521 2.940 4.044 4.557 4.663 4.697 4.701 4.703 4.703
46.51 57.56 22.45 8.25 3.33 2.89 2.87 2.86 2.86 2.86 2.86
Similarly, for a step input using Pad& approximation about s = 0 to start the iteration, the values in Table IV are obtained for 2nd-order models of the transient part. Again, after a poor initial estimate, good reduced models are quickly generated and the optimal model (to 3 decimal places) is obtained after 9 iterations, giving the final model as 1.25 (- 1.096s+3.654) ~R(s)= ~ + ~~ ~_._ _.~. ._~ s2 +3.015s+4.703 s s which yields a very good relutiue integral-square-error value of 4.1%. It is interesting to notice again that the optimal expansion points for this step model are complex conjugate pairs at 1.508 f 1.559j (repeated). VII. Conchions
A novel method for calculating optimal reduced-order models in transfer function form has been given. It is based upon a new multipoint Pad& approximation method, applied iteratively, which allows complex expansion points in the iteration Vol 330, No. I, pp. 79-93. Prmted m Great Britam
1993
91
T. N. Lucas without the need for any complex arithmetic. This allows quick computation of the reduced models when compared to current optimal model techniques. The method as presented produces optimal impulse and step response models, but can be extended to ramp and higher degree polynomial inputs by considering the transfer function of the transient part of the response only. It should be noted that the method can be easily adapted to find suboptimal models in the sense of finding a reduced numerator, given a reduced denominator polynomial, which yields the least integral-square-error in Eq. (9) for the given set of reduced system poles. The method is not iterative in this case as only the first half of the conditions in Eq. (10) will apply, i.e. G(s) = R(s) at the negative values of the poles of R(s). The problem then becomes one of solving Eq. (7) for d only. The question of guaranteed convergence of the method to the optimal model is left open to further mathematical investigation. However, the experience of applying the technique to a number of different types of systems suggests that it is very robust. This is probably due to the quadratic convergence property in the vicinity of the optimal model.- An obvious extension of the method is to discrete-time systems and work is proceeding in this direction.
References
(1) R. Genesio and M. Milanese, “A note on the derivation and use of reduced order models”, IEEE Trans. Autom. Control, pp. 118-I 22, 1976.
(2) A. Bultheel and M. Van Bare& “Pade techniques (3) (4) (5) (6) (7) (8) (9) (10)
(11) (12) (13)
92
for model reduction in linear system theory: a survey”, J. Comput. Appl. Math., Vol. 14, pp. 401438, 1986. M. S. Mahmoud and M. G. Singh, “Large Scale Systems Modelling”, Pergamon Press, Oxford, 198 1. T. N. Lucas and A. M. Davidson, “Frequency-domain reduction of linear systems using Schwarz approximation”, ht. J. Control, Vol. 37, pp. 1167-l 178, 1983. C. P. Therapos, “Balanced minimal realisation of SISO systems”, Electron. Lett., Vol. 19, pp. 424426, 1983. C. P. Therapos and J. E. Diamessis, “Sampling method for linear system reduction”, J. Franklin Inst., Vol. 317, pp. 3599371, 1984. T. N. Lucas, “Continued-fraction expansion about two or more points: a flexible approach to linear system reduction”, J. Franklin Inst., Vol. 321, pp. 49960, 1986. C. Hwang and M. Y. Chen, “A multipoint continued-fraction expansion for linear system reduction”, IEEE Trans. Autom. Control, Vol. AC-3 1, pp. 648-65 1, 1986. Y. Shamash, “Linear system reduction using Pade approximation to allow retention of dominant poles”, ht. J. Control, Vol. 21, pp. 2577272, 1975. M. F. Hutton and B. Friedland, “Routh approximations for reducing order of linear, time-invariant sytems”, IEEE Trans. Autom. Control, Vol. AC-20, pp. 3299337, 1975. of transfer functions by the T. C. Chen, C. Y. Chang and K. W. Han, “Reduction Stability Equation method”, J. Franklin Inst., Vol. 308, pp. 389404, 1979. W. Krajewski, A. Lepschy, G. A. Mian and U. Viaro, “On model reduction by Lzoptimal pole retention”, J. Franklin Inst., Vol. 327, pp. 61-70, 1990. “New reduction technique by step S. S. Lamba, R. Gorez and B. Bandyopadhyay, error minimization for multivariable systems”, ht. J. Syst. Sci., Vol. 19, pp. 9991009, 1988.
Optimal Model Reduction (14) C. P. Therapos, “Low-order modelling via constrained least squares minimisation”, Electron. Lett., Vol. 24, pp. 5499550, 1988. Int. J. Control, Vol. 32, pp. 141~747, (15) R. Luus, “Optimization in model reduction”, 1980. (16) R. Vepa, “On reducing the order of transfer functions”, J. Dyn. Syst. Meas. Control, Vol. 108, pp. 27&272, 1986. (17) A. E. Bryson and A. Carrier, “Second order algorithm for optimal model order reduction”, J. Guid. Control Dyn., Vol. 13, pp. 8877892, 1990. (18) T. N. Lucas, “Finding roots by deflated polynomial approximation”, J. Franklin Inst., Vol. 327, pp. 819-830, 1990. (19) I. M. Longman, “Best rational function approximation for Laplace transform inversion”, SIAM J. Math. Anal., Vol. 5, pp. 574-580, 1974.
Vol 330. No. I. pp. 79-03, Prmtcd in Chat Bnlam
1993
93