Many classes of nonlinear systems that are described by analytic
nonlinearities can be represented equivalently as quadratic-bilinear systems
(by means of lifting). The proposed method combines the Loewner framework with
the Volterra series theory and constructs reduced quadratic-bilinear systems
from input-output time-domain data that approximate the original nonlinear
system.

# A framework for fitting quadratic-bilinear systems with applications to
models of electrical circuits

Dimitrios S. Karachalios<EMAIL_ADDRESS>Ion Victor Gosea
<EMAIL_ADDRESS>Athanasios C. Antoulas<EMAIL_ADDRESS>
###### Abstract

In this contribution, we propose a data-driven procedure to fit quadratic-
bilinear surrogate models from data. Although the dynamics characterizing the
original model are strongly nonlinear, we rely on lifting techniques to embed
the original model into a quadratic-bilinear format. Here, data represent
generalized transfer function values. This method is an extension of methods
that do bilinear, or quadratic inference, separately. It is based on first
fitting a linear model with the classical Loewner framework, and then on
inferring the best supplementing nonlinear operators, in a least-squares
sense. The application scope of this method is given by electrical circuits
with nonlinear components (such as diodes). We propose various test cases to
illustrate the performance of the method.

###### keywords:

Data-Driven Methods, Non-Intrusive Modeling, Model Reduction, Nonlinear
Dynamics, Quadratic-Bilinear Systems, System Identification, Lift and Learn
Approach.

##  1 Introduction

System Identification (SI) and data-driven Model Order Reduction (MOR) are two
fairly established methodologies that aim at discovering/computing robust
surrogates of dynamical systems from data. This is performed without having an
exact access to the system’s structure or matrices scaling the various terms.
In the case of SI, the aim of discovery of known classes of dynamical systems
that constitute an appropriate mathematical formalism capable of describing
dynamical phenomena. We refer the reader to [21] and to [19] for more details
on various methodologies. In the case of MOR, the need for approximating the
underlying dynamical system is dictated mostly by the increased dimension
under hand (the number of internal variables). Conventional MOR methods are
indeed intrusive, in the sense that they require an explicit formulation of
the dynamical system to be reduced (in terms of matrices or various
operators). We refer the reader to [7] and to [1], for more details. However,
data-driven MOR methods are generally non-intrussive, since they require only
data (snapshots of the states, input-output measurements, etc.), and not the
full/exact description of the model. Methods that fall into this category
range from DMD (dynamical mode decomposition), OpInf (operator inference), to
LF (the Loewner framework). Such data-driven methods can be used for
identifying surrogate models without having exact access to the original
operators. Starting with [14], many algorithms have been developed to identify
linear dynamical systems in the state-space realization:

$\small\left\\{\begin{aligned}
\dot{{\textbf{x}}}(t)&={\textbf{A}}{\textbf{x}}(t)+{\textbf{B}}u(t),\ \
y(t)={\textbf{C}}{\textbf{x}}(t),\end{aligned}\right.$ (1)

where x is the state variable of dimension $n$, while the system matrices are
given by ${\textbf{A}}\in{\mathbb{R}}^{n\times
n},~{}{\textbf{B}},{\textbf{C}}^{T}\in{\mathbb{R}}^{n\times 1}$. We refer to ,
for more details on various methodologies. In recent years, the ideas for
developing methods for linear systems have steadily expanded to fit certain
classes of nonlinear systems (such as polynomial). We consider nonlinear
system can be written in the following state-space representation

$\small\left\\{\begin{aligned}
\dot{{\textbf{x}}}(t)&={\textbf{A}}{\textbf{x}}(t)+{\textbf{f}}({\textbf{x}}(t),{\textbf{u}}(t))+{\textbf{B}}u(t),\
\ y(t)={\textbf{C}}{\textbf{x}}(t),\end{aligned}\right.$ (2)

where
${\textbf{f}}:{\mathbb{R}}^{n}\times{\mathbb{R}}\rightarrow{\mathbb{R}}^{n}$
is the nonlinear operator that will be approximated (using Carleman
linearization [9]) or lifted (McCormick relaxation [13]) to a polynomial
structure. In such scenarios, the scaling and superposition principles that
hold for linear models, do not hold here anymore, making the analysis fairly
difficult. Additionally, many useful phenomena of nonlinear nature can not be
accurately approximated by means of linearization methods (that could be
performed only locally). Two sub-class of polynomial models that belong to the
class described in Eq. (2) are the quadratic and the bilinear ones. The
operators that we will be concerned with in this study are given as follows:

  * •

Quadratic case:
${\textbf{f}}_{q}({\textbf{x}}(t),u(t))={\textbf{Q}}({\textbf{x}}(t)\otimes{\textbf{x}}(t))$
and

  * •

Bilinear case:
${\textbf{f}}_{b}({\textbf{x}}(t),u(t))={\textbf{N}}{\textbf{x}}(t)u(t)$,

or with linear combinations of these two classes that appear in the case of
quadratic-bilinear (QB) systems. Here, denotes the Kronecker product. For
cases when the non-linear operator of the original system is not directly
written as $\alpha f_{q}+\beta f_{b}$, we can employ lifting techniques to
embed the original nonlinear dynamics into the required format (without any
approximation whatsoever). This is performed by using specifically tailored
lifting transformations. More specifically, auxiliary variables and equations
are introduced in order to reformulate the equations in desired form. This
allows to apply conventional MOR methods to more general nonlinear systems.
Specific lifting transformations were discussed in [13, 8, 18].

A viable alternative is to employ data-driven methods is the Loewner framework
(LF), where the construction of low-order models can be achieved directly from
data. It is to be noted that LF has been recently extended to fit certain
classes of nonlinear systems from data, such as bilinear systems in [3], and
quadratic-bilinear (QB) systems in [11, 2]. However, in this methods, data
used in the computation process can not be easily inferred from practical
experiments (it is challenging to obtain in practice). Another data-driven
method that has emerged in recent years is OpInf, which uses time-domain state
measurements (snapshots of the state variable), and then fits a particular
nonlinear model (quadratic or quadratic-bilinear) by computing the appropriate
matrices. Details of the operator inference method can be found in [20] and
also in more recent works such as [6], [5].

In this contribution, we propose a data-driven procedure that can be used to
infer quadratic-bilinear surrogate models from data. This can be viewed as an
extension of the methods in [15] and in [12], which dealt with bilinear, and
quadratic inference, separately. We propose a specific application for testing
the method, i.e., electrical circuits with nonlinear components (such as
diodes). One aspect that distinguishes the OpInf framework from the extended
Loewner-based frameworks that we propose in this contribution, is that the
former requires measurements of the whole state variable. Our framework
requires only input-output measurements (transfer function measurements of
higher-order generalized transfer functions). The paper is structured as
follows; after the introduction, Section 2 introduces the class of QB systems
with its generalized transfer functions. Here, we also briefly introduce the
classical Loewner framework. Next, the newly-proposed method is introduced in
Section 3. Then, in Section 4, we first go through a simple example of a
nonlinear circuit to show various reformulation of its structure. Section 5
includes a detailed numerical study for applying the method to nonlinear
ladder circuit, while Section 6 gives the conclusions.

##  2 Quadratic-bilinear systems and the Loewner framework

###  2.1 State-space format and properties of QB systems

We analyze in what follows dynamical systems as in (2), with quadratic-
bilinear (QB) nonlinearities for which f is given as:
${\textbf{f}}({\textbf{x}}(t),u(t))={\textbf{Q}}({\textbf{x}}(t)\otimes{\textbf{x}}(t))+{\textbf{N}}{\textbf{x}}(t)u(t)$.
More precisely, let the state-space representation of such a system be given
as:

$\small\left\\{\begin{aligned}
{\textbf{E}}\dot{{\textbf{x}}}(t)&={\textbf{A}}{\textbf{x}}(t)+{\textbf{Q}}({\textbf{x}}(t)\otimes{\textbf{x}}(t))+{\textbf{N}}{\textbf{x}}(t)u(t)+{\textbf{B}}u(t),\\\
y(t)&={\textbf{C}}{\textbf{x}}(t),\end{aligned}\right.,$ (3)

where ${\textbf{x}}(0)={\textbf{x}}_{0}=\textbf{0}$ and the matrix
${\textbf{E}}\in{\mathbb{R}}^{n\times n}$ is non-singular,
${\textbf{A}}\in{\mathbb{R}}^{n\times
n},~{}{\textbf{Q}}\in{\mathbb{R}}^{n\times
n^{2}},~{}{\textbf{N}}\in{\mathbb{R}}^{n\times
n},~{}{\textbf{B}}\in{\mathbb{R}}^{n\times 1}$ and
${\textbf{C}}\in{\mathbb{R}}^{1\times n}$. Moreover, assume that Q satisfies
the property
${\textbf{Q}}({\textbf{v}}\otimes{\textbf{w}})={\textbf{Q}}({\textbf{w}}\otimes{\textbf{v}})$,
i.e., it is represented in a ”symmetrizable format”.

The first two generalized symmetric transfer functions functions of a QB
system as in (3) (sometimes referred to as Volterra kernels in the frequency
domain) are:

$\displaystyle\small H_{1}(s_{1})$
$\displaystyle={\textbf{C}}_{\ell}\boldsymbol{\Phi}(s_{1}){\textbf{B}},$ (4)
$\displaystyle H_{2}(s_{1},s_{2})$
$\displaystyle={\textbf{C}}_{\ell}\boldsymbol{\Phi}(s_{1}+s_{2}){\textbf{Q}}(\boldsymbol{\Phi}(s_{1}){\textbf{B}}\otimes\boldsymbol{\Phi}(s_{2}){\textbf{B}})$
$\displaystyle+\frac{1}{2}{\textbf{C}}\boldsymbol{\Phi}(s_{1}+s_{2}){\textbf{N}}(\boldsymbol{\Phi}(s_{1}){\textbf{B}}+\boldsymbol{\Phi}(s_{2}){\textbf{B}}).$

For more details on deriving such functions we refer the reader to [8] and to
[11, 2]. An important property of these particular functions (sometimes called
symmetric transfer functions), is that their samples can be inferred from the
spectrum of the observed output when using a purely oscillatory control input.

###  2.2 The Loewner framework

We start with a brief summary of the Loenwer framework (LF) for fitting linear
systems as in (10). For more details, we refer the reader to [4]. The starting
point for LF is to collect measurements corresponding to the (first) transfer
function, which can be inferred in practice from the first harmonic. The data
are first partitioned into two disjoint subsets, as:

$\displaystyle\begin{split}{\textrm{r}ight\
data}&:(\lambda_{j};w_{j}),~{}j=1,\ldots,k,~{}{\textrm{a}nd},\\\
{\textrm{l}eft\ data}&:(\mu_{i};v_{i}),~{}i=1,\ldots,k,\end{split}$ (5)

find the function ${\textbf{H}}(s)$, such that the following interpolation
conditions are (approximately) fulfilled:

${\textbf{H}}(\mu_{i})=v_{i},~{}~{}~{}{\textbf{H}}(\lambda_{j})=w_{j}.$ (6)

The Loewner matrix ${\mathbb{L}}\in{\mathbb{C}}^{k\times k}$ and the shifted
Loewner matrix ${\mathbb{L}}_{s}\in{\mathbb{C}}^{k\times k}$ are defined as
follows

${\mathbb{L}}_{(i,j)}=\frac{v_{i}-w_{j}}{\mu_{i}-\lambda_{j}},\
{\mathbb{L}}_{s(i,j)}=\frac{\mu_{i}v_{i}-\lambda_{j}w_{j}}{\mu_{i}-\lambda_{j}},$
(7)

while the data vectors ${\mathbb{V}},{\mathbb{W}}^{T}\in{\mathbb{R}}^{k}$ are
introduced as

${\mathbb{V}}_{(i)}=v_{i},\ \
{\mathbb{W}}_{(j)}=w_{j},~{}\text{for}~{}i,j=1,\ldots,k.$ (8)

The Loewner model is hence constructed as follows:

$\displaystyle{\textbf{E}}=-{\mathbb{L}},~{}~{}{\textbf{A}}=-{\mathbb{L}}_{s},~{}~{}{\textbf{B}}={\mathbb{V}},~{}~{}{\textbf{C}}={\mathbb{W}}.$

Provided that enough data is available, the pencil
$({\mathbb{L}}_{s},\,{\mathbb{L}})$ is often singular. In these cases, a
singular value decomposition (SVD) of Loewner matrices is needed to compute
projection matrices ${\textbf{X}}_{r},{\textbf{Y}}_{r}\in{\mathbb{C}}^{k\times
r}$. Here, $r<n$ represents the truncation index.

Then, the system matrices corresponding to a projected Loewner model of
dimension $r$ can be computed using matrices ${\textbf{X}}_{r}$ and
${\textbf{Y}}_{r}$, as:

$\small\hat{{\textbf{E}}}=-{\textbf{X}}_{r}^{*}{\mathbb{L}}{\textbf{Y}}_{r},\
\ \hat{{\textbf{A}}}=-{\textbf{X}}_{r}^{*}{\mathbb{L}}_{s}{\textbf{Y}}_{r},\ \
\hat{{\textbf{B}}}={\textbf{X}}_{r}^{*}{\mathbb{V}},\ \
\hat{{\textbf{C}}}={\mathbb{W}}{\textbf{Y}}_{r},$ (9)

and hence, directly find a state-space realization corresponding to the
reduced-order system of equations

$\left\\{\begin{aligned}
{\hat{\textbf{E}}}\dot{{\hat{\textbf{x}}}}(t)&={\hat{\textbf{A}}}{\textbf{x}}(t)+{\hat{\textbf{B}}}u(t),\
\ \hat{y}(t)={\hat{\textbf{C}}}{\hat{\textbf{x}}}(t).\end{aligned}\right.$
(10)

More implementation details and properties on the LF procedure can be found in
[4] and in [15].

##  3 The proposed hybrid method based on the Loewner framework and LS solves

The idea is to recover all the operators from measurements. The Loewner
framework is capable of recovering the linear part by providing a fitted
realization of dimension $r$ as in (9), denoted with
$({\hat{\textbf{E}}},{\hat{\textbf{A}}},{\hat{\textbf{B}}},{\hat{\textbf{C}}})$.
Now, based on these matrices, introduce
${\hat{\mathbf{\Phi}}}(s)=(s{\hat{\textbf{E}}}-{\hat{\textbf{A}}})^{-1}\in{\mathbb{C}}^{r\times
r}$.

The sampling domain is denoted with
$\Omega=\\{(\zeta_{1}^{(i)},\zeta_{2}^{(i)})\in{\mathbb{C}}^{2}|1\leq i\leq
K\\}$. For a particular pair of sampling points
$(\zeta_{1}^{(i)},\zeta_{2}^{(i)})\in{\mathbb{C}}^{2}$, we define the
following quantities:

$\displaystyle\begin{split}{\hat{\cal
O}}(\zeta_{1}^{(i)},\zeta_{2}^{(i)})&:={\hat{\textbf{C}}}{\hat{\mathbf{\Phi}}}(\zeta_{1}^{(i)}+\zeta_{2}^{(i)})\in{\mathbb{C}}^{1\times
r},\\\ {\hat{\cal
R}}_{q}(\zeta_{1}^{(i)},\zeta_{2}^{(i)})&:=({\hat{\mathbf{\Phi}}}(\zeta_{1}^{(i)}){\hat{\textbf{B}}}\otimes{\hat{\mathbf{\Phi}}}(\zeta_{2}^{(i)}){\hat{\textbf{B}}})\in{\mathbb{C}}^{r^{2}\times
1},\\\ {\hat{\cal
R}}_{b}(\zeta_{1}^{(i)},\zeta_{2}^{(i)})&:=({\hat{\mathbf{\Phi}}}(\zeta_{1}^{(i)}){\hat{\textbf{B}}}+{\hat{\mathbf{\Phi}}}(\zeta_{2}^{(i)}){\textbf{B}})\in{\mathbb{C}}^{r\times
1}.\end{split}$ (11)

It is to be noted that the vectors introduced in (11) are computed solely in
terms of the matrices
$({\hat{\textbf{E}}},{\hat{\textbf{A}}},{\hat{\textbf{B}}},{\hat{\textbf{C}}})$
corresponding to the data-driven Loewner surrogate model. Let
${\textbf{v}}\in{\mathbb{C}}^{K}$ be the vector of data measurements, i.e.,
containing samples of the second symmetric transfer function
$H_{2}(s_{1},s_{2})$ evaluated on the $\Omega$ grid. More precisely, let
${\textbf{v}}_{i}=H_{2}(\zeta_{1}^{(i)},\zeta_{2}^{(i)})$. Now, since we would
like to fit a reduced-order QB model to interpolate the 2D data, it follows
that
$H_{2}(\zeta_{1}^{(i)},\zeta_{2}^{(i)})={\hat{H}}_{2}(\zeta_{1}^{(i)},\zeta_{2}^{(i)})$.
We can write:

$\displaystyle\underbrace{{\hat{H}}_{2}(\zeta_{1}^{(i)},\zeta_{2}^{(i)})}_{{\textbf{v}}_{i}\in{\mathbb{C}}}=\underbrace{{\hat{\textbf{C}}}{\hat{\mathbf{\Phi}}}(s_{1}+s_{2})}_{{\hat{\cal
O}}(\zeta_{1}^{(i)},\zeta_{2}^{(i)})\in{\mathbb{C}}^{1\times
r}}{\hat{\textbf{Q}}}\underbrace{({\hat{\mathbf{\Phi}}}(\zeta_{1}^{(i)}){\hat{\textbf{B}}}\otimes{\hat{\mathbf{\Phi}}}(s_{2}){\hat{\textbf{B}}})}_{{\hat{\cal
R}}_{q}(\zeta_{1}^{(i)},\zeta_{2}^{(i)})\in{\mathbb{C}}^{r^{2}\times 1}}$ (12)
$\displaystyle+(1/2)\underbrace{{\textbf{C}}\boldsymbol{\Phi}(\zeta_{1}^{(i)}+\zeta_{2}^{(i)})}_{{\hat{\cal
O}}(\zeta_{1}^{(i)},\zeta_{2}^{(i)})\in{\mathbb{C}}^{1\times
r}}{\hat{\textbf{N}}}\underbrace{({\hat{\mathbf{\Phi}}}(\zeta_{1}^{(i)}){\hat{\textbf{B}}}+{\hat{\mathbf{\Phi}}}(\zeta_{2}^{(i)}){\textbf{B}})}_{{\hat{\cal
R}}_{b}(\zeta_{1}^{(i)},\zeta_{2}^{(i)})\in{\mathbb{C}}^{r\times 1}},$ (13)

and hence it follows that:

$\displaystyle{\hat{\textbf{H}}}_{2}(\zeta_{1}^{(i)},\zeta_{2}^{(i)})$
$\displaystyle={\hat{\cal
O}}(\zeta_{1}^{(i)},\zeta_{2}^{(i)}){\hat{\textbf{Q}}}{\hat{\cal
R}}_{q}(\zeta_{1}^{(i)},\zeta_{2}^{(i)})$ (14) $\displaystyle+(1/2){\hat{\cal
O}}(\zeta_{1}^{(i)},\zeta_{2}^{(i)}){\hat{\textbf{N}}}{\hat{\cal
R}}_{b}(\zeta_{1}^{(i)},\zeta_{2}^{(i)})$ (15)

###### Definition 3.1

Given a matrix ${\textbf{X}}\in{\mathbb{C}}^{m\times n}$, we denote with
$\text{vec}({\textbf{X}})$ the vector $(mn)\times 1$ computed as follows:

$\text{vec}({\textbf{X}})=\left[\begin{array}[]{ccc}{\textbf{X}}(1,:)&\cdots&{\textbf{X}}(m,:)\end{array}\right]^{T}\in{\mathbb{C}}^{mn},$
(16)

where the MATLAB notation ${\textbf{X}}(k,:)\in{\mathbb{C}}^{1\times n}$ was
used to refer to the kth row of X.

The vectorization procedure adapted to the data-driven problem shown here is
presented in (18). Let us denote with ${\cal
T}\in{\mathbb{C}}^{K\times(r^{3}+r^{2})}$, the matrix for which the ith row is
given by ${\cal T}(i,:)=\left[\begin{matrix}{\cal T}_{q}(i,:)&{\cal
T}_{b}(i,:)\end{matrix}\right]\in{\mathbb{C}}^{1\times(r^{3}+r^{2})}$

$\displaystyle\begin{cases}{\cal
T}_{q}(i,:)=\left[\begin{array}[]{c}{\hat{\cal
O}}(\zeta_{1}^{(i)},\zeta_{2}^{(i)})\otimes{\hat{\cal
R}}_{q}^{T}(\zeta_{1}^{(i)},\zeta_{2}^{(i)})\end{array}\right]\in{\mathbb{C}}^{K\times
r^{3}},\\\ {\cal T}_{b}(i,:)=\left[\begin{array}[]{c}{\hat{\cal
O}}(\zeta_{1}^{(i)},\zeta_{2}^{(i)})\otimes{\hat{\cal
R}}_{b}^{T}(\zeta_{1}^{(i)},\zeta_{2}^{(i)})\end{array}\right]\in{\mathbb{C}}^{K\times
r^{2}}.\end{cases}$ (17)

$\displaystyle{\textbf{v}}_{i}={\hat{H}}_{2}(\zeta_{1}^{(i)},\zeta_{2}^{(i)})$
$\displaystyle={\hat{\cal
O}}(\zeta_{1}^{(i)},\zeta_{2}^{(i)}){\hat{\textbf{Q}}}{\hat{\cal
R}}_{q}(\zeta_{1}^{(i)},\zeta_{2}^{(i)})+(1/2){\hat{\cal
O}}(\zeta_{1}^{(i)},\zeta_{2}^{(i)}){\hat{\textbf{N}}}{\hat{\cal
R}}_{b}(\zeta_{1}^{(i)},\zeta_{2}^{(i)})$ (18) $\displaystyle=\left({\hat{\cal
O}}(\zeta_{1}^{(i)},\zeta_{2}^{(i)})\otimes{\hat{\cal
R}}_{q}^{T}(\zeta_{1}^{(i)},\zeta_{2}^{(i)})\right)vec({\hat{\textbf{Q}}})+(1/2)\left({\hat{\cal
O}}(\zeta_{1}^{(i)},\zeta_{2}^{(i)})\otimes{\hat{\cal
R}}_{b}^{T}(\zeta_{1}^{(i)},\zeta_{2}^{(i)})\right)vec({\hat{\textbf{N}}})$
$\displaystyle=\underbrace{\left[\begin{array}[]{c|c}\underbrace{{\hat{\cal
O}}(\zeta_{1}^{(i)},\zeta_{2}^{(i)})\otimes{\hat{\cal
R}}_{q}^{T}(\zeta_{1}^{(i)},\zeta_{2}^{(i)})}_{={\cal
T}_{q}(i,:)}&\underbrace{{\hat{\cal
O}}(\zeta_{1}^{(i)},\zeta_{2}^{(i)})\otimes{\hat{\cal
R}}_{b}^{T}(\zeta_{1}^{(i)},\zeta_{2}^{(i)})}_{={\cal
T}_{b}(i,:)}\end{array}\right]}_{={\cal T}(i,:)\ \in\
{\mathbb{C}}^{1\times(r^{3}+r^{2})}}\cdot\underbrace{\left[\begin{array}[]{c}\text{vec}({\hat{\textbf{Q}}})\\\
\hline\cr(1/2)\text{vec}({\hat{\textbf{N}}})\end{array}\right]}_{={\textbf{z}}\
\in\ {\mathbb{C}}^{(r^{3}+r^{2})\times 1}}$

Now, from (18), by varying the index $i$ such as $1\leq i\leq K$, it follows
that we can put together a linear system of equations in $r^{3}+r^{2}$
unknowns as follows:

${\cal T}{\textbf{z}}={\textbf{v}},$ (19)

where
${\textbf{z}}=\left[\begin{matrix}\text{vec}({\hat{\textbf{Q}}})^{T}&(1/2)\text{vec}({\hat{\textbf{N}}})^{T}\end{matrix}\right]^{T}$
is the vector of variables which contains the entries of the vectorized
operators of the surrogate reduced-order QB system. In order to ensure an
over-determined linear system of equations, we clearly need to have enough
data measurements, i.e., the condition $K\geq r^{3}+r^{2}$ needs to hold true.
Then, we can employ a direct solution of system (19), e.g., by means of the
Moore-Penrose pseudo-inverse or by using Gaussian elimination. If the matrix
${\cal T}$ isindeed of full column rank, then it does not matter what
procedure is chosen (for the direct solve). However, in most cases, the matrix
${\cal T}$ is not of full column rank, and hence direct solves need to be
carefully dealt with (by introducing regularization techniques). In what
follows, we will use a truncated singular value decomposition (tSVD) approach.
This is an attractive and powerful method since it uses the optimal rank-k
approximation of the SVD (in the 2 norm). Such approach has been already used
for applying OpInf [5], together with the Tikhonov regularization scheme [20]
and the tQR (truncated QR decomposition) approach.

###### Remark 3.1

Finding a viable solution for the least-squares problem stated in (19) is not
a straightforward task and can be computationally challenging task because of
the ill-conditioning of matrix ${\cal T}$ (the column rank of matrix ${\cal
T}$ is sometimes much smaller than $r^{3}+r^{2}$). Another challenge is the
computational cost which grows cubically in the order $r$ of the reduced-order
system.

##  4 A simple model of a non-linear circuit

We consider a simple circuit constructed from two blocks connected in series.
Each block contains a capacitor in parallel with a diode. This is modeled as a
dynamical system in two variables given by the voltage drops on each block.
The input is given by the current through the circuit, while the observed
output is the sum of the two variables. More precisely, we can write the
differential equations characterizing the dynamics, as follows:

$\left\\{\begin{aligned}
C_{1}\frac{dV_{1}(t)}{dt}&=I(t)-I_{r_{1}}\left(e^{\frac{1}{V_{t_{1}}}V_{1}(t)}-1\right),\\\
C_{2}\frac{dV_{2}(t)}{dt}&=I(t)-I_{r_{2}}\left(e^{\frac{1}{V_{t_{2}}}V_{2}(t)}-1\right),\end{aligned}\right.$
(20)

while the output is $y(t)=V_{1}(t)+V_{2}(t)$. The capacitance values are
denoted with $C_{i},1\leq i\leq 2$, while other constants are denoted with
$I_{r_{i}}$ and $V_{t_{i}}$. Next, we introduce
$x_{1}(t)=V_{1}(t)/V_{t_{1}},~{}x_{2}(t)=V_{2}(t)/V_{t_{2}}$, and let
$a=\frac{1}{C_{1}V_{t_{1}}},~{}b=\frac{1}{C_{2}V_{t_{2}}},~{}c=I_{r_{1}},~{}d=I_{r_{2}}$.
Using all of these, the nonlinear system in (20) is rewritten as:

$\left\\{\begin{aligned}
\dot{x}_{1}(t)&=aI(t)-ac\left(e^{x_{1}(t)}-1\right),\\\
\dot{x}_{2}(t)&=bI(t)-bd\left(e^{x_{2}(t)}-1\right),\\\
y(t)&=x_{1}(t)V_{t_{1}}+x_{2}(t)V_{t_{2}}.\end{aligned}\right.$ (21)

###  4.1 Analysis based on Taylor series truncation

The methods discussed here are inexact, i.e., based on approximation (e.g., on
truncating the Taylor series associated to the nonlinearity). First, by using
a truncated Taylor series given by the formula $e^{x_{i}(t)}\approx
1+x_{i}(t)$, and substitute it in (21), it follows that the following linear
dynamical system is obtained:

$\left\\{\begin{aligned} \dot{x}_{1}(t)&=aI(t)-acx_{1}(t),\\\
\dot{x}_{2}(t)&=bI(t)-bdx_{2}(t),\\\
y(t)&=x_{1}(t)V_{t_{1}}+x_{2}(t)V_{t_{2}}+R\cdot I(t).\end{aligned}\right.$
(22)

Then, identify the following linear realization:

$\displaystyle{\textbf{A}}_{L}=\left[\begin{matrix}-ac&0\\\
0&-bd\end{matrix}\right],\ \ {\textbf{B}}_{L}=\left[\begin{matrix}a\\\
b\end{matrix}\right],\ \ \
{\textbf{C}}_{L}=\left[\begin{matrix}V_{t_{1}}&V_{t_{2}}\end{matrix}\right].$

Next, we discuss Carleman’s linearization proposed in [9]. It is a method used
to embed a nonlinear system of differential equations of finite dimension into
a system of bilinear differential equations of infinite dimension. By
truncating the hence obtained bilinear system at finite orders, one finds a
systematic way of achieving an arbitrary-order approximation. By following
this procedure, introduce the augmented new state by stopping at quadratic
terms (and ignoring all other higher powers), as

${\textbf{x}}^{C}=\left[\begin{matrix}{\textbf{x}}\\\
{\textbf{x}}\otimes{\textbf{x}}\end{matrix}\right]=\left[\begin{matrix}x_{1}&x_{2}&x_{1}^{2}&x_{1}x_{2}&x_{2}x_{1}&x_{2}^{2}\end{matrix}\right]^{T},$
(23)

and appropriately compute the derivatives of the last four entries as, for
example, of:

$\displaystyle\begin{split}\dot{x}^{C}_{3}(t)&=\frac{d}{dt}x_{1}^{2}(t)=2x_{1}(t)\dot{x}_{1}(t)\\\
&=2a{\textbf{x}}_{1}(t)I(t)-2acx_{1}^{2}(t)-acx_{1}^{3}(t).\end{split}$ (24)

Finally, the approximation steps follows, i.e., we neglect the powers in
$x_{i}(t)$ higher than 3, and write:

$\displaystyle\begin{split}\dot{x}^{C}_{3}(t)&\approx
2ax_{1}(t)I(t)-2acx_{1}^{2}(t)\\\ &\Rightarrow\dot{x}^{C}_{3}(t)\approx
2ax^{C}_{1}(t)I(t)-2acx_{2}^{C}(t),\end{split}$ (25)

which includes a linear term, i.e., $x_{2}^{C}(t)$, and a bilinear term:
$x^{C}_{1}(t)I(t)$. The same procedure is applied for the other entries of the
derivative of the new state vector $\dot{{\textbf{x}}}^{C}(t)$. Hence, an
approximate bilinear systems is derived:

$\left\\{\begin{aligned}
\dot{x}^{C}(t)&={\textbf{A}}{\textbf{x}}^{C}(t)+{\textbf{N}}{\textbf{x}}^{C}(t)u(t)+{\textbf{B}}u(t),\
\ y(t)={\textbf{C}}{\textbf{x}}^{C}(t).\end{aligned}\right.$

###  4.2 Polynomial lifting

The first step towards implementing a ”lifting approach” is to identify the
quantities that depend non-linearly on the original states, and to consider
them as new artificial states. For example, we introduce the auxiliary
variables

$x_{3}(t):=e^{x_{1}(t)}-1,\ \ \ \ x_{4}(t):=e^{x_{2}(t)}-1.$ (26)

In this way, we also enforce zero initial conditions. The augmented system is
hence written as:

$\hskip 5.69054pt\left\\{\begin{aligned} \dot{x}_{1}(t)&=aI(t)-acx_{3}(t),\ \
\dot{x}_{2}(t)=bI(t)-bdx_{4}(t),\\\
\dot{x}_{3}(t)&=-acx_{3}(t)-acx_{3}^{2}(t)+aI(t)+ax_{3}(t)I(t)\\\
\dot{x}_{4}(t)&=-bdx_{4}(t)-bdx_{4}^{2}(t)+bI(t)+bx_{4}(t)I(t)\\\
y(t)&=x_{1}(t)V_{t_{1}}+x_{2}(t)V_{t_{2}}+R\cdot I(t).\end{aligned}\right.$
(27)

The system in (27) can be written in the general QB form:

$\left\\{\begin{aligned}
\dot{{\textbf{x}}}(t)&={\textbf{A}}{\textbf{x}}(t)+{\textbf{Q}}({\textbf{x}}(t)\otimes{\textbf{x}}(t))+{\textbf{N}}{\textbf{x}}(t)I(t)+{\textbf{B}}I(t),\\\
y(t)&={\textbf{C}}{\textbf{x}}(t),\end{aligned}\right.\\\ $

where the matrices are defined next:

A $\displaystyle=\left[\begin{array}[]{cccc}0&0&-ac&0\\\ 0&0&0&-bd\\\
0&0&-ac&0\\\
0&0&0&-bd\end{array}\right],{\textbf{N}}=\left[\begin{array}[]{cccc}0&0&0&0\\\
0&0&0&0\\\ 0&0&a&0\\\
0&0&0&b\end{array}\right],~{}{\textbf{B}}=\left[\begin{array}[]{c}a\\\ b\\\
a\\\ b\end{array}\right],$ C
$\displaystyle=\left[\begin{array}[]{cccc}V_{t_{1}}\\\ V_{t_{2}}\\\ 0\\\
0\end{array}\right]^{T},{\textbf{Q}}=\left[\begin{array}[]{c|cccc}\textbf{0}_{2\times
12}&\textbf{0}_{2\times 4}\\\ \hline\cr\textbf{0}_{2\times
12}&{\textbf{Q}}_{22}\\\
\end{array}\right],~{}{\textbf{Q}}_{22}=\left[\begin{array}[]{cccc}-ac&0&0&0\\\
0&0&0&-bd\end{array}\right]$

The first two symmetric transfer functions, where
$\boldsymbol{\Phi}(s)=(s{\textbf{E}}-{\textbf{A}})^{-1}\in\mathbb{C}^{n\times
n}$, can be written as:

$\displaystyle\begin{split}H_{1}(s)&={\textbf{C}}\boldsymbol{\Phi}(s){\textbf{B}}=\frac{V_{t_{1}}}{C_{1}V_{t_{1}}s+I_{r_{1}}}+\frac{V_{t_{2}}}{C_{2}V_{t_{2}}s+I_{r_{2}}},\\\
H_{2}(s,s)&={\textbf{C}}\boldsymbol{\Phi}(2s){\textbf{Q}}\left[\boldsymbol{\Phi}(s){\textbf{B}}\otimes\boldsymbol{\Phi}(s){\textbf{B}}\right]+{\textbf{C}}\boldsymbol{\Phi}(2s){\textbf{N}}\boldsymbol{\Phi}(s){\textbf{B}}\\\
&=-\frac{{I}_{r_{1}}{V_{t_{1}}}}{2\,{\left({I}_{r_{1}}+C_{1}\,{V_{t_{1}}}\,s\right)}^{2}\,\left({I}_{r_{1}}+2\,C_{1}\,V_{t_{1}}\,s\right)}\\\
&-\frac{{I}_{r_{2}}V_{t_{2}}}{2\,\left({I}_{r_{2}}+C_{2}\,V_{t_{2}}\,s\right)^{2}\,\left({I}_{r_{2}}+2\,C_{2}\,V_{t_{2}}\,s\right)}\end{split}$

The parameters can be recovered from LF by considering the 2nd kernel as a
univariate rational function, but the pole residue form needs a special
treatment and is usually quite challenging. Similar studies have been proposed
in [17, 12] for inferring bilinear or quadratic systems respectively, where an
improvement towards bilinear identification was shown in [16] and in the
current study.

##  5 Numerical experiments: a nonlinear RC Ladder circuit

We analyze a nonlinear RC-ladder electronic circuit first introduced in [10].
Various variants of this model were also mentioned in other MOR works,i.e.,
[13] and [8]. This nonlinear first-order system models a resistor-capacitor
network that exhibits a nonlinear behaviour caused by the nonlinear resistors
consisting of a parallel connected resistor with a diode.

As presented in [10], the underlying model is given by a SISO system of the
form:

$\displaystyle\dot{x}(t)=\left[\begin{matrix}-g(x_{1}(t))-g(x_{1}(t)-x_{2}(t))\\\
g(x_{1}(t)-x_{2}(t))-g(x_{2}(t)-x_{3}(t))\\\ \vdots\\\
g(x_{k-1}(t)-x_{k}(t))-g(x_{k}(t)-x_{x+1}(t))\\\ \vdots\\\
g(x_{N-1}(t)-x_{N}(t))\end{matrix}\right]+\left[\begin{matrix}u(t)\\\ 0\\\
\vdots\\\ 0\\\ \vdots\\\ 0\end{matrix}\right],$ (28)

with $y(t)=x_{1}(t)$, where the mapping $g$ is given by
$g:\mathbb{R}\to\mathbb{R}$ defined as $g(x_{i})=g_{D}(x_{i})+x_{i}$, which
combines the effect of a diode and a resistor. The non-linearity $g_{D}$
models a diode as a nonlinear resistor, based on the classical Shockley model:

$g_{D}(x_{i})=i_{S}(\exp(u_{P}x_{i})-1),$ (29)

with material parameters $i_{S}>0$ and $u_{P}>0$. For this benchmark, the
parameters are selected as follows: $i_{S}$ = 1 and $u_{P}=40$ as in [10]. By
substituting these values into (29), we get that
$g_{D}(x_{i})=\exp(40x_{i})-1$, and hence it follows that
$g(x_{i})=\exp(40x_{i})+x_{i}-1$. In what follows, we will apply the proposed
methods for:

###  5.1 Bilinear treatment via Carleman’s approach

The original nonlinear system is transformed into a bilinear system by means
of Carleman linearization, as originally shown in [10] and later in [8]. The
matrices are:

$\displaystyle{\textbf{A}}=\left[\begin{array}[]{cc}{\textbf{A}}_{1}&\frac{1}{2}{\textbf{A}}_{2}\\\
\textbf{0}&{\textbf{A}}_{1}\otimes{\textbf{I}}+{\textbf{I}}\otimes{\textbf{A}}_{1}\end{array}\right],~{}{\textbf{x}}=\left[\begin{array}[]{c}{\textbf{v}}\\\
{\textbf{v}}\otimes{\textbf{v}}\end{array}\right],$ (30)
$\displaystyle{\textbf{N}}=\left[\begin{array}[]{cc}\textbf{0}&\textbf{0}\\\
{\textbf{B}}\otimes{\textbf{I}}+{\textbf{I}}\otimes{\textbf{B}}&\textbf{0}\end{array}\right],~{}{\textbf{B}}=\left[\begin{array}[]{c}{\textbf{b}}\\\
\textbf{0}\end{array}\right],~{}{\textbf{C}}=\left[\begin{array}[]{c}{\textbf{c}}\\\
\textbf{0}\end{array}\right]^{T}.$

Consequently, the resulting bilinear system has dimension $n^{2}+n$ with $n$
as the number of circuit blocks of the original system. More details on the
structure of the involved matrices in Eq. (30) can be found in [8].

###  5.2 Lifting to Quadratic-Bilinear form

Analogues with the example in Sec. 4.2, the original RC-ladder model can be
lifted to an equivalent quadratic bilinear model. The introduced additional
state variables ${\textbf{x}}_{1}=\upsilon_{1}$ and
$x_{i}=\upsilon_{i}-\upsilon_{i+1}$ followed by introducing the additional
state variables $z_{1}=e^{-40\upsilon_{1}}-1$ and $z_{i}=e^{40x_{i}}$ can
transform equivalently the original system Eq. (28) to a quadratic-bilinear
form Eq. (3) with dimension $2n$ [8].

In Fig. 1 the original nonlinear system along with the equivalent quadratic-
bilinear and the approximated bilinear are depicted. The numerical difference
between the original and the quadratic-bilinear has reached machine precision
where the bilinear for this amplitude starts to differ significantly.

Figure 1: The original RC-ladder model with $n=50$ circuit blocks. The lifted
QB model is equivalent with the original nonlinear. The bilinear model of
dimension $n^{2}+n$ offers good approximation only for relative small input
amplitudes.

The aim is to achieve MOR by measuring the first symmetric Volterra kernels
from input-output time-domain simulations. As we infer the operators from the
$2$nd kernel a double tone input is considered. The scheme for kernel
separation and harmonic indexing remains the same as in [16]. By simulating
the original nonlinear model in the time domain under the excitation of a
double-tone harmonic input, an accurate separation of kernels can be achieved
(the Fourier transform is indeed accurate). Here, as we want to illustrate the
efficacy of the proposed method by inferring the operators from the $2$nd
Kernel, we assume a perfect measurement setup.

Step 1: The first harmonic can be measured under excitation with a single tone
input. Then, measurements of the first kernel $H_{1}(s_{1})$ (e.g., the
magnitude and phase) can be derived. The Loewner framework constructs a low-
order rational interpolant and identifies the minimal linear sub-system of
order $r$. In Fig. 2 and on the left pane, the Loewner singular value decay
offers the criterion for reduction. The order $r=10$ was chosen, as the $11$th
singular value is close to machine precision. Hence, it is an indicator for
recovering the original linear dynamics. Denote with
$\Sigma_{\texttt{lin}}:(\hat{{\textbf{A}}},\hat{{\textbf{B}}},\hat{{\textbf{C}}})$
the linear reduced system of order $r=10$.

In Fig. 2 and on the right pane, the approximation results are depicted.

Figure 2: Left pane: The Loewner singular value decay (used to decide on the
reduced order). Right pane: The reduced-order Loewner interpolant computed
from the $1$st kernel that can reach machine precision approximation.

Step 2: By having access to the reduced operators of the linear system
$\hat{\Sigma}_{\texttt{lin}}$, we can infer the remaining nonlinear quadratic
and bilinear operators from Eq. (18). The second harmonic can be measured with
a double tone excitation. Thus, measurements of the second kernel
$H_{2}(s_{1},s_{2})$ over the whole complex domain of definition can be
collected repetitively. It is important here to mention that the amount of
measurements is related to the reduced dimension $r$, where at least
$K\geq(r^{3}+r^{2})$ measurements ensures enough data, for the solution of Eq.
(19). By enforcing quadratic symmetries (the matrix $\hat{{\textbf{Q}}}$ is
set to satisfy the property
$\hat{{\textbf{Q}}}({\textbf{w}}\otimes{\textbf{v}})=\hat{{\textbf{Q}}}({\textbf{v}}\otimes{\textbf{w}})$),
the complexity can be further reduced. Using some algebraic adjustments,
another simplification can be performed by replacing the symmetric Kronecker
product with the asymmetric one ($\otimes^{\prime}$ as in [6]). Solving for
the vector z yields a reduced-order QB system
$\hat{\Sigma}_{QB}=(\hat{{\textbf{A}}},\hat{{\textbf{Q}}},\hat{{\textbf{N}}},\hat{{\textbf{B}}},\hat{{\textbf{C}}})$.

Figure 3: Left pane: The original 2nd kernel along with the measurements
(green dots). Right pane: The 2nd kernel of the reduced system.
$\|H_{2}(s_{1},s_{2})-\hat{H}_{2}(s_{1},s_{2})\|_{\infty}\sim 10^{-8}$.

Finally, in Fig. 4 the time domain solution is depicted and is compared to the
original nonlinear response.

Figure 4: Left pane: The original nonlinear RC-ladder along with the
approximant that is constructed with the new LoewQB method. Right pane: The
absolute error over a dense grid.

The proposed method successfully constructs a reduced quadratic-bilinear model
from input-output time domain data that accurately approximates the response
of the original nonlinear system. The simulations are performed with the
multi-step Runge Kutta scheme (e.g., using ODE45 in Matlab) with input
$u(t)=0.01e^{-t}$, and the maximum error is $\|(y(t)-y_{r}(t))\|_{\infty}\sim
10^{-6}$.

##  6 Conclusion and future endeavours

A non-intrusive data-driven method that constructs nonlinear models with the
quadratic-bilinear structure from input-output time-domain data was presented.
The new method was based on the Loewner and Volterra frameworks. The second
symmetric kernel can be inferred as the $2$nd harmonic of the transformed
output and hence, the problem of estimating the nonlinear operators can be
resolved into a linear LS system. Although this LS system may be under-
determined, it contains enough information for estimating the nonlinear
operators with the Moore-Penrose pseudo-inverse (truncated SVD), thus
providing good approximations. The use of higher harmonics (e.g., $3$rd
kernel, etc.) that results in a nonlinear optimization problem along with the
use of other regularization techniques as the one proposed by Tikhonov, will
be the topics of future research endeavors in connection also with the
analysis on the sensitivity of noise that real data contain.

## References

  * [1] A. C. Antoulas, C. A. Beattie, and S. Gugercin.  Interpolatory Methods for Model Reduction.  SIAM, Philadelphia, 2020.
  * [2] A. C. Antoulas, I. V. Gosea, and M. Heinkenschloss.  On the Loewner framework for model reduction of Burgers’ equation.  In R. King, editor, Active Flow and Combustion Control, pages 255–270. Springer, Cham, Switzerland, 2019.
  * [3] A. C. Antoulas, I. V. Gosea, and A. C. Ionita.  Model reduction of bilinear systems in the Loewner framework.  SIAM Journal on Scientific Computing, 38(5):B889–B916, 2016.
  * [4] A. C. Antoulas, S. Lefteriu, and A. C. Ionita.  A tutorial introduction to the Loewner framework for model reduction.  In Model Reduction and Approximation, chapter 8, pages 335–376. SIAM, 2017.
  * [5] P. Benner, P. Goyal, J. Heiland, and I. Pontes Duff.  Operator inference and physics-informed learning of low-dimensional models for incompressible flows.  e-prints 2010.06701, arXiv, 2020.  math.DS.
  * [6] P. Benner, P. Goyal, B. Kramer, B. Peherstorfer, and K. Willcox.  Operator inference for non-intrusive model reduction of systems with non-polynomial nonlinear terms.  Comp. Methods in App. Mechanics and Engineering, 372, 2020.
  * [7] P. Benner, M. Ohlberger, A. Cohen, and K. Willcox, editors.  Model Reduction and Approximation: Theory and Algorithms.  Computational Science & Engineering. Society for Industrial and Applied Mathematics, Philadelphia, PA, 2017.
  * [8] T. Breiten.  Interpolatory Methods for Model Reduction of Large-Scale Dynamical Systems.  Dissertation, Department of Mathematics, Otto-von-Guericke University, Magdeburg, Germany, 2013.
  * [9] T. Carleman.  Application de la théories des équations intégrales linéaires aux systèmes d’équations différentielles non linéaires.  Acta. Math., 59:63–87, 1932.
  * [10] Y. Chen.  Model reduction for nonlinear systems.  Master’s thesis, M.I.T, Boston, USA, 1999.
  * [11] 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.
  * [12] I. V. Gosea, D. S. Karachalios, and A. C. Antoulas.  Learning reduced-order models of quadratic control systems from input-output data.  e-print 2012.02075, arXiv, 2020.  math.OC, presented in the European Control Conference, June 29 - July 2, 2021, Rotterdam, The Netherlands.
  * [13] C. Gu.  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.
  * [14] B. L. Ho and R. E. Kalman.  Effective construction of linear state-variable models from input/output functions.  at - Automatisierungstechnik, 14(1-12):545–548, 1966.
  * [15] D. S. Karachalios, I. V. Gosea, and A. C. Antoulas.  The Loewner framework for system identification and reduction.  In P. Benner, S. Grivet-Talocia, A. Quarteroni, G. Rozza, W. H. A. Schilders, and L. M. Silveira, editors, Model Reduction Handbook: Volume I: System- and Data-Driven Methods and Algorithms, volume 1 of Handbook on Model Reduction. De Gruyter, 2020.  accepted December 2019.
  * [16] D. S. Karachalios, I. V. Gosea, and A. C. Antoulas.  On bilinear time-domain identification and reduction in the Loewner framework.  In Model Reduction of Complex Dynamical Systems, volume 171 of International Series of Numerical Mathematics, pages 3–30. Birkhäuser, Cham, 2021.
  * [17] Dimitrios S. Karachalios, Ion Victor Gosea, and Athanasios C. Antoulas.  A bilinear identification-modeling framework from time domain data.  PAMM, 19(1):e201900246, 2019.
  * [18] B. Kramer and K. E. Willcox.  Nonlinear model order reduction via lifting transformations and proper orthogonal decomposition.  AIAA Journal, 57(6):2297–2307, 2019.
  * [19] L. Ljung.  System Identification: Theory for the User.  Prentice Hall, Upper Saddle River, New Jersey, 1987 (2nd edition in 1999), 1999.
  * [20] 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.
  * [21] P. van Overschee and B L. de Moor.  Subspace Identification for Linear Systems: Theory — Implementation — Applications.  Springer, 1996.