Extending Lagrangian and Hamiltonian Neural Networks with Differentiable Contact Models
AA Differentiable Contact Model to Extend Lagrangian and Hamiltonian NeuralNetworks for Modeling Hybrid Dynamics
Yaofeng Desmond Zhong Biswadip Dey Amit Chakraborty Abstract
The incorporation of appropriate inductive biasplays a critical role in learning dynamics fromdata. A growing body of work has been explor-ing ways to enforce energy conservation in thelearned dynamics by incorporating Lagrangian orHamiltonian dynamics into the design of the neu-ral network architecture. However, these existingapproaches are based on differential equations,which does not allow discontinuity in the states,and thereby limits the class of systems one canlearn. Real systems, such as legged robots androbotic manipulators, involve contacts and colli-sions, which introduce discontinuities in the states.In this paper, we introduce a differentiable con-tact model, which can capture contact mechanics,both frictionless and frictional, as well as bothelastic and inelastic. This model can also accom-modate inequality constraints, such as limits onthe joint angles. The proposed contact model ex-tends the scope of Lagrangian and Hamiltonianneural networks by allowing simultaneous learn-ing of contact properties and system properties.We demonstrate this framework on a series ofchallenging 2D and 3D physical systems withdifferent coefficients of restitution and friction.
In a large class of real-world problems, the underlyingphysical systems evolve in a piecewise-continuous manner.For example, we can play tennis since tennis balls collidewith the ground and rackets with high elasticity, but followsmooth trajectories governed by the laws of physics in be-tween those collisions. Our ability to walk and run dependsheavily on the frictional contacts between our shoes and theground; icy roads are slippery since the contact surface isalmost frictionless. Robotics systems, such as legged robotsand manipulators, also rely on contacts to perform tasks. Siemens Corporation, Princeton, New Jersey, USA. Corre-spondence to: Y. D. Zhong
These scenarios highlight the relevance and importance ofcontacts and collisions, which exist everywhere, albeit withdifferent properties.The last few years have witnessed a growing interest in learn-ing the dynamics governing motions of physical objects.Deep Lagrangian Networks (Lutter et al., 2019), Hamilto-nian Neural Networks (HNN) (Greydanus et al., 2019), Sym-plectic ODE-Net (SymODEN) (Zhong et al., 2020a), Sym-plectic Recurrent Neural Networks (SRNN) (Chen et al.,2020) and Lagrangian Neural Networks (LNN) (Cranmeret al., 2020) incorporate physics prior in the form of La-grangian/Hamiltonian dynamics to learn physical motionswith limited data and improved generalization. To learn thedynamics of energy-conserved systems, these physics pri-ors encode energy conservation into the architecture of theneural network. However, since the underlying dynamicsare represented by an ordinary differential equation (ODE),these previous works (except SRNN) cannot deal with sys-tems that involve collisions since collisions introduce a sud-den change in velocity, making the trajectories non-smooth.To capture the elastic collision of a billiard ball, SRNN man-ually reverses the momentum of the ball orthogonal to thecontact surface. However, this specialized technique cannotbe applied to frictional contact, inelastic contact, or objectswhich can rotate. The restriction on contact free systemsmakes it impossible to apply these energy-based approachesto physical systems with contacts.Another line of work, such as Neural Physics Engine(NPE) (Chang et al., 2016) and Interaction Networks (IN)(Battaglia et al., 2016), use relational priors to learn mo-tions of physical systems. These models are not ODE-basedand can handle collisions and contacts. A recent work byYang et al. (2020) directly learns physical constraints andcontact with an iterative neural projection (INP) algorithm.However, only frictionless 2D contacts have been testedon NPE, IN and INP. The learning of frictional contacts aswell as 3D contact is underexplored in the literature. In thiswork, we introduce a contact model, which can handle fric-tional contacts both with or without elasticity. The proposedcontact model solves two convex optimization problemsfor calculating the jump in velocity during contact. Weleverage recent progress on differentiating through convex a r X i v : . [ c s . R O ] F e b Differentiable Contact Model to Extend Lagrangian and Hamiltonian Neural Networks for Modeling Hybrid Dynamics optimization problems (Agrawal et al., 2019) make our con-tact model differentiable. We demonstrate the performanceof our contact model in learning coefficients of restitutionand friction of different 2D and 3D contacts.The three main contributions of our work are:(1) We propose a differentiable contact model, which cancapture elastic or inelastic, frictional or frictionless contactsand collisions as well as joint limits.(2) We demonstrate simultaneous learning of system andcontact properties of various physical systems by com-bining the proposed contact model and Constrained La-grangian/Hamiltonian neural networks. This way, we ex-tend the applicability of previous Lagrangian/Hamiltonian-inspired deep learning methods to more realistic systemswith contact and collisions.(3) The learned contact properties, i.e., coefficients of resti-tution and friction, are interpretable and matches the groundtruth with high accuracy.
Consider a rigid body system with holonomic constraints(equality constraints, see Sec. 2.3) but free of collisionand contact. The configuration of the system at time t can be described by a set of coordinates x ( t ) =( x ( t ) , x ( t ) , ..., x D ( t )) . Then the time evolution of therigid body system can be expressed as the following second-order ordinary differential equation (ODE), ¨ x = h ( x , ˙ x ; p s ) , (1)where p s denote system properties, which may includeinertia of objects M ( x ) and potential energy V ( x ) . Thevector-valued function h can be derived from the laws ofphysics, e.g, Lagrangian/Hamiltonian dynamics. This ODEcan also be written as the following first-order ODE (cid:18) ˙ x ˙ v (cid:19) = (cid:18) vh ( x , v ; p s ) (cid:19) = g ( x , v ; p s ) , (2)where v := ˙ x . There are two popular choices for the coor-dinates x – the generalized coordinates, and the Cartesiancoordinates. The generalized coordinates are usually chosenas a set of independent coordinates which implicitly en-forces holonomic constraints in the system. The Cartesiancoordinates are in general not independent of each other,so that the holonomic constraints in the system must beenforced explicitly in the dynamics (2). Although the actualform of g is usually derived from Lagrangian mechanics,here we demonstrate our results based on the Lagrangianand Hamiltonian formalism with Cartesian coordinates. Theproposed contact model is independent of the set of coordi-nates and form of g . We provide the form of g used in thiswork in the Supplementary Material. Algorithm 1:
Rigid Body Dynamics with Contact
Input : t , t , ..., t N Sequence of time points ( x , v ) Initial condition at t p s = ( M ( x ) , V ( x )) System properties p c = ( µµµ, e P ) Contact properties g ( x , v ; p s ) First-order system dynamics Initialize output trajectories T = { ( x , v ) } . for i = 0 → N − do ( x i +1 , v i +1 ) ← ODESolve ( g , ( x i , v i ) , t i , t i +1 ) Get active contacts c a (collision detection) if exists active contacts then ∆ v ← ContactModel ( x i +1 , v i +1 , c a , p s , p e ) v i +1 ← v i +1 + ∆ v T ← T ∪ { ( x i +1 , v i +1 ) } In robotics tasks, the above assumption of no collision andcontact no longer holds. For example, legged robots movearound through repeated collisions/contacts between therobot legs and the ground, and robot arms grasp objectsby making frictional contact with them. The difficulty ofmodeling these phenomena is that they essentially make thedynamics discontinuous. For example, when a ball hits theground, its velocity changes from pointing downward topointing upward in an infinitesimally small period of time,which can be modelled as an instantaneous change in ve-locity. In general, contacts, collisions and joint limits canall be modelled in this way. Algorithm 1 shows the generalprocedure of modelling rigid body with contacts, where ajump in velocity is calculated by a contact model when thereexists active contacts. From a simulation perspective, withknown system properties, contact properties (coefficients offriction and restitution) and vector field g , the trajectory ofthe system can be simulated by Algorithm 1. From a learn-ing perspective, we frame the problem as learning unknownsystem and contact properties from a given set of trajecto-ries given a model prior of vector field g . In this case, wecan parametrize the unknown system and contact proper-ties ( p s , p c ) by neural networks and learnable parameters,predict trajectories by Algorithm 1 and minimize the dif-ference between the predicted and true trajectories by backpropagation. This training scheme requires all operations inthe forward pass (Algorithm 1) to be differentiable. Thereare two key parts in the forward pass, the ODE solver mod-ule and the contact model. Operations in the ODE solverare in general differentiable and Neural ODE (Chen et al.,2018) provides a framework of back propagating throughODE solvers with constant memory cost. Neural ODE isleveraged by Symplectic ODE-Net (Zhong et al., 2020a)and Constrained Lagrangian/Hamiltonian Neural Network(CLNN/CHNN) (Finzi et al., 2020) to learn unknown sys- Differentiable Contact Model to Extend Lagrangian and Hamiltonian Neural Networks for Modeling Hybrid Dynamics tem properties in rigid body dynamics without contacts. Inthis work, we provide a differentiable contact model so thatwe can extend these previous works to learn rigid bodydynamics with contacts.
Holonomic constraints are equality constraints which canbe collected into a column vector Φ( x ) ∈ R E with equality Φ( x ) = . Differentiating this constraint w.r.t. time, wehave ˙Φ = ( D x Φ) ˙ x = ( D x Φ) v = J E ( x ) · v = , (3)where we denote the equality constraint Jacobian J E ( x ) := D x Φ ∈ R E × D . Eqn. (3) implies that holonomic constraintsrequire the velocity v to be always in the null space of equal-ity constraint Jacobian J E ( x ) . We will use this property toderive impulses caused by equality constraints. Contacts in general can be expressed as inequalities Φ C ( x ) ≥ . A ball bouncing on the ground, for exam-ple, requires the whole ball to be above the ground. Whenthe equality holds for a contact, we refer to the contact as anactive contact, otherwise, an inactive contact. If there existsactive contacts, contact impulses will cause an instantaneousvelocity change. In practice, the set of active contacts iscalculated by a collision detection (CD) module.A conceptual contact can contribute to one or more dimen-sions in the contact space, corresponding to one or moredimensions of contact impulse. Take Fig. 1 as an example.Mass 2 at the end of the pendulum would experience a con-tact impulse f C = ( f n , f t ) in the two dimensional contactspace - f n is the component normal to the contact surface,and f t is the friction impulse tangential to the contact sur-face. For 3D systems, the contact space is three dimensionalwith two tangential components. Assume that the contactspace for all active contacts is C dimensional, then we definecontact Jacobian J C ( x ) ∈ R C × D , which maps velocities v in the coordinate space to v C in the contact space, v C = J C ( x ) · v . (4) The elasticity of a collision can be captured by the coeffi-cient of restitution (COR). According to Newton’s hypoth-esis (Newton, 1999), COR is defined as the ratio of thenormal relative velocity after the collision to that before thecollision, ranging from 0 to 1. This definition of COR cancause unrealistic energy increases when the contact is fric-tional and the COR is close to 1 (Kane & Levinson, 1985;Stronge, 1991). Alternatively, Poisson (Poisson, 1817) di-vides the collision into two phases. The former, referred Space
Contact Space Impulses During Contact
Frictional Contact
Figure 1.
A ball collide with a pendulum. The equality constraintimpulse f E ensures that equality constraints are always satisfied. to as the compression phase, start with the first contact ofthe bodies and stops at the greatest compression. The lat-ter, referred to as the restitution phase, start right after thecompression phase till the separation of bodies. Accordingto Poisson’s hypothesis, the COR is defined as the ratio ofthe normal contact impulse in the restitution phase to thatin the compression phase. Poisson’s hypothesis is favoredin simulation since it will not lead to unrealistic energy in-crease. For a detailed comparison of different hypotheses,please refer to (Djerassi, 2009a;b). In this paper, we defineCOR e P in accordance with Poisson’s hypothesis. In this section, we introduce a differentiable contact modelfor learning rigid body dynamics with contacts. For sim-plicity, we present the model by referring to x and v asposition and velocity in Cartesian space, but the derivationis valid for any other choice of coordinate system. We listthe notations used in the paper in Supplementary Material. When there are active contacts, we construct the contactJacobian J C for active contacts. For brevity of notation,we drop explicit dependence on x from now onward. FromNewton’s second law, the change of momentum duringcontact equals the impulses, which can be described as Mv + = Mv − + J TC f C + J TE f E , (5)where v − and v + denote the Cartesian space velocity beforeand after the instantaneous velocity change, M is the inertiamatrix, J TC maps contact impulses in the contact space f C tocontact impulses in Cartesian space and J TE maps equalityconstraint impulses f E to equality constraint impulses inCartesian space. The impulses f C and f E should not beconfused with forces. An impulse is an integral of forceover time, which contributes to the change in momentum.The equality constraint impulses f E is caused by contactimpulse f C . See Fig. 1 for an intuitive example. Theirdependence can be revealed from the fact that the velocityin Cartesian space at any time is in the null space of J E (Sec.2.3.) We can left multiply the above equation by J E M − Differentiable Contact Model to Extend Lagrangian and Hamiltonian Neural Networks for Modeling Hybrid Dynamics and solve for f E , f E = − ( J E M − J TE ) − J E M − J TC f C . (6)Thus, from Eqn. (5) and (6), we can express instantaneousvelocity change as v + = v − + ˆM − J TC f C , (7)where ˆM − = M − − M − J TE ( J E M − J TE ) − J E M − . (8) ˆM can be interpreted as the inertia that incorporates equalityconstraints.In order to solve for contact impulses, we left multiply Eqn.(7) by J C to project the instantaneous velocity change intothe contact space: v + C = v − C + Af C , (9)where A = J C ˆM − J TC , which can be interpreted as the in-verse inertia in the contact space. Our contact model solvescontact impulses in two phases - the compression phase andthe restitution phase, both of which can be described by Eqn.(9). Here we express the instantaneous velocity change intwo phases as follows v c + C = v c − C + Af cC , (10) v r + C = v c + C − v ∗ C + Af rC , (11)where f cC and f rC are the contact impulses during the com-pression phase and the restitution phase and need to besolved by the contact model. The target velocity v ∗ C isincluded in the restitution phase to compensate existing pen-etration in the simulation. Please see Sec. 3.4 and Supple-mentary Material for details on compensating penetration. In this work, we consider two types of contact - frictionalcontact and limit constraint. For any conceptual frictionalcontact i in 3D contact space, the contact force f i ∈ R must lie in the friction cone, µ i f i,n ≥ (cid:113) f i,t + f i,t , ∀ i, (12)where µ i is the coefficient of friction for conceptual contact i . The normal forces must be non-negative, since objectscan only push but not pull others: f i,n ≥ , ∀ i. (13)For any limit constraint such as limit in joint angle or dis-tance, the contact space is essentially one dimensional andthe constraints on contact forces are only f i,n ≥ . Mathe-matically this is equivalent to setting up a 3D contact spacelike that in frictional contact and letting µ i = 0 in Eqn. (12). The idea behind solving contact impulses during compres-sion phase is the maximum dissipation principle (Stewart,2000), which says the compression impulses should maxi-mize the rate of energy dissipation. Equivalently, the com-pression impulses are those that minimizes the kinetic en-ergy at the end of the compression phase. This yields thefollowing optimization problem Minimize f cC , v c + C
12 ( v c + C ) T A − v c + C (14)subject to (10) , (12) , (13) . Substitute Eqn. (10) into the objective, we get the followingequivalent optimization problem, from which we can solvefor contact impulses in compression phase f cC ,Minimize f cC
12 ( f cC ) T Af cC + ( f cC ) T v c − C (15)subject to (12) , (13) . Similarly, we can set up an optimization problem to solvefor the contact force in restitution phase f rC . Accordingto Poisson’s hypothesis, the normal components in f rC arethose in f cC scaled by coefficient of restitution, so that weset up the following constraint: f ri,n ≥ e P,i · f ci,n , ∀ i (16)We have inequality instead of equality here since we wouldlike to compensate existing penetration in the simulation.Since we simulate the rigid body system in discrete timesteps, almost every time when a collision is detected, pene-tration has occurred between objects involved in that colli-sion. If the collision is perfectly inelastic, i.e., COR e P = 0 ,then normal impulse during restitution phase would be zero,which fails to fix existing unrealistic penetration. By settingup the constraint as in Eqn. (16), a larger normal impulseis allowed to fix existing penetration. Thus, with Eqn. (11)and (16), we take penetration compensation into accountand set up the following optimization problem from theprinciple of maximum dissipation to solve for f rC .Minimize f rC
12 ( f rC ) T Af rC + ( f rC ) T ( v c + C − v ∗ C ) (17)subject to (12) , (13) , (16) . From the principle of maximum dissipation, the equalityin Eqn. (16) would hold for the solution when COR islarge and penetration is small, thus respecting Poisson’shypothesis when penetration can be fixed naturally. Strictly speaking, this form is not correct because A is invert-ible only if there exists no equality constraint in the system. Whenequality constraints do exist, a pseudo-inverse of A should be usedhere, and the form (15) can still be derived. Differentiable Contact Model to Extend Lagrangian and Hamiltonian Neural Networks for Modeling Hybrid Dynamics (a) (b) (c) (d) (e)
Figure 2.
Simulated systems with contact. From left to right: (a) bouncing point masses, (b) bouncing disks, (c) chained pendulums withground, (d) gyroscope with a wall, and (e) rope. The black lines in bouncing disks show the orientations of disks.
Solving optimization problems (15) and (17) for contact im-pulses allow us to calculate instantaneous velocity changeto perform simulation (Algorithm 1). Moreover, we wouldlike to back-propogating through the contact model to learnunknown properties. In fact, our proposed contact modelis differentiable, thanks to recent progress on differentiat-ing through convex optimization problems. Both problems(15) and (17) are convex optimization problems with con-vex quadratic objectives as well as linear constraints andsecond order cone constraints. (We proof that A is positivesemi-definite in the supplementary materials.) We can thenexpress our problems using disciplined parametrized pro-gramming and set up these two problems as differentiablelayers using CvxpyLayers (Agrawal et al., 2019). Thus, ourmodel is differentiable and can be used in learning tasks. We simulate five different systems with contacts (Fig. 2),and propose 8 tasks based on these systems with differentcontact properties to test our contact model (Table 1.)
Bouncing point masses.
This system is often referred toas bouncing balls in previous works (Battaglia et al., 2016;Chang et al., 2016). We call it bouncing point masses in-stead since each object is essentially a circle with a pointmass at the center and cannot rotate like real balls. For n ob-jects bouncing in the box, there exists n ( n − / possiblecontacts between objects and n possible contacts betweenobjects and walls. These contacts cannot be all active simul-taneously. We set up two tasks of 5 bouncing point masseswith different configurations, which will be referred to asBP5-e and BP5. BP5-e is a homogeneous setting, where themasses and radii are the same for all objects, and contactproperties ( e P = 1 and µ = 0 ) are the same for all contacts.This task conserves energy since no energy is lost duringthe collisions and during collision-free periods. BP5 is aheterogeneous setting where the masses and radii are differ-ent for different points and contact properties are differentfor different contacts. Bouncing disks.
A real 2D disk has mass spread over thecircle and thus can rotate, especially when frictional contactsare involved. Thus we extend the bouncing point massessystem to bouncing disks, where all the disks can rotate. Weuse the extended bodies representation introduced in (Finziet al., 2020) to embed the motion of disks in Cartesian coor-dinates. The idea is to use the motion of 3 points - the centerof mass as well as the unit vectors aligned with two principleaxes - to describe the motion of a disk. Since the relativeposition of these 3 points are fixed, this representation willintroduce 3 equality constraints for each disk. A contactimpulse will be distributed properly into these 3 points ina way that obeys the law of physics. Please refer to (Finziet al., 2020) for more details on this representation. We sim-ulate 5 heterogeneous bouncing disks with heterogeneouscontact properties. This task is referred to as BD5.
Chained pendulums with ground.
The 2-pendulum col-liding with ground has been used to study and analyze con-tact models more than three decades ago (Kane & Levin-son, 1985). Until recently, some works (Finzi et al., 2020;Zhong et al., 2020b) have studied learning dynamics ofN-pendulums without contacts. Here we simulate a 3-pendulum system above the ground where the lowest pen-dulum can collide with the ground. The masses are locatedat the joints and the sizes of the joints are different. Wefollow the convention to assume that pendulums cannot col-lide with each other. We propose two tasks: CP3-e with e P = 1 and µ = 0 , where energy is conserved, and CP3with e P = 0 and µ = 0 . . Gyroscope with a wall.
Gyroscope is a 3D system thatexhibits complex dynamics such as precession and nutation.In order to test our contact model in 3D space, we extend thegyroscope system by putting a wall near it so that collisionscan happen. The motion of the gyroscope is embedded inCartesian coordinates using the extended bodies represen-tation(Finzi et al., 2020). This representation introduces6 equality constraints. As the gyroscope is attached to aball joint, one more equality constraint is introduced. Wepropose two tasks: Gyro-e with e P = 1 and µ = 0 , whichconserves energy, and Gyro with e P = 0 . and µ = 0 . . Rope.
Our contact model can also capture limits in joint
Differentiable Contact Model to Extend Lagrangian and Hamiltonian Neural Networks for Modeling Hybrid Dynamics
Table 1.
Benchmark tasks. The columns D , E , max( C ) denote dimension of the dynamics, number of equality constraints, and themaximum number of contacts that could be simultaneously active , respectively.Name System D E max( C ) Space Same e P , µ forall contacts Linearpotential ConserveenergyBP5-e Bouncing point masses
10 0 ∼
2D Y Y YBP5 Bouncing point masses
10 0 ∼
2D N Y NCP3-e Chained pendulums with ground
2D Y Y YCP3 Chained pendulums with ground
2D Y Y NBD5 Bouncing disks
30 15 ∼
2D N Y NRope Rope
20 0 ∼
2D Y N NGyro-e Gyroscope with a wall
12 7 1
3D Y Y YGyro Gyroscope with a wall
12 7 1
3D Y Y N
BP5-e BP5 CP3-e CP3 BD5 Rope Gyro-e Gyro −5 −3 −1 T r a j . r e l a t i v e e rr o r CM-CD-CLNN CM-CD-CHNN CMr-CD-CLNN CMr-CD-CHNN MLP-CD-CLNN IN-CP-CLNN IN-CP-SP
Figure 3.
Trajectory relative error (log scale). Each error is averaged over 100 test trajectories of length 5. angles and distances, which we show in this rope system.The motion of the rope is described by 10 equally spacedpoints along the rope. The distance between adjacent pointsare not fixed as in the chained pendulums system. Instead,the rope can be stretched. The stretch is modelled by elasticsprings connecting each pair of adjacent points. We set themaximum stretch and minimum stretch to be 1.2 and 0.8,respectively, which implies that two adjacent points are notallowed to be pushed or pulled by more than of theirdistance at rest. We also assume that two adjacent segmentscannot be bent over a predefined angle (0.2rad). The abovestretch and bending constraint can be handled by our contactmodel with e = 0 and µ = 0 . During simulation of the rope,a total of 19 “contacts" can be active at the same time, whichmakes it nontrivial to solve for contact impulses. The forceof the spring is modelled via the potential energy, whichresults in a potential energy function that is not linear inthe location of points. This is the only system tested in thiswork that has a nonlinear potential energy function. Thissetup is similar to the rope proposed in (Yang et al., 2020),and differs from string proposed in (Battaglia et al., 2016),as the latter impose no bending constraint. In each task, we jointly learn system and contact propertiesfrom trajectory data by extending CLNN/CHNN with theproposed contact model.
Data:
For each task, the training set is generated by ran- domly sampling 800 collision-free initial conditions andthen simulating the dynamics for 100 times step. Since forsome systems, there are very few data points in a trajectorythat involves collision, we select a small chunk containing 5consecutive time steps from each simulated trajectory suchthat the final training set contains 800 trajectories of length5, where half of the trajectories contain collisions and theother half are collision-free. We also make sure that theinitial state of these selected chunks is collision-free. Thetest set is generated in a similar way with 100 trajectories.
Architecture and Training details:
In the experiments,we assume the system properties, i.e., object inertia andpotential energy, as well as contact properties, i.e., coeffi-cients of friction and restitution, are unknown and need tobe learned from trajectory data. The system properties areparametrized as in CLNN and CHNN (Finzi et al., 2020).As for contact properties, each coefficient of friction is posi-tive, so it is parametrized by a learnable parameter passedthrough a
ReLu function. Each coefficient of restitution isin the interval of [0 , , so it is parametrized by a learnableparameter passed through a hard sigmoid function. Thepredicted trajectories are generated by running Algorithm1 with parametrized system and contact properties. We useRK4 as the ODE solver in Neural ODE. We compute the L loss of the difference between predicted and true trajectories.The gradients are computed by differentiating through Al-gorithm 1, and neural network parameters are updated usingthe AdamW optimizer (Kingma & Ba, 2014; Loshchilov &Hutter, 2019) with a learning rate of 0.001. Differentiable Contact Model to Extend Lagrangian and Hamiltonian Neural Networks for Modeling Hybrid Dynamics −5 −2 BP5-e
CP3-e
BD5
Gyro-e
Time (s) −5 −2 BP5
Time (s)CP3
Time (s)Rope
Time (s)Gyro R e l a t i v e e rr o r i n s t a t e s CM-CD-CLNN CM-CD-CHNN CMr-CD-CLNN CMr-CD-CHNN MLP-CD-CLNN IN-CP-CLNN IN-CP-SP
Figure 4.
Relative error (log scale) along long test trajectories (50 times steps). Each curve is averaged over 100 test trajectories. Verticaldashed lines show the trajectory length during training.
Models:
We implement two slightly different versions ofthe contact model. The first version, referred to as CM, setup optimization problems exactly as stated in (15) and (17).The second version, referred to as CMr, adds a diagonalpositive regularization matrix R to A in (15) and (17), suchthat ( A + R ) is always positive definite, which ensuresa unique global minimum in each problem. These twoversions are combined with CLNN and CHNN to set up thefollowing four neural network models: ( i ) CM-CD-CLNN,( ii ) CM-CD-CHNN, ( iii ) CMr-CD-CLNN, and ( iv ) CMr-CD-CHNN. The “CD" in model names emphasizes that weassume that a collision detection module is given. Baselines:
We also set up three baselines. In the first base-line, MLP-CD-CLNN, we calculate the instantaneous ve-locity change from a multi-layer perceptron (MLP) insteadof our contact model. Our second baseline, IN-CP-CLNN,calculate velocity change from an interaction network (IN)(Battaglia et al., 2016) without collision detection, since INhas the ability to learn collisions and contact. IN requirestrue system and contact properties as input. Here the “CP"in the model name emphasizes true contact properties aregiven and the system properties learned by CLNN are fedinto IN. Our last baseline, IN-CP-SP, is an original interac-tion network which has shown strong ability in predicting2D rigid body trajectories without equality constraints, buthaven’t been tested on systems with equality constraints or3D systems. The name emphasizes that true system and con-tact properties are given. For baseline models, we transformeach trajectory into multiple one-step pairs, as done in IN.
Prediction:
We report the average relative L error over thetest trajectories of 7 models on 8 tasks in Fig. 3. In all 2Dtasks, our models beat baseline models. The performancedifference between CLNN and CHNN is minor since their architectures are similar. In most tasks, CM outperformsCMr. In the BP5-E task, CM beats CMr by 2 orders ofmagnitude. In two 3D gyroscope tasks, CM performs badly,while CMr performs the best. The reason is discussed be-low. IN does not perform well even in BP5 tasks since ourtraining set (3.2k one-step pairs) is much smaller than thedataset (1M one-step pairs) used in the IN paper. We alsoreport average relative L errors along test trajectories of 50time steps in Fig. 4 in order to show each model’s abilityin long term prediction. We observe that our contact modelCM and CMr outperform baselines in all 2D tasks and CMroutperforms baselines in 3D tasks. Interpretability:
Table 2 shows the learned contact proper-ties by our 4 models in 6 tasks where the contact propertiesare the same for all contacts. For 4 2D tasks, CM can learncontact properties accurately which explains its good perfor-mance in prediction. CMr is an approximate model and doesnot infer contact properties as accurately as CM. For 3Dtasks, CM fails to infer contact properties, which explainsthe large prediction error in Figure 3. This can probablybe attributed to the fact that the positive semi-definitenessof matrix A can introduce some errors in gradient compu-tation. By using a positive semi-definite matrix ( A + R ) ,CMr avoids this problem and infers contact properties thatare close to the ground truth. Energy:
The prior of Lagrangian/Hamiltonian dynamicsconserve energy along each trajectory, which is one of thereason that Lagrangian/Hamiltonian-based neural networkmodels perform better in prediction and generalization. Anelastic frictionless contact ( e P = 1 , µ = 0 ) also conservesenergy. Fig. 5 illustrates how the total energy changes overtime for the predicted trajectories of 7 models on three tasksthat conserves energy. Models using CM conserves energyin BP5-e and CP3-e tasks since CM learns contact prop-erties perfectly (Table 2). This shows that the proposedcontact model can uncover the energy conserving aspect Differentiable Contact Model to Extend Lagrangian and Hamiltonian Neural Networks for Modeling Hybrid Dynamics
Table 2.
Learned contact properties from our 4 neural network models on 6 tasks that has unique contact properties for all contacts. Boldnumbers are the best learned contact properties in each task across 4 models.BP5-e CP3-e CP3 Rope Gyro-e Gyro µ e P µ e P µ e P µ e P µ e P µ e P Ground Truth 0.000 1.000 0.000 1.000 0.500 0.000 0.000 0.000 0.000 1.000 0.100 0.800CM-CD-CLNN
CMr-CD-CHNN
Time (s)
BP5-e
Time (s) −1.9−1.8−1.7−1.6
CP3-e
Time (s)
Gyro-e E n e r g y a l o n g t r a j . True CM-CD-CLNN CM-CD-CHNN CMr-CD-CLNN CMr-CD-CHNN MLP-CD-CLNN IN-CP-CLNN IN-CP-SP
Figure 5.
Energy of the predicted trajectories of all 7 models on a sampled test initial condition from BP5-e, CP3-e and Gyro-e tasks. Thetrue energy in each task is represented by the horizontal black line in the middle, which is constant along the trajecotry. even though energy conservation has not been not enforcedexplicitly. Models using CMr loses energy each time there isa collision in CP3-e and Gyro-e, since they learn nonzero co-efficients of friction (Table 2). The baseline models performthe worst in terms of energy conservation.
Contact Model:
Our contact model shares similarity withthe contact model of MuJoCo (Todorov, 2011; 2014). Thereare three differences: (1) our contact model handles elasticcontacts while MuJoCo only focuses on inelastic contacts;(2) MuJoCo solves the convex optimization problem witha generalization of Projected Gauss-Seidel method whilewe leverage the open source scs solver (O’donoghue et al.,2016) to solve the optimization problem; (3) the dynamics inMuJoCo is described using generalized coordinates, whilewe use Cartesian coordinates. This is because previouswork has shown it is easier to learn system properties underCartesian coordinates (Finzi et al., 2020). Another categoryof contact models solve contact impulses by solving a linearcomplementarity problem (LCP) (Anitescu & Potra, 1997).Recently, a number of works (de Avila Belbute-Peres et al.,2018; Degrave et al., 2019; Liang & Lin, 2020) has proposeddifferentiable LCPs to solve contact impulses. However,none of these works attempted to learn contact properties.
Neural ODEs with Jumps:
Neural ODE solvers and theirvariants (Chen et al., 2018; Dupont et al., 2019; Massaroliet al., 2020) have emerged as effective tools for learning dynamics from data. Although these solvers and their ap-plications in learning physics-based models primarily focuson systems with continuous trajectories, a small body ofwork has been investigating ways to learn dynamic modelsgoverning discontinuous state trajectories. Jia & Benson(2019) use a stochastic process to introduce discontinuitiesfor learning the dynamics governing piecewise-continuousstate trajectories. Gwak et al. (2020) leverage a separateneural network to learn discrete, deterministic jumps thatcan happen only at pre-specified time-steps. Herrera et al.(2021) uses a Neural ODE with jumps to model the condi-tional expectation of a stochastic process given its previousobservations.
In this work, we have introduced a differentiable contactmodel, which can capture contact dynamics with differentproperties. Our contact model extends CHNN/CLNN toenable learning of hybrid dynamics in rigid body systemsand offer interpretability about system and contact prop-erties. Future works will incorporate model-based controland explore learning interpretable safe control policies forrobotics applications.
References
Agrawal, A., Amos, B., Barratt, S., Boyd, S., Diamond,S., and Kolter, J. Z. Differentiable convex optimizationlayers. In
Advances in Neural Information Processing
Differentiable Contact Model to Extend Lagrangian and Hamiltonian Neural Networks for Modeling Hybrid Dynamics
Systems , pp. 9562–9574, 2019.Anitescu, M. and Potra, F. A. Formulating dynamic multi-rigid-body contact problems with friction as solvablelinear complementarity problems.
Nonlinear Dynamics ,14(3):231–247, 1997.Battaglia, P., Pascanu, R., Lai, M., Jimenez Rezende, D., andkavukcuoglu, k. Interaction networks for learning aboutobjects, relations and physics. In
Advances in NeuralInformation Processing Systems , pp. 4502–4510, 2016.Chang, M. B., Ullman, T., Torralba, A., and Tenenbaum,J. B. A compositional object-based approach to learningphysical dynamics. arXiv preprint arXiv:1612.00341 ,2016.Chen, R. T., Rubanova, Y., Bettencourt, J., and Duvenaud,D. K. Neural ordinary differential equations. In
Advancesin Neural Information Processing Systems , volume 31,pp. 6571–6583, 2018.Chen, Z., Zhang, J., Arjovsky, M., and Bottou, L. Symplec-tic recurrent neural networks. In
International Conferenceon Learning Representations , 2020.Cranmer, M., Greydanus, S., Hoyer, S., Battaglia, P.,Spergel, D., and Ho, S. Lagrangian neural networks.In
ICLR 2020 Workshop on Integration of Deep NeuralModels and Differential Equations , 2020.de Avila Belbute-Peres, F., Smith, K., Allen, K., Tenenbaum,J., and Kolter, J. Z. End-to-end differentiable physics forlearning and control.
Advances in neural informationprocessing systems , 31:7178–7189, 2018.Degrave, J., Hermans, M., Dambre, J., et al. A differentiablephysics engine for deep learning in robotics.
Frontiers inneurorobotics , 13:6, 2019.Djerassi, S. Collision with friction; part a: Newton’s hy-pothesis.
Multibody System Dynamics , 21(1):37, 2009a.Djerassi, S. Collision with friction; part b: Poisson’s andstornge’s hypotheses.
Multibody System Dynamics , 21(1):55, 2009b.Dupont, E., Doucet, A., and Teh, Y. W. Augmented neuralodes. In
Advances in Neural Information ProcessingSystems , volume 32, pp. 3140–3150, 2019.Engin, M., Wang, L., Zhou, L., and Liu, X. Deepkspd:Learning kernel-matrix-based spd representation for fine-grained image recognition. In
Proceedings of the Euro-pean Conference on Computer Vision (ECCV) , pp. 612–627, 2018. Finzi, M., Wang, K. A., and Wilson, A. G. Simplifyinghamiltonian and lagrangian neural networks via explicitconstraints.
Advances in Neural Information ProcessingSystems , 33, 2020.Greydanus, S., Dzamba, M., and Yosinski, J. Hamiltonianneural networks. In
Advances in Neural InformationProcessing Systems , pp. 15379–15389, 2019.Gwak, D., Sim, G., Poli, M., Massaroli, S., Choo, J., andChoi, E. Neural ordinary differential equations for inter-vention modeling, 2020.Herrera, C., Krach, F., and Teichmann, J. Neural jump ordi-nary differential equations. In
International Conferenceon Learning Representations , 2021.Jia, J. and Benson, A. R. Neural Jump Stochastic Differ-ential Equations. In
Advances in Neural InformationProcessing Systems , volume 32, pp. 9847–9858, 2019.Kane, T. R. and Levinson, D. A.
Dynamics, theory andapplications . McGraw Hill, 1985.Kingma, D. P. and Ba, J. Adam: A method for stochasticoptimization. arXiv preprint arXiv:1412.6980 , 2014.Liang, J. and Lin, M. C. Differentiable Physics Simulation.In
ICLR 2020 Workshop on Integration of Deep NeuralModels and Differential Equations , 2020.Loshchilov, I. and Hutter, F. Decoupled weight decay reg-ularization. In
International Conference on LearningRepresentations , 2019.Lutter, M., Ritter, C., and Peters, J. Deep lagrangian net-works: Using physics as model prior for deep learning. In
International Conference on Learning Representations ,2019.Massaroli, S., Poli, M., Park, J., Yamashita, A., and Asama,H. Dissecting Neural ODEs. In
Advances in NeuralInformation Processing Systems , volume 33, 2020.Newton, I.
The Principia: mathematical principles of natu-ral philosophy . Univ of California Press, 1999.O’donoghue, B., Chu, E., Parikh, N., and Boyd, S. Conicoptimization via operator splitting and homogeneous self-dual embedding.
Journal of Optimization Theory andApplications , 169(3):1042–1068, 2016.Poisson, S. Mechanics, vol. ii.
Trans. HH Harte, Longman,London , 1817.Stewart, D. E. Rigid-body dynamics with friction and im-pact.
SIAM review , 42(1):3–39, 2000.Stronge, W. J. Friction in collisions: Resolution of a paradox.
Journal of Applied Physics , 69(2):610–612, 1991.
Differentiable Contact Model to Extend Lagrangian and Hamiltonian Neural Networks for Modeling Hybrid Dynamics
Todorov, E. A convex, smooth and invertible contact modelfor trajectory optimization. In , pp. 1071–1076.IEEE, 2011.Todorov, E. Convex and analytically-invertible dynamicswith contacts and constraints: Theory and implementationin mujoco. In , pp. 6054–6061. IEEE,2014.Walter, S. F. and Lehmann, L. Algorithmic differentiationof linear algebra functions with application in optimumexperimental design (extended version). arXiv preprintarXiv:1001.1654 , 2010.Yang, S., He, X., and Zhu, B. Learning physical constraintswith neural projections. In
Advances in Neural Informa-tion Processing Systems , 2020.Zhong, Y. D., Dey, B., and Chakraborty, A. Symplecticode-net: Learning Hamiltonian dynamics with control. In
International Conference on Learning Representations ,2020a.Zhong, Y. D., Dey, B., and Chakraborty, A. Benchmarkingenergy-conserving neural networks for learning dynamicsfrom data. arXiv preprint arXiv:2012.02334 , 2020b.
Differentiable Contact Model to Extend Lagrangian and Hamiltonian NeuralNetworks for Modeling Hybrid Dynamics (Supplementary Material)
A Notation M ( x ) Inertia matrix V ( x ) Potential energy µµµ coefficients of friction e P coefficients of restitution p s = ( M ( x ) , V ( x )) System properties p c = ( µµµ, e P ) Contact properties f C Contact Impulses f E Equality Constraint Impulses J C ( x ) Contact Jacobian J E ( x ) Equality Constraint Jacobian v − / v + velocities (in Cartesian space) before/after a general impulse v − C / v + C velocities (in contact space) before/after a general impulse v c − C / v c + C velocities (in contact space) before/after the compression phase v r + C velocities (in contact space) after the restitution phase B Functional Form of the Constrained Lagrangian and Hamiltonian Dynamics
In this section, we present the functional form of system dynamics that we use in the experiments. Instead of usinggeneralized coordinates, we use Cartesian coordinates. This is because the inertia matrix under Cartesian coordinates areconstant and independent of the coordinates, which makes the learning of inertia easier, as pointed out in Finzi et al. (2020).
B.1 Constrained Lagrangian Dynamics
The first-order dynamics can be obtained from Finzi et al. (2020), which is (cid:18) ˙ x ˙ v (cid:19) = g ( x , v ; p s ) = (cid:18) vM − J TE [ J E M − J TE ] − [ J E M − ∇ x V − ( D x ( J E · v )) · v ] − M − ∇ x V (cid:19) . (S.1) B.2 Constrained Hamiltonian Dynamics
The Hamiltonian dynamics deal with position x ∈ R D and momentum p x = Mx instead of ( x , v ) . The derivation is not asstraightforward as in the Lagrangian case. We denote z = ( x , p x ) . The Hamiltonian equals the total energy of the systemand can be written as H ( x , p x ) = 12 p T x M − p x + V ( x ) , (S.2)For the E holonomic constraints Φ( x ) = , we can get another E constraints on position and momentum, i.e., ˙Φ( x , p x ) = 0 ,and collect these E constraints in a vector Ψ( z ) = (Φ , ˙Φ) . Then the Hamiltonian dynamics in z can be written as thefollowing differential equations ˙ z = J ∇ z H − J ( D z Ψ) T [( D z Ψ) J ( D z Ψ) T ] − ( D z Ψ) J ∇ z H, (S.3)where J is a symplectic matrix J = (cid:20) D − I D (cid:21) , (S.4) Differentiable Contact Model to Extend Lagrangian and Hamiltonian Neural Networks for Modeling Hybrid Dynamics and I D is the D × D identity matrix. In order to convert the ODE into a set of ODE in ( x , v ) , we introduce the matrix ˜ M − = (cid:20) I D
00 M − (cid:21) , (S.5)then we obtain the first order ODE (cid:18) ˙ x ˙ v (cid:19) = (cid:20) I D
00 M − (cid:21) (cid:18) ˙ x ˙ p x (cid:19) = ˜ M − J ∇ z H − ˜ M − J ( D z Ψ) T [( D z Ψ) J ( D z Ψ) T ] − ( D z Ψ) J ∇ z H (S.6) C Proof of Positive Semi-definiteness of A By definition, we have A = J C ˆM − J TC ∈ R C × C , where ˆM − = M − − M − J TE ( J E M − J TE ) − J E M − . (S.7)For any real physical system, the inertia matrix M is symmetric and positive definite. Thus, its inverse exists and can bedecomposed using Cholesky decomposition M = LL T . We can then express the inverse inertia that incorporates equalityconstraints as ˆM − = L ( I − P ) L T , where P is a projection matrix P = L T J TE ( J E M − J TE ) − J E L , (S.8)which satisfies P = P . A property of projection matrices is that the eigenvalues can only take two values: 1 or 0. Byeigen-decomposition, P and I − P can be written as P = (cid:0) V V (cid:1) (cid:18) (cid:19) (cid:18) V T V T (cid:19) = V V T , (S.9) I − P = (cid:0) V V (cid:1) (cid:18) I 00 I (cid:19) (cid:18) V T V T (cid:19) − (cid:0) V V (cid:1) (cid:18) (cid:19) (cid:18) V T V T (cid:19) = V V T , (S.10)where V ∈ R D × ( D − E ) . So we can decompose A into A = A Td A d , where A d = V T L T J TC . Then for any vector c ∈ R C ,we have c T Ac = ( A d c ) T A d c ≥ , (S.11)which proves that A is positive semi-definite. If the system does not have equality constraints, A has full rank and can bedecomposed using Cholesky decomposition. D Solving Contact Impulses using CvxpyLaers
In this section, we show an implementation of setting up the differentiable optimization problem in the compression phaseusing CvxpyLayers and PyTorch. We then show how we use this implementation in CM and CMr. import torch import cvxpy as cp from cvxpylayers.torch import CvxpyLayer def solve_compression_impulse( A_d: torch.Tensor, v_: torch.Tensor, mu: torch.Tensor, n_cld: int, d: int, ): C = v_.shape[0] f = cp.Variable((C, 1)) A_d_p = cp.Parameter(A_d.shape) v_p = cp.Parameter((C, 1)) mu_p = cp.Parameter((mu.shape[0], 1)) Differentiable Contact Model to Extend Lagrangian and Hamiltonian Neural Networks for Modeling Hybrid Dynamics objective = cp.Minimize(0.5 * cp.sum_squares(A_d_p @ f) + cp.sum(cp.multiply(f, v_p))) constraints = [cp.SOC(cp.multiply(mu_p[i], f[i*d]), f[i*d+1:i*d+d]) for i in range(n_cld)] \ + [f[i*d] >= 0 for i in range(n_cld)] problem = cp.Problem(objective, constraints) cvxpylayer = CvxpyLayer(problem, parameters=[A_d_p, v_p, mu_p], variables=[f]) impulse, = cvxpylayer(A_d, v_, mu) return impulse Listing 1.
Implementation of solving compression phase impulse using CvxpyLayers A_d impulse = solve_compression_impulse(A_d, v_, mu, n_cld, d) Listing 2. pseudocode of solving compression phase impulse in CM ... R = torch.eye(A.shape[0]).type_as(A)*1e-2 A_d = torch.cholesky(A+R, upper=True) impulse = solve_compression_impulse(A_d, v_, mu, n_cld, d) Listing 3. pseudocode of solving compression phase impulse in CMr
The only difference between CM and CMr is how we construct the matrix A d . In CM, A d is constructed as described inSection C, while in CMr, A d is constructed by adding a regularization and performing Cholesky decomposition. E Differentiable Symeig operation with non-distinct eigenvalues
The implementation of CM requires an operation of eigen decomposition of symmetric matrices. The operation shouldbe differentiable so that back-propagation can be performed. The current implementation of torch.symeig operation isbased on the general Daleckii-Krein formula (Walter & Lehmann, 2010; Engin et al., 2018) which assumes the eigenvaluesto be distinct. Therefore it does not support backward gradient calculation with non-distinct eigenvalues. If non-distincteigenvalues exists, some gradient values will be NaN. However, in our problem, the matrix P has repeated eigenvalues, andhence we cannot use the default torch.symeig to calculate gradients. Here, we derive a simple alternative backward gradientcalculation of this operation, which supports non-distinct eigenvalues.For a symmetric matrix X ∈ R n × n , the operation can be written as e , V = symeig( X ) , (S.12)where e ∈ R n is the collection of the eigenvalues and V ∈ R n × n is the collection of eigenvectors. In the backward pass,we need to calculate d X , the derivative of the input matrix, from d e and d V , the derivative of the output eigenvalues andeigenvectors. By definition, we have the eigen decomposition X = V ΣΣΣ V T , (S.13)where Σ = diag( e ) . Taking the derivative of this eigen decomposition, we have d X = V d ΣΣΣ V T + d V ΣΣΣ V T + V ΣΣΣ d V T (S.14) = V diag( d e ) V T + d V ΣΣΣ V T + V ΣΣΣ d V T (S.15)The first term is the contribution from the variation of eigenvalues and the latter two terms are the contribution from thevariation of the eigenvectors. It is easy to check that d X is also a symmetric matrix and the above gradient calculation isvalid for non-distinct eigenvalues. Differentiable Contact Model to Extend Lagrangian and Hamiltonian Neural Networks for Modeling Hybrid Dynamics
F Penetration Compensation
To compensate for an existing penetration during restitution phase in the simulation, we use the target velocity v ∗ C ∈ R C and the optimization problem (17). In this section, we discuss how to calculate the target velocity v ∗ C so that it does notviolate the equality constraints of the system. The calculation of v ∗ C might be nontrivial; however, in the backward pass,the gradients of v ∗ C are not required for learning contact properties. Thus, in practice, we do not calculate the backwardgradients for every calculation introduced in this section.To choose the target velocity v ∗ C , we first come up with a desired velocity v dC ∈ R C . For each direction normal to contactsurfaces, the component in v dC is calculated as the depth of penetration divided by integration time interval. For eachtangential dimension, the component in v dC is set to zero. This choice of v dC will fix penetration in the next time step. Thedownside is that for totally inelastic contacts, in the next few time steps, the bodies in collision might separate (because therelative velocity normal to the contact surface is greater than zero), which make the contact looks like partially elastic. Thisphenomenon can be avoided by using more than one time step to compensate the penetration, i.e., by setting the componentsin v dC to be a fraction of the depth of the penetration, as shown in the figure below. The reason that we cannot use v dC f i x i n t i m e s t e p f i x i n t i m e s t e p t=0 N o c o m p e n s a t i o n t=0.96 t=1.00 t=1.04 t=1.08 t=1.12 t=1.16 Figure 6.
Different v dC for compensation in a bouncing point mass with gravity and COR=0. First row : fixing penetration in 1 time step.The circle bounces off the ground after touching the ground.
Second row : fixing penetration in 4 time steps. The penetration is fixed andthe point mass doesn’t bounce up.
Third row : no penetration compensation. The penetration are not fixed over time. as the target velocity is that v dC might violate equality constraints. Take a ball collide with a pendulum as an example(Fig. 1 in main paper), v dC would violate the equality constraint of the pendulum, i.e., the velocity of object 2 can only beperpendicular to the pendulum. Thus, we need to transform v dC into target velocity v ∗ C that satisfies the equality constraints.The idea of obtaining a target velocity v ∗ C is to project the desired velocity v dC into Cartesian space, make corrections tosatisfy equality constraints and then project it back to the contact space. In Section 3.1, we showed how to project impulsesbetween contact space and Cartesian space. However, strictly speaking, the projection defined would introduce a scalingif one, say, project a vector from contact space to Cartesian space and back to contact space. This is not a problem forsolving contact impulses, but it will be problematic if we have this scaling in calculating target velocity v ∗ C . To fix this issue,we need to introduce the pseudoinverse of J C . Let’s assume we have some form of pseudoinverse J + C (we will define itlater.) Then the target velocity in Cartesian space is the sum of the desired velocity v dC projected into Cartesian space and acorrection term. v ∗ = J + C v dC + J TE v dE (S.16)The equality constraints require J E v ∗ = , from which we can solve for v dE = − ( J E J TE ) − J E J + C v dC , then we project thetarget velocity in Cartesian space v ∗ into contact space and get v ∗ C = J C v ∗ = J C ( I − J C ( J E J TE ) − J E ) J + C v dC (S.17)The form of the pseudoinverse J + C we use here is dependent on the shape of J C . We define the pseudoinverse of J C as J + C = (cid:40) J TC ( J C J TC ) − , if C ≤ D ( J TC J C ) − J TC , if C > D (S.18)
Differentiable Contact Model to Extend Lagrangian and Hamiltonian Neural Networks for Modeling Hybrid Dynamics
When C ≤ D , the dimension of contact space is smaller than that of Cartesian space, we can verify that projecting a velocityfrom contact space to Cartesian space and back to contact space equals the original velocity, i.e., J C J + C = I C . When C > D , the dimension of contact space is greater than that of Cartesian space, we can verify that projecting a velocity fromCartesian space to contact space and back to Cartesian space equals the original velocity, i.e., J + C J C = I D . G Sequences of Prediction Trajectories
In this section, we show prediction sequences of a sample trajectory of 7 models for all 8 tasks. The sequence length are 70time steps, while our training time steps are only 4 time steps. These sequences shows the ability of long term prediction fordifferent models. In all 2D tasks, our 4 models (the middle four rows in each figure) give prediction sequences that is closeto the ground truth. In the Gyro tasks, our CMr models give prediction sequences that is close to the ground truth. The 3baseline models give physically incorrect predictions in most of the tasks. A blank in a figure means some coordinate datacontain value Inf or NaN, and thus the scene cannot be generated. T r u e C M - C D - C L NN C M - C D - C H NN C M r - C D - C L NN C M r - C D - C H NN M L P - C D - C L NN I N - C P - C L NN t=0.00 I N - C P - S P t=0.05 t=0.10 t=0.15 t=0.20 t=0.25 t=0.30 t=0.35 t=0.40 t=0.45 t=0.50 t=0.55 t=0.60 t=0.65 t=0.70 Figure 7.
Prediction sequences of a test trajectory in the BM5-e task.
Differentiable Contact Model to Extend Lagrangian and Hamiltonian Neural Networks for Modeling Hybrid Dynamics T r u e C M - C D - C L NN C M - C D - C H NN C M r - C D - C L NN C M r - C D - C H NN M L P - C D - C L NN I N - C P - C L NN t=0.00 I N - C P - S P t=0.05 t=0.10 t=0.15 t=0.20 t=0.25 t=0.30 t=0.35 t=0.40 t=0.45 t=0.50 t=0.55 t=0.60 t=0.65 t=0.70 Figure 8.
Prediction sequences of a test trajectory in the BM5 task. T r u e C M - C D - C L NN C M - C D - C H NN C M r - C D - C L NN C M r - C D - C H NN M L P - C D - C L NN I N - C P - C L NN t=0.00 I N - C P - S P t=0.05 t=0.10 t=0.15 t=0.20 t=0.25 t=0.30 t=0.35 t=0.40 t=0.45 t=0.50 t=0.55 t=0.60 t=0.65 t=0.70 Figure 9.
Prediction sequences of a test trajectory in the CP3-e task.
Differentiable Contact Model to Extend Lagrangian and Hamiltonian Neural Networks for Modeling Hybrid Dynamics T r u e C M - C D - C L NN C M - C D - C H NN C M r - C D - C L NN C M r - C D - C H NN M L P - C D - C L NN I N - C P - C L NN t=0.00 I N - C P - S P t=0.05 t=0.10 t=0.15 t=0.20 t=0.25 t=0.30 t=0.35 t=0.40 t=0.45 t=0.50 t=0.55 t=0.60 t=0.65 t=0.70 Figure 10.
Prediction sequences of a test trajectory in the CP3 task. T r u e C M - C D - C L NN C M - C D - C H NN C M r - C D - C L NN C M r - C D - C H NN M L P - C D - C L NN I N - C P - C L NN t=0.00 I N - C P - S P t=0.05 t=0.10 t=0.15 t=0.20 t=0.25 t=0.30 t=0.35 t=0.40 t=0.45 t=0.50 t=0.55 t=0.60 t=0.65 t=0.70 Figure 11.
Prediction sequences of a test trajectory in the BD5 task.
Differentiable Contact Model to Extend Lagrangian and Hamiltonian Neural Networks for Modeling Hybrid Dynamics T r u e C M - C D - C L NN C M - C D - C H NN C M r - C D - C L NN C M r - C D - C H NN M L P - C D - C L NN I N - C P - C L NN t=0.000 I N - C P - S P t=0.025 t=0.050 t=0.075 t=0.100 t=0.125 t=0.150 t=0.175 t=0.200 t=0.225 t=0.250 t=0.275 t=0.300 t=0.325 t=0.350 Figure 12.
Prediction sequences of a test trajectory in the Rope task. T r u e C M - C D - C L NN C M - C D - C H NN C M r - C D - C L NN C M r - C D - C H NN M L P - C D - C L NN I N - C P - C L NN t=0.0 I N - C P - S P t=0.1 t=0.2 t=0.3 t=0.4 t=0.5 t=0.6 t=0.7 t=0.8 t=0.9 t=1.0 t=1.1 t=1.2 t=1.3 t=1.4 Figure 13.
Prediction sequences of a test trajectory in the Gyro-e task.