$\mathcal{H}_{2}$-optimal Model Reduction of Linear Quadratic Output Systems in Finite Frequency Range

This study focuses on a specific class of nonlinear dynamical systems with weak nonlinearity. These systems have linear time-invariant (LTI) state equations but feature quadratic nonlinear terms in the output equation, referred to as linear quadratic output (LQO) systems [1]. They naturally arise in scenarios where it is necessary to observe quantities involving products of state components in either the time or frequency domain. For example, they are used in applications that quantify energy or power, such as the internal energy functional of a system [2] or the objective cost function in optimal quadratic control problems [3]. Additionally, they are used to measure the deviation of a state’s coordinates from a reference point, such as the root mean squared displacement of spatial coordinates around an excitation point, or in stochastic modeling for calculating the variance of a random variable [4].

Consider a LQO system described by the following state and output equations:

wherein $x(t)\in\mathbb{R}^{n\times 1}$ represents the state vector, $u(t)\in\mathbb{R}^{m\times 1}$ represents inputs, and $y(t)\in\mathbb{R}^{p\times 1}$ represents outputs. $(A,B,C,M_{1},\cdots,M_{p})$ is the state-space realization of $G$ with $A\in\mathbb{R}^{n\times n}$, $B\in\mathbb{R}^{n\times m}$, $C\in\mathbb{R}^{p\times n}$, and $M_{i}\in\mathbb{R}^{n\times n}$. The state equation in (1) is identical to that of standard LTI systems. However, the output equation in (1) introduces a nonlinearity in the form of the quadratic function of states $x(t)^{T}M_{i}x(t)$.

Let us denote $y_{1}(t)=Cx(t)$ and $y_{2,i}=x(t)^{T}M_{i}x(t)$. The input-output mapping between $u(t)$ and $y_{1}(t)$ is represented by the following transfer function:

Additionally, the input-output mapping between $u(t)$ and $y_{2,i}(t)$ is represented by the following multivariate transfer function:

cf. [5].

To ensure high fidelity in the mathematical modeling of complex physical phenomena, it is often necessary to use dynamical systems with very high orders, sometimes exceeding several thousand. This high order $n$ makes it computationally difficult or even prohibitive to simulate and analyze the model (1). Therefore, it is important to approximate (1) with a reduced-order model (ROM) of much lower order $k$ (where $k\ll n$). This approach simplifies the simulation and analysis processes. The process of creating such a ROM while preserving the key features of the original model is known as model order reduction (MOR); refer to [6, 7, 8, 9, 10] for a deeper understanding of the topic.

Let us denote the $k^{th}$-order ROM of $G$ as $G_{k}$, characterized by the following state and output equations:

wherein $A_{k}=W_{k}^{T}AV_{k}\in\mathbb{R}^{k\times k}$, $B_{k}=W_{k}^{T}B\in\mathbb{R}^{k\times m}$, $C_{k}=CV_{k}\in\mathbb{R}^{p\times k}$, and $M_{k,i}=V_{k}^{T}M_{i}V_{k}\in\mathbb{R}^{k\times k}$, satisfying the Petrov-Galerkin projection condition $W_{k}^{T}V_{k}=I$. The projection matrices $V_{k}\in\mathbb{R}^{n\times k}$ and $W_{k}\in\mathbb{R}^{n\times k}$ project $G$ onto a reduced subspace to obtain the ROM $G_{k}$. Various MOR methods differ in how they construct $V_{k}$ and $W_{k}$. The choice of $V_{k}$ and $W_{k}$ depends on the specific characteristics of $G$ that need to be preserved in $G_{k}$.

Let $y_{k,1}(t)=C_{k}x_{k}(t)$ and $y_{k,2,i}=x_{k}(t)^{T}M_{k,i}x_{k}(t)$. The input-output relationships between $u(t)$ and $y_{k,1}(t)$, as well as $u(t)$ and $y_{k,2,i}(t)$, are described by the following transfer functions:

Throughout this paper, it is assumed that both $A$ and $A_{k}$ are Hurwitz.

The Balanced Truncation (BT) method, introduced in 1981, is a prominent technique for MOR [11]. This approach retains states that significantly effect energy transfer between inputs and outputs, discarding those with minimal impact as indicated by their Hankel singular values. A notable advantage of BT is its ability to estimate errors in advance of creating the ROM [12]. Moreover, BT preserves the stability of the original system. Originally proposed for LTI systems, BT’s application has broadened to include descriptor systems [13], second-order systems [14], linear time-varying systems [15], parametric systems [16], nonlinear systems [17], and bilinear systems [18]. Additionally, BT has been tailored to maintain specific system properties such as positive realness [19], bounded realness [20], passivity [21], and special structural characteristics [22]. For a detailed overview of the various BT algorithms, refer to the survey [23]. BT has been adapted for LQO systems in [24, 25, 26]. Of these algorithms, only the one proposed in [26] preserves the LQO structure in the ROM.

The $\mathcal{H}_{2}$-optimal MOR problem for standard LTI systems has been thoroughly explored in the literature. Wilson’s conditions, which are necessary conditions for achieving a local optimum of the $\mathcal{H}_{2}$ norm of error, are outlined in [27]. The iterative rational Krylov algorithm (IRKA) was the first to apply interpolation theory to meet these conditions [28]. Other algorithms using Sylvester equations and projection have been developed and enhanced for improved robustness in [29] and [30]. Recently, the $\mathcal{H}_{2}$-optimal MOR problem for LQO systems has been addressed in [31], where a Sylvester equation-based algorithm has been proposed to achieve a local optimum upon convergence.

Many MOR problems are inherently frequency-limited, with certain frequency ranges being more important. For example, when creating a ROM for a notch filter, it is crucial to minimize the approximation error near the notch frequency [32]. Similarly, for closed-loop stability, the ROM of a plant must accurately capture the system’s behavior in the crossover frequency region [33, 34]. In interconnected power systems, low-frequency oscillations are essential for small-signal stability studies, so the ROM should accurately represent the behavior within the frequency range of inter-area and inter-plant oscillations [35, 36]. This need has led to frequency-limited MOR, which focuses on achieving high accuracy within specific frequency intervals rather than the entire spectrum [37].

The frequency-limited MOR problem aims to create a ROM that ensures that

are small when $\nu$, $\nu_{1}$, and $\nu_{2}$ are within the specified frequency range of $[0,\omega]$ rad/sec.

BT typically offers an accurate approximation of the original model across the entire frequency spectrum. In [37], BT was adapted to address the frequency-limited MOR problem, resulting in the frequency-limited BT (FLBT) algorithm. However, FLBT does not preserve the stability and a priori error bounds that BT does. The computational aspects of FLBT and efficient methods for handling large-scale systems are discussed in [38]. Additionally, FLBT has been expanded to cover a wider range of systems, including descriptor systems [39], second-order systems [40], and bilinear systems [41]. FLBT has been recently generalized for LQO systems in [42].

The frequency-limited $\mathcal{H}_{2}$ norm for LTI systems is defined in [43], with Gramian-based conditions outlined for achieving a local optimum. Generalizations of the iterative rational Krylov algorithm (IRKA) have resulted in the development of the frequency-limited IRKA (FLIRKA) in [44, 45, 46], which approximately fulfills these conditions. This paper examines the frequency-limited $\mathcal{H}_{2}$-optimal MOR problem for LQO systems.

This research work makes several important contributions. First, it introduces the frequency-limited $\mathcal{H}_{2}$ norm ($\mathcal{H}_{2,\omega}$ norm) for LQO systems and shows how to compute it using frequency-limited system Gramians defined in [42]. Second, it derives the necessary conditions for achieving a local optimum of $||G-G_{k}||_{\mathcal{H}_{2,\omega}}^{2}$. Third, it compares these conditions with those for standard $\mathcal{H}_{2}$-optimal model order reduction [31], highlighting that Petrov-Galerkin projection generally cannot achieve a local optimum in the frequency-limited scenario. Fourth, it proposes a stationary point algorithm based on Petrov-Galerkin projection, which meets three of the four necessary conditions for optimality upon convergence. The paper includes a numerical example demonstrating the algorithm’s accuracy within the specified frequency range and showing it outperforms existing methods.

In this section, we will briefly explore two key MOR algorithms relevant to LQO systems within the context of the problem at hand. The first is the FLBT [42], and the second is the $\mathcal{H}_{2}$-optimal MOR method [31].

In this section, we define the frequency-limited $\mathcal{H}_{2}$ norm and establish its connection to the frequency-limited observability Gramian. We then derive the necessary conditions for achieving a local optimum of the (squared) frequency-limited $\mathcal{H}_{2}$ norm of the error. Building on these optimality conditions, we present a projection-based iterative algorithm that meets three of the four optimality conditions. The challenge of satisfying the fourth optimality condition within the projection framework is also discussed. Finally, the computational aspects of the proposed algorithm are briefly discussed.