On the asymptotic accuracy of reduced-order models

Popular model reduction methods can easily be adapted to retain the asymptotic response to inputs with rational transform. To this purpose, the forced response of the high-order system is decomposed into a transient and a steady-state component. Then, the reduced-order model is obtained by combining the unaltered steady-state component with an approximation of the transient component. Examples show that forcing the reduced-order model to retain the steady-state component does not compromise the transient accuracy.


INTRODUCTION
The model reduction problem has received considerable attention since the dawn of system and control theory (see [1][2][3] for recent surveys) and the subject has reached a high level of maturity, even if important new impulses to research are still coming from emerging areas [4] and clever ideas [5][6][7][8]. In fact, despite the increase of computing capabilities witnessed in recent years, reduced-order models are being widely used for both analysis and design purposes; suffice it here to recall the control of flexible mechanical structures [9], the simulation of integrated circuits [10], the study of power grid networks [11] and climate modelling [12]. Usually, reduced-order models are obtained by referring to the minimisation of a norm or a semi-norm of the approximation error in the response to suitable inputs. For instance, the classic Padé technique and the so-called moment-matching methods set to zero the error transform at specific frequencies, which, for stable systems, entails zeroing asymptotically the error in the response to harmonic input signals at these frequencies. These techniques, however, do not ensure stability retention, at least in their original versions. On the other hand, popular stability-preserving methods, such as balanced truncation [13], Hankel-norm approximation [14] and L 2 model reduction [1], do not ensure the retention of the asymptotic response to any canonical input (except possibly for the step, at the expense of exact optimality, by resorting to singular perturbation approximation [15] [16]), even if their transient accuracy is usually very good compared to alternative techniques. Moreover, fairly robust (and easily accessible) algorithms exist for their implementation [17][18][19][20], which makes them applicable to the Manuscript  reduction of very large-scale systems despite their nonnegligible numerical complexity (related to the solution of high-dimensional Lyapunov equations in the case balanced truncation or the iterative solution of interpolation problems in the case of L 2 approximation). This note suggests a simple way to ensure the retention of the asymptotic response to predetermined persistent inputs with rational transform. To this purpose, the transform of the related original forced response is decomposed into the sum of a system component and an input component, according to the terminology in [21], where the model reduction problem is not considered. For stable systems, the first component accounts for the system's transient behaviour, whereas the second accounts for its asymptotic behaviour. Then, the approximation procedure is applied only to the system component whereas the input component is reproduced exactly. In a sense, a similar approach has been followed in [22,23] for approximating unstable systems via stable/antistable decomposition. The present contribution is limited to single-input single-output (SISO) systems. Its extension to multi-input multi-output (MIMO) system is not straightforward and would require considering either transfer matrices or state-space descriptions.

DECOMPOSITION OF THE FORCED RESPONSE
The suggested reduction technique is based on the additive decomposition of the forced response presented in [21] and briefly outlined next. Denote the strictly-proper rational transfer function of the original system by W (s) = B(s)/A(s), where A(s) and B(s) are coprime polynomials, and the proper rational Laplace transform of the input u(t) by U(s) = N(s)/D(s), where N(s) and D(s) are also coprime polynomials. Throughout, it is further assumed that B(s) and D(s) are coprime as well as N(s) and A(s). In this way, no cancellation occurs in the expression of the Laplace transform of the forced response y f (t) If polynomials A(s) and D(s) at the denominator of (1) have common roots, they can be expressed as A(s) = A(s)C A (s) and D(s) =D(s)C D (s), where C A (s) is a factor containing all the roots of A(s) (with their multiplicities) that are also roots of D(s) and C D (s) is a factor containing all the roots of D(s) (with their multiplicities) that are also roots of A(s). Denoting the product of these two factors by C(s) = C A (s)C D (s), the three pairs (Ā(s),D(s)), (Ā(s),C(s)) and (C(s),D(s)) are all formed by coprime polynomials, so that Y f (s) can be expressed uniquely as the sum of three strictly-proper rational functions as: where , as may easily be proved by resorting to the Bezout identity (see, e.g., [24, p. 204]). Since the inverse Laplace transform y w (t) of the first addendum on the right-hand side of (2), whose poles are those of W (s) not in common with U(s), is a combination of system modes, it is called system component in [21]. Similarly, y u (t) = LT −1 [Y u (s)] is a combination of the modes of u(t) that are not in common with w(t) and is therefore called input component. The poles of the third addendum Y c (s) in (2) are those common to W (s) and U(s), but their multiplicity is the sum of the multiplicities of the same poles in W (s) and U(s); therefore, y c (t) = LT −1 [Y c (s)] can reasonably be called resonant component. A characterisation of these components in terms of the solution of either a homogeneous or a nonhomogeneous differential equation is provided in [21] together with some suggestions on how to extend the definitions to non-rational input and system transforms.
In most cases, A(s) and D(s) have no common factors so that the forced-response is given by the sum of the system and input components only, i.e., This simplifying assumption will be made in the sequel since the consideration of the more general case would entail a substantial increase in notation without a corresponding gain in insight. Although the decomposition (3) holds for stable and unstable systems, as well as for both persistent and vanishing inputs, it is particularly meaningful for the purpose of model reduction when u(t), and thus y u (t), is persistent and w(t), and thus y w (t), tends asymptotically to zero, i.e., the system is BIBO stable. In this case, y u (t) corresponds to the steady-state (or asymptotic) response and y w (t) to the transient response. In the sequel, reference is made to such a situation.

REDUCTION PROCEDURE
According to the previous considerations, it is assumed next that the original system is BIBO stable and that all of the poles of U(s) are in the closed right half-plane so that A(s) and D(s) are necessarily coprime (C(s) = 1) and all of the modes of u(t) are persistent. The reduction procedure consists of the following steps.
according to a suitable reduction criterion (e.g., balance truncation, Hankel-norm approximation or L 2 reduction).
2) Form the reduced-order forced response transform as where Y a (s) = N a (s)/A a (s) is an auxiliary strictly-proper stable rational function whose task is to guarantee that step 3 admits a solution.
3) Determine the strictly-proper reduced-order transfer function W r (s) as Remark 1: The auxiliary term Y a (s) is needed since, otherwise, the reduction procedure would admit a solution only in the case of impulsive inputs (see Remark 3).

Remark 2:
The poles of Y a (s) must lie in the open left half-plane not to influence the asymptotic response. Also, these poles should lie far away from the poles of Y w r (s) not to affect appreciably the transient behaviour at the dominant frequencies of the reduced model.

Remark 3:
The order of Y a (s), i.e., the degree of its denominator, depends on the order of the component Y u (s) that must be retained and, thus, on the order of the input. To clarify this point, consider equation (5) which, in view of steps 1 and 2, leads to the polynomial identity: whose degree, given the properness of U(s) and the strict properness of Y w r (s), Y u (s) and Y a (s), is It follows that, by equating the coefficients of the equal powers of s on both sides of (6), a system of N + 1 equations is obtained. In order Remark 4: A consequence of Remark 3 is that the order of the reduced model W r (s) is greater than the order of the simplified transient term Y w r (s) by an amount equal to deg[D(s)]. For steps, ramps and sinusoidal inputs, deg[D(s)] ≤ 2 which is usually quite negligible compared to the order of the approximating model of a very high-dimensional system [25].

Remark 5:
The aforementioned approximation procedure ensures that the reduced-order model is strictly proper like the original system, whereas the balanced truncation procedure and its DC-gain-preserving variant, illustrated, e.g., in [16], lead, in general, to exactly-proper reduced-order models.
Remark 6: Clearly, the steady-state response is exactly preserved only for the selected input signal which is assumed to be "stationary," i.e., its time-and frequencydomain characteristics keep unchanged. However, the frequency response of the reduced models is close the original frequency response in the neighbourhood of the input frequency (which is zero for the singularity input functions such as the impulse, step, and ramp functions) as the Bode diagrams for the examples in the next section shows.

EXAMPLES
In the following examples it is assumed that the asymptotic response to be retained is either the one in the step response (first example) or the one in the ramp response (second and third examples). In all cases, the reducedorder transient response to such inputs is determined by minimising the Hankel norm of the approximation error. Of course, other simplification criteria could be chosen to reduce the transient term. This choice depends crucially on the availability of efficient and reliable reduction algorithms (not all of which are applicable to systems of very large dimension [1] like the one in the following Subsection 4.3). The results of the suggested approximation procedure are then compared with those afforded by the balanced truncation and Hankel-norm approximation methods applied directly to the original system without consideration of the asymptotic response to a ramp input. Note, however, that the reduced model has been adjusted by truncating the original balanced realisation in such a way that the steady-state value in the response to a step input is retained exactly. This adjustment obviously entails that the reduced model is no longer a perfect truncation of the original balanced realisation. The second and the third of the following three examples are taken, rather arbitrarily, from the collection of high-order benchmark examples [25]. The first one, instead, considers an original system of rather low dimension just to show in detail how the method operates. However, according to Remark 4 of Section 3 the suggested approach makes more sense when the orders of both the original system and the reduced model are substantially greater than the order of the input whose asymptotic response must be retained.

Elucidative example
Assume that a second-order model of the system described by the 4th-order stable transfer function is required, and that the reduced model must retain the steady-state value of the step response (equal to the lowfrequency or DC gain, which is 1 in this case). Since no resonant component is present (C(s) = 1), the Laplace transform of the forced response of (7) to a unit step input decomposes as , which is impossible since r 1 (s) should be equal to B r (s), whose degree is 1, and r 2 (s) should be equal to q(s), whose degree is 2. As shown in Remark 3, Y a (s) must be of the same order as U(s), i.e., 1, and strictly proper, which in the present case implies that its numerator must be a constant. From (4) it follows that also the order of Y w r (s) must be 1. The first-order Hankelnorm approximation of the original transient term Y w (s) turns out to be Y w r (s) = −1.046/(s + 1.119). Step 2) of the procedure entails the determination of the overall forced response (4) of the reduced model to the chosen input (unit step), which in turn requires choosing the (stable) far-off pole p a of the auxiliary term Y a (s). This choice is not critical, being it enough that the pole of Y a (s) lies sufficiently to the left of the pole of Y w r (s) which is at Step 3) of the procedure entails forming the reducedorder transfer function by diving Y f r (s) by the input transform according to (5), which leads to W r (s) = −0.2095s + 33.57 (s + 1.119)(s + 30) . Fig. 1 shows the step responses of: (i) the original system with transfer function W (s) given by (7) (OM), (ii) the second-order model whose transfer function W H (s) is obtained by minimising the Hankel norm of the difference between W (s) and the reduced transfer function (which does not preserve the steady-state value) (HN), and (iii) the second-order model with transfer function W r (s) given by (8) (AR). Fig. 2 shows the magnitude Bode plots for the same models.

Hospital building
This example, whose original equations can be found in and [25], describes the dynamics of a hospital building with 8 floors, each having 3 degrees of freedom (horizontal and vertical displacements and rotation). The order of  1 (t)). The output of interest coincides with the derivative of the first coordinate (i.e., the 25th state x 25 (t)). Fig. 3 depicts the original model (OM) response to a unit step input as well as the step responses of the 6thorder models obtained using Hankel-norm approximation (HN), balanced truncation (with and without DC gain adjustment) (BY and BN, respectively) and the method (AR) outlined in Section 3 The last model has been obtained from the 4th-order approximation of the transient component in the ramp response; of course, the approximation would be better if a 6th-order approximation were used for this purpose, but the resulting reduced-order transfer function would then be of order 8.
Figs. 4 and 5 compare, respectively, ramp responses and Bode diagrams. Observe that, although the ramp response of the balanced-truncation model adjusted to retain the DC gain seems to coincide asymptotically with the original ramp response, it slightly differs from it for t → ∞. Clearly, the frequency responses of the reduced models obtained by minimising the Hankel norm of the error (HN) and by retaining the asymptotic ramp response (AR) are much closer to the original one at high frequencies since they are strictly proper like the original system.

1006th-order system
A 1006th-order model has been described in [26] and detailed in [25]. This model, which has been considered also in [5] to show the remarkable efficiency of a reduction method based on orthogonal polynomial expansions in the time domain, does not correspond to a real system but represents an interesting "purely theoretical text exam-  ple" [26] because its magnitude Bode plot is quite jagged.
Although the ramp response of model BY seems to coincide with the original ramp response, it actually exhibits a small asymptotic error equal to 0.074. Fig. 7 shows the original Bode diagrams together with the Bode diagrams of all of the aforementioned 9th-order models. Clearly, forcing the retention of the asymptotic response (t → ∞) to canonical inputs (steps and ramps) leads to a better approximation in the low frequency range (ω → 0) as is the case with the classic Padé approximation. Note, finally, that, the magnitude plot of the gain-adjusted model BY diverges from the original magnitude plot for ω → ∞ because its transfer function is not strictly proper. Fig. 6 depicts the original model response to a ramp input (OM) and the ramp responses of the 9th-order models obtained with the Hankel-norm approximation method (HN), the balanced truncation method (BY) (modified to match the DC gain) and the method (retaining the asymp- totic ramp response) (AR) outlined in Section 3 The latter has been formed from the 7th-order approximation of the transient component of the ramp response.

CONCLUSIONS
It has been shown how model reduction techniques can be adapted to reproduce the asymptotic response to inputs with rational transform. The suggested procedure is not computationally demanding and entails only a small increment of the reduced model order with respect to the order of the approximating transient term. Benchmark examples have shown that the responses of the reduced-order models obtained in this way compare favourably with those afforded by conventional reduction techniques that do not take care ab initio of the asymptotic response.