Tutorial: Reducing a heat equation using balanced truncation
============================================================

.. include:: jupyter_init.txt

Heat equation
-------------

We consider the following one-dimensional heat equation over :math:`(0, 1)` with
two inputs :math:`u_1, u_2` and three outputs :math:`y_1, y_2, y_2`:

.. math::

    \begin{align}
        \partial_t T(\xi, t) & = \partial_{\xi \xi} T(\xi, t) + u_1(t),
        & 0 < \xi < 1,\ t > 0, \\
        -\partial_\xi T(0, t) & = -T(0, t) + u_2(t),
        & t > 0, \\
        \partial_\xi T(1, t) & = -T(1, t),
        & t > 0, \\
        y_1(t) & = T(0, t),
        & t > 0, \\
        y_2(t) & = T(0.5, t),
        & t > 0, \\
        y_3(t) & = T(1, t),
        & t > 0.
    \end{align}

In the following, we will create a discretized |Model| and reduce it using the
balanced truncation method to approximate the mapping from inputs
:math:`u = (u_1, u_2)` to outputs :math:`y = (y_1, y_2, y_3)`.


Discretized model
-----------------

We need to construct a linear time-invariant (LTI) system

.. math::

    \begin{align}
        E \dot{x}(t) & = A x(t) + B u(t), \\
        y(t) & = C x(t) + D u(t).
    \end{align}

In pyMOR, these models are captured by |LTIModels| from the
:mod:`pymor.models.iosys` module.

There are many ways of building an |LTIModel|.
Here, we show how to build one from custom matrices instead of using a
discretizer as in :doc:`tutorial_builtin_discretizer` and the
:meth:`~pymor.models.basic.InstationaryModel.to_lti` of |InstationaryModel|.
In particular, we will use the
:meth:`~pymor.models.iosys.LTIModel.from_matrices` method of |LTIModel|, which
instantiates an |LTIModel| from NumPy or SciPy matrices.

First, we do the necessary imports.

.. jupyter-execute::

    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

Next, we can assemble the matrices based on a centered finite difference
approximation:

.. jupyter-execute::

    k = 50
    n = 2 * k + 1

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

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

    C = np.zeros((3, n))
    C[0, 0] = C[1, k] = C[2, -1] = 1

Then, we can create an |LTIModel|:

.. jupyter-execute::

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

We can get the internal representation of the |LTIModel| `fom`

.. jupyter-execute::

    fom

From this, we see that the matrices were wrapped in |NumpyMatrixOperators|,
while default values were chosen for :math:`D` and :math:`E` matrices
(respectively, zero and identity). The operators in an |LTIModel| can be
accessed directly, e.g., `fom.A`.

We can also see some basic information from `fom`'s string representation

.. jupyter-execute::

    print(fom)

To visualize the behavior of the `fom`, we can draw its magnitude plot, i.e.,
a visualization of the mapping :math:`\omega \mapsto H(i \omega)`, where
:math:`H(s) = C (s E - A)^{-1} B + D` is the transfer function of the system.

.. jupyter-execute::

    w = np.logspace(-2, 8, 50)
    fom.mag_plot(w)
    plt.grid()

Plotting the Hankel singular values shows us how well an LTI system can be
approximated by a reduced-order model

.. jupyter-execute::

    hsv = fom.hsv()
    fig, ax = plt.subplots()
    ax.semilogy(range(1, len(hsv) + 1), hsv, '.-')
    ax.set_title('Hankel singular values')
    ax.grid()

As expected for a heat equation, the Hankel singular values decay rapidly.


Running balanced truncation
---------------------------

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. In
particular, there exist invertible transformation matrices
:math:`T, S \in \mathbb{R}^{n \times n}` such that the equivalent full-order
model with
:math:`\widetilde{E} = S^T E T = I`,
:math:`\widetilde{A} = S^T A T`,
:math:`\widetilde{B} = S^T B`,
:math:`\widetilde{C} = C T`
has Gramians :math:`\widetilde{P}` and :math:`\widetilde{Q}`, i.e., solutions to
Lyapunov equations

.. math::

    \begin{align}
        \widetilde{A} \widetilde{P}
        + \widetilde{P} \widetilde{A}^T
        + \widetilde{B} \widetilde{B}^T
        & = 0, \\
        \widetilde{A}^T \widetilde{Q}
        + \widetilde{Q} \widetilde{A}
        + \widetilde{C}^T \widetilde{C}
        & = 0, \\
    \end{align}

such that
:math:`\widetilde{P} = \widetilde{Q} = \Sigma = \operatorname{diag}(\sigma_i)`,
where :math:`\sigma_i` are the Hankel singular values.
Then, taking as basis matrices :math:`V, W \in \mathbb{R}^{n \times r}` the
first :math:`r` columns of :math:`T` and :math:`S` (possibly after
orthonormalization), gives a reduced-order model

.. math::

    \begin{align}
        \widehat{E} \dot{\widehat{x}}(t)
        & = \widehat{A} \widehat{x}(t) + \widehat{B} u(t), \\
        \widehat{y}(t)
        & = \widehat{C} \widehat{x}(t) + D u(t),
    \end{align}

with
:math:`\widehat{E} = W^T E V`,
:math:`\widehat{A} = W^T A V`,
:math:`\widehat{B} = W^T B`,
:math:`\widehat{C} = C V`,
which satisfies the :math:`\mathcal{H}_\infty` (i.e., induced
:math:`\mathcal{L}_2`) error bound

.. math::

    \sup_{u \neq 0} \frac{\lVert y - \widehat{y} \rVert_{\mathcal{L}_2}}{\lVert u \rVert_{\mathcal{L}_2}}
    \leqslant 2 \sum_{i = r + 1}^n \sigma_i.

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

.. math::

    \sup_{u \neq 0} \frac{\lVert y - \widehat{y} \rVert_{\mathcal{L}_2}}{\lVert u \rVert_{\mathcal{L}_2}}
    \geqslant \sigma_{r + 1}.

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

.. jupyter-execute::

    bt = BTReductor(fom)

Calling its :meth:`~pymor.reductors.bt.GenericBTReductor.reduce` method runs the
balanced truncation algorithm. This reductor additionally has an `error_bounds`
method which can compute the a priori :math:`\mathcal{H}_\infty` error bounds
based on the Hankel singular values:

.. jupyter-execute::

    error_bounds = bt.error_bounds()
    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')
    ax.grid()

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

.. jupyter-execute::

    rom = bt.reduce(10)

Instead, or in addition, a tolerance for the :math:`\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

.. jupyter-execute::

    fig, ax = plt.subplots()
    fom.mag_plot(w, ax=ax, label='FOM')
    rom.mag_plot(w, ax=ax, linestyle='--', label='ROM')
    ax.legend()
    ax.grid()

and plot the magnitude plot of the error system

.. jupyter-execute::

    (fom - rom).mag_plot(w)
    plt.grid()

We can compute the relative errors in :math:`\mathcal{H}_\infty` or
:math:`\mathcal{H}_2` (or Hankel) norm

.. jupyter-execute::

    print(f'Relative Hinf error: {(fom - rom).hinf_norm() / fom.hinf_norm():.3e}')
    print(f'Relative H2 error:   {(fom - rom).h2_norm() / fom.h2_norm():.3e}')

.. note::

    To compute the :math:`\mathcal{H}_\infty` norms, pyMOR uses the dense solver
    from Slycot, and therefore all of the operators have to be converted to
    dense matrices. For large systems, this may be very expensive.

Download the code:
:jupyter-download:script:`tutorial_bt`
:jupyter-download:notebook:`tutorial_bt`