System‐theoretic model order reduction with pyMOR

This paper shows recent developments in pyMOR, in particular the addition of system‐theoretic methods. All methods are implemented using pyMOR's abstract interfaces, which allows the application to partial differential equation (PDE) models implemented with third‐party libraries. We demonstrate this by applying balanced truncation to a PDE model discretized in FEniCS.


Introduction
pyMOR is a free software library for building model order reduction (MOR) applications in the Python programming language. All of its algorithms are based on abstract interfaces which allow seamless integration with external partial differential equation (PDE) solver packages (see [1]). Its initial focus was on reduced basis methods for linear and nonlinear parameterized PDEs. We present here recent developments in system-theoretic MOR methods and related algorithms.
In Section 2, we give a brief overview of linear systems and balanced truncation (BT). Next, we describe pyMOR's abstract interfaces and the implementation of the methods in Section 3. We demonstrate BT applied to a PDE model discretized using FEniCS in Section 4.
2 System-theoretic model order reduction We consider continuous-time, linear time-invariant systems with E, A ∈ R n×n (E invertible), B ∈ R n×m , C ∈ R p×n , input u(t) ∈ R m , state x(t) ∈ R n , and output y(t) ∈ R p . The goal of MOR is to find a reduced-order model (ROM) withÊ,Â ∈ R r×r ,B ∈ R r×m ,Ĉ ∈ R p×r , for some r n such that y −ŷ is small for every input u in some function space. BT is a projection-based method whereÊ,Â,B,Ĉ are obtained by Petrov-Galerkin projection of the original system matrices on appropriate ansatz and test spaces. The main computational cost in finding these spaces lies in solving two largescale Lyapunov equations: AP E T + EP A T + BB T = 0 and A T QE + E T QA + C T C = 0. Low-rank approximations for P and Q are computed using the low-rank alternating directions implicit (LR-ADI) method. The method consists of solving a sequence of linear systems 3 Generic interfaces for model order reduction pyMOR's PDE solver integration is based on three abstract interfaces: VectorArrays, Operators, and Models. Each VectorArray is an element of a VectorSpace, Operator is a mapping (possibly nonlinear or parametric) between two VectorSpaces, and a Model is a structured collection of Operators. The interfaces are completely agnostic of the solver's internal data representation. All reduction algorithms in pyMOR are implemented in terms of these interfaces.
In particular, in pyMOR 0.5, LTISystem is a model representing (1), where E, A, B, C are Operators. In case of the example in Section 4, these operators come from a discretization in FEniCS. Via linear combination and the Operator's apply_inverse method, the LR-ADI algorithm instructs FEniCS to solve (p i E + A)x i = b i .

Numerical example
As an example, we consider the heat equation over a cross section of a heat sink (see [2] for the problem description). Finite element discretization, using FEniCS, leads to a system of the form (1) with n = 12 296, m = 1, and p = 1. The top-middle plot in Figure 1 shows the state of the full-order model at t = 10 when the input is u(t) = sin( π 3 t) 2 . The average run time is about 8 seconds on a laptop (model: Vaio VJS132C11L; processor: Intel c Core TM i7-8550U CPU @ 1.80 GHz, 4 cores, 8 logical processors, hyper-threading activated; RAM: 8 GB; cache: 8 MB).
We applied BT (using pyMOR 0.5) to get a ROM of order 10, with a run time of about 7 seconds. Based on Hankel singular values, we find the relative H ∞ -error lies in the interval [7.27 × 10 −5 , 1.23 × 10 −4 ], shown also in the left plot of Figure 1. The relative H 2 -error can be computed and its value is 7.37 × 10 −3 . In the bottom row of Figure 1, we see the state and output errors of the ROM. As expected, we see that the output is approximated better than the state. The average run time of the simulation is about 15 milliseconds, which is roughly 500 times faster then for the full-order model.
The entire script to generate these results can be found in [2].

Conclusion and outlook
We demonstrated BT applied to PDE models discretized in the PDE solver package FEniCS. This is possible due to the LR-ADI method's implementation using pyMOR's abstract interfaces. The example used FEniCS, but other solver packages are also possible. Currently also deal.II, DUNE, and NGSolve are supported. Custom (domain specific) solvers can be easily integrated with pyMOR.
In pyMOR 0.5, further variants of BT for first-order systems are available, together with a few methods for second-order systems. Furthermore, interpolatory methods for first-order and second-order systems are also included.
Some PDE solver packages do not support complex numbers (e.g., FEniCS). In the example, we considered a system with symmetric E and A, where real shifts p i were sufficient. For general systems, we plan to add support to pyMOR for solving linear systems with complex shifts based on real linear algebra to further enable applicability of system-theoretic methods to models defined with such PDE solver packages.