# Port-Hamiltonian Dynamic Mode Decomposition

Riccardo Morandin† and Jonas Nicodemus⋆ and Benjamin Unger⋆ † Institute of
Mathematics MA 4-5, Technical University Berlin, Straße des 17. Juni 136,
10623 Berlin, Germany<EMAIL_ADDRESS>⋆ Stuttgart Center for
Simulation Science (SC SimTech), University of Stuttgart, Universitätsstr. 32,
70569 Stuttgart, Germany<EMAIL_ADDRESS>stuttgart.de

###### Abstract.

We present a novel physics-informed system identification method to construct
a passive linear time-invariant system. In more detail, for a given quadratic
energy functional, measurements of the input, state, and output of a system in
the time domain, we find a realization that approximates the data well while
guaranteeing that the energy functional satisfies a dissipation inequality. To
this end, we use the framework of port-Hamiltonian (pH) systems and modify the
dynamic mode decomposition, respectively operator inference, to be feasible
for continuous-time pH systems. We propose an iterative numerical method to
solve the corresponding least-squares minimization problem. We construct an
effective initialization of the algorithm by studying the least-squares
problem in a weighted norm, for which we present the analytical minimum-norm
solution. The efficiency of the proposed method is demonstrated with several
numerical examples.

Keywords: dynamic mode decomposition, port-Hamiltonian systems, system
identification, dissipation inequality, passivity, knowledge-driven
realization

AMS subject classification: 37J06,37M99,65P10,93A30,93B30,93C05

##  1\. Introduction

Incorporating prior knowledge into modern learning architectures becomes
increasingly important in several applications. Such knowledge-driven or
physics-informed approaches [26, 39] exploit expert knowledge during the
learning process, either by optimizing only over a suitable set of candidate
functions, or by using physical information in the cost functional. In our
work, we deal with data from physical systems with the goal of identifying a
linear dynamical system, which we refer to as a realization that approximates
the data as well as possible with respect to the Frobenius norm. In more
detail, assume that we have measurements
$(u(t_{i}),x(t_{i}),y(t_{i}))\in\mathbb{R}^{m}\times\mathbb{R}^{n}\times\mathbb{R}^{m}$
at time instances $t_{i}$ for $i=0,\ldots,M$. Then, we want to determine
matrices $A\in\mathbb{R}^{n\times n}$, $B\in\mathbb{R}^{n\times m}$,
$C\in\mathbb{R}^{m\times n}$, and $D\in\mathbb{R}^{m\times m}$ such that the
data can be approximately recovered by the linear time-invariant system

$\displaystyle\dot{x}$ $\displaystyle=Ax+Bu,$ (1.1) $\displaystyle y$
$\displaystyle=Cx+Du.$

Since we assume that the data is based on a physical process, we want to
ensure that the realization satisfies a dissipation inequality, i.e., that the
rate of change of the energy associated with the system is bounded by the
externally supplied energy. We thus incorporate physical knowledge by
prescribing the (quadratic) energy functional

$\mathcal{H}(x)\vcentcolon=\tfrac{1}{2}x^{\mathsf{T}}Hx$ (1.2)

with symmetric positive definite matrix $H\in\mathbb{R}^{n\times n}$. Our main
goal is then to determine the matrices in (1.1) such that any solution of
(1.1) satisfies the dissipation inequality

$\tfrac{\mathrm{d}}{\mathrm{d}t}\mathcal{H}(x(t))\leq y(t)^{\mathsf{T}}u(t)$
(1.3)

for any $t$. One of the main advantages of requiring the learned model to
satisfy a dissipation inequality is that whenever the model is coupled with
another passive model via a power-conserving or dissipative interconnection,
then the coupled model is also passive. Moreover, since the Hamiltonian also
serves as a Lyapunov function, we are guaranteed that the identified system is
stable (independent of the underlying physical system and the quality of the
measurements). Our framework can thus be used to guarantee the physical
behavior of coupled first-principle and purely data-inferred dynamical
systems. To achieve this goal, we use the framework of _port-Hamiltonian_ (pH)
systems [24, 50] and modify the _dynamic mode decomposition_ (DMD) [41, 49,
28] and _operator inference_ (OI) [35] accordingly. Our main contributions are
the following:

  1. (i)

Since DMD is designed to compute a discrete-time dynamical system, we follow
[27, 30] and present a definition of a discrete-time pH system in section 2.2,
which is motivated from the structure-preserving time-discretization of a
continuous-time pH system. The corresponding modified DMD optimization problem
is formulated in 3.1.

  2. (ii)

Although 3.1 is convex and solvable, a closed-form solution formula is not
immediately available. Instead, we propose an iterative method (Algorithm 1)
combining the result of a skew-symmetric Procrustes problem [15] with a
projected fast gradient method for a positive semidefinite Procrustes problem
[17].

  3. (iii)

For an efficient initialization of Algorithm 1, we consider a weighted
Frobenius norm, where we weight the problem according to the relevant
information in the data, see the forthcoming Section 3.3. The analytic
minimum-norm solution of the weighted problem is then presented in Theorem 3.7
and used, up to some modification, as initialization for our iterative method.

###  1.1. Literature review

The construction of a realization of the form (1.1) from data is a well-
studied subject with many popular approaches. We mention the eigensystem
realization algorithm [25], the Loewner framework [29], vector fitting [19],
neural ordinary differential equations [13], OI [35], and DMD [28].
Introducing expert knowledge to these approaches is not a novel idea and may
even be the key idea to quantifying the error between the realization and the
true physical system, cf. [22, 23]. Exploiting expert knowledge in the
identification of linear systems is reported in various applications and
methods. For instance, the specific structure of mechanical systems is
exploited in a vector fitting framework in [51]. The authors of [42, 43]
exploit the fact that certain wave-type phenomena can be represented with
delay equations to construct accurate linear time-invariant surrogate models
from data. In the context of the Navier-Stokes equation, the specific
structure of the (semi-discretized) equations is exploited in an OI approach
in [7]. Extensions of OI to Hamiltonian and Lagrangian dynamics are reported
in [48, 47]. Similar ideas are discussed in the context of DMD in [4] by
restricting the discrete-time iteration matrix to specific manifolds. In
essence, our method relies on a similar idea but with a specific structure of
the iteration matrix not discussed in [4].

As detailed above, we use the framework of pH systems to achieve our goal of
identifying a system that satisfies a dissipation inequality. In contrast to
our approach, most existing results on learning linear time-invariant pH
systems work in the frequency domain. Using rational interpolation, a pH
system is constructed in [3] within the Loewner framework by interpolating the
transfer function at the spectral zeros. Since these are typically not known a
priori, the authors of [8] propose first identifying an (unstructured) system
and then computing the spectral zeros from this system. A parameterization of
pH systems is used in [45, 44, 46] to approximately mimic an
$\mathcal{H}_{\infty}$-type cost functional. Note that these ideas can be
transferred to time-domain data using similar ideas as, for instance, in [34].
Instead of directly identifying a pH structure (as we do in our contribution),
the authors of [14] propose to first identify an unstructured model and then
find the nearest pH systems. Let us emphasize that in contrast to the other
methods, we assume knowledge of the energy functional, and instead of finding
an arbitrary pH system, our goal is to identify a system such that the
dissipation inequality for this energy functional is guaranteed.

We mention that many of the identification methods discussed above are closely
related to intrusive structure-preserving model order reduction methods for pH
systems. Popular approaches include the effort-and flow-constraint reduction
method [37], tangential interpolation [52, 36, 18], generalized balancing [11,
10], $\mathcal{H}_{2}$-and $\mathcal{H}_{\infty}$-optimal approaches [40, 45,
32], and spectral factorization [12]. For further methods, we refer to the
references cited in [31, Rem. 8.2].

###  1.2. Organization of the manuscript

After this introduction, we recall basic results on DMD and OI (cf. section
2.1) and review the theory of linear time-invariant pH systems in section 2.2.
In particular, we introduce a definition for a discrete-time pH system,
following the more general discussion in [27, 30]. The _port-Hamiltonian
dynamic mode decomposition_ (pHDMD) problem is formulated in section 3.1,
yielding a convex minimization problem. An iterative algorithm to solve the
pHDMD problem is presented in section 3.2, whereas we discuss a clever
initialization in section 3.3 by solving the pHDMD problem in a weighted norm.
The algorithm is then applied to several numerical examples in section 4.

###  1.3. Notation

We use the symbols $\mathbb{N}$, $\mathbb{R}$, $\mathbb{R}^{n}$, and
$\mathbb{R}^{n\times m}$ to denote the positive integers, the real numbers,
the set of columns vectors with $n\in\mathbb{N}$ real entries, and the set of
$n\times m$ real matrices, respectively. For a matrix $A\in\mathbb{R}^{n\times
m}$ we denote its transpose with $A^{\mathsf{T}}$ and its Moore-Penrose
pseudoinverse with $A^{\dagger}$. A matrix $A\in\mathbb{R}^{n\times n}$ is
called symmetric (skew-symmetric), if $A=A^{\mathsf{T}}$
($A=-A^{\mathsf{T}}$). The decomposition of a squared matrix
$A\in\mathbb{R}^{n\times n}$ into its skew-symmetric resp. symmetric part is
denoted as

$\operatorname{sym}(A)=\tfrac{1}{2}(A+A^{\mathsf{T}}),\qquad\mathrm{resp.}\qquad\operatorname{skew}(A)=\tfrac{1}{2}(A-A^{\mathsf{T}}).$

The sets of all $n\times n$ symmetric positive definite and symmetric positive
semi-definite matrices with real entries are denoted with
$\mathcal{S}^{n}_{\succ}$ and $\mathcal{S}^{n}_{\succeq}$. The projection
$\mathcal{P}_{\succeq}(A)$ of a matrix $A\in\mathbb{R}^{n\times n}$ onto
$\mathcal{S}^{n}_{\succeq}$ is given by

$\mathcal{P}_{\succeq}(A)=\Xi\operatorname{diag}(\Lambda_{+},0,0)\Xi^{\mathsf{T}},$

where
$\operatorname{sym}(A)=\Xi\operatorname{diag}(\Lambda_{+},\Lambda_{-},0)\Xi^{\mathsf{T}}$
is the ordered eigendecomposition of the symmetric part of $A$. Here,
$\Lambda_{+}$ and $\Lambda_{-}$ are the diagonal matrices containing the
positive and negative eigenvalues, respectively. The set of nonsingular
matrices of size $n\times n$ is denoted with $\mathrm{GL}_{n}(\mathbb{R})$.
The Stiefel manifold of $n\times r$ dimensional matrices with real entries is
denoted by

$\mathrm{St}(n,r)\vcentcolon=\left\\{U\in\mathbb{R}^{n\times r}\mid
U^{\mathsf{T}}U=I_{r}\right\\},$

where $I_{r}$ denotes the $r\times r$ identity matrix. Furthermore, given
$\Omega\in\mathcal{S}^{n}_{\succeq}$ we denote the weighted Frobenius semi-
norm as

$\left\|A\right\|_{\Omega}\vcentcolon=\sqrt{\operatorname{trace}{(A^{\mathsf{T}}\Omega
A)}}.$ (1.4)

##  2\. Preliminaries

###  2.1. Dynamic Mode Decomposition and Operator Inference

Assume data triples
$(u_{i},x_{i},y_{i})\in\mathbb{R}^{m}\times\mathbb{R}^{n}\times\mathbb{R}^{p}$
(DMD) or data quadruples
$(u_{i},x_{i},\dot{x}_{i},y_{i})\in\mathbb{R}^{m}\times\mathbb{R}^{n}\times\mathbb{R}^{n}\times\mathbb{R}^{p}$
(OI) for $i=0,\ldots,M$ of a dynamical system available, which may be obtained
from measurements of real phenomena or the simulation of a model. The goal of
(input-output) DMD [38, 2] or OI [35] is to find matrices
$\widetilde{A},\widetilde{B},\widetilde{C},\widetilde{D}$ of suitable size
minimizing

$\min_{\widetilde{A},\widetilde{B},\widetilde{C},\widetilde{D}}\sum_{i=0}^{M-1}\|\Delta
x_{i}-\widetilde{A}x_{i}-\widetilde{B}u_{i}\|_{2}^{2}+\|y_{i}-\widetilde{C}x_{i}-\widetilde{D}u_{i}\|_{2}^{2},$
(2.1)

where $\Delta x_{i}\vcentcolon=\dot{x}_{i}$ if we assume derivative
information of the state to be available, and $\Delta
x_{i}\vcentcolon=x_{i+1}$, otherwise. The goal of solving (2.1) is to
determine a linear system of the form

$\displaystyle\Delta\widetilde{x}$
$\displaystyle=\widetilde{A}\widetilde{x}+\widetilde{B}u,$ (2.2)
$\displaystyle\widetilde{y}$
$\displaystyle=\widetilde{C}\widetilde{x}+\widetilde{D}u,$

that best approximates the data (in the sense of (2.1)), where $\Delta$
denotes the differentiation operator with respect to time if we are in the
continuous time setting described by OI, and $\Delta$ denotes the forward
shift in the discrete-time setting used in DMD. This is conveniently achieved
by introducing the matrices

$\displaystyle Z_{0}\vcentcolon=\begin{bmatrix}x_{0}&\ldots&x_{M-1}\\\
u_{0}&\ldots&u_{M-1}\end{bmatrix}\in\mathbb{R}^{(n+m)\times
M}\qquad\text{and}\qquad Z_{1}\vcentcolon=\begin{bmatrix}\Delta
x_{0}&\ldots&\Delta x_{M-1}\\\
y_{0}&\ldots&y_{M-1}\end{bmatrix}\in\mathbb{R}^{(n+p)\times M}$

and studying the equivalent problem

$\min_{\mathcal{A}\in\mathbb{R}^{(n+p)\times(n+m)}}\left\|Z_{1}-\mathcal{A}Z_{0}\right\|_{\mathrm{F}},$
(2.3)

whose solution, conformably partitioned, yields the matrices
$\widetilde{A},\widetilde{B},\widetilde{C}$, and $\widetilde{D}$. The minimum
norm solution $\mathcal{A}^{\star}$ of (2.3) is given by the Moore-Penrose
pseudo-inverse of the data, i.e., $\mathcal{A}^{\star}=Z_{1}Z_{0}^{\dagger}$.
We emphasize that this is the unique solution of (2.3) if and only if
$r\vcentcolon=\operatorname{rank}(Z_{0})=n+m$. The pseudo-inverse can be
computed efficiently using the _singular value decomposition_ (SVD). In more
detail, set $r\vcentcolon=\operatorname{rank}(Z_{0})$ and let $Z_{0}=V\Sigma
W^{\mathsf{T}}$ with $V\in\mathrm{St}(n+m,r)$,
$\Sigma\in\mathrm{GL}_{r}(\mathbb{R})$, $W\in\mathrm{St}(M,r)$ denote the
skinny SVD of $Z_{0}$. Then, $Z_{0}^{\dagger}=W\Sigma^{-1}V$. Note that from a
numerical point of view, we truncate singular values below a given tolerance
to ensure that the inverse $\Sigma^{-1}$, respectively the associated linear
system, can be computed accurately. Such a truncation is equivalent to adding
a regularization term to the minimization problem (2.3), see [9] for further
details.

######  Remark 2.1.

If no measurements of the derivative are available, then a classical finite
difference approach of the form

$\displaystyle\Delta
x_{i}\approx\dot{\hat{x}}_{i}\vcentcolon=\frac{x_{i+1}-x_{i}}{\delta_{i}},$
(2.4)

may be used as an approximation, where $\delta_{i}$ is the time step between
the snapshots. This is a common approach in the literature, see, e.g., [35],
for which a convergence result for $\delta_{i}\to 0$ is available [35, Thm.
1]. Throughout this manuscript, we do not assume $\dot{x}_{i}$ available, and
hence work with the approximation (2.4), instead.

The method can be further extended to include dimensionality reduction, which
is then called _reduced OI_ [35]. There, a low dimensional basis for the state
space must is constructed. Afterwards, the data is projected onto this low-
dimensional basis. Finally, (2.1) is solved for the reduced data, which
results in a _reduced order model_. The low dimensional basis generation can
be achieved, for instance, by applying the SVD to the state data, i.e., let
$X=V\Sigma W$ denote the skinny SVD of $X$. Then the reduced basis $\Phi$ of
size $r$ can be defined as the first $r$ columns of $V$.

###  2.2. Port-Hamiltonian system

As outlined in the introduction, pH systems are a promising modeling tool that
generalize the notion of Hamiltonian systems to allow for interaction with the
environment. Since the aim of DMD is to produce a linear system, we review the
pH framework only for linear dynamics and refer the reader for a more general
definition to [50, 31].

######  Definition 2.2 (Linear time-invariant pH system).

Assume that we have given a quadratic function
$\mathcal{H}(x)=\tfrac{1}{2}x^{\mathsf{T}}\widehat{H}x$, called the
_Hamiltonian_ , with $\widehat{H}\in\mathcal{S}^{n}_{\succ}$ and a suitable
factorization $\widehat{H}=E^{\mathsf{T}}Q$ with $E,Q\in\mathbb{R}^{n\times
n}$. Then the system

$\begin{bmatrix}E\dot{x}\\\ y\end{bmatrix}=\begin{bmatrix}A&B\\\
C&D\end{bmatrix}\begin{bmatrix}Qx\\\
u\end{bmatrix},\qquad\begin{bmatrix}A&B\\\
-C&-D\end{bmatrix}=\mathcal{J}-\mathcal{R},$ (2.5)

with $\mathcal{J}=-\mathcal{J}^{\top}\in\mathbb{R}^{(n+m)\times(n+m)}$ and
$\mathcal{R}\in\mathcal{S}^{n+m}_{\succeq}$ is called a _pH system_.

In Definition 2.2, the matrix $\mathcal{J}$ represents the conservation of
energy, while $\mathcal{R}$ describes dissipation. Note that the specific
structure requires the input and output dimension to match and directly
implies that pH systems are stable and passive, and the Hamiltonian serves as
a Lyapunov function. In more detail, any pH system together with its
Hamiltonian satisfies the dissipation inequality (1.3). Conversely, any
passive linear time-invariant system, i.e., any system for which a quadratic
Hamiltonian exists such that the system with this Hamiltonian satisfies the
dissipation inequality (1.3), has a pH representation, cf. [5]. Every pH
system (2.5) can be equivalently written as

$\displaystyle E\dot{x}$ $\displaystyle=(J-R)Qx+(G-P)u,$ (2.6a) $\displaystyle
y$ $\displaystyle=(G+P)^{\mathsf{T}}Qx+(S-N)u,$ (2.6b)

with $J=-J^{\top}$, $N=-N^{\top}$, and

$\mathcal{R}=\begin{bmatrix}R&P\\\
P^{\mathsf{T}}&S\end{bmatrix}\in\mathcal{S}^{n+m}_{\succeq}.$

For further details on the properties that are encoded within this structure
we refer to [6]. Multiplication of (2.6a) from the left with $Q^{\mathsf{T}}$,
see [31, Sec. 4.3], and introducing the matrices

$\displaystyle\widehat{J}$ $\displaystyle\vcentcolon=Q^{\mathsf{T}}JQ,$
$\displaystyle\widehat{R}$ $\displaystyle\vcentcolon=Q^{\mathsf{T}}RQ,$
$\displaystyle\widehat{G}$ $\displaystyle\vcentcolon=Q^{\mathsf{T}}G,$
$\displaystyle\widehat{P}$ $\displaystyle\vcentcolon=Q^{\mathsf{T}}P,$
$\displaystyle\widehat{S}$ $\displaystyle\vcentcolon=S,$
$\displaystyle\widehat{N}$ $\displaystyle\vcentcolon=N$

allows us to rewrite (2.6) as

$\displaystyle\widehat{H}\dot{x}$
$\displaystyle=(\widehat{J}-\widehat{R})x+(\widehat{G}-\widehat{P})u,$ (2.7a)
$\displaystyle y$
$\displaystyle=(\widehat{G}+\widehat{P})^{\mathsf{T}}x+(\widehat{S}-\widehat{N})u,$
(2.7b)

which is linear in the system matrices. Due to the congruence transformation,
we immediately conclude $\widehat{J}=-\widehat{J}^{\mathsf{T}}$,
$\widehat{N}=-\widehat{N}^{\mathsf{T}}$ and

$\widehat{\mathcal{R}}=\begin{bmatrix}\widehat{R}&\widehat{P}\\\
\widehat{P}^{\mathsf{T}}&\widehat{S}\end{bmatrix}\in\mathcal{S}^{n+m}_{\succeq}.$

######  Remark 2.3.

Let us emphasize that from a mathematical perspective, the choice of the
energy functional yielding to a pH system is generally not unique. Indeed, any
solution of the corresponding Kalman-Yakubovich-Popov inequality can be used
as the Hessian of the energy; see [5] for further details. Furthermore, recent
results detail that the choice of the energy functional characterizes how easy
it is to approximate the system [12].

To understand pH systems in the framework of DMD, we need to find a discrete-
time analogue of (2.7). Since symplectic Gauss-Legendre collocation methods
are able to preserve the underlying Dirac structure of the pH system [27, 30],
and DMD is able to recover the original dynamics if the DMD approximation is
based on a discretization with a 1-stage Runge-Kutta method [20], we use the
implicit midpoint rule for the time-discretization. More precisely, for a
constant step size $\delta_{t}>0$, the continuous dynamics (2.7a) are replaced
with the discrete dynamics

$\widehat{H}\tfrac{x_{i+1}-x_{i}}{\delta_{t}}=(\widehat{J}-\widehat{R})\tfrac{x_{i+1}+x_{i}}{2}+(\widehat{G}-\widehat{P})\tfrac{u_{i+1}+u_{i}}{2},$

with $x_{i}\approx x(i\delta_{t})$, and $u_{i}\vcentcolon=u(i\delta_{t})$.
Note that with respect to our goal to generate a discrete-time pH system from
data, we have replaced $u((i+\tfrac{1}{2})\delta_{t})$ from the implicit
midpoint rule with the approximation

$u((i+\tfrac{1}{2})\delta_{t})\approx\tfrac{u_{i+1}+u_{i}}{2}.$

It is natural to also replace the continuous output equation (2.6b) with its
discrete counterpart using $y_{i}\approx y(i\delta_{t})$. To preserve port-
Hamiltonian structure, we take the average of consecutive output equations
arriving at

$\tfrac{y_{i+1}+y_{i}}{2}=(\widehat{G}+\widehat{P})^{\mathsf{T}}\tfrac{x_{i+1}+x_{i}}{2}+(\widehat{S}-\widehat{N})\tfrac{u_{i+1}+u_{i}}{2}.$

The previous analysis motivates the following definition.

######  Definition 2.4 (discrete-time pH system).

Consider sequences of states $x_{i}\in\mathbb{R}^{n}$, inputs
$u_{i}\in\mathbb{R}^{m}$, outputs $y_{i}\in\mathbb{R}^{m}$, and a constant
time step $\delta_{t}>0$. We call a system of the form

$\displaystyle\widehat{H}\tfrac{x_{i+1}-x_{i}}{\delta_{t}}$
$\displaystyle=(\widehat{J}-\widehat{R})\tfrac{x_{i+1}+x_{i}}{2}+(\widehat{G}-\widehat{P})\tfrac{u_{i+1}+u_{i}}{2},$
(2.8a) $\displaystyle\tfrac{y_{i+1}+y_{i}}{2}$
$\displaystyle=(\widehat{G}+\widehat{P})^{\mathsf{T}}\tfrac{x_{i+1}+x_{i}}{2}+(\widehat{S}-\widehat{N})\tfrac{u_{i+1}+u_{i}}{2},$
(2.8b)

a _discrete-time pH system_ with quadratic Hamiltonian
$\mathcal{H}(x)\vcentcolon=\tfrac{1}{2}x^{\mathsf{T}}\widehat{H}x$ if and only
if $\widehat{J}=-\widehat{J}^{\mathsf{T}}$,
$\widehat{N}=-\widehat{N}^{\mathsf{T}}$, and
$\widehat{\mathcal{R}}=\left[\begin{smallmatrix}\widehat{R}&\widehat{P}\\\
\widehat{P}^{\mathsf{T}}&\widehat{S}\end{smallmatrix}\right]\in\mathcal{S}^{n+m}_{\succeq}$.

######  Remark 2.5.

The discrete-time pH system (2.8) is not in the standard form of a discrete-
time dynamical system. Nevertheless, assuming $\delta_{t}$ to be sufficiently
small, defining the matrices $\widehat{A}\vcentcolon=\widehat{J}-\widehat{R}$,
$\widehat{B}\vcentcolon=\widehat{G}-\widehat{P}$,
$\widehat{C}\vcentcolon=(\widehat{G}+\widehat{P})^{\mathsf{T}}$, and
$\widehat{D}\vcentcolon=\widehat{S}-\widehat{N}$, and assuming a consistent
initialization of the output, i.e., $y_{0}=\widehat{C}x_{0}+\widehat{D}u_{0}$,
then (2.8) can be rewritten as

$\displaystyle x_{i+1}$
$\displaystyle=(\tfrac{1}{\delta_{t}}\widehat{H}-\tfrac{1}{2}\widehat{A})^{-1}(\tfrac{1}{\delta_{t}}\widehat{H}+\tfrac{1}{2}\widehat{A})x_{i}+\tfrac{1}{2}(\tfrac{1}{\delta_{t}}\widehat{H}-\tfrac{1}{2}\widehat{A})^{-1}\widehat{B}(u_{i+1}+u_{i}),$
(2.9) $\displaystyle y_{i}$ $\displaystyle=\widehat{C}x_{i}+\widehat{D}u_{i}.$

While the equation (2.9) can be useful, since it represents $x_{i+1}$
explicitly in terms of $x_{i}$ and of the output variables, the pH structure
of the original system is not evident. Because of that, we prefer working with
the system (2.8) instead.

Let us emphasize that (2.8) can be equivalently written as

$\begin{bmatrix}\widehat{H}\tfrac{x_{i+1}-x_{i}}{\delta_{t}}\\\
-\tfrac{y_{i+1}+y_{i}}{2}\end{bmatrix}=(\widehat{\mathcal{J}}-\widehat{\mathcal{R}})\begin{bmatrix}\tfrac{x_{i+1}+x_{i}}{2}\\\
\tfrac{u_{i+1}+u_{i}}{2}\end{bmatrix},$ (2.10)

where $\widehat{\mathcal{J}}=-\widehat{\mathcal{J}}^{\mathsf{T}}$ and
$\widehat{\mathcal{R}}\in\mathcal{S}^{n+m}_{\succ}$. Defining
$\Delta\mathcal{H}_{i}\vcentcolon=\tfrac{1}{\delta_{t}}(\mathcal{H}(x_{i+1})-\mathcal{H}(x_{i}))$
and observing

$\Delta\mathcal{H}_{i}=(\tfrac{x_{i+1}+x_{i}}{2})^{\mathsf{T}}\widehat{H}(\tfrac{x_{i+1}-x_{i}}{\delta_{t}})$

immediately yields the following discrete version of the dissipation
inequality (1.3).

######  Lemma 2.6.

(Discrete-time dissipation inequality) Any discrete-time pH system in the form
of (2.8) satisfies the discrete-time dissipation inequality

$\Delta\mathcal{H}_{i}\leq(\tfrac{y_{i+1}+y_{i}}{2})^{\mathsf{T}}(\tfrac{u_{i+1}+u_{i}}{2}).$
(2.11)

######  Remark 2.7.

To simplify our presentation, we use a constant time step $\delta_{t}$
throughout this manuscript. Nevertheless, it is straightforward to use a
variable time step in all what follows.

##  3\. Port-Hamiltonian Dynamic Mode Decomposition

In this section, we discuss a variant of DMD that allows to construct a
continuous-time pH system from discrete-time data.

###  3.1. The port-Hamiltonian Dynamic Mode Decomposition problem

If the data $(x_{i},u_{i},y_{i})$ at hand are obtained from a physical system,
we expect the real system to have a pH representation and thus want to reflect
that in our time-discrete realization. Following our concise definition of a
discrete-time pH system (2.10), we are thus interested in solving the
following problem.

######  Problem 3.1 (Port-Hamiltonian Dynamic Mode Decomposition).

Consider data points
$(x_{i},u_{i},y_{i})\in\mathbb{R}^{n}\times\mathbb{R}^{m}\times\mathbb{R}^{m}$
for $i=0,1,\ldots,M$ and define the matrices

$\displaystyle\dot{\hat{X}}$
$\displaystyle\vcentcolon=\tfrac{1}{\delta_{t}}\begin{bmatrix}x_{1}-x_{0}&\ldots&x_{M}-x_{M-1}\end{bmatrix}\in\mathbb{R}^{n\times
M},$ $\displaystyle\hat{X}$
$\displaystyle\vcentcolon=\tfrac{1}{2}\begin{bmatrix}x_{1}+x_{0}&\ldots&x_{M}+x_{M-1}\end{bmatrix}\in\mathbb{R}^{n\times
M},$ $\displaystyle\hat{U}$
$\displaystyle\vcentcolon=\tfrac{1}{2}\begin{bmatrix}u_{1}+u_{0}&\ldots&u_{M}+u_{M-1}\end{bmatrix}\in\mathbb{R}^{m\times
M},$ $\displaystyle\hat{Y}$
$\displaystyle\vcentcolon=\tfrac{1}{2}\begin{bmatrix}y_{1}+y_{0}&\ldots&y_{M}+y_{M-1}\end{bmatrix}\in\mathbb{R}^{m\times
M}.$

Given any matrix $\widehat{H}\in\mathcal{S}^{n}_{\succ}$ (defining the
Hamiltonian), and a reduced basis $\Phi\in\mathbb{R}^{n\times\tilde{r}}$, find
matrices
$\widetilde{\mathcal{J}},\widetilde{\mathcal{R}}\in\mathbb{R}^{\tilde{n}\times\tilde{n}}$
with $\tilde{n}\vcentcolon=\tilde{r}+m$ that solve

$\min\left\|\widetilde{\mathcal{Z}}-(\widetilde{\mathcal{J}}-\widetilde{\mathcal{R}})\widetilde{\mathcal{T}}\right\|_{\mathrm{F}}\qquad\text{such
that
$\widetilde{\mathcal{J}}=-\widetilde{\mathcal{J}}^{\top}\in\mathbb{R}^{\tilde{n}\times\tilde{n}}$
and $\widetilde{\mathcal{R}}\in\mathcal{S}^{\tilde{n}}_{\succeq}$},$ (3.1)

where
$\widetilde{\mathcal{Z}}\vcentcolon=\begin{bmatrix}\Phi^{\mathsf{T}}H\Phi\Phi^{\mathsf{T}}\dot{\hat{X}}\\\
-\hat{Y}\end{bmatrix}$ and
$\widetilde{\mathcal{T}}\vcentcolon=\begin{bmatrix}\Phi^{\mathsf{T}}\hat{X}\\\
\hat{U}\end{bmatrix}$. If no dimension reduction is applied, then we set
$\tilde{r}=n$ and $\Phi=I_{n}$.

######  Remark 3.2.

Note that the structural properties of $\widetilde{\mathcal{J}}$ and
$\widetilde{\mathcal{R}}$ are preserved under the transformation with $\Phi$,
since it is a congruence transformation.

Using standard arguments, it is easy to establish that the minimization
problem (3.1) is convex and solvable. Moreover, let us emphasize that by
solving 3.1, we directly learn the matrices of a continuous-time pH system
without requiring further postprocessing. Moreover, 3.1 includes two important
special cases. First, suppose we know a-priori that no dissipation is
involved. In that case, we can set $\widetilde{\mathcal{R}}=0$, and the
minimization problem (3.1) reduces to a skew-symmetric Procrustes problem,
which can be solved analytically. We refer to [15, 4] and the forthcoming
Section 3.2. Second, if the system has no input and output, i.e., we have
$m=0$, then the task reduces to the identification of a dissipative
Hamiltonian system, i.e., we want to solve

$\min\left\|\widehat{H}\dot{\hat{X}}-\begin{bmatrix}\widetilde{J}-\widetilde{R}\end{bmatrix}\hat{X}\right\|_{\mathrm{F}}\qquad\text{such
that
$\widetilde{J}=-\widetilde{J}^{\top}\in\mathbb{R}^{\tilde{n}\times\tilde{n}}$
and $\widetilde{R}\in\mathcal{S}^{\tilde{n}}_{\succeq}$}.$ (3.2)

######  Remark 3.3.

If instead of prescribing the Hamiltonian and thus the matrix $\widehat{H}$,
one also wants to determine $\widehat{H}\in\mathcal{S}^{n}_{\succ}$ from the
data, then we observe that (3.2) is not a suitable formulation, since in this
case (3.2) is not solvable. To see this, notice that $\mathcal{S}^{n}_{\succ}$
is an open set. Hence, for any $\varepsilon>0$, set
$\widetilde{J}=\widetilde{R}=0\in\mathbb{R}^{n\times n}$ and
$\widehat{H}=\smash{\tfrac{\varepsilon}{2\|\dot{\hat{X}}\|_{\mathrm{F}}}I_{n}}$.
Then $\widehat{H}\in\mathcal{S}^{n}_{\succ}$,
$\smash{\widetilde{J}}=-\smash{\widetilde{J}^{\mathsf{T}}}$,
$\smash{\widetilde{R}}\in\mathcal{S}^{n}_{\succeq}$ and we have

$\left\|\widehat{H}\dot{\hat{X}}-\begin{bmatrix}\widetilde{J}-\widetilde{R}\end{bmatrix}\hat{X}\right\|_{\mathrm{F}}=\bigg{\|}\tfrac{\varepsilon}{2\|\dot{\hat{X}}\|_{\mathrm{F}}}\dot{\hat{X}}\bigg{\|}_{\mathrm{F}}=\tfrac{\varepsilon}{2}<\varepsilon,$

implying that the infimum of (3.2) is zero with infimizer
$H=\widetilde{J}=\widetilde{R}=0$.

###  3.2. An iterative algorithm for the pHDMD problem

Having established the existence of a solution for 3.1, we will now derive a
numerical algorithm to solve the minimization problem (3.1). In [16], a _fast
gradient method_ (FGM), cf. [33, p. 90], for the nearest stable matrix to a
given matrix is proposed, which is formulated similarly as 3.1. The FGM is an
optimal first-order method for convex optimization, which means no other
first-order method can converge faster while using the same first-order
information. However, we notice that if we already have a guess for
$\widetilde{\mathcal{R}}$, then 3.1 simplifies to the skew-symmetric
Procrustes problem

$\min\left\|\mathcal{Z}_{1}-\widetilde{\mathcal{J}}\mathcal{T}\right\|_{\mathrm{F}}\qquad\text{such
that
$\widetilde{\mathcal{J}}=-\widetilde{\mathcal{J}}^{\top}\in\mathbb{R}^{\tilde{n}\times\tilde{n}}$},$
(3.3)

where
$\mathcal{Z}_{1}\vcentcolon=\mathcal{Z}+\widetilde{\mathcal{R}}\mathcal{T}$.
Fortunately, the solution of (3.3) can be computed analytically as detailed in
the following theorem taken as a special case of [15, Lem. 2.1], so we propose
to include this knowledge into the iterative algorithm to achieve
acceleration.

######  Theorem 3.4.

Let $\mathcal{T},\mathcal{Z}_{1}\in\mathbb{R}^{\tilde{n}\times M}$ and let
$V\Sigma W^{\mathsf{T}}=\mathcal{T}$ denote the SVD of $\mathcal{T}$ with

$\Sigma=\begin{bmatrix}\Sigma_{1}&0\\\ 0&0\end{bmatrix},\qquad\text{with}\
\Sigma_{1}=\operatorname{diag}(\sigma_{1},\ldots,\sigma_{r})\in\mathbb{R}^{r\times
r}\ \text{and}\ r=\operatorname{rank}(T).$

Define $\Phi=[\phi_{ij}]\in\mathbb{R}^{r\times r}$ via
$\phi_{ij}=\frac{1}{\sigma_{i}^{2}+\sigma_{j}^{2}}\qquad\text{for}\
i,j=1,\ldots,r$. Then

$\widetilde{\mathcal{J}}=V^{\mathsf{T}}\begin{bmatrix}\Phi\odot\left(2\operatorname{skew}(Z_{1}\Sigma_{1})\right)&-\Sigma_{1}^{-1}Z_{3}^{\mathsf{T}}\\\
Z_{3}\Sigma_{1}^{-1}&\widetilde{\mathcal{J}}_{4}\end{bmatrix}V$ (3.4)

is a solution of (3.3) for any skew-symmetric matrix
$\widetilde{\mathcal{J}}_{4}\in\mathbb{R}^{(\tilde{n}-r)\times(\tilde{n}-r)}$,
where

$\displaystyle Z_{1}$
$\displaystyle=\begin{bmatrix}I_{r}&0\end{bmatrix}V\mathcal{Z}_{1}W^{\mathsf{T}}\begin{bmatrix}I_{r}&0\end{bmatrix}^{\mathsf{T}}\in\mathbb{R}^{r\times
r},$ $\displaystyle Z_{3}$
$\displaystyle=\begin{bmatrix}0&I_{n-r}\end{bmatrix}V\mathcal{Z}_{1}W^{\mathsf{T}}\begin{bmatrix}I_{r}&0\end{bmatrix}^{\mathsf{T}}\in\mathbb{R}^{(\tilde{n}-r)\times
r}.$

The solution is unique if, and only if, $\operatorname{rank}(T)=\tilde{n}$.

###### Proof.

The proof follows along the proof of [15, Lem. 2.1]. ∎

On the other hand if $\widetilde{\mathcal{J}}$ is given, 3.1 simplifies to a
symmetric positive definite Procrustes problem

$\min\left\|\mathcal{Z}_{2}-\widetilde{\mathcal{R}}\mathcal{T}\right\|_{\mathrm{F}}\qquad\text{s.t.
$\widetilde{\mathcal{R}}\in\mathcal{S}^{\tilde{n}}_{\succeq}$},$ (3.5)

where $\mathcal{Z}_{2}=\widetilde{\mathcal{J}}\mathcal{T}-\mathcal{Z}$.
Algorithmic solutions for this problem are available [17], for instance a FGM.
By modifying this algorithm to allow a $\mathcal{Z}_{2}$ depending on the
optimal solution $\widetilde{\mathcal{J}}$ of the skew-symmetric Procrustes
problem (3.3), we arrive at Algorithm 1.

Algorithm 1 Semi-analytical fast gradient method for (3.1)

1:Input: Data matrices
$\mathcal{Z},\mathcal{T}\in\mathbb{R}^{\tilde{n}\times\tilde{M}}$, initial
guess for the dissipative part
$\widetilde{\mathcal{R}}^{(0)}\in\mathcal{S}^{n}_{\succeq}$

2:Output: Matrices
$\widetilde{\mathcal{J}}=-\widetilde{\mathcal{J}}^{\mathsf{T}}\in\mathbb{R}^{\tilde{n}\times\tilde{n}}$
and $\widetilde{\mathcal{R}}\in\mathcal{S}^{\tilde{n}}_{\succeq}$ that
minimize (3.1)

3:

4:$L=\sigma^{2}_{1}(\mathcal{T})$, $q=\frac{\sigma^{2}_{r}}{L}$;

5:$k=0$; $\alpha_{1}\in(0,1)$;

6:$Q=\widetilde{\mathcal{R}}^{(0)}$;

7:while not converged do

8: $\mathcal{Z}_{1}=\mathcal{Z}+\widetilde{\mathcal{R}}^{(k)}\mathcal{T}$

9: Solve (3.3) for $\mathcal{Z}_{1},\mathcal{T}$ according to Theorem 3.4 to
obtain $\widetilde{\mathcal{J}}^{(k+1)}$

10: $\mathcal{Z}_{2}=\widetilde{\mathcal{J}}^{(k+1)}\mathcal{T}-\mathcal{Z}$

11: $\nabla=Q\mathcal{T}\mathcal{T}^{\mathsf{T}}-\mathcal{Z}_{2}\mathcal{T}$

12:
$\widetilde{\mathcal{R}}^{(k+1)}=\mathcal{P}_{\succeq}(Q-\frac{1}{L}\nabla)$

13:
$\alpha_{k+1}=\tfrac{1}{2}(q-\alpha_{k}^{2}+\sqrt{(q-\alpha_{k}^{2})^{2}+4\alpha_{k}^{2}}),\quad\beta_{k}=\frac{\alpha_{k}(1-\alpha_{k})}{\alpha_{k}^{2}+\alpha_{k+1}}$

14:
$Q=\widetilde{\mathcal{R}}^{(k+1)}+\beta_{k}(\widetilde{\mathcal{R}}^{(k+1)}-\widetilde{\mathcal{R}}^{(k)})$

15: $k=k+1$

16:end while

17:$\widetilde{\mathcal{J}}=\widetilde{\mathcal{J}}^{(k)}$,
$\widetilde{\mathcal{R}}=\widetilde{\mathcal{R}}^{(k)}$

In more detail, within each iteration step of Algorithm 1, we first compute
the solution of the skew-symmetric Procrustes problem (8 and 9) and then
compute the gradient with respect to the matrix $\widetilde{\mathcal{R}}$ in
10 and 11, ignoring for the moment that we need the new iterate to be
symmetric positive definite. This is achieved by projecting onto the cone of
symmetric positive definite matrices in 12. The update is then computed as a
linear combination of the current and previous iterate with the fast gradient
coefficients (cf. 13 and 14).

######  Remark 3.5.

In general, we cannot guarantee convergence of Algorithm 1. Nevertheless, we
can use the standard safety strategy for the fast-gradient, as, for instance,
reported in [16], by using a reinitialization with a standard gradient step
and a backtracking line search. By doing so, classical convergence results can
be obtained since our objective function is convex and solvable.

We notice in our numerical experiments that the performance of the algorithm
strongly depends on the initialization and may need many iterations to
converge if a poor initialization is used. We thus study a particular
initialization strategy in the next subsection by analyzing 3.1 in a weighted
norm.

###  3.3. A weighted pHDMD problem

A different but related problem arises when the Frobenius norm of the 3.1 is
replaced with the weighted Frobenius seminorm introduced in (1.4), where
$\mathcal{T}^{\mathsf{T}}\mathcal{T}$ is used as the semi-definite weighting
matrix.

######  Problem 3.6 (Weighted Input-Output port-Hamiltonian Dynamic Mode
Decomposition).

For given data
$\mathcal{Z},\mathcal{T}\in\mathbb{R}^{\tilde{n}\times\tilde{M}}$, solve the
optimization problem

$\min\left\|\mathcal{T}^{\mathsf{T}}(\mathcal{Z}-(\widetilde{\mathcal{J}}-\widetilde{\mathcal{R}})\mathcal{T})\right\|_{\mathrm{F}}\qquad\text{such
that
$\widetilde{\mathcal{J}}=-\widetilde{\mathcal{J}}^{\top}\in\mathbb{R}^{\tilde{n}\times\tilde{n}}$
and $\widetilde{\mathcal{R}}\in\mathcal{S}^{\tilde{n}}_{\succeq}$}.$ (3.6)

Two remarks are in order, first we immediately notice that by solving the
weighted problem we solve the original problem, but weighted according to the
relevant information in the data. Second, using the definition of the data
matrices from 3.1, we observe that the $(i,i)$ entry of the matrix
$\mathcal{T}^{\mathsf{T}}(\mathcal{Z}-(\widetilde{\mathcal{J}}-\widetilde{\mathcal{R}})\mathcal{T})$
reads

$\displaystyle
e_{i}^{\mathsf{T}}(\mathcal{T}^{\mathsf{T}}(\mathcal{Z}-(\widetilde{\mathcal{J}}-\widetilde{\mathcal{R}})\mathcal{T})e_{i}$
$\displaystyle=e_{i}^{\mathsf{T}}(\hat{X})^{\mathsf{T}}\widehat{H}\dot{\hat{X}}e_{i}-e_{i}^{\mathsf{T}}\hat{U}^{\mathsf{T}}\hat{Y}e_{i}+e_{i}^{\mathsf{T}}\begin{bmatrix}(\hat{X})^{\mathsf{T}}&\hat{U}^{\mathsf{T}}\end{bmatrix}\widetilde{\mathcal{R}}\begin{bmatrix}\hat{X}\\\
\hat{U}\end{bmatrix}e_{i}$
$\displaystyle=\Delta\mathcal{H}_{i}-(\tfrac{y_{i+1}+y_{i}}{2})^{\mathsf{T}}(\tfrac{u_{i+1}+u_{i}}{2})+e_{i}^{\mathsf{T}}\begin{bmatrix}(\hat{X})^{\mathsf{T}}&\hat{U}^{\mathsf{T}}\end{bmatrix}\widetilde{\mathcal{R}}\begin{bmatrix}\hat{X}\\\
\hat{U}\end{bmatrix}e_{i},$

which resembles the power-balance equation corresponding to the discrete-time
dissipation inequality (2.11). In particular, the weighted problem (3.6) uses
the dissipation inequality as part of the cost-functional to determine the
dissipative component, for which we need a good initialization.

######  Theorem 3.7.

Let $\mathcal{T},\mathcal{Z}\in\mathbb{R}^{\tilde{n}\times\tilde{M}}$ and let
$V_{1}\Sigma_{1}W_{1}^{\mathsf{T}}=\mathcal{T}$ denote the skinny SVD of
$\mathcal{T}$, i.e., $V_{1}\in\mathrm{St}(\tilde{n},r)$,
$\Sigma_{1}\in\mathrm{GL}_{r}(\mathbb{R})$, and
$W_{1}\in\mathrm{St}(\tilde{M},r)$, where
$r\vcentcolon=\operatorname{rank}(\mathcal{T})$. Moreover, let $W_{2}$
complement $W_{1}$ to an orthogonal matrix, i.e.,
$\begin{bmatrix}W_{1}&W_{2}\end{bmatrix}\in\mathrm{St}(\tilde{M},\tilde{M})$.
Define

$\displaystyle\widetilde{\mathcal{Z}}_{1}$
$\displaystyle\vcentcolon=\Sigma_{1}V_{1}^{\mathsf{T}}\mathcal{Z}W_{1},$
$\displaystyle\widetilde{\mathcal{Z}}_{2}$
$\displaystyle\vcentcolon=\Sigma_{1}V_{1}^{\mathsf{T}}\mathcal{Z}W_{2},$

and

$\displaystyle\widetilde{\mathcal{J}}^{\star}\vcentcolon=V_{1}\Sigma_{1}^{-1}\operatorname{skew}(\widetilde{\mathcal{Z}}_{1})\Sigma_{1}^{-1}V_{1}^{\mathsf{T}}\qquad\text{and}\qquad\widetilde{\mathcal{R}}^{\star}\vcentcolon=V_{1}\Sigma_{1}^{-1}\mathcal{P}_{\succeq}(-\widetilde{\mathcal{Z}}_{1})\Sigma_{1}^{-1}V_{1}^{\mathsf{T}}.$
(3.7)

Then, $\widetilde{\mathcal{J}}^{\star}$ and $\widetilde{\mathcal{R}}^{\star}$
are the unique minimum-norm minimizers of (3.6) with

$\left\|\mathcal{T}^{\mathsf{T}}\mathcal{Z}-\mathcal{T}^{\mathsf{T}}(\widetilde{\mathcal{J}}^{\star}-\widetilde{\mathcal{R}}^{\star})\mathcal{T}\right\|_{\mathrm{F}}^{2}=\left\|\widetilde{\mathcal{Z}}_{2}\right\|_{\mathrm{F}}^{2}+\left\|\Lambda_{+}\right\|_{\mathrm{F}}^{2},$
(3.8)

where $\Lambda_{+}$ is the diagonal matrix which contains the positive
eigenvalues of $\operatorname{sym}(\widetilde{\mathcal{Z}}_{1})$

######  Remark 3.8.

The two terms on the right-hand side of (3.8) can be interpreted as follows.
The first term is only present if we have too much data to fit, i.e., if
$\mathcal{T}$ has more columns than rows. Then, in general, no linear system
can perfectly capture the data, and the corresponding error contribution is
$\left\|\smash{\widetilde{\mathcal{Z}}_{2}}\right\|_{\mathrm{F}}^{2}$. On the
other hand, due to the specific pH structure, not every linear system can be
written as a pH system. This potential deviation corresponds to the second
error term given by $\left\|\Lambda_{+}\right\|_{\mathrm{F}}^{2}$.

Before we present the proof, we first illustrate Theorem 3.7 with an academic
toy example and need some further preliminary results that provide the best-
fit in the Frobenius norm within the class of skew-symmetric and symmetric
positive semi-definite matrices.

######  Example 3.9.

Consider $\mathcal{Z}=\left[\begin{smallmatrix}-1&2\\\
\phantom{-}2&-1/2\end{smallmatrix}\right]$ and
$\mathcal{T}=\left[\begin{smallmatrix}1&0\\\ 0&2\end{smallmatrix}\right]$. We
immediately notice that we can choose $V_{1}=W_{1}=I_{2}$ and
$\Sigma_{1}=\mathcal{T}$. We thus obtain
$\widetilde{\mathcal{Z}}_{1}=\Sigma_{1}V_{1}^{\mathsf{T}}\mathcal{Z}W_{1}=\left[\begin{smallmatrix}-1&\phantom{-}2\\\
\phantom{-}4&-1\end{smallmatrix}\right]$ and thus

$\operatorname{sym}(\widetilde{\mathcal{Z}}_{1})=\begin{bmatrix}-1&\phantom{-}3\\\
\phantom{-}3&-1\end{bmatrix}\qquad\text{and}\qquad\operatorname{skew}(\widetilde{\mathcal{Z}}_{1})=\begin{bmatrix}0&-1\\\
1&\phantom{-}0\end{bmatrix}.$

Then, the ordered eigendecomposition of
$\widetilde{\mathcal{Z}}_{\mathrm{sym}}$ is given by,

$\widetilde{\mathcal{Z}}_{\mathrm{sym}}=\Xi\Lambda\Xi^{\mathsf{T}}=\begin{bmatrix}\frac{\sqrt{2}}{2}&\phantom{-}\frac{\sqrt{2}}{2}\\\
\frac{\sqrt{2}}{2}&-\frac{\sqrt{2}}{2}\end{bmatrix}\begin{bmatrix}2&\phantom{-}0\\\
0&-4\end{bmatrix}\begin{bmatrix}\frac{\sqrt{2}}{2}&\phantom{-}\frac{\sqrt{2}}{2}\\\
\frac{\sqrt{2}}{2}&-\frac{\sqrt{2}}{2}\end{bmatrix}.$

Thus, the unique minimizers are given by

$\widetilde{\mathcal{J}}^{\star}=\begin{bmatrix}0&-\frac{1}{2}\\\
\frac{1}{2}&\phantom{-}0\end{bmatrix}\qquad\text{and}\qquad\widetilde{\mathcal{R}}^{\star}=\begin{bmatrix}\phantom{-}2&-1\\\
-1&\phantom{-}\tfrac{1}{2}\end{bmatrix}.$ (3.9)

The minimal value of the optimization is given by

$\left\|\mathcal{T}^{\mathsf{T}}\mathcal{Z}-\mathcal{T}^{\mathsf{T}}(\widetilde{\mathcal{J}}^{\star}-\widetilde{\mathcal{R}}^{\star})\mathcal{T}\right\|_{\mathrm{F}}=\left\|\begin{bmatrix}1&1\\\
1&1\end{bmatrix}\right\|_{\mathrm{F}}=2=\left\|\Lambda_{+}\right\|_{\mathrm{F}}.$

For the proof of Theorem 3.7 we need the following technical result, taken
from [16, Lem. 7], see also [21].

######  Lemma 3.10.

Let $Z\in\mathbb{R}^{r\times r}$, then the minimization problem

$\min\left\|Z-(J-R)\right\|_{\mathrm{F}}\qquad\text{such that
$J=-J^{\mathsf{T}}$ and $R\in\mathcal{S}^{r}_{\succeq}$}.$ (3.10)

is solved by $J^{\star}\vcentcolon=\operatorname{skew}(Z)$ and
$R^{\star}\vcentcolon=\mathcal{P}_{\succeq}(-Z)$.

###### Proof of Theorem 3.7.

Let $\mathcal{T}=V\Sigma W^{\mathsf{T}}$ denote the singular value
decomposition of $\mathcal{T}$, with partitioning

$\displaystyle V$ $\displaystyle=\begin{bmatrix}V_{1}&V_{2}\end{bmatrix},$
$\displaystyle\Sigma$ $\displaystyle=\operatorname{diag}(\Sigma_{1},0),$
$\displaystyle W$ $\displaystyle=\begin{bmatrix}W_{1}&W_{2}\end{bmatrix}$
(3.11a) such that $\Sigma_{1}\in\mathrm{GL}_{r}(\mathbb{R})$, where
$r=\operatorname{rank}(\mathcal{T})$. In particular, we obtain
$\mathcal{T}=V_{1}\Sigma_{1}W_{1}^{\mathsf{T}}$. Let
$\widetilde{\mathcal{J}}=-\widetilde{\mathcal{J}}^{\mathsf{T}}\in\mathbb{R}^{\tilde{n}\times\tilde{n}}$
and $\widetilde{\mathcal{R}}\in\mathcal{S}^{\tilde{n}}_{\succeq}$. Define
$\displaystyle\begin{bmatrix}\widetilde{\mathcal{J}}_{11}&-\widetilde{\mathcal{J}}_{21}^{\mathsf{T}}\\\
\widetilde{\mathcal{J}}_{21}&\phantom{-}\widetilde{\mathcal{J}}_{22}\end{bmatrix}$
$\displaystyle\vcentcolon=V^{\mathsf{T}}\widetilde{\mathcal{J}}V,$
$\displaystyle\begin{bmatrix}\widetilde{\mathcal{R}}_{11}&\widetilde{\mathcal{R}}_{21}^{\mathsf{T}}\\\
\widetilde{\mathcal{R}}_{21}&\widetilde{\mathcal{R}}_{22}\end{bmatrix}$
$\displaystyle\vcentcolon=V^{\mathsf{T}}\widetilde{\mathcal{R}}V$
$\displaystyle\begin{bmatrix}\mathcal{Z}_{11}&\mathcal{Z}_{12}\\\
\mathcal{Z}_{21}&\mathcal{Z}_{22}\end{bmatrix}$
$\displaystyle\vcentcolon=V^{\mathsf{T}}\mathcal{Z}W$ (3.11b)

Note that throughout the proof we work with congruence transformations of
$\widetilde{\mathcal{J}}$ and $\widetilde{\mathcal{R}}$, which preserve the
skew-symmetry and the positive semi-definiteness. We then obtain

$\displaystyle\left\|\mathcal{T}^{\mathsf{T}}\mathcal{Z}-\mathcal{T}^{\mathsf{T}}(\widetilde{\mathcal{J}}-\widetilde{\mathcal{R}})\mathcal{T}\right\|_{\mathrm{F}}^{2}$
$\displaystyle=\left\|\Sigma V^{\mathsf{T}}\mathcal{Z}W-\Sigma
V^{\mathsf{T}}(\widetilde{\mathcal{J}}-\widetilde{\mathcal{R}})V\Sigma\right\|_{\mathrm{F}}^{2}$
$\displaystyle=\left\|\begin{bmatrix}\Sigma_{1}\mathcal{Z}_{11}&\Sigma_{1}\mathcal{Z}_{12}\\\
0&0\end{bmatrix}-\begin{bmatrix}\Sigma_{1}(\widetilde{\mathcal{J}}_{11}-\widetilde{\mathcal{R}}_{11})\Sigma_{1}&0\\\
0&0\end{bmatrix}\right\|_{\mathrm{F}}^{2}$
$\displaystyle=\left\|\widetilde{\mathcal{Z}}_{1}-\Sigma_{1}(\widetilde{\mathcal{J}}_{11}-\widetilde{\mathcal{R}}_{11})\Sigma_{1}\right\|_{\mathrm{F}}^{2}+\left\|\widetilde{\mathcal{Z}}_{2}\right\|_{\mathrm{F}}^{2}.$

We immediately notice that
$\widetilde{\mathcal{J}}_{21},\widetilde{\mathcal{J}}_{22},\widetilde{\mathcal{R}}_{21}$,
and $\widetilde{\mathcal{R}}_{22}$ do not influence the objective function and
can thus be chosen arbitrarily (provided that
$\widetilde{\mathcal{J}}_{22}=-\widetilde{\mathcal{J}}_{22}^{\mathsf{T}}$ and
$\widetilde{\mathcal{R}}\in\mathcal{S}^{\tilde{n}}_{\succeq}$). For our
further construction we set them to zero, in agreement with (3.7). It thus
suffices to minimize over all skew-symmetric matrices
$\widetilde{\mathcal{J}}_{11}\in\mathbb{R}^{r\times r}$ and all
$\widetilde{\mathcal{R}}_{11}\in\mathcal{S}^{r}_{\succeq}$. Thus, using Lemma
3.10, we obtain

$\displaystyle\min_{\begin{subarray}{c}\widetilde{\mathcal{J}}=-\widetilde{\mathcal{J}}^{\mathsf{T}},\\\
\widetilde{\mathcal{R}}\in\mathcal{S}^{\tilde{n}}_{\succeq}\end{subarray}}\left\|\mathcal{T}^{\mathsf{T}}\mathcal{Z}-\mathcal{T}^{\mathsf{T}}(\widetilde{\mathcal{J}}-\widetilde{\mathcal{R}})\mathcal{T}\right\|_{\mathrm{F}}^{2}$
$\displaystyle=\min_{\begin{subarray}{c}\widetilde{\mathcal{J}}_{11}=-\widetilde{\mathcal{J}}_{11}^{\mathsf{T}},\\\
\widetilde{R}_{11}\in\mathcal{S}^{r}_{\succeq}\end{subarray}}\left\|\widetilde{\mathcal{Z}}_{1}-\Sigma_{1}(\widetilde{\mathcal{J}}_{11}-\widetilde{\mathcal{R}}_{11})\Sigma_{1}\right\|_{\mathrm{F}}^{2}+\left\|\widetilde{\mathcal{Z}}_{2}\right\|_{\mathrm{F}}^{2}$
$\displaystyle=\min_{\widetilde{\mathcal{R}}_{11}\in\mathcal{S}^{r}_{\succeq}}\left\|\widetilde{\mathcal{Z}}_{1}-\Sigma_{1}(\widetilde{\mathcal{J}}_{11}^{\star}-\widetilde{\mathcal{R}}_{11})\Sigma_{1}\right\|_{\mathrm{F}}^{2}+\left\|\widetilde{\mathcal{Z}}_{2}\right\|_{\mathrm{F}}^{2}$
$\displaystyle=\min_{\widetilde{\mathcal{R}}_{11}\in\mathcal{S}^{r}_{\succeq}}\left\|\operatorname{sym}(\widetilde{\mathcal{Z}}_{1})+\Sigma_{1}\widetilde{\mathcal{R}}_{11}\Sigma_{1}\right\|_{\mathrm{F}}^{2}+\left\|\widetilde{\mathcal{Z}}_{2}\right\|_{\mathrm{F}}^{2}$
$\displaystyle=\left\|\operatorname{sym}(\widetilde{\mathcal{Z}}_{1})+\Sigma_{1}\widetilde{\mathcal{R}}_{11}^{\star}\Sigma_{1}\right\|_{\mathrm{F}}+\left\|\widetilde{\mathcal{Z}}_{2}\right\|_{\mathrm{F}}^{2}=\left\|\Lambda_{+}\right\|_{\mathrm{F}}^{2}+\left\|\widetilde{\mathcal{Z}}_{2}\right\|_{\mathrm{F}}^{2}.$

It remains to show that $\widetilde{\mathcal{J}}^{\star}$ and
$\widetilde{\mathcal{R}}^{\star}$ are the minimizers with minimal norm. To
this end, let
$\widetilde{\mathcal{J}}=-\widetilde{\mathcal{J}}^{\mathsf{T}}\in\mathbb{R}^{\tilde{n}\times\tilde{n}}$
and $\widetilde{\mathcal{R}}\in\mathcal{S}^{\tilde{n}}_{\succeq}$ be further
minimizers of (3.6). Then

$\left\|\widetilde{\mathcal{J}}\right\|_{\mathrm{F}}=\left\|W^{\mathsf{T}}\widetilde{\mathcal{J}}W\right\|_{\mathrm{F}}=\left\|\begin{bmatrix}\widetilde{\mathcal{J}}_{11}&-\widetilde{\mathcal{J}}_{12}^{\mathsf{T}}\\\
\widetilde{\mathcal{J}}_{21}&\phantom{-}\widetilde{\mathcal{J}}_{22}\end{bmatrix}\right\|_{\mathrm{F}}.$

Lemma 3.10 implies
$\widetilde{\mathcal{J}}_{11}=\widetilde{\mathcal{J}}_{11}^{\star}$, and thus

$\left\|\widetilde{\mathcal{J}}\right\|_{\mathrm{F}}=\left\|\begin{bmatrix}\widetilde{\mathcal{J}}_{11}^{\star}&-\widetilde{\mathcal{J}}_{12}^{\mathsf{T}}\\\
\widetilde{\mathcal{J}}_{21}&\phantom{-}\widetilde{\mathcal{J}}_{22}\end{bmatrix}\right\|_{\mathrm{F}}\geq\left\|\begin{bmatrix}\widetilde{\mathcal{J}}_{11}^{\star}&0\\\
0&0\end{bmatrix}\right\|_{\mathrm{F}}=\left\|W\widetilde{\mathcal{J}}^{\star}W\right\|_{\mathrm{F}}=\left\|\widetilde{\mathcal{J}}^{\star}\right\|_{\mathrm{F}}.$

A similar argument shows
$\left\|\widetilde{\mathcal{R}}\right\|_{\mathrm{F}}\geq\left\|\widetilde{\mathcal{R}}^{\star}\right\|_{\mathrm{F}}$,
which completes the proof. ∎

###  3.4. Relation between the optimization problems

In this subsection we discuss the relation between the original pHDMD
optimization 3.1 and the weighted 3.6 discussed in the previous section. We
immediately obtain the following result, which showcases that whenever the
data is sufficiently rich, then the cost functional of the projected
optimization problem (3.6) provides an upper bound for the original
minimization problem (3.1).

######  Lemma 3.11.

For given $\mathcal{Z},\mathcal{T}\in\mathbb{R}^{\tilde{n}\times\tilde{M}}$
with $\operatorname{rank}(\mathcal{T})=\tilde{n}$ there exists a constant
$c>0$ such that for every
$\widetilde{\mathcal{J}},\widetilde{\mathcal{R}}\in\mathbb{R}^{\tilde{n}\times\tilde{n}}$
we have

$\left\|\mathcal{Z}-(\widetilde{\mathcal{J}}-\widetilde{\mathcal{R}})\mathcal{T}\right\|_{\mathrm{F}}\leq
c\left\|\mathcal{T}^{\mathsf{T}}\mathcal{Z}-\mathcal{T}^{\mathsf{T}}(\widetilde{\mathcal{J}}-\widetilde{\mathcal{R}})\mathcal{T}\right\|_{\mathrm{F}}.$
(3.12)

###### Proof.

Let $\mathcal{T}=V\Sigma W^{\mathsf{T}}$ denote the singular value
decomposition of $\mathcal{T}$. We then have
$\mathcal{T}^{\dagger}=W\Sigma^{\dagger}V^{\mathsf{T}}$. Define
$c\vcentcolon=\left\|\smash{\mathcal{T}^{\dagger}}\right\|_{\mathrm{F}}=\left\|\smash{\Sigma^{\dagger}}\right\|_{\mathrm{F}}$.
Using $\operatorname{rank}{\mathcal{T}}=\tilde{n}$, we conclude
$\Sigma^{\dagger}\Sigma=I_{\tilde{n}}$. Let
$\widetilde{\mathcal{J}},\widetilde{\mathcal{R}}\in\mathbb{R}^{\tilde{n}\times\tilde{n}}$.
Then

$\displaystyle\left\|\mathcal{Z}-(\widetilde{\mathcal{J}}-\widetilde{\mathcal{R}})\mathcal{T}\right\|_{\mathrm{F}}$
$\displaystyle=\left\|VV^{\mathsf{T}}(\mathcal{Z}-(\widetilde{\mathcal{J}}-\widetilde{\mathcal{R}})\mathcal{T})\right\|_{\mathrm{F}}=\left\|(\mathcal{T}^{\mathsf{T}})^{\dagger}\mathcal{T}^{\mathsf{T}}(\mathcal{Z}-(\widetilde{\mathcal{J}}-\widetilde{\mathcal{R}})\mathcal{T})\right\|_{\mathrm{F}}$
$\displaystyle\leq\left\|(\mathcal{T}^{\mathsf{T}})^{\dagger}\right\|_{\mathrm{F}}\left\|\mathcal{T}^{\mathsf{T}}\mathcal{Z}-\mathcal{T}^{\mathsf{T}}(\widetilde{\mathcal{J}}-\widetilde{\mathcal{R}})\mathcal{T}\right\|_{\mathrm{F}},$

which completes the proof. ∎

If $\mathcal{T}$ has not full row rank, i.e.,
$\operatorname{rank}(\mathcal{T})<\tilde{n}$, then we cannot expect to obtain
a similar result, in particular if we use the minimum norm-minimizers from
Theorem 3.7. The main reason for this behavior is that in this case, using the
notation as in (3.11), we project out the data corresponding to
$\mathcal{Z}_{21}$ and $\mathcal{Z}_{22}$. While the latter corresponds,
similarly as $\mathcal{Z}_{12}$ to too much data, the contribution
$\mathcal{Z}_{21}$ corresponds to the components
$\widetilde{\mathcal{J}}_{21}-\widetilde{\mathcal{R}}_{21}$, which are set
zero in Theorem 3.7. In more detail, using the notation as in (3.11), we
obtain for
$\widetilde{\mathcal{J}}=-\widetilde{\mathcal{J}}^{\mathsf{T}}\in\mathbb{R}^{\tilde{n}\times\tilde{n}}$
and $\widetilde{\mathcal{R}}\in\mathcal{S}^{\tilde{n}}_{\succeq}$

$\displaystyle\left\|\mathcal{Z}-(\widetilde{\mathcal{J}}-\widetilde{\mathcal{R}})\mathcal{T}\right\|_{\mathrm{F}}$
$\displaystyle=\left\|\begin{bmatrix}V_{1}^{\mathsf{T}}\\\
V_{2}^{\mathsf{T}}\end{bmatrix}\mathcal{Z}\begin{bmatrix}W_{1}&W_{2}\end{bmatrix}-\begin{bmatrix}V_{1}^{\mathsf{T}}\\\
V_{2}^{\mathsf{T}}\end{bmatrix}(\widetilde{\mathcal{J}}-\widetilde{\mathcal{R}})\begin{bmatrix}V_{1}&V_{2}\end{bmatrix}\begin{bmatrix}\Sigma_{1}&0\\\
0&0\end{bmatrix}\right\|_{\mathrm{F}}$
$\displaystyle=\left\|\begin{bmatrix}\mathcal{Z}_{11}&\mathcal{Z}_{12}\\\
\mathcal{Z}_{21}&\mathcal{Z}_{22}\end{bmatrix}-\begin{bmatrix}\widetilde{\mathcal{J}}_{11}-\widetilde{\mathcal{R}}_{11}&-\widetilde{\mathcal{J}}_{21}^{\mathsf{T}}-\widetilde{\mathcal{R}}_{12}^{\mathsf{T}}\\\
\widetilde{\mathcal{J}}_{21}-\widetilde{\mathcal{R}}_{21}&\phantom{-}\widetilde{\mathcal{J}}_{22}-\widetilde{\mathcal{R}}_{22}\end{bmatrix}\begin{bmatrix}\Sigma_{1}&0\\\
0&0\end{bmatrix}\right\|_{\mathrm{F}}$
$\displaystyle=\left\|\begin{bmatrix}\mathcal{Z}_{11}-(\widetilde{\mathcal{J}}_{11}-\widetilde{\mathcal{R}}_{11})\Sigma_{1}&\mathcal{Z}_{12}\\\
\mathcal{Z}_{21}-(\widetilde{\mathcal{J}}_{21}-\widetilde{\mathcal{R}}_{21})\Sigma_{1}&\mathcal{Z}_{22}\end{bmatrix}\right\|_{\mathrm{F}}.$

Thus, whenever we find
$\widetilde{\mathcal{J}}_{21},\widetilde{\mathcal{R}}_{21}$ such that
$\mathcal{Z}_{21}\Sigma_{1}^{-1}=\widetilde{\mathcal{J}}_{21}-\widetilde{\mathcal{R}}_{21}$
we only have to ensure that $\widetilde{\mathcal{J}}$ remains skew-symmetric
and $\widetilde{\mathcal{R}}$ remains symmetric positive semi-definite. A
simple way to achieve this is via the choice

$\widetilde{\mathcal{J}}_{21}^{\star}\vcentcolon=\mathcal{Z}_{21}\Sigma_{1}^{-1}\qquad\text{and}\qquad\widetilde{\mathcal{R}}_{21}^{\star}=0.$
(3.14)

In view of Algorithm 1, where we only require an initialization for
$\widetilde{\mathcal{R}}$, we can directly use the result of Theorem 3.7.
Nevertheless, in our numerical experiments, we observe that the initialization
with $\widetilde{\mathcal{J}}^{\star}$ and $\widetilde{\mathcal{R}}^{\star}$
with the modification from (3.14) already yields promising results such that
they can be used even without a further application of Algorithm 1. However,
in general these are not the optimal choices as the next example illustrates.

######  Example 3.12.

Consider again Example 3.9 with the minimizers of the weighted problem
$\widetilde{\mathcal{J}}^{\star}$ and $\widetilde{\mathcal{R}}^{\star}$ as
presented in (3.9). We obtain
$\left\|\mathcal{Z}-(\widetilde{\mathcal{J}}^{\star}-\widetilde{\mathcal{R}}^{\star})\mathcal{T}\right\|_{\mathrm{F}}=\sqrt{\frac{5}{2}}$
for the original pHDMD minimization problem. Nevertheless, for

$\widetilde{\mathcal{J}}^{(1)}=\begin{bmatrix}0&-\frac{1}{5}\\\
\frac{1}{5}&\phantom{-}0\end{bmatrix}$

we obtain
$\left\|\mathcal{Z}-(\widetilde{\mathcal{J}}^{(1)}-\widetilde{\mathcal{R}}^{\star})\mathcal{T}\right\|_{\mathrm{F}}=\sqrt{\frac{200}{105}}<\sqrt{\frac{5}{2}}$,
detailing that $\widetilde{\mathcal{J}}^{\star}$ is not the optimal choice.

##  4\. Numerical experiments

In this section, we demonstrate the theoretical discussion on two exemplary
port-Hamiltonian systems. The first one is a Mass-Spring-Damper system taken
from [18] and the second one is a linear poroelastic network model; for
details see [1]. For our numerical experiments, we stop Algorithm 1 whenever

$\tfrac{\left\|\widetilde{\mathcal{J}}^{(k+1)}-\widetilde{\mathcal{J}}^{(k)}\right\|_{\mathrm{F}}}{\left\|\widetilde{\mathcal{J}}^{(k+1)}\right\|_{\mathrm{F}}}+\tfrac{\left\|\widetilde{\mathcal{R}}^{(k+1)}-\widetilde{\mathcal{R}}^{(k)}\right\|_{\mathrm{F}}}{\left\|\widetilde{\mathcal{R}}^{(k+1)}\right\|_{\mathrm{F}}}\leq\varepsilon$

with prescribed tolerance $\varepsilon$. In our experiments, we use
$\varepsilon\vcentcolon=1\mathrm{e-}10$. To report the progress during our
iterative method, we introduce the relative values of the cost functional of
3.1 and 3.6, which we denote by

$f(\widetilde{\mathcal{J}}^{(k)},\widetilde{\mathcal{R}}^{(k)})=\tfrac{\left\|\mathcal{Z}-(\widetilde{\mathcal{J}}^{(k)}-\widetilde{\mathcal{R}}^{(k)})\mathcal{T}\right\|_{\mathrm{F}}}{\left\|\mathcal{Z}\right\|_{\mathrm{F}}}\quad\text{and}\quad
f_{\mathcal{T}}(\widetilde{\mathcal{J}}^{(k)},\widetilde{\mathcal{R}}^{(k)})=\tfrac{\left\|\mathcal{T}^{\mathsf{T}}\mathcal{Z}-\mathcal{T}^{\mathsf{T}}(\widetilde{\mathcal{J}}^{(k)}-\widetilde{\mathcal{R}}^{(k)})\mathcal{T}\right\|_{\mathrm{F}}}{\left\|\mathcal{T}^{\mathsf{T}}\mathcal{Z}\right\|_{\mathrm{F}}},$

respectively. Moreover, since we are working with academic toy examples, we
report the $\mathcal{H}_{2}$ and $\mathcal{H}_{\infty}$ errors for the
identified systems, defined as

$\left\|\mathcal{G}\right\|_{\mathcal{H}_{2}}=\left(\frac{1}{2\pi}\int_{-\infty}^{\infty}\left\|\mathcal{G}(\imath\omega)\right\|_{\mathrm{F}}^{2}\mathrm{d}\omega\right)^{\frac{1}{2}}\quad\text{and}\quad\left\|\mathcal{G}\right\|_{\mathcal{H}_{\infty}}=\sup_{\omega\in\mathbb{R}}\left\|\mathcal{G}(\imath\omega)\right\|_{\mathrm{F}},$

with $\mathcal{G}$ is the transfer function of the error system. Note that
these errors require access to the original linear time-invariant system,
which of course is not available in practical applications.

######  Remark 4.1.

To have a fair comparison of our method with OI, we use also the implicit
midpoint information as defined in 3.1 for OI in all our numerical examples.

The code and data used to generate the subsequent results are accessible via
doi:10.5281/zenodo.6497497 under MIT Common License.

###  4.1. SISO Mass-Spring-Damper system

$m_{1}$$m_{2}$$m_{\frac{n}{2}}$$\scriptstyle\cdots$$k_{1}$$c_{1}$$k_{2}$$c_{2}$$k_{\frac{n}{2}}$$c_{\frac{n}{2}}$$u$
Figure 1. Illustration of the Mass-Spring-Damper system.

In our first experiment, we want to identify the Mass-Spring-Damper system,
visualized in Figure 1 with masses $m_{i}$, spring constants $k_{i}$ and
damping constants $c_{i}\geq 0$ for $i=1,\ldots,\frac{n}{2}$. A minimal
realization of the Mass-Spring-Damper system, as pH system (2.7) for the order
$n=6$, which corresponds to three masses, three springs and three dampers is
given by

$\displaystyle\widehat{H}=\left[\begin{smallmatrix}k_{1}&0&-k_{1}&0&0&0\\\
0&\frac{1}{m_{1}}&0&0&0&0\\\ -k_{1}&0&k_{1}+k_{2}&0&-k_{2}&0\\\
0&0&0&\frac{1}{m_{2}}&0&0\\\ 0&0&-k_{2}&0&k_{2}+k_{3}&0\\\
0&0&0&0&0&\frac{1}{m_{3}}\end{smallmatrix}\right],\qquad\widehat{J}=\left[\begin{smallmatrix}0&\frac{k_{1}}{m_{1}}&0&-\frac{k_{1}}{m_{2}}&0&0\\\
-\frac{k_{1}}{m_{1}}&0&\frac{k_{1}}{m_{1}}&0&0&0\\\
0&-\frac{k_{1}}{m_{1}}&0&\frac{k_{1}+k_{2}}{m_{2}}&0&-\frac{k_{2}}{m_{3}}\\\
\frac{k_{1}}{m_{2}}&0&-\frac{k_{1}+k_{2}}{m_{2}}&0&\frac{k_{2}}{m_{2}}&0\\\
0&0&0&-\frac{k_{2}}{m_{2}}&0&\frac{k_{2}+k_{3}}{m_{3}}\\\
0&0&\frac{k_{2}}{m_{3}}&0&-\frac{k_{2}+k_{3}}{m_{3}}&0\end{smallmatrix}\right],$
$\displaystyle\widehat{R}=\operatorname{diag}(0,\tfrac{c_{1}}{m_{1}^{2}},0,\tfrac{c_{2}}{m_{2}^{2}},0,\tfrac{c_{3}}{m_{3}^{2}}),\qquad\widehat{G}^{\mathsf{T}}=\begin{bmatrix}0&\frac{1}{m_{1}}&0&0&0&0\\\
\end{bmatrix},$

and $P=0$, $S=0$, $N=0$. The parameter of the masses, springs and damper are
chosen as $m_{i}=4$, $k_{i}=4$, and $c_{i}=1$ for $i=1,\ldots,3$. The
generation of the training data $(x_{i},u_{i},y_{i})$ for $i=0,1,\ldots,M$
takes place via the implicit midpoint rule. We collect $M=100$ snapshots by
running the simulation for 4 $\mathrm{s}\text{/}$, i.e., we set
$\delta_{t}=\tfrac{1}{25}$. A suitable training input is
$u(t)=\exp(-\tfrac{t}{2})\sin(t^{2})$, since it has increasing frequency to
excite the model.

$\displaystyle{0}$$\displaystyle{1}$$\displaystyle{2}$$\displaystyle{3}$$\displaystyle{4}$Time
(s)$\displaystyle{-0.25}$$\displaystyle{0.00}$$\displaystyle{0.25}$$\displaystyle{0.50}$Input$\displaystyle
u$ (a) Training input $u$.

$\displaystyle{0}$$\displaystyle{2}$$\displaystyle{4}$$\displaystyle{6}$$\displaystyle{8}$$\displaystyle{10}$Time
(s)$\displaystyle{-1}$$\displaystyle{0}$$\displaystyle{1}$Input$\displaystyle
u$ (b) Testing input $u$.

Figure 2. Training and testing input for the SISO Mass-Spring-Damper system.

With the collected data in place, Algorithm 1 with $\tilde{r}=n$ and
$\Phi=I_{n}$ is performed using the initialization from Theorem 3.7. The
algorithm terminates after the first iteration. The evolution of the
optimization is shown in Table 1, where the values for $f_{\mathcal{T}}$, $f$,
and the relative $\mathcal{H}_{2}$ and relative $\mathcal{H}_{\infty}$ errors
are displayed in each iteration.

Table 1. Evolution of Algorithm 1 applied for the SISO Mass-Spring-Damper system in terms of $f_{\mathcal{T}}$, $f$, and the relative $\mathcal{H}_{2}$ and $\mathcal{H}_{\infty}$ error over the iterations. | $f_{\mathcal{T}}$ | $f$ | $\mathcal{H}_{2}$ | $\mathcal{H}_{\infty}$
---|---|---|---|---
$\widetilde{\mathcal{J}}^{(0)},\widetilde{\mathcal{R}}^{(0)}$ | $2.61\text{\times}{10}^{-15}$ | $1.66\text{\times}{10}^{-14}$ | $2.54\text{\times}{10}^{-08}$ | $1.73\text{\times}{10}^{-11}$
$\widetilde{\mathcal{J}}^{(1)},\widetilde{\mathcal{R}}^{(1)}$ | $2.73\text{\times}{10}^{-14}$ | $1.38\text{\times}{10}^{-14}$ | $2.89\text{\times}{10}^{-08}$ | $1.74\text{\times}{10}^{-11}$

We emphasize that the initialization yields a fairly good approximation of the
minimizer, which is only slightly improved with Algorithm 1. Moreover, since
in this scenario the rank of $\mathcal{T}=\tilde{n}$, according to Lemma 3.11
the optimal value of 3.6 provides an upper bound to 3.1,

$\displaystyle\left\|\mathcal{Z}-(\widetilde{\mathcal{J}}^{(0)}-\widetilde{\mathcal{R}}^{(0)})\mathcal{T}\right\|_{\mathrm{F}}\leq
c\left\|\mathcal{T}\mathcal{Z}-\mathcal{T}(\widetilde{\mathcal{J}}^{(0)}-\widetilde{\mathcal{R}}^{(0)})\mathcal{T}\right\|_{\mathrm{F}}=$1.31\text{\times}{10}^{-11}$$
(4.1)

with $c=$2.49\text{\times}{10}^{03}$$. The original and identified models are
simulated for the testing input, shown in Figure 2(b), for 10
$\mathrm{s}\text{/}$, i.e., on a longer time horizon than during the training.
We compare the result of our method with the result of the standard DMD
approach, where DMD identifies a discrete-time LTI system. The resulting
outputs $\widetilde{y}_{\mathrm{DMD}}$ resp. $\widetilde{y}_{\mathrm{pHDMD}}$
and their absolute error to the original output trajectory $y$ are presented
in Figure 3.

$\displaystyle{0.0}$$\displaystyle{2.5}$$\displaystyle{5.0}$$\displaystyle{7.5}$$\displaystyle{10.0}$Time$\displaystyle{-1.5}$$\displaystyle{-1.0}$$\displaystyle{-0.5}$$\displaystyle{0.0}$$\displaystyle{0.5}$Testing
output$\displaystyle
y$$\displaystyle\widetilde{y}_{\mathrm{DMD}}$$\displaystyle\widetilde{y}_{\mathrm{pHDMD}}$
(a) Output trajectories.

$\displaystyle{0.0}$$\displaystyle{2.5}$$\displaystyle{5.0}$$\displaystyle{7.5}$$\displaystyle{10.0}$Time$\displaystyle{10^{-13}}$$\displaystyle{10^{-7}}$$\displaystyle{10^{-1}}$Absolute
errorDMDpHDMD (b) Absolute error of the output trajectories.

Figure 3. Output trajectories of the original and identified systems and the
corresponding absolute error for the SISO Mass-Spring-Damper system.
($\delta_{t}=$4.00\text{\times}{10}^{-02}$$)

We notice that DMD identifies an unstable system with five (out of 6) unstable
eigenvalues. This can not happen with our method, which guarantees that the
identified system is stable. However, if we decrease the stepsize, e.g., to
$\delta_{t}=$1\text{\times}{10}^{-4}$$, yielding a total of $M=40000$ data
points, DMD is also able to identify a stable system. Nevertheless, the error
of DMD is still significantly larger than pHDMD; see Figure 4(a). Let us
emphasize that the results of pHDMD are, at least in parts, attributed to the
generation of the data with the implicit midpoint rule. If we use a different
time discretization scheme instead, for instance, the Runge-Kutta 45 (RK45)
method, then the approximation quality of pHDMD reduces several orders of
magnitudes; see Figure 4(b). The reason for this behavior is that the data
generated by RK45 does not satisfy the discrete dissipation inequality (2.11).
As before, applying DMD to data generated by evaluating the RK45 solution on a
grid with the larger step size $\delta_{t}=$4.00\text{\times}{10}^{-02}$$
yields an unstable system (cf. Figure 4(b)). If we use RK45 with the small
step size, then DMD produces a slightly better approximation than pHDMD.

$\displaystyle{0.0}$$\displaystyle{2.5}$$\displaystyle{5.0}$$\displaystyle{7.5}$$\displaystyle{10.0}$Time$\displaystyle{10^{-14}}$$\displaystyle{10^{-9}}$$\displaystyle{10^{-4}}$Absolute
errorDMDpHDMD (a) Implicit midpoint with
$\delta_{t}=$1\text{\times}{10}^{-4}$$.

$\displaystyle{0.0}$$\displaystyle{2.5}$$\displaystyle{5.0}$$\displaystyle{7.5}$$\displaystyle{10.0}$Time$\displaystyle{10^{-12}}$$\displaystyle{10^{-6}}$$\displaystyle{10^{0}}$Absolute
errorDMDpHDMD (b) RK45 with $\delta_{t}=$4\text{\times}{10}^{-2}$$.

Figure 4. Absolute error for pHDMD and DMD using data generated by the
implicit midpoint rule and RK45.

###  4.2. SISO Mass-Spring-Damper system with noisy data

We repeat the experiment from Section 4.1 but add noise to the data. In
particular, we add Gaussian noise with a standard deviation
$s=$1\text{\times}{10}^{-4}$$ to the training data ($x_{i}$, $y_{i}$) and
compare the results of our method with the results of OI, which also
identifies a continuous time LTI system. The comparison is displayed in Figure
5.

$\displaystyle{0.0}$$\displaystyle{2.5}$$\displaystyle{5.0}$$\displaystyle{7.5}$$\displaystyle{10.0}$Time$\displaystyle{-0.50}$$\displaystyle{-0.25}$$\displaystyle{0.00}$$\displaystyle{0.25}$Testing
output$\displaystyle
y$$\displaystyle\widetilde{y}_{\mathrm{OI}}$$\displaystyle\widetilde{y}_{\mathrm{pHDMD}}$
(a) Output trajectories.

$\displaystyle{0}$$\displaystyle{5}$$\displaystyle{10}$$\displaystyle{15}$$\displaystyle{20}$Time$\displaystyle{10^{-4}}$$\displaystyle{10^{0}}$$\displaystyle{10^{4}}$Absolute
errorOIpHDMD (b) Absolute error of the output trajectories.

Figure 5. Output trajectories of the original and identified systems and the
corresponding absolute error for the SISO Mass-Spring-Damper system with noisy
data (standard deviation $s$ = $1\text{\times}{10}^{-4}$).

We observe that the error of OI increases over time while the error of pHDMD
remains small. This is again because OI identifies an unstable system. In
contrast, our method guarantees that the identified system is stable. We want
to emphasize that the same argument holds for passivity

###  4.3. MIMO Mass-Spring-Damper system

In this subsection, we apply the proposed dimensionality reduction to the
Mass-Spring-Damper system with increased order $n=100$. We also add a second
input $u_{2}$ (and consequently also a second output), which is applied to
$m_{2}$ the same way as $u_{1}$ to $m_{1}$ in Figure 1. We decrease
$\delta_{t}$ to generate sufficiently rich data. As training input we use
$u_{1}(t)=\exp(-\tfrac{t}{2})\sin(t^{2})$ and
$u_{2}(t)=\exp(-\tfrac{t}{2})\cos(t^{2})$. In Figure 6, we compare our method
with reduced OI and _proper orthogonal decomposition_ (POD) by plotting the
relative $\mathcal{H}_{2}$ error for different values of $\tilde{r}$.
Moreover, we perturb the data by adding Gaussian noise with standard deviation
$s$.

$\displaystyle{20}$$\displaystyle{40}$$\displaystyle{60}$$\displaystyle{80}$$\displaystyle{100}$Reduced
order$\displaystyle{10^{-7}}$$\displaystyle{10^{-6}}$$\displaystyle{10^{-5}}$$\displaystyle{10^{-4}}$$\displaystyle{10^{-3}}$$\displaystyle{10^{-2}}$$\displaystyle{10^{-1}}$$\displaystyle{10^{0}}$$\displaystyle\mathcal{H}_{2}$
errorPODOIpHDMDPOD ($\displaystyle s=10^{-4}$)OI ($\displaystyle
s=10^{-4}$)pHDMD ($\displaystyle s=10^{-4}$)POD ($\displaystyle s=10^{-6}$)OI
($\displaystyle s=10^{-6}$)pHDMD ($\displaystyle s=10^{-6}$) Figure 6.
$\mathcal{H}_{2}$-norm of the error system over reduced orders $\tilde{r}$ for
different noise levels with standard deviation $s$.

We observe that the $\mathcal{H}_{2}$ error of the non-intrusive reduced
models is fairly close to the $\mathcal{H}_{2}$ error of the reduced POD model
for both pHDMD and OI whenever the noise is sufficiently small.

###  4.4. MIMO port-Hamiltonian Poroelastic Network Model

In our final experiment, we apply our algorithm to a linear poroelasticity
problem and use the pH formulation discussed in [1], see also [12]. The model
has an order of $n=980$, with two inputs and two outputs ($m=2$). We generate
$M=10000$ snapshots with the same settings as before and obtain
$\operatorname{rank}(\mathcal{T})=100<\tilde{n}=982$. Hence, we cannot expect
to identify the complete system dynamics. Nevertheless, Algorithm 1 yields
matrices $\widetilde{\mathcal{J}}^{\star}$, $\widetilde{\mathcal{R}}^{\star}$
with
$f(\widetilde{\mathcal{J}}^{\star},\widetilde{\mathcal{R}}^{\star})=$7.11\text{\times}{10}^{-09}$$,
i.e., 3.1 is solved fairly accurately. To verify the approximation quality, we
simulate the identified system with the testing input Figure 2(b), then we
approximately recover the true output as reported in Figure 7(a) for the first
0.1 $\mathrm{s}\text{/}$ with the corresponding absolute error displayed in
Figure 7(b).

$\displaystyle{0.00}$$\displaystyle{0.02}$$\displaystyle{0.04}$$\displaystyle{0.06}$$\displaystyle{0.08}$$\displaystyle{0.10}$Time$\displaystyle{-2}$$\displaystyle{-1}$$\displaystyle{0}$$\displaystyle{1}$Testing
output$\displaystyle y$$\displaystyle\widetilde{y}_{\mathrm{pHDMD}}$ (a) Time-
discretized output $y_{1}$ for the poroelastic network model.

$\displaystyle{0.000}$$\displaystyle{0.025}$$\displaystyle{0.050}$$\displaystyle{0.075}$$\displaystyle{0.100}$Time$\displaystyle{10^{-4}}$$\displaystyle{10^{-3}}$$\displaystyle{10^{-2}}$Absolute
errorpHDMD (b) Absolute error of $y_{1}$ for the poroelastic network model.

Figure 7. Time-discretized output of the original and identified system and
the corresponding absolute error for the poroelastic network model.

For the full simulation until 10 $\mathrm{s}\text{/}$ the identified system
results in a relative $L^{2}$ error between the output of the original system
and the output of the identified system of $2.95\text{\times}{10}^{-02}$ and a
relative $L^{\infty}$ error of $2.37\text{\times}{10}^{-02}$. In comparison,
if we use the training input, then the relative $L^{2}$ and $L^{\infty}$
errors are $5.39\text{\times}{10}^{-08}$ and $3.47\text{\times}{10}^{-07}$,
respectively.

##  5\. Conclusions

We have developed a physics-informed system identification and dimensionality
reduction algorithm that uses time-domain samples of the input, state, and
output to infer a linear continuous-time dynamical system. Our algorithm is
physics-informed in the sense that for a prescribed energy functional, the
identified system is guaranteed to satisfy a dissipation inequality along any
solution of the system. Hence, the identified system is guaranteed to be
stable and passive independently of the data used for the identification. To
achieve this goal, we present a generalization of dynamic mode decomposition
and operator inference to port-Hamiltonian systems. The resulting method is an
iterative algorithm based on a fast gradient method. For the initialization,
we study a weighted problem, where we use the dominant information of the data
as a weighting factor. For the weighted problem, we derive the analytical
solution and detail in the numerical examples that the initialization is close
to the optimum.

### Acknowledgments

The work of R. Morandin is funded by the Deutsche Forschungsgemeinschaft (DFG)
within the CRC/Transregio 154 _Mathematical Modelling, Simulation and
Optimization using the Example of Gas Networks_ and the Werner-Von-Siemens
Centre for Industry and Science within the project _Maintenance, Repair &
Overhaul_. J. Nicodemus and B. Unger acknowledge funding from the DFG under
Germany’s Excellence Strategy – EXC 2075 – 390740016 and are thankful for
support by the Stuttgart Center for Simulation Science (SimTech). The authors
like to thank the anonymous referees for valuable comments that significantly
improved the manuscript.

## References

  * [1] R. Altmann, V. Mehrmann, and B. Unger.  Port-Hamiltonian formulations of poroelastic network models.  Math. Comput. Model. Dyn. Sys., 27(1):429–452, 2021.
  * [2] J. Annoni, P. Gebraad, and P. Seiler.  Wind farm flow modeling using an input-output reduced-order model.  In Am. Control Conf., Boston, MA, USA, 2016. IEEE.
  * [3] A. C. Antoulas, S. Lefteriu, and A. C. Ionita.  Chapter 8: A tutorial introduction to the Loewner framework for model reduction.  In P. Benner, A. Cohen, M. Ohlberger, and K. Willcox, editors, Model Reduction and Approximation, pages 335–376. SIAM, Philadelphia, PA, USA, 2017.
  * [4] P. J. Baddoo, B. Herrmann, B. J. McKeon, J. N. Kutz, and S. L. Brunton.  Physics-informed dynamic mode decomposition (piDMD).  ArXiv e-print 2112.04307, 2021.
  * [5] C. Beattie, V. Mehrmann, and P. Van Dooren.  Robust port-Hamiltonian representations of passive systems.  Automatica J. IFAC, 100:182–186, 2019.
  * [6] C. Beattie, V. Mehrmann, H. Xu, and H. Zwart.  Port-Hamiltonian descriptor systems.  Math. Control Signals Systems, 30(17):1–27, 2018.
  * [7] P. Benner, P. Goyal, J. Heiland, and I. Pontes Duff.  Operator inference and physics-informed learning of low-dimensional models for incompressible flows.  Electron. Trans. Numer. Anal., 56:28–51, 2022.
  * [8] P. Benner, P. Goyal, and P. Van Dooren.  Identification of port-Hamiltonian systems from frequency response data.  Systems Control Lett., 143:104741, 2020.
  * [9] P. Benner, C. Himpe, and T. Mitchell.  On reduced input-output dynamic mode decomposition.  Adv. Comput. Math., pages 1–18, 2018.
  * [10] P. Borja, J. M. A. Scherpen, and K. Fujimoto.  Extended balancing of continuous LTI systems: a structure-preserving approach.  IEEE Trans. Automat. Control, 2021.
  * [11] T. Breiten, R. Morandin, and P. Schulze.  Error bounds for port-Hamiltonian model and controller reduction based on system balancing.  Comput. Math. Appl., 116:100–115, 2022.
  * [12] T. Breiten and B. Unger.  Passivity preserving model reduction via spectral factorization.  Automatica J. IFAC, 142:110368, 2022.
  * [13] R. T. Q. Chen, Y. Rubanova, J. Bettencourt, and D. Duvenaud.  Neural ordinary differential equations.  In 32nd Conference on Neural Information Processing Systems (NeurIPS 2018), 2018.
  * [14] K. Cherifi, V. Mehrmann, and K. Hariche.  Numerical methods to compute a minimal realization of a port-Hamiltonian system.  ArXiv e-print 1903.07042, 2019.
  * [15] Y.-B. Deng, X.-Y. Hu, and L. Zhang.  Least squares solution of $BXA^{T}=T$ over symmetric, skew-symmetric, and positive semidefinite $X$.  SIAM J. Matrix Anal. Appl., 25(2):486–494, 2003.
  * [16] N. Gillis and P. Sharma.  On computing the distance to stability for matrices using linear dissipative hamiltonian systems.  Automatica J. IFAC, 85:113–121, 2017.
  * [17] N. Gillis and P. Sharma.  A semi-analytical approach for the positive semidefinite procrustes problem.  Linear Algebra Appl., 540:112–137, 2018.
  * [18] S. Gugercin, Rostyslav V. Polyuga, C. Beattie, and A. van der Schaft.  Structure-preserving tangential interpolation for model reduction of port-Hamiltonian systems.  Automatica J. IFAC, 48(9):1963–1974, 2012.
  * [19] B. Gustavsen and A. Semlyen.  Rational approximation of frequency domain responses by vector fitting.  IEEE Trans. Power Deliv., 14(3):1052–1061, 1999.
  * [20] J. Heiland and B. Unger.  Identification of linear time-invariant systems with dynamic mode decomposition.  Mathematics, 10(3):418, 2022.
  * [21] Nicholas J Higham.  Matrix nearness problems and applications.  1988\.
  * [22] B. Hillebrecht and B. Unger.  Certified machine learning: A posteriori error estimation for physics-informed neural networks.  In 2022 International Joint Conference on Neural Networks (IJCNN), pages 1–8, 2022.
  * [23] B. Hillebrecht and B. Unger.  Certified machine learning: Rigorous a posteriori error bounds for PDE defined PINNs.  ArXiv e-print 2210.03426, 2022.
  * [24] B. Jacob and H. Zwart.  Linear port-Hamiltonian systems on infinite-dimensional spaces.  Operator Theory: Advances and Applications. Birkhäuser, Basel, 2012\.
  * [25] J.-N Juang and R. S. Pappa.  An eigensystem realization algorithm for modal parameter identification and model reduction.  J. Guidance Control Dynam., 8(5):620–627, 1985.
  * [26] G. E. Karniadakis, I. G. Kevrekidis, L. Lu, P. Perdikaris, S. Wang, and L. Yang.  Physics-informed machine learning.  Nature Reviews Physics, 3:422–440, 2021.
  * [27] P. Kotyczka and L. Lefèvre.  Discrete-time port-Hamiltonian systems based on Gauss-Legendre collocation.  IFAC-PapersOnLine, 51(3):125–130, 2018.
  * [28] J. Kutz, S. Brunton, B. Brunton, and J. Proctor.  Dynamic Mode Decomposition.  SIAM, Philadelphia, PA, 2016.
  * [29] A. J. Mayo and A. C. Antoulas.  A framework for the solution of the generalized realization problem.  Linear Algebra Appl., 425(2-3):634–662, 2007.
  * [30] V. Mehrmann and R. Morandin.  Structure-preserving discretization for port-Hamiltonian descriptor systems.  In 58th IEEE Conference on Decision and Control (CDC), Nice, France, pages 6863–6868, 2019.
  * [31] V. Mehrmann and B. Unger.  Control of port-Hamiltonian differential-algebraic systems and applications.  ArXiv e-print 2201.06590, 2022.
  * [32] T. Moser and B. Lohmann.  A new Riemannian framework for efficient $\mathcal{H}_{2}$-optimal model reduction of port-Hamiltonian systems.  In 2020 59th IEEE Conference on Decision and Control (CDC), pages 5043–5049, 2020.
  * [33] Y. Nesterov.  Introductory lectures on convex optimization: A basic course, volume 87.  Springer Science & Business Media, 2003.
  * [34] B. Peherstorfer, S. Gugercin, and K. E. Willcox.  Data-driven reduced model construction with time-domain Loewner models.  SIAM J. Sci. Comput., 39(5):2152–2178, 2017.
  * [35] 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.
  * [36] R. V. Polyuga and A. van der Schaft.  Structure-preserving moment matching for port-Hamiltonian systems: Arnoldi and Lanczos.  IEEE Trans. Automat. Control, 56(6):1458–1462, 2011.
  * [37] R. V. Polyuga and A. van der Schaft.  Effort- and flow-constraint reduction methods for structure preserving model reduction of port-Hamiltonian systems.  Systems Control Lett., 61(3):412–421, 2012.
  * [38] J. L. Proctor, S. L. Brunton, and J. N. Kutz.  Dynamic mode decomposition with control.  SIAM J. Appl. Dyn. Syst., 15(1):142–161, 2014.
  * [39] M. Raissi, P. Perdikaris, and G.E. Karniadakis.  Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations.  J. Comput. Phys., 378:686–707, 2019.
  * [40] K. Sato and H. Sato.  Structure-Preserving $H^{2}$ Optimal Model Reduction Based on the Riemannian Trust-Region Method.  IEEE Trans. Automat. Control, 63(2):505–512, 2018.
  * [41] P. J. Schmid.  Dynamic mode decomposition of numerical and experimental data.  J. Fluid Mech., 656:5–28, 2010.
  * [42] P. Schulze and B. Unger.  Data-driven interpolation of dynamical systems with delay.  Systems Control Lett., 97:125–131, 2016.
  * [43] P. Schulze, B. Unger, C. Beattie, and S. Gugercin.  Data-driven structured realization.  Linear Algebra Appl., 537:250–286, 2018.
  * [44] P. Schwerdtner.  Port-Hamiltonian system identification from noisy frequency response data.  ArXiv e-print 2106.11355, 2021.
  * [45] P. Schwerdtner and M. Voigt.  SOBMOR: Structured optimization-based model order reduction.  ArXiv e-print 2011.07567, 2020.
  * [46] P. Schwerdtner and M. Voigt.  Adaptive sampling for structure-preserving model order reduction of port-Hamiltonian systems.  IFAC-PapersOnLine, 54(19):143–148, 2021.
  * [47] H. Sharma and B. Kramer.  Preserving Lagrangian structure in data-driven reduced-order modeling of large-scale mechanical systems.  ArXiv e-print 2203.06361, 2022.
  * [48] H. Sharma, Z. Wang, and B. Kramer.  Hamiltonian operator inference: Physics-preserving learning of reduced-order models for canonical Hamiltonian systems.  Phys. D, 432:133122, 2022.
  * [49] Jonathan H. Tu, Clarence W. Rowley, Dirk M. Luchtenburg, Steven L. Brunton, and J. Nathan Kutz.  On dynamic mode decomposition: Theory and applications.  J. Comput. Dyn., 1(2):391–421, 2014.
  * [50] A. van der Schaft and D. Jeltsema.  Port-Hamiltonian systems theory: An introductory overview.  Foundations and Trends in Systems and Control, 1(2-3):173–378, 2014\.
  * [51] S. W. R. Werner, I. V. Gosea, and S. Gugercin.  Structured vector fitting framework for mechanical systems.  IFAC-PapersOnLine, 55(20):163–168, 2022.
  * [52] T. Wolf, B. Lohmann, R. Eid, and P. Kotyczka.  Passivity and structure preserving order reduction of linear port-Hamiltonian systems using Krylov subspaces.  Eur. J. Control, 16(4):401–406, 2010.