# Approximating a flexible beam model in the Loewner framework

Alexander Zuyev<EMAIL_ADDRESS>Ion Victor Gosea  gosea@mpi-
magdeburg.mpg.de

###### Abstract

The paper develops the Loewner approach for data-based modeling of a linear
distributed-parameter system. This approach is applied to a controlled
flexible beam model coupled with a spring-mass system. The original dynamical
system is described by the Euler-Bernoulli partial differential equation with
the interface conditions due to the oscillations of the lumped part. The
transfer function of this model is computed analytically, and its sampled
values are then used for the data-driven design of a reduced model. A family
of approximate realizations of the corresponding input-output map is
constructed within the Loewner framework. It is shown that the proposed
finite-dimensional approximations are able to capture the key properties of
the original dynamics over a given range of observed frequencies. The
robustness of the method to noisy data is also investigated.

###### keywords:

Data-driven modeling, flexible structure, Euler–Bernoulli beam, distributed
parameter systems, Loewner framework, linear systems, model reduction, noisy
data, input-output map, transfer function.

##  1 Introduction

Model reduction techniques can be employed to replace a large-scale system
with a complex structure (characterized by multidimensional systems of
ordinary differential equations and/or partial differential equations), with a
much simpler and smaller dynamical system (characterized by few equations with
well-understood dynamics). In the last decades, there have been many
methodologies proposed in this direction; we refer the reader to [2, 32, 24,
9, 3] for more details.

A viable alternative to using classical model reduction approaches based on
single or double-sided projections (that usually require explicit access to a
large-scale model) is to use instead data-driven methods. These latter do not
require explicit access to the large-scale model’s structure or matrices. We
mention here the Loewner framework (LF) [27], Vector fitting (VF) [19], or the
AAA algorithm [28]. When using these, low-order models can be constructed
directly from data in the frequency domain (samples of the transfer function).
Such methods can be viewed as rational approximation tools by means of
interpolation (LF), least-squares fit (VF), or a mixed approach (AAA). Other
data-driven methods that has emerged in recent years are dynamic mode
decomposition (DMD) and operator inference (OpInf), which use time-domain
snapshots of the state variables and then fit a particular structured model by
computing the appropriate matrices (in reduced coordinates). Details on DMD
can be found in [24], while details on OpInf can be found in [29, 6].

We consider here the problem of data-driven rational approximation by means of
fitting a linear time-invariant (LTI) dynamical system to a set of
measurements (in the frequency domain). The fitted LTI system is characterized
in the state-space by the following equations:

$\begin{cases}{\mathbf{E}}\dot{{\mathbf{x}}}(t)={\mathbf{A}}{\mathbf{x}}(t)+{\mathbf{B}}u(t),\\\
y(t)\ \hskip
4.2679pt={\mathbf{C}}{\mathbf{x}}(t)+{\mathbf{D}}u(t),\end{cases}$ (1)

where $u(t)\in\mathbb{R}$ is the input, $y(t)\in\mathbb{R}$ is the output,
${\mathbf{x}}(t)\in\mathbb{R}^{n}$ is the state vector, and the system
matrices are ${\mathbf{A}},{\mathbf{E}}\in{\mathbb{R}}^{n\times
n},\leavevmode\nobreak\ {\mathbf{B}},{\mathbf{C}}^{T}\in{\mathbb{R}}^{n\times
1}$. The transfer function of (1) is given by
$H(s)={\mathbf{C}}(s{\mathbf{E}}-{\mathbf{A}})^{-1}{\mathbf{B}}+{\mathbf{D}}$.
We refer to [2] for more details on various methodologies especially tailored
to the reduction of linear systems. In recent years, such methods for linear
systems have been steadily extended to particular classes of structured linear
systems [5, 30], or even to nonlinear structured systems in [16] (without
preservation of structure) or in [7] (with preservation of structure). The
structures treated include distributed parameters, delay terms or integro-
differential equations. In a data-driven setup, the structure-preserving
approach is [33] was proposed. Note that the analytical representations of
transfer functions have been obtained only for particular classes of
distributed parameter systems. We refer to [11, 1] for surveys of results in
this area. In the former tutorial article, the authors provide the derivation
of a variety of (irrational) transfer functions for systems described by
partial-differential equations. It is also shown that the choice of boundary
conditions have an influence on the dynamics and on the locations of poles and
zeros. In most practical situations, it is desirable to approximate the
irrational transfer function by a rational one, for the purpose of controller
design.

The Loewner framework approach was shown to be extremely powerful in data-
based control problems for wide classes of finite-dimensional control systems,
whose transfer functions are rational [4, 3]. However, the efficiency of the
Loewner framework for distributed-parameter systems (characterized by
irrational transfer functions) still remains to be verified, and the present
paper aims at filling this gap. Preliminary analysis was provided for linear
time-delay systems in [34, 26], fractional-order systems in [10], or for
control purposes in [17, 31]. A recent overview was provided in [22],
including amongst others, rational approximation of the Bessel function, of a
hyperbolic sine, and of a vibrating beam model from [11]. Hence, the LF was
studied in the context of approximating infinite-dimensional models of
vibrating beams (with finite-dimensional ones). However, as far as the authors
are aware, this is the first contribution that also takes into account the
effects of perturbed data, i.e., under additive Gaussian noise.

##  2 Vibrating beam with attached mass

Consider the Euler–Bernoulli equation describing the transverse vibrations of
a flexible beam of length $l$:

$\ddot{w}(x,t)+\frac{EI}{\rho}w^{\prime\prime\prime\prime}(x,t)+d\,\dot{w}^{\prime\prime\prime\prime}(x,t)=\frac{1}{\rho}\sum\limits_{j=1}^{k}\psi_{j}^{\prime\prime}(x)u_{j},$
(2)

where $w(x,t)$ is the beam deflection at point $x\in[0,l]$ and time $t$, $E$
is the Young’s modulus, $I$ is the area moment of inertia of the cross-
section, $\rho$ is the mass per unit length of the beam, and $d\geq 0$ is the
structural damping coefficient. We denote the derivative with respect to time
by a dot, while the prime denotes the spatial derivative (i.e., with respect
to $x$). We assume that a mass-spring system (shaker) is attached to the beam
at point $x=l_{0}$, so that equation (2) holds for $x\in[0,l_{0})$ and
$x\in(l_{0},l]$, and the interface condition is imposed at $x=l_{0}$:

$(m\ddot{w}+\varkappa
w)\Big{|}_{x=l_{0}}=(EIw^{\prime\prime})^{\prime}\Big{|}_{x=l_{0}-0}-(EIw^{\prime\prime})^{\prime}\Big{|}_{x=l_{0}+0}+u_{0}.$
(3)

The beam is hinged at both ends, which is formalized by the boundary
conditions

$w\Big{|}_{x=0}=w\Big{|}_{x=l}=0,\quad
w^{\prime\prime}\Big{|}_{x=0}=w^{\prime\prime}\Big{|}_{x=l}=0.$ (4)

System (2)–(4) is controlled by the force $u_{0}$ applied to the shaker at
$x=l_{0}$ and $k$ piezo actuators, whose actions $u_{j}$ are characterized in
terms of shape functions $\psi_{j}(x)$, $j=1,...,k$. It is also assumed that
$p$ piezo sensors are located at the points $x=l_{i}$, $i=1,..,p$, i.e. the
system outputs are

$y_{i}(t)=\left.\frac{\partial^{2}w(x,t)}{\partial
x^{2}}\right|_{x=l_{i}},\quad i=1,...,p.$ (5)

The above mathematical model has been presented in [21, 20] for the case
without damping; here we take into account the structural damping by
introducing the parameter $d$ in (2).

##  3 Computation of the transfer function

Let $u_{0}(t)$, $u_{1}(t)$, …, $u_{k}(t)$ $(t\geq 0)$ be inputs of the control
system (2)–(5) with zero initial data, denote the Laplace transform of the
inputs and outputs by

$U_{j}(s)=\int_{0}^{+\infty}u_{j}(t)e^{-st}dt,\quad j=0,1,...,k,$

and

$Y_{i}(s)=\int_{0}^{+\infty}y_{i}(t)e^{-st}dt,\quad i=1,...,p,$

respectively.

After introducing the Laplace transform of $w(x,t)$ with respect to $t$:
$W_{1}^{s}(x)=\int_{0}^{+\infty}w(x,t)e^{-st}dt$, we obtain from (2) the
following system of ordinary differential equations for
$W^{s}(x)=\left(W_{1}^{s},W_{2}^{s},W_{3}^{s},W_{4}^{s}\right)^{T}$:

$\frac{d}{dx}W^{s}(x)=AW^{s}(x)+\Phi^{s}(x),A=\left(\begin{array}[]{cccc}0&1&0&0\\\
0&0&1&0\\\ 0&0&0&1\\\ -4\gamma^{4}&0&0&0\\\ \end{array}\right),$ (6)

with

$\Phi^{s}(x)=\left(\begin{array}[]{c}0\\\ 0\\\ 0\\\
\phi^{s}(x)\end{array}\right),\;\gamma=\frac{\alpha\sqrt{2s}}{2},$ (7)
$\alpha^{4}=\frac{\rho}{EI+\rho ds}>0,\,\phi^{s}(x)=\frac{1}{EI+\rho
ds}\sum_{j=1}^{k}\psi_{j}(x)U_{j}(s),$ (8)

and the $s$ variable is treated as a parameter in (6). The general solution of
(6) is represented with the matrix exponential as follows:

$W^{s}(x)=\left\\{\begin{array}[]{ll}e^{xA}\bar{W}^{0}+\int_{0}^{x}e^{(x-y)A}\Phi^{s}(y)dy,&x\in[0,l_{0}],\\\
e^{(x-l)A}\bar{W}^{l}-\int_{x}^{l}e^{(x-y)A}\Phi^{s}(y)dy,&x\in(l_{0},l],\end{array}\right.$
(9)

where

$e^{xA}=\left(\begin{array}[]{cccc}z_{1}(x)&z_{2}(x)&z_{3}(x)&z_{4}(x)\\\
-4\gamma^{4}z_{4}(x)&z_{1}(x)&z_{2}(x)&z_{3}(x)\\\
-4\gamma^{4}z_{3}(x)&-4\gamma^{4}z_{4}(x)&z_{1}(x)&z_{2}(x)\\\
-4\gamma^{4}z_{2}(x)&-4\gamma^{4}z_{3}(x)&-4\gamma^{4}z_{4}(x)&z_{1}(x)\\\
\end{array}\right),$ $\displaystyle z_{1}(x)$ $\displaystyle=\cosh(\gamma
x)\cos(\gamma x),$ (10) $\displaystyle z_{2}(x)$
$\displaystyle=\frac{\cosh(\gamma x)\sin(\gamma x)+\sinh(\gamma x)\cos(\gamma
x)}{2\gamma},$ $\displaystyle z_{3}(x)$ $\displaystyle=\frac{\sinh(\gamma
x)\sin(\gamma x)}{2\gamma^{2}},$ $\displaystyle z_{4}(x)$
$\displaystyle=\frac{\cosh(\gamma x)\sin(\gamma x)-\sinh(\gamma x)\cos(\gamma
x)}{4\gamma^{3}}.$

Formula (9) represents the solutions of (6) in terms of their boundary values
$\bar{W}^{0}$ and $\bar{W}^{l}$ at $x=0$ and $x=l$, respectively, From the
boundary conditions (4), we conclude that

$\bar{W}^{0}=(0,\bar{W}^{0}_{2},0,\bar{W}^{0}_{4})^{T},\;\bar{W}^{l}=(0,\bar{W}^{l}_{2},0,\bar{W}^{l}_{4})^{T}.$
(11)

To eliminate the parameters $\bar{W}^{0}_{2}$, $\bar{W}^{0}_{4}$,
$\bar{W}^{l}_{2}$, $\bar{W}^{l}_{4}$, we exploit the property that $W^{s}(x)$
is of class $C^{2}[0,l]$ together with the interface condition (3). As a
result, we get the following linear algebraic system with respect to
$\bar{W}=(\bar{W}^{0}_{2},\bar{W}^{0}_{4},\bar{W}^{l}_{2},\bar{W}^{l}_{4})^{T}$:

$M\bar{W}=R,$ (12)

where

$\displaystyle{\scriptsize
M=\begin{pmatrix}-z_{2}(l_{0})&-z_{4}(l_{0})&z_{2}(l_{0}-l)&z_{4}(l_{0}-l)\\\
-z_{1}(l_{0})&-z_{3}(l_{0})&z_{1}(l_{0}-l)&z_{3}(l_{0}-l)\\\
4\gamma^{4}z_{4}(l_{0})&-z_{2}(l_{0})&-4\gamma^{4}z_{4}(l_{0}-l)&z_{2}(l_{0}-l)\\\
\beta z_{2}(l_{0})+4\gamma^{4}z_{3}(l_{0})&\beta
z_{4}(l_{0})-z_{1}(l_{0})&-4\gamma^{4}z_{3}(l_{0}-l)&z_{1}(l_{0}-l)\end{pmatrix},R=\begin{pmatrix}\int_{0}^{l}z_{4}(l_{0}-y)\phi^{s}(y)dy\\\
\int_{0}^{l}z_{3}(l_{0}-y)\phi^{s}(y)dy\\\
\int_{0}^{l}z_{2}(l_{0}-y)\phi^{s}(y)dy\\\
\frac{U_{0}}{EI}+\int_{0}^{l}z_{1}(l_{0}-y)\phi^{s}(y)dy-\beta\int_{0}^{l_{0}}z_{4}(l_{0}-y)\phi^{s}(y)dy\end{pmatrix},}$
(13)

and

$\beta=\frac{s^{2}m+\varkappa}{EI}.$

Thus, the vector-valued function $W^{s}(x)$ is defined by (9) and (11) with
$\bar{W}=M^{-1}R$, and the components of
$R=(R_{1}^{s},R_{2}^{s},R_{3}^{s},R_{4}^{s})^{T}$ are linear combinations of
$U_{0}$, …, $U_{k}$:

$R_{i}^{s}=\sum_{j=0}^{k}r_{ij}^{s}U_{j},\;\;i=1,2,3,4,$ (14)

where the coefficient matrix $r^{s}=(r_{ij}^{s})$ can be obtained from (8),
(13). Then the function $W_{3}^{s}(x)$, corresponding to the second
$x$-derivative of $w(x,t)$, is expressed as:

$\displaystyle
W_{3}^{s}(x)=-4\gamma^{4}z_{4}(x)\bar{W}^{0}_{2}+z_{2}(x)\bar{W}_{4}^{0}$
$\displaystyle+\sum_{j=1}^{k}\frac{U_{j}}{EI+\rho
ds}\int_{0}^{x}K_{j}(x,y)dy\;\;\text{for}\;x\in[0,l_{0}],$ $\displaystyle
W_{3}^{s}(x)=-4\gamma^{4}z_{4}(x-l)\bar{W}^{l}_{2}+z_{2}(x-l)\bar{W}_{4}^{l}$
(15) $\displaystyle-\sum_{j=1}^{k}\frac{U_{j}}{EI+\rho
ds}\int_{x}^{l}K_{j}(x,y)dy\;\text{for}\;x\in(l_{0},l],$

where $K_{j}(x,y)=z_{2}(x-y)\psi_{j}^{\prime\prime}(x)$ and
$(\bar{W}^{0}_{2},\bar{W}^{0}_{4},\bar{W}^{l}_{2},\bar{W}^{l}_{4})^{T}=M^{-1}r^{s}(U_{0},U_{1},...,U_{k})^{T}$
because of (14). Thus, at each $x\in[0,l]$, the above formulas define
$W_{3}^{s}(x)$ as a linear combination of $U_{0}$, $U_{1}$,…, $U_{k}$:

$W^{s}_{3}(x)=\sum_{j=0}^{k}h_{j}^{s}(x)U_{j},$ (16)

with the coefficients $h_{j}^{s}(x)$ collected from (15). Recalling that the
output of the considered system is given by (5), we summarize the computation
of the transfer function in the following lemma.

###### Lemma 1

The transfer function of the multi-input multi-output control system (2)–(5)
is presented in the form:

$H(s)=\begin{pmatrix}h_{0}^{s}(l_{1})&h_{1}^{s}(l_{1})&...&h_{k}^{s}(l_{1})\\\
h_{0}^{s}(l_{2})&h_{1}^{s}(l_{2})&...&h_{k}^{s}(l_{2})\\\
\vdots&\vdots&\ddots&\vdots\\\
h_{0}^{s}(l_{p})&h_{1}^{s}(l_{p})&...&h_{k}^{s}(l_{p})\\\ \end{pmatrix},$

where $h_{j}^{s}(x)$ are taken from (16).

###  3.1 Single-input single-output (SISO) case

Let the system be controlled by the shaker force $u_{0}$ only and the scalar
output signal $y_{1}(t)$ be available. In this particular case, the scalar
transfer function $H(s)$ is such that $Y_{1}(s)=H(s)U_{0}(s)$. Lemma 1 implies
the following result in the considered SISO case.

###### Lemma 2

Assume that $k=0$, $p=1$, and $l_{1}\leq l_{0}$. Then the transfer function of
the control system (2)–(5) is

$H(s)=\frac{1}{EI}\left(-4\gamma^{4}z_{4}(l_{j})M^{-1}_{14}+z_{2}(l_{j})M^{-1}_{24}\right).$
(17)

Here $M^{-1}_{ik}$ are elements of $M^{-1}$ (the matrix $M$ is given in (13))
and $z_{i}(x)$ are defined in (10).

For further numerical simulations, we consider the beam actuated by the shaker
force only ($k=0$) with single output ($p=1$) and take the following realistic
mechanical parameters [21] (see also [13]):

$l=1.905\,\text{m},\;l_{0}=1.4\,\text{m},\;\rho_{0}=2700\,\text{kg}/\text{m}^{3},$
$S=2.25\cdot 10^{-4}\text{m}^{2},\;\rho=\rho_{0}S,\;E=6.9\cdot
10^{10}\,\text{Pa},$ $I=1.6875\cdot
10^{-10}\,\text{m}^{4},\;m=0.1\,\text{kg},\;\varkappa=7\,\text{N/mm},$
$l_{1}=732.5\,\text{mm}.$ (18)

##  4 The Loewner framework for fitting linear time-invariant systems

In what follows, we provide a brief summary of the LF for fitting linear
dynamical systems as in (1), from data. The starting point for LF is having
access to measurements corresponding to the transfer function of the
underlying dynamical process, which can be inferred in practice by means of
experimental or model-based procedures. The data set is given by:

${\cal D}=\\{(\omega_{\ell};H(\omega_{\ell}))\ |\ \ell=1,\ldots,2k\\},$ (19)

by means of sampling $H:{\mathbb{C}}\rightarrow{\mathbb{C}}$ is an analytic
function (not necessarily rational) on a particular (complex) grid of points
$\omega_{\ell}$’s. It is to be noted that data sets with an odd number of
measurements can also be accommodated in the LF. The first step is to
partition the data set in (19) into two disjoint subsets, as follows:

$\displaystyle\begin{split}{\textrm{r}ight\ data}&:\ {\cal
D}_{R}=\\{(\lambda_{j};w_{j})\ |\ j=1,\ldots,k\\},\leavevmode\nobreak\
{\textrm{a}nd},\\\ {\textrm{l}eft\ data}&:\ {\cal D}_{L}=\\{(\mu_{i};v_{i})\
|\ \leavevmode\nobreak\ i=1,\ldots,k\\},\end{split}$ (20)

For simplicity, all points are assumed distinct and also
$\mu_{i}\neq\lambda_{j}$, for all $1\leq i,j\leq k$; extensions to Hermite
interpolation were proposed in [27].

A typical approach for splitting the data, commonly used in the LF
publications, is the “alternate splitting scheme”, described as follows. The
left and right sample nodes (and points) are chosen so that they are
interlacing each other. More precisely, for $1\leq i\leq k$, we can write:

$\displaystyle\begin{split}\mu_{i}&=\omega_{2i-1},\ \
\lambda_{i}=\omega_{2i},\\\ v_{i}=H(\mu_{i})=H(&\omega_{2i-1}),\ \
w_{i}=H(\lambda_{i})=H(\omega_{2i}).\end{split}$ (21)

We refer the reader to [22] for a more comprehensive account of data
partitioning strategies in the Loewner framework; there, the “half-hal”
splitting is mentioned together with approaches that split the data based on
the magnitude of the data samples.

The goal is to find a rational function denoted with $\tilde{H}(s)$, such that
the following interpolation conditions are (approximately) fulfilled:

$\tilde{H}(\mu_{i})=v_{i},\leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \tilde{H}(\lambda_{j})=w_{j}.$ (22)

In order to accomplish this scope, we first arrange the elements of the
original data set ${\cal D}$, partitioned as in (20) in matrix format. Hence,
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}},$
(23)

while the data vectors ${\mathbb{V}},{\mathbb{W}}^{T}\in{\mathbb{R}}^{k}$ are
given by:

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

The Loewner model is hence constructed as follows:

$\displaystyle{\mathbf{E}}=-{\mathbb{L}},\leavevmode\nobreak\
\leavevmode\nobreak\ {\mathbf{A}}=-{{{\mathbb{L}}_{s}}},\leavevmode\nobreak\
\leavevmode\nobreak\ {\mathbf{B}}={\mathbb{V}},\leavevmode\nobreak\
\leavevmode\nobreak\ {\mathbf{C}}={\mathbb{W}}.$

The following Sylvester equations are satisfied by the Loewner and shifted
Loewner matrices, as shown in [4] (here,
${\mathbbm{1}}_{q}=\left[\begin{matrix}1&\cdots&1\end{matrix}\right]^{T}\in{\mathbb{C}}^{q}$):

$\begin{cases}{\mathbf{M}}{\mathbb{L}}-{\mathbb{L}}\boldsymbol{\Lambda}={\mathbb{V}}{\mathbbm{1}}_{k}^{T}-{\mathbbm{1}}_{q}{\mathbb{W}},\\\
{\mathbf{M}}{{{\mathbb{L}}_{s}}}-{{{\mathbb{L}}_{s}}}\boldsymbol{\Lambda}={\mathbf{M}}{\mathbb{V}}{\mathbbm{1}}_{k}^{T}-{\mathbbm{1}}_{q}{\mathbb{W}}\boldsymbol{\Lambda},\end{cases}\vspace{-2mm}$
(25)

where ${\mathbf{M}}=\text{diag}(\mu_{1},\cdots,\mu_{q})$ and
$\boldsymbol{\Lambda}=\text{diag}(\lambda_{1},\cdots,\lambda_{k})$. The
following relations expressing the shifted Loewner matrix to the Loewner
matrix, in two distinct ways, hold:

${{{\mathbb{L}}_{s}}}={\mathbb{L}}\boldsymbol{\Lambda}+{\mathbb{V}}{\mathbbm{1}}_{k}^{T}={\mathbf{M}}{\mathbb{L}}+{\mathbbm{1}}_{q}{\mathbb{W}}.$
(26)

Hence, the explicit computation of large Loewner matrices can be avoided by
means of computing (approximated, low-rank) solutions of the Sylvester
equations in (25). This can be accomplished, e.g., by means of using optimized
and robust numerical tools such as [8].

Provided that enough data are available, the pencil
$({{{\mathbb{L}}_{s}}},\,{\mathbb{L}})$ is often singular. For example, if the
data ${\cal D}$ in 19 were generated from a rational function $H(s)$
corresponding to a minimal LTI system of dimension $n$, i.e.,
$H(s)={\mathbf{C}}(s{\mathbf{I}}_{n}-{\mathbf{A}})^{-1}{\mathbf{B}}$, then
this corresponds to the case $k>n$. Then, in the perfect setup (no noisy
data), one would encounter $k-n$ zero singular values when performing the
singular value decomposition (SVD) of the pencil
$\zeta{\mathbb{L}}-{{{\mathbb{L}}_{s}}}$, where $\zeta\in{\mathbb{C}}$ is
chosen to be different than the eigenvalues of matrix ${\mathbf{A}}$. In such
cases, an SVD of augmented Loewner matrices is computed, and the dominating
part is selected as:

$\left[{\mathbb{L}},\leavevmode\nobreak\
{{{\mathbb{L}}_{s}}}\right]\approx{\mathbf{Y}}\widehat{\Sigma}_{{r}}\tilde{{\mathbf{X}}}^{*},\leavevmode\nobreak\
\left[\begin{array}[]{c}{\mathbb{L}}\\\
{{{\mathbb{L}}_{s}}}\end{array}\right]\approx{\tilde{\mathbf{Y}}}\Sigma_{{r}}{\mathbf{X}}^{*},$
(27)

where $\widehat{\Sigma}_{{r}}$, $\Sigma_{{r}}$
$\in{\mathbb{R}}^{{{r}}\times{r}}$,
${\mathbf{Y}}\in{\mathbb{C}}^{k\times{r}}$,${\mathbf{X}}\in{\mathbb{C}}^{k\times{r}}$,
$\tilde{{\mathbf{Y}}}\in{\mathbb{C}}^{2k\times{r}}$,
$\tilde{{\mathbf{X}}}\in{\mathbb{C}}^{r\times{2k}}$ and
$({\mathbf{X}})^{*}\in{\mathbb{C}}^{r\times k}$ denotes the conjugate-
transpose of matrix ${\mathbf{X}}$. This is performed in order to find
projection matrices ${\mathbf{X}}_{r},{\mathbf{Y}}_{r}\in{\mathbb{C}}^{k\times
r}$, as described in [4]. Here, $r<n$ represents the truncation index. Then,
the system matrices corresponding to a projected Loewner model of dimension
$r$ can be computed using the truncated singular vector matrices
${\mathbf{X}}_{r}$ and ${\mathbf{Y}}_{r}$:

$\displaystyle\begin{split}\hat{{\mathbf{E}}}&=-{\mathbf{X}}_{r}^{*}{\mathbb{L}}{\mathbf{Y}}_{r},\
\
\hat{{\mathbf{A}}}=-{\mathbf{X}}_{r}^{*}{{{\mathbb{L}}_{s}}}{\mathbf{Y}}_{r},\\\
\hat{{\mathbf{B}}}&={\mathbf{X}}_{r}^{*}{\mathbb{V}},\ \
\hat{{\mathbf{C}}}={\mathbb{W}}{\mathbf{Y}}_{r},\end{split}$ (28)

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

$\begin{cases}{\hat{\textbf{E}}}\dot{{\hat{\textbf{x}}}}(t)={\hat{\textbf{A}}}{\mathbf{x}}(t)+{\hat{\textbf{B}}}u(t),\\\
\hat{y}(t)\ \hskip
5.12149pt={\hat{\textbf{C}}}{\hat{\textbf{x}}}(t).\end{cases}.$ (29)

The transfer function of the reduced Loewner model in (29) is written as
${\hat{H}}(s)={\hat{\textbf{C}}}(s{\hat{\textbf{E}}}-{\hat{\textbf{A}}})^{-1}{\hat{\textbf{B}}}$,
and it provides a good approximant to the original transfer function $H(s)$.
Then, ${\hat{H}}(s)$ may be expanded in a pole/zero or pole/residue format.
These values represent system invariants and can be related to the inherent
dynamics. It is noted that the state-space realization is not unique, and that
is why an extra step is required. More implementation details and properties
of the LF procedure can be found in [4, 22].

##  5 Numerical examples

###  5.1 Analysis on the unperturbed data

In this section, we present two numerical test cases based on sampling the
transfer function explicitly derived in Section 3. We consider the SISO system
with the choice of physical parameters (18). Then the scalar transfer function
$H(s)$ in Lemma 2 is sampled at the purely imaginary grid points
$s=\omega_{\ell}$, $\ell=1,\ldots,k$ with $k=1000$, which are equally
distributed in the range of physical frequencies from $0Hz$ to $250Hz$. The
data partitioning scheme (into the left and right disjoint subsets) chosen
here is the “alternate” one, previously described in Section 4. We also note
that complex conjugated data is added to the process to enforce real-valued
models; more details on how this is achieved can be found in [4].

Two values of the structural damping parameter $d$ in (2) are considered for
the numerical simulations: $d_{1}=0.0249$ (“large” damping) and $d_{2}=0.001$
(“small” damping). The case $d=d_{1}$ is depicted in Figs. 1–4, for which we
are fitting a Loewner model of order $r=20$ (with a rational transfer
function). Additionally, the case $d=d_{2}$ is presented in Figs. 5–8; for
this case, we are fitting a Loewner model of order $r=27$. It should be noted
that the Loewner framework does not automatically impose stability; post-
processing methods can be applied whenever unstable poles appear, as described
in [15].

We observe that the Loewner model approximates the original transfer function
with good accuracy in the whole range of test frequencies (the approximation
error of order $10^{-6}$ in Fig. 3, and of order $10^{-5}$ in Fig. 7). It is
also clear that these approximate models preserve the stability property (the
maximal real part of $\lambda_{i}$ in Fig. 4 is
${\textrm{m}ax}({\textrm{R}e}\,\lambda_{i})=-5.2780\cdot 10^{-1}<0$, while the
maximal real part of $\lambda_{i}$ in Fig. 8 is
${\textrm{m}ax}({\textrm{R}e}\,\lambda_{i})=-2.1187\cdot 10^{-2}<0$).

Although the transfer function of the considered distributed parameter
function is not rational, the decay of singular values associated with the
Loewner matrices in Figs. 1 and 5 seems to indicate the opportunity to enforce
rational approximation. More precisely, in both cases, a plateau (flat portion
of the graph) is observed after a steep decay. In Figs. 1 and 5, we only
depict the first $50$ singular values (out of $1000$, which is the dimension
of the Loewner matrices ${\mathbb{L}}$ and ${{{\mathbb{L}}_{s}}}$). In the
first test case, the decay is faster than in the second one. The dimensions of
the reduced-order models were chosen in accordance with this phenomenon (i.e.
$r=20$ for the first and $r=27$ for the second). However, such a clear and
steep decay as seen in Figs. 1 and 5, is seldom noticed in experimental data.
In such scenarios, the data can be perturbed, i.e., by means of noise. We
treat this case below.

Figure 1: The decay of the singular values for the augmented Loewner matrices
in 27. Figure 2: The original data (transfer function samples) vs. the Loewner
model fit. Figure 3: The approximation error. Figure 4: The poles of the
fitted Loewner model, as eigenvalues of pencil
$({\hat{\textbf{A}}},{\hat{\textbf{E}}})$. Figure 5: The decay of the singular
values for the augmented Loewner matrices in 27. Figure 6: The original data
(transfer function samples) vs. the Loewner model fit. Figure 7: The
approximation error. Figure 8: The poles of the fitted Loewner model, as
eigenvalues of pencil $({\hat{\textbf{A}}},{\hat{\textbf{E}}})$.

###  5.2 Analysis on the perturbed data (by means of artificial additive
Gaussian noise)

In this subsection, we analyze the robustness of the LF when applied to
perturbed (noisy) data. A preliminary analysis of such endeavors was reported
in [25, 12, 23] and in [35, 14]. Similarly to the approaches in these
publications, we will include additive Gaussian noise into the measurements of
the transfer function $H(s)$ from Lemma 2. More precisely, for all $k=1000$
previous grid points $s=\omega_{\ell}$, the new data are, for
$\ell=1,\ldots,k$ and $\nu=1,\ldots,4$:

$H(\omega_{\ell})[1+\epsilon^{(\nu)}(\alpha_{\ell}+\imath\beta_{\ell})].$ (30)

Here, the “noise power” $\epsilon^{(\nu)}>0$ is chosen so that
$\epsilon^{(\nu)}=10^{-\nu}$, while the values $\alpha_{\ell}$ and
$\beta_{\ell}$ are drawn from the standard normal distribution. For
simplicity, we use the term “noise level $\nu$” (in the subsequent text and
plots). Also, for the sake of brevity, we will restrict our analysis to the
case of a smaller damping coefficient: $d=0.001$.

The first numerical experiment that is shown here is concerned with the decay
of singular values for Loewner matrices under the influence of noise. As
previously pointed out in [25, 18], the effect of noise is generally reflected
in the “flattening” of the singular value curve. This is especially valid for
the “alternate” splitting, and not as much for the “half-half” splitting (as
shown in [18]). This is precisely the phenomenon observed in Fig. 9.
Additionally, there are a number of dominant singular values that are less
prone to perturbations; depending on the noise levels, it is clear that this
number decreases with the increase of the noise power. For example, when
$\nu=2$ there are 22 singular values that seem to be stagnant. This is hence
another effect of the noise in data for LF; deciding the order of the fitted
model becomes a more challenging task, and one needs to be careful to avoid
over-fitting.

Figure 9: The decay of singular values for the augmented Loewner matrices for
different ”noise levels/powers”.

Next we fix the noise level at $\nu=2$ and record the poles of the fitted
Loewner model on the noisy data. We compare those with the original poles, and
the results are depicted in Fig. 10. As expected, most of the (dominant) poles
seem to be matching well, with the remark that the noisy data introduces
spurious poles.

Figure 10: The poles of the fitted Loewner models (with and without noise).

However, as shown in Fig. 11, the effects of the noise (for this level) do not
seem to be drastic; the response of the Loewner model fitted to the noisy data
faithfully follows the original (unperturbed) data.

Figure 11: The original data vs. the Loewner model fit (on the noisy data for
level $\nu=2$) $\&$ the approximation error.

It is to be noted that for levels of noise higher than the one used in the
previous experiment, i.e., for $\nu=1$, the results are significantly less
accurate. In that case, for some perturbed data sets, the first two dominant
peaks are partially or completely missed (as shown in Fig. 12). Additionally,
it was noticed that the asymptotic stability of the reduced-order Loewner
model was sometimes lost (however not in the experiment reported here). This
behavior is mostly due to the high level of the noise signal added here. As
reported in [25], for high signal-to-noise scenarios, it is very challenging
to extract all the meaningful information from the perturbed data. A more
thorough numerical analysis based on the so-called stabilization diagrams will
be left for future research endeavors.

Figure 12: The original data vs. the Loewner model fit (on the noisy data for
level $\nu=1$) $\&$ the approximation error.

##  6 Conclusion and outlook

The contribution of this work is twofold. First, an analytic construction of
the transfer function for an infinite-dimensional flexible structure with the
Euler–Bernoulli beam and a spring-mass system has been proposed. Second, we
have applied the Loewner framework (LF) for data-driven modeling of the
considered class of flexible structures on the basis of transfer function
measurements. The presented case study of the singular value decay for the
Loewner pencil clearly indicates a possible choice of the dimension of an
acceptable reduced-order model, depending on the damping parameter $d$. Note
that our results are not limited to deterministic measurements since the
performed study evaluates the effects of noise in the LF for data sets with
Gaussian perturbations as well. A thorough numerical analysis of the perturbed
(noisy) data has been carried out. The preliminary results show the robustness
of the method to low and moderate levels of noise, but also point out some
challenges for the cases in which the data are perturbed with higher noise
levels. For future research endeavors, we intend to tackle the following open
issues:

  * •

Analyze data from real experimental measurements to supplement the information
attained from simulation data.

  * •

Take into account and quantify the effects of measurement noise in the LF for
such data sets (using pseudospectra theory [14]).

  * •

Study the decay of the singular values depending on the influence of the
damping parameter and on data splitting strategies. Extend this study to
flexible structures with a different asymptotic distribution of the
eigenvalues, e.g., the Timoshenko beam with attached rigid bodies [36].

  * •

Examine the applicability of the Loewner-based reduced-order systems for the
control design of the original infinite-dimensional plant.

## References

  * [1] B. Altıner, A. Delibaşı, and B. Erol.  Modeling and control of flexible link manipulators for unmodeled dynamics effect.  Proceedings of the Institution of Mechanical Engineers, Part I: Journal of Systems and Control Engineering, 233(3):245–263, 2019.
  * [2] A. C. Antoulas.  Approximation of large-scale dynamical systems.  SIAM, Philadelphia, 2005.
  * [3] A. C. Antoulas, C. A. Beattie, and S. Gugercin.  Interpolatory Methods for Model Reduction.  SIAM, Philadelphia, 2020.
  * [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] C. Beattie and S. Gugercin.  Interpolatory projection methods for structure-preserving model reduction.  Systems & Control Letters, 58(3):225–232, 2009.
  * [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, S. Gugercin, and S. W. R. Werner.  Structure-preserving interpolation for model reduction of parametric bilinear systems.  Automatica, 132:109799, 2021.
  * [8] P. Benner, M. Köhler, and J. Saak.  Matrix equations, sparse solvers: M-M.E.S.S.-2.0.1 – philosophy, features and application for (parametric) model order reduction.  In Model Reduction of Complex Dynamical Systems, volume 171 of International Series of Numerical Mathematics, pages 369–392. Birkhäuser, Cham, 2021.
  * [9] P. Benner, M. Ohlberger, A. Cohen, and K. Willcox.  Model Reduction and Approximation.  Society for Industrial and Applied Mathematics, Philadelphia, PA, 2017\.
  * [10] D. Casagrande, W. Krajewski, and U. Viaro.  The integer–order approximation of fractional–order systems in the Loewner framework.  IFAC-PapersOnLine, 52(3):43–48, 2019.
  * [11] R. Curtain and K. Morris.  Transfer functions of distributed parameter systems: A tutorial.  Automatica, 45(5):1101–1116, 2009.
  * [12] Z. Drmač and B. Peherstorfer.  Learning low-dimensional dynamical-system models from noisy frequency-response data with Loewner rational interpolation.  In Realization and Model Reduction of Dynamical Systems, pages 39–57. Springer, 2022.
  * [13] C. Dullinger, A. Schirrer, and M. Kozek.  Advanced control education: optimal & robust MIMO control of a flexible beam setup.  IFAC Proceedings Volumes, 47(3):9019–9025, 2014.
  * [14] M. Embree and A. C. Ionita.  Pseudospectra of Loewner matrix pencils.  In Realization and Model Reduction of Dynamical Systems, pages 59–78. Springer, 2022.
  * [15] I. V. Gosea and A. C. Antoulas.  Stability preserving post-processing methods applied in the Loewner framework.  In 2016 IEEE 20th Workshop on Signal and Power Integrity (SPI), pages 1–4, 2016.
  * [16] I. V. Gosea and I. Pontes Duff.  An iterative realization-free approach for model reduction of bilinear systems via Hermitian interpolation.  In 2022 European Control Conference (ECC), pages 584–589. IEEE, 2022.
  * [17] I. V. Gosea, C. Poussot-Vassal, and A. C. Antoulas.  On Loewner data-driven control for infinite-dimensional systems.  In 2021 European Control Conference (ECC), pages 93–99. IEEE, 2021\.
  * [18] I. V. Gosea, Q. Zhang, and A. C. Antoulas.  Data-driven modeling from noisy measurements.  PAMM, 20(S1):e202000358, 2021.
  * [19] B. Gustavsen and A. Semlyen.  Rational approximation of frequency domain responses by vector fitting.  IEEETransPD, 14(3):1052–1061, 1999.
  * [20] J. Kalosha and A. Zuyev.  Asymptotic stabilization of a flexible beam with an attached mass.  Ukrainian Mathematical Journal, 73:1537–1550, 2022.
  * [21] J. Kalosha, A. Zuyev, and P. Benner.  On the eigenvalue distribution for a beam with attached masses.  In Stabilization of Distributed Parameter Systems: Design Methods and Applications, pages 43–56. Springer, 2021.
  * [22] D. S. Karachalios, I. V. Gosea, and A. C. Antoulas.  The Loewner framework for system identification and reduction.  In Model Order Reduction: Volume I: System-and Data-Driven Methods and Algorithms, pages 181–228. De Gruyter, 2021.
  * [23] P. Kergus and I. V. Gosea.  Data-driven approximation and reduction from noisy data in matrix pencil frameworks.  arXiv preprint arXiv:2202.09568, to appear in the MTNS22 proceedings volume, 2022.
  * [24] J. N. Kutz, S. L. Brunton, B. W. Brunton, and J. L. Proctor.  Dynamic mode decomposition: data-driven modeling of complex systems.  SIAM, 2016.
  * [25] S. Lefteriu, A. C. Ionita, and A. C. Antoulas.  Modeling systems based on noisy frequency and time domain measurements.  In Perspectives in Mathematical System Theory, Control, and Signal Processing, pages 365–378. Springer, 2010.
  * [26] P. Lietaert and K. Meerbergen.  Comparing Loewner and Krylov based model order reduction for time delay systems.  In 2018 European Control Conference (ECC), pages 545–550, 2018\.
  * [27] A. J. Mayo and A. C. Antoulas.  A framework for the solution of the generalized realization problem.  Linear Algebra and Its Applications, 425(2-3):634–662, 2007.
  * [28] Y. Nakatsukasa, O. Sete, and L. N. Trefethen.  The AAA algorithm for rational approximation.  SIAM Journal on Scientific Computing, 40(3):A1494–A1522, 2018.
  * [29] 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.
  * [30] I. Pontes-Duff, C. Poussot-Vassal, and C. Seren.  Realization independent single time-delay dynamical model interpolation and $\mathcal{H}_{2}$-optimal approximation.  In 2015 IEEE Conference on Decision and Control (CDC), pages 4662–4667. IEEE, 2015.
  * [31] C. Poussot-Vassal, P. Kergus, and P. Vuillemin.  Interpolation-based irrational model control design and stability analysis.  In Realization and Model Reduction of Dynamical Systems, pages 353–371. Springer, 2022.
  * [32] A. Quarteroni and G. Rozza.  Reduced Order Methods for Modeling and Computational Reduction, volume 9 of MS&A – Modeling, Simulation and Applications.  Springer International Publishing, Cham, Switzerland, 2014.
  * [33] P. Schulze and B. Unger.  Data-driven interpolation of dynamical systems with delay.  Systems & Control Letters, 97:125–131, 2016.
  * [34] P. Schulze, B. Unger, C. Beattie, and S. Gugercin.  Data-driven structured realization.  Linear Algebra and its Applications, 537:250–286, 2018.
  * [35] Q. Zhang, I. V. Gosea, and A. C. Antoulas.  Factorization of the Loewner matrix pencil and its consequences.  arXiv preprint arXiv:2103.09674, 2021.
  * [36] Alexander Zuyev and Oliver Sawodny.  Stabilization and observability of a rotating Timoshenko beam model.  Mathematical Problems in Engineering, 2007.