Discrete-time port-Hamiltonian systems based on Gauss-Legendre collocation⁎

6th IFAC Workshop on Lagrangian and Hamiltonian Methods or 6th IFAC Workshop on Lagrangian and Hamiltonian Methods or Nonlinear Control on Lagrangian and Hamiltonian Methods or 6th IFAC Workshop Nonlinear Control Available online at www.sciencedirect.com Universidad Técnicaon Federico Santa María 6th IFAC Workshop Lagrangian and Hamiltonian Methods or Nonlinear Control Universidad Técnica Federico Santa María Valparaíso, Chile, May 1-4, 2018 Universidad Técnica Federico Santa María Nonlinear Control Valparaíso, Chile, May 1-4, 2018 Universidad Técnica Federico Santa María Valparaíso, Chile, May 1-4, 2018 Valparaíso, Chile, May 1-4, 2018

ScienceDirect

IFAC PapersOnLine 51-3 (2018) 125–130

Discrete-time port-Hamiltonian systems Discrete-time port-Hamiltonian systems Discrete-time port-Hamiltonian systems  based on Gauss-Legendre collocation Discrete-time port-Hamiltonian systems based on Gauss-Legendre collocation based on Gauss-Legendre collocation  based on Gauss-Legendre collocation Paul Kotyczka ∗∗ , Laurent Lef` evre ∗∗ ∗∗

Paul Kotyczka ∗∗ , Laurent Lef` evre ∗∗ Paul Kotyczka , Laurent Lef` evre ∗∗ ∗ ∗∗ ∗ Paul Kotyczka , Laurent Lef` evre University of Munich, Department of Mechanical ∗ Technical Technical University of Munich, Department of Mechanical ∗ ∗ Engineering, Chair of Automatic Control, Boltzmannstraße 15, Technical University of Munich, Department of Mechanical Engineering, Chair of Automatic Control, Boltzmannstraße 15, 85748 85748 ∗ Technical University of Munich, Department of Mechanical Garching, Germany (e-mail: [email protected]) Engineering, Chair of Automatic Control, Boltzmannstraße 15, 85748 Garching, Germany (e-mail: [email protected]) ∗∗ Engineering, ChairAlpes, ofGermany Automatic Control, Boltzmannstraße 15, 85748 Univ. Grenoble LCIS, rue eel´ Garching, (e-mail: [email protected]) ∗∗ Grenoble Alpes, LCIS, 50 50 rue Barth´ Barth´ l´eemy my de de Laff´ Laff´eemas, mas, 26902 26902 ∗∗ ∗∗ Univ. Garching, Germany (e-mail: [email protected]) Univ. Grenoble Alpes, LCIS, 50 rue Barth´ e l´ e my de Laff´ e mas, 26902 Valence, France (e-mail: [email protected]) Valence, France (e-mail: [email protected]) ∗∗ Univ. Grenoble Alpes, LCIS,[email protected]) 50 rue Barth´el´emy de Laff´emas, 26902 Valence, France (e-mail: Valence, France (e-mail: [email protected]) Abstract: We introduce Abstract: We introduce a a family family of of discrete-time discrete-time lossless lossless input-state-output input-state-output port-Hamiltonian port-Hamiltonian systems based on numerical time integration with symplectic collocation For Abstract: We introduce a family of discrete-time lossless input-state-output port-Hamiltonian systems based on numerical time integration with lossless symplectic collocation schemes. schemes. For systems systems Abstract: We introduce a family of discrete-time input-state-output port-Hamiltonian with non-zero input, symplecticity extends to the conservation of a discrete energy systems based on numerical time integration with symplectic collocation schemes. For balance, systems with non-zero input, symplecticity extends to the conservation of a discrete energy balance, systems based on numerical timeDirac integration with symplectic collocation schemes. For systems with non-zero input, symplecticity extends to the conservation ofGauss-Legendre a discrete energy balance, based on which a discrete-time structure is defined. Using collocation, based on whichinput, a discrete-time Dirac structure is defined. UsingofGauss-Legendre collocation, with non-zero symplecticity extends to the conservation a discrete energy balance, based on which a quadrature discrete-time Dirac structure is defined. Using Gauss-Legendre collocation, the corresponding formula allows to quantify the discretization error for the supplied the corresponding formula to quantify the discretization error for the supplied based on which a quadrature discrete-time Dirac allows structure is defined. Using Gauss-Legendre collocation, energy. On a linear example, backward error analysis and numerical experiments are performed the corresponding quadrature formula allows to quantify the discretization error for the supplied energy. On a linear example, backward error analysis and numerical experiments are performed the corresponding quadrature formula allows to quantify the discretization error for the supplied in order to illustrate the accuracy of the resulting structure-preserving integration schemes. energy. On a linear example, backward error analysis and numerical experiments are performed in orderOn to illustrate the accuracy of theerror resulting structure-preserving integration energy. linear example, backward analysis and numerical experiments areschemes. performed in order to aillustrate the accuracy of the resulting structure-preserving integration schemes. © order 2018, IFAC (International Federationofofthe Automatic Control) Hosting by Elsevier Ltd. All rights reserved. in to illustrate the accuracy resulting structure-preserving integration schemes. Keywords: Keywords: Port-Hamiltonian Port-Hamiltonian systems, systems, Dirac Dirac structures, structures, discrete-time discrete-time systems, systems, geometric geometric Keywords: integration, Port-Hamiltonian systems, Dirac structures, discrete-time systems, geometric numerical symplectic methods. numerical symplectic methods. Keywords: integration, Port-Hamiltonian systems, Dirac structures, discrete-time systems, geometric numerical integration, symplectic methods. numerical integration, symplectic methods. 1. INTRODUCTION on on discrete discrete manifolds manifolds (spaces (spaces that that locally locally look look like like discretidiscreti1. 1. INTRODUCTION INTRODUCTION zation grids or the set of floating-point numbers). Objects on discrete manifolds (spaces that locally look like discretization grids or the set of floating-point numbers). Objects The geometric geometric integration integration of ordinary ordinary differential differential equatiequati- on 1. INTRODUCTION discrete manifolds (spaces thatgeometry locallynumbers). look like discretiand operations from differential are adapted to zation grids or the set of floating-point Objects The of and operations from differential geometry are adapted to The integrationand of ordinary differential ons, geometric see e. e. g. g. Leimkuhler Leimkuhler Reich (2004), (2004), Hairerequatiet al. al. zation grids setting or the set ofdiscrete-time floating-point numbers). Objects and operations from differential geometry are adapted to the discrete and Dirac structures are ons, see and Reich Hairer et the discrete setting and discrete-time Dirac structures are The geometric integration of ordinary differential equations, see e. g. Leimkuhler and Reich (2004), Hairer et al. (2006) is an important approach to perform long-time siand operations from differential geometry are adapted to defined. In the discrete setting, the chain rule is not valid, the discrete setting and discrete-time Dirac structures are (2006) approach to perform long-time siIn the discrete setting, the chain rule is not valid, ons, seeis g. important Leimkuhler and Reich (2004), Hairer et al. (2006) is e.an an important approach toSymplectic perform long-time si- defined. mulations of Hamiltonian systems. integration the discrete setting andchange discrete-time Dirac structures are which means that the of energy over a sampling defined. In the discrete setting, the chain rule is not valid, mulations of Hamiltonian systems. Symplectic integration means the change energy over isa not sampling (2006) is an approach toSymplectic perform long-time si- which mulations of important Hamiltonian systems. integration conserves not only the the symplectic symplectic form in the the (mechanical) (mechanical) defined. Inonly thethat discrete setting,byof the chain rule valid, which means that the change of over a sampling interval is approximated aaenergy product of the discrete conserves not only form in interval is only approximated by product of the discrete mulations of Hamiltonian systems. Symplectic integration conserves not only the symplectic form in the (mechanical) phase space, but also invariants of motion (Casimirs, first which means that the change of energy over a sampling interval is only approximated by a product of the discrete gradient  H(x) and the increment ∆x of the state. phase space, but of motion (Casimirs, first xx H(x) and the increment ∆x of the state. conserves not onlyalso the invariants symplectic inare the (mechanical) phase space, but also invariants ofform motion (Casimirs, first gradient integrals). Symplectic integrators that derived based interval approximated by a product of the discrete gradientisxxonly H(x) and the increment ∆x of the state. integrals). Symplectic integrators that are derived based give a new definition of discrete-time Dirac structures phase space, but also invariants of that motion first We integrals). Symplectic integrators are(Casimirs, derived based on discrete versions of Hamilton’s principle are called gradient  H(x) and the increment ∆x of the We give a xnew definition of discrete-time Dirac state. structures on discrete versions of Hamilton’s principle are called We give a new definition ofPH discrete-time Diracisstructures lossless discrete-time systems, based integrals). are based on discreteSymplectic versions ofintegrators Hamilton’s principle variational integrators. They “workthat very wellderived forare bothcalled con- and and lossless discrete-time PH systems, which which isstructures based on on variational integrators. They “work very well for both conWe give a new definition of discrete-time Dirac 1 symplectic collocation methods for the geometric numerical on discrete versions of Hamilton’s principle are called and lossless discrete-time PH systems, which is based on variational integrators. They “work very well for both conservative and dissipative or forced mechanical systems” 1. symplectic collocation methods for the geometric numerical servative and dissipative or systems” 1 and lossless ofdiscrete-time PH systems, which with isnumerical based on 1 .. integration Gaussvariational integrators. “workmechanical very well for both consymplectic methods for collocation the geometric servative and dissipativeThey or forced forced mechanical systems” integration collocation of ODEs. ODEs. Continuous Continuous collocation with Gauss1 The port-Hamiltonian (PH) approach, see Duindam et al. symplectic collocation methods for the geometric numerical integration of ODEs. Continuous collocation with Gaussnodes leads aa family discrete-time represenservative and dissipative or forced mechanical systems” The port-Hamiltonian (PH) approach, see Duindam et al.. Legendre nodes leads to to family of ofcollocation discrete-time represenThe port-Hamiltonian (PH)for approach, see simulation Duindam etand al. Legendre (2009), is very appealing modeling, integration of ODEs. Continuous with Gausstations/numerical integrators for systems. The correLegendre nodes leads to a family ofPH discrete-time represen(2009), is very appealing for modeling, simulation and tations/numerical integrators for PH systems. The correThe port-Hamiltonian (PH)for approach, seePH Duindam etand al. Legendre (2009), is complex very appealing modeling, simulation control of multiphysics systems. systems genodes leads to a family ofPH discrete-time represensponding quadrature formulas and backward error analysis tations/numerical integrators for systems. The correcontrol of complex multiphysics systems. PH systems gesponding quadrature formulas and backward error analysis (2009), isthe very appealing for modeling, simulation and neralize Hamiltonian system representation by the control of complex multiphysics systems. PH systems ge- tations/numerical integrators for PH systems. The correallow for quantitative statements on the approximation sponding quadrature formulas and backward error analysis neralize the Hamiltonian system representation by the for quadrature quantitative statements on the approximation control ofthe complex multiphysics systems. PH systems ge- allow additional definition of ports, which are pairs of power neralize Hamiltonian system representation by the sponding formulas and backward error analysis allow for quantitative statements on the approximation error both of and energy. additional definition of ports, which are pairs of by power both of the the solution solution and the the supplied supplied energy. neralize the Hamiltonian system representation the error additional definition of ports, which are pairs of system power variables that characterize energy exchange in the allow for quantitative statements on the approximation error both of the solution and the supplied energy. variables that characterize energy exchange in the system The paper is organized as follows. In Section II, additional definition of The ports, which are pairs of system power and over its boundary. numerical integration of PH variables that characterize energy exchange in the error both of the solution the supplied energy. The paper is organized as and follows. In Section II, we we introintroand over its boundary. The numerical integration of PH The paper is organized asoffollows. In Section II, we introduce the considered class finite-dimensional PH systems variables that characterize energy exchange in the Besides system systems has to account for this energy exchange. and over its boundary. The numerical integration of PH duce the considered class of finite-dimensional PH systems systems has to account for this energy exchange. Besides The paper is organized as follows. In Section II, we introand the underlying Dirac structure. Section III contains as and over its boundary. The numerical integration of PH duce the considered class of finite-dimensional PH systems the error of numerical itself, the of systems to account for solution this energy exchange. the underlying Dirac structure. Section IIIPH contains as the errorhas of the the numerical solution itself, the error errorBesides of the the and duce the considered class of finite-dimensional systems main results the definitions of discrete-time Dirac structusystems has to account for this energy exchange. Besides and the underlying Dirac structure. Section III contains as the error of the numerical solution itself, the error of the energy transmitted over the discrete ports is of fundamenmain results the definitions of discrete-time Dirac structuenergy transmitted over thesolution discreteitself, ports the is oferror fundamenthe underlying Dirac structure. Section III contains as res and PH systems based on the collocation method. We the error of inthe numerical of the and main results the definitions of discrete-time Dirac structutal interest the simulation of PH systems. energy transmitted over the discrete ports is of fundamenres and PH systems based on the collocationDirac method. We tal interest in the simulation of PH systems. main results the definitions of discrete-time structures and PH systems based on the collocation method. We on the Gauss-Legendre collocation. Section energy transmitted over the discrete ports is of fundamen- focus tal interest in the simulation of PH systems. focus onPH the case case of of based Gauss-Legendre collocation. Section Most existing works on formulation and on the collocation method. We IV contains aacase numerical example to illustrate the statetal interest in the simulation of PH systems. focus on thesystems of Gauss-Legendre collocation. Section Most existing works on the the discrete-time discrete-time formulation of of res IV contains numerical example to illustrate the statePH systems make use of a discrete gradient, defined from Most existing works on the discrete-time formulation of focus on the case of Gauss-Legendre collocation. Section ments on the errors of the numerical solution and the IV contains a numerical example to illustrate the statePH systems make use of a discrete gradient, defined from on the errors of the numerical solutionthe and the Most existing works on the discrete-time formulation of ments a finite-differences point (G¨ o ren-S¨ u mer and Yal¸ ccιn PH systems make use of of a view discrete gradient, defined from IV contains a numerical example to Section illustrate stateenergy increment. In the concluding V, we point ments on the errors of the numerical solution and the a finite-differences point of view (G¨ o ren-S¨ u mer and Yal¸ ιn energy increment. In the concluding Section V, we point PH systems make use of a discrete gradient, defined from (2008), Aoues et al. (2013), Falaize H´ lie (2016)). Taa finite-differences point of view (G¨ oand ren-S¨ ueemer and Yal¸ cιn energy ments on the errors of the numerical solution andpoint the increment. In the concluding Section V, we out perspectives to extend and generalize our results. (2008), Aoues et al. (2013), Falaize and H´ lie (2016)). Taperspectives toIn extend and generalize our V, results. a finite-differences point of view (G¨ oand ren-S¨ ueof mer and Yal¸ cιn out (2008), Aoues et al. (2013), Falaize H´ liePH (2016)). Talasila et al. (2006) give a generic definition dynamics energy increment. the concluding Section we point out perspectives to extend and generalize our results. lasila et Aoues al. (2006) give a generic definition dynamics (2008), et al. (2013), Falaize and H´eof liePH (2016)). Ta- out perspectives to extend and generalize our results. lasila et al. (2006) give a generic definition of PH dynamics 2. 2. FINITE-DIMENSIONAL FINITE-DIMENSIONAL PH PH SYSTEMS SYSTEMS lasila et al. (2006) give a generic definition of PH dynamics  P. Kotyczka received financial support as a part-time post-doctoral 2. FINITE-DIMENSIONAL PH SYSTEMS  P. Kotyczka received financial support as a part-time post-doctoral  We the PH researcher (03/17–08/17) from the DFG-ANR funded project INFI2. FINITE-DIMENSIONAL PH SYSTEMS P. Kotyczka received financial support as a part-time post-doctoral We consider consider the class class of of lossless lossless finite-dimensional finite-dimensional PH researcher (03/17–08/17) from the DFG-ANR funded project INFI systems (see e. g. the books Duindam et al. (2009) or We consider the class of lossless finite-dimensional PH DHEM (noo(03/17–08/17) ANR-16-CE92-0028) and byas a apart-time fellowresearcher from the DFG-ANR fundedvisiting project INFIP. Kotyczka received financial support part-time post-doctoral systems (see e. g. the books Duindam et al. (2009) or DHEM (noo ANR-16-CE92-0028) and by a part-time visiting fellowWe consider the class of an lossless finite-dimensional PH ship of Grenoble INP in summer term 2017. The work makes also researcher (03/17–08/17) from theand DFG-ANR funded project INFIvan der Schaft (2017) for overview) systems (see e. g. the books Duindam et al. (2009) or DHEM (n ANR-16-CE92-0028) by a part-time visiting fellowvan der Schaft (2017) for an overview) ship of Grenoble INP in summer term 2017. The work makes also systems (see e. g. the books Duindam et al. (2009) or part theo ANR-16-CE92-0028) project KOin4750/1-1, funded the Research ship of Grenoble INP summerand term TheGerman work makes also DHEM (n by2017. a by part-time visiting fellowvan der Schaft (2017) for an overview) part of the project KO 4750/1-1, funded by the German Research ˙˙ x(t) = J (x(t))∇H(x(t)) + G(x(t))u(t) x(t) = J (x(t))∇H(x(t)) + G(x(t))u(t) Foundation (DFG). ship of Grenoble INP in4750/1-1, summer term 2017. TheGerman work makes also van der Schaft (2017) for an overview) part the project KO funded by the Research (1) Foundation (DFG). T ˙ x(t) = J (x(t))∇H(x(t)) + G(x(t))u(t) 1 See Lew et (1) al. (2004), Section 2. funded by the German Research Foundation (DFG). part the et project KO 4750/1-1, T (x(t))∇H(x(t)) y(t) = G 1 SeeofLew al. (2004), Section 2. (1) y(t) =J G(x(t))∇H(x(t)) ˙ x(t) + G(x(t))u(t) T (x(t))∇H(x(t)) 1 T 1 Foundation See Lew et(DFG). al. (2004), Section 2. y(t) = GT (x(t))∇H(x(t)) (1) 1 See Lew et al. (2004), Section 2. y(t) = G (x(t))∇H(x(t)) 2405-8963 © © 2018 2018, IFAC IFAC (International Federation of Automatic Control) Copyright 125 Hosting by Elsevier Ltd. All rights reserved. Copyright © under 2018 IFAC 125 Control. Peer review responsibility of International Federation of Automatic Copyright © 2018 IFAC 125 10.1016/j.ifacol.2018.06.035 Copyright © 2018 IFAC 125

IFAC LHMNC 2018 126 Valparaíso, Chile, May 1-4, 2018

Paul Kotyczka et al. / IFAC PapersOnLine 51-3 (2018) 125–130

with state vector x ∈ Rn , in- and collocated output u, y ∈ Rm . The Hamiltonian H : Rn → R is bounded from below with a strict minimum in x∗ , which is the equilibrium state for u ≡ 0. By skew-symmetry of the interconnection matrix J = −J T and the definition of the collocated output, the differential energy balance ˙ H(x(t)) = y T (t)u(t) (2) holds, or in integral form,  H(x(t2 ))−H(x(t1 )) =

t2

T

y (s)u(s) ds,

t1

∀t1 ≤ t2 . (3)

This energy balance is a structural or geometric property, i. e. it holds independently of H(x). To verify this, define the vectors of flows and efforts ˙ f (t) := −x(t), e(t) := ∇H(x(t)), (4) T ˙ T ˙ which due to H = (∇H) x = −e f are called powerconjugated. The differential energy balance (2) can be written as the power balance equation on the bond space F × E of power-conjugated variables, with F = Rn × Rm  (f , y) and E = Rn × Rm  (e, u): eT (t)f (t) + y T (t)u(t) = 0.

(5)

By −f (t) = J (x(t))e(t) + G(x(t))u(t)

(6) y(t) = GT (x(t))e(t), the bond variables are constrained to a subspace, on which in particular (5) holds. This subspace is called a Dirac structure. For more details on Dirac structures and PH systems, see e. g. van der Schaft (2017), Chapter 6. 3. DISCRETE-TIME PH SYSTEMS BASED ON COLLOCATION In this main section of the paper, we define a class of discrete time port-Hamiltonian systems that arise from the definition of a discrete-time Dirac structure. The latter is obtained by applying the collocation method to the class of PH systems (1) and by defining in an appropriate manner discrete flow and effort vectors for every sampling interval. It turns out that the symplectic Gauss-Legendre collocation method, if applied to an open PH system, is structurepreserving in the sense of the Dirac structure. Moreover, the associated quadrature formula allows for a quantification of the error in the discrete energy balance. 3.1 Collocation method

k−1

ts+1 tk−1 s tk0

tk+1 tk+1 0 1

Ik tk1

tk2

tks−1 tks tks+1

...

t

h

Fig. 1. Illustration of the k-th sampling interval Ik with interior collocation points tki = tk0 + ci h, i = 1, . . . , s. x ˜˙ (tki ) = −f ki , −f ki

i = 1, . . . s,

= (J (x)∇H(x))|x=˜x(tk ) + G(˜ x(tki ))u(tki ).

(7)

i

Notation: Arguments in latin letters (t or s under the integral), refer to time functions evaluated on I k . Greek letters (τ or σ) refer to the same function, mapped to the normalized interval [0, 1]. 3.2 Approximation of flow and state variables Based on the evaluations of the state differential equation in the collocation points ci , denoted −f ki ∈ Rn and called the discrete flow coordinates (with negative sign), the interpolation formula s  x ˜˙ (tk0 + τ h) =: −f˜ (tk0 + τ h) = − f ki i (τ ), (8) i=1

with i the i-th Lagrange interpolation polynomial s  τ − cj i (τ ) = , τ ∈ [0, 1], c − cj j=1 i

(9)

j=i

˙ gives a polynomial approximation of x(t) on I k . The flow coordinates are merged in the discrete-time flow vector T  f k := (f k1 )T . . . (f ks )T ∈ Rsn , (10) ˜ (tk + τ h), based on which the the numerical solution x τ ∈ [0, 1] is obtained by integration:  τ ˜ (tk0 ) + h ˜ (tk0 + τ h) = x x ˜˙ (tk0 + σh) dσ x 0

˜ (tk0 ) − h =x xki

˜ (tki ), x

s  

f kj

j=1



τ

0

 j (σ) dσ .

(11)

The values := i = 1, . . . , s + 1 of the numerical solution inside and at the end of the sampling interval then are computed as s  xki = xk0 − h aij f kj , i = 1, . . . , s, (12) j=1

We consider equidistant sampling intervals I k = [tk0 , tks+1 ] = [(k − 1)h, kh], k ∈ N for the time t with tks+1 = tk0 + h, see Fig. 1. With t = ((k − 1) + τ )h, the sampling intervals are parametrized by the normalized time τ ∈ [0, 1]. The polynomial approximations of the system variables will be denoted with a tilde. As described in Section II.1.2 of Hairer et al. (2006), the numerical approximation of ˜ (t) ∈ Rn the solution x(t) of (1) is given by the vector x of collocation polynomials of degree s. Assume first the ˜ (tk0 ) = x(tk0 ) to be known. The initial value xk0 := x ˜ (t) is then the vector of continuous numerical solution x polynomials whose time derivative matches the right hand side of the differential equation in the s collocation points tki := tk0 + ci h with 0 ≤ ci ≤ 1: 126

xks+1 = xk0 − h with 2 aij =



s 

0

(13)

j=1

ci

j (σ) dσ,

bj f kj ,

bj =



1

j (σ) dσ,

i, j = 1, . . . , s.

0

(14) In continuous collocation methods, the numerical solution at the start tk+1 = tks+1 of the subsequent interval is 0 k+1 initialized by x0 = xks+1 . 2 Note that these values are, together with c , the coefficients of i the Butcher tableau for the Runge-Kutta (RK) interpretation of the collocation method, see Hairer et al. (2006), Theorem II.1.4.

IFAC LHMNC 2018 Valparaíso, Chile, May 1-4, 2018

Paul Kotyczka et al. / IFAC PapersOnLine 51-3 (2018) 125–130

3.3 Effort approximation and discrete structure equation The definition of the discrete flow coordinates f k1 , . . . f ks in (7) requires to evaluate the effort vector ∇H(x(t)), the input u(t) and the interconnection and input matrices J (x(t)), G(x(t)) in the flow collocation points ci , (15) eki := ∇H(x)|xk , uki := u(tk0 + ci h), i

and

J ki := J (xki ), Gki := G(xki ) (16) for i = 1, . . . , s. The discrete-time counterpart of the first equation in (6) is J ki

−f ki = J ki eki + Gki uki ,

−(J ki )T .

(17)

= With ∈ R and ∈ R the with discrete effort and discrete input coordinates according to (15), the polynomial approximations of the effort and the input vector are s s   ˜ (tk0 +τ h) = ˜(tk0 +τ h) = eki i (τ ), u uki i (τ ). (18) e eki

n

i=1

uki

m

i=1

In accordance with the approximate flows, we define the discrete-time effort vector  T (19) ek = (ek1 )T . . . (eks )T ∈ Rsn and the discrete-time input vector  T (20) uk = (uk1 )T . . . (uks )T ∈ Rsm . Defining the block-diagonal matrices J k = −(J k )T = blockdiag(J k1 , . . . , J ks ), (21) Gk = blockdiag(Gk1 , . . . , Gks ), Equation (17), i. e. the structure equation on the sampling interval I k can be rewritten as (22) −f k = J k ek + Gk uk . Remark 1. Defining the discrete effort and input coordinates based on different collocation points d1 , . . . , dr is also conceivable. In this case, the terms of the right hand side of (17) must be replaced by interpolations between the effort collocation points according to (18). In this paper, we restrict ourselves to identical collocation points. 3.4 Discrete-time energy balance We seek to express a discrete-time counterpart of the integral energy balance equation (3) on the time interval I k . To this end, we integrate  ∂H(x)  k ˙ x ˜˙ (tk0 + τ h) H(˜ x(t0 + τ h)) = ∂x  k x=˜ x(t0 +τ h) T k −˜ e (t0 + τ h)f˜ (tk0

+ τ h) = over the normalized time interval [0, 1], ˜ k := H(˜ ∆H x(tks+1 )) − H(˜ x(tk0 ))  1 ˜T (tk0 + σh)f˜ (tk0 + σh) dσ. e = −h

(23)

(24)

0

Substituting the polynomial approximations of efforts (18) and flows (8), we obtain the bilinear form ˜ k = −h(ek )T M f k ∆H (25)

with the symmetric matrix M = M T ∈ Rsn×sn that consists of the n × n submatrices  1 [M ]ij = mij I n , mij = mji = i (σ)j (σ) dσ. (26) 0

127

127

Remark 2. We can understand the term M ek as a generalization of the discrete gradient and −hf k as a vector expressing a generalized increment of the numerical solution in the integration step. 3.5 Discrete-time Dirac structure Substituting the relation (22) in the right hand side of the discrete energy balance (25), we obtain −h(ek )T M f k = h(ek )T M J k ek + h(ek )T M Gk uk . (27) At this stage, we want to recover a discrete-time energy balance equation with the structure of (5). To this end, the first term on the right hand side must vanish: ! h(ek )T M J k ek = 0, or written out (recall mij = mji ),   k   e1 m11 I n . . . m1s I n J k1  k T   ..   ..  . .. . h (e1 ) . . . (eks )T  ...  .  . . .  ms1 I n . . . mss I n J ks eks s  s  ! =h (eki )T mij J kj ekj = 0. (28) i=1 j=1 k J j , we have

By skew-symmetry of (ekj )T mjj J kj ekj = 0 for all j = 1, . . . , s. The remaining elements of the sum cancel to zero, if (eki )T mij J kj ekj = −(ekj )T mji J ki eki holds for all i = j. With (eki )T mij J kj ekj = −((eki )T mij (J kj )T ekj )T = −(ekj )T mji J kj eki , the requirement translates to

(29) eTj mji J ki ei = eTj mji J kj ei , which is true if either of the following conditions holds. (C1) (C2)

mij = 0 for all i = j, J ki = J kj = const. for all i, j = 1, . . . , s.

While (C1) is a condition on the discretization method, the constant interconnection structure according to (C2) is a system property 3 . In both cases (C1) and (C2), the discrete energy balance (27) boils down to h(ek )T M f k + h(ek )T M Gk uk = 0. The definition of the discrete output vector y k := h(Gk )T M ek yields (using M = M T )

(30) (31)

(32) h(M ek )T f k + (y k )T uk = 0 as the discrete-time counterpart of the (structural) power balance equation (5). Unlike the interpretation of (instantaneous) power flows in (5), the terms of (32) represent the (accumulated) energy flows to the storage elements and over the system boundary during the sampling interval I k . We are ready to define the discrete-time Dirac structure which is associated to the sketched family of time discretization schemes. Proposition 3. (Discrete-time Dirac structure). Given the s collocation points 0 ≤ ci ≤ 1, i = 1, . . . , s. The system of equations (22), (31), i. e. −f k = (J k M −1 )M ek + Gk uk y k = h(Gk )T M ek

(33)

3 The Darboux-Lie theorem states that every Poisson system can be transformed locally to a canonical Hamiltonian system with a constant skew-symmetric matrix (Hairer et al. (2006), Section VII.3).

IFAC LHMNC 2018 128 Valparaíso, Chile, May 1-4, 2018

Paul Kotyczka et al. / IFAC PapersOnLine 51-3 (2018) 125–130

with discrete flow, effort and input vectors f k , ek , uk according to (10), (19), (20), the block matrices J k = −(J k )T , Gk according to (21), and the symmetric block matrix M = M T according to (26), represents a discretetime Dirac structure on the time interval I k = [(k−1)h, kh], if either condition (C1) or (C2) is satisfied. Proof. Rewrite (33) in kernel representation    k M ek hf +E =0 (34) F uk yk with  k −1  G J M . (35) F = I s(n+m) , E = h −GT 0 According to van der Schaft (2017), Prop. 6.6.6, the subspace of discrete bond variables (hf k , M ek ) and (y k , uk ) on which (34) holds, is a Dirac structure, if (i) EF T + F E T = 0 and (ii) the matrix [F E] has full row rank s(n + m). The latter is obvious from the identity matrix F in (35). Straightforward computations show that (i) is true if and only if J k M or, equivalently, M J k is skew-symmetric. As shown above, the latter holds if the continuous-time system has a constant interconnection matrix (C2). Otherwise, the choice of collocation points must guarantee (C1).  3.6 Discrete-time port-Hamiltonian system The discrete-time Dirac structure will be complemented by dynamics and constitutive equations. Definition 4. (Discrete-time PH system). Equations (33), together with the discrete dynamics (36a) xk0 = xk−1 s+1 , s  xki = xk0 − h aij f kj , i = 1, . . . , s, (36b) xks+1 = xk0 − h

j=1 s 

bj f kj ,

(36c)

j=1

with the Runge-Kutta coefficients aij and bj according to (14) and the discrete constitutive equations i = 1, . . . , s, (37) eki = ∇H(x)|x=xk , i

define a discrete-time PH system. Using (25) and (32), the numerical increment of the energy H(x) on the sampling interval I k = [(k − 1)h, kh] is given by ˜ k = (y k )T uk . (38) ∆H

Remark 5. Equations (36b) and (36c) can be written in more compact form using the Kronecker product (39a) xki = xk0 − h(aTi ⊗ I n )f k , i = 1, . . . , s, xks+1 = xk0 − h(bT ⊗ I n )f k ,

(39b)

with aTi = [ai1 . . . ais ], bT = [b1 . . . bs ]. Remark 6. The structure of the discrete-time PH models is independent of the sampling time h, which, however, determines the approximation quality.

numerical integration to PH systems in a way that it coincides with a symplectic scheme for the autonomous case (u ≡ 0). Moreover, we want to quantify and minimize the approximation error for the energy increment in each sampling interval. As we do not want to restrict the integration scheme to constant interconnection matrices, we seek to satisfy condition (C1) by the choice of collocation points. To have  1 i (σ)j (σ) dσ = 0 ∀i = j, (40) mij = 0

the interpolation polynomials {1 (τ ), . . . , s (τ )} must form a system of orthogonal functions. This is the case, if we take the collocation points c1 , . . . , cs as the zeros of the shifted Legendre polynomials 4 ds (τ s (τ − 1)s ) (41) dτ s in normalized time. With the resulting interpolation polynomials, the equalities  1  1 1 bi = i (σ) dσ = 2i (σ) dσ = mii (42) h 0 0 hold, where mii are the factors in the block diagonal matrix M and the bi are used as weights in the GaussLegendre quadrature formula  1 s  f (σ) dσ ≈ bi f (ci ) (43) 0

i=1

for a function f (τ ) on [0, 1] with nodes ci , i = 1, . . . , s.

This quadrature formula, see Table A.1 with the values of ci and bi for s = 1, 2, 3, is exact for polynomials f (τ ) up to order 2s − 1. With 2s the order of the quadrature formula, the approximation error of the integral is of order 2s + 1. Proposition 7. (Approximation error). Consider a discretetime PH system according to Definition 4 with quadratic energy H(x) = 21 xT Qx, Q = QT > 0 and Gauss-Legendre collocation points. The local error of the numerical solu˜ (t) is of order 2s+1. The relative error of the supplied tion x ˜ k in each sampling interval is of order 2s. energy ∆H Proof. The RK method based on Gauss-Legendre collocation has the same order 2s as the quadrature formula, see Hairer et al. (2006), Section II.1.3. Hence, the local ˜ (tk0 ) = x(tk0 ) approximation error is of order 2s+1 and for x ˜ x(tk0 + h) − x(tk0 + h) ≤ ck h2s+1

(44)

˜ x(tk0 + h) − x(tk0 + h)Q ≤ ckQ h2s+1

(45)

k

with some c > 0. Equivalently,

1 T 2 x Qx.

By triangle with the energy norm xQ := inequality and mean value theorem, we get for ∆H k = 0     ˜ k − ∆H k  ˜ ∆H x(tk0 + h)Q − x(tk0 + h)Q      = ∆H k  ˙ x(tk + τ ∗ h)) hH(˜ 0 ckQ ckQ h2s+1 =  h2s (46) ≤  ˙ x(tk + τ ∗ h)) ˙ x(tk + τ ∗ h)) H(˜ hH(˜ 0

3.7 Gauss-Legendre collocation In the previous definition, the collocation points c1 , . . . , cs are degrees of freedom to parametrize the numerical integration scheme. We want to extend the notion of geometric 128

0

with τ ∗ ∈ [0, 1], which represents the relative energy error on I k .  4

See Hairer et al. (2006), Section II.1.3.

IFAC LHMNC 2018 Valparaíso, Chile, May 1-4, 2018

Paul Kotyczka et al. / IFAC PapersOnLine 51-3 (2018) 125–130

129

Remark 8. The order of the local approximation error can be proven by backward error analysis 5 . See Appendix B for its application to the example of Section 4 with s = 1. 4. EXAMPLE To illustrate our definition of discrete-time PH systems based on the collocation method, and to verify the statement on the accuracy of the supplied energy, we consider as the simplest PH system, a linear oscillator with state x(t) ∈ R2 , excited by the input u(t) ∈ R, ˙ x(t) = J ∇H(x(t)) + gu(t) (47a) y(t) = g T ∇H(x(t)).

(47b)

1 T 2 x x,

The quadratic Hamiltonian is H(x) = and input matrices are     0 1 0 J= , g= . −1 0 1

the structure (48)

The system is integrated numerically as described in Definition 4, which defines a discrete-time PH system. Table A.2 shows the coefficients of the associated RK method based on Gauss-Legendre collocation. 4.1 Local error

For illustration, we consider s = 1. The resulting RK method with the coefficients a11 = c1 = 12 and b1 = 1 is the midpoint rule, and can be written as h h (49a) k1 = J (xk0 + k1 ) + gu(tk0 + ), 2 2 = xk0 + hk1 . (49b) xk+1 0 By linearity of the system, it is possible to explicitly determine the slope k1 as h h h k1 = (I − J )−1 J xk0 + (I − J )−1 gu(tk0 + ). (50) 2 2 2 For h < 2, the eigenvalues of h2 J lie inside the unit disc. The inverse can then be written as a convergent infinite geometric series of matrices and the flow of the numerical method, starting with initial value xk0 at time tk0 is ∞ ∞  h i k  h i h  k k ( J ) J x0 + ( J ) gu(tk0 + ) . ψ h (x0 ) = x0 + h 2 2 2 i=0 i=0 (51) In Appendix B, we verify for s = 1, using backward error analysis, that the error between the numerical and the exact solution for s = 1 is of order O(h2s+1 ) = O(h3 ). 4.2 Global energy error We illustrate the approximation order for the total supplied energy in a time interval [0, T ] by the following numerical experiment. Consider system (47), (48) with initial value x0 = [0 − 1]T and the squared sine input pulse (0 ≤ tA ≤ tB ≤ T )  sin2 ( t − tA ), tA ≤ t ≤ tB , u(t) = (52) tB − t A 0, else. 5

See e. g. Hairer et al. (2006), Chapter IX, or Leimkuhler and Reich (2004), Chapter 5.

129

Fig. 2. Relative energy error during the input pulse. The exact solution with tA = 0.5π and tB = 1.5π can be computed as     1 sin2 (t) − 5 sin(2t) + 1 sin(t) xI = − , , xII = cos(t) sin(2t) − 5 cos(t) 3   7 sin(t) xIII = − , (53) 3 cos(t) where the indices I, II and III denote the time intervals before, during and after the input pulse. The total supplied energy ∆Htot during phase II is 1 7 20 ∆Htot = (( )2 − 1) = . (54) 2 3 9 The solution x(t) in the phase plane is lifted to a circle with increased radius. The results from backward error analysis can be exploited to quantify the error in the numerical approximation of the energy, which is supplied to the system over its control ¯ (see port (u, y). The change of the modified Hamiltonian H Appendix B) along the solutions of the modified system is ˜ tot along the numerical an approximation of the change ∆H solution. Comparing this increment with its exact value ∆Htot yields an expression to quantify the order of the ¯ = global energy error. With the modified Hamiltonian H h2 h2 ¯ = (1 − 12 )g, (1 − 12 )H and the modified input vector g we can compute  T  T ¯ ∂H ¯ +g ¯˙ ds = ¯ tot = ¯ u) ds (J ∇H H ∆H ∂x 0 0 (55) h4 h2 ), = ∆Htot − ∆Htot ( − 6 144 which is an indication  for the order 2s of the global energy  ∆H˜ tot −∆Htot  error  . Note the slopes of the curves for s = 1 ∆Htot and s = 2 in Fig. 2, where the error is plotted, based on the results of numerical integration, over the step size h. 5. CONCLUSION We presented a new definition of lossless discrete-time PH systems, which is based on geometric numerical integration using continuous collocation methods. We expressed discrete-time relations between power variables and gave conditions under which these relations characterize a discrete-time Dirac structure. Adding discrete constitutive equations and discrete-time update equations (dynamics)

IFAC LHMNC 2018 130 Valparaíso, Chile, May 1-4, 2018

Paul Kotyczka et al. / IFAC PapersOnLine 51-3 (2018) 125–130

Table A.1. Nodes and weights for GaussLegendre quadrature on the interval [0, 1]

results in a very general definition of discrete-time PH systems. Choosing Gauss-Legendre collocation points, the tight relation between integration method and quadrature formula allows for quantifying the order of the error in the discrete-time energy balance.

s 1

Ongoing work is on the extension of the results towards partitioned and discontinuous methods in order to provide a general framework for the structure-preserving numerical integration of PH systems. In such a framework, integrating dissipation and considering PH systems described by DAEs are important aspects. Combined with structurepreserving spatial discretization, see e. g. Kotyczka et al. (2018), the presented approach contributes to the full discretization of distributed-parameter PH systems.

2 3

c1 1 2√

1 3 − 2 √6 1 15 − 2 10

1 2

Appendix A. GAUSS-LEGENDRE COLLOCATION Tables A.1 and A.2 show nodes and weights related to Gauss-Legendre collocation. Appendix B. BACKWARD ERROR ANALYSIS A modified differential equation is constructed, whose flow is compared with the numerical solution based on a Taylor series expansion in the time step h. We consider the example system (47) and the case s = 1. We start with the free system, i. e. u ≡ 0. Denote f := J ∇H and let f¯ := f + f 1 h + f 2 h2 , (B.1) 130

c3

s

b1

b2

b3

− √

−

1 2

− 1 2 4 9

−

− √ 15 1 + 2 10

1 1 2 5 18

3 1 + 2 6 1 2

3

− 5 18

Table A.2. Butcher tables for Gauss-Legendre collocation as Runge-Kutta method (s = 1, 2)

REFERENCES Aoues, S., Eberard, D., and Marquis-Favre, W. (2013). Canonical interconnection of discrete linear portHamiltonian systems. In 52nd IEEE Conference on Decision and Control, 3166–3171. Brewer, J. (1978). Kronecker products and matrix calculus in system theory. IEEE Transactions on Circuits and Systems, 25(9), 772–781. Duindam, V., Macchelli, A., Stramigioli, S., and Bruyninckx, H. (2009). Modeling and Control of Complex Physical Systems: The Port-Hamiltonian Approach. Springer Science & Business Media. Falaize, A. and H´elie, T. (2016). Passive guaranteed simulation of analog audio circuits: A port-Hamiltonian approach. Applied Sciences, 6(10), 273. G¨ oren-S¨ umer, L. and Yal¸cιn, Y. (2008). Gradient based discrete-time modeling and control of Hamiltonian systems. IFAC Proceedings Volumes, 41(2), 212–217. Hairer, E., Lubich, C., and Wanner, G. (2006). Geometric numerical integration: structure-preserving algorithms for ordinary differential equations, volume 31. Springer Science & Business Media. Kotyczka, P., Maschke, B., and Lef`evre, L. (2018). Weak form of Stokes-Dirac structures and geometric discretization of port-Hamiltonian systems. Journal of Computational Physics, 361, 442–476. Leimkuhler, B. and Reich, S. (2004). Simulating Hamiltonian dynamics, volume 14. Cambridge University Press. Lew, A., Marsden, J.E., Ortiz, M., and West, M. (2004). An overview of variational integrators. Talasila, V., Clemente-Gallardo, J., and van der Schaft, A. (2006). Discrete port-Hamiltonian systems. Systems & Control Letters, 55(6), 478 – 486. van der Schaft, A. (2017). L2-Gain and Passivity Techniques in Nonlinear Control. Springer, 3rd edition.

c2

1 2 1

√ 3 1 − 2 √6 3 1 + 2 6

1 4√

1 3 + 4 6 1 2

√ 3 1 − 4 6 1 4 1 2

with f 1 and f 2 to be determined. The flow of this modified vector field at time h, starting from xk0 is written as h2 ¯˙ k h3 ¨¯ k ϕh,f¯(xk0 ) = xk0 +hf¯ (xk0 )+ f (x0 )+ f (x0 )+. . . (B.2) 2 6 with (using Vetter’s notation for the Jacobian 6 )  ¯  ¯ ∂f ¯ ¯ ¨¯ = ∂ ¯˙ = ∂ f f¯ , f f f , etc. (B.3) f T T ∂x ∂x ∂xT By matching the coefficient functions in front of the powers of h in (B.2) and (51), we obtain 1 f 1 = 0, (B.4) f 2 = − J ∇H(x), 12 such that the modified vector field can be expressed in ¯ terms of a modified Hamiltonian H, 2 ¯ ¯ = (1 − h )H. (B.5) f¯ = J ∇H, H 12 By this choice, also the 4th order coefficients of (B.2) and (51) coincide, which means that the error between both flows at time h is of order O(h5 ). For u = 0, we modify the input vector field accordingly: h2 ¯ = (1 − )g. (B.6) g 12 For the flow ϕh,f¯+¯gu , we first add the terms

h2 k h3 k ¯ u˙ 0 + g ¯u g ¨0 + . . . (B.7) 2 6 to (B.2). Second, u and its time derivatives of order ν = 0, 1, 2, . . . must be expressed in terms of the values in the center of the sampling interval (subscript “1”), (ν)k (ν) h (ν+1)k u 0 = u k1 − u 1 + O(h2 ). (B.8) 2 The difference of the coefficients in the Taylor series of the numerical solution and the modified flow yields the error h3 ψ h (xk0 ) − ϕh,f¯+¯gu (xk0 ) = − (J g u˙ k1 + g u ¨k1 ) + O(h4 ). 12 (B.9) The error between the flows of the original vector field ¯ ¯ u is also of (47a) and the modified version J ∇H(x) +g order O(h3 ), which concludes the demonstration of second order accuracy of the considered integration scheme. h¯ g uk0 +

6

See Brewer (1978).