Tutorial: Reducing an LTI system using balanced truncation

Run this tutorial

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

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

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, 0] = C[1, k] = C[2, -1] = 1

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{split}\begin{align} E \dot{x}(t) & = A x(t) + B u(t), \\ y(t) & = C x(t) + D u(t), \end{align}\end{split}\]

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^{\operatorname{T}} E T = I\), \(\tilde{A} = S^{\operatorname{T}} A T\), \(\tilde{B} = S^{\operatorname{T}} B\), \(\tilde{C} = C T\) has Gramians \(\tilde{P}\) and \(\tilde{Q}\) satisfying \(\tilde{P} = \tilde{Q} = \Sigma = \operatorname{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, \quad 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{split}\begin{align} \hat{E} \dot{\hat{x}}(t) & = \hat{A} \hat{x}(t) + \hat{B} u(t), \\ \hat{y}(t) & = \hat{C} \hat{x}(t) + D u(t), \end{align}\end{split}\]

with \(\hat{E} = W^{\operatorname{T}} E V\), \(\hat{A} = W^{\operatorname{T}} A V\), \(\hat{B} = W^{\operatorname{T}} B\), \(\hat{C} = C V\).

It is known that the reduced-order model is asymptotically stable if \(\sigma_r > \sigma_{r + 1}\). Furthermore, it satisfies the \(\mathcal{H}_\infty\) error bound

\[\lVert H - \hat{H} \rVert_{\mathcal{H}_\infty} \leqslant 2 \sum_{i = r + 1}^n \sigma_i.\]

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

\[\lVert H - \hat{H} \rVert_{\mathcal{H}_\infty} \geqslant \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}_\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)

Instead, or in addition, a tolerance for the \(\mathcal{H}_\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 = np.logspace(-2, 8, 300)
fig, ax = plt.subplots()
fom.mag_plot(w, ax=ax, label='FOM')
rom.mag_plot(w, ax=ax, linestyle='--', label='ROM')
_ = ax.legend()

as well as Bode plots

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

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

\[\begin{split}\begin{align} \begin{bmatrix} E & 0 \\ 0 & \hat{E} \end{bmatrix} \begin{bmatrix} \dot{x}(t) \\ \dot{\hat{x}}(t) \end{bmatrix} & = \begin{bmatrix} A & 0 \\ 0 & \hat{A} \end{bmatrix} \begin{bmatrix} x(t) \\ \hat{x}(t) \end{bmatrix} + \begin{bmatrix} B \\ \hat{B} \end{bmatrix} u(t), \\ y(t) - \hat{y}(t) & = \begin{bmatrix} C & -\hat{C} \end{bmatrix} \begin{bmatrix} x(t) \\ \hat{x}(t) \end{bmatrix}. \end{align}\end{split}\]

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

and its Bode plot

_ = err.bode_plot(w)

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:   3.702e-05
Relative H2 error:     6.399e-04
Relative Hankel error: 7.629e-05

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