remarkRemark Stability Domains for Quadratic-Bilinear Reduced-Order ModelsB.
Kramer

# Stability Domains for Quadratic-Bilinear Reduced-Order Models††thanks:
Submitted to the editors .

Boris Kramer Department of Mechanical and Aerospace Engineering, University of
California San Diego, CA (, kramer.ucsd.edu<EMAIL_ADDRESS>
###### Abstract

We propose a computational approach to estimate the stability domain of
quadratic-bilinear reduced-order models (ROMs), which are low-dimensional
approximations of large-scale dynamical systems. For nonlinear ROMs, it is not
only important to show that the origin is locally asymptotically stable, but
also to quantify if the operative range of the ROM is included in the region
of convergence. While accuracy and structure preservation remain the main
focus of development for nonlinear ROMs, computational methods that go beyond
the existing highly conservative analytical results have been lacking thus
far. In this work, for a given quadratic Lyapunov function, we first derive an
analytical estimate of the stability domain, which is rather conservative but
can be evaluated efficiently. With the goal to enlarge this estimate, we
provide an optimal ellipsoidal estimate of the stability domain by solving a
convex optimization problem. This provides us with valuable information about
stability properties of the ROM, an important aspect of predictive simulation.
We do not assume a specific ROM method, so a particular appeal is that the
approach is applicable to quadratic-bilinear models obtained via data-driven
approaches, where ROM stability properties cannot—per definition—be derived
from the full-order model. Numerical results for a LQG-balanced ROM of
Burgers’ equation, a proper orthogonal decomposition ROM of FitzHugh-Nagumo,
and a non-intrusive ROM of Burgers’ equation demonstrate the scalability and
quantitative advantages of the proposed approach. The optimization-based
estimates of the stability domain are found to be up to four orders of
magnitude less conservative than analytical estimates.

###### keywords:

Stability domain; domain of attraction; reduced-order models; quadratic-
bilinear dynamical systems

34D20, 34D35, 37E99, 37M22.

##  1 Introduction

Reduced-order modeling provides a mathematical framework for efficient
simulation of complex systems, where large-scale nonlinear dynamical systems
are approximated on low-dimensional manifolds. High-dimensional systems,
arise, e.g., from semi-discretization of partial differential equation (PDE),
and their dimension (degrees of freedom) can be in the order of thousands and
hundreds of thousands. Reduced-order models (ROMs) of those systems are
important in the context of prediction, control, design, and optimization,
see, e.g., [3, 5, 30]. While the development and analysis of ROMs for linear
systems has matured in recent years, ROMs for nonlinear systems—as
expected—face significantly different challenges. For general nonlinear
systems, proper orthogonal decomposition (POD) [18] and the reduced basis
method [17] are the most commonly used model reduction methods. If we restrict
ourselves to the class of quadratic-bilinear (QB) systems, ROM methods have
been developed in the intrusive setting (where governing equations are
available) [4, 6, 10] and in the non-intrusive setting [14, 32, 29, 7], where
the model has to be learned from data. Notably, QB models are a less
restrictive class than the name suggests: many nonlinear dynamical systems can
be transformed into QB form via variable transformations and the introduction
of auxiliary variables, see [24, 20, 31, 15, 4, 21, 32, 16, 29]. Consequently,
reduced-order modeling and analysis developed for QB systems apply to a large
class of nonlinear systems.

The characterization of stability domains is important for open and closed
loop simulation. In the nonlinear control and dynamical systems community,
several approaches exist for estimating the stability domain of (very low-
dimensional) polynomial systems, and in particular, quadratic systems. The
authors in [26] use polynomial recasting of nonlinear dynamical systems
together with sum-of-squares decomposition to prove stability of nonlinear
dynamical systems. The result is illustrated on numerical examples with at
most two DoFs. In [25], a fast sampling approach to estimate the domain of
attraction is used. While this direct approach works well in lower dimensions
(examples up to third order), scalability issues for higher dimensions remain.
In [34] it is shown that the arc length function is a maximal Lyapunov
function (particularly, it is a Lyapunov function inside the domain of
attraction, and approaches infinity on the boundary). A rather expensive
computational procedure for approximating the arc length function is proposed,
and illustrated on numerical examples with at most three DoFs. The authors in
[11] parametrize the Lyapunov function used in the stability proof, and
construct the domain of attraction as a union of (potentially infinitely many)
Lyapunov functions. While results in two and three DoF systems are shown, the
authors point out that scalability will be a major issue. For quadratic
systems, [13] present a conservative ellipsoidal estimate of the domain of
attraction. The method enlarges the ellipsoidal shape of the Lyapunov function
by evaluating the $2^{n-1}$ corners of a polytope in $\mathbb{R}^{n}$. In
[22], an analytical method for estimating the domain of attraction of
polynomial systems is proposed via parametrized quadratic Lyapunov functions.
The method, however, provides only conservative analytical estimates of the
stability domain. Examples of at most two DoFs are presented. The authors in
[23] present a conservative stability analysis for linear systems with
quadratic state feedback controllers, where switched feedback is shown to
improve the transient system response. The considered systems are only of two
DoFs. The domain of attraction can virtually have any shape, and while
ellipsoidal, circular, polytope solutions are common, they are also
conservative. The authors in [28] propose a method that enlarges an initial
guess of the stability domain by using fuzzy polynomials together with sum-of-
squares techniques. However, only examples with two DoFs are presented. The
emergence of sum-of-squares techniques for convex programming has also led to
improved algorithms for stability analysis for polynomial systems, mostly with
quadratic Lyapunov functions, see [12, 2] for more details. As evidenced by
the above literature, stability domain computations have largely focused on
very low-dimensional ($\leq 4$ DoFs) polynomial systems, and scalability
remains a concern.

This work enables quantitative analysis of stability domains for projection-
based and fully data-driven ROMs via a computationally tractable optimization
algorithm. The analysis of accurate ROMs of high-dimensional systems is a
first step towards stability analysis for complex high-dimensional systems.
Our focus is solely on scalable approaches that work well in those (higher,
yet not large-scale) ROM dimensions of $\mathcal{O}(10)$. In particular, we
propose to use a convex optimization-based approach from [33] to compute
estimates of the stability domain for quadratic-bilinear systems. Our
contributions are thus threefold: First, we derive a new analytical estimate
of the stability domain, which is rather conservative but can be evaluated
efficiently, even in the case of high-dimensional systems. Second, with the
goal to enlarge this estimate, we show that our efficient implementation of
the convex optimization-based approach from [33] scales well for
$\mathcal{O}(10)$-dimensional ROMs, which goes well beyond the four
dimensional systems analyzed in the literature. This allows us to certify ROM
simulations with respect to stability domains centered around equilibria.
Third, we demonstrate the methods’ flexibility on several different
(projection-based and fully data-driven) ROMs with up to 21 degrees of
freedom. This provides us with new insight into the stability characteristics
of those ROM techniques. We do not assume a specific ROM method, so a
particular appeal is that the approach is applicable to quadratic-bilinear
models obtained via data-driven approaches, where ROM stability properties
cannot—per definition—be derived from the full-order model.

This paper is organized as follows. Section 2 presents stability definitions
and some necessary background material. Section 3 derives an analytical
estimate for the stability domain, and Section 4 presents an optimization-
based approach to obtain larger, less-conservative estimates. Section 5
illustrates our findings on three test problems with different ROMs for semi-
discretized PDE systems. Section 6 offers conclusions and an outlook to future
work.

##  2 Background: Stability analysis for quadratic-bilinear systems

We present necessary background material and define the problem under
consideration in this paper. Section 2.1 introduces the specific model form of
quadratic-bilinear (QB) systems, Section 2.2 discusses equilibrium solutions,
and Section 2.3 defines the domain of attraction (also called stability
domain, region of attraction, or basin of attraction).

###  2.1 Quadratic-bilinear systems

Consider a quadratic-bilinear (QB) system of the form

(1)
$\displaystyle\mathbf{E}\dot{\mathbf{x}}=\mathbf{A}\mathbf{x}+\mathbf{H}(\mathbf{x}\otimes\mathbf{x})+\sum_{i=1}^{n}\mathbf{N}_{i}\mathbf{x}\mathbf{u}_{i}+\mathbf{B}\mathbf{u}$

with state $\mathbf{x}=\mathbf{x}(t)\in\mathbb{R}^{n}$, initial condition
$\mathbf{x}(0)=\mathbf{0}$, input $\mathbf{u}=\mathbf{u}(t)\in\mathbb{R}^{m}$,
matrices $\mathbf{E}\in\mathbb{R}^{n\times n},\
\mathbf{A}\in\mathbb{R}^{n\times n},\ \mathbf{H}\in\mathbb{R}^{n\times
n^{2}},\ \mathbf{B}\in\mathbb{R}^{n\times m}$, and
$\mathbf{N}_{i}\in\mathbb{R}^{n\times n}$ for $i=1,2,\ldots,m$. Here, the
matrix $\mathbf{E}$ is assumed non-singular, and the symbol $\otimes$ denotes
the standard Kronecker product. Without loss of generality, we assume the
matrix $\mathbf{H}$ is symmetric in that
$\mathbf{H}(\mathbf{x}_{1}\otimes\mathbf{x}_{2})=\mathbf{H}(\mathbf{x}_{2}\otimes\mathbf{x}_{1})$
for any $\mathbf{x}_{1},\mathbf{x}_{2}\in\mathbb{R}^{n}$, which can be
enforced without changing the dynamics, see, e.g., [4].

To ease notation, the material in the following two sections on stability
domains is presented for a generic $n$-dimensional system. In Section 5 we
introduce the reduced-order modeling context formally, which results in low-
order models for which these methods can be applied.

###  2.2 Stability of equilibrium solutions and domain of attraction

To investigate the stability of equilibrium solutions for QB systems, consider
the system

(2)
$\displaystyle\mathbf{E}\dot{\mathbf{x}}=\mathbf{A}\mathbf{x}+\mathbf{H}(\mathbf{x}\otimes\mathbf{x}),$

and invertible matrix $\mathbf{E}$. The above systems includes the autonomous
case, i.e, equation (1) with $\mathbf{u}\equiv\mathbf{0}$ and also the case
with state-dependent feedback $\mathbf{u}=\mathcal{K}\mathbf{x}$ (with proper
redefinition of the matrices $\mathbf{A},\mathbf{H}$). It is clear that
$\mathbf{x}_{e}=\mathbf{0}$ is an equilibrium. Moreover, it is sufficient to
study the zero equilibrium of the QB system. To see this, note that all
equilibrium solutions are given by $\mathbf{x}_{e}\neq\mathbf{0}$ that satisfy

(3)
$\displaystyle\mathbf{0}=\mathbf{A}\mathbf{x}_{e}+\mathbf{H}(\mathbf{x}_{e}\otimes\mathbf{x}_{e}).$

Introducing $\mathbf{z}=\mathbf{x}-\mathbf{x}_{e}$ and inserting into equation
(2) gives

(4) $\displaystyle\mathbf{E}\dot{\mathbf{z}}$
$\displaystyle=\mathbf{A}\mathbf{z}+\mathbf{H}(\mathbf{z}\otimes\mathbf{z})+2\mathbf{H}(\mathbf{x}_{e}\otimes\mathbf{z})+\underbrace{\mathbf{A}\mathbf{x}_{e}+\mathbf{H}(\mathbf{x}_{e}\otimes\mathbf{x}_{e})}_{=\mathbf{0}}$
(5)
$\displaystyle=[\mathbf{A}+2\mathbf{H}(\mathbf{I}\otimes\mathbf{x}_{e})]\mathbf{z}+\mathbf{H}(\mathbf{z}\otimes\mathbf{z}).$

The resulting system is again quadratic and has zero as an equilibrium.
Therefore, without loss of generality we study the zero equilibrium of the QB
system, and note that the theory carries over to nonzero equilibria after the
model is shifted as shown above.

###  2.3 Domain of attraction

For nonlinear systems, stability of equilibrium point requires a local
characterization–whereas for linear systems, stability of the system matrix
automatically guarantees global stability. The domain of attraction (DA), also
called stability domain, is the set of all initial conditions that result in
bounded trajectories that converge to the equilibrium solution. Let
$\phi(\mathbf{x}_{0},t)$ be a solution to (2) for a given initial condition
$\mathbf{x}_{0}$. The domain of attraction for the equilibrium
$\mathbf{x}_{e}$ is defined as

(6)
$\displaystyle\mathcal{A}(\mathbf{x}_{e}):=\left\\{\mathbf{x}_{0}:\lim_{t\rightarrow\infty}\phi(\mathbf{x}_{0},t)=\mathbf{x}_{e}\right\\}.$

The domain of attraction can have complicated geometric structure. Computing
the domain of attraction analytically is impossible even in very simple cases.
A common solution to this problem is to lower-bound the DA by a set
$\mathcal{D}\subseteq\mathcal{A}$ of simple geometric shapes, such as
ellipsoids or polyhedrons. The task then becomes to make this lower bound as
tight as possible.

We next state the well-known theorem that relates local stability to a
suitable Lyapunov function, and characterizes the domain of attraction.

######  Theorem 2.1.

[33, 11] If there exists a Lyapunov function
$v(\cdot):\mathbb{R}^{n}\mapsto\mathbb{R}^{+}$ such that

$v(\mathbf{x})>0,\qquad\dot{v}(\mathbf{x})<0$

for all $\mathbf{x}$ in a neighborhood of zero, then the zero solution is
locally asymptotically stable. Moreover,

(7) $\mathcal{D}(\rho)=\\{\mathbf{x}:v(\mathbf{x})\leq\rho^{2},\
\dot{v}(\mathbf{x})<0\\}.$

is an estimate of the domain of attraction for the zero equilibrium.

This theorem presents the basis for our computational algorithms following in
the next section. Finding the largest $\rho$ that satisfies (7) can then be
shown to require the solution of a convex optimization problem.

##  3 Analytical estimates of the stability domain

Lyapunov functions are the standard tool for stability analysis of nonlinear
systems, and hence for approximating the DA. The choice of the Lyapunov
function is non-trivial and influences the DA approximation as the level sets
of the Lyapunov function are needed in DA estimates. In stability analysis for
quadratic systems, Lyapunov functions such as polyhedral functions [1],
polynomials of higher degree [11], and arc length function approximations [34]
have been used. Here, we present an analytical approach to approximate the
domain of attraction. Quadratic Lyapunov functions lead to ellipsoidal
estimates of the DA [13, 22] which are easy to compute, and provide an initial
conservative estimate of the DA.

A Lyapunov matrix is a positive definite matrix $\mathbf{P}$ that satisfies

(8)
$\mathbf{A}^{\top}\mathbf{P}\mathbf{E}+\mathbf{E}^{\top}\mathbf{P}\mathbf{A}+\mathbf{Q}=\mathbf{0}$

for some positive definite matrix
$\mathbf{Q}=\mathbf{Q}_{f}^{\top}\mathbf{Q}_{f}$. Given a Lyapunov matrix, we
define a quadratic, nonnegative Lyapunov function via

(9)
$v(\mathbf{x})=\mathbf{x}^{\top}\mathbf{E}^{\top}\mathbf{P}\mathbf{E}\mathbf{x},$

with a derivative along trajectories as

(10)
$\displaystyle\dot{v}(\mathbf{x})=\dot{\mathbf{x}}^{\top}\mathbf{E}^{\top}\mathbf{P}\mathbf{x}+\mathbf{x}^{\top}\mathbf{P}\mathbf{E}\dot{\mathbf{x}}.$

######  Proposition 3.1.

Let $\mathbf{A}$ be Hurwitz and $\mathbf{P}$ be a Lyapunov matrix, i.e., a
solution to equation (8). Let
$v(\mathbf{x})=\mathbf{x}^{\top}\mathbf{E}^{\top}\mathbf{P}\mathbf{E}\mathbf{x}$
be a Lyapunov function. Then $\mathbf{x}_{e}=\mathbf{0}$ is a locally stable
equilibrium and $\mathcal{D}(\rho)\subseteq\mathcal{A}(\mathbf{0})$ with
$\rho=\frac{\sigma^{2}_{\text{min}}(\mathbf{Q}_{f})}{2\|\mathbf{H}\|_{2}\sqrt{\|\mathbf{P}\|_{2}}}$
is an estimate of the domain of attraction.

######  Proof 3.2.

Since $\mathbf{P}$ is symmetric positive definite, and $\mathbf{E}$ is
positive definite, the Lyapunov function
$v(\mathbf{x})=\mathbf{x}^{\top}\mathbf{E}^{\top}\mathbf{P}\mathbf{E}\mathbf{x}$
is positive. We next consider the region where its derivative is negative. We
have

(11) $\displaystyle\dot{v}(\mathbf{x})$
$\displaystyle=\dot{\mathbf{x}}^{\top}\mathbf{E}^{\top}\mathbf{P}\mathbf{E}\mathbf{x}+\mathbf{x}^{\top}\mathbf{E}^{\top}\mathbf{P}\mathbf{E}\dot{\mathbf{x}},$
(12)
$\displaystyle=[\mathbf{A}\mathbf{x}+\mathbf{H}(\mathbf{x}\otimes\mathbf{x})]^{\top}\mathbf{P}\mathbf{E}\mathbf{x}+\mathbf{x}^{\top}\mathbf{E}^{\top}\mathbf{P}[\mathbf{A}\mathbf{x}+\mathbf{H}(\mathbf{x}\otimes\mathbf{x})],$
(13)
$\displaystyle=\mathbf{x}^{\top}[\mathbf{A}^{\top}\mathbf{P}\mathbf{E}+\mathbf{E}^{\top}\mathbf{P}\mathbf{A}]\mathbf{x}+(\mathbf{x}^{\top}\otimes\mathbf{x}^{\top})\mathbf{H}^{\top}\mathbf{P}\mathbf{E}\mathbf{x}+\mathbf{x}^{\top}\mathbf{E}^{\top}\mathbf{P}\mathbf{H}(\mathbf{x}\otimes\mathbf{x}),$
(14)
$\displaystyle=-\mathbf{x}^{\top}\mathbf{Q}_{f}^{\top}\mathbf{Q}_{f}\mathbf{x}+2\mathbf{x}^{\top}\mathbf{E}^{\top}\mathbf{P}\mathbf{H}(\mathbf{x}\otimes\mathbf{x}),$
(15)
$\displaystyle\leq-\sigma^{2}_{\text{min}}(\mathbf{Q}_{f})\|\mathbf{x}\|_{2}^{2}+2\|\mathbf{x}\|_{2}^{3}\|\mathbf{E}\|_{2}\|\mathbf{P}\|_{2}\|\mathbf{H}\|_{2}.$

where $\sigma_{\text{min}}(\mathbf{Q}_{f})$ is the smallest non-zero singular
value of $\mathbf{Q}_{f}$. Thus we have that

(16) $\dot{v}(\mathbf{x})<0\ \ \Leftrightarrow\ \
\|\mathbf{x}\|_{2}<\frac{\sigma^{2}_{\text{min}}(\mathbf{Q}_{f})}{2\|\mathbf{E}\|_{2}\|\mathbf{P}\|_{2}\|\mathbf{H}\|_{2}}.$

This shows that the zero equilibrium is locally asymptotically stable. To get
an estimate of the domain of attraction (7), we need to find the radial upper
bound on the Lyapunov function for all $\mathbf{x}$ where
$\dot{v}(\mathbf{x})<0$. In particular,

(17)
$v(\mathbf{x})=\mathbf{x}^{\top}\mathbf{E}^{\top}\mathbf{P}\mathbf{E}\mathbf{x}\leq\|\mathbf{x}\|_{2}^{2}\|\mathbf{P}\|_{2}\|\mathbf{E}\|_{2}^{2}<\frac{\sigma^{4}_{\text{min}}(\mathbf{Q}_{f})}{4\|\mathbf{P}\|_{2}\|\mathbf{H}\|^{2}_{2}}=\rho^{2}$

which yields the claimed result after taking a square root.

######  Remark 3.3.

Several observations are in order. First, for $\mathbf{H}=\mathbf{0}$ we see
that $\rho=\infty$, so the stability region is $\mathbb{R}^{n}$, which is
consistent with the assumption that $\mathbf{A}$ is Hurwitz. Second, note that
$\rho$ is always nonzero as $\sigma_{\text{min}}(\mathbf{Q}_{f})$ is the
smallest non-zero singular value of $\mathbf{Q}_{f}$, and the Lyapunov matrix
is positive definite for the right-hand-side
$\mathbf{Q}_{f}^{\top}\mathbf{Q}_{f}$. Therefore $\mathbf{x}_{e}=\mathbf{0}$
is guaranteed to be a stable solution of the QB system as long as $\mathbf{A}$
is Hurwitz. However, the size of the corresponding stability domain (indicated
by $\rho$) intuitively is inversely proportional to $\|\mathbf{H}\|_{2}$. For
weak nonlinear terms, the estimate of the stability domain is larger, and for
strong nonlinear terms, this estimate shrinks. Third, as with most analytical
approaches, the result in Proposition 3.1 is rather conservative.
Computational approaches can enlarge the estimate of the stability domain,
which we discuss next.

##  4 Estimating the stability domain via optimization

We present an optimization-based approach to enlarge the estimate of the
stability domain. For the purpose of this analysis we rewrite the QB system in
the form

(18)
$\mathbf{E}\dot{\mathbf{x}}=\mathbf{A}\mathbf{x}+\sum_{i=1}^{n}x_{i}\mathbf{K}_{i}\mathbf{x},$

with $\mathbf{x}(0)=\mathbf{0}$ and where $x_{i}$ is the $i$th component of
the vector $\mathbf{x}$, $\mathbf{K}_{i}\in\mathbb{R}^{n\times n}$ are
matrices such that $\mathbf{K}_{i}=\mathbf{H}(:,(i-1)n+1:in)$ for
$i=1,2,\ldots n$. We observe that the quadratic part of the right-hand side is
invariant under skew-symmetric matrix additions, namely

$\sum_{i=1}^{n}x_{i}\mathbf{K}_{i}\mathbf{x}=\sum_{i=1}^{n}x_{i}[\mathbf{K}_{i}+\mathbf{S}_{i}]\mathbf{x}$

for any skew-symmetric matrices $\mathbf{S}_{i}$. Therefore, the nonlinearity
can be parametrized with the vector

(19)
$\boldsymbol{\mu}=\text{vec}([\mathbf{S}_{1},\mathbf{S}_{2},\ldots,\mathbf{S}_{n}])\in\mathbb{R}^{d_{n}},\qquad
d_{n}=n^{2}(n-1)/2.$

where vec() stacks the columns of the matrix into a column vector. We then
define

(20)
$\mathbf{F}(\mathbf{x},\boldsymbol{\mu})=\sum_{i=1}^{n}x_{i}[\mathbf{K}_{i}+\mathbf{S}_{i}]=\sum_{i=1}^{n}x_{i}\mathbf{M}_{i}(\boldsymbol{\mu}),$

and note that
$\mathbf{F}(\mathbf{x},\boldsymbol{\mu}_{1})\mathbf{x}=\mathbf{F}(\mathbf{x},\boldsymbol{\mu}_{2})\mathbf{x}$
for $\boldsymbol{\mu}_{1}\neq\boldsymbol{\mu}_{2}$ but
$\mathbf{F}(\mathbf{x},\boldsymbol{\mu}_{1})\neq\mathbf{F}(\mathbf{x},\boldsymbol{\mu}_{2})$.
The parametrized dynamics of the quadratic system can be written as

(21)
$\mathbf{E}\dot{\mathbf{x}}=\mathbf{A}\mathbf{x}+\mathbf{F}(\mathbf{x},\boldsymbol{\mu})\mathbf{x}.$

This parametrized formulation is the basis for the optimization routine to
estimate the stability domain.

In order to compute a less conservative estimate of the stability domain, we
define a program to maximize the stability radius $\rho$, as in [33]. The
optimization problem becomes

(22) $\displaystyle{\rho^{*}}^{2}$
$\displaystyle=\inf_{\mathbf{x}\in\mathbb{R}^{n}}\
\mathbf{x}^{\top}\mathbf{E}^{\top}\mathbf{P}\mathbf{E}\mathbf{x},$ s.t.
$\displaystyle\quad\dot{v}(\mathbf{x})=0.$

The derivative of $v(\cdot)$ along trajectories of equation (21) is

(23) $\displaystyle\dot{v}(\mathbf{x})$
$\displaystyle=\dot{\mathbf{x}}^{\top}\mathbf{E}^{\top}\mathbf{P}\mathbf{E}\mathbf{x}+\mathbf{x}^{\top}\mathbf{E}^{\top}\mathbf{P}\mathbf{E}\dot{\mathbf{x}}$
(24)
$\displaystyle=[\mathbf{A}\mathbf{x}+\mathbf{F}(\mathbf{x},\boldsymbol{\mu})\mathbf{x}]^{\top}\mathbf{P}\mathbf{E}\mathbf{x}+\mathbf{x}^{\top}\mathbf{E}^{\top}\mathbf{P}[\mathbf{A}\mathbf{x}+\mathbf{F}(\mathbf{x},\boldsymbol{\mu})\mathbf{x}]$
(25)
$\displaystyle=\mathbf{x}^{\top}\left[-\mathbf{Q}+\sum_{i=1}^{n}x_{i}[\mathbf{M}_{i}(\boldsymbol{\mu})^{\top}\mathbf{P}\mathbf{E}+\mathbf{E}^{\top}\mathbf{P}\mathbf{M}_{i}(\boldsymbol{\mu})]\right]\mathbf{x}$

for the positive definite matrix
$-(\mathbf{A}^{\top}\mathbf{P}\mathbf{E}+\mathbf{E}^{\top}\mathbf{P}\mathbf{A})=\mathbf{Q}=\mathbf{Q}_{f}^{\top}\mathbf{Q}_{f}$.
Let $\mathbf{P}=\mathbf{P}_{f}^{\top}\mathbf{P}_{f}$, and define the matrices

$\mathbf{G}(\boldsymbol{\mu})=\left[\mathbf{M}_{1}(\boldsymbol{\mu})^{\top}\mathbf{P}\mathbf{E}+\mathbf{E}^{\top}\mathbf{P}\mathbf{M}_{1}(\boldsymbol{\mu})|\ldots|\mathbf{M}_{n}(\boldsymbol{\mu})^{\top}\mathbf{P}\mathbf{E}+\mathbf{E}^{\top}\mathbf{P}\mathbf{M}_{n}(\boldsymbol{\mu})\right]^{\top}\in\mathbb{R}^{n^{2}\times
n}$

as well as

(26)
$\mathbf{J}(\boldsymbol{\mu})=(\mathbf{P}_{f}^{\top}\otimes\mathbf{Q}_{f}^{\top})^{-1}\mathbf{G}(\boldsymbol{\mu})\mathbf{Q}_{f}^{-1}\in\mathbb{R}^{n^{2}\times
n}$

which is affine linear in the parameters $\boldsymbol{\mu}$, a property that
is important for optimization later111If the Lyapunov equation (8) is solved
with low-rank methods, see e.g., the survey [8], then
$\mathbf{P}_{f}^{\top}\in\mathbb{R}^{n\times n_{r}}$ with $n_{r}\ll n$, in
which case the inverse in (26) becomes a pseudo-inverse.. Taken together, we
have

(27)
$\dot{v}(\mathbf{x})=\mathbf{x}^{\top}\mathbf{Q}_{f}^{\top}\left[-\mathbf{I}_{n}+(\mathbf{x}^{\top}\mathbf{P}_{f}^{\top}\otimes\mathbf{I}_{n})\mathbf{J}(\boldsymbol{\mu})\right]\mathbf{Q}_{f}\mathbf{x}.$

The next theorem then addresses the solution of the optimization problem (22).

######  Theorem 4.1.

[33] For a given Lyapunov matrix $\mathbf{P}$, a guaranteed stability domain
$\mathcal{D}(\rho)\subseteq\mathcal{A}(\mathbf{0})$ is given through

(28) $\rho<\rho^{*}=\frac{1}{\alpha^{*}},$

where

(29)
$\alpha^{*}=\min_{\boldsymbol{\mu}\in\mathbb{R}^{d_{n}}}\|\mathbf{J}(\boldsymbol{\mu})\|_{2}.$

Moreover, the estimate is optimal, i.e., $\rho^{*}$ is the maximum achievable
value of $\rho$ in equation (7).

We observe that the optimization problem (29) is convex, as it involves the
two-norm computation of an matrix that is affine linear in $\boldsymbol{\mu}$.
We can therefore employ special convex optimization solvers, which makes this
problem formulation appealing, despite the polynomial growth (see (19)) of the
optimization variables with respect to the dimension of the model.

##  5 Numerical results

We present three different test cases, spanning different PDE models, as well
as different techniques to obtain the ROMs. This shows the wide applicability
of this framework, as it is agnostic to the way the ROM was obtained. We
consider Burgers’ equation in Section 5.1 for which we obtain ROMs through
LQG-balanced truncation, and in Section 5.2 we use POD projection-based model
reduction to obtain ROMs for the FitzHugh-Nagumo system. Section 5.3 then
illustrates our framework on a ROM that was learned purely from data, where
that data was generated from yet a different configuration of Burgers’
equation.

To solve the optimization problem (29), we use Matlab’s fmincon optimizer with
relative tolerances TolX=0.1, tolFun=0.001, MaxIter=1000, random initial
condition $\boldsymbol{\mu}_{0}$, and we bound the entries
$\boldsymbol{\mu}^{(i)}$ of the vector $\boldsymbol{\mu}$ within the
optimization as $-10^{4}\leq\boldsymbol{\mu}^{(i)}\leq 10^{4}$.

###  5.1 LQG-balanced ROM for Burgers’ equation

####  5.1.1 Burgers’ equation and discretization

We consider the one-dimensional Burgers’ equation following the setup in [9].
The PDE model is

(30) $\dot{z}(\xi,t)=\epsilon
z_{\xi\xi}(\xi,t)-\frac{1}{2}(z^{2}(\xi,t))_{\xi}+\sum_{k=1}^{m}\chi_{[(k-1)/m,k/m]}(\xi)u_{k}(t)$

for $t>0$, $\xi\in[0,1]$ is the spatial variable and
$z(\xi,0)=z_{0}(\xi)=0.5\sin(2\pi\xi)^{2}$ for $\xi\in[0,0.5]$ and zero
otherwise. The notation
$z_{\xi\xi}(\xi,t):=\frac{\partial^{2}}{\partial\xi^{2}}z(\xi,t)$ denotes a
second order spatial derivative; similarly, $z_{\xi}(\xi,t)$ denotes a first
spatial derivative. Moreover, $z\in H_{\text{per}}^{1}(0,1)$, which means that
the system has periodic boundary conditions. Here, $u_{i}(t),i=1,2,\ldots,m$
are the input functions. The function $\chi_{[a,b]}(x)$ denotes the
characteristic function on $[a,b]$. Practically, this means that the spatial
domain is subdivided into $m$ intervals of equal length, and each control is
applied in one of the corresponding intervals. When discretized with linear
finite elements, the $N$-dimensional semi-discretized system has the form

(31) $\displaystyle\mathbf{E}_{N}\dot{\mathbf{x}}_{N}$ $\displaystyle=$
$\displaystyle\mathbf{A}_{N}\mathbf{x}+\mathbf{H}_{N}(\mathbf{x}_{N}\otimes\mathbf{x}_{N})+\mathbf{B}_{N}\mathbf{u}\qquad\mathbf{x}_{N}(0)=\mathbf{x}_{N,0}$
(32) $\displaystyle\mathbf{y}_{N}$ $\displaystyle=$
$\displaystyle\mathbf{C}_{N}\mathbf{x}_{N}.$

The output matrix $\mathbf{C}_{N}$ produces an observation $\mathbf{y}_{N}$ of
the system. Here, we observe the entire state, so
$\mathbf{C}_{N}=\mathbf{I}_{N}$. We choose $m=3$ inputs and the viscosity is
set to $\epsilon=10^{-3}$ to make the nonlinear quadratic term dominant. The
model dimension is $N=101$, which is still not tractable for stability domain
computation.

####  5.1.2 LQG-balanced reduced-order model

The ROM for Burgers’ equation is obtained through the LQG-balancing framework
[19], which also allows for finding a proper energy function. For linear time-
invariant systems, the LQG-balancing framework finds a coordinate
transformation (and subsequent reduction) so that the LQG solutions to the
algebraic Riccati equations of the ROM are equal and diagonal. This yields a
ROM with favorable control-theoretic properties. We outline the main steps of
this well-known method here, and refer to [19, 3] for details. First, we
compute the solutions to the LQG algebraic Riccati equations based on the
linearized system matrices, i.e.,

(33)
$\displaystyle\mathbf{A}_{N}\mathbf{P}_{N}+\mathbf{P}_{N}\mathbf{A}_{N}^{\top}-\mathbf{P}_{N}\mathbf{C}_{N}^{\top}\mathbf{C}_{N}\mathbf{P}_{N}+\mathbf{B}_{N}\mathbf{B}_{N}^{\top}$
$\displaystyle=\mathbf{0},$
$\displaystyle\mathbf{A}_{N}^{\top}\mathbf{Q}_{N}+\mathbf{Q}_{N}\mathbf{A}_{N}-\mathbf{Q}_{N}\mathbf{B}_{N}\mathbf{B}_{N}^{\top}\mathbf{Q}_{N}+\mathbf{C}_{N}^{\top}\mathbf{C}_{N}$
$\displaystyle=\mathbf{0}.$

Since $\mathbf{P}_{N}$ and $\mathbf{Q}_{N}$ are symmetric positive definite
they can be used to define a Lyapunov function in (9), in a similar way as the
symmetric positive definite solutions to the linear Lyapunov equation are
used. The transformation requires the Cholesky factors
$\mathbf{P}_{N}=\mathbf{R}_{N}\mathbf{R}_{N}^{\top}$ and
$\mathbf{Q}_{N}=\mathbf{L}_{N}\mathbf{L}_{N}^{\top}$, and then the singular
value decomposition

$\mathbf{U}_{N}\boldsymbol{\Sigma}_{N}\mathbf{W}_{N}=\mathbf{L}_{N}^{\top}\mathbf{R}_{N}.$

Based on the singular value decay, we choose a truncation order $n\ll N$, and
define $\boldsymbol{\Sigma}_{n}=\boldsymbol{\Sigma}_{N}(1:n,1:n)$. The LQG-BT
projection matrices are
$\mathbf{T}_{n}=\mathbf{R}_{N}(:,1:n)\mathbf{W}_{N}(:,1:n)\boldsymbol{\Sigma}_{n}^{-\frac{1}{2}}$
and
$\mathbf{T}_{n}^{-1}=\boldsymbol{\Sigma}_{n}^{-\frac{1}{2}}\mathbf{U}_{N}(:,1:n)^{\top}\mathbf{L}_{N}(:,1:n)^{\top}$,
which yields the ROM matrices in $n$-dimensions:

(34)
$\mathbf{A}=\mathbf{T}_{n}^{-1}\mathbf{A}_{N}\mathbf{T}_{n},\qquad\mathbf{B}=\mathbf{T}_{n}^{-1}\mathbf{B}_{N},\quad\mathbf{C}=\mathbf{C}_{N}\mathbf{T}_{n},\quad\mathbf{H}=\mathbf{T}_{n}^{-1}\mathbf{H}_{N}(\mathbf{T}_{n}\otimes\mathbf{T}_{n}).$

Note that the ROM is a quadratic model of dimension $n$ and we only used the
matrices of the linearized system to compute the transformations. Per
definition of the LQG balancing transformation, the matrices in reduced
dimensions satisfy

(35)
$\displaystyle\mathbf{A}\boldsymbol{\Sigma}_{n}+\boldsymbol{\Sigma}_{n}\mathbf{A}^{\top}-\boldsymbol{\Sigma}_{n}\mathbf{C}^{\top}\mathbf{C}\boldsymbol{\Sigma}_{n}+\mathbf{B}\mathbf{B}^{\top}$
$\displaystyle=\mathbf{0},$
$\displaystyle\mathbf{A}^{\top}\boldsymbol{\Sigma}_{n}+\boldsymbol{\Sigma}_{n}\mathbf{A}-\boldsymbol{\Sigma}_{n}\mathbf{B}\mathbf{B}^{\top}\boldsymbol{\Sigma}_{n}+\mathbf{C}^{\top}\mathbf{C}$
$\displaystyle=\mathbf{0}.$

Figure 1, left, shows the singular values for the matrix
$\mathbf{L}_{N}^{\top}\mathbf{R}_{N}$, which are typically used to decide
about the truncation of the system. The singular values decay rather slowly at
first, but have a steep drop at $n=20$, after which they become machine zero.
Therefore, increasing $n$ beyond twenty should not be expected to yield better
results ROMs. We compute ROMs of dimensions $n=3,5,7,\ldots,21$, for which we
analyze the stability domain in the next section.

####  5.1.3 Stability domain computation

We have seen that $\boldsymbol{\Sigma}_{n}$ satisfies the LQG Riccati
equations (35), i.e., it is symmetric positive definite (trivially so, as a
diagonal matrix). Thus, we can choose $\mathbf{P}=\boldsymbol{\Sigma}_{n}$ to
define a Lyapunov function
$v(\mathbf{x})=\mathbf{x}^{\top}\boldsymbol{\Sigma}_{n}\mathbf{x}$, and
moreover $\mathbf{P}_{f}=\boldsymbol{\Sigma}_{n}^{1/2}$ and
$\mathbf{Q}_{f}=\mathbf{C}$. We then solve the optimization problem (22) to
compute the optimal estimate $\rho^{*}$. Figure 1, right, shows the stability
radii $\rho$ obtained from the analytical estimate in Proposition 3.1, and the
optimal estimate $\rho^{*}$ from Theorem 4.1.

We observe that the optimized stability domain is several orders of magnitude
larger than the conservative analytical estimate from Proposition 3.1.
Furthermore, we observe that while the estimate $\rho^{*}$ initially reduces
significantly until $n=13$, the stability region again increases until $n=21$.
This hints at the fact that as the ROM increases in fidelity, it also becomes
more stable. This observation is not possible by considering only the
analytical estimate $\rho$, and is therefore a major appeal to using the
optimized estimate. The analytic estimate from Proposition 3.1 for the full-
order model is $\rho_{\text{FOM}}=9.81\times 10^{-4}$. Note, that the ROM-
based analytic estimate appears to converge to zero, therefore not approaching
the FOM analytical estimate. We observe that the optimization-based result for
the ROM is much closer to $\rho_{\text{FOM}}$. Due to the high dimension of
the state space, the optimization-based estimate $\rho^{*}_{\text{FOM}}$ was
not computationally tractable.

We note that we tried the SMRSOFT toolbox
https://www.eee.hku.hk/~chesi/y_smrsoft.htm, but for $n\in\\{2,3\\}$ it did
not yield a lower bound other than zero, and for higher-order systems it
became very cumbersome to implement due to the toolbox not accepting matrix-
vector multiplications in right-hand-side. The toolbox seemingly was written
for low-order systems, and our work is specifically aimed at high-order
systems.

Figure 1: Burgers’ equation: (Left) Singular values of the matrix
$\mathbf{L}_{N}^{\top}\mathbf{R}_{N}$ in LQG balancing; (Right) The value of
$\rho$ for the stability domain estimate for the analytic method from
Proposition 3.1 and the optimal estimate $\rho^{*}$ from Theorem 4.1.

###  5.2 Proper orthogonal decomposition ROM for FitzHugh-Nagumo equation

####  5.2.1 FitzHugh-Nagumo equation and discretization

This section illustrates our nonlinear model reduction approach on the
FitzHugh-Nagumo system, which is a model for the activation and deactivation
of a spiking neuron. The original FitzHugh-Nagumo is a cubic system and here
we consider a lifted quadratic bilinear model, see [4, 21] for more details on
the model. The lifted QB system then reads as

$\displaystyle\epsilon\dot{v}(\xi,t)$
$\displaystyle=\epsilon^{2}v_{\xi\xi}(\xi,t)-v(\xi,t)^{3}+0.1v(\xi,t)^{2}-0.1v(\xi,t)-w(\xi,t)+c,$
$\displaystyle\dot{w}(\xi,t)$ $\displaystyle=hv(\xi,t)-\gamma w(\xi,t)+c$
$\displaystyle\dot{z}(\xi,t)$
$\displaystyle=2[\epsilon^{2}vv_{\xi\xi}(\xi,t)-z^{2}(\xi,t)+0.1z(\xi,t)v(\xi,t)-0.1z(\xi,t)-w(\xi,t)v(\xi,t)+cv(\xi,t)].$

where $\xi\in[0,L]$ is the spatial variable and the time horizon of interest
is $t\in[0,t_{f}]$. The states of the system are voltage $v(\xi,t)$ and
recovery of voltage $w(\xi,t)$. The initial conditions are specified as
$v(\xi,0)=w(\xi,0)=z(\xi,0)=0$ for $\xi\in[0,L]$, and the boundary conditions
are $v_{\xi}(0,t)=u(t)$, where $u(t)=5\times 10^{4}\ t^{3}\exp(-15t)$ so the
system is excited through the boundary; $v_{\xi}(L,t)=0$ for $t\geq 0$;
$z_{\xi}(L,t)=2v(L,t),z_{s}(0,t)=2v(0,t)u(t)$. In the problem setup we
consider, the parameters are given by $L=0.1$, $c=0.05$, $\gamma=2$, $h=0.5$,
and $\epsilon=0.015$. The PDE model is semi-discretized in space by using
finite differences, resulting in the finite-dimensional QB system

$\mathbf{E}_{N}\dot{\mathbf{x}}_{N}=\mathbf{A}_{N}\mathbf{x}_{N}+\mathbf{B}_{N}\mathbf{u}+\mathbf{H}_{N}(\mathbf{x}_{N}\otimes\mathbf{x}_{N})+\sum_{k=1}^{2}\mathbf{N}_{N,k}\mathbf{x}_{N}\mathbf{u}_{k},$

where $\mathbf{E}_{N}=\epsilon\mathbf{I}_{N}$ is diagonal,
$\mathbf{A}_{N},\mathbf{N}_{N,1},\mathbf{N}_{N,2}\in\mathbb{R}^{N\times N}$
and $\mathbf{H}_{N}\in\mathbb{R}^{N\times N^{2}}$. The input matrix is
$\mathbf{B}_{N}\in\mathbb{R}^{N\times 2}$, with the second column of
$\mathbf{B}_{N}$ being copies of $c$ (the constant in the FHN PDE) and the
first column of $\mathbf{B}_{N}$ having a 1 at the first entry. Thus, the
input $\mathbf{u}=[u(t),1]$. Here, each variable is discretized with $200$
degrees of freedom, i.e., the overall dimension of the QB model is $N=600$.

####  5.2.2 Proper orthogonal decomposition reduced-order model

We generate a ROM via the method of proper orthogonal decomposition [18] and
follow directly the implementation in [21]. We simulate the system for
$t_{f}=12s$ with Matlab’s ode15s solver with ’RelTol’,1e-8,’AbsTol’,1e-10
tolerances, and collect simulated data every 0.1s, for a total of 120
snapshots. For each of the three discretized variables we compute a POD basis
of order $n$, where we use the same $n$ for each variable. The POD basis is
the optimal basis to represent the snapshot set in the $\ell_{2}$ sense. The
projection matrix is then assembled as a block-diagonal matrix with the POD
basis for each variable as blocks. Thus, the ROM is the $3n$ dimensional ROM

$\mathbf{E}\dot{\mathbf{x}}=\mathbf{A}\mathbf{x}+\mathbf{B}\mathbf{u}+\mathbf{H}(\mathbf{x}\otimes\mathbf{x})+\sum_{k=1}^{2}\mathbf{N}_{k}\mathbf{x}\mathbf{u}_{k}.$

Figure 2, left, shows the decay of the singular values of the POD snapshot
matrices for each variable. We observe that until $n=12$ the singular values
drop significantly, after which a slow decay appears. Therefore, for $n\leq
12$, increasing $n$ should lead to improved ROM accuracy.

####  5.2.3 Stability domain computation

For the Lyapunov matrix, we choose $\mathbf{Q}=\mathbf{I}_{n}$ and compute the
corresponding solution $\mathbf{P}$ of the Lyapunov equation (8). We again
compute the stability radii $\rho$ obtained from the analytical estimate in
Proposition 3.1, and the optimal estimate $\rho^{*}$ from Theorem (4.1).
Figure 2, right, shows the obtained results, where the horizontal axis plots
the reduced dimension $3n$, so for the ROM of dimension 21, each variable was
approximated with seven POD modes. Similar to the previous example, we observe
that the optimal computed estimate is roughly three orders of magnitude larger
than the analytic one. Moreover, this time the size of the stability domain
seems to stay approximately the same after six modes are used. We also compute
the analytic estimate from Proposition 3.1 for the full-order model, and
obtain $\rho_{\text{FOM}}=2.27\times 10^{-6}$. We observe that the ROM
stability radius approximates the FOM stability radius well, within the same
order of magnitude. This, together with a measure of accuracy of predictions
could give confidence in having a good reduced-order model strategy. Due to
the high dimension of the state space, the optimization-based estimate
$\rho^{*}_{\text{FOM}}$ was not computationally tractable.

Figure 2: FitzHugh-Nagumo equation: (Left) Singular values of the POD snapshot
matrix for each variable. (Right) The value of $\rho$ for the stability domain
estimate for the analytic method from Proposition 3.1 and the optimal estimate
from Theorem 4.1.

###  5.3 Non-intrusive ROM for Burgers’ equation

In this example, we demonstrate one of the major appeals of this method,
namely that we can apply it to non-intrusive ROMs. The model presented here is
learned from data of Burgers’ equation using the operator inference method
following the setup in [27, Sec 4.2]. This setup deviates from Section 5.1 in
the discretization parameters and implementation of the control term.

####  5.3.1 Learned reduced-order model

We are given a simulator for the following Burgers’ equation PDE:

(36) $\dot{z}(\xi,t)=\epsilon
z_{\xi\xi}(\xi,t)-\frac{1}{2}(z^{2}(\xi,t))_{\xi}$

for $t>0$, $z(\cdot,0)=0$, and the input enters through the Dirichlet boundary
condition, $z(t,0)=u(t)$ and $z(t,1)=-u(t)$. We generate non-intrusive ROMs
via the operator inference222For a Matlab code for this model, and the
operator inference approach, see https://github.com/elizqian/operator-
inference. A scalable Python implementation can be found at
https://github.com/Willcox-Research-Group/rom-operator-inference-Python3.
approach from [27]. The data are generated from a finite difference solver
with equidistant grid of $N=128$. The viscosity is set to $\epsilon=10^{-1}$
and the model is simulated for $t\in[0,t_{f}]$ and $t_{f}=1$s. The one-
dimensional input $u(t)$ is generated via Matlab’s rand command with setting
rng default for reproducibility. We collect snapshots of the state every
$dt=10^{-4}$ steps and store the snapshots and the corresponding inputs as

$\mathbf{X}_{N}=[\mathbf{x}_{N}^{(0)}\
\dots\mathbf{x}_{N}^{(K)}]\in\mathbb{R}^{N\times
K},\qquad\mathbf{\mathbf{U}}=[\mathbf{u}_{0},\ \dots,\
\mathbf{u}_{K}]\in\mathbb{R}^{m\times K},$

where $\mathbf{x}_{N}^{(i)}=\mathbf{x}_{N}(t_{i})$ with
$0=t_{0}<t_{1}<\dots<t_{K}=t_{f}$. We next compute a low-dimensional POD
subspace in which to optimally represent the snapshot data, i.e., we compute
the singular value decomposition of the snapshot matrix as

$\mathbf{X}_{N}=\mathbf{V}_{N}\boldsymbol{\Sigma}_{N}\mathbf{W}_{N}^{\top},$

where $\mathbf{V}_{N}\in\mathbb{R}^{N\times K}$,
$\boldsymbol{\Sigma}\in\mathbb{R}^{K\times K}$ and
$\mathbf{W}_{N}\in\mathbb{R}^{K\times K}$. We then obtain the $n\ll
N$-dimensional POD basis as $\mathbf{V}_{n}=\mathbf{V}_{N}(:,1:n)$. The high-
dimensional snapshot data are projected onto the POD subspace spanned by the
columns of $\mathbf{V}_{n}$, which yields the projected data

$\mathbf{X}=\mathbf{V}^{\top}_{n}\mathbf{X}_{N}=[\mathbf{x}^{(0)}\quad\dots\quad\mathbf{x}^{(K)}]\in\mathbb{R}^{n\times
K},\qquad\dot{\mathbf{X}}=[\dot{\mathbf{x}}^{(0)}\quad\dot{\mathbf{x}}^{(1)}\quad\dots\quad\dot{\mathbf{x}}^{(K)}]\in\mathbb{R}^{n\times
K}.$

The columns of the time-derivative data matrix $\dot{\mathbf{X}}$ are computed
with a fourth-order implicit Runge-Kutta backward differencing method.

In order to learn the ROM, the operator inference framework solves a least
squares problem to find the reduced operators that yield the ROM that best
matches the projected snapshot data in a minimum residual sense. To learn the
Burgers’ ROM, we solve

$\min_{\mathbf{A}\in\mathbb{R}^{n\times n},\mathbf{H}\in\mathbb{R}^{n\times
n^{2}},\mathbf{B}\in\mathbb{R}^{n\times
m}}\left\|\mathbf{X}_{N}^{\top}\mathbf{A}^{\top}+(\mathbf{X}_{N}\otimes\mathbf{X}_{N})^{\top}\mathbf{H}^{\top}+\mathbf{U}^{\top}\mathbf{B}^{\top}-\dot{\mathbf{X}}_{N}^{\top}\right\|^{2}_{2}.$

This allows us to compute the ROM operators $\mathbf{A}$, $\mathbf{H}$, and
$\mathbf{B}$ without needing explicit access to the original high-dimensional
operators $\mathbf{A}_{N}$, $\mathbf{H}_{N}$, $\mathbf{B}_{N}$, which
constitutes a fully non-intrusive method. Finally, the learned ROM takes the
form

(37)
$\dot{\mathbf{x}}=\mathbf{A}\mathbf{x}+\mathbf{H}(\mathbf{x}\otimes\mathbf{x})+\mathbf{B}\mathbf{u}.$

Figure 3, left, shows the relative state reconstruction error
$\|\mathbf{X}_{N}-\mathbf{V}_{N}\mathbf{X}\|/\|\mathbf{X}_{N}\|$ for the
operator inference ROM. For $n=1,2,3$ we did not obtain stable models. The
state reconstruction error decays monotonically, and operator inference
produces accurate ROM simulations.

####  5.3.2 Stability domain computation

For the Lyapunov matrix, we choose $\mathbf{Q}=\mathbf{I}_{n}$ and compute the
corresponding solution $\mathbf{P}$ of the Lyapunov equation (8). We again
compute the stability radii $\rho$ obtained from the analytical estimate in
Proposition 3.1, and the optimal estimate $\rho^{*}$ from Theorem
(4.1).333Since this is a non-intrusive setting, we cannot compute the analytic
stability radius for the full-order model, as we only have model data, but not
the FOM operators. Figure 3, right, shows the numerical result. First, we note
that for $n>12$, the ROM simulations continued to increase in accuracy, but
the $\mathbf{A}$ matrix had few eigenvalues that were positive (in the order
$\mathcal{O}(10^{-8})$), so the approach to compute a quadratic Lyapunov
function via solution of the Lyapunov matrix 8 is not applicable. Other
approaches that work in this setting are needed, which we will discuss as part
of future work. Similar to the previous examples, we observe from Figure 3,
right, that the analytical estimate $\rho$ is two to four orders of magnitude
smaller than the estimate $\rho^{*}$ computed via optimization. Moreover, we
observe that while $\rho,\rho^{*}$ initially decrease, after $n=8$ they
stagnate and then increase, indicating improved stability properties of the
larger-dimensional ROMs.

Figure 3: Learned ROM for Burgers’ equation: (Left) Relative state
reconstruction error
$\|\mathbf{X}_{N}-\mathbf{V}_{N}\mathbf{X}\|/\|\mathbf{X}_{N}\|$. (Right) The
value of $\rho$ for the stability domain estimate for the analytic method from
Proposition 3.1 and the optimal estimate $\rho^{*}$ from Theorem 4.1.

##  6 Conclusions and future directions

For nonlinear reduced-order models (ROMs), computing the stability domain of
equilibrium points is important for both open and closed-loop applications. We
presented a framework to compute the optimal stability domain estimates for
quadratic-bilinear ROMs for a given quadratic Lyapunov function. Quadratic-
bilinear ROMs represent a large class of nonlinear systems, as many of those
systems can be recast into QB form via variable transformations and the
addition of extra variables. Our numerical findings on three different ROM
test problems show that the classical analytical estimates of the stability
domain are overly conservative—up to four orders of magnitude—compared to the
stability radii computed via the suggested convex optimization problem. The
numerical results also demonstrate various ways to pick quadratic Lyapunov
functions, which can be informed by the model reduction process itself. This
work motivates several directions of future research. First, as seen in
Section 5.3, this approach has its limitations when the linear system matrix
is unstable, since the Lyapunov equation does not have a solution. Alternative
approaches in this case are needed that can take into account the quadratic
influence on the system. Moreover, Section 5.3 suggests an interesting
direction of future research. The operator inference problem is unconstrained,
however, the stability domain computation could be integrated as an additional
constraint for the learning problem, yielding learned ROMs with favorable
stability properties. Lastly, in the projection-based setting, theoretical
results relating the ROM stability radius to the high-fidelity stability
radius would be desirable. For large-scale systems, a proper choice of norm
would like need to be made to avoid degeneration of the $l_{2}$ norm.

## References

  * [1] F. Amato, F. Calabrese, C. Cosentino, and A. Merola.  Stability analysis of nonlinear quadratic systems via polyhedral Lyapunov functions.  Automatica, 47(3):614–617, 2011.
  * [2] J. Anderson and A. Papachristodoulou.  Advances in computational Lyapunov analysis using sum-of-squares programming.  Discrete & Continuous Dynamical Systems-B, 20(8):2361, 2015.
  * [3] A. C. Antoulas.  Approximation of Large-Scale Dynamical Systems.  Advances in Design and Control. Society for Industrial and Applied Mathematics, Philadelphia, PA, USA, 2005.
  * [4] P. Benner and T. Breiten.  Two-sided projection methods for nonlinear model order reduction.  SIAM Journal on Scientific Computing, 37(2):B239–B260, 2015.
  * [5] P. Benner, A. Cohen, M. Ohlberger, and K. Willcox, editors.  Model Reduction and Approximation: Theory and Algorithms.  Computational Science & Engineering. SIAM Publications, Philadelphia, PA, 2017.
  * [6] P. Benner, P. Goyal, and S. Gugercin.  H2-quasi-optimal model order reduction for quadratic-bilinear control systems.  SIAM Journal on Matrix Analysis and Applications, 39(2):983–1032, 2018.
  * [7] P. Benner, P. Goyal, B. Kramer, P. Peherstorfer, and K. Willcox.  Operator inference for non-intrusive model reduction of systems with non-polynomial nonlinear terms.  Computer Methods in Applied Mechanics and Engineering, 372:113433, 2020.
  * [8] P. Benner and J. Saak.  Numerical solution of large and sparse continuous time algebraic matrix riccati and lyapunov equations: a state of the art survey.  GAMM-Mitteilungen, 36(1):32–52, 2013.
  * [9] J. Borggaard and L. Zietsman.  The Quadratic-Quadratic Regulator problem: Approximating feedback controls for quadratic-in-state nonlinear systems.  In 2020 American Control Conference (ACC), pages 818–823. IEEE, 2020.
  * [10] X. Cao, J. Maubach, S. Weiland, and W. Schilders.  A novel Krylov method for model order reduction of quadratic bilinear systems.  In 2018 IEEE Conference on Decision and Control (CDC), pages 3217–3222. IEEE, 2018.
  * [11] G. Chesi.  Estimating the domain of attraction via union of continuous families of Lyapunov estimates.  Systems & Control Letters, 56(4):326–333, 2007.
  * [12] G. Chesi.  Domain of attraction: analysis and control via SOS programming, volume 415.  Springer Science & Business Media, 2011.
  * [13] R. Genesio and A. Tesi.  Stability analysis of quadratic systems.  In Nonlinear Control Systems Design 1989, pages 195–199. Elsevier, 1990.
  * [14] I. V. Gosea and A. C. Antoulas.  Data-driven model order reduction of quadratic-bilinear systems.  Numerical Linear Algebra with Applications, 25(6):e2200, 2018.
  * [15] C. Gu.  QLMOR: A projection-based nonlinear model order reduction approach using quadratic-linear representation of nonlinear systems.  IEEE Transactions on Computer-Aided Design of Integrated Circuits and Systems, 30(9):1307–1320, 2011.
  * [16] L. Guillot, B. Cochelin, and C. Vergez.  A Taylor series-based continuation method for solutions of dynamical systems.  Nonlinear Dynamics, 98(4):2827–2845, 2019.
  * [17] J. S. Hesthaven, G. Rozza, B. Stamm, et al.  Certified reduced basis methods for parametrized partial differential equations, volume 590.  Springer.
  * [18] P. Holmes, J. L. Lumley, and G. Berkooz.  Turbulence, Coherent Structures, Dynamical Systems and Symmetry.  Cambridge Monographs on Mechanics. Cambridge University Press, 1996.
  * [19] E. Jonckheere and L. Silverman.  A new set of invariants for linear systems–application to reduced order compensator design.  IEEE Transactions on Automatic Control, 28(10):953–964, 1983\.
  * [20] E. H. Kerner.  Universal formats for nonlinear ordinary differential systems.  Journal of Mathematical Physics, 22(7):1366–1371, 1981.
  * [21] B. Kramer and K. Willcox.  Nonlinear model order reduction via lifting transformations and proper orthogonal decomposition.  AIAA Journal, 57(6):2297–2307, 2019.
  * [22] A. Levin.  An analytical method of estimating the domain of attraction for polynomial differential equations.  IEEE Transactions on Automatic Control, 39(12):2471–2475, 1994\.
  * [23] J.-S. Liu and S.-L. Chen.  On global stability of quadratic state feedback controlled linear systems.  Systems & control letters, 21(5):371–379, 1993.
  * [24] G. P. McCormick.  Computability of global solutions to factorable nonconvex programs: Part I–convex underestimating problems.  Mathematical Programming, 10(1):147–175, 1976.
  * [25] E. Najafi, R. Babuška, and G. A. Lopes.  A fast sampling method for estimating the domain of attraction.  Nonlinear dynamics, 86(2):823–834, 2016.
  * [26] A. Papachristodoulou and S. Prajna.  Analysis of non-polynomial systems using the sum of squares decomposition.  In Positive polynomials in control, pages 23–43. Springer, 2005\.
  * [27] B. Peherstorfer and K. Willcox.  Data-driven operator inference for nonintrusive projection-based model reduction.  Computer Methods in Applied Mechanics and Engineering, 306:196–215, 2016.
  * [28] J. L. Pitarch, A. Sala, and C. V. Arino.  Closed-form estimates of the domain of attraction for nonlinear systems via fuzzy-polynomial models.  IEEE transactions on cybernetics, 44(4):526–538, 2013.
  * [29] E. Qian, B. Kramer, B. Peherstorfer, and K. Willcox.  Lift & Learn: Physics-informed machine learning for large-scale nonlinear dynamical systems.  Physica D: Nonlinear Phenomena, 406:132401, 2020.
  * [30] A. Quarteroni and G. Rozza, editors.  Reduced Order Methods for Modeling and Computational Reduction, volume 9 of MS & A.  Springer, Berlin, Heidelberg, New York, 2014.
  * [31] M. A. Savageau and E. O. Voit.  Recasting nonlinear differential equations as S-systems: a canonical nonlinear form.  Mathematical Biosciences, 87(1):83–115, 1987.
  * [32] R. Swischuk, B. Kramer, C. Huang, and K. Willcox.  Learning physics-based reduced-order models for a single-injector combustion process.  AIAA Journal, 58:6:2658–2672, 2020.
  * [33] A. Tesi, F. Villoresi, and R. Genesio.  On the stability domain estimation via a quadratic Lyapunov function: convexity and optimality properties for polynomial systems.  IEEE Transactions on Automatic Control, 41(11):1650–1657, 1996\.
  * [34] M. Zarei, A. Kalhor, and D. Brake.  Arc length based maximal Lyapunov functions and domains of attraction estimation for polynomial nonlinear systems.  Automatica, 90:164–171, 2018.