Model Boundary Approximation Method as a Unifying Framework for Balanced Truncation and Singular Perturbation Approximation

TL;DR

This paper unifies two major model reduction techniques, Balanced Truncation and Singular Perturbation, within a framework derived from the Model Boundary Approximation Method, enabling interpolation and potential extension to nonlinear systems.

Contribution

It demonstrates that these techniques are boundary cases on a model manifold and introduces a method to interpolate between them, suggesting extensions to nonlinear systems.

Findings

Balanced Truncation and Singular Perturbation are boundary points on a model manifold.

MBAM allows interpolation between the two reduction techniques.

Potential extension of these methods to nonlinear systems.

Abstract

We show that two widely accepted model reduction techniques, Balanced Truncation and Balanced Singular Perturbation Approximation, can be derived as limiting approximations of a carefully constructed parameterization of Linear Time Invariant (LTI) systems by employing the Model Boundary Approximation Method (MBAM), a recent development in the Physics literature. This unifying framework of these popular model reduction techniques shows that Balanced Truncation and Balanced Singular Perturbation Approximation each correspond to a particular boundary point on a manifold, the "model manifold," which is associated with the specific choice of model parameterization and initial condition, and is embedded in a sample space of measured outputs, which can be chosen arbitrarily, provided that the number of samples exceeds the number of parameters. We also show that MBAM provides a novel way to interpolate between Balanced Truncation and Balanced Singular Perturbation Approximation, by exploring the set of approximations on the boundary of the manifold between the elements that correspond to the two model reduction techniques; this allows for alternative approximations of a given system to be found that may be better under certain conditions. The work herein suggests similar types of approximations may be obtainable in topologically similar places (i.e. on certain boundaries) on the model manifold of nonlinear systems if analogous parameterizations can be achieved, therefore extending these widely accepted model reduction techniques to nonlinear systems.

Equations

$$\begin{array}{rl}\dot{x}(t)& =Ax(t)+Bu(t)\\ y(t)& =Cx(t)+Du(t)\end{array}$$

$$G(s)=C(sI-A{)}^{-1}B+D$$

$$\begin{array}{rl}\dot{\bar{x}}(t)& =\bar{A}\bar{x}(t)+\bar{B}u(t)\\ y(t)& =\bar{C}\bar{x}(t)+\bar{D}u(t)\end{array}$$

$$\begin{array}{rl}{\bar{A}}^{T}X+X\bar{A}& =-{\bar{C}}^T\bar{C}\text{ and}\\ \bar{A}X+X{\bar{A}}^T& =-\bar{B}{\bar{B}}^T\end{array}$$

$$\begin{array}{rl}\bar{A}& =\left[\begin{array}{c}{\bar{A}}_{11}{\bar{A}}_{12}\\ {\bar{A}}_{21}{\bar{A}}_{22}\end{array}\right],\\ \bar{B}& =\left[\begin{array}{c}{\bar{B}}_1\\ {\bar{B}}_2\end{array}\right],\\ \bar{C}& =\left[\begin{array}{c}{\bar{C}}_1{\bar{C}}_2\end{array}\right],\text{ and}\\ \bar{x}& =\left[\begin{array}{c}{\bar{x}}_1\\ {\bar{x}}_2\end{array}\right]\end{array}$$

$$\begin{array}{rl}{\dot{\bar{x}}}_1(t)& ={\bar{A}}_{11}{\bar{x}}_1(t)+{\bar{B}}_{1}u(t)\\ y(t)& ={\bar{C}}_1{\bar{x}}_1(t)+Du(t).\end{array}$$

$${\bar{G}}_r(s)={\bar{C}}_1(sI-{\bar{A}}_{11}{)}^{-1}{\bar{B}}_1+D.$$

$$\Vert G-{\bar{G}}_r{\Vert }_{\infty }\le 2{\theta }_{n-k+1}.$$

$$\Vert G-{\bar{G}}_r{\Vert }_{\infty }\le 2\sum _{i=n-k+1}^n{\theta }_i.$$

$$\begin{array}{rl}{\dot{\hat{x}}}_1(t)& ={\hat{A}}_{11}{\bar{x}}_1(t)+{\hat{B}}_{1}u(t)\\ y(t)& ={\hat{C}}_1{\hat{x}}_1(t)+\hat{D}u(t)\end{array}$$

$$\begin{array}{rl}\hat{A}& ={\bar{A}}_{11}-{\bar{A}}_{12}{\bar{A}}_{22}^{-1}{\bar{A}}_{21},\\ \hat{B}& ={\bar{B}}_1-{\bar{A}}_{12}{\bar{A}}_{22}^{-1}{\bar{B}}_2,\\ \hat{C}& ={\bar{C}}_1-{\bar{C}}_2{\bar{A}}_{22}^{-1}{\bar{A}}_{21},\text{ and}\\ \hat{D}& =\bar{D}-{\bar{C}}_2{\bar{A}}_{22}^{-1}{\bar{B}}_2.\end{array}$$

$${\hat{G}}_r(s)=\hat{C}(sI-\hat{A}{)}^{-1}\hat{B}+\hat{D}.$$

$$\Vert G-{\hat{G}}_r{\Vert }_{\infty }\le 2{\theta }_{n-k+1}.$$

$$\Vert G-{\hat{G}}_r{\Vert }_{\infty }\le 2\sum _{i=n-k+1}^n{\theta }_i.$$

$$\dot{\mathbf{x}}$$

$$0$$

$$\mathbf{y}=\mathbf{h}(\mathbf{x},\mathbf{z},\mathbf{p},\mathbf{u},t).$$

$$\frac{d}{dt}\frac{\partial \mathbf{x}}{\partial \mathbf{p}}$$

$$0$$

$$\frac{\partial \mathbf{y}}{\partial \mathbf{p}}$$

$$\frac{d^2{p}^i}{d{\tau }^2}=\sum _{j,k}{\Gamma }_{jk}^i\frac{d{p}^j}{d\tau }\frac{d{p}^k}{d\tau };{\Gamma }_{jk}^i=\sum _{l,m}={\left({\mathbf{I}}^{-1}\right)}^{il}\frac{\partial y_m}{\partial p_l}\frac{{\partial }^2{y}_m}{\partial p_j\partial p_k}$$

$$\dot{x}(t)=\frac{-{\rho }_{1}x(t)}{{\rho }_2+x(t)}$$

$$\dot{x}(t)$$

$$\begin{array}{rl}\left[\begin{array}{c}{\dot{x}}_1(t)\\ {\dot{x}}_2(t)\end{array}\right]& =\left[\begin{array}{cc}-{\rho }_1& 0\\ 0& -{\rho }_2\end{array}\right]\left[\begin{array}{c}{x}_1(t)\\ x_2(t)\end{array}\right]\\ y& =\left[\begin{array}{cc}1& 1\end{array}\right]\left[\begin{array}{c}{x}_1(t)\\ x_2(t)\end{array}\right]\end{array}$$

$$y(t)=e^{-{\rho }_{1}t}{x}_1(0)+e^{-{\rho }_{2}t}{x}_2(0).$$

$$Y(s)=G(s)U(s)$$

$$G(s)=C(sI-A{)}^{-1}B+D.$$

$$G(s)=\left[\begin{array}{ccc}\frac{{\nu }_{11}(s)}{\delta (s)}& \cdots & \frac{{\nu }_{1m}(s)}{\delta (s)}\\ ⋮& \ddots & ⋮\\ \frac{{\nu }_{p1}(s)}{\delta (s)}& \cdots & \frac{{\nu }_{pm}(s)}{\delta (s)}\end{array}\right].$$

$$G(s)=\frac{\nu (s)}{\delta (s)}$$

$$G(s)=D+\sum _{i=1}^N\frac{K_i}{s-{\lambda }_i}$$

$$n^2+nm+np+pm$$

Full text

I Introduction

Modern systems theorists are studying and engineering systems that are larger and more complex than ever before [2]. Typical examples of these complex systems include economic networks [3], biological systems [4, 5, 6], and the Internet [7]. Due to the unprecedented size of these systems, simplified models are necessary to reason about them effectively [8, 9]. Specifically, we detail four important motivations for building simplified approximations of large scale complicated systems, and then provide an overview of how this work relates singular perturbation and balanced truncation methods.

I-A Motivation

When attempting to learn a system from limited data, large models cannot be identified, making simplified models necessary. First principles models typically have many parameters that must be tuned correctly for the model to reflect the behavior of a real system. Examples are everywhere, from agronomy and biochemical reaction networks to swarms of autonomous flying robots or power systems. Using data to learn the correct values of parameters is the purview of system identification, and a rich theory has developed quantifying when data is informative enough to accurately estimate parameter values [10, 11, 12, 13]. Typically, however, there is much less data than needed to learn all the parameters in a first-principles model, so simplifying the model to yield one with fewer parameters can help identify the system from data.

Second, the need for simplified models arises when designing controllers for complex systems. The complexity of an optimal controller often mirrors that of the system being controlled; therefore a complex system may often suggest the need for complicated controllers. Nevertheless, when engineering such complicated systems is unreasonable, designing controllers for simplified approximations can lead to acceptable trade-offs between complexity and performance degradation.

Third, simplified models can be an important link between macro-scale and micro-scale models. Generative models, for example, often detail micro-scale phenomena, such as consumer-choice models or models of a single neuron or molecular organization, etc., and then hypothesize that macro-scale behavior, such as consumer demand or regions of the brain or material properties, etc., is the aggregation of a large number of micro-scale instances[14, 15, 16]. The complexity of modeling a macro-scale model composed of thousands or millions of micro-scale instances, however, can be unwieldy, and such models almost never exhibit behavioral complexity commensurate with the descriptive complexity of the model as an aggregate of many micro-scale models. Thus, a systematic technique for developing simplified descriptions of macro-scale models that resemble, in a principled way, the aggregate of micro-scale phenomena can be critical in such applications.

Fourth, understanding the resilience and vulnerability of large-scale critical infrastructures demands techniques for modeling the attack surface of complex cyber-physical-human systems. The attack surface is typically a simplified model of the system that highlights the exposed variables and dynamics strongly affected by or observed from them [17, 18, 19, 20, 21, 22]. Thus, the techniques discussed here offer the possibility of modeling the attack surface of large-scale cyber-physical-human systems and contribute to a science of system security.

I-B Overview

Two important model reduction techniques for linear systems include Balanced Truncation and singular perturbation methods. Each of these approaches focuses on a particular aspect of the system to preserve.

Balanced Truncation (BT) was first proposed in [23] and has been explored for continuous and discrete time [24]. A clear presentation is provided in [25]. The basic idea is that a change of basis is used to order the states equally observable and controllable, and to order them from most to least controllable and observable, and then the least controllable and observable states are truncated; see [9] for a survey of projection-based reduction techniques, including balanced truncation, for parameterized dynamical systems. A complete explanation of BT is presented in Section II-A.

Perturbation theory is a well studied area and has a rich background in linear operator theory and the controls field [26, 27, 28], and it has seen renewed interest in recent years [29, 30, 31]. It is commonly applied in the context of well-separated time scales. In this case the ratio of time-scales identifies an explicit “small” parameter in which a series expansion can be computed. The theory has also been applied to balanced realizations, which we will refer to as Balanced Singular Perturbation Approximation (BSPA) [32, 33, 34, 35]. In [33], an alternative parameterization of the transfer function is used to provide a different unification of BT and BSPA using the “generalized singular perturbational approximation.” This parameterization is quite different than the one proposed here and does not lead to the insight that BT and BSPA occur on the boundary of a particular manifold representing the original model class. A complete explanation of BSPA is presented in Section II-B. The Manifold Boundary Approximation Method, introduced in [1], is an algorithm that, given a parameterized model class, systematically generates a “representative” subset of the class parameterized by fewer parameters than the original set of models. The process is enabled by the fact that a parameterized model can be interpreted as a mapping between a parameter space and prediction space. The set of all possible models generated by varying the parameters over all values typically corresponds to a manifold with the model mapping acting as a coordinate chart. To make this abstract idea more concrete, consider a model class with $N$ parameters that is sampled at $M$ different points (times or frequencies). In this case, the model defines a mapping from $\mathbb{R}^{N}$ to $\mathbb{R}^{M}$ whose image is an $N$ dimensional submanifold of $\mathbb{R}^{M}$. We illustrate this conceptually for several simple models in this manuscript. Each point on the manifold corresponds to the sample of a different set of parameters. For practical problems, explicitly constructing the entire manifold would be impractical, but such a manifold can be explored using computational differential geometry. Systematic study of model manifolds from a variety of classes has revealed that they are typically bounded and that the boundary consists of a hierarchical cell complex, i.e., a hierarchy of faces, edges, corners, etc., similar to a high-dimensional polygon. Each boundary cell is a manifold that corresponds to a model class with fewer parameters. The “best” approximation of this model can be found by locating the closest boundary point with the desired level of complexity, where closest depends on your metric of choice. MBAM will be formalized in Section III.

In this paper we explore the problem of system approximation and identify an intrinsic structure to the problem by proposing an appropriate parameterization for LTI systems that enables us to show that BT and BSPA can be derived by applying MBAM to it. Given this insight, and since MBAM does not depend on linearity, only on a given parameterization, we conjecture that MBAM may provide a framework for the systematic model reduction of more general (nonlinear) systems. In [36], Andersson et al. provide bounds for model reduction on systems that include uncertainties or nonlinearities, modeled as a delta block, assuming the delta block can be modeled using integral quadratic constraints (IQC). The authors also show that BT and BSPA appear as natural limits of the proposed IQC model, but do not provide any insight as to how to interpolate between the two. While the model in [36] is more general than the model considered herein, we believe that the ideas proposed in this work can be applied to a broader classes of models, independent of the IQC assumption.

In Section II we review Balanced Truncation and Balanced Singular Perturbation Approximation. We present MBAM in Section III. In Section IV, we present a parameterization of LTI systems, and then in Section V we use it to derive BT and BSPA using MBAM. In Section VI we present a simple example that illustrates the result and how MBAM can give insight into model reduction of linear systems, allowing a way to interpolate between BT and BSPA. A condensed version of the work is given in [29]. Further contributions of this paper include: 1) the inclusion of the complete proofs; 2) the addition of several illustrative examples; 3) a clearer, more in-depth presentation of MBAM; and 4) an extensive discussion of the parameterization of LTI systems.

II Model Reduction of LTI Systems

In this work we consider LTI systems described by,

[TABLE]

where $x(t)\in\mathbb{R}^{n}$, $y(t)\in\mathbb{R}^{p}$, and $u(t)\in\mathbb{R}^{m}$. We will assume the system described by (1) is minimal and stable, and note that it generates a unique transfer function given by:

[TABLE]

which is a $p\times m$ matrix of rational functions in the complex variable $s\in\mathbb{C}$. Note that while (1) uniquely specifies (2), there are many equivalent realizations of the form (1) for any given input-output map (2).

Model reduction is an important, well–studied problem. Although there are many approaches to model reduction, this work focuses on two types of model reduction for LTI systems: BT and BSPA. We will review these two methods and their accompanying literature.

II-A Balanced Truncation

Consider a stable system given by a bounded map (2), and let some (arbitrary) minimal realization of it, (1), be known. In general, the description (1) may not reveal any particular insight about the system, compared with any other realization of the system, and all such minimal realizations can be found by changing coordinates of the state variable, i.e. by rewriting (1) in terms of a new basis of the state space obtained by multiplying $x$ by any invertible matrix, $\bar{x}=Tx$, and considering the system in terms of $\bar{x}$; we call $T$ a state transformation.

Nevertheless, while an arbitrary realization may be no more useful than any other, some particular realizations are more useful for particular types of analysis than others. For example, modal realizations change basis of the state space to reveal the eigenstructure of each dynamic mode of the system. Any such “particular” realization is called canonical when there is a well-defined process for constructing it, i.e. when such a realization always exists, and when it is unique for every system $G$. Thus, a canonical realization is in one-to-one correspondence with any given input-output map (2).

It turns out that there is a canonical realization of LTI systems that has a property we call balanced. That is to say, although there may be multiple realizations with the property of being balanced, we can specify a fixed procedure that will repeatably construct such a realization for any system $G$. We will represent this realization by:

[TABLE]

where $\bar{A}$ satisfies the Lyapunov equations:

[TABLE]

with $X=diag(\theta_{1},\dots,\theta_{n})$ and the $\theta_{i}$’s are the Hankel singular values (HSVs) of the system [25]. These Lyapunov equations are solved by the controllability and observability Gramians; since they are solved by the same diagonal matrix, $X$, we call the realization balanced, that is, each state in this basis is as observable as it is controllable, quantified by the corresponding $\theta_{i}$. The HSVs are ordered from largest to the smallest in magnitude, and consequently the states are ordered from the most controllable/observable to the least.

Once balanced, one may partition the states into two sets, a highly controllable/observable set and a low controllable/observable set. Partitioning the state matrices commensurate with this partition of the states then yields:

[TABLE]

where $\bar{A}_{11}\in\mathbb{R}^{(n-k)\times(n-k)}$, $\bar{A}_{22}\in\mathbb{R}^{k\times k}$, and the rest of the blocks are the appropriate dimensions. BT of $k$ states projects the system to become:

[TABLE]

Note that stability is preserved under truncation [25]. The transfer function for the truncated system becomes

[TABLE]

Using the H-infinity norm, if $\theta_{n-k+1}<\theta_{n-k}$ and $\theta_{n-k+1}=\cdots=\theta_{n}$,

[TABLE]

If $\theta_{n-k+1}<\theta_{n-k}$ and $\theta_{n-k+1}>\cdots>\theta_{n}$,

[TABLE]

Proofs of these error bound can be found in [37, 38, 39, 25]. BT has been extended to several different classes of systems including time-varying, multidimensional, and uncertain systems [40, 41, 42], and many structure preserving model reduction techniques, with varying definitions, have been considered as extensions of BT [43, 44, 45, 46, 47]. Extensions have also been developed for nonlinear systems [48, 49, 50].

II-B Balanced Singular Perturbation Approximation

Given a partitioned balanced realization characterized as in (5), the reduced system given by BSPA is

[TABLE]

where

[TABLE]

The matrix $\hat{A}$ is the Schur complement of $\bar{A}$, which we will denote by $\bar{A}/\bar{A}_{22}$. Note, similarly to BT, stability is preserved under the BSPA process [35].

Define the transfer function of the BSPA system to be

[TABLE]

Then, identically to BT, if $\theta_{n-k+1}<\theta_{n-k}$ and

$\theta_{n-k+1}=\cdots=\theta_{n}$,

[TABLE]

If $\theta_{n-k+1}<\theta_{n-k}$ and $\theta_{n-k+1}>\cdots>\theta_{n}$,

[TABLE]

These bounds are proven in [35].

In many respects, BT and BSPA are complementary types of approximations. They are both derived from the same block partition of a balanced realization. They share the same error bounds, illustrated in (6)-(7) and (9)-(10). It is well known that BT typically provides better approximations at high frequencies while BSPA works well at low frequencies [35]. In what follows, we will see how these similarities allow us to unify both approximations as limiting approximations of a balanced system.

III Manifold Boundary Approximation Method

The Manifold Boundary Approximation Method, i.e. MBAM, was originally described in [1], but we present an overview here for completeness. The idea is that a parameterized set of models can be viewed as a mapping between its parameter space and a prediction space, sometimes called “data space”, which is the space of sampled, measured outputs. As such, it is natural to interpret a model class, or a continuously parameterized set of models, as a manifold embedded in the space of possible predictions (i.e. embedded in data space). We refer to this manifold as the model manifold, denoted by $\mathcal{M}(\mathrm{D})$.

Control-oriented models, whether linear or nonlinear, are often written as differential algebraic equations, or DAEs:

[TABLE]

where $\mathbf{x}$ is the vector of (differential) state variables, $\mathbf{z}$ are the algebraic variables, $\mathbf{p}$ are parameters, $\mathbf{u}$ are inputs (typically assumed to be known in estimation studies) and t is the (scalar) time variable. For example, [51] illustrates how DAEs with differentiation index zero can be used to represent interconnections of systems. Likewise, in general we assume the system measurement vector is of the form:

[TABLE]

Note that we begin our work with state space models because typically first-principles modeling techniques result in large sets of such systems of equations, often with more parameters than one can accurately estimate from available data, thus motivating the need for a reduction method. These equations lead to the following definition:

Definition 1**.**

A model mapping of the set of systems in (11)-(13) sampled at fixed times is defined as the map $\mathcal{M}~{}:~{}\mathrm{D}~{}\subset~{}{\cal P}\rightarrow{\cal D}$, where ${\cal D}=\mathbb{R}^{M\cdot p}$ is the data or prediction space, with $\bar{N}<M$.

In this work, we assume that $f$, $g$, and $h$ in (11), (12), and (13) are such that $\mathcal{M}$ is smooth. Furthermore, note that in Definition 1, we have assumed that $\bar{N}<M$, with no bound on $M$. It is possible for the data or prediction space to be infinite dimensional. In many cases, however, we sample the output at specific times, making the data space finite dimensional. Note that while other methods, such as [48], sample system data to build an empirical controllability Gramian, the sampling here embeds the model manifold in an appropriately large data-space. Moreover, the topological properties of the manifold are largely invariant to details of the sampling [52]. Thus, the model manifold becomes the central object to represent the parameterized model class describing the system.

We use the notation $\mathcal{M}(\mathrm{D})$ to denote the image of $\mathrm{D}$, or the set of points in $\mathbb{R}^{M\cdot p}$ that are mapped from $\mathrm{D}$; this set of points is a manifold embedded in $\cal D$, called the model manifold. Each point on the manifold, then, corresponds to a particular choice of parameter, $p\in\cal P$, which specifies a given system in the model class. For a given model manifold $\mathcal{M}(\mathrm{D})$, we denote the closure of the manifold as $\overline{\mathcal{M}(\mathrm{D})}$.

Definition 2**.**

Given a model mapping $\mathcal{M}:\mathrm{D}\subset{\cal P}\rightarrow{\cal D}$, a point $d\in\overline{\mathcal{M}(\mathrm{D})}$ is interior if there exists an open neighborhood of $\overline{\mathcal{M}(\mathrm{D})}$ centered at $d$. If $d$ is not interior then it is a boundary point.

The set of boundary points of the manifold closure defines the boundary of the manifold, denoted by $\partial\mathcal{M}(\mathrm{D})$.

Definition 3**.**

Given a model mapping $\mathcal{M}:\mathrm{D}\subset{\cal P}\rightarrow{\cal D}$, a model mapping $\mathcal{\tilde{M}}$ is a $k^{th}$-order manifold boundary approximation of $\mathcal{M}$ if:

$\mathcal{\tilde{M}}:\mathrm{\tilde{D}}\subset\tilde{{\cal P}}=\mathbb{R}^{\bar{N}-k}\rightarrow{\cal D}$* and* 2. 2.

$\mathcal{\tilde{M}}(\mathrm{\tilde{D}})\subset\partial\mathcal{M}(\mathrm{D})$.

That is to say, a $k^{th}$ order model boundary approximation is itself a manifold defined as being on the boundary of the model manifold and having $k$ less parameters than the original model manifold. A model boundary approximation method, then, is a systematic process for choosing a $k^{th}$-order approximation to the model mapping on the boundary of the model manifold.

Note that in the case where the output is sampled, one needs to collect more samples than parameters to ensure that the model manifold is embedded in a large enough space. In general, though, as long as there are enough sample points, the topological features of the resulting model manifold are invariant to the specific times when the samples are taken as long as they are decided randomly–see [52] for details. Once a model manifold is constructed, however, we then turn our attention to the specific details of choosing a $k^{th}$-order approximation.

To develop such an MBAM method, we consider the sensitivity of such models to their parameters by taking partial derivatives as follows:

[TABLE]

Note that these equations are linear in terms of sensitivities, but the matrices involved do vary with time. We see that (14)-(16) are derived by differentiating (11)-(13) with respect to the parameters and applying the chain rule to account for the implicit dependence of the dynamic, algebraic, and measured variables on the parameters. We also use the second order sensitivities in the sequel, the equations for which can be derived similarly. The problem of calculating parametric sensitivities for dynamical systems is well-known, with a long history [53]. However, deriving expressions for the first and second order sensitivities by hand can be tedious and error prone (particularly for large models). We therefore use automatic symbolic and numerical differentiation [54, 55] to simplify the process.

Previous work suggests that model manifolds for real, physical systems tend to be “thin,” in the data space where they are embedded; that is to say, measured variables tend to be considerably less sensitive to some parameter combinations than others [56, 57, 1]. MBAM exploits this “sloppy” character of such model classes by approximating this “thin” manifold by its boundary. We identify the boundary by numerically computing a geodesic on the model manifold and using the results from this calculation to identify an approximate model. This “boundary” model will have fewer parameters than the original, and it approximates the original model manifold the same way a long edge of a ribbon approximates the ribbon.

Geodesics are the analogs of straight lines generalized to curved surfaces, and we compute them numerically as the solution to a second order ordinary differential equation in parameter space (while utilizing quantities from the data space):

[TABLE]

where $\Gamma$ are the Christoffel symbols [58, 59, 57], containing curvature information about the mapping between parameter space and data space. These are expressed in terms of the first and second order parametric sensitivities (14)-(16) and $\mathbf{I}$ is the Fisher Information Matrix (FIM) for the measurement process. The geodesic is parameterized by $\tau$, which is proportional to the arc length of the geodesic as measured on the model manifold, i.e., in data space. Solving (17) gives a parameterized curve $\mathbf{p}(\tau)$ in parameter space that is used to reveal a limiting behavior in the model as we demonstrate below.

Equation (17) is an ordinary differential equation that we solve as an initial value problem. Here, we take the model’s nominal, or “true” parameter values as the starting point of the geodesic. The initial “velocity” is taken to be the least sensitive parameter direction as measured by the FIM. In this way, the MBAM procedure can be summarized as a five step algorithm. Here, we describe the algorithm and give a summary in Fig. 1.

Step a: The least sensitive parameter combination is identified from an eigenvalue decomposition of the FIM which becomes the initial geodesic velocity, $\partial\mathbf{p}/\partial\tau$ The geodesic acceleration $\partial^{2}p/\partial\tau^{2}$ is given by (17). If $\sum_{i}(\partial p^{i}/\partial\tau)(\partial^{2}p^{i}/\partial\tau^{2})<0$, then reverse the initial velocity: $\partial\mathbf{p}/\partial\tau\rightarrow-\partial\mathbf{p}/\partial\tau$. This heuristic resolves the ambiguous direction associated with the eigenvalue.

Step b: A geodesic on the model manifold is constructed by numerically solving (17) using the nominal or ”true” parameter values and the velocity $\partial\mathbf{p}/\partial\tau$ calculated in Step a as initial conditions. The limiting behavior of this curve identifies the boundary of the model manifold.

Step c: Having found the edge of the model manifold, the corresponding lower-dimensional model class is identified as an approximation to the original model class. By inspecting which values of the parameter vector become infinite, we identify the boundary as a limiting approximation in the model class. We evaluate this limit to construct the approximate model class parameterized by fewer parameters.

Step d: We choose parameter values of this approximate model class by fitting the approximate model to the behavior of the original model using, for example, least squares regression.

Step e: The procedure is repeated until the reduced model is unable to approximate quantities of interest within the desired tolerance.

Note that there can be many model boundary approximation methods, depending on which boundary one chooses as an approximation of the model manifold. A common approach is to choose a metric and then select the boundary that is closest to a given point in the data space. In the original presentation of MBAM [1], an information distance was imposed on the data space defined by the FIM, which defines a Riemannian metric on the parameter space, [58, 59]. From this Riemannian metric, computational differential geometry can be used to identify candidate boundary approximations [1]. Nonetheless, MBAM is a topological operation in general and is agnostic to the actual choice of metric [52].

III-A Michaelis Menten Example

The MBAM concept is best illustrated by example. Consider the dynamics of the Michaelis Menten Reaction (MMR),

[TABLE]

where $\rho_{i}$ indicates a parameter in the dynamic equations characterizing the family of models under consideration [60, 61, 62]. The parameters $\rho_{i}$ are only physically relevant for positive values; we therefore restrict our attention to the domain $\mathrm{D}$ as the positive quadrant of the parameter space, that is, where $\rho_{i}\geq 0\ \forall i$. Assuming a fixed initial condition, $x_{0}=1$, and varying $\rho_{1}$ and $\rho_{2}$ among all possible elements of the domain, $\mathrm{D}$, we observe some sample time series for this model as shown in Fig. 2(a).

Consider now three observations of this system at time points, $(t_{1},t_{2},t_{3})$, indicated by the vertical dotted lines in Fig. 2(a). To visualize the sensitivity of the model to variations in the parameters, it is common to construct an “error function” in parameter space as in Fig. 2(b). An alternative approach is to recognize that the space of all possible experimental data at these three time points forms a three-dimensional “data space.” All possible model predictions, for different values of $\rho_{1}$ and $\rho_{2}$, correspond to a two-dimensional subset of this space. This two-dimensional surface is the model manifold characterized as in Definition 1. The manifold for the MMR model illustrated in Fig. 2(c) is bounded, a feature shown to be common among model manifolds with many more parameters than can be estimated from available data [56, 57]. The existence of these boundaries is the crucial element that enables the manifold boundary approximation, where each boundary corresponds to a different model reduction of the original model in (18). The “top” boundary reached by the blue geodesic curve, corresponds to the limit where $\rho_{1},\rho_{2}\rightarrow\infty$, and the ratio $\rho_{1}/\rho_{2}$ becomes the new parameter. The “bottom” boundary reached by the green geodesic curve corresponds to the limit where $\rho_{2}\rightarrow 0$, leaving $\rho_{1}$ as the sole parameter. These two different limiting approximations each lead to two distinct reduced models of the form:

[TABLE]

respectively. A well–known interpretation of the MMR in biochemistry is that it describes a saturable reaction. That is, at low concentrations (small values of $x$), the reaction rate is nearly first order (i.e., linear in $x$), while at high concentrations the rate saturates into a zeroth order (i.e., constant with respect to $x$). Notice how these two limits are naturally identified as the boundaries of the model manifold in Fig. 2(c) and together bound the space of predictions that can be reached by the model.

To be clear, the use of MBAM to approximate the MMR system here reduced the set of systems from a two dimensional set to a one dimensional set, but it did not reduce the dynamical order of the system, which started as a first order system and remained, in each limiting case, a first order system. This is because parameter reduction does not necessarily lead to model-order reduction. The next example will demonstrate a case where parameter reduction does reduce the order of the model. Although the dynamical order was not reduced in the example above, the reduced models are simpler. While the original model exhibited nonlinear dynamics, the approximate models in (19) and (20) are linear and constant, respectively. Note that code for computing each step of the MBAM process can be found at [63].

III-B Second-Order LTI Example

As a second example, consider the set of systems described by:

[TABLE]

where $t,\rho_{1},\rho_{2}\in\mathbb{R}^{+}$. For any choice of $\rho_{1}$ and $\rho_{2}$, the behavior of the corresponding system is the set of functions

[TABLE]

When we fix the initial condition, say to x(0)=\left[\begin{array}[]{cc}1&1\end{array}\right]^{T}, this set is completely parameterized by the parameters $\rho_{1}$ and $\rho_{2}$. One way to characterize the behavior of this set of systems is to associate each choice of parameters with an observation vector given by v=\left[\begin{array}[]{ccc}y(t_{1})&y(t_{2})&y(t_{3})\end{array}\right]^{T}, where the observation times $t_{1}$, $t_{2}$, and $t_{3}$ are fixed. Spanning over all choices of parameters $\rho_{1}$ and $\rho_{2}$, the observation vectors $v(\rho_{1},\rho_{2})$ then sweep out a two-dimensional manifold embedded in $\mathbb{R}^{3}$, shown in Fig. 3(b)-3(c). This is the model manifold for this set of systems, and we note that although its size and shape somewhat depend on the observation times $t_{1}$, $t_{2}$, and $t_{3}$, its characteristic topological features (number of cusps and edges and the relationships among them) are determined solely by the relationship between the parameters, $\rho_{1}$ and $\rho_{2}$, and observations, $y(t)$. That is to say, regardless of how complicated the underlying dynamics are, an $\bar{N}$-parameter system generates an $\bar{N}$-dimensional manifold embedded in a larger data space (a fact that therefore requires the number of observation points to be larger than the number of parameters), and the relationship between these cusps and edges is fixed, regardless of when we make the observations or whether we make more than $n$ observations, etc. [52].

This invariance of the key topological features of the model manifold allows us to consider a graphical representation among these features, shown in Fig. 4.

This graph encodes the relationship between the model manifold surface in Fig. 3 and its boundaries, and between its boundaries and their endpoints. Each of these objects are themselves a manifold of lower dimension, and this hierarchy of lower dimensional manifolds forms a partial ordering that reveals different ways the fully parameterized model can be approximated by representations of lower complexity.

Thus, for example, each point on the model manifold, $M$, in Fig. 3 is specified by choosing particular values for both of the parameters $\rho_{1}$ and $\rho_{2}$; the fact that the model manifold is two-dimensional results from the fact that two parameters must be specified. When one of these parameters is set to one of its limiting values, zero or infinity, the resulting set of models parameterized by the other parameter form a boundary on the model manifold, a boundary of dimension one (since only one parameter is left free). Thus one edge is formed when one of the parameters (without loss of generality, assume it is $\rho_{1}$) is fixed at zero, yielding the remaining behaviors $y(t)=1+e^{-\rho_{2}t}$ represented by the one-dimensional manifold, the edge $e_{AB}$. Another edge is formed when $\rho_{2}$ approaches infinity, yielding $y(t)=e^{-\rho_{1}t}$ and represented by the edge $e_{BC}$. Finally, the last boundary of $M$ is formed when $\rho_{1}$ is fixed to the same value as $\rho_{2}$, yielding $y(t)=2e^{-\rho_{2}t}$, and is represented by the edge $e_{AC}$.

The three points where these edges intersect form zero-dimensional manifolds and represent systems where both parameters are fixed at one of the limiting values. At point $A$, $\rho_{1}=\rho_{2}=0$ and the output $y(t)\equiv 2$. At point $B$, one of the parameters ($\rho_{1}$ in our case above) is zero while the other approaches infinity, with $y(t)\equiv 1$. At point $C$ both parameters approach infinity and $y(t)\equiv 0$.

The partial ordering among these models of decreasing complexity is captured by the Hasse diagram, [52], shown in Fig. 4; this diagram captures the relationship between model complexity and different types of approximations. So, for example, $e_{AB}$, with behavior $y(t)=1+e^{-\rho_{2}t}$, is an approximation that considers that the unmodeled mode of the system is much slower than the other; it models the fast dynamics of the system. On the other hand, $e_{BC}$, with behavior $y(t)=e^{-\rho_{1}t}$, is an approximation that considers that the unmodeled mode of the system is much faster than the other and has already reached equilibrium; it models the slow dynamics of the system. Finally, $e_{AC}$ describes the situation where both decay rates are comparable; it models the dynamics of a symmetric system.

Each of these three types of approximations capture a different idealized behavior of every model in the set of systems parameterized by $\rho_{1}$ and $\rho_{2}$. Likewise, each of these approximations are themselves approximated at their endpoints by one of three extreme behaviors $A$, $B$, or $C$. These behaviors are themselves stationary, describing the regimes of behavior characterizing these systems. In this example, $A$ represents a “high” value, $B$ a “medium” value, and $C$ a “low” value, but in other contexts they might represent different phenotypes or behavior regimes such as “healthy” and “sick,” etc.

These examples all illustrate the underlying mechanism driving MBAM techniques–that of looking for approximations on the boundary of the model manifold. Current research is exploring efficient algorithms for conducting such a search, as well as application of the ideas to large scale systems, [52], but this work reports how both BT and BSPA, and a host of hybrid techniques, become special cases of MBAM methods when applied to LTI systems.

IV Parameterizations of LTI Systems

The previous examples show how MBAM can reduce the number of parameters used to describe a set of models. The reduced parameter descriptions for the first example, given by (19) and (20), parameterize a subset of the original model class, but this subset is chosen to be “representative” of the original class in the sense that it is a set of models on the boundary of the original set. We can parameterize this set of boundary models with fewer parameters than the original set, hence MBAM is a parameter-reduction technique.

Parameter reduction is not necessarily model order reduction, however. The example in Section III-A shows that MBAM can reduce the number of parameters without changing the order of the system, while the example in Section III-B shows that depending on the parameterization the order can be reduced. The interest of this work is to use MBAM to perform model order reduction on general LTI systems. This requires a careful examination of the parameters of such systems, so that reduced-order models generated by BT or BSPA can be shown to be instances of reduced-order models generated by MBAM. In this section we present and classify the parameters of general LTI systems, and the following sections will leverage this parameterization to prove this relationship between these two model order reduction techniques and provide a framework to interpolate between them.

IV-A Parameterizations of Transfer Functions

Consider the LTI system in (1), which has $m$ inputs and $p$ outputs. The input-output behavior of such a system is characterized by

[TABLE]

where $U(s)$ and $Y(s)$ are the Laplace transforms of the input and output signals, respectively, and $G(s)$ is given by

[TABLE]

In this context, the operator $G$ is the transfer function of the LTI system, and shows what the output of the system will be given any input. Thus $G$ is a $p\times m$ matrix, where each entry is a proper polynomial of the Laplace variable $s\in\mathbb{C}$, which we will write as

[TABLE]

Without loss of generality, we may consider that each entry has a common denominator $\delta(s)$, with degree $N$, and the degree of every polynomial $\nu_{ij}(s)$ is equal to $N$. The matrix $G$ can be rewritten as

[TABLE]

where $\nu(s)=[\nu_{ij}(s)]$. We assume that $\delta(s)$ is monic (the leading coefficient is one), since if it were not, in each entry, both $\delta(s)$ and $\nu_{ij}(s)$ could be divided by the leading coefficient in $\delta(s)$ without changing the behavior of $G(s)$. Thus, each $\nu_{ij}(s)$ is determined by the $N+1$ coefficients; therefore $\nu(s)$ is completely characterized by $pm(N+1)$ coefficients. Adding in the $N$ coefficients of $\delta(s)$ (not $N+1$ because $\delta(s)$ is monic), it follows that $G(s)$ is completely characterized by $Npm+N+pm$ coefficients (see Fig. 5(a)).

These coefficients provide a way to minimally parameterize the set of all $p\times m$ LTI operators of degree $N$. The parameterization is minimal in that if any coefficient is perturbed, then this perturbation will uniquely change the behavior of the LTI system, and no coordinated changes among the other parameters can produce the same change in the input-output behavior. Moreover, all LTI systems of size $p\times m$ with degree $N$ have such a parameterization. In this context we conclude that a $p\times m$ transfer function $G(s)$ with degree $N$ is uniquely characterized by $Npm+N+pm$ parameters.

Another way to minimally parameterize the set of $p\times m$ degree $N$ transfer functions is to employ a partial fraction expansion on each element $\frac{\nu_{ij}}{\delta(s)}$. For systems where every pole has a geometric multiplicity of one, the transfer function then becomes

[TABLE]

where $\lambda_{i}\in\mathbb{C}$ are the roots of $\delta(s)$, and $D,K_{i}\in\mathbb{R}^{p\times m}$. The matrix $D$ is a matrix of the leading coefficients for each $\nu_{ij}$ and is equal to $D$ in (1) for any state space realization of $G(s)$; indeed, every state realization of $G(s)$ will have the same $D$ matrix. The $K_{i}$ matrices contain the scalar numerators corresponding to $s-\lambda_{i}$ in the partial fraction expansion. In this representation (called Gilbert’s realization [64]), the parameters are $D$ ($pm$ entries), $\lambda_{i}$’s ($N$ elements of the spectrum), and $K_{i}$ ($pm$ entries for each $i=1,...,N$). Thus we see that the number of parameters remains the same as the previous parameterization: $Npm+N+pm$ (see Fig. 5(b)). Nevertheless, this particular parameterization is often more descriptive of the actual system dynamics: not only are the input-output dynamics clearly defined, but as we will see in Section IV-B, each $\lambda_{i}$ is a mode of the system and the $D$ matrix represents how the inputs directly affect the outputs.

Example 1**.**

Consider the representation in (24) where $K_{i}$ are all rank one matrices; this scenario arises when the poles of the system have algebraic multiplicity of one (i.e., they are distinct). For a particular $i$, instead of $K_{i}$ being parameterized by $pm$ parameters, here many of the parameters of $K_{i}$ are dependent on the other parameters, in order to ensure that the rank one constraint is met. Let $c_{i}\in\mathbb{R}^{p}$ and $b_{i}\in\mathbb{R}^{m}$ for all $i$, where $c_{i}b_{i}^{T}=K_{i}$. We can assume that the first entry of $c_{i}$ is one, if not, we could divide $c_{i}$ by its first entry and multiply $b_{i}$ by the same, thus leaving $K_{i}$ unchanged. Under this circumstance, each $K_{i}$ now only requires $p+m-1$ entries to be specified in order to completely characterize it. Additionally, now any perturbations in these entries will uniquely change $K_{i}$. Therefore, these $p+m-1$ entries are the parameters of $K_{i}$. Thus for this example, the total number of parameters for $G(s)$ is $Np+Nm+pm$ (see Fig. 5(c)).

We note that both parametrizations of the set of all $p\times m$ transfer functions of degree $N$, (23) and (24), have the same number of parameters. Generically, almost any element of this set is uniquely identified by a choice of these parameters, and sufficiently rich input-output data enables the identification of the parameters specifying the system generating the data. As a result, in either case, we call these $Npm+N+pm$ parameters identifiable. This definition of parameter identifiability is consistent with the definition of structurally identifiable [65].

IV-B Parameterizations of State Realizations

In the previous section we saw that the set of all transfer functions of a given size and degree was efficiently parameterized by a particular number of identifiable parameters, and we demonstrated two distinct ways of doing this. In this section we extend these ideas to state realizations.

Recall that the standard form for state realizations is given in $\eqref{eq:ss}$, and that we have assumed minimality. The set of all state realizations with $m$ inputs and $p$ outputs of order $n$ is parameterized by four matrices: $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}$. Counting each entry in each matrix as a parameter yields a total of

[TABLE]

parameters for the set.

In the generic case, $A$ is diagonalizable with no repeated eigenvalues. This implies that there exists a transformation $T$ such that $x=T\tilde{x}$ and $T^{-1}AT=\tilde{A}$ is diagonal. Likewise, $T^{-1}B=\tilde{B}$, $CT=\tilde{C}$, and $D=\tilde{D}$. In this case, by (21),

[TABLE]

Each term $\tilde{c}_{i}\tilde{b}_{i}^{T}$ is a rank one matrix (being the outer product of two vectors), thus it matches the transfer function parameterization given in (24), satisfying the constraint given in Example 1 that each $K_{i}$ is rank one. Example 1 shows that there are $Nm+Np+pm$ parameters in the set of such transfer functions, where $N$ is the degree of $\delta(s)$. Since $A$ has no repeated eigenvalues, it follows from (26) that $N=n$. This implies that the set of corresponding transfer functions has $nm+np+pm$ parameters. Therefore, by (25), the set of corresponding state realizations has $n^{2}$ more parameters.

This discrepancy in parameters is not surprising, since it is well-known that any transfer function will have infinitely many corresponding state realizations. Further, it is also well-known that any two state realizations of a transfer function are related by an $n\times n$ state transformation matrix $T$ as described prior to (26). Therefore, once a $G$ has been specified with $nm+np+pm$ parameters, $n^{2}$ additional parameters are required to specify an instance of $(A,B,C,D)$. We say that these $n^{2}$ parameters are structural, since they do not affect the system dynamics and are not identifiable from input-output data, but they do affect the internal structure of the system.

The matrices given in (1) are one way to parameterize the set of state realizations. In this parameterization, it is clear from (24) that $D$ can be identified directly from knowing $G$, so its entries are also identifiable parameters. That is to say, every state space realization of $G$ has the same $D$ matrix. However, the rest of the parameters cannot be identified simply by knowing $G$: one must also know $n^{2}$ additional parameters (fixing the structure, in a sense, by fixing the basis of the state space). For instance, if we fix $G$, there are many admissible choices for $B$ and $C$, but it could be the case that if $B$ were known, $C$ could be determined from data. It also is the case that the reverse is true: fixing $C$ would allow $B$ to be determined from data. Therefore, we refer to the entries of $B$ and $C$ as conditionally identifiable: once a certain number of them have been fixed, the rest can be determined from data (see Fig. 7(a)). An illustrative example is also shown in Fig. 6.

Conditionally identifiable parameters are not ideal when performing MBAM, because it is not clear whether or how the structure or dynamics will be affected. Therefore, it is important to find a state space parameterization which has only identifiable and structural parameters.

IV-C Parameterization Using the Balanced Realization

When performing MBAM on a state realization, it is important to know which parameters will affect the dynamics of the system (identifiable) and which will affect the implementation of the dynamics (structural), since either may be fixed, depending on the application. Parameterizations that partition the parameters into identifiable and structural without the need for conditionally identifiable parameters are useful when applying MBAM so that unintended parameters are not eliminated. Unique canonical realizations allow one to accomplish this goal. All parameters from a unique canonical realization can be identified from data, since the canonical realization fixes the structure. Then, any realization can be parameterized using the parameters of the canonical realization as the identifiable parameters and the transformation matrix $T$ between the two realizations as the structural parameters.

We illustrate this point by employing the balanced realization. Consider a minimal, stable system as in (1). It is well known that there exists a state transformation from these matrices to an input–output equivalent balanced realization $(\bar{A},\bar{B},\bar{C},\bar{D})$ as in (3), where $X$ is the diagonal matrix of HSVs ($\theta_{i}$’s) satisfying (4). Notice the following simple statement is true for balanced realizations.

Lemma 1**.**

Given a balanced realization $(\bar{A},\bar{B},\bar{C},\bar{D})$, $diag(\bar{B}\bar{B}^{T})=diag(\bar{C}^{T}\bar{C})$.

Proof.

Since $(\bar{A},\bar{B},\bar{C},\bar{D})$ is balanced, the observability and controllability Gramians are equal with the HSVs, $\theta_{1},\dots,\theta_{n}$, on the diagonal. The diagonals of the Lyapunov equations give

[TABLE]

where the subscript $ii$ indicates the $i$th diagonal entry of the matrix and $a_{ii}$ is the $i$th diagonal entry of the $\bar{A}$ matrix. Therefore $diag(\bar{B}\bar{B}^{T})=diag(\bar{C}^{T}\bar{C})$. ∎

Note that this lemma is also a result of Theorem 1 in [67].

Let the common diagonal entries of $\bar{B}$ and $\bar{C}$ be denoted by $r_{1}^{2},\dots r_{n}^{2}$. We can then write the $\bar{B}$ matrix as

[TABLE]

where the $\beta_{i}$’s are a collection of normalized column vectors in $\mathbb{R}^{m}$ satisfying $\beta_{i}^{T}\beta_{i}=1$, for all $i=1,\dots,n$. Denoting the $n\times m$ matrix whose rows correspond to $\beta_{i}^{T}$ as $\beta^{T}$ and introducing $R=diag(r_{i},\dots,r_{n})$, it follows that $\bar{B}=R\beta^{T}$. Clearly by construction $diag(\bar{B}\bar{B}^{T})=(r_{1}^{2},\dots,r_{n}^{2})$.

Similarly, we can write $\bar{C}$ as

[TABLE]

where the $\gamma_{i}$’s are a collection of normalized column vectors in $\mathbb{R}^{p}$ satisfying $\gamma_{i}^{T}\gamma_{i}=1$, for all $i=1,\dots,n$. We write $\bar{C}=\gamma R$; note that $diag(\bar{C}^{T}\bar{C})=(r_{1}^{2},\dots,r_{n}^{2})=diag(\bar{B}\bar{B}^{T})$, consistent with Lemma 1.

Plugging $r_{i}^{2}$ into (27) gives

[TABLE]

for the diagonal elements of the $\bar{A}$ matrix. From the off-diagonals of the Lyapunov equations we find,

[TABLE]

for $i\neq j$, where

[TABLE]

Leveraging these properties, we see that a balanced realization is specified by $(\theta,\beta,\gamma,R,D)$ as follows:

[TABLE]

where $\bar{A}$ is defined in (28) and (29). We know that $\theta$ and $R$ both contain $n$ parameters and that $D$ contains $pm$ parameters. While $\beta$ is an $n\times m$ matrix, since $\beta^{T}_{i}\beta_{i}=1$, specifying $n-1$ entries in $\beta_{i}$ fixes the magnitude of the final entry. Thus $\beta$ carries only $n(m-1)$ parameters. Likewise, $\gamma$ carries $n(p-1)$. It follows that any realization can be specified by these parameters along with an additional transformation matrix $T$ that dictates the change of basis from the balanced realization.

Using this balanced parameterization, we again have $n^{2}+nm+np+pm$ parameters (see Fig. 7(b)), but now each parameter is clearly labeled as either identifiable or structural. Since any realization can be uniquely described by $T$, its transformation matrix to the balanced realization, we know that there are $n^{2}$ structural parameters. As discussed, the remaining $nm+np+pm$ parameters match the same number as $G(s)$ in Example 1, and therefore are all identifiable from data. As a result this parameterization enables one, when using MBAM to reduce the number of parameters, to be deliberate as to how the structure and/or dynamics are affected.

V Unification of Model Reduction Techniques

Given the parameterization of LTI systems of the previous section, we are equipped to apply the algorithm summarized in Figure 1. In our experience solving the geodesic equation on multiple LTI systems, we empirically observed that the boundary approximations always take one of two forms when reducing the order of the system. This empirical observation suggests that it may not be necessary to sample the model predictions and solve Eq. 17. In other words, it is not necessary to explicitly construct the data space to find the mathematical form of the boundary approximation, all of this work can be done directly in parameter space.

In this section, we make this empirical observation rigorous through a sequence of theorems demonstrating that balanced truncation and balanced singular perturbation are each manifold boundary approximations of LTI systems. Note that in doing this, we have emphasized the conceptual distinction between a boundary approximation (see Definition 3) and the algorithm for finding them that requires solving a geodesic (see Figure 1).

V-A Balanced Truncation from MBAM

For the first theorem of this section we will restrict ourselves to considering the HSVs as the parameters, holding $R$, $\beta$, $\gamma$ and $D$ fixed.

Theorem 1**.**

Consider a balanced realization, as in (30), where the Hankel singular values ($\theta_{i}$’s) are parameters and the matrices $R$, $\beta$, $\gamma$ and $D$ are fixed. Balanced Truncation of $k$ states is equivalent to $k$ iterations of MBAM taking the relevant $\theta_{i}$’s $\rightarrow 0$.

Proof.

Consider the equation for $\dot{x}_{n}(t)$ from (30):

[TABLE]

Multiplying through by $\theta_{n}$ gives

[TABLE]

Performing an MBAM approximation by taking the limit $\theta_{n}\rightarrow 0$ gives that $\bar{x}_{n}=0$. Plugging this back into the dynamics of the rest of the system, i.e. $\dot{\bar{x}}_{i}(t),\ i<n$, gives

[TABLE]

which is equivalent to BT of one state. Clearly also having $x_{n}=0$ gives that

[TABLE]

Iterating this MBAM approximation $k-1$ more times, always choosing the smallest HSV, completes the proof. ∎

Theorem 2**.**

Consider a balanced realization, as in (30), where $r_{1},\dots,r_{n}$ are parameters and the Hankel singular values, $\beta$, $\gamma$, and $D$ are fixed. Balanced Truncation of $k$ states is equivalent to $k$ iterations of MBAM taking the relevant $r_{i}$’s $\rightarrow 0$.

Proof.

Consider the equation for $\dot{\bar{x}}_{n}(t)$ in (31). Performing an MBAM approximation by taking the limit $r_{n}\rightarrow 0$ gives that $\dot{\bar{x}}_{n}=0$. Plugging $r_{n}=0$ into the dynamics of the rest of the system, i.e. $\dot{\bar{x}}_{i}(t),\ i<n$, gives

[TABLE]

Also, when $r_{n}=0$,

[TABLE]

This gives BT of one state. Iterating this limit $k$ times, always choosing the $r_{i}$ with the largest subscript, completes the proof. ∎

Note that Theorems 1 and 2 indicate two different paths along the model manifold that both converge to the same point on the boundary. These paths are found by taking different limits of different parameters, but the common limit point is the Balanced Truncation approximation for the system.

V-B Singular Perturbation Approximation from MBAM

In a similar manner BSPA can be derived by applying MBAM to the balanced parameterization.

Theorem 3**.**

Consider a balanced realization, as in (30), where $r_{1},\dots,r_{n}$ are parameters and the Hankel singular values, $\beta$, $\gamma$, and $D$ are fixed. Balanced Singular Perturbation Approximation of $k$ states is equivalent to $k$ iterations of MBAM taking the relevant $r_{i}$’s $\rightarrow\infty$.

Proof.

We will prove this by induction, starting with the one state case. Dividing (31) by $r_{n}^{2}$ gives

[TABLE]

Letting $r_{n}\rightarrow\infty$ then yields

[TABLE]

driving

[TABLE]

Dividing (31) by $r_{n}$ gives

[TABLE]

Taking the limit as $r_{n}\rightarrow\infty$ with $r_{n}\bar{x}_{n}$ remaining finite, since, by (32), $\bar{x}_{n}\rightarrow 0$, gives

[TABLE]

Furthermore, in the remaining equations for $\dot{\bar{x}}_{i}(t)$, $i~{}=~{}1,\dots,n-1$, and $y_{i}(t)$, $i=1,\dots,p$, we find that $r_{n}$ and $\bar{x}_{n}$ always appear in the combination $r_{n}\bar{x}_{n}$. Therefore, this limit is a well-defined boundary approximation for this parameterization. Plugging (33) into the rest of the system, i.e. $\dot{\bar{x}}_{i}(t),\ i<n$, gives

[TABLE]

which is the system in (8) for $k=1$, that is, BSPA of one state. By similarly plugging (33) into $\bar{C}$ and $\bar{D}$, $\hat{C}$ and $\hat{D}$ from (8) for $k=1$ appear.

Now for the inductive step, assume BSPA of $k$ states is equivalent to $k$ MBAM approximations, giving the system in (8), denoted by $(\hat{A},\hat{B},\hat{C},\hat{D})$. Let the system matrices that result from BSPA of $k+1$ states be denoted by $(\hat{A}^{k+1},\hat{B}^{k+1},\hat{C}^{k+1},\hat{D}^{k+1})$. Partition $\bar{A}$ such that

[TABLE]

where $\bar{A}_{k+1,k+1}\in\mathbb{R}^{k+1\times k+1}$ and $\bar{A}_{k,k}\in\mathbb{R}^{k\times k}$; therefore, using the Schur complement notation, $\hat{A}^{k+1}=\bar{A}/\bar{A}_{k+1,k+1}$ and $\hat{A}=\bar{A}/\bar{A}_{k,k}$. By the Quotient Formula for Schur complements (Theorem 1.4 in [68]),

[TABLE]

where

[TABLE]

with $\bar{A}_{k,k}\in\mathbb{R}^{k\times k}$ and $\bar{a}_{k+1,k+1}\in\mathbb{R}$. This means that $\hat{A}^{k+1}$ resulting from BSPA of $k+1$ states (the left hand side of (35)) is equivalent to the analogous part in (34) (the right hand side of (35)).

Now we will show the same holds for $\hat{B}^{k+1}$, which is slightly more laborious since we cannot appeal to Schur complement properties. Partitioning $\bar{A}_{1,k}$ as

[TABLE]

with $\bar{a}_{k+1,k}$ defined in (36) allows us to partition $\bar{A}_{1,k+1}$ as

[TABLE]

where $\bar{a}_{k+1}\in\mathbb{R}^{n-(k+1)}$. Using the above, (8), and (36), and partitioning $\bar{B}_{k+1}$ as

[TABLE]

with $\bar{B}_{k}\in\mathbb{R}^{n-k\times m}$ and $\bar{b}_{k+1}\in\mathbb{R}^{m}$, we have

[TABLE]

where $\bar{B}_{1}\in\mathbb{R}^{n-(k+1)\times m}$. By (36), the block matrix inversion formula, matrix multiplication, and some rearranging of terms,

[TABLE]

where

[TABLE]

Partition $\hat{A}$ and $\hat{B}$ as

[TABLE]

with $\hat{a}_{n-k,n-k}\in\mathbb{R}$ and $\hat{b}_{n-k}\in\mathbb{R}^{m}$. Note that $c$ is equal to $\hat{a}_{n-k,n-k}$, $(\bar{a}_{k+1}-\tilde{A}_{1,k}\bar{A}_{k,k}^{-1}\bar{a}_{k,k+1})=\hat{a}_{1,n-k}$, and $d~{}=~{}\hat{b}^{T}_{n-k}$. Also note that $\bar{B}_{1}-\bar{A}_{1,k}\bar{A}_{k,k}^{-1}\bar{B}_{k}=\hat{B}_{1}$. Therefore (37) becomes

[TABLE]

Similar expressions for $\hat{C}^{k+1}$ and $\hat{D}^{k+1}$ can be found in an analogous way. Therefore (35), (38), and the analogous expressions for $\hat{C}^{k+1}$ and $\hat{D}^{k+1}$ can be obtained, similar to (35), by performing an iteration of MBAM on $(\hat{A},\hat{B},\hat{C},\hat{D})$ by taking $r_{n-k}\rightarrow~{}\infty$, thus completing the inductive step. Therefore, by induction, it holds for all $k$. ∎

The resulting BSPA system is Hurwitz and is still a balanced realization with the $n-1$ HSVs, the same as the largest $n-1$ HSVs of the original system [32].

Note: In (34) if $\theta_{n}$, which is clearly no longer a HSV, is equal to its original value, then the reduced system is BSPA. However if $\theta_{n}$ is set to zero, then the resulting reduction is BT. Therefore, Theorem 3 offers a whole new class of reduced systems ranging between BSPA and BT, and given a metric or goal, the optimal reduction can be found.

VI Illustrative Example

Consider the balanced parameterization for state realizations, presented in Section IV-C, with $x\in\mathbb{R}^{2}$,

[TABLE]

where $\theta_{1}=r_{1}=1$ and $\theta_{2}$ and $r_{2}$ are the free parameters (while enforcing $\theta_{2}\in[0,\theta_{2}]$), and we set $d=0$ and $\beta_{1}~{}=~{}\beta_{2}=\gamma_{1}=\gamma_{2}=1$. We can calculate the frequency responses of the different systems by varying the free parameters $\theta_{2}$ and $r_{2}$. See Fig. 8(a) for some sample frequency responses. We then use the frequencies $(s_{1},s_{2},s_{3})$ to create the manifold in the data space. The model manifold, in Figures 8(b)-8(d), has two boundaries. The cyan boundary, at the top of the manifold, corresponds to the parameterizations where $r_{2}\rightarrow\infty$ (we use one hundred since it is two orders of magnitude larger than one). The red boundary corresponds to parameterizations where $\theta_{1}=\theta_{2}$. The point at which the two boundaries meet at the bottom of the manifold, the magenta dot, is BT, corresponding to the parameterization where $\theta_{2}$ and/or $r_{2}\rightarrow 0$.

Consider the system given by the parameter values $(\theta_{2},r_{2})=(0.7,8)$, depicted on the manifold by the black triangle. It is clear that BT (the magenta dot) is not going to be a good approximation of the system by the distance between the two systems on the manifold. However, BSPA denoted by the green “plus” symbol is a very close approximation. Consider another system given by the parameter values $(\theta_{2},r_{2})=(0.01,0.8)$, depicted on the manifold by the black circle. The reduced model given by BSPA, shown as a green “x”, is still fairly close to the true system; however, it is clear that BT (the magenta dot) is much closer.

Although, as was explained in Section II, BT and BSPA have the same a priori $H_{\infty}$ error bound, BT typically gives better results at high frequencies while BSPA excels at low frequencies [35]. The optimal approximation requires one to identify a metric customized to a context of interest. As a concrete example, returning to the system given by the parameter values $(\theta_{2},r_{2})=(0.01,0.8)$ (the green circle), the cyan curve connecting the black “x” to the magenta dot represents a family of candidate reduced models that interpolates between BT and BSPA. The optimal reduced model relative to the three frequencies sampled, using the two-norm as the metric as measured in the data-space, is the point on this cyan curve that is closest to the green dot. It is identified numerically as $(\theta_{2},r_{2})=(5.3958\times 10^{-4},0.8)$, depicted by the purple diamond. Clearly this point is much closer to the original system than either BSPA or BT and therefore a better approximation for the chosen metric of the two-norm.

Recall that the model manifold, in Figures 8(b)-8(d), is constructed from the sampled frequencies $(s_{1},s_{2},s_{3})$. Different samplings will give different manifolds, and therefore could provide better reduced models for different samplings. If a system designer had a certain frequency of importance $s^{*}$, a manifold could be constructed using sampled frequencies $(\hat{s}_{1},\hat{s}_{2},\hat{s}_{3})$ around $s^{*}$, and find the best approximation of the system by finding the closest boundary point to the original system.

VII Conclusion

This paper demonstrates how both Balanced Truncation and Singular Perturbation Approximations can be viewed as a type of manifold boundary approximation. The key idea unifying these different techniques is to choose a canonical parameterization for any given system that partitions its parameters into identifiable and structural sets, eliminating conditionally identifiable parameters. Application of the Manifold Boundary Approximation Method, an information geometry technique that requires neither linearity nor time invariance, then illustrates that each technique finds approximations on different boundaries of the system’s model manifold embedded in an appropriate data space. Depending on the choice of metric, this manifold can give insight into which approximation is the best for a given application and give alternative approximations interpolating between BT and BSPA.

As we saw in Section III, nothing about the MBAM approximation requires a linear model class, indeed, the Michaelis-Menten Reaction example was nonlinear. This fact, combined with the promising results for linear systems, suggests that by building an analogous parameterization for certain classes of nonlinear systems, similar types of approximations may be obtainable. This framework, then, focuses attention on obtaining good parameterizations of nonlinear systems in order to recover approximations with desirable qualities similar to those of BT and BSPA, questions that can be explored in future work. The work in [48, 49, 50] explores extending balanced truncation to nonlinear systems, providing one possible reduced model. Alternatively, after the cost of constructing a parameterization, the MBAM approach, as illustrated in Section VI, would provide a whole class of possible reduced models, extending BT and BSPA to nonlinear systems and providing a spectrum of reduced-order systems interpolating between the two.