Structure-Preserving Model-Reduction of Dissipative Hamiltonian Systems

TL;DR

This paper introduces a structure-preserving model reduction technique for dissipative Hamiltonian systems that maintains energy dissipation and stability by coupling the original system with a heat bath, enabling accurate and stable simulations.

Contribution

It presents a novel reduction method that preserves symplectic structure and dissipation in Hamiltonian systems, improving stability and accuracy over traditional methods.

Findings

Method preserves energy dissipation and stability.

Numerical simulations demonstrate accuracy and stability.

Applicable to dissipative wave equations and electric circuits.

Abstract

Reduced basis methods are popular for approximately solving large and complex systems of differential equations. However, conventional reduced basis methods do not generally preserve conservation laws and symmetries of the full order model. Here, we present an approach for reduced model construction, that preserves the symplectic symmetry of dissipative Hamiltonian systems. The method constructs a closed reduced Hamiltonian system by coupling the full model with a canonical heat bath. This allows the reduced system to be integrated with a symplectic integrator, resulting in a correct dissipation of energy, preservation of the total energy and, ultimately, in the stability of the solution. Accuracy and stability of the method are illustrated through the numerical simulation of the dissipative wave equation and a port-Hamiltonian model of an electric circuit.

Figures

Equations

$$\dot{z}(t)$$

$$z(0)$$

$${\mathbb{J}}_{2n}=\left(\begin{array}{cc}{0}_n& I_n\\ -I_n& 0_n\end{array}\right),$$

$$\Omega (x,y)=\Lambda (\alpha (x),\alpha (y)),\text{for all }x,y\in {\mathbb{R}}^{2n}.$$

$$A^T{\mathbb{J}}_{2n}A={\mathbb{J}}_{2k},$$

$$A^{+}={\mathbb{J}}_{2k}^T{A}^T{\mathbb{J}}_{2n}.$$

$$p_{n+1/2}$$

$$q_{n+1}$$

$$p_{n+1}$$

$$\dot{z}={\mathbb{J}}_{2n}{K}^{T}f(t),$$

$$z(0)=z_0,$$

$$f(t)+{\int }_0^t\chi (t-s)\cdot f(s)ds=Kz.$$

$$\ddot{q}+r\dot{q}+kq=0$$

$$\dot{q}(t)=f(t),\dot{p}(t)=-kq(t),f(t)+{\int }_0^{t}rf(s)ds=p(t).$$

$$\dot{x}={\mathcal{J}}_{2n}\frac{\delta H_{\text{ex}}}{\delta x}.$$

$$\text{Im}(\xi \hat{\chi }(\xi ))\ge 0,\forall \xi =\omega +i\eta ,\eta \ge 0,$$

$$\hat{\chi }(\xi )={\int }_0^{\infty }{e}^{i\xi t}\chi (t)dt.$$

$$\dot{z}(t)={\mathbb{J}}_{2n}{K}^{T}f(t),$$

$${\partial }_t\varphi (t,x)=\theta (t,x),$$

$${\partial }_t\theta (t,x)={\partial }_x^2\varphi (t,x)+\sqrt{2}{\delta }_0(x)\cdot \sqrt{\chi }f(t),$$

$$z(0)=z_0,\varphi (0,\cdot )=0,\theta (0,\cdot )=0.$$

$$f(t)+\sqrt{2}\cdot \sqrt{\chi }\varphi (t,0)=Kz(t).$$

$$\varphi (t,x)=\frac{\sqrt{2}}{2}{\int }_0^{t-|x|}\sqrt{\chi }f(s)ds,\theta (t,x)=\frac{\sqrt{2}}{2}\cdot \sqrt{\chi }f(t-|x|).$$

$$H_{\text{ex}}(z,\varphi ,\theta )=\frac{1}{2}\left(\Vert Kz-\varphi (t,0){\Vert }_2^2+\Vert \theta (t){\Vert }_{{\mathcal{H}}^{2n}}^2+\Vert {\partial }_x\varphi (t){\Vert }_{{\mathcal{H}}^{2n}}^2\right)$$

$${\int }_0^{t_n}f(s)ds\approx \frac{\Delta t}{2}f(t_n)+\frac{\Delta t}{2}f(t_{n-1})+{\int }_0^{t_{n-1}}f(s)ds,$$

$${\int }_0^{t_n}f(s)ds\approx \sum _{i=0}^k{\omega }_{i}f(t_{n-i})+{\int }_0^{t_{n-k}}f(s)ds$$

$$z=Ay,$$

$$A\dot{y}={\mathbb{J}}_{2n}{\nabla }_{z}H(Ay).$$

$$\dot{y}=A^{+}{\mathbb{J}}_{2n}(A^{+}{)}^T{\nabla }_{y}H(Ay).$$

$$\dot{y}(t)$$

$$y(0)$$

$$t^{k+1}:=\underset{t\in [0,T]}{\mathrm{argmax}}\Vert z(t)-A_{2k}{A_{2k}}^{+}z(t)\Vert .$$

$$A_{2k+2}=\{e_1,\dots ,e_k,e_{k+1},f_1,\dots ,f_k,{\mathbb{J}}_{2n}^T{e}_{k+1}\}.$$

$${\omega }^{k+1}:=\underset{\omega \in \Omega }{\mathrm{argmax}}|H(z(\omega ))-H(A{A}^{+}z(\omega ))|.$$

Full text

1 Introduction

The need for increased accuracy has led to more complex models and the use of large systems of partial differential equations in engineering and science. As a consequence, direct numerical methods for solving PDEs have become computationally demanding and, at times impractical.

During the past decade, reduced basis methods have emerged as a powerful approach to reduce the cost of evaluating large systems of partial differential equations Ito:1998up ; Ito:1998ch ; Ito:2001ev . These methods construct a low-dimensional linear subspace, the reduced space, that approximately represents the solution to the system of differential equations. The projection of the original system onto the reduced space then allows the exploration of the solution with a significantly reduced computational complexity hesthaven2015certified ; quarteroni2015reduced .

Hamiltonian systems are an important class of systems that appear in engineering and science. In these systems, preserving system energy is essential to obtain a correct numerical solution. Therefore, the development of model reduction techniques that preserve the symplectic symmetry is crucial. However, the classical reduced basis methods do not generally preserve the conservation laws and intrinsic symmetries of Hamiltonian systems Amsallem:2014ef ; prajna2003pod . This often results in an unstable or a qualitatively wrong solution.

It is demonstrated in Maboudi:2016 ; Lall:2003iy ; Carlberg:2014ky ; Peng:2014di , that if the basis for the reduced space is not chosen carefully, the symplectic symmetry of Hamiltonian and Lagrangian systems will be destroyed by model reduction. To resolve this issue, in Maboudi:2016 ; Peng:2014di ; Lall:2003iy a reduced order configuration space is constructed that inherits symmetries of the full configuration space. By using a proper time integrator scheme, the symmetries are preserved in the reduced system. A greedy-type algorithem is developed in Maboudi:2016 for construction of a basis for such a reduced configuration space.

Most models in engineering appear as a dissipative perturbation of a Hamiltonian system. In these systems, conservation of energy is taken as a fundamental principle of the system dynamics, while dissipative forces, e.g. friction, can change the energy of the system vanderSchaft:2014:PST:2693645.2693646 . As the energy is no longer preserved for such systems, existing methods can no longer be applied directly peng2016geometric .

For dissipative and forced Hamiltonian system, Peng et al. peng2016geometric suggest a symplectic model reduction method that preserves the Hamiltonian and the dissipative structure of the original system. However, since this method uses a symplectic integrator for a non-conservative system, there is no guarantee that the evolution of the energy is translated correctly to the reduced system.

In the context of network modeling and circuit simulation, considerable work has been done in the development of structure preserving, and in particular energy preserving, model reduction techniques. Model reduction for port-Hamiltonian systems are given in Polyuga:2010gj ; beattie2011structure ; chaturantabut2016structure and the references therein. These methods use a Krylov or a Proper Orthogonal Decomposition (POD) approach to construct a reduced port-Hamiltonian system that preserves the passivity, and, thus, the stability of the original system. However, these methods do not generally guarantee the correct distribution of the energy among the energy consuming and energy storing units. Furthermore, over long time integration, accumulation of local errors might produce an erroneous solution.

In this paper, we present the Reduced Dissipative Hamiltonian (RDH) method as a structure preserving model reduction approach for dissipative Hamiltonian systems. A key difference between this method and the other existing methods is that the RDH enables the reduced system to be integrated using a symplectic integrator. By considering a canonical heat bath, also known as hidden strings Figotin:2006jy ; Figotin:2005 , the reduced system is extended to a closed and conservative system. Therefore, a symplectic time integrator can be used to guarantee conservation of the system energy and the correct dissipation of energy. Furthermore, the hidden strings assure that the local errors in the dissipation of energy do not accumulate, resulting in a correct evolution of the system energy.

What remains of this paper is organized as follow. Section 2 covers the required background on Hamiltonian systems, dissipative Hamiltonian systems, and the Hamiltonian extension. In Section 3 we discussion the symplectic model reduction for Hamiltonian systems and introduce the Reduced Dissipative Hamiltonian method. Accuracy, stability and efficiency of the RDH method is discussed in Section 4, and illustrated through simulation of the dissipative wave equation and a linear port-Hamiltonian system of an electrical circuit. We offer conclusive remarks in Section 5.

2 Dissipative Hamiltonian Systems

We first summarize needed concepts around the geometry of a symplectic linear vector space and then discuss Hamiltonian systems subject to dissipative forces. Finally, we introduce the Hamiltonian extension of dissipative Hamiltonian systems.

2.1 Hamiltonian Systems

Suppose that $(\mathbb{R}^{2n},\Omega)$ is a symplectic linear vector space, where $\mathbb{R}^{2n}$ is a configuration space and $\Omega:\mathbb{R}^{2n}\times\mathbb{R}^{2n}\to\mathbb{R}$ is a closed, skew-symmetric and non-degenerate 2-form on $\mathbb{R}^{2n}$. Given a smooth Hamiltonian function $H:\mathbb{R}^{2n}\to\mathbb{R}$, Hamilton’s equations of evolution are given as

[TABLE]

where $z\in\mathbb{R}^{2n}$ is the state coordinates and $\mathbb{J}_{2n}$ is a $2n\times 2n$ matrix such that $\Omega(x,y)=x^{T}\mathbb{J}_{2n}y$, for all $x,y\in\mathbb{R}^{2n}$ Marsden:2010:IMS:1965128 . By using the Symplectic Gram-Schmidt de2006symplectic method, one can construct a coordinate system in which $\mathbb{J}_{2n}$ takes the form

[TABLE]

where $I_{n}$ and $0_{n}$ are the identity matrix and the zero matrix of size $n$, respectively. A main feature of Hamiltonian systems is the conservation of the Hamiltonian along the integral curves.

Theorem 2.1

Marsden:2010:IMS:1965128 * Consider the flow $\phi_{t}:\mathbb{R}\times\mathbb{R}^{2n}\to\mathbb{R}^{2n}$ of the Hamiltonian system (1). Then $H\circ\phi_{t}=H$.*

In many physical problems, $H$ represents the system energy and is bounded from below. Here, to avoid difficulties with well-posedness of Hamiltonian systems, we often assume that $H$ is a quadratic Hamiltonian, i.e., it takes the form $H(z)=\frac{1}{2}z^{T}K^{T}Kz$, where $K$ is a full rank $2n\times 2n$ matrix. This assumption leads to a linear system of evolution (1). We emphasise that the Hamiltonian extension in Section 2.2 and the Reduced Dissipative Hamiltonian method in Section 3.3 can be naturally extended to Hamiltonians of the form $H(z)=\frac{1}{2}z^{T}K^{T}Kz+g(z)$, where $g:\mathbb{R}^{2n}\to\mathbb{R}$ is an arbitrary function of $z$.

Under general coordinate transformations, the equations of evolution may not take the form in (1). Let $(\mathbb{R}^{2n},\Omega)$ and $(\mathbb{R}^{2k},\Lambda)$ be two symplectic linear vector spaces. A linear transformation $\alpha:\mathbb{R}^{2n}\to\mathbb{R}^{2k}$ is called a symplectic transformation Marsden:2010:IMS:1965128 if

[TABLE]

In matrix notation, $A\in\mathbb{R}^{2n\times 2k}$ is called a symplectic matrix if

[TABLE]

where the superscript $T$ represents the transpose operator. As the symplectic 2-form is preserved under a symplectic transformation, the form of the equations of evolution remains invariant through a symplectic coordinate transformation Marsden:2010:IMS:1965128 . The symplectic inverse of $A$, is a pseudo-inverse given by

[TABLE]

It is shown in Peng:2014di that $A^{+}A=I_{2k}$ and that $(A^{+})^{T}$ is a symplectic matrix. Furthermore, one easily checks that $AA^{+}$ is idempotent, i.e. it is a projection operator onto the column span of $A$.

It is known that symplectic matrices are usually ill-conditioned Karow:2006cf . Under some conditions on a symplectic space da2003introduction , one can construct a symplectic basis which is also ortho-normal, and thus norm bounded. A basis which is both symplectic and orthogonal is called a ortho-symplectic basis. We refer the reader to da2003introduction for conditions on existence and construction of an ortho-symplectic basis.

It is natural to expect the numerical integrator that solves (1) also satisfy the conservation law in Theorem 2.1. However, common numerical integrators, e.g. the Runge-Kutta method, do not generally preserve the Hamiltonian. Symplectic numerical integrators are a class of numerical integrators for Hamiltonian systems that preserves the symplectic structure and ensure stability during long-time integration. The Strömer-Verlet time-stepping scheme

[TABLE]

where $p$ and $q$ are the canonical coordinates $z=(q^{T},p^{T})^{T}$, is an example of such numerical integrators. More information on the construction and application of symplectic integrators can be found in Hairer:1250576 .

2.2 Dissipative Hamiltonian Systems and Hamiltonian Extensions

Many systems in engineering and science appear as a perturbation of a Hamiltonian system, where the perturbation can be regarded as dissipation. In these systems, the energy tends to decrease over time, and thus, the conservation law in Theorem 2.1 does not hold. Therefore, it is common to take the conservation of energy as a fundamental principle and consider the dissipative system coupled with a heat bath that absorbs the dissipated energy of the original system.

To account for dissipation in a quadratic Hamiltonian $H(z)=\frac{1}{2}z^{T}K^{T}Kz$, we rewrite (1) as a time dispersive and dissipative (TDD) Figotin:2006jy system

[TABLE]

where $f$ is the solution to the Volterra integral equation corduneanu1991integral

[TABLE]

Here $\chi:\mathbb{R}^{+}\to\mathbb{R}^{2n\times 2n}$ is a bounded matrix valued function with respect to the Frobenius norm and is called the general susceptibility. Note that the integral term in (8) accounts to the accumulation of the dissipation, whereas $\chi(s)=0$ implies (7) is equivalent to (1). Furthermore, under suitable assumptions on $K$, both (1) and (7) are well-posed Figotin:2006jy .

Example 1

Consider the dynamics of the damped harmonic oscillator

[TABLE]

where $k$ is the Hooke’s constant and $r$ is the spring’s damping factor. Note that without a damping term, (9) is a Hamiltonian system. The TDD formulation for the damped harmonic oscillator takes the form

[TABLE]

Here $(q,p)$ are the canonical coordinates and the susceptibility is the constant function $r$.

It is shown in Figotin:2006jy ; Figotin:2005 that under natural assumptions on the linear susceptibility $\chi(t)$ (see below), one can couple a TDD system of the form (7) with a canonical heat bath where the dissipated energy is captured in the heat bath in a canonical sense. In other words, one can construct a Hilbert space $\mathcal{H}$ and an isometric injection $I:\mathbb{R}^{2n}\to\mathbb{R}^{2n}\times\mathcal{H}^{2n}$ where the solution $z$ to (7) is the projection of $x$ onto $\mathbb{R}^{2n}$, and $x$ is the solution to

[TABLE]

Here $H_{\text{ex}}:\mathbb{R}^{2n}\times\mathcal{H}^{2n}\to\mathbb{R}$ is an extended quadratic Hamiltonian function and $\mathcal{J}_{2n}$ is the symplectic operator defined on $\mathbb{R}^{2n}\times\mathcal{H}^{2n}$ respectively.

Theorem 2.2

Suppose that $K$ is full rank and $\chi(t)$ is symmetric. Then there is a quadratic extension to (7) of the form (11), if

[TABLE]

where $\hat{\chi}$ is the Fourier-Laplace transform of $\chi$

[TABLE]

Proof. Here we prove the theorem for the case where $\chi$ is a constant symmetric matrix, where condition (12) corresponds to $\chi$ being positive semi-definite. We refer the reader to Figotin:2006jy for the proof of the general case. Consider the Hamiltonian system

[TABLE]

together with the initial condition

[TABLE]

Here $\theta$ and $\phi$ are vector valued functions in $\mathcal{H}^{2n}$, $\delta_{0}(s)$ is the Dirac’s delta function, $\sqrt{\chi}$ is the matrix square root of $\chi$ and $f$ is the solution to the equation

[TABLE]

To show that the Hamiltonian system (14a)-(14b) is an extension to (7) in the sense discussed above, it is enough to show that the solution $f$ to equation (16) also satisfies (8). Equations (14b) and (14c) are equations for a vibrating string, and can be solved analytically

[TABLE]

Then, we recover (8) by substituting (17) into (16). The extended Hamiltonian $H_{\text{ex}}$ for the system (14a)-(14b) takes the quadratic from

[TABLE]

where $\|\cdot\|_{2}$ is the Euclidean norm on $\mathbb{R}^{2n}$ and $\|\cdot\|_{\mathcal{H}^{2n}}$ is the induced norm from the inner product on $\mathcal{H}^{2n}$.

Equations (14b) and (14c) are called the hidden strings. The dissipation of energy in the original system (7) is carried away, as vibrations, along the added strings making the extended system conservative. The Hamiltonian extension of the damped harmonic oscillator in Example 1 is exactly the Lamb model lamb:1900 which is a harmonic oscillator coupled with a vibrating string, and the tension in the string causes linear dissipation in the dynamics of the harmonic oscillator.

Note that the time integration of (14a)-(14b) involves the integration of $f$ in (17). In general, the history of $f(t)$ must be stored and may cause storage limitation in long-time integration. However, we are interested solely in finding $z(t)$ which depends on $f$ at time $t$, and $\phi(t,0)$, i.e. the integral of the history of $f$. So by carefully choosing a quadrature rule that uses the same quadrature nodes as the time integrator we can avoid storing the history of $f$. For example for the trapezoidal rule, we recover the recursive relation

[TABLE]

where $\Delta t$ is the time step. The recursive relation in (19) suggests that storing the value of the integral term together with the state of $f$ in the previous time step suffices to evaluate the integral for the new time step. For other interpolation based quadrature rules, we can construct similar recursive rules of the form

[TABLE]

for some quadrature weights $\omega_{i}$, $i=1,\dots,k$ with $k\ll n$. Thus, time integration of (14a)-(14b), only requires storage of $k$ evaluations of $f$.

3 Model Order Reduction

In this section we first explain the main results of Maboudi:2016 ; Peng:2014di regarding model reduction of Hamiltonian systems and, subsequently we introduce the Reduced Dissipative Hamiltonian method.

3.1 Symplectic Model Order Reduction

Consider a Hamiltonian system of the form (1) together with a quadratic Hamiltonian of the form $H(z)=\frac{1}{2}z^{T}K^{T}Kz$. In this paper we focus on reducing the complexity of the numerical evaluation of (1) with respect to time $t$. Nevertheless, one can extend the results of this paper to a Hamiltonian system that depends on a set of physical or geometrical parameters belonging to a compact subset of a Euclidean space.

The main idea behind model order reduction is that the solution manifold, $\mathcal{M}=\{z(t):t\in[0,T]\}$ can be approximated by a low dimensional linear subspace hesthaven2015certified ; quarteroni2015reduced . A basis for such a subspace is called a reduced basis, and its span is referred to as the reduced space hesthaven2015certified .

Suppose that a reduced symplectic basis $A\in\mathbb{R}^{2n\times 2k}$ is provided with $k\ll n$. The approximate solution to (1) in this basis is expressed as

[TABLE]

where $y\in\mathbb{R}^{2k}$ are the coordinates of the approximation with respect to the basis $A$. Substituting (21) into (1) yields

[TABLE]

Multiplying both sides with the symplectic inverse of $A$ and using the chain rule we obtain

[TABLE]

As $(A^{+})^{T}$ is a symplectic matrix it implies that $A^{+}\mathbb{J}_{2n}(A^{+})^{T}=\mathbb{J}_{2k}$. By defining the reduced Hamiltonian $\tilde{H}:\mathbb{R}^{2k}\to\mathbb{R}$, as $\tilde{H}(y)=H(Ay)$, we recover the reduced system

[TABLE]

Equation (24) is called the symplectic Galerkin projection Peng:2014di of the Hamiltonian system (1). Conventional model reduction techniques, e.g. Galerkin or Petrov-Galerkin methods hesthaven2015certified ; quarteroni2015reduced , do not yield a Hamiltonian reduced system and the reduced system does not necessarily preserve the conservation law in Theorem 2.1 which results in a qualitatively wrong and often unstable solution Peng:2014di . On the other hand the reduced system, obtained by the symplectic Galerkin projection, is a Hamiltonian system, and the system energy is therefore preserved over time Peng:2014di . The following theorem guarantees the boundedness of the solution of the reduced system obtained by the symplectic Galerkin projection for quadratic Hamiltonians.

Theorem 3.1

Peng:2014di * Let $S$ be a bounded open subset of $\mathbb{R}^{2n}$ that contain $z_{0}$. Furthermore, assume that $H(z_{0})<H(z)$ or $H(z_{0})>H(z)$ for all $z\in\partial S$, the boundary of $S$. For a given symplectic basis $A$, provided $z_{0}$ is in the range of $A$, then both the original system and the reduced system obtained by the symplectic Galerkin projection are bounded.*

In the next section we introduced the greedy method to construct a symplectic reduced basis Maboudi:2016 .

3.2 The Greedy Approach to Symplectic Basis Generation

Greedy generation of the reduced basis is an iterative procedure which, in each iteration, adds the two best possible basis vectors to the symplectic basis to enhance overall accuracy. Suppose that at the $k$-th generic step of the greedy basis selection, an ortho-symplectic basis $A_{2k}=\{e_{1},\dots,e_{k},f_{1},\dots,f_{k}\}$ of size $2k$ is provided. The first step in an iteration of the greedy basis selection comprises finding $t^{k+1}\in[0,T]$ such $z(t^{k+1})$ is worst approximated by the current reduced space. In other words

[TABLE]

Suppose that $e_{k+1}$ is the vector obtained by $\mathbb{J}_{2n}$-orthogonalization (the symplectic Gram-Schmidt process Salam2014 ) of $z(t^{k+1})$ with respect to $A_{2k}$. Then the enriched basis $A_{2k+2}$ takes the form

[TABLE]

It is easily checked that $A_{2k+2}$ is ortho-symplectic. Furthermore, it is shown in Maboudi:2016 that under natural assumptions on the solution manifold $\mathcal{M}$, the greedy method converges exponentially fast.

For parametric problems, evaluation of the projection error in (25) for the entire parameter space is computationally demanding. One can use the loss in the Hamiltonian function as a cheap surrogate to the projection error, i.e.

[TABLE]

Here $\Omega$ is a closed and bounded set of parameters for the original Hamiltonian system. Note that $\omega^{k+1}$ can be identified prior to time integration since the Hamiltonian function is constant in time, i.e. $H(z)=H(z_{0})$ and $H(AA^{+}z)=H(AA^{+}z_{0})$. We summarize the greedy method for symplectic basis generation in Algorithm 1 and refer the reader to Maboudi:2016 for further details.

3.3 The Reduced Dissipative Hamiltonian Method

Since the symplectic model reduction in Section 3.2 is based on the conservation law in Theorem 2.1, it can no longer be applied to dissipative Hamiltonian systems. Instead in the Reduced Dissipative Hamiltonian method, we considers a Hamiltonian extension to a dissipative Hamiltonian system to construct a closed system. A symplectic model reduction can then be naturally applied to conserve the total energy.

Consider a dissipative Hamiltonian system of the form (7) with a quadratic Hamiltonian, $H(z)=z^{T}K^{T}Kz$. Since $K^{T}K$ is symmetric and positive definite, it has a unique Cholesky factorization $K^{T}K=L^{T}L$ where $L$ is upper triangular strang09 . So we can write

[TABLE]

Further, suppose that the solution $z(t)$ lies on a low-dimensional symplectic subspace such that $z=Ay$, where $A$ is an ortho-symplectic matrix of the size $2n\times 2k$ and $y$ is the expansion coefficients of $z$ in the basis of $A$. The Hamiltonian $H$ then takes the form

[TABLE]

Since $A^{T}L^{T}LA$ is symmetric and positive definite, it too has a unique Cholesky factorization $\tilde{L}^{T}\tilde{L}$ where $\tilde{L}=A^{T}LA$ is an upper triangular matrix of size $2k\times 2k$. Writing (7) in terms of the reduced coordinates $y$ and the Hamiltonian (28) reads

[TABLE]

together with the complementary equation

[TABLE]

Multiplying (30) with $A^{+}$ and (31) with $A^{T}$ yields

[TABLE]

where we use the fact that $A^{+}\mathbb{J}_{2n}=\mathbb{J}_{2k}A^{T}$. If we define $f=A\tilde{f}$, $\phi=A\tilde{\phi}$, $\theta=A\tilde{\theta}$ and the reduced susceptibility as $\tilde{\chi}=A^{T}\chi A$ we recover the reduced Hamiltonian system

[TABLE]

together with the auxiliary equation

[TABLE]

Equations (34a)-(34c) is a Hamiltonian system on the symplectic linear vector space $\mathbb{R}^{2k}\times\mathcal{H}^{2k}$ and contributes to the reduced TDD system

[TABLE]

Therefore, the system energy will be conserved along integral curves of (34a)-(34c).

We point out that the transformation that connects (14a)-(14c) to (34a)-(34c) is given by

[TABLE]

This is a symplectic transformation, since $\mathcal{A}^{T}\mathcal{J}_{2n}\mathcal{A}=\mathcal{J}_{2k}$. Furthermore, the dissipation of energy in the reduced system only depends on the reduced susceptibility. Thus, the choice of $A$ should be independent of the hidden strings $(\phi,\theta)$. In other words, if the reduced space is chosen to be a symplectic subspace, then the actions of model reduction and Hamiltonian extension commute. We summarize the algorithm for model reduction of dissipative Hamiltonian systems in Algorithm 2.

Note that Algorithm 2 does not depend on the choice of the method to construct an ortho-symplectic basis $A$. Thus, for basis generation, the greedy approach introduced in section 3.2 or an SVD-based method, e.g. the cotangent lift Peng:2014di , can be applied.

The main advantage of the RDH method compared to the existing methods is that it enables the reduced system to be integrated using a symplectic integrator. The reduced system constructed using the RDH is a closed Hamiltonian system, therefore the conservation law in Theorem 2.1 holds and a symplectic integrator guarantees that the total energy is preserved in the reduced system. Alternative methods, e.g. peng2016geometric ; Polyuga:2010gj ; beattie2011structure , either integrate the reduced system with a non-symplectic integrator, or do not construct a closed reduced system which result in accumulation of local errors or unstable solution during long time integration, respectively Hairer:1250576 .

4 Numerical Results

In the following we illustrate the performance of the method through the reduced order model of the dissipative wave equation and a port-Hamiltonian model for a dissipative circuit.

4.1 Dissipative wave equation

Consider the dissipative linear wave equation

[TABLE]

where $x$ belongs to a one-dimensional torus of length $L$ and $r:[0,1]\to[0,1]$ is a positive semi-definite real valued function.

We discretize the torus into $N_{\Delta x}$ equidistant points and define $\Delta x=L/N_{\Delta x}$, $x_{i}=i\Delta x$, $q_{i}=q(t,x_{i})$ and $p_{i}=p(t,x_{i})$ for $i=1,\dots,N_{\Delta x}$. The discretization of $r$ corresponds to a diagonal and semi-positive definite matrix $r_{\Delta}$. Furthermore, we discretize (38) using a standard central finite differences schemes to obtain

[TABLE]

where $z=(q_{1},\dots,q_{N_{\Delta x}},p_{1},\dots,p_{N_{\Delta x}})$ and $K$ and $R$ are given as

[TABLE]

with $D_{x}^{T}D_{x}=D_{xx}$ as the central finite differences matrix operator. Writing (39) in a TDD formulation yields

[TABLE]

Since $R$ is not time dependent, it commutes with the integration operator. The Hamiltonian extension of (41), then takes the form (14a)-(14c).

The initial condition used is given by

[TABLE]

where $h(s)$ is the cubic spline function

[TABLE]

For the numerical time integration of the extended Hamiltonian system, the Strömer-Verlet time stepping scheme (6) is used. In each time step, the system of linear equations (16) is solved to recover $z$. System parameters are summarized below.

[TABLE]

The first numerical experiment corresponds to an inhomogeneous dissipative media. Here, $r_{\Delta}$ is a diagonal matrix with diagonal elements $r_{i}:=0.1+0.9(i/N_{\Delta x})$, for $i=1,\dots,N_{\Delta x}$.

Figure 1.(a) shows the solution of the original dissipative wave equation (38) at $t\in\{0,2.5,5,7.5\}$. For a nonzero $r_{\Delta}$ the solution will converge to $(q(t=\infty,x),p(t=\infty,x))=(\rho,0)$ where $\rho$ is the center of mass of $q_{0}$.

We construct the RDH reduced system according to the Algorithm 2. The performance of the method is then compared to the POD and the method proposed in peng2016geometric , referred to as the PSD.

Figure 1.(b) illustrates the decay of the singular values of the snapshot matrix hesthaven2015certified , for the POD, PSD, and the RDH methods. Note that the snapshots for the PSD and the RDH are different since they have different canonical representations. The fast decay of the eigenvalues in all methods is a strong indicator for the existence of a low dimensional reduced system. The reduced bases are then constructed using 20, 40 and 60 number of modes.

The $L^{2}$-error between the full system and the RDH, the PSD, and the POD methods are presented in Figure 1.(c). We notice that the symplectic methods provide a more accurate solution when compared to the POD method. In fact, the POD method does not yield a stable reduced system. Furthermore, it is seen that enriching the PSD reduced basis does not yield a significant improvement in the accuracy of the reduced system. This happens as the PSD method, numerically integrate a non-conservative system with a symplectic integrator. This results in an incorrect evolution of the energy and eventually, in a qualitatively wrong numerical solution.

On the other hand, we notice that the RDH method with 40 modes provides a significantly more accurate solution compared to the PSD method with 60 modes. The RDH method provides a conservative reduced system where the dissipated energy is absorbed by the hidden strings and the conservation of the energy is then guaranteed by using a symplectic integrator. Therefore, we observe remarkable increase in the accuracy by enriching the RDH reduced basis.

Figure 1.(d) shows the conservation of the energy in the different methods. The conservation law expressed in Theorem 2.1 is destroyed through the POD model reduction and as a consequence we observe blow-up of the system energy. The symplectic methods preserves the energy significantly better. As discussed above, enriching the PSD basis does not significantly improve the preservation of energy. On the contrary, the RDH provides a substantial improvement in the accuracy of the energy.

In Figure 1.(e) we show the transfer of the energy from the TDD system to the hidden strings, for the full system and the RDH reduced system. We notice that the RDH method preserves the total energy of the extended Hamiltonian system. Furthermore, the transfer of energy to the hidden strings in the full model is correctly translated in the reduced system.

The second numerical experiment is the dissipative wave equation (38) in a near-zero dissipation regime. The numerical setting is taken to be identical to the previous numerical experiment, but with the difference that $r_{i}=10^{-5}$, for $i=1,\dots,N_{\Delta x}$. Figure 1.(f) shows the $L^{2}$-error between the solution to the reduced system and the full system, for the POD, the PSD, and the RDH methods. We notice that the POD does not yield a stable reduced system as the symplectic structure is lost via model reduction. Furthermore we notice that error for the PSD and the RDH coincide as the two methods become identical as $\|\chi\|_{\infty}\to 0$.

4.2 The sine-Gordon equation

Consider the one-dimensional dissipative nonlinear wave equation

[TABLE]

defined on a domain of length $L$, which is known as the sine-Gordon equation. In the absence of dissipation, $r(x)=0$, the kink solution to (44) is given as

[TABLE]

where $|v|<1$ is the wave speed. In the presence of dissipation, where $r(x)\geq 0$, the traveling wave de-accelerates and stops. The TDD formulation for (44) takes the form

[TABLE]

where $z$, $K$, $R$ and $f$ are defined similar to (41), with $r_{\Delta}=rI_{n}$, and

[TABLE]

and $z_{\text{bd}}$ is the term corresponding to the Dirichlet boundary condition. Note that the extended Hamiltonian $H_{ex}$ takes the form

[TABLE]

where $G(z)$ is a potential for $g(z)$ given as $G(q_{i})=1-\sin(q_{i})$, for $i=1,\dots,N_{\Delta_{x}}$. System parameters are summarized below

[TABLE]

The RDH reduced system is constructed following Algorithm 2. To reduce the complexity of the nonlinear term, we used the symplectic discrete empirical interpolation method (SDEIM) Maboudi:2016 . The performance of the method is then compared to the PSD and the POD where the SDEIM proposed in Peng:2014di and the classical DEIM Chaturantabut:2010cz is applied to reduce the complexity of the nonlinear term, respectively.

Figure 2.(a) shows the The $L^{2}$-error between the full system and the RDH, the PSD, and the POD methods. Although the Hamiltonian system of the sine-Gordon equation is nonlinear, the errors for the different methods show a similar behavior as those in Section (4.1). We observe that the POD does not yield a stable reduce system while the symplectic methods provide a high accuracy. Furthermore, we notice that enriching the PSD basis does not significantly enhance the accuracy of the method.

The evolution of error in the kinetic energy $K(p)=\|p\|_{2}^{2}/2$ is illustrated in Figure 2.(b). We see that the POD does not conserve the evolution of the kinetic energy. The RDH method conserves the kinetic energy of the system with a higher accuracy than the PSD method. Furthermore, the accuracy of the RDH method is better scaled under enrichment of the reduced basis, compared to the PSD method.

It is observed in Figure 1 that the symplectic treatment of the nonlinear terms is essential in correct model reduction of Hamiltonian systems. In addition, the SDEIM can be combined with the RDH method to construct a reduced Hamiltonian system that can be integrated using a symplectic integrator. Thus, the combination preserves the system energy and the symplectic symmetry of Hamiltonian systems.

4.3 Port-Hamiltonian Systems

Port-Hamiltonian systems are popular in network modeling and electrical engineering. In the framework of port-Hamiltonian modelling, energy conservation and Hamiltonian structure is the fundamental principle of the dynamics of the system. Ports in the system network then allows the exchange of energy with the environment in the form of sources, capacitors, and dissipations vanderSchaft:2014:PST:2693645.2693646 . Port-Hamiltonian systems can be viewed as a forced and dissipative Hamiltonian system.

Consider the $n$-dimensional linear ladder network in Figure 3. Here $C_{i}$, $L_{i}$ and $R_{i}$, $i=1,\dots,n$, are the capacitance, inductance and resistance of the corresponding capacitors, inductors, and resistors, respectively, and $R_{n+1}$ is the load capacitor. The port-Hamiltonian model of the linear ladder network takes the form

[TABLE]

Here $x=(c_{1},\phi_{1},\dots,c_{n},\phi_{n})^{T}$ where $c_{i}$ and $\phi_{i}$, for $i=1,\dots,n$, are the charge and the flux of $C_{i}$ and $L_{i}$ respectively. $Q$ and $R$ are given as

[TABLE]

$u=(1,0,\dots,0)^{T}$ is a constant input current and $J_{2n}$ is a skew-symmetric $2n\times 2n$ matrix with -1 and 1 on the superdiagonal and subdiagonal, respectively.

The energy associated with a port-Hamiltonian system of the form (49) at time $t$, is given as $H(x(t))=\frac{1}{2}x^{T}Q^{T}Qx$. Since $J_{2n}$ is skew symmetric we have that that $\frac{d}{dt}H(x)=u^{T}Q^{T}Qx-x^{T}Q^{T}QRQ^{T}Qx\leq u^{T}Q^{T}Qx$ which is referred to as the passivity of the system (49) vanderSchaft:1996es ; Willems:1972ek .

Since $J_{2n}$ is full rank, one can always find a coordinate transformation $x=T\tilde{x}$ such that $T^{-1}J_{2n}T^{-T}=\mathbb{J}_{2n}$. The dissipative Hamiltonian formulation of (49) takes the form

[TABLE]

where $\tilde{Q}=QT$, $\tilde{R}=T^{-1}RT^{-T}Q^{T}Q$ and $\tilde{u}=T^{-1}u$. Note that in this case, $\tilde{R}$ is symmetric and semi-positive definite since $T$ is orthogonal and $R$ is diagonal. The TDD formulation of (51) takes the form

[TABLE]

The extended Hamiltonian formulation (14a)-(14c) with a quadratic Hamiltonian $H_{\text{ex}}$ can be carried out following Section 2.2. We note that due to the input $\tilde{u}$, the Hamiltonian $H_{\text{ex}}$ is time dependent. In fact $\frac{d}{dt}H_{\text{ex}}=\tilde{u}^{T}\tilde{Q}^{T}\tilde{Q}\tilde{x}$. If we define $\overset{\circ}{H}_{\text{ex}}:\mathbb{R}^{2n}\times\mathcal{H}^{2n}\times\mathbb{R}^{2}\to\mathbb{R}$ as

[TABLE]

it is easily checked that $\frac{d}{dt}\overset{\circ}{H}_{\text{ex}}=0$ Hairer:1250576 . However for the time integration of the Hamiltonian system related to $\overset{\circ}{H}_{\text{ex}}$ we can apply a symplectic integrator directly on (53), since the evolution of $\tilde{x}$, $\phi$ and $\theta$ does not explicitly depend on $e$. Thus, the passivity of (49) will be preserved through a symplectic time integration of (52).

Using an ortho-symplectic reduced basis $A$, the Reduced Dissipative Hamiltonian method can be applied to (52) to construct a reduced system of the form (34a)-(34c) together with the extended Hamiltonian $\tilde{H}_{\text{ex}}$. We note that $\frac{d}{dt}\tilde{H}_{\text{ex}}=(A^{+}\tilde{u})^{T}A^{T}\tilde{Q}^{T}\tilde{Q}Ay$, ensuring that the reduced system is passive. Furthermore, the dissipative Hamiltonian structure of the reduced system indicates that the reduced system also carries a port-Hamiltonian structure.

We consider a 100-dimensional ($n=50$) port-Hamiltonian system for the ladder network discussed above. We take $C_{i}=1$, $L_{i}=1$, $R_{i}=0.2$ for $i=1,\dots,50$, and $R_{51}=0.4$. We construct the RDH reduced system following Algorithm 2.

The solution of the RDH method is compared to the main results of Polyuga:2010gj , where a passivity-preserving model reduction is developed using a moment matching method at infinity. The charge in $C_{1}$ is chosen to be the single out put for the moment matching method.

Reduced bases of size $2k=10$, $2k=20$ and $2k=30$ are constructed with the RDH and the moment matching method. Figure 4.(c) shows the error in the charge in $C_{1}$ for the two methods. We observe that although the moment matching method is bounded over long-time integration, the RDH method provides a significantly more accurate solution. In the moment matching method, the passivity of the reduced system implies that the energy of the system will be bounded by the input energy. However, there is no guarantee that the dissipation of energy in the reduced system mimics the one of the original system. On the other hand, the RDH method allows a correct dissipation of energy through the hidden strings and the symplectic time integration in the RDH method guarantees that the total energy is preserved.

Over short-time integration, we notice that the moment matching method with 10 modes provides a more accurate solution than the RDH with 10 modes. Furthermore, the moment matching method with 20 and 30 modes provide a comparable accuracy to the RDH method with 20 and 30 modes. However, the RDH method maintains the high accuracy during long-time integration, while the moment matching method loses up to 3 orders of magnitude in the accuracy, independent of the number of modes.

Figure 4.(a) and Figure 4.(b) show the average temporal error in the charge and flux of the capacitors and inductors, respectively. The RDH method provides a significantly better accuracy compared to the moment matching method. This is because the charge of $C_{1}$ is specified as the output of interest in the moment matching method and so it is expected that that method provides low accuracy for computing other outputs. On the other hand, the RDH method not only provides high accuracy in computing the charge for $C_{1}$ but also high accuracy for all components of the system.

5 Conclusion

In this paper we present the Reduced Dissipative Hamiltonian method. The method preserves the symplectic structure of dissipative Hamiltonian systems and guarantees the correct dissipation of energy through time integration. The RDH method couples the reduced system with a canonical heat bath such that the reduced system forms a closed system.

The main advantage of the RDH method compared to the existing methods is that it enables the reduced system to be integrated using a symplectic integrator which naturally preserves the Hamiltonian structure and the symplectic symmetry of the Hamiltonian systems. Applying a symplectic integrator to a non-conservative system or using a non-symplectic integrator for the reduced system can cause accumulation of local errors or wrong qualitative solution over long-time integration, respectively.

The numerical simulations illustrate that the RDH method preserves the system energy with significantly higher accuracy than other methods. Furthermore, it is shown that the hidden strings assure that the dissipation of energy in the reduce system mimics the dissipation of energy in the full system. This ensures that the local error do not accumulate over long-time integration.