Run this tutorial

Click here to run this tutorial on binder:# Tutorial: Model order reduction for port-Hamiltonian systems¶

In the first section, we introduce the class of
port-Hamiltonian systems and their relationship to two other system-theoretic properties called
*passivity* and *positive realness*. After introducing a toy example in
the second section, we
look into structure-preserving model reduction schemes for port-Hamiltonian systems
in the third section.

## Port-Hamiltonian LTI systems¶

Port-Hamiltonian systems have several favorable properties for modeling, control and simulation, for example, composability and stability. Furthermore, they adhere to a power balance equation. Port-Hamiltonian systems are especially suited for network-based modeling and problems involving multi-physics phenomena. We refer to [MU23] for a general introduction to port-Hamiltonian descriptor systems and their applications.

We say a LTI system is *port-Hamiltonian* if it can be expressed as

with \(H := Q^T E\), and if the structure matrix

and the dissipation matrix

satisfy \(H = H^T \succ 0\), \(\Gamma^T = -\Gamma\), and \(\mathcal{W} = \mathcal{W}^T \succcurlyeq 0\).

The quadratic (energy) function \(\mathcal{H}(x) := \tfrac{1}{2} x^T H x\), typically called Hamiltonian, corresponds to the energy stored in the system. In applications, \(E\) and/or \(Q\) often are identity matrices.

It is known that if the LTI system is minimal and stable, the following are equivalent:

The system is passive.

The system is port-Hamiltonian.

The system is positive real.

See for example [BU22] for more details.

In pyMOR, there exists a `PHLTIModel`

class. Currently, pyMOR only supports
port-Hamiltonian systems with nonsingular E. `PHLTIModel`

inherits from
`LTIModel`

, so `PHLTIModel`

can be used with all reductors that expect
an `LTIModel`

. For model reduction, it is often desirable to preserve the
port-Hamiltonian structure, i.e., to compute a ROM that is also port-Hamiltonian.

If desired, a passive `LTIModel`

can be converted into a `PHLTIModel`

using
the `from_passive_LTIModel`

method.
Consequentely, one option to preserve port-Hamiltonian structure is to use a reductor
that preserves passivity (but returns a ROM of type `LTIModel`

) and convert the
ROM into a `PHLTIModel`

in a post-processing step.

## A toy problem: Mass-spring-damper chain¶

As a toy problem, we use a mass-spring-damper chain, which can be formulated as a port-Hamiltonian system (see [GPBvdS12]):

Here, the spring constants are denoted by \(k_i\) and the damping constants by
\(c_i\), \(i=1,\dots,n/2\).
The inputs \(u_1\) and \(u_2\) are the external forces on the first two
masses \(m_1\) and \(m_2\). The system outputs \(y_1\) and \(y_2\)
correspond to the velocities of the first two masses \(m_1\) and \(m_2\).
The toy problem is included in pyMOR in the `pymor.models.examples`

module as
`msd_example`

.

## Structure-preserving model order reduction¶

pyMOR provides three reductors which can be used for model order reduction while preserving the port-Hamiltonian structure:

pH-IRKA (

`PHIRKAReductor`

) [GPBvdS12],PRBT (

`PRBTReductor`

) [DP84, GA04, HJS84],passivity preserving model reduction via spectral factorization (

`SpectralFactorReductor`

) [BU22].

In this section, we apply all three reductors on our toy example and compare their performance. All three reductors are described in [BU22] in more detail.

Note: Currently, the `PRBTReductor`

and
`SpectralFactorReductor`

reductors require
the symmetric part of \(D\) (i.e., the \(S\) matrix in the port-Hamiltonian system)
to be nonsingular. The MSD example has a zero \(D\) matrix. Therefore,
we have to add a small regularization feedthrough term, i.e., we replace \(D\) with
\(D+\varepsilon I_m\). This is a limitation of the current implementation since the
numerical solution of the KYP-LMI is obtained by solving a related Riccati
equation, for instance

which is only possible if \(D + D^\top\) is nonsingular.
For `from_passive_LTIModel`

, \(D + D^\top\) must be
nonsingular for the same reasons.

```
import numpy as np
from pymor.models.examples import msd_example
fom = msd_example(50, 2)
# tolerance for solving the Riccati equation instead of KYP-LMI
# by introducing a regularization feedthrough term D
# (required for PRBTReductor and SpectralFactorReductor reductors)
S = fom.S.matrix.copy()
S += np.eye(S.shape[0]) * 1e-12
fom = fom.with_(S=fom.S.with_(matrix=S),
solver_options={'ricc_pos_lrcf': 'slycot'})
```

The `ricc_pos_lrcf`

solver option refers to the solver used for the underlying
Riccati equation relevant for `PRBTReductor`

and
`SpectralFactorReductor`

. Possible choices are
`scipy`

or `slycot`

(if installed). Currently, we recommend `slycot`

, since `scipy`

gets into trouble if the associated Hamiltonian pencil has eigenvalues close to the
imaginary axis.

### pH-IRKA¶

The pH-IRKA reductor `PHIRKAReductor`

directly returns
a ROM of type `PHLTIModel`

. pH-IRKA works similar to the standard IRKA reductor
`IRKAReductor`

but with fewer degrees of freedom to preserve
the port-Hamiltonian structure. In more detail, the IRKA fixed-point iteration is performed,
but the left projection matrix is chosen as \(W = QV\), which then automatically yields a
reduced pH system with \(\hat{Q} = I_r\).

```
from pymor.reductors.ph.ph_irka import PHIRKAReductor
reductor = PHIRKAReductor(fom)
rom1 = reductor.reduce(10)
print(f'rom1 is of type {type(rom1).__qualname__}.')
```

```
rom1 is of type PHLTIModel.
```

### Positive-real balanced truncation (PRBT)¶

Positive-real balanced truncation (PRBT) works analogously to the standard balanced truncation
method described in Tutorial: Reducing an LTI system using balanced truncation, but uses positive real controllability
and observability Gramians instead. PRBT preserves passivity but returns a ROM
of type `LTIModel`

. Thus, we convert the ROM into a `PHLTIModel`

in a
post-processing step. Note that PRBT can be used with any passive `LTIModel`

FOM.

```
from pymor.models.iosys import PHLTIModel
from pymor.reductors.bt import PRBTReductor
reductor = PRBTReductor(fom)
rom2 = reductor.reduce(10)
rom2 = rom2.with_(solver_options={'ricc_pos_lrcf': 'slycot'})
rom2 = PHLTIModel.from_passive_LTIModel(rom2)
print(f'rom2 is of type {type(rom2).__qualname__}.')
```

```
rom2 is of type PHLTIModel.
```

### Passivity preserving model reduction via spectral factorization¶

The `SpectralFactorReductor`

method
is a wrapper reductor for another generic reductor. The method extracts a
spectral factor from the FOM (this is only possible if the system is passive),
which subsequentely is reduced by a reductor specified by the user.
A spectral factor is a standard `LTIModel`

, and hence any LTI reduction can be used.
For our example, we use the `IRKAReductor`

as the inner reductor.
If the inner reductor returns a stable ROM, passivity is preserved.
The spectral factor method and PRBT are related since the computation of
the optimal spectral factor for model reduction depends on the computation of
the positive-real observability Gramian. The spectral factor method can be used with
any passive `LTIModel`

FOM. Again, we convert the ROM of type `LTIModel`

into a
`PHLTIModel`

in a post-processing step.

```
from pymor.reductors.spectral_factor import SpectralFactorReductor
from pymor.reductors.h2 import IRKAReductor
reductor = SpectralFactorReductor(fom)
rom3 = reductor.reduce(
lambda spectral_factor, mu : IRKAReductor(spectral_factor, mu).reduce(10)
)
rom3 = rom3.with_(solver_options={'ricc_pos_lrcf': 'slycot'})
rom3 = PHLTIModel.from_passive_LTIModel(rom3)
print(f'rom3 is of type {type(rom3).__qualname__}.')
```

```
rom3 is of type PHLTIModel.
```

### Comparison¶

Let us compare the \(\mathcal{H}_2\) errors of the three methods:

```
err1 = fom - rom1
err2 = fom - rom2
err3 = fom - rom3
print(f'pHIRKA - Relative H2 error: {err1.h2_norm() / fom.h2_norm():.3e}')
print(f'PRBT - Relative H2 error: {err2.h2_norm() / fom.h2_norm():.3e}')
print(f'spectral_factor - Relative H2 error: {err3.h2_norm() / fom.h2_norm():.3e}')
```

```
pHIRKA - Relative H2 error: 2.430e-01
PRBT - Relative H2 error: 1.031e-02
spectral_factor - Relative H2 error: 3.239e-03
```

We can plot a magnitude plot of the three error systems:

```
import matplotlib.pyplot as plt
w = (1e-4, 1e3)
fig, ax = plt.subplots()
err1.transfer_function.mag_plot(w, ax=ax, label='pHIRKA')
err2.transfer_function.mag_plot(w, ax=ax, linestyle='--', label='PRBT')
err3.transfer_function.mag_plot(w, ax=ax, label='spectral_factor')
_ = ax.legend()
```

Download the code:
`tutorial_ph.md`

,
`tutorial_ph.ipynb`

.