State-response decomposition for model reduction

ABSTRACT The decomposition of the overall forced response into a steady-state and a transient component can be exploited to find reduced-order models that retain both long-term and short-term characteristics of the original system behaviour. To this purpose, the state-space expressions of the aforementioned components in response to inputs with rational transform are derived. The reduced-order model is then obtained by considering separately the asymptotic and transient terms. The suggested approach is tested on benchmark examples of very high dimension.


Introduction
Recently, model reduction methods based on the decomposition of the forced response into a transient and a steady-state component (Dorato, Lepschy, and Viaro, 1994) have been proposed with the purpose of generating reduced-order models that retain, at least approximately, some important characteristics of both the transient component, for example, its impulseresponse energies and moments, and the asymptotic component in response to particular inputs. Also the interpolation-based moment-matching methods (Antoulas, Beattie, and Gugercin, 2010;Astolfi, 2010;Bultheel and De Moor, 2000) retain the steady-state response to particular inputs but do not even ensure, in their standard version, the preservation of stability. To combine transient features with the retention of at least the static or direct-current (DC) gain and, thus, of the steady-state value in the step response, a variant of the balancedtruncation method (Moore, 1981) based on perturbation theory (Kokotovic, Khail, and Reilly, 1986) has been suggested, which can be implemented using popular control packages (Chaturvedi, 2009). No such (limited) extensions are as yet available for other optimal reduction techniques, such as the Hankel-norm (Glover, 1984) and L 2 -norm (Antoulas, 2005;Krajewski and Viaro, 2009;Xu and Zeng, 2011;Zeng and Lu, 2015) approximation methods. Not even the newly proposed metaheuristic model reduction methods (Desai and Prasad, 2013;Ramawat and Kumar, 2015;Rana, Prasad, and Singh, 2014; CONTACT Daniele Casagrande daniele.casagrande@uniud.it Sikander and Prasad, 2015) explicitly address the problem of asymptotic accuracy. The reconciliation of optimal or suboptimal approximation methods with specifications on the asymptotic behaviour could be done in the domain of Laplace transforms (transfer functions). Unfortunately, the transferfunction approach is not suited for the reduction of very high-dimensional systems given in state-space form because the conversion from state-space representation to transfer-function is not numerically robust and cannot be safely applied to systems whose order is greater than a few tens. On the other hand, the test inputs often differ from singularity inputs and can involve, to more advantage, combinations of (possibly complex) exponentials, for example, sinusoids.
To overcome these problems, this paper first expresses the components of the forced response to inputs with arbitrary rational transform in terms of the original state-space representation (Section 2): since, in the time domain, these inputs are linear combinations of functions of the form t β e γ t , with β ∈ N and γ ∈ C, the state response to this kind of functions is preliminarily determined (Proposition 2.1). The reduced-order model is finally obtained by adding an approximation of the transient component to the original steady-state component (Section 3). Essentially, the contribution of this paper consists in the direct time-domain decomposition of the state response and in the separate treatment of its two components. The suggested approach is tested on a pair of high-order benchmark examples taken from the relevant literature (Chahlaoui and Van Dooren, 2005) (Section 4). Possible extensions of the method, with particular reference to its applicability to unstable system, are briefly outlined in Section 5.

Decomposition of the forced response
In the following, it is assumed that the original system is linear time-invariant strictly proper and asymptotically stable. Denoting its (minimal) state-space representation by the triplet = {A, B, C}, its forced state response to a causal input u(t) is given by As is known, if the Laplace transform U(s) of u(t) is real rational and strictly proper with poles in the closed right half-plane, u(t) can be expressed as a linear combination of real non-decreasing exponentials, possibly multiplied by powers of t and/or, in the presence of pairs of conjugate complex poles, sinusoidal functions. If U(s) is exactly proper (Bernstein, 2009), that is, proper but not strictly, an impulsive term is also present, but this case is excluded in the sequel for simplicity. Attention will be limited here to such persistent inputs because this is precisely the kind of signals (including singularity inputs and sinusoids) usually employed to test the system behaviour. Note that, since the system is bounded-input bounded-output stable, no resonance phenomena (interaction between input and system modes Dorato et al., 1994) can occur. By denoting the (possibly complex) poles of a strictly proper U(s) by q i , the corresponding input u(t) can also be written as a combination of functions that are not necessarily real as where n d is the number of distinct poles of u(t), μ i denotes the multiplicity of pole q i with Re(q i ) ≥ 0 and Q i,j ∈ C. The combination coefficients Q ij in (2) are real whenever q i is so, and complex otherwise. Also, since the coefficients of U(s) are assumed to be real, if q i is a complex pole of U(s), then its conjugateq i is also a pole of U(s), and the combination coefficient of every function t j e q i t associated with q i is conjugate to the coefficient of the homologous function t j eq i t associated withq i . Under these conditions, in fact, the overall input (2) is real for all t. In the sequel, reference is made to this expression of the input because it contains functions of one type only, namely t j e q i t . By linearity, the forced state response to the general input (2) is given by a linear combination of the responses to these basic signals. The following result, whose conceptually simple proof can be found in the Appendix, provides a useful tool for finding the forced response to an input of the general form (2).

Proposition 2.1:
The forced state response of the asymptotically stable system = {A, B, C} to the causal input is given by Remark 2.1: The invertibility of the matrix q i I − A is guaranteed because the eigenvalues of the asymptotically stable matrix A lie in the open left half-plane while q i lies in the closed right half-plane.

Remark 2.2:
The first term on the right-hand side of relation (4) contains the exponential function e At , which depends on the system modes only. Therefore, this term represents the transient component x tr (t) of the forced response to the considered input. Instead, the addenda under the summation symbol depend only on the input modes associated with the pole q i of the transform U q i ,j (s) = j!/(s − q i ) j+1 of (3) and together form the asymptotic (or steady-state) component x ss (t) of the response.

Remark 2.3:
As Equation (4) shows, in general the asymptotic component y ss (t) of the overall forced output response, which is related to the overall forced state response via is not simply proportional to the input but is a linear combination of the input and its first derivatives.
It is easy to derive from (4) the forced responses to the singularity inputs, which are reported next for convenience. Precisely, the response to u 0,j (t) = t j , t ≥ 0, turns out to be In particular, the step and ramp responses are, respectively, The response to can be obtained by combining the responses to e jωt and e −jωt according to the conjugate coefficients 1/2j and −1/2j . After some trivial manipulations, this response turns out to be whose transient and steady-state components can be identified easily.

Model reduction
Model reduction is, by its very nature, an input-output problem. Therefore, in the following, attention focuses on the forced output response y f (t), simply related to x f (t) via (5) According to the approach outlined in Section 1, the suggested reduction procedure can be presented as follows. This procedure is schematically represented in Figure 1.  (5) that holds for the overall forced responses. The second matrix B tr of the triplet tr , instead, depends on both A and B. For instance, in the case of the ramp input, it is given by B tr = A −2 B (see Equation (8) and in that of the sinusoidal input (9) by B tr = (ω 2 I + A 2 ) −1 Bω (see Equation (10)).

Remark 3.2:
Step (ii) of Procedure 3.1 leads to an approximate realization a tr = {A a , B a tr , C a } of tr . Since the dimension of a tr is much smaller than that of the original system, the conversion from its state-space representation to its s-domain representation does not pose any numerical problem.
Step (iii) of Procedure 3.1 is not so straightforward as it might seem at first sight. To illustrate this point, let us refer to the s-domain representation of the reduced model, and denote by Y a tr (s) the Laplace transform of the reduced-order transient component, which is completely determined after Step (ii), and by Y ss (s) the Laplace transform of the (known) original asymptotic component to be preserved. The transfer function W a (s) of the reducedorder model should satisfy the equation However, since the denominator d a w (s) of the strictly proper W a (s) must coincide with the denominator d a tr (s) of Y a tr (s), the only unknowns in (11) are the coefficients of the numerator n a w (s) of W a (s) whose number is equal to deg[d a tr (s)], whereas the number of equations that may be formed by equating the numerator coefficients of the equal powers of s on both sides of (11)

Examples
This section shows the results of the application of the model reduction technique based on the suggested state-response decomposition to a pair of very high-order examples for which the transfer-function approach would be unreliable or even impracticable since round-off errors and very small changes in the input data produce wildly different s-domain representations. The time responses of the reduced-order models that retain the steady-state component of the original response to selected inputs are compared with those obtained, without explicit consideration of the asymptotic behaviour, by means of the popular Hankel-norm and balanced-truncation methods that can be implemented using readily available computer programs. Since the purpose of this paper is just to show how reduction techniques can be adapted to the case in which the asymptotic behaviour need be reproduced and not to propose a new approximation algorithm, other reduction techniques are not considered even if they might lead to better approximations. Of course, the improvement in the steady-state behaviour of the suggested approach is limited to the inputs of interest.

Example 1
Consider first the 348th-order model of a clamped beam illustrated in Chahlaoui and Van Dooren (2005). Figure 2 shows the ramp responses of: (i) the original system, (ii) the sixth-order model retaining the original steady-state response to the ramp input, (iii) the sixth-order model obtained using Hankel-norm approximation, and (iv) the sixth-order model obtained using balanced truncation (with DC gain adjustment, thus retaining the steady-state response to a step input).
The transient component of the model that retains the steady-state component has been approximated according to the Hankel-norm criterion. Figure 3 compares the responses to the sinusoidal input u(t) = sin 5t. In order to highlight the differences, attention is limited to the interval 190 ≤ t ≤ 200 s.

Example 2
Consider now the 48th-order model that describes the dynamics of a hospital building. This model, which is illustrated in Chahlaoui and Van Dooren (2005), has also been considered in Antoulas, Sorensen, and Gugercin (2001). Figure 4 compares again the ramp response of the original model with the ramp responses of the sixth-order models obtained by means of: (i) the proposed technique ensuring the retention of the asymptotic component, (ii) the Hankel-norm approximation, and (iii) the balancedtruncation method with DC gain adjustment. Observe, in particular, that the ramp response of the balancedtruncation model does not coincide asymptotically with the original response, as shown by the smaller picture inside the figure. Figure 5 compares instead the responses to the input u(t) = 5 sin 10t. In this case, attention is limited to the interval 196 ≤ t ≤ 200 s to make the differences more visible.

Conclusions
It has been shown how standard reduction techniques can be adapted to retain the steady-state response to a numerous class of inputs, namely those expressed by a linear combination of exponentials, possibly multiplied by powers of t. This result is achieved by increasing a little the order of the reduced model with respect to the order of the models obtained from the impulse response only. This increase is negligible when the dimension of the  original system is high, typically, many tens or even hundreds, as is the case in many practical applications (see, e.g. Chahlaoui and Van Dooren, 2005;Gawronski, 1998). In these cases, the resort to state-space representations is mandatory for numerical robustness. To this purpose, the components of the original forced response have preliminarily been expressed in this form.
Examples have shown that the approximation of the responses to the chosen test inputs improves significantly not only at steady state but also during the transients.
The same approach could be adopted to approximate unstable original systems by first separating the stable and unstable parts and then applying the reduction procedure to the first one only. The reduced model will finally be obtained by combining the reduced stable component with the original unstable one. In this case too, care must be taken to ensure that the reduced transferfunction denominator is the same as the denominator of the reduced stable transient term. This can again be obtained by including an auxiliary term as in (12).
The approach could also be extended to the case where the input denominator has a factor in common with the denominator of the original transfer function, which would require considering a resonant or interaction term (Dorato et al., 1994) besides the transient and asymptotic terms.

Disclosure statement
No potential conflict of interest was reported by the authors.
which is equivalent to (A2) for j = 0. Now, suppose that (A2) holds for a generic j. By applying the integration-by-parts rule, we obtain which cancels the second addendum in (A4). In addition, the first addendum can be included in the summation by simply extending it to k = −1. As a consequence, Finally, by setting h = k + 1, we obtain which proves the claim for j+1.