Operator splitting for coupled linear port-Hamiltonian systems

Abstract

Operator splitting methods tailored to coupled linear port-Hamiltonian systems are developed. We present algorithms that are able to exploit scalar coupling, as well as multirate potential of these coupled systems. The obtained algorithms preserve the dissipative structure of the overall system and are convergent of second order. Numerical results for coupled mass-spring-damper chains illustrate the computational efficiency of the splitting methods compared to a straight-forward application of the implicit midpoint rule to the overall system.

[BUW]organization=Bergische Universität Wuppertal, IMACM, addressline=Gaußstraße 20, city=Wuppertal, postcode=D-42119, country=Germany

1 Introduction

Many (complex) physical systems can be modeled as linear port-Hamiltonian systems (PHS) of ordinary differential equations (ODEs) [1]

$$\dot{x}(t) =\left[J-R\right]\nabla H(x)+Bu(t),\quad x(t_{0})=x_{0},$$ (1a)

$$y(t) =B^{\top}\nabla H(x),$$ (1b)

with state variable $x\colon[t_{0},t_{\mathrm{end}}]\to\mathbb{R}^{n}$, quadratic Hamiltonian $H(x)=\tfrac{1}{2}x^{\top}Qx$, $J=-J^{\top}\in\mathbb{R}^{n\times n}$ skew-symmetric, $R\in\mathbb{R}^{n\times n}$ positive semi-definite, input $u\colon[t_{0},t_{\mathrm{end}}]\to\mathbb{R}^{d}$ and output $y\colon[t_{0},t_{\mathrm{end}}]\to\mathbb{R}^{d}$, as well as the port matrix $B\in\mathbb{R}^{n\times d}$. Supposing sufficient regularity of the input, there exists an unique solution that satisfies the energy balance condition

$$\tfrac{\mathrm{d}}{\mathrm{d}t}H(x(t))=-\nabla H(x(t))^{\top}R(x(t))\nabla H(x(t))+y(t)^{\top}u(t)\leq y(t)^{\top}u(t).$$ (2)

The dissipation inequality (2) can also be stated in integral form, i.e., $H(x(t_{0}+h))-H(x(t_{0}))\leq\int\nolimits_{t_{0}}^{t_{0}+h}y(\tau)^{\top}u(\tau)\;\mathrm{d}\tau$. We demand numerical solutions of (1) that preserve the energy balance (2) at a discrete level. This can be achieved, for example, by using discrete gradient methods [2]. Furthermore, we are aiming at deriving efficient numerical integration schemes which exploit the fine structure of the PHS. In this paper, we consider systems (1) with finer structure

$$\dot{x}(t) =\left[\!\begin{pmatrix}J_{1}&-\tilde{J}^{\top}\\ \tilde{J}&J_{2}\end{pmatrix}-\begin{pmatrix}R_{1}&0\\ 0&R_{2}\end{pmatrix}\!\right]\begin{pmatrix}Q_{1}&0\\ 0&Q_{2}\end{pmatrix}x+Bu(t),\quad x(t_{0})=x_{0},$$ (3a)

$$y(t) =B^{\top}\nabla H(x),$$ (3b)

where $J_{1}\in\mathbb{R}^{n_{1}\times n_{1}},\;J_{2}\in\mathbb{R}^{n_{2}\times n_{2}}$ are skew-symmetric, $R_{1}\in\mathbb{R}^{n_{1}\times n_{1}},\;R_{2}\in\mathbb{R}^{n_{2}\times n_{2}}$ are positive semi-definite matrices, and $Q_{1}\in\mathbb{R}^{n_{1}\times n_{1}},\;Q_{2}\in\mathbb{R}^{n_{2}\times n_{2}}$, $n=n_{1}+n_{2}$. This structure typically arises by coupling linear PHS via port matrices, e.g. in the modeling of electric circuits [3, Remark 7], and then shifting the coupling in the off-diagonal blocks of the skew-symmetric matrix.

Operator splitting techniques [4] enable the exploitation of the particular structure by additively splitting the right-hand side of (3). Recently, operator splitting techniques have been applied successfully to linear port-Hamiltonian ODE systems by splitting the PHS into its energy-conservative and energy-dissipative parts [5]. In this work, we will consider different splitting approaches that split the overall system (3) into its subsystems. Nonetheless, we want to emphasize that the resulting subsystems in our approach can be split further into their conservative and dissipative parts by applying the approaches introduced in [5] in a hierarchical manner.

The paper is organized as follows. In Section 2 we give a brief overview on operator splitting techniques and dissipativity-preserving numerical integration schemes for linear PHS. Section 3 is devoted to the special case where the matrix $\tilde{J}$ introduces a scalar coupling. Often, the port-Hamiltonian subsystems have transient behavior that is characterized by different time scales. In Section 4, we will develop multiple time stepping schemes that are integrating the subsystems using different step sizes. The proposed splitting methods are applied to coupled mass-spring-damper chains in Section 5. The paper closes with a summary and an outlook for future research.

2 Operator Splitting for Linear PHS

Operator splitting [4] is a framework to efficiently solve additively partitioned initial-value problems (IVPs) of ODEs

$$\dot{x}(t)=f(x,t)=f^{\{1\}}(x,t)+f^{\{2\}}(x,t),\quad x(t_{0})=x_{0},$$ (4)

where the right-hand side is split into two different parts w.r.t., for example, stiffness, nonlinearity, dynamical behavior and computational cost. Rewriting (4) in the homogeneous form

$$\begin{pmatrix}\dot{x}\\ \dot{s}\end{pmatrix}=\begin{pmatrix}f(x,s)\\ 1\end{pmatrix}=\underbrace{\begin{pmatrix}f^{\{1\}}(x,s)\\ 1\end{pmatrix}}_{\eqqcolon\tilde{f}^{\{1\}}(x,s)}+\underbrace{\begin{pmatrix}f^{\{2\}}(x,s)\\ 0\end{pmatrix}}_{\eqqcolon f^{\{2\}}(x,s)},\quad\begin{pmatrix}x(t_{0})\\ s(t_{0})\end{pmatrix}=\begin{pmatrix}x_{0}\\ 0\end{pmatrix},$$ (5)

operator splitting techniques approximate the exact flow $\varphi_{t}(x_{0})$ of the dynamical system (5) by composing the flows $\varphi_{t}^{\{1\}}(x_{0})$ and $\varphi_{t}^{\{2\}}(x_{0})$ of the subsystems $\dot{x}=\tilde{f}^{\{1\}}(x,s)$ and $\dot{x}=\tilde{f}^{\{2\}}(x,s)$, respectively [5]. Splitting methods of order $p>2$ require negative weights [6], which violates the dissipativity inequality and thus introduces an order barrier. The paper aims at developing efficient splitting methods of order $p=2$ that are tailored to the specific structure of the linear PHS (3). For this purpose, we consider the Strang splitting [7]

$$\Psi_{h}=\varphi_{h/2}^{\{1\}}\circ\varphi_{h}^{\{2\}}\circ\varphi_{h/2}^{\{1\}}.$$ (6)

Assuming that both subsystems are linear PHS, splitting methods satisfy the dissipativity inequality at a discrete level since the composition of dissipativity preserving maps is again dissipativity preserving.

Instead of using the exact flows $\varphi_{t}^{\{1\}}(x_{0})$ and $\varphi_{t}^{\{2\}}(x_{0})$, we can replace them by numerical approximations $\Phi_{h}^{\{1\}}(x_{0})$ and $\Phi_{h}^{\{2\}}(x_{0})$, respectively. If these are obtained by discrete gradient methods of order $p=2$, it does not affect the convergence order, as well as the energy balance of the splitting method. For linear PHS (1), the implicit midpoint rule defines a discrete gradient method of order $p=2$,

$$\begin{split}x_{1}&=x_{0}+h[J-R]Q\tfrac{x_{0}+x_{1}}{2}+hBu_{1},\\ y_{1}&=B^{\top}Q\tfrac{x_{0}+x_{1}}{2}.\end{split}$$ (7a)

with $$u_{1}=u(\tfrac{t_{0}+t_{1}}{2}) Q\tfrac{x_{0}+x_{1}}{2} \nabla H(x)=Qx H(x)=\tfrac{1}{2}x^{\top}Qx (x_{1}-x_{0})^{\top}Q\tfrac{x_{0}+x_{1}}{2}=H(x_{1})-H(x_{0})$$. In particular, the midpoint is a discrete approximation of the gradient for the quadratic Hamiltonian . As a discrete gradient, it satisfies , which implies the energy balance [8]

$$\frac{H(x_{1})-H(x_{0})}{h}=\tfrac{(x_{1}-x_{0})^{\top}}{h}Q\tfrac{x_{1}+x_{0}}{2}=\left([J-R]Q\tfrac{x_{0}+x_{1}}{2}+Bu_{1}\right)^{\top}Q\tfrac{x_{1}+x_{0}}{2}=-\left(\tfrac{x_{0}+x_{1}}{2}\right)^{\top}Q^{\top}R^{\top}Q\tfrac{x_{0}+x_{1}}{2}+u_{1}^{\top}y_{1}\leq u_{1}^{\top}y_{1}.$$

By defining $$A\coloneqq(J-R)Q$$, the update of the state variable can be written as

$$\left(I-\tfrac{h}{2}A\right)x_{1} =\left(I+\tfrac{h}{2}A\right)x_{0}+hBu_{1} \Leftrightarrow x_{1} =\left(I-\tfrac{h}{2}A\right)^{-1}\left(I+\tfrac{h}{2}A\right)x_{0}+h\left(I-\tfrac{h}{2}A\right)^{-1}Bu_{1},$$ (7b)

emphasizing its close connection to the Cayley transform. The next two sections will discuss different splittings of the partitioned linear PHS (3) exploiting the fine structure of the system, resulting in an efficient computational process.

3 Coupled Linear PHS with Scalar Coupling

In this section, we consider the linear PHS (3) as an additively partitioned ODE system (4) with splitting

$$f^{\{1\}}(x,t) =\underbrace{\begin{pmatrix}0&-\tilde{J}^{\top}\\ \tilde{J}&0\end{pmatrix}}_{A_{1}}\nabla H(x)+Bu(t), f^{\{2\}}(x,t) =\underbrace{\begin{pmatrix}J_{1}-R_{1}&0\\ 0&J_{2}-R_{2}\end{pmatrix}}_{A_{2}}\nabla H(x).$$ (8)

This splitting preserves the port-Hamiltonian structure in each subsystem, allowing the preservation of the energy balance (2) of the system using appropriate numerical integration schemes. Moreover, we can consider the subsystem $\dot{x}=f^{\{2\}}(x,t)$ as the internal evolution equation of the systems and $\dot{x}=f^{\{1\}}(x,t)$ as the coupling part.

Computational complexity. Considering the implicit midpoint rule and the computation of the update (7b) via an LU-decomposition, the splitting approach is generally more complex than the straight-forward application of the implicit midpoint rule to the overall system. The situation changes if we make further assumptions on the coupling.

Definition 1.

The system (8) is scalar coupled, if the coupling matrix $A_{1}Q$ has just one non-zero above, as well as below the diagonal.

Scalar coupling occurs in many use cases of linear PHS. We can exploit this structure to reduce the computational cost to compute (7b) for the second subsystem. As we have seen in (7b), the implicit midpoint rule demands the solution of a linear system of equations. Scalar coupling enables us to compute the transformation directly in a numerical stable way with no need of a matrix decomposition. Let ${\tilde{j}}_{1}$ be the non-zero value below and $\tilde{j}_{2}$ the non-zero value above the diagonal of $A_{1}Q$. Then, it holds

$$\begin{pmatrix}1&&&\ldots&&&0\\ &\ddots&&&&&\\ &&1&\ldots&-{\tilde{j}}_{2}&&\\ \vdots&&\vdots&\ddots&\vdots&&\vdots\\ &&-{\tilde{j}}_{1}&\ldots&1&&\\ &&&&&\ddots&\\ 0&&&\ldots&&&1\end{pmatrix}^{-1}\begin{pmatrix}1&&&\ldots&&&0\\ &\ddots&&&&&\\ &&1&\ldots&{\tilde{j}}_{2}&&\\ \vdots&&\vdots&\ddots&\vdots&&\vdots\\ &&{\tilde{j}}_{1}&\ldots&1&&\\ &&&&&\ddots&\\ 0&&&\ldots&&&1\end{pmatrix}= \begin{pmatrix}1&&&\ldots&&&0\\ &\ddots&&&&&\\ &&\frac{1+\tilde{j}_{1}\tilde{j}_{2}}{1-\tilde{j}_{1}\tilde{j}_{2}}&\ldots&\frac{2\tilde{j}_{2}}{1-\tilde{j}_{1}\tilde{j}_{2}}&&\\ \vdots&&\vdots&\ddots&\vdots&&\vdots\\ &&\frac{2\tilde{j}_{1}}{1-\tilde{j}_{1}\tilde{j}_{2}}&\ldots&\frac{1+\tilde{j}_{1}\tilde{j}_{2}}{1-\tilde{j}_{1}\tilde{j}_{2}}&&\\ &&&&&\ddots&\\ 0&&&\ldots&&&1\end{pmatrix}.$$

Applying the implicit midpoint rule (7a) to the overall system (1) demands the computation of $(I-\tfrac{h}{2}A)$ and $(I+\tfrac{h}{2}A)$, and the LU decomposition of the first matrix, resulting in setup costs of $\mathcal{O}(n^{3})$. Moreover, one time step of the implicit midpoint rule requires $\mathcal{O}(n^{2}+nd)$ operations. In contrast, the splitting method (6) with splitting (8) and assuming scalar coupling has setup costs of size $\mathcal{O}(n_{1}^{3}+n_{2}^{3})$ and demands $\mathcal{O}(nd+n_{1}^{2}+n_{2}^{2})$ operations per time step. Since it usually holds $d\ll n$, the splitting approach becomes more efficient. Particularly, a straight-forward computation shows that the splitting method requires less operations, provided that $d\leq n$ and $n_{1},n_{2}\geq 2$.

4 Coupled Linear PHS with Multirate Potential

In this section, we consider linear PHS (3), and we introduce the splitting approach

$$\dot{x} =f^{\{1\}}(x,t)+f^{\{2\}}(x,t) f^{\{1\}}(x,t) =\underbrace{\begin{pmatrix}0&-\tilde{J}^{\top}\\ \tilde{J}&J_{2}-R_{2}\end{pmatrix}}_{M_{1}}\nabla H(x)+Bu(t),\quad f^{\{2\}}(x,t)=\underbrace{\begin{pmatrix}J_{1}-R_{1}&0\\ 0&0\end{pmatrix}}_{M_{2}}\nabla H(x),$$ (9)

where $f^{\{1\}}$ is characterized by slow dynamics and $f^{\{2\}}$ by fast dynamics. Moreover, we assume $n_{1}\ll n_{2}$. To exploit this sort of multirate behavior, we aim at using a large macro-step size $h$ for the slow subsystem

$$\dot{x}=M_{1}\nabla H(x)+Bu(t),$$ (10a)

and a small micro-step size $$h/m,\ m\in\mathbb{N},\ m\gg 1,$$ for the fast subsystem

$$\dot{x}=M_{2}\nabla H(x).$$ (10b)

This can be achieved by using the impulse method, a multirate generalization to the Strang splitting [9]. It reads

$$\Psi_{h}=\varphi_{h/2}^{\{1\}}\circ\left(\varphi_{h/m}^{\{2\}}\right)^{m}\circ\varphi_{h/2}^{\{1\}}$$ (11)

and computes the numerical approximation $x_{1}\approx x(h)$ by composing (approximated) flows of the subsystems (10a) and (10b). Since the integrator (11) is a composition of flows of the subsystems (10a) and (10b) that are port-Hamiltonian and all weights are positive, the numerical approximation satisfies the dissipativity inequality.

Computational complexity. Due to the simple structure of $M_{2}$, the fast subsystem (10b) simplifies to

$$\dot{x}_{1}=\tilde{M}_{2}\nabla_{1}H(x),$$ (12)

with $\tilde{M}_{2}\coloneqq J_{1}-R_{1}$, $x_{1}$ denoting the first $n_{1}$ components and $\nabla_{1}H(x)$ denoting the first $n_{1}$ components of $\nabla H(x)$.

Neglecting the setup costs, computing one time step of the impulse method (11) with macro-step size $h$ demands $\mathcal{O}(n^{2}+mn_{1}^{2})$ operations. In contrast, computing $m$ steps of the implicit midpoint rule (7a) with step size $h/m$ to resolve the fast dynamics accurately enough demands $\mathcal{O}(mn^{2})$ operations. If the system (9) exhibits multirate potential so that $n_{1}\ll n_{2}$ and $m\gg 1$, the multiple time stepping approach requires less operations.

Remark 2.

An alternative approach would shift the coupling part into the fast subsystem. However, this would imply that the fast subsystem is of full dimension $n=n_{1}+n_{2}$ so that the $m\gg 1$ evaluations of the fast flow $\varphi^{\{2\}}$ would become more expensive, resulting in a less efficient computational process.

5 Numerical Results

A frequently used benchmark problem for linear PHS is given by a mass-spring-damper chain (MSD-chain) [10]. We consider two coupled MSD-chains, consisting of $n_{1}$ and $n_{2}$ mass-spring-damper elements, respectively. The two chains are coupled via a spring with stiffness $K_{co}>0$, see Fig. 1. The first chain consists of homogeneous masses $m_{1}>0$, mass-less springs $K_{1}>0$, and damping $r_{1}>0$. For the second chain, we have homogeneous parameters $m_{2}>0$, $K_{2}>0$, and $r_{2}>0$. Introducing positions $q_{i,1},\ldots,q_{i,n_{i}}$ and momenta $p_{i,1},\ldots,p_{i,n_{i}}$, $i=1,2$, we can formulate the equations of motion as a PHS of the form (3) by setting $q=q_{1,1}$ as the coupling variables, resulting in [3]

$$J_{1} =\begin{pmatrix}0&-1&&&&-1\\ 1&0&&&&0\\ &&\ddots&&&\vdots\\ &&&0&-1&\\ &&&1&0&0\\ 1&0&\cdots&&0&0\end{pmatrix}\in\mathbb{R}^{(2n_{1}+1)\times(2n_{1}+1)}, J_{2} =\begin{pmatrix}0&-1&&&&\\ 1&0&&&&\\ &&&\ddots&&\\ &&&&0&-1\\ &&&&1&0\\ \end{pmatrix}\in\mathbb{R}^{2n_{2}\times 2n_{2}},$$

$$R_{1} =\mathrm{diag}(r_{1},0,r_{1},0,\ldots,r_{1},0,0)\in\mathbb{R}^{(2n_{1}+1)\times(2n_{1}+1)}, R_{2} =\mathrm{diag}(r_{2},0,r_{2},0,\ldots,r_{2},0)\in\mathbb{R}^{2n_{2}\times 2n_{2}},$$

$$\tilde{J} =\begin{pmatrix}0&\cdots&0&1\\ 0&\cdots&0&0\\ \vdots&&\vdots&\vdots\\ 0&\cdots&0&0\end{pmatrix}\in\mathbb{R}^{2n_{2}\times(2n_{1}+1)}, Q =\begin{pmatrix}Q_{1}&0&0\\ 0&K_{co}&0\\ 0&0&Q_{2}\end{pmatrix}\in\mathbb{R}^{(2n+1)\times(2n+1)},$$

where the two matrices $Q_{1},Q_{2}$ are given by

$$Q_{i}=\begin{pmatrix}\tfrac{1}{m_{i}}&0&0&&&&&&\\ 0&K_{i}&0&-K_{i}&&&&&\\ 0&0&\tfrac{1}{m_{i}}&0&0&&&&\\ &-K_{i}&0&2K_{i}&0&-K_{i}&&&\\ &&\ddots&&\ddots&&\ddots&&\\ &&&\ddots&&\ddots&&0&\\ &&&&\ddots&&2K_{i}&0&-K_{i}\\ &&&&&0&0&\tfrac{1}{m_{i}}&0\\ &&&&&&-K_{i}&0&2K_{i}\\ \end{pmatrix}\in\mathbb{R}^{2n_{i}\times 2n_{i}},\quad i=1,2.$$

The state variables are $x=(p_{1,1},q_{1,1},\ldots,p_{1,n_{1}},q_{1,n_{1}},q_{1,1}-q,p_{2,1},q_{2,1},\ldots,p_{2,n_{2}},q_{2,n_{2}})^{\top}\in\mathbb{R}^{2n+1}$.

Figure 1: Coupled mass-spring-damper chain, consisting of two chains with $n_{1}$ and $n_{2}$ masses, respectively.

Figure 2: Coupled mass-spring-damper chain. Global error at $t_{\mathrm{end}}=2$ vs. the total number of operations for the implicit midpoint rule ($\circ$) and the Strang splitting ($\triangle$).

Figure 3: Coupled mass-spring-damper chain. Global error at $t_{\mathrm{end}}=2$ vs. the total number of operations for the implicit midpoint rule ($\circ$) and the impulse method with multirate factor $m=10$ ($\triangle$).

5.1 The singlerate case

In a first simulation, we choose $n_{1}=n_{2}=25$, $K_{1}=K_{2}=50$, $K_{co}=50$, $m_{1}=m_{2}=0.3$, $r_{1}=r_{2}=0.1$ and as initial value at $t=0$ we choose zero, apart from $x_{6}(0)=q_{1,3}(0)=0.1$. We perform numerical simulations on the time interval $t\in[0,2]$ using the implicit midpoint rule without any splitting approach (7a), as well as the Strang splitting (6) with splitting (8) and approximation of the flows via the implicit midpoint rule. For both approaches, we investigate the global error at $t_{\mathrm{end}}=2$ vs. the total number of operations for different step sizes $h=2^{-k},\ k=7,\ldots,11$. The results are depicted in Fig. 2, emphasizing the efficiency of the splitting method.

5.2 The multirate case

In a second example, we choose $n_{1}=5$, $n_{2}=45$, $K_{1}=100$, $K_{2}=10$, $K_{co}=10$, $m_{1}=0.1$, $m_{2}=0.4$, and $r_{1}=r_{2}=0.1$. Consequently, the masses of the first chain are characterized by faster dynamics than the masses of the second chain, introducing multirate potential that motivates the use of multiple time stepping techniques like the impulse method (11). Again, we choose zero as the initial value, apart from $x_{6}(0)=q_{1,3}(0)=0.1$ and we perform numerical simulations on the time interval $t\in[0,2]$ using the implicit midpoint rule (7a), as well as the impulse method (11) with splitting (9) and approximation of the flows via the implicit midpoint rule. Analogously to the previous example, we compare the efficiency of the splitting approach compared to a straight-forward application of the implicit midpoint rule. As multirate factor, we choose $m=10$ and perform numerical simulations using the implicit midpoint rule (7a) and the impulse method (11) with (macro) step sizes $h^{k}$, $k=9,\ldots,13$. The results in Fig. 3 highlight that the multiple time stepping approach results in a more efficient computational process.

6 Conclusion and Outlook

In this work, we developed splitting methods tailored to coupled linear port-Hamiltonian systems. The resulting algorithms are dissipativity-preserving and convergent of second order. The investigations also include multiple time stepping techniques, enabling the use of different step sizes for the subsystems. Numerical results emphasize the efficiency of the proposed splitting methods compared to a straight-forward application of the implicit midpoint rule.

Future research will deal with the derivation of more advanced decomposition algorithms for linear PHS, particularly by using the (constant) Hessian $\nabla^{2}H(x)$ of the Hamiltonian [11], enabling the derivation of fourth-order decomposition algorithms without using negative weights [12, 13] and thus satisfying the energy balance at a discrete level. We will moreover investigate the proposed splitting methods within the generalized additive Runge–Kutta (GARK) framework [14, 15] that allows to couple certain steps which may result, for example, in improved stability properties.