Discrete-Time Models for Implicit Port-Hamiltonian Systems
Fernando Castaños, Hannah Michalska, Dmitry Gromov, Vincent Hayward
aa r X i v : . [ c s . S Y ] J a n DISCRETE-TIME MODELS FORIMPLICIT PORT-HAMILTONIAN SYSTEMS
FERNANDO CASTA ˜NOS ∗ , HANNAH MICHALSKA † , DMITRY GROMOV ‡ , AND
VINCENTHAYWARD § Abstract.
Implicit representations of finite-dimensional port-Hamiltonian systems are studiedfrom the perspective of their use in numerical simulation and control design. Implicit representationsarise when a system is modeled in Cartesian coordinates and when the system constraints are appliedin the form of additional algebraic equations (the system model is in a differential-algebraic equationform). Such representations lend themselves better to sample-data approximations. Given an implicitrepresentation of a port-Hamiltonian system it is shown how to construct a sampled-data model thatpreserves the port-Hamiltonian (PH) structure under sample and hold.
Key words. port-Hamiltonian systems, nonlinear implicit systems, symplectic integration,sampled-data systems
1. Introduction.
The class of Hamiltonian systems has a prominent role inmany disciplines. It was recently extended in [6] to include open systems, i.e. systemsthat interact with the environment via a set of inputs and outputs (called ports ),giving rise to port-Hamiltonian (PH) systems. Such extended models immediatelyreveal the passive properties of the underlying systems, making them particularlywell suited for designing passivity-based control (PBC) laws. Two types of modelrepresentations of Hamiltonian systems are in widespread use: the explicit repre-sentation stated in the form of an ordinary differential equation (ODE) on an ab-stract manifold [28, 29, 25, 26] and the implicit representation stated in the form of adifferential-algebraic equation (DAE) usually evolving in a Euclidean space [4]. Whileexplicit representations are usually preferred in the context of analytical mechanics;see [1, 8, 19], the implicit DAEs models lend themselves better for numerical compu-tations as they lead to simpler expressions for the Hamiltonian functions. The formalrelations between the two representations and their equivalence can be establishedif the system’s configuration space is regarded as an embedded submanifold of theEuclidean space; see [4] for a full development.Of principal interest here will be the construction of sampled-data (discrete-time)models of PH systems that best approximate their continuous-time counterparts.Sampled-data models are important in digital control implementations and permitfor simpler design of PBC laws directly in discrete time. In this context, the notionof “best approximation” deserves clarification.For linear systems, exact sampled-data models can be constructed by requiringthat the solutions of the sampled-data and continuous-time systems coincide at thediscretization points. Point-wise model matching is usually impossible in the case ofnonlinear systems short of explicit derivation of analytical expressions for their solu-tions. For general dynamic systems, it is thus the choice of an integration methodfor generating an approximate solution (like Euler or Runge-Kutta) whose precision ∗ Departamento de Control Automatico, Cinvestav del IPN, Av. IPN 2508, Col. San PedroZacatenco, C.P. 07360, M´exico D.F., M´exico ( [email protected]) † McGill Centre for Intelligent Machines, 3480 University Street, Montreal Quebec, Canada H3A2A7 ( [email protected]) ‡ Faculty of Applied Mathematics, St.Petersburg State University, University pr., 35, Petrod-vorets, 198 904 St.Petersburg, Russia ( [email protected]) § UPMC Paris 6, UMR7222, Institut de Syst`emes et de Robotique, Paris, France1
Fernando Casta˜nos, Hannah Michalska, Dmitry Gromov and Vincent Hayward relates to the compromise between complexity and order of approximation of thegiven continuous-time system by its sampled-data counterpart. In the special case ofHamiltonian systems, structure preservation is usually the main criterion for choosinga numerical method; see, e.g., [11] for more details on structure-preserving numericalschemes. Structure-preserving methods guarantee accuracy in long-horizon simula-tions.It is known that autonomous Hamiltonian systems conserve two quantities: theHamiltonian function H (i.e., the energy or storage function) and a certain two-form ω , called the symplectic form. Numerical integration algorithms can eitherconserve H [14, 15] or ω [27], but not both. Conservation of ω is usually preferredover conservation of H as the symplectic form unambiguously defines the class ofHamiltonian vector fields; see Theorem 2.3 and Remark 1. Contribution.
Adopting a symplectic-form preserving approach, a sampled-data model for a given continuous-time PH system in developed here. The approachemploys implicit representations of PH systems as the latter lead to simpler expres-sions for the Hamiltonian functions. Specifically, the Hamiltonian functions arisingfrom implicit representations lend themselves well to the application of flow-splittingnumerical integration methods; see [27] for a splitting method that applies to au-tonomous Hamiltonian systems (Hamiltonian systems without ports). Additionally,the discrete-time modeling approach presented here preserves passivity of PH systems(Theorem 4.2). The passive structure is preserved in the sense that the discrete modelis the exact representation of another, possibly non affine, continuous-time PH sys-tem which, up to an approximation error of order two with respect to the samplinginterval, has the same storage function H and the same output function y . Paper structure.
Section 2 introduces the implicit representations of PH sys-tems. This section recalls the results from [4] and serves mainly to introduce thenotation and to state the main assumptions (1 and 2). Also, the definition of PH sys-tems is extended to the case of non affine control systems. Splitting and symplecticintegration methods for autonomous systems are recalled in §
3. Discrete-time mod-els for implicit PH systems based on vector splitting methods are developed in §
2. Implicit representations of port-Hamiltonian systems.
We begin thissection by defining the configuration and phase spaces of Hamiltonian systems inimplicit form. Conservative properties are also discussed. Similarly to the case ofexplicit representations of Hamiltonian systems, both the Hamiltonian function andthe symplectic two-form of the implicit representation are invariant under the flowof the system. Following [6], input and output ports are added to the implicit rep-resentation, insuring passivity of the resulting PH system. Finally, we extend thedefinition of PH to the the case of non affine control systems. This will be neededwhen performing backward error analysis later in § The class of systems considered here is re-stricted to mechanical finite-dimensional systems with scleronomic constraints thatevolve in continuous time. Systems of this type typically consist of M rigid bodiesheld together as one structure by the action of constraint forces. The position of eachof the rigid bodies can be unambiguously described in terms of the position of its cen-ter of mass and its orientation in space. The configuration space of a spatial system as iscrete-Time Models for Implicit port-Hamiltonian Systems
3a whole can thus be viewed as a subset, G , of a Cartesian space of dimension n = 6 M that is expressed in terms of smooth, independent constraint functions g : R n → R k ,with k ≤ n , as the level set G := g − (0) = { r ∈ R n | g ( r ) = 0 } (2.1)where g − (0) is the inverse image of 0 ∈ R k . Functional independence of the con-straints is expressed in terms of the following rank condition . Assumption 1.
The rank of g is equal to k at all points of the set g − (0) . Recalling that the rank of a mapping is the rank of its tangent map (push-forward by g ), then in any coordinates, the condition requires that the Jacobian G := { ∂g i /∂r j } ij has full rank, i.e., rank( G ) = k at every point of g − (0). The value0 ∈ R k is then called a regular value of g , the level set g − (0) is a regular level set of g , and g is said to be a defining map for g − (0); see [16], p. 113-114. It hence followsthat G can be given the differentiable structure of a closed embedded submanifold of R n ; see [16], Corollary 5.24, p. 114. The Constant-Rank Level Set Theorem ([16],Theorem 5.22, p. 113), further specifies its dimension as o = n − k .Because G is an embedded submanifold, it can also be regarded as an abstract o -dimensional manifold with local coordinates q (called generalized coordinates inmechanics). It is then also possible to exhibit an injective inclusion map: ı : G ֒ → R n (embedding), with rank( ı ) = dim G , which is a homeomorphism onto the image ı ( G ) ⊂ R n and such that ı ( q ) = r and g ◦ ı ≡
0. This inclusion serves as a mapbetween local coordinates q and global Cartesian coordinates r . Let T ∗ G be the cotangent bundle of G and let (cid:8) q i , ˆ p i (cid:9) be local coordinates. A system is said to be Hamiltonian if its trajectoriesare integral curves of the Hamiltonian vector field D ˆ H : T ∗ G → T ( T ∗ G ), D ˆ H = ∂ ˆ H∂ ˆ p i ∂∂q i − ∂ ˆ H∂q i ∂∂ ˆ p i , where ˆ H : T ∗ G → R is the sum of the system’s kinetic and potential energies ˆ K : T ∗ G → R and ˆ V : G → R , respectively; see [1] for more details and a coordinate-freedefinition.The implicit model for a Hamiltonian system is defined as follows. Let T ∗ R n be the cotangent bundle of R n and let (cid:8) r i , p i (cid:9) be global coordinates. The implicitHamiltonian vector field takes the form X H,g = D H − λ j ∂g j ∂r i ∂∂p i , g = 0 , (2.2)with D H = ∂H∂p i ∂∂r i − ∂H∂r i ∂∂p i (2.3)the unconstrained part of the Hamiltonian vector field and H : T ∗ R n → R is againthe sum of the kinetic and potential energies, but expressed in global coordinates.That is, K : T ∗ R n → R and V : R n → R ; see [11, 27, 3] for details on the derivationof this equation. The vector field unravels as the DAE system˙ r = ∇ p H ( r, p ) , ˙ p = −∇ r H ( r, p ) − G ( r ) ⊤ λ (2.4a)0 = g ( r ) . (2.4b) Fernando Casta˜nos, Hannah Michalska, Dmitry Gromov and Vincent Hayward
By applying X H,g to both sides of the constraint equations g j = 0, one obtainsthe so-called hidden constraints f j := X H,g ( g j ) = ∂H∂p i ∂g j ∂r i = 0 . (2.5)Thus, the system evolves on the closed submanifold L G = { ( r, p ) ∈ T ∗ R n | g ( r ) = 0 , f ( r, p ) = 0 } (2.6)and we have X H,g : L G → T ( L G ). A more rigorous construction of the phase space L G is given in [4], where it is regarded as an embedding of T ∗ G in T ∗ R n . The mapˆ p p as well as the formal relation between ˆ H and H can also be found in [4].The Lagrange multipliers λ j are defined implicitly by (2.2). Precisely, applicationof X H,g to the hidden constraints makes λ j appear: X H,g ( f l ) = D H ( f l ) − λ j ∂g j ∂r i ∂f l ∂p i = 0 . (2.7)Thus, if the matrix (cid:26) ∂g j ∂r i ∂f l ∂p i (cid:27) jl = (cid:26) ∂g j ∂r i ∂ H∂p i ∂p m ∂g l ∂r m (cid:27) jl (2.8)is non-singular on L G , then there are unique λ j satisfying (2.7) and forcing the integralcurves to stay on L G . Assumption 2.
The Hessian matrix (cid:8) ∂ H ( r, p ) /∂p i ∂p j (cid:9) ij is positive definite forall ( r, p ) ∈ T ∗ R n so H ( r, p ) is convex in p . This assumption, satisfied by most mechanical systems, together with Assump-tion 1 ensures that the matrix in (2.8) is invertible and the system is well defined. Inmechanical systems, λ is the covector of constraint forces that insure satisfaction ofthe constraints during the evolution of the system. It is well known that the flowsgenerated by Hamiltonian vector fields in explicit form preserve the Hamiltonian func-tion, i.e. the total energy of the system is conserved during the evolution of the sys-tem [32, 23, 24]. As the explicit and implicit system representations are equivalent,the conservation also holds for the implicit Hamiltonian vector field (2.2), as is easilyconfirmed by computing the Lie derivative of H along the flow generated by X H,g [11].Indeed, L X H,g H = X H,g ( H ) = D H ( H ) − λ j ∂g j ∂r i ∂H∂p i = λ j f j , (2.9)where the third equality holds because D H ( H ) = 0 and because of the definition of f given in (2.5). Eq. (2.9) shows that L X H,g H (cid:12)(cid:12) L G = 0 (cf. Eq. (2.6)), so H remainsconstant along the system trajectories.It is also well known that, besides the 0-form H , Hamiltonian flows preserve acertain 2-form called the symplectic form. In the case of implicit representations thesymplectic form is given by the formula ω := d r i ∧ d p i , (2.10) See [1] for a coordinate-free definition.iscrete-Time Models for Implicit port-Hamiltonian Systems T ( T ∗ R n ), with Einstein’s summation convention implied. Definition 2.1.
A differentiable mapping φ : T ∗ R n → T ∗ R n is called symplectic if φ ∗ ω = ω . (2.11)Here, φ ∗ denotes the pull-back map associated with φ , defined by Definition 2.2.
Let φ : M → M be a smooth map of manifolds and let p ∈ M . The pull-back map φ ∗ : T ∗ φ ( p ) M → T ∗ p M associated with φ is a dual mapto the push-forward map φ ∗ and is characterized by h φ ∗ ξ, X i = h ξ, φ ∗ X i for ξ ∈ T ∗ φ ( p ) M , X ∈ T p M . The application of this definition permits to re-write (2.11) in the equivalent form ω ( φ ∗ ξ, φ ∗ η ) = ω ( ξ, η ) for all ξ, η ∈ T ( T ∗ R n ) . (2.12)The conservation of ω along X H,g can be established by showing that the Liederivative L X H,g ω is equal to zero, the demonstration bearing similarity to that of theconservation of H ; see [4] for the explicit derivation. Here, we cite the result Theorem 2.3. [11, 17] Let H be twice continuously differentiable. The flow φ t : L G → L G of X H,g governed by (2.2) is a symplectic transformation on L G , i.e., φ ∗ t ω = ω for every t for which φ t is defined. Remark 1.
The converse statement, that every symplectic flow φ t solves Hamil-ton’s equations for some H , is also true, so symplecticity is a characteristic propertyof Hamiltonian systems [1]. This does not translate to the case of energy conservation,i.e., while every Hamiltonian system conserves energy, not every energy-conservingsystem is Hamiltonian. In the presence of external forces and dissi-pation it is convenient to represent (2.2) as an input-output system equipped with apair of port variables ( u, y ), giving rise to a
PH system ; see [20, 32, 23] for the originaldefinition as stated with respect to Hamiltonian systems in explicit form. Extend-ing on this definition the Hamiltonian systems in implicit form, a port-Hamiltoniansystem is described in terms of the vector field X H,u,g : L G × ( R m ) ∗ → T ( L G ): X H,u,g = D H + (cid:18) u l U li − λ j ∂g j ∂r i (cid:19) ∂∂p i , g = 0 , (2.13)with u ∈ ( R m ) ∗ defined as the controlled or input variable, y ∈ R m is defined as theoutput variable that satisfies y l = U li ∂H∂p i (2.14)and where U li are maps from R n to R .The vector field (2.13) and the output unravel as the DAE system˙ r = ∇ p H ( r, p ) , ˙ p = −∇ r H ( r, p ) − G ( r ) ⊤ λ + U ( r ) uy = U ( r ) ⊤ ∇ p H ( r, p )0 = g ( r ) . Fernando Casta˜nos, Hannah Michalska, Dmitry Gromov and Vincent Hayward
By analogy with the results described in Sec. 2.2, one can determine the Lagrangemultipliers λ explicitly. The constraints f a = 0 imply that X H,u,g ( f a ) = D H ( f a ) + u l U li ∂f a ∂p i − λ j ∂g j ∂r i ∂f a ∂p i = 0 , (2.15)from where it follows that, as long as (2.8) is non-singular, there are unique λ j (ingeneral dependent on u as well as on r and p ) such that X H,u,g ( f a ) = 0 and such thatthe integral curve stays on L G .It can be readily seen that an implicit PH system described by (2.13) no longerpreserves H . The Lie derivative of H is now L X H,u,g ( H ) = X H,u,g ( H ) = D H ( H ) + u l U li ∂H∂p i − λ j ∂g j ∂r i ∂H∂p i = u l y l − λ j f j , which establishes the power balance L X H,u,g H (cid:12)(cid:12) L G = u l y l . (2.16)Since the product u l y l is equal to the rate of change in energy, we say that ( u, y ) is apower-conjugated pair of port variables. If, in addition, the restriction of H to L G isbounded from below, i.e., if the image of L G under H is bounded from below, then(2.13) is called passive, or more precisely, lossless . Boundedness of H can be easilyassessed using the following proposition. Proposition 2.4. [4] If the potential energy V is lower semi-continuous and G iscompact, then H restricted to L G is bounded from below (hence, the vector field (2.13) describes a lossless system). With the inclusion of the control variable u , it can no longer be expected thatthe flow of (2.13) be symplectic. Indeed, it is not hard to see that the Lie derivativeof ω along X H,u,g is in general different from zero.
Proposition 2.5. [4] The Lie derivative of ω (2.10) along (2.13) and restrictedto L G satisfies L X H,u,g ω (cid:12)(cid:12) L G = d r i ∧ d( u l U li ) . (2.17) Example: A double planar pendulum.
Let us recall an example from [4], onwhich we will elaborate when discussing sampled-data models.Consider the model of a double planar pendulum shown in Fig. 2.1 that comprisesa pair of point masses m a and m b whose coordinate positions are r a = ( r a x , r a y ) and r b = ( r b x , r b y ), respectively. The massless bars are of fixed lengths l a and l b whichgives rise to the two holonomic constraints: g ( r ) = k r a k − l a = 0 , g ( r ) = k r δ k − l b = 0 (2.18)with r := ( r a , r b ) ∈ R n , n = 4, k = 2, r δ := r b − r a . The rank of the constraintJacobian is full as rank G ( r ) = rank (cid:18) r a x r a y − r δ x − r δ y r δ x r δ y (cid:19) = k (2.19)for all r ∈ G . Therefore, 0 is a regular value of g and G is an embedded submanifoldof R . iscrete-Time Models for Implicit port-Hamiltonian Systems xy Fig. 2.1 . A double planar pendulum, a simple Hamiltonian system.
The total energy is H ( r, p ) = 12 p ⊤ M − p + ¯ g ( m a r a y + m a r a y ) , (2.20)where M := (cid:18) m a I n n n m b I n (cid:19) and ¯ g is the standard gravity.Substituting (2.18) and (2.20) in (2.4) gives˙ r a = m − a p a , ˙ r b = m − b p b (2.21a) ˙ p a x ˙ p a y ˙ p b x ˙ p b y = − gm a gm b − r a x − r δ x r a y − r δ y r δ x r δ y (cid:18) λ λ (cid:19) (2.21b)which, together with (2.18), constitutes a set of DAEs describing the motion of thedouble pendulum in implicit form. The multipliers λ and λ are the magnitudes ofthe tension along the two bars.Now, assume that the double pendulum is actuated by application of torques u and u to the joints that correspond to the angles q and q , respectively. Theresulting linear forces are then U u and U u with U := { U i } i = − r a y r a x l a and U := { U i } i = r δ y − r δ x − r δ y r δ x l b − U . (2.22)The manifold defined by (2.18) is compact and the potential energy is continu-ous, which confirms that the double pendulum is passive with passive outputs y l = U li ∂H∂p i = U li ˙ r i .An explicit model for the double pendulum as an ODE can be derived by choosingthe generalized coordinates on G as q ∈ ( − π, π ) and q ∈ ( − π, π ), as motivated byFig. 2.1. The associated embedding r = ı ( q ) satisfying g ◦ ı ≡ r a x r a y r b x r b y = l a cos q l a sin q l a cos q + l b cos q t l a sin q + l b sin q t , q t := q + q . (2.23) Fernando Casta˜nos, Hannah Michalska, Dmitry Gromov and Vincent Hayward
In local coordinates, the total energy isˆ H ( q, ˆ p ) = 12 ˆ p ⊤ ˆ M ( q ) − ˆ p + ¯ g (cid:0) m t l a sin q + m b l b sin q t (cid:1) . (2.24)with ˆ M ( q ) = (cid:18) m t l a + m b l b + 2 m b l a l b cos q m b l b + m b l a l b cos q m b l b + m b l a l b cos q m b l b (cid:19) The motion of the system is described by˙ q = ˆ M ( q ) − ˆ p , ˙ˆ p = −∇ q V ( q ) − ∇ q (cid:18)
12 ˆ p ⊤ ˆ M ( q ) − ˆ p (cid:19) + u (2.25)(see [4] for more details).Two representations for the same system were derived, one in the form of anODE (2.25) and the other as a DAE (2.21). The main point of this example can besummarized in the following remark. Remark 2.
The implicit DAE representation of the Hamiltonian systems rendersa separable Hamiltonian function (2.20) , i.e. the kinetic energy does not depend on r and the inertia matrix M is a constant diagonal matrix. Moreover, the potentialenergy is linear , which results in a constant potential energy gradient. On the otherhand it is easily verified that the explicit Hamiltonian representation does not sharesuch fortunate properties. The explicit Hamiltonian is no longer separable and theinertia matrix depends on the generalized coordinates q . The a cost of a simpler Hamiltonian function is the higher dimensional implicitmodel representation and the appearance of the Lagrange multipliers. As will soon beseen, however, implicit representations are particularly advantageous for the purposesof system discretization.
It was assumed in (2.13) thatthe control variables enter the vector field affinely. Many physical systems exhibitthis property so, from a modeling point of view, this assumption is not too restrictive.However, for the purpose of performing backward error analysis in §
4, we will need toconsider PH systems for which the control might enter in a non affine way. Motivatedby the fact that L X H,u,g ω (cid:12)(cid:12) L G = 0 when u ≡
0, we propose the following extendeddefinition of a PH system.
Definition 2.6.
A controlled vector field X : L G × ( R m ) ∗ → T ( L G ) is said tobe port-Hamiltonian if u ≡ implies that the generated flow is symplectic. Passivity of a non affine port-Hamiltonian system can be established by redefiningthe passive output.
Lemma 2.7.
A (not necessarily affine) smooth PH system described by a vectorfield X can always be decomposed as X = X + u l Z l , where X : L G → T ( L G ) is aHamiltonian vector field and Z l : L G × ( R m ) ∗ → T ( L G ) are the input vector fields.Hence, X satisfies the power balance X ( H ) = u l y l for some real-valued function H and real-valued output functions y l = Z l ( H ) . (The output functions may now dependdirectly on u as well as on r and p .)Proof . Following [18], we first show that a smooth control vector field can be splitinto a drift and a set of vector fields having u factored out. Let us define the drift X as X ( x ) = X ( x,
0) and let us define the vector fields W l by the equations W l ( α ) = ∂∂u l (cid:16) X ( α ) (cid:17) for all α ∈ C ∞ ( L G , R ) . iscrete-Time Models for Implicit port-Hamiltonian Systems u l W l ( x, θu )( α ) = dd θ (cid:16) X ( x, θu ) (cid:17) ( α ) . Upon integration on both sides of the equation we arrive at u l Z W l ( x, θu )( α )d θ = X ( x, θu )( α ) (cid:12)(cid:12)(cid:12) = X ( x, u )( α ) − X ( x )( α ) . Therefore, we have X = X + u l Z l , (2.26)where the input vector fields Z l are defined by Z l ( x, u )( α ) = Z W l ( x, θu )( α )d θ for all α ∈ C ∞ ( L G , R ) . It follows from the hypothesis that X generates a symplectic flow, so it is aHamiltonian vector field and satisfies X ( H ) = 0 for some real-valued function H .Applying X to H shows that X ( H ) = X ( H ) + u l Z l ( H ) = u l y l with y l = Z l ( H ). Remark 3.
For an affine PH system, the formulae of this lemma recover theoutput functions (2.14) with Z l = U lj ∂∂p j , that is, y l = Z l ( H ) = U li ∂H∂p i .
3. Sampled-data models for autonomous Hamiltonian systems.
Com-puting a sampled-data model of a dynamical system basically amounts to computingan approximate solution of the differential equations during a small interval of time.This problem has been studied extensively in the literature of numerical analysis, fromwhich we borrow some results and terminology. In numerical analysis, a sampled-datamodel is the central component of an integration method or a numerical integrator ,up to the point that these terms are used interchangeably.Mathematical models for sampled-data systems arise in diverse circumstances. Inthe direct approach to digital control, i.e., as opposed to the emulation of continuouscontrol laws, the design of the controller is performed in discrete time, the designerworking directly over a sampled-data model. When designing directly in discrete time,the controller can be directly implemented on a digital device. Also, it is possible toexploit the advantages of switched controls or, e.g., multirate control techniques [21].Computing the sampled-data model for a given nonlinear system relies on thecomputation of a solution φ t of the corresponding ODE or DAE, which is in generalimpossible to do analytically, so one has to settle for an approximate solution.When simulating the behavior of dynamic systems, a discrete-time model of thecontinuous system is also used for computing a numerical solution to the initial-valueproblem. Many different integration methods (or methods for short) can be found inthe literature. Let us first focus on integration methods for autonomous systems andfurther restrict our attention to one-step methods defined by a transformation ψ h : x α x α +1 , where the constant step-size h is regarded as a parameter of the method . For a giveninitial condition x in the phase space, ψ h is applied recursively to generate a discrete For a more general method, the value of the x α +1 need not depend only on x α , but may alsodepend on the previous values x α − , x α − ,. . . (a multistep method). Also, the value of h need notbe constant in general. Fernando Casta˜nos, Hannah Michalska, Dmitry Gromov and Vincent Hayward flow x , x , x , . . . that approximates the true flow φ h ( x ) , φ h ( x ), φ h ( x ) , . . . ofa given vector field X at time instants h, h, h, . . . In this sense, the map ψ h is adiscrete-time approximation of φ h (or a sampled-data model of φ h ). Definition 3.1.
A one-step method has order s if the local error satisfies ψ h ( x ) − φ h ( x ) = O ( h s +1 ) as h → uniformly in x . A one-step method is said to be consistent if s ≥ . Let us now discuss some important properties of numerical integrators, like orderand symmetry.
If a sampled-data model approximates the discretetime behavior of a Hamiltonian system, one could hope for ψ h to inherit its funda-mental qualitative properties: energy conservation and symplecticity. Unfortunately,it is not possible to preserve H and ω simultaneously, unless ψ h agrees with the exactflow φ h up to a reparametrization of time [7]. For this reason, one has to choose ei-ther in favor of one or the other invariant . Energy conserving methods have receivedsome attention [9, 31, 10, 12, 13], but in light of Remark 1, most of the literaturefocuses on symplectic integration algorithms (see [11, 17, 5] and references therein).A comparison between both approaches is carried out in [30] for the rigid body.A theoretical advantage of constructing a symplectic one-step method is that,even though ψ h only approximates φ h up to the s ’th order, it coincides exactly (ifone disregards convergence issues) with the flow of another Hamiltonian system, amodified Hamiltonian system described by a modified differential equation . Theorem 3.2. [11, p. 352] A symplectic method ψ h : L G → L G for the con-strained Hamiltonian system (2.4) has a modified equation that is locally of the form ˙ r = + ∇ p ˜ H ( r, p ) , ˙ p = −∇ r ˜ H ( r, p ) − G ( r ) ⊤ ˜ λ (3.2a)0 = g ( r ) (3.2b) with ˜ H = H + hH + h H + . . . Furthermore, ∂H j ( r, p ) ∂p i ∂g l ( r ) ∂r i = 0 , for all ( r, p ) ∈ L G , all l = 1 , . . . , k and all j . Note that the actual value of the Legendre multipliers ˜ λ differ in general from those obtained for the original system (2.4). In other words, for an initial condition x , ψ h ( x ) is equal to the solution of (3.2)at time t = h . Note that (3.1) provides information about the difference betweenthe actual flow φ h of X and the approximate discrete flow ψ h . This is the kindof information that forward error analysis aims at. While certainly useful as anindicator of the quality of the approximation, Eq. (3.1) only evaluates the behaviorof the approximate flow on the first iteration, but says nothing about its long timebehavior. From (3.1) alone we cannot infer anything about the error x α − φ αh ( x )when α is large, so we do not know if errors accumulate or if they average out to zero. We use big-O notation when quantifying approximation errors, i.e., for a given pair of functions e ( h ), e ( h ), we write e ( h ) = O ( e ( h )) as h → a as shorthand for lim sup h → a k e ( h ) kk e ( h ) k < ∞ . For particular Hamiltonians there might be other invariants, such as momentum or angularmomentum, but in general there need not be.iscrete-Time Models for Implicit port-Hamiltonian Systems ψ h is symplectic, then thereexists a modified continuous system whose flow coincides exactly with the discreteflow generated by ψ h . The modified system (3.2) preserves the Hamiltonian structureof the original system (2.4) and it is ‘close’ to it in the sense that ˜ H = H + O ( h s )for a method of order s . In other words, a symplectic integration method preservesthe original 2-form ω and a different (but close) Hamiltonian function. This propertyguaranties that the good behavior of the integration scheme is maintained duringmany iterations, giving a global nature to the local property (3.1). This observationis at the center of backward error analysis [11]. A practical advantage of symplectic schemes is thatthey lend themselves well to the application of splitting methods. To illustrate theidea, consider again the unconstrained or, otherwise, explicit Hamiltonian vector field D ˆ H on an abstract manifold T ∗ G . If the Hamiltonian function is separable, i.e., if itcan be written as ˆ H ( q, ˆ p ) = ˆ H a ( q ) + ˆ H b (ˆ p ), then the vector field can be spilt into twoHamiltonian vector fields D ˆ H a = − ∂ ˆ H a ∂q i ∂∂ ˆ p i and D ˆ H b = ∂ ˆ H b ∂ ˆ p i ∂∂q i with D ˆ H = D ˆ H a + D ˆ H b . Notice that, taken separately, each vector field can be triviallyintegrated exactly. For ( q α , ˆ p α ) ∈ T ∗ G (we use the subindex α to refer to an elementin a sequence, not a particular coordinate) we have (cid:18) ˙ q ˙ˆ p (cid:19) = (cid:18) −∇ q ˆ H a ( q ) (cid:19) = ⇒ (cid:18) q α +1 ˆ p α +1 (cid:19) = (cid:18) q α ˆ p α − h · ∇ q ˆ H a ( q α ) (cid:19) = φ a ,h (cid:18) q α ˆ p α (cid:19) and (cid:18) ˙ q ˙ˆ p (cid:19) = (cid:18) + ∇ p ˆ H b (ˆ p )0 (cid:19) = ⇒ (cid:18) q α +1 ˆ p α +1 (cid:19) = (cid:18) q α + h · ∇ ˆ p ˆ H b (ˆ p α )ˆ p α (cid:19) = φ b ,h (cid:18) q α ˆ p α (cid:19) . A first-order symplectic method can be easily constructed by performing the compo-sition ψ h = φ b ,h ◦ φ a ,h . (3.3)Indeed, the maps φ b ,h and φ a ,h are symplectic because they are the exact flows ofHamiltonian vector fields. Since the composition of two symplectic maps is againsymplectic, ψ h is symplectic.Many simple Hamiltonian systems with phase space T ∗ G are not governed byseparable Hamiltonians, so splitting methods cannot be applied directly. However,the Hamiltonian function of many mechanical systems becomes separable if the phasespace is embedded in T ∗ R n (see, e.g., Remark 2). An interesting symplectic methodthat is particularly well suited for this class of systems was proposed in [27]. Roughlyspeaking, the idea is to compute a symplectic method for the unconstrained Hamilto-nian vector field D H : T ∗ R n → T ( T ∗ R n ), i.e., a symplectic map ψ H,h approximatingthe solution of the ODE (notice the absence of the constraint equations)˙ r = + ∇ p H ( r, p )˙ p = −∇ r H ( r, p )at time t = h .2 Fernando Casta˜nos, Hannah Michalska, Dmitry Gromov and Vincent Hayward If H is separable, ψ H,h can be readily found. The method Ψ
H,g,h for the originalconstrained Hamiltonian vector field X H,g is then constructed by taking the image of ψ H,h and applying a correction term that ensures that the value of Ψ
H,g,h belongs to L G , so that the constraints are satisfied. The correction is done in a careful way sothat the resulting map is still symplectic (see also [17]). Depending on the accuracyof ψ H,h , the resulting Ψ
H,g,h can be of first or second order (see § For each t for which the solution is defined, theflow φ t ( x ) of an autonomous differential equation defines a transformation on thephase space. It follows from the group property of the flow [2] that the inverse ofthe transformation can be obtained simply by reversing time, that is, φ − t ( x ) = φ − t ( x ) = x . Needless to say, this property does not hold in general for a discreteapproximation ψ h , which motivates the following definition. Definition 3.3.
The adjoint method ψ ∗ h of a method ψ h is the inverse map ofthe original method with reversed time step − h , i.e., ψ ∗ h := ( ψ − h ) − . In other words, x = ψ ∗ h ( x ) is implicitly defined by ψ − h ( x ) = x . A method forwhich ψ ∗ h = ψ h is called symmetric . From a theoretical point of view, an approximate discrete-time flow should besymmetric because actual continuous flows are. But symmetry is important from apractical point of view too. It has been proved in [34] that all symmetric methodsare of even order, a fact that can be exploited to construct high-order methods fromsimple lower-order methods. For example, one can take a first-order non-symmetricmethod, compute its adjoint and construct a symmetric methodΨ h = ψ h ◦ ψ ∗ h . (3.4)We know that Ψ h is at least first order, but since Ψ h is symmetric, we also know thatthe order has to be even, so we conclude that the method is actually of second order.The scheme (3.4) works particularly well for splitting methods. Take, e.g., theintegration scheme (3.3). The maps φ b ,h and φ a ,h are symmetric (because they areexact solutions of a differential equation), but their composition is not symmetric ingeneral. To remedy this, one can compute the adjoint method ψ ∗ h = ( φ b , − h ◦ φ a , − h ) − = φ − , − h ◦ φ − , − h = φ a ,h ◦ φ b ,h (3.5)and, using (3.3), (3.5) and (3.4), constructΨ h = φ b , h ◦ φ a , h ◦ φ a , h ◦ φ b , h = φ b , h ◦ φ a ,h ◦ φ b , h , which is a second-order symmetric method. Consider avector-valued function F and a vector field X , both defined on L G . If F and X areanalytic, then the composition of F and the generated flow φ t ( x ) can be expandedin a Taylor series around t = 0, F ◦ φ t ( x ) = exp( tX ) F ( x ) := ∞ X i =0 t i i ! X i ( F )( x ) , iscrete-Time Models for Implicit port-Hamiltonian Systems ZOH Continuous PH System Sampler
Sampled-data PH system
Fig. 4.1 . Sampled-data PH system with sampling period h . The zero-order hold produces piece-wise constant inputs u l ( t ) ≡ u lα for t ∈ [ αh, αh + h ) , which is fed to the continuous-time PH system.The output is then sampled to generate the discrete-time output sequences { y lα } . where X ( F ) = F , X ( F ) = X ( X ( F )), X ( F ) = X ( X ( F )) etc (see [22, 33] fordetails). In particular, if F is taken as the identity function Id, one obtains the flow φ t ( x ) = exp( tX ) Id( x ). Since an s -order method ψ h for X coincides with the flow ofa modified vector field ˜ X = X + O ( h s ) [11, p. 340], it is also possible to expand ψ h ina Taylor series, ψ h ( x ) = exp( h ˜ X ) Id( x ). This exponential notation is a convenientway to express the relationship between a vector field and the flow generated by it,as well as to analyze the composition of flows.
4. A splitting method for implicit port-Hamiltonian systems.
Supposethat there is a sequence of commands { u lα } α ∈ N . Each command in the sequence ar-rives at the discrete instants of time α = 0 , h, h, . . . where h is a positive real number—such commands could be generated, e.g., by a computer program. Suppose furtherthat a zero-order hold transforms this sequence into piece-wise constant controls u l ( t ) ≡ u lα for t ∈ [ αh, αh + h ) , (4.1)which are fed into a PH system. Let φ t ( x , u ( · )) be the integral curve of the non-autonomous vector field (2.13) with control (4.1) and passing through x ∈ L G at t = 0. Let y l ( t ) = (cid:18) U li ∂H∂p i (cid:19) ◦ φ t ( x , u ( · ))be the corresponding outputs and let { y lα } with y lα := y l ( αh ) be the sequence ob-tained by sampling them at discrete instants of time αh (see Fig. 4.1). We call theresulting system a sampled-data port-Hamiltonian system .The goal here is to develop a method for derivation of discrete-time (or sampled-data) models for PH systems given by implicit vector fields. The underlying ideais to split the PH vector field into two components: the vector field describing anunconstrained system with state-space equal to the whole T ∗ R n , and a vector fieldcontaining the Lagrange multipliers, the one that maintains the trajectories on thesubmanifold L G . Splitting the vector field simplifies the computation of the sampled-data model by decomposing the problem into two simpler subproblems.Now we extend the results of [27] to the PH case. We show that, with a straightfor-ward modification, the method presented in [27], originally intended as an integrationscheme for autonomous Hamiltonian systems, can be used to compute sampled-datamodels that preserve the main properties of a PH system.Consider again the implicit vector field X H,u,g = D H + (cid:18) u l U li − λ j ∂g j ∂r i (cid:19) ∂∂p i , g = 0 , y l = U li ∂H∂p i , (4.2)4 Fernando Casta˜nos, Hannah Michalska, Dmitry Gromov and Vincent Hayward defined on L G , with D H as in (2.3) and with piece-wise constant controls (4.1). Sup-pose that a method ψ H,u,h : T ∗ R n → T ∗ R n of order s ≥ X H,u = D H + u l U li ∂∂p i has been computed. Again, in many cases H isseparable so a high-order and symmetric method with a PH modified vector field canbe easily found. The controls u l are constant during each sampling interval, whichfurther simplifies the task of finding ψ H,u,h .Let us define the mapΦ Λ ,h (cid:18) r α p α (cid:19) := (cid:18) r α p α − hG ( r α ) ⊤ Λ (cid:19) . (4.3)Loosely, this is an approximation of the integral curve of the remnant vector field − λ j ∂g j ∂r i ∂∂p i evaluated at t = h and subject to g = 0. More precisely, for arbitraryfunctions λ j of r and p , we have that D λ j g j := ∂ ( λ j g j ) ∂p i ∂∂r i − ∂ ( λ j g j ) ∂r i ∂∂p i = g j ∂λ j ∂p i ∂∂r i − (cid:18) λ j ∂g j ∂r i + g j ∂λ j ∂r i (cid:19) ∂∂p i and, for r ∈ G , the vector field reduces to D λ j g j = λ j D g j = − λ j ∂g j ∂r i ∂∂p i . (4.4)In other words, when restricted to G , the vector field − λ j ∂g j ∂r i ∂∂p i is Hamiltonian (henceit generates a symplectic flow). Lemma 4.1. [27] Let g ( r α ) = 0 . Then, the map (4.3) is a first-order symplecticmethod for D Λ j g j . That is, for r α ∈ G , Φ Λ ,h (cid:18) r α p α (cid:19) = exp( h ˜ D Λ j g j ) Id (cid:18) r α p α (cid:19) , ˜ D Λ j g j = D ˜Λ j g j with ˜Λ a modified or perturbed version of Λ . A method for (4.2) can be obtained from the symmetric compositionΨ
H,u,g,h = Φ µ, h ◦ ψ H,u,h ◦ Φ ν, h . (4.5)For each ( r α , p α ) ∈ L G , the values of µ and ν are determined implicitly by theconstraints g ( r α +1 ) = 0 and f ( r α +1 , p α +1 ) = 0 (i.e., by ( r α +1 , p α +1 ) ∈ L G ). In thisway, Ψ H,u,g,h defines a transformation on L G .The transformation Ψ H,u,g,h produces an approximate discrete flow for a givencommand sequence { u lα } . From this flow, an approximate output sequence can beobtained by evaluating the output function y l = U li ∂H∂p i at each discrete time αh . Theorem 4.2.
Consider the implicit method Ψ H,u,g,h (4.5) and let ˜ X H,u be themodified vector field of ψ H,u,h . (i) The method preserves the constraints g j = 0 , f j = 0 and is of order ¯ s =min( s, , where s is the order of ψ H,u,h . (ii) The method is symmetric if ψ H,u,h is symmetric. (iii) If ψ H,u,h is symplectic for u lα ≡ (i.e., if ˜ X H,u is port-Hamiltonian), thenthe modified vector field ˜ X H,u,g : L G → T ( L G ) is also port-Hamiltonian with Hamil-tonian and output functions ˜ H = H + O ( h ¯ s ) and ˜ y l = y l + O ( h ¯ s ) . (4.6) iscrete-Time Models for Implicit port-Hamiltonian Systems Proof . The method preserves the constraints by construction. The proof aboutthe order of the method follows the same lines as the one given in [27] except that,since we are dealing with port-Hamiltonian vector fields, Lie brackets have to beused instead of Poisson brackets. We will compute ˜ X H,u,g , the modified vector fieldgenerating Ψ
H,u,g,h , and show that it agrees with X H,u,g up to the first or secondorder, depending on whether Ψ
H,u,g,h is, respectively, first or second order.Let us consider the case s = 1. Using the exponential notation and Lemma 4.1,the composition (4.5) takes the formΨ H,u,g,h = exp (cid:18) h D ˜ ν j g j (cid:19) exp (cid:16) h ˜ X H,u (cid:17) exp (cid:18) h D ˜ µ j g j (cid:19) Id , where ˜ X H,u = X H,u + O ( h ). Applying the Baker-Campbell-Hausdorff (BCH) for-mula [33] to the product of the first two factors and truncating after the first termgives exp (cid:18) h D ˜ ν j g j (cid:19) exp (cid:16) h ˜ X H,u (cid:17) = exp (cid:16) h ˜ X ′ (cid:17) (4.7)with ˜ X ′ = ˜ X H,u + D ˜ ν j g j + O ( h ). Applying BCH again to include the third factorgives exp (cid:16) h ˜ X ′ (cid:17) exp (cid:18) h D ˜ µ j g j (cid:19) = exp (cid:16) h ˜ X H,u,g (cid:17) (4.8)with the modified vector field ˜ X H,u,g = ˜ X H,u + (cid:0) D ˜ ν j g j + D ˜ µ j g j (cid:1) + O ( h ). Using (4.4)and ˜ X H,u = X H,u + O ( h ), we can write the modified vector field as˜ X H,u,g = X H,u + ˜ ν j + ˜ µ j D g j + O ( h ) . (4.9)The hidden constraints f l = 0 imply that˜ X H,u,g ( f l ) = X H,u ( f l ) + ˜ ν j + ˜ µ j D g j ( f l ) + O ( h ) = 0 . (4.10)It follows from (4.10), (4.4) and (2.15), that the Lagrange multipliers λ j and the‘modified Lagrange multipliers’ ˜ ν j and ˜ µ j are related by the equation (˜ ν j + ˜ µ j ) / λ j + O ( h ), which when substituted back in (4.9) gives the desired result:˜ X H,u,g = X H,u + λ j D g j + O ( h ) = X H,u,g + O ( h ) . For s = 2 we follow the same procedure, but we truncate the BCH formula afterthe second term. For the expression (4.7), the intermediate vector field is˜ X ′ = ˜ X H,u + 12 D ˜ ν j g j + h h D ˜ ν j g j , ˜ X H,u i + O ( h )where [ · , · ] is the standard Lie bracket. Using the initial assumption ˜ X H,u = X H,u + O ( h ), we can write ˜ X ′ as˜ X ′ = X H,u + 12 D ˜ ν j g j + h (cid:2) D ˜ ν j g j , X H,u (cid:3) + O ( h ) . Fernando Casta˜nos, Hannah Michalska, Dmitry Gromov and Vincent Hayward
Regarding (4.8), the modified vector field for the complete scheme is˜ X H,u,g = X H,u + 12 (cid:0) D ˜ ν j g j + D ˜ µ j g j (cid:1) + h (cid:2) D ˜ ν j g j , X H,u (cid:3) + h (cid:20) X H,u + 12 D ˜ ν j g j , D ˜ µ j g j (cid:21) + O ( h ) . Using (4.4) and the skew symmetry and bilinearity of the Lie bracket, the vector fieldcan be equivalently written as˜ X H,u,g = X H,u + ˜ ν j + ˜ µ j D g j + h (cid:2) (˜ ν j − ˜ µ j ) D g j , X H,u (cid:3) + h (cid:2) ˜ ν j D g j , ˜ µ j D g j (cid:3) + O ( h ) . (4.11)In order to extract information from the equation ˜ X H,u,g ( g l ) = 0, let us first openthe brackets in (4.11) and write˜ X H,u,g = X H,u + ˜ ν j + ˜ µ j D g j + h ν j − ˜ µ j ) D g j X H,u − h X H,u (cid:0) (˜ ν j − ˜ µ j ) D g j (cid:1) + h ν j D g j (cid:0) ˜ µ j D g j (cid:1) − h µ j D g j (cid:0) ˜ ν j D g j (cid:1) + O ( h ) . Taking into account that f l = X H,u ( g l ) = 0 and D g j ( g l ) ≡
0, we have that˜ X H,u,g ( g l ) = h ν j − ˜ µ j ) D g j f l + O ( h ) = 0 , from which we can see that modified Lagrange multipliers satisfy the order relation˜ ν j − ˜ µ j = O ( h ) . (4.12)By substituting (4.12) back in (4.11) we can verify that the commutators are actuallysecond order, that is,˜ X H,u,g = X H,u + ˜ ν j + ˜ µ j D g j + h (cid:2) ˜ ν j D g j , (˜ ν j + O ( h )) D g j (cid:3) + O ( h )= X H,u + ˜ ν j + ˜ µ j D g j + O ( h ) . (4.13)From ˜ X H,u,g ( f l ) = 0 and (2.15) we conclude that, when s = 2, (˜ ν j + ˜ µ j ) / λ j + O ( h ), so the desired result follows: ˜ X H,u,g = X H,u + λ j D g j + O ( h ) = X H,u,g + O ( h ).For statement (ii), notice that, when restricted to L G , the method (4.5) can bedescribed by the implicit equations (cid:18) r α +1 p α +1 + hG ( r α +1 ) ⊤ µ (cid:19) = ψ H,u,h (cid:18) r α p α − hG ( r α ) ⊤ ν (cid:19) (4.14a) g ( r α +1 ) = g ( r α ) (4.14b) f ( r α +1 , p α +1 ) = f ( r α , p α ) , (4.14c)where r α , p α are the independent variables and r α +1 , p α +1 are the dependent variables.The vectors ν, µ are (also dependent) dummy variables that can be discarded after r α +1 , p α +1 have been found. iscrete-Time Models for Implicit port-Hamiltonian Systems h by − h ), equation (4.14a) be-comes (cid:18) r α +1 p α +1 − hG ( r α +1 ) ⊤ µ (cid:19) = ψ H,u, − h (cid:18) r α p α + hG ( r α ) ⊤ ν (cid:19) . Recall that ψ H,u, − h = ψ − H,u,h if ψ H,u,h is symmetric. Therefore, when restricted to L G , the reverse-time method is (cid:18) r α p α + hG ( r α ) ⊤ ν (cid:19) = ψ H,u,h (cid:18) r α +1 p α +1 − hG ( r α + 1) ⊤ µ (cid:19) (4.15a) g ( r α +1 ) = g ( r α ) (4.15b) f ( r α +1 , p α +1 ) = f ( r α , p α ) , (4.15c)which is the same as (4.14), but with r α , p α and ν interchanged with r α +1 , p α +1 and µ , respectively. This implies that, if we input r α +1 , p α +1 as independent variables, werecover r α , p α as the dependent variables, that is: Ψ H,u,g, − h is the inverse mapping ofΨ H,u,g,h . (In general, the vectors ν, µ obtained using (4.14) will be different from thoseobtained using (4.15), but this is inconsequential since these are dummy variables.)In statement (iii), the fact that ˜ X H,u,g is PH follows directly from Lemma 4.1,Definition 2.6 and the fact that the composition of symplectic maps is again symplec-tic. In other words, Ψ
H,u,g,h is symplectic when u lα ≡
0, so˜ X H,u,g = X + u l Z l = D ˜ H + u l Z l + ˜ λ j D g j (4.16)for some Hamiltonian function ˜ H and some input vector fields Z l . Since the methodis of order ¯ s , we have˜ X H,u,g = D H + u l U li ∂H∂p i + λ j D g j + O ( h ¯ s ) . (4.17)By setting u l = 0 and recalling that ˜ λ j = λ j + O ( h ¯ s ), it follows from (4.16) and (4.17)that D ˜ H = D H + O ( h ¯ s ), which implies that ˜ H = H + O ( h ¯ s ), and, in turn, that Z l = U li ∂∂p i + O ( h ¯ s ). The modified output functions are thus ˜ y l = Z l ( ˜ H ), that is,˜ y l = U li ∂∂p i (cid:0) H + O ( h ¯ s ) (cid:1) + O ( h ¯ s ) = y l + O ( h ¯ s ).Let us now turn to the problem of energy balance under sample and hold. Thepower balance (2.16) implies that H α +1 − H α = Z αh + hαh u l ( t ) y l ( t )d t , (4.18)where we have defined the sampled Hamiltonian H α := H ( x α ). A usual way to im-prove the transient behavior of the system is to add damping by means of a continuouscontrol law [23] u l ( t ) = − K lj y j ( t ) , (4.19)with { K lj } a symmetric and positive semi-definite matrix. With the control law(4.19), the power balance (4.18) results in the dissipation inequality H α +1 − H α ≤ H α decreases monotonically and, if the right conditions aremet, the system converges to a state of minimal energy.8 Fernando Casta˜nos, Hannah Michalska, Dmitry Gromov and Vincent Hayward
Suppose that the output is being sampled and that the input is being held atintervals of length h . The control sequence is then given by u lα = − K lj y jα (4.20)and the power balance (4.18) takes the form H α +1 − H α = u lα Z αh + hαh y l ( t )d t = u lα Z h y l ( αh + τ )d τ . Applying Taylor’s theorem to the integral term gives H α +1 − H α = X l u lα (cid:0) y lα h + O ( h ) (cid:1) = − hK lj y lα y jα + X l u lα O ( h ) , so H α decreases when h is small enough and the norm of y α is large enough.Since the approximate sampled-data model (4.5) is also PH (cf. item (iii) ofTheorem 4.2), it satisfies (again, after applying Taylor’s theorem)˜ H α +1 − ˜ H α = X l u lα (cid:0) ˜ y lα h + O ( h ) (cid:1) (4.21)for some ˜ H and ˜ y . According to (4.6), the energy balance (4.21) takes the form H α +1 − H α + O ( h ¯ s ) = X l u lα (cid:0)(cid:0) y lα + O ( h ¯ s ) (cid:1) h + O ( h ) (cid:1) . The same control sequence (4.20) produces H α +1 − H α = − hK lj y lα y jα + X l u lα O ( h ) + O ( h ¯ s ) . Thus, for ¯ s = 2, the qualitative behavior of the approximated sampled data model isthe same as the exact one: H α decreases when h is small enough and the norm of y α is large enough. Example: A double planar pendulum (continued).
Let us compute asampled-data model for the double pendulum described in the previous examples.The first step is to compute a sample-data model for the simple unconstrained PHsystem X H,u = D H + u l U li ∂∂p i , where H is given by (2.20) and U li by (2.22). Theunconstrained and unactuated Hamiltonian vector field D H describes a pair of masseswith initial positions r a and r b and initial momenta p a p b , simply falling underthe influence of gravity. The exact flow generated by D H , the drift, denoted by( r α +1 , p α +1 ) = φ H,h ( r α , p α ), is then given by r a x α +1 = r a x α + hm a p a x α , r a y α +1 = r a y α + hm a p a y α − ¯ g h r b x α +1 = r b x α + hm b p b x α , r b y α +1 = r b y α + hm b p b y α − ¯ g h p a x α +1 = p a x α , p a y α +1 = p a y α − m a ¯ ghp b x α +1 = p b x α , p b y α +1 = p b y α − m b ¯ gh . iscrete-Time Models for Implicit port-Hamiltonian Systems Value Description l a = 0 . l b = 0 . m a = 0 . m b = 0 . g = 9 . / s] Acceleration due to gravity
Table 4.1
Parameters for the double pendulum.
The exact flow generated by u l U li ∂∂p i , the control vector field without drift, is denotedby ( r α +1 , p α +1 ) = φ u,h ( r α , p α ). It is given by r iα +1 = r iα , i ∈ { a x , a y , b x , b y } and p a x α +1 = p a x α − hl a ( u − u ) r a y α , p a y α +1 = p a y α + hl a ( u − u ) r a x α p b x α +1 = p b x α − hl b u r δ y α , p b y α +1 = p b y α + hl b u r δ x α . From § X H,u is ψ H,u,h = φ H, h ◦ φ u,h ◦ φ H, h . (4.22)Notice that φ u,h = Id when u ≡
0, so ψ H,u,h = φ H, h ◦ φ H, h = φ H,h , which is asymplectic map because it is the exact solution of a Hamiltonian system. Therefore,˜ X H,u is PH, and, from Theorem 4.2, it follows that the implicit method (4.5), withΨ
H,u,h as in (4.22) is the exact solution of a PH system ˜ X H,u,g with Hamiltonianfunction ˜ H = H + O ( h ) and output function ˜ y = y + O ( h ) (i.e., ¯ s = s = 2).The sampled-data model was tested using the parameters shown in Table 4.1.For illustration purposes, we chose a damping control u α = − . · y α and simulatedthe closed-loop system using the sampled-data model (4.5). Figure 4.2 shows thediscrete-time series of H for different values of h . It can be seen that the time seriesconverge and, as expected, the value of H α decreases monotonically when h is smallenough (in this case, less or equal to 30 [ms]). For comparison purposes, we haveincluded the evolution of H that is obtained by simulating (with Matlab’s Simulink)the explicit model developed in [4] in series with a sampler and a zero-order hold.
5. Conclusions.
We have extended the second-order integration method pre-sented in [27]. The original method applies to autonomous Hamiltonian systems and,being symplectic, preserves the Hamiltonian structure of the continuous-time system.The extended method can be applied to port-Hamiltonian systems, which are Hamil-tonian systems equipped with input–output pairs. The extended method preservesthe port-Hamiltonian structure. Affinity in the controls is lost by the method but,fortunately, the passivity properties of the continuous-time system can be recoveredby a suitable redefinition of the output. Interestingly, the relation between the originaland the new output is also of second order.The integration method can be used with the purposes of numerical simulationor with the purpose of deriving discrete-time models to be used in the design ofdiscrete-time control laws.0
Fernando Casta˜nos, Hannah Michalska, Dmitry Gromov and Vincent Hayward
Discrete flowsSimulink solution for explicit vector fieldImplicit discrete-time model h = 0.0171 h = 0.0555 h = 0.0005 h = 0.0016 h = 0.0029 h = 0.0053 Time [s] H a m i l t o n i a n [ J ] h = 0.0308 Fig. 4.2 . Results of the numerical experiment. The Hamiltonian function is plotted againsttime. The explicit model (simulated using Matlab’s module Simulink) is compared with the implicitmodel (simulated using a Matlab script). As expected, H is monotonically decreasing when h issmall enough and the time series converges as h goes to zero. REFERENCES[1]
Vladimir I. Arnold , Mathematical methods of classical mechanics , Springer-Verlag, NewYork, 1989.[2] ,
Ordinary Differential Equations , Springer-Verlag, Berlin, 1992.[3]
Vladimir I. Arnold, Valery Kozlov, and Anatoly I. Neishtadt , Mathematical aspects ofclassical and celestial mechanics , Springer-Verlag, 2006.[4]
Fernando Casta˜nos, Dmitry Gromov, Vincent Hayward, and Hannah Michalska , Im-plicit and explicit representations of continuous-time port-Hamiltonian systems , Systemsand Control Lett., 62 (2013), pp. 324 – 330.[5]
P. J. Channell and F. R. Neri , Symplectic integrators , in Integration Algorithms and Clas-sical Mechanics, Jerrold E. Marsden, George W. Patrick, and William F. Shadwick, eds.,AMS, Rhode Island, 1996, pp. 45 – 57.[6]
Morten Dalsmo and Arjan J. van der Schaft , On representations and integrability ofmathematical structures in energy-conserving physical systems , SIAM J. Control Optim.,37 (1999), pp. 54 – 91.[7]
Zhong Ge and Jerrold E. Marsden , Lie-Poisson Hamilton-Jacobi theory and Lie-Poissonintegrators , Physical Review Letters A, 133 (1988), pp. 134 – 139.[8]
Herbert Goldstein, Charles P. Poole, and John L. Safko , Classical Mechanics , AddisonWesley, 2002.[9]
Oscar Gonz´alez , Design and analysis of conserving integrators for nonlinear Hamiltoniansystems with symmetry , PhD thesis, Stanford University, 1996.[10]
Oscar Gonz´alez , Mechanical systems subject to holonomic constraints: Differentialalgebraicformulations and conservative integration , Physica D, 132 (1999), pp. 165 – 174.[11]
Ernst Hairer, Christian Lubich, and Gerhard Wanner , Geometric numerical integration:structure-preserving algorithms for ordinary differential equations , Springer-Verlag, 2006.[12]
Robert A. LaBudde and Donald Greenspan , Energy and momentum conserving methods ofarbitrary order for the numerical integration of equations of motion. I. motion of a singleparticle , Numer. Math., 25 (1976), pp. 323 – 346.[13] ,
Energy and momentum conserving methods of arbitrary order for the numerical in-tegration of equations of motion. II. motion of a system of particles , Numer. Math., 26(1976), pp. 1 – 16.[14]
Dina Shona Laila and Alessandro Astolfi , On the construction of discrete-time modelfor port-controlled Hamiltonian systems with applications , Systems and Control Lett., 55(2006), pp. 673 – 680.[15] ,
Direct discrete-time design for sampled-data Hamiltonian control systems , in La-iscrete-Time Models for Implicit port-Hamiltonian Systems grangian and Hamiltonian Methods for Nonlinear Control 2006, Francesco Bullo and KenjiFujimoto, eds., Springer, 2007, pp. 87 – 98.[16] John M. Lee , Introduction to Smooth Manifolds , Springer-Verlag, New York, 2003.[17]
Benedict Leimkuhler and Sebastian Reich , Simulating Hamiltonian Dynamics , CambrigeUniversity Press, Cambridge, UK, 2004.[18]
Wei Lin , Feedback stabilization of general nonlinear control systems: A passive system ap-proach , Systems and Control Lett., 25 (1995), pp. 41 – 52.[19]
Jerrold E. Marsden and Tudor S. Ratiu , Introduction to Mechanics and Symmetry ,Springer-Verlag, New York, 1999.[20]
Bernhard Maschke, Arjan J. van der Schaft, and Peter C. Breedveld , An intrinsicHamiltonian formulation of network dynamics: Non-standard Poisson structures and gy-rators , Journal of the Franklin Institute, 329 (1992), pp. 923 – 966.[21]
S. Monaco and Doroth´ee Normand-Cyrot , Advanced tools for nonlinear sampled-data sys-tems’ analysis and control , European Journal of Control, 13 (2007), pp. 221 – 241.[22]
Peter J. Olver , Applications of Lie groups to differential equations , Springer-Verlag, NewYork, 1993.[23]
Romeo Ortega, Arjan J. van der Schaft, Iven Mareels, and Bernhard Maschke , Putting energy back in control , IEEE Control Syst. Mag., (2001), pp. 18–33.[24]
Romeo Ortega, Arjan J. van der Schaft, Bernhard Maschke, and Gerardo Escobar , Interconnection and damping assignment passivity-based control of port-controlled Hamil-tonian systems , Automatica, 38 (2002), pp. 585 – 596.[25]
Sebastian Reich , On a geometrical interpretation of differential-algebraic equations , CircuitsSyst. Signal Process., 9 (1990), pp. 367 – 382.[26] ,
On an existence and uniqueness theory for nonlinear differential-algebraic equations ,Circuits Syst. Signal Process., 10 (1991), pp. 343 – 359.[27] ,
Symplectic integration of constrained Hamiltonian systems by composition methods ,SIAM J. Numer. Anal., 33 (1996), pp. 475 – 491.[28]
Werner C. Rheinboldt , Differential-algebraic systems as differential equations on manifolds ,Mathematics of computation, 43 (1984), pp. 473–482.[29] ,
On the existence and uniqueness of solutions of nonlinear semi-implicit differential-algebraic equations , Nonlinear Analysis. Theory, Methods and Applications, 16 (1991),pp. 647 – 661.[30]
J. C. Simo, N. Tarnow, and K. K. Wong , Exact energy-momentum conserving algorithmsand symplectic schemes for nonlinear dynamics , Int. J. Numerical Methods in Engineering,100 (1992), pp. 63 – 116.[31]
J. C. Simo and K. K. Wong , Unconditionally stable algorithms for rigid body dynamics thatexactly preserve energy and momentum , Int. J. Numerical Methods in Engineering, 31(1991), pp. 19 – 52.[32]
Arjan J. van der Schaft , L -Gain and Passivity Techniques in Nonlinear Control , Springer-Verlag, London, 2000.[33] Veeravalli S. Varadarajan , Lie Groups, Lie Algebras, and Their Representations , Springer-Verlag, New York, 1984.[34]