Relative Error-based Time-limited $\mathcal{H}_{2}$ Model Order Reduction via Oblique Projection

Computer-aided simulation and analysis have been a major driver in scientific innovation and progress for the last several decades. The dynamical systems and processes, which can be natural or artificial, are mathematically modeled to be used in a computer program that can simulate and analyze that system for endless possible inputs, settings, and scenarios. During the mathematical modeling process, the behavior of the system or process is generally described by partial differential equations (PDEs). To utilize the mathematical tools available in dynamical system theory, the PDEs are converted into ordinary differential equations (ODEs), which are then expressed as a state-space model [1, 2].

Let $u(t)\in\mathbb{R}^{m\times 1}$ be the inputs, $y(t)\in\mathbb{R}^{p\times 1}$ be the outputs, and $x(t)\in\mathbb{R}^{n\times 1}$ be the states of the state-space model of the dynamical system $H$. Then the state-space equations of $H$ can be written as

wherein $A\in\mathbb{R}^{n\times n}$, $B\in\mathbb{R}^{n\times m}$, $C\in\mathbb{R}^{p\times n}$, and $D\in\mathbb{R}^{p\times m}$. The state-space realization $(A,B,C,D)$ has the following equivalence with the transfer function of $H$

Note that the simulation of $H$ requires the solution of $n$ differential equations. The computing power of modern-day computers has been exponentially increasing following Moore’s law, and memory resources have been increasing tremendously [2]. However, the complexity of modern-day systems and processes has also been increasing tremendously following the scientific and technological advancements enabled by high-speed chip manufacturing. Therefore, the value of $n$ in most modern-day systems and processes is normally several thousand. The simulation and analysis of such a high-order model remain a computational challenge, which has been a consistent motivation for developing efficient model order reduction (MOR) algorithms [3, 4, 5]. MOR is a process of obtaining a reduced-order model (ROM) that has similar properties and behavior but can be simulated cheaply without exhausting the available computational resources. The ROM can then be used as a surrogate for the original model with tolerable numerical error.

Let us denote the ROM as $\hat{H}$. Further, let us denote the states and outputs of $\hat{H}$ as $\hat{x}(t)$ and $\hat{y}(t)$, respectively. Then the state-space equations of $\hat{H}$ can be written as

wherein $\hat{A}\in\mathbb{R}^{r\times r}$, $\hat{B}\in\mathbb{R}^{r\times m}$, $\hat{C}\in\mathbb{R}^{p\times r}$, $\hat{D}\in\mathbb{R}^{p\times m}$, and $r\ll n$. The state-space realization $(\hat{A},\hat{B},\hat{C},\hat{D})$ has the following equivalence with the transfer function of $\hat{H}$

Ideally, the MOR algorithm should have three main properties:

1. 1.

   $\hat{H}$ should closely mimic $H$ in the frequency and/or time domains.

2. 2.

   $\hat{H}$ should preserve important system properties like stability, passivity, contractivity, etc.

3. 3.

   The computational cost required to obtain $\hat{H}$ should be way less than that required to simulate $H$.

No MOR algorithm has all these properties; nevertheless, the significance of a MOR algorithm is generally accessed against these three standards. Mathematically, two possible ways to express the MOR problem are the following:

1. 1.

   $H=\hat{H}+\Delta_{add}$.

2. 2.

   $H=\hat{H}(I+\Delta_{rel})$.

Most of the MOR algorithms use additive error $\Delta_{add}=H-\hat{H}$ as their performance criterion, i.e., they minimize $\Delta_{add}$. The quality of approximation in these algorithms is accessed by computing system norms (like $\mathcal{H}_{2}$ norm and $\mathcal{H}_{\infty}$ norm) of $\Delta_{add}$. If the output responses of $H$ and $\hat{H}$ are inherently small, a small output response of $\Delta_{add}$ can be a misleading notion of accuracy. The relative error $\Delta_{rel}$ is a better approximation criterion as it reveals percentage error. However, $\Delta_{rel}$ is a less used approximation criterion because most of the algorithms that use this criterion are computationally expensive due to the inherent nature of the problem, which will be discussed later. Various performance specifications and specific properties to be preserved lead to various MOR families. Within the family, several algorithms with different discourses are available, each having some distinct features. In the sequel, a brief literature survey of the important MOR algorithms is presented.

Balanced truncation (BT) is one of the most famous MOR techniques due to its useful theoretical properties, like stability preservation and supreme accuracy [6]. BT has led to several important MOR algorithms constituting a rich family of algorithms known as balancing- or gramian-based MOR techniques in the literature. The ROM produced by BT is often suboptimal in the $\mathcal{H}_{\infty}$ norm, i.e., a suboptimum of $||\Delta_{add}(s)||_{\mathcal{H}_{\infty}}$. An apriori upper bound on $||\Delta_{add}(s)||_{\mathcal{H}_{\infty}}$ for the ROM generated by BT was proposed in [7]. Most of the algorithms in the BT family offer a similar apriori upper bound, which is a distinct and important feature of balancing-based MOR algorithms. BT is a computationally expensive algorithm as it requires the solution of two large-scale Lyapunov equations, which limits its applicability to medium-scale systems. To overcome this problem, several numerical algorithms are reported that compute low-rank solutions of these Lyapunov equations, which are cheaper to compute [8, 9, 10, 11]. By using these low-rank solutions, the applicability of BT is extended to large-scale systems [12]. Some other numerical methods to avoid ill-conditioning in the BT algorithm are also reported, like [13, 14]. BT has been extended to preserve other system properties like passivity and contractivity in [15, 16, 17] and [18], respectively. Further, its applicability has also been extended to more general classes of linear and nonlinear systems, cf.[19, 20, 21, 22]. A closely related yet different MOR method known as optimal Hankel norm approximation (OHNA) is proposed in [23]. The ROM produced by OHNA is optimal in the Hankel norm, and an apriori upper bound on $||\Delta_{add}(s)||_{\mathcal{H}_{\infty}}$, which is half of that in the case of BT, also holds [24]. An important extension to BT, which is relevant to the problem considered in this paper, is the balanced stochastic truncation (BST) [25]. It minimizes $||\Delta_{rel}(s)||_{\mathcal{H}_{\infty}}$, and an apriori upper bound on $||\Delta_{rel}(s)||_{\mathcal{H}_{\infty}}$ also exists [26]. If the original model is a stable minimum phase system, BST guarantees that the ROM preserves this property.

In a large-scale setting, the computation of the $\mathcal{H}_{\infty}$ norm is expensive. Moreover, the MOR algorithms that tend to reduce the $\mathcal{H}_{\infty}$ norm of error are generally also expensive [27, 28, 29]. The $\mathcal{H}_{2}$ norm, on the other hand, can be computed cheaply due to its relation with the system gramians, for which several efficient low-rank algorithms are available. Interestingly, the MOR algorithms that seek to ensure less error in the $\mathcal{H}_{2}$ norm happen to be the most computationally efficient ones in the available literature [30]. The selection of the $\mathcal{H}_{2}$ norm as a performance specification in this paper is motivated by this factor. A locally optimal solution in the $\mathcal{H}_{2}$ norm is also possible, and several efficient algorithms are available that achieve the local optimum upon convergence. The necessary conditions for a local optimum of $||\Delta_{add}(s)||_{\mathcal{H}_{2}}^{2}$ (known as Wilson’s conditions) are related to interpolation conditions, i.e., $\hat{H}(s)$ interpolates $H(s)$ at selected points in the $s$-plane, cf. [31, 32, 33, 34]. The iterative rational Krylov algorithm (IRKA) is the pioneer and most famous interpolation algorithm for single-input single-output (SISO) systems that can achieve the local optimum of $||\Delta_{add}(s)||_{\mathcal{H}_{2}}^{2}$ upon convergence [32]. IRKA is generalized for multi-input multi-output (MIMO) systems in [33]. A more general algorithm that does not assume that $H(s)$ and $\hat{H}(s)$ have simple poles is presented in [34] that can capture the local optimum of $||\Delta_{add}(s)||_{\mathcal{H}_{2}}^{2}$ upon convergence. This algorithm is based on the Sylvester equations and uses the connection between the Sylvester equations and projection established in [35]. Its numerical properties are further improved in [36], and a more robust algorithm is presented. To speed up convergence in IRKA, some trust region-based techniques are proposed in [37]. Some other $\mathcal{H}_{2}$-optimal algorithms based on manifold theory are also reported, like [38, 39, 40, 41]. To guarantee stability, some suboptimal solutions to $||\Delta_{add}(s)||_{\mathcal{H}_{2}}^{2}$ are reported, like [42, 43, 44].

The simulation of a model does not encompass the entire time horizon since no practical system is run over the entire time range. Thus it is reasonable to ensure good accuracy in the actual time range of operation, which is often limited to a small interval. For instance, low-frequency oscillations in interconnected power systems only last for 15 seconds, after which these are successfully damped by the power system stabilizers and damping controllers. The first 15 seconds thus become important in small-signal stability analysis of interconnected power systems [45, 46, 47]. Similarly, in the time-limited optimal control problem, the behavior of the plant in the desired time interval is important [48]. This motivates time-limited MOR, wherein it is focused to ensure supreme accuracy in the desired time interval rather than trying to achieve good accuracy over the entire time range. BT was generalized for the time-limited MOR problem in [49] to obtain a time-limited BT (TLBT) algorithm. TLBT does not retain the stability preservation and apriori error bound properties of BT. However, an $\mathcal{H}_{2}$ norm-based upper bound on $\Delta_{add}$ is proposed in [50, 51] for TLBT. Some ad hoc approaches to guarantee stability in TLBT are also reported, like [52]; however, these are quite inaccurate and computationally expensive. TLBT, like BT, is computationally expensive in a large-scale setting and requires solutions of two large-scale Lyapunov equations. In [53], these Lyapunov equations are replaced with their low-rank solutions to extend the applicability of TLBT to large-scale systems. TLBT is further generalized to preserve more properties like passivity and to extend its applicability to a more general class of linear and nonlinear systems in the literature, cf. [54, 55, 56, 57, 58]. BST is extended to the time-limited MOR scenario in [59]; however, an upper bound on $\Delta_{rel}$ does not exist for time-limited BST (TLBST).

In [60], a definition for the time-limited $\mathcal{H}_{2}$ norm is presented, and the necessary conditions for a local optimum of $||\Delta_{add}(s)||_{\mathcal{H}_{2,\tau}}^{2}$ (wherein $\mathcal{H}_{2,\tau}$ represents the time-limited $\mathcal{H}_{2}$ norm) are derived. Then, IRKA is heuristically extended to time-limited IRKA (TLIRKA) for achieving a local optimum of $||\Delta_{add}(s)||_{\mathcal{H}_{2,\tau}}^{2}$. However, TLIRKA does not satisfy the necessary conditions for a local optimum. In [61], interpolation-based necessary conditions for a local optimum of $||\Delta_{add}(s)||_{\mathcal{H}_{2,\tau}}^{2}$ are derived, and a nonlinear optimization-based algorithm is proposed to construct the local optimum. This algorithm, however, is computationally expensive, and its applicability is limited to the order of a few hundred. Another nonlinear optimization-based algorithm for this problem is reported in [62], which is very similar to that given in [61]. The equivalence between gramian-based conditions and interpolation conditions is established in [63, 64]. A stability-preserving interpolation algorithm that achieves a subset of the necessary conditions is proposed in [65]. To the best of the authors’ knowledge, there is no algorithm in the literature that seeks to reduce $||\Delta_{rel}(s)||_{\mathcal{H}_{2,\tau}}$, which has motivated the results obtained in this paper.

In this paper, necessary conditions for the local optimum of $||\Delta_{rel}(s)||_{\mathcal{H}_{2,\tau}}^{2}$ are derived. It is shown that it is inherently not possible to achieve a local optimum of $||\Delta_{rel}(s)||_{\mathcal{H}_{2,\tau}}^{2}$ in an oblique projection framework. Nevertheless, inspired by these necessary conditions, a computationally efficient oblique projection-based algorithm is proposed that seeks to obtain the local optimum of $||\Delta_{rel}(s)||_{\mathcal{H}_{2,\tau}}^{2}$ and offers high fidelity. Unlike TLBST, the proposed algorithm does not require the solutions of large-scale Lyapunov and Ricatti equations. The proposed algorithm is numerically compared to TLBT, TLBST, and TLIRKA by considering benchmark examples. The numerical results confirm that the proposed algorithm offers a small relative error within the desired time interval.

Let the original high-order system be an $n$-th order $m\times m$ square linear time-invariant system $H$ and $(A,B,C,D)$ be its state-space realization. Further, let us denote the impulse response of $H$ as $h(t)=Ce^{At}B$. The time-limited controllability gramian $P=\int_{0}^{t_{d}}e^{A\tau}BB^{T}e^{A^{T}\tau}d\tau$ and time-limited observability gramian $Q=\int_{0}^{t_{d}}e^{A^{T}\tau}C^{T}Ce^{A\tau}d\tau$ of the state-space realization $(A,B,C)$ within the desired time interval $[0,t_{d}]$ sec solve the following Lyapunov equations

cf. [49].

In this section, first, a state-space realization and an expression for the $\mathcal{H}_{2,\tau}$ norm of $\Delta_{rel}(s)$ are obtained. Second, necessary conditions for the local optimum of $||\Delta_{rel}(s)||_{\mathcal{H}_{2,\tau}}^{2}$ are derived. Third, an oblique projection algorithm is presented, which seeks to satisfy these conditions. Fourth, the algorithmic and computational aspects of the proposed algorithm are discussed.