Run this tutorial

Click here to run this tutorial on binder:
Please note that starting the notebook server may take a couple of minutes.

Tutorial: Reducing an LTI system using balanced truncation

Here we briefly describe the balanced truncation method, for asymptotically stable LTI systems with an invertible $E$ matrix, and demonstrate it on the heat equation example from Tutorial: Linear time-invariant systems. First, we import necessary packages, including BTReductor.

import matplotlib.pyplot as plt
import numpy as np
import scipy.sparse as sps
from pymor.models.iosys import LTIModel
from pymor.reductors.bt import BTReductor

plt.rcParams['axes.grid'] = True

Then we build the matrices

import numpy as np
import scipy.sparse as sps

k = 50
n = 2 * k + 1

E = sps.eye(n, format='lil')
E[0, 0] = E[-1, -1] = 0.5
E = E.tocsc()

d0 = n * [-2 * (n - 1)**2]
d1 = (n - 1) * [(n - 1)**2]
A = sps.diags([d1, d0, d1], [-1, 0, 1], format='lil')
A[0, 0] = A[-1, -1] = -n * (n - 1)
A = A.tocsc()

B = np.zeros((n, 2))
B[:, 0] = 1
B[0, 0] = B[-1, 0] = 0.5
B[0, 1] = n - 1

C = np.zeros((3, n))
C[0, :n//3] = C[1, n//3:2*n//3] = C[2, 2*n//3:] = 1
C /= C.sum(axis=1)[:, np.newaxis]

and form the full-order model.

fom = LTIModel.from_matrices(A, B, C, E=E)

Balanced truncation

As the name suggests, the balanced truncation method consists of finding a balanced realization of the full-order LTI system and truncating it to obtain a reduced-order model.

The balancing part is based on the fact that a single LTI system has many realizations. For example, starting from a realization

$$\begin{array}{r}\begin{array}{rl}E\dot{x}(t)& =Ax(t)+Bu(t),\\ y(t)& =Cx(t)+Du(t),\end{array}\end{array}$$

another realization can be obtained by replacing $x(t)$ with $T\tilde{x}(t)$ or by pre-multiplying the differential equation with an invertible matrix. In particular, there exist invertible transformation matrices $T,S\in {\mathbb{R}}^{n\times n}$ such that the realization with $\tilde{E}=S^{\mathrm{T}}ET=I$, $\tilde{A}=S^{\mathrm{T}}AT$, $\tilde{B}=S^{\mathrm{T}}B$, $\tilde{C}=CT$ has Gramians $\tilde{P}$ and $\tilde{Q}$ satisfying $\tilde{P}=\tilde{Q}=\mathrm{\Sigma }=\mathrm{diag}({\sigma }_i)$, where ${\sigma }_i$ are the Hankel singular values (see Tutorial: Linear time-invariant systems for more details). Such a realization is called balanced.

The truncation part is based on the controllability and observability energies. The controllability energy $E_c(x_0)$ is the minimum energy (squared ${\mathcal{L}}_2$ norm of the input) necessary to steer the system from the zero state to $x_0$. The observability energy $E_o(x_0)$ is the energy of the output (squared ${\mathcal{L}}_2$ norm of the output) for a system starting at the state $x_0$ and with zero input. It can be shown for the balanced realization (and same for any other realization) that, if $\tilde{P}$ is invertible, then

$$E_c(x_0)=x_0{\tilde{P}}^{-1}{x}_0,E_o(x_0)=x_0\tilde{Q}{x}_0.$$

Therefore, states corresponding to small Hankel singular values are more difficult to reach (they have a large controllability energy) and are difficult to observe (they produce a small observability energy). In this sense, it is then reasonable to truncate these states. This can be achieved by taking as basis matrices $V,W\in {\mathbb{R}}^{n\times r}$ the first $r$ columns of $T$ and $S$, possibly after orthonormalization, giving a reduced-order model

$$\begin{array}{r}\begin{array}{rl}\hat{E}\dot{\hat{x}}(t)& =\hat{A}\hat{x}(t)+\hat{B}u(t),\\ \hat{y}(t)& =\hat{C}\hat{x}(t)+Du(t),\end{array}\end{array}$$

with $\hat{E}=W^{\mathrm{T}}EV$, $\hat{A}=W^{\mathrm{T}}AV$, $\hat{B}=W^{\mathrm{T}}B$, $\hat{C}=CV$.

It is known that the reduced-order model is asymptotically stable if ${\sigma }_r>{\sigma }_{r+1}$. Furthermore, it satisfies the ${\mathcal{H}}_{\mathrm{\infty }}$ error bound

$$\Vert H-\hat{H}{\Vert }_{{\mathcal{H}}_{\mathrm{\infty }}}⩽2\sum _{i=r+1}^n{\sigma }_i.$$

Note that any reduced-order model (not only from balanced truncation) satisfies the lower bound

$$\Vert H-\hat{H}{\Vert }_{{\mathcal{H}}_{\mathrm{\infty }}}⩾{\sigma }_{r+1}.$$

Balanced truncation in pyMOR

To run balanced truncation in pyMOR, we first need the reductor object

bt = BTReductor(fom)

Calling its reduce method runs the balanced truncation algorithm. This reductor additionally has an error_bounds method which can compute the a priori ${\mathcal{H}}_{\mathrm{\infty }}$ error bounds based on the Hankel singular values:

error_bounds = bt.error_bounds()
hsv = fom.hsv()
fig, ax = plt.subplots()
ax.semilogy(range(1, len(error_bounds) + 1), error_bounds, '.-')
ax.semilogy(range(1, len(hsv)), hsv[1:], '.-')
ax.set_xlabel('Reduced order')
_ = ax.set_title(r'Upper and lower $\mathcal{H}_\infty$ error bounds')

To get a reduced-order model of order 10, we call the reduce method with the appropriate argument:

rom = bt.reduce(10)

Log Output

00:01 gram_schmidt: Orthonormalizing vector 2 again

00:01 gram_schmidt: Orthonormalizing vector 3 again

00:01 gram_schmidt: Orthonormalizing vector 4 again

00:01 gram_schmidt: Orthonormalizing vector 5 again

00:01 gram_schmidt: Orthonormalizing vector 6 again

00:01 gram_schmidt: Orthonormalizing vector 7 again

00:01 gram_schmidt: Orthonormalizing vector 8 again

00:01 gram_schmidt: Orthonormalizing vector 9 again

00:01 gram_schmidt: Orthonormalizing vector 3 again

00:01 gram_schmidt: Orthonormalizing vector 4 again

00:01 gram_schmidt: Orthonormalizing vector 5 again

00:01 gram_schmidt: Orthonormalizing vector 6 again

00:01 gram_schmidt: Orthonormalizing vector 7 again

00:01 gram_schmidt: Orthonormalizing vector 8 again

00:01 gram_schmidt: Orthonormalizing vector 9 again

00:01 LTIPGReductor: Operator projection ...

00:01 LTIPGReductor: Building ROM ...

Instead, or in addition, a tolerance for the ${\mathcal{H}}_{\mathrm{\infty }}$ error can be specified, as well as the projection algorithm (by default, the balancing-free square root method is used). The used Petrov-Galerkin bases are stored in bt.V and bt.W.

We can compare the magnitude plots between the full-order and reduced-order models

w = (1e-2, 1e3)
fig, ax = plt.subplots()
fom.transfer_function.mag_plot(w, ax=ax, label='FOM')
rom.transfer_function.mag_plot(w, ax=ax, linestyle='--', label='ROM')
_ = ax.legend()

as well as Bode plots

fig, axs = plt.subplots(6, 2, figsize=(8, 10), sharex=True, constrained_layout=True)
fom.transfer_function.bode_plot(w, ax=axs)
_ = rom.transfer_function.bode_plot(w, ax=axs, linestyle='--')

Also, we can plot the magnitude plot of the error system, which is again an LTI system.

$$\begin{array}{r}\begin{array}{rl}\left[\begin{array}{cc}E& 0\\ 0& \hat{E}\end{array}\right]\left[\begin{array}{c}\dot{x}(t)\\ \dot{\hat{x}}(t)\end{array}\right]& =\left[\begin{array}{cc}A& 0\\ 0& \hat{A}\end{array}\right]\left[\begin{array}{c}x(t)\\ \hat{x}(t)\end{array}\right]+\left[\begin{array}{c}B\\ \hat{B}\end{array}\right]u(t),\\ y(t)-\hat{y}(t)& =\left[\begin{array}{cc}C& -\hat{C}\end{array}\right]\left[\begin{array}{c}x(t)\\ \hat{x}(t)\end{array}\right].\end{array}\end{array}$$

err = fom - rom
_ = err.transfer_function.mag_plot(w)

and its Bode plot

fig, axs = plt.subplots(6, 2, figsize=(8, 10), sharex=True, constrained_layout=True)
_ = err.transfer_function.bode_plot(w, ax=axs)

Finally, we can compute the relative errors in different system norms.

print(f'Relative Hinf error:   {err.hinf_norm() / fom.hinf_norm():.3e}')
print(f'Relative H2 error:     {err.h2_norm() / fom.h2_norm():.3e}')
print(f'Relative Hankel error: {err.hankel_norm() / fom.hankel_norm():.3e}')

Relative Hinf error:   1.236e-06
Relative H2 error:     6.797e-05

Relative Hankel error: 2.377e-06

Download the code: tutorial_bt.md, tutorial_bt.ipynb.