Next Article in Journal
End-to-End Continuous/Discontinuous Feature Fusion Method with Attention for Rolling Bearing Fault Diagnosis
Previous Article in Journal
A Comprehensive Collection and Analysis Model for the Drone Forensics Field
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Feedback Integrators for Mechanical Systems with Holonomic Constraints

by
Dong Eui Chang
1,*,
Matthew Perlmutter
2 and
Joris Vankerschaver
3,4
1
School of Electrical Engineering, Korea Advanced Institute of Science and Technology, Daejeon 34141, Korea
2
Department of Mathematics, Universidade Federal de Minas Gerais, Belo Horizonte 31270-901, Brazil
3
Center for Biosystems and Biotech Data Science, Ghent University Global Campus, Incheon 21985, Korea
4
Department of Applied Mathematics, Computer Science and Statistics, Ghent University, 9000 Ghent, Belgium
*
Author to whom correspondence should be addressed.
Sensors 2022, 22(17), 6487; https://doi.org/10.3390/s22176487
Submission received: 20 July 2022 / Revised: 24 August 2022 / Accepted: 26 August 2022 / Published: 29 August 2022
(This article belongs to the Section Sensors and Robotics)

Abstract

:
The feedback integrators method is improved, via the celebrated Dirac formula, to integrate the equations of motion for mechanical systems with holonomic constraints so as to produce numerical trajectories that remain in the constraint set and preserve the values of quantities, such as energy, that are theoretically known to be conserved. A feedback integrator is concretely implemented in conjunction with the first-order Euler scheme on the spherical pendulum system and its excellent performance is demonstrated in comparison with the RATTLE method, the Lie–Trotter splitting method, and the Strang splitting method.

1. Introduction

The method of feedback integrators was proposed in [1,2] to numerically integrate the equations of motion for dynamical systems in order to preserve their domain manifolds and first integrals. The method is summarized as follows. Suppose that there is a given invariant set Λ of a continuous-time dynamical system on a manifold P, where Λ can be the intersection of level sets of the first integrals of the system. Embed P into some Euclidean space, extend the system to the ambient Euclidean space, and modify it outside Λ to turn Λ into a local attractor of the resulting dynamical system in the ambient space. Then, trajectories originating from points in Λ generated by integration of the modified dynamics remain in Λ theoretically and near Λ numerically, irrespective of the choice of numerical integration schemes. It is rigorously shown [1,3,4] that the discrete-time dynamical system derived from any one-step numerical integrator with uniform step size h for the modified continuous-time system has an attractor Λ h that contains Λ in its interior and converges to Λ as h 0 + . In this procedure, the set of equations of motion of the modified dynamical system is called a feedback integrator. Feedback integrators can be implemented by any usual integration scheme such as Euler, Runge–Kutta, or Matlab ode45 in one single global Cartesian coordinate system for the ambient Euclidean space and do not require projecting numerical trajectories to a certain set or solving algebraic equations during integration.
Here, we propose a way to apply feedback integrators to mechanical systems with holonomic constraints, which is not addressed in [1,2]. Consider a symplectic manifold P, a symplectic submanifold S of P, and a Hamiltonian function H on P, where it is often the case that P = T * R n for some n and S is the cotangent bundle of a set of holonomic constraints in R n . The Hamiltonian function H defines a Hamiltonian vector field X H on P and the restriction H | S of H to S defines a Hamiltonian vector field, denoted X H | S , on S. In general, X H does not coincide with X H | S on S; so, we employ the celebrated Dirac formula to extend the vector field X H | S from S to P, such that the dynamical system extends from S to P. The manifold S is an invariant set of the extended dynamical system on P. Thus, we can apply the usual method of feedback integrators to the extended system on P such that S is numerically well-preserved.
This paper is organized as follows. The Dirac formula is first reviewed for the sake of completeness; then, the procedure for implementing feedback integrators for mechanical systems with holonomic constraints is presented. A simulation study is carried out on the spherical pendulum system, which is a mechanical system with a holonomic constraint, to demonstrate the excellent performance of feedback integrators in preserving the constraint set, the total energy, and the vertical component of the angular momentum vector, in comparison with the RATTLE method, the Lie–Trotter method, and the Strang method. Refer to [5] for more information on the three methods. We conclude with a small-scale simulation of the planar pendulum to show that the feedback method generally gives rise to integrators that are computationally more efficient as well.

Related Work

The numerical integration of mechanical systems with holonomic constraints has been an area of active research interest over the past decades. In this section, we briefly compare our approach with existing methods from the literature. Holonomic constraints can be viewed as index-1 differential-algebraic equations (DAEs) [6,7]. These systems can be integrated by direct discretization or by reformulating the DAE as an equivalent set of ODEs to which a standard numerical integration scheme may be applied. The resulting integration scheme typically involves the solution of a nonlinear equation representing the constraint. Another approach starts from the Hamiltonian or variational nature of mechanical systems to come up with discretizations that preserve the symplectic structure and the constraint manifold [5,8,9,10]. Such discretizations, the RATTLE algorithm [11] and its generalizations [12,13,14] chiefly among them, typically exhibit superior long-term integration properties compared with standard, nonsymplectic integration algorithms, and have found wide application in the numerical integration of mechanical and control systems [15,16,17]. This approach has also been extended to the case of classical field theories with constraints [18,19], or to systems with dissipation [20]. What all these methods have in common is that they seek to enforce the constraint equations directly, which requires special-purpose integration algorithms. By contrast, the feedback integrator method described in this paper modifies the equations of motion directly, to approximately conserve the constraint equations and other integrals of motion. As we pointed out before, the advantage of this approach is that it allows for standard, off-the-shelf numerical integrators to be used, such as the Euler or Runge–Kutta method.

2. Main Results

We first review the Dirac formula [21,22], explain how to construct feedback integrators for mechanical systems with holonomic constraints, and design a feedback integrator for the planar and spherical pendulum systems to demonstrate its excellent integration performance and versatility in comparison with the RATTLE method, the Lie–Trotter method, and the Strang method.

2.1. Review of the Dirac Formula

In this section, we recall the construction of Dirac for the decomposition of the Hamiltonian vector field along a symplectic submanifold S P , where ( P , ω ) is a symplectic manifold and S is symplectic, such that for all z S we have T z S T z S ω = T z P . Suppose that we are given a Hamiltonian function H on P. We will work semiglobally as follows. Suppose that S locally is expressed as the zero level set of a function f : P R 2 k that has 0 as a regular value. Then, we have dim S = 2 n 2 k where dim P = 2 n and we denote f = ( f 1 , , f 2 k ) . Since S is symplectic, we can restrict the Hamiltonian H to S and, pulling back the symplectic form to S, we have
ι S * ω ( z ) ( X H | S ( z ) , v z ) = d ( H | S ) ( z ) · v z
for v z T z S .
Proposition 1.
Each of the X f i | S is a section of T S ω . Furthermore, these sections are linearly independent and, at any point, span this distribution.
Proof. 
We have that v T z S if and only if d f i ( z ) · v = 0 for all i { 1 , 2 k } . Now, let v T z S . We then have, for each i, ω ( z ) ( X f i , v ) = d f i ( z ) · v = 0 ; therefore, X f i | S is a section of T S ω . From the independence of the f i , a dimension count shows that the linearly independent X f i ( z ) span T z S ω , completing the proof. □
Now, define for each z P a ( 2 k × 2 k ) matrix by
C i j ( z ) = { f i , f j } ( z ) .
We then have the following.
Proposition 2.
Fix z S . The matrix C i j ( z ) is invertible.
Proof. 
Fix z S . We know that T z S is symplectic and, therefore, T z S ω is too. By the previous proposition, we know that the X f i ( z ) ’s form a basis of the symplectic space T z S ω . Therefore, the entry of the matrix C i j ( z ) is precisely the symplectic form evaluated on the vectors X f i ( z ) , X f j ( z ) , which is, thus, a non-degenerate matrix since the restriction of ω to T z S ω is non-degenerate. □
We then have the following version of the Dirac formula for the Hamiltonian vector field.
Theorem 1.
With the definition of C i j ( z ) given as above, denoting its inverse by C i j ( z ) for z S , the following formula holds:
X H | S ( z ) = X H ( z ) i , j = 1 2 k C i j ( z ) { H , f i } ( z ) X f j ( z ) .
Proof. 
Since S is symplectic, we have that for all z S , T z P = T z S T z S ω . We will show that the projection π z : T z P T z S relative to this decomposition is given by the right-hand side of (1). This is equivalent to showing that I π z : T z P T z S ω is given by
( I π z ) X H ( z ) = i , j = 1 2 k { H , f i } ( z ) C i j ( z ) X f j ( z ) .
To establish this, first observe that the right-hand side lies in T z S ω by Proposition 1. Next, if H = f where { 1 , 2 k } , then the right-hand side of (2) is
i , j = 1 2 k C i ( z ) C i j ( z ) X f j ( z ) = j = 1 2 k δ j ( z ) X f j ( z ) = X f ( z )
which shows that the operator is the identity on the subspace T z S ω . Next, suppose that X H ( z ) T z S . Then, we have d f i ( z ) · X H ( z ) = 0 for all i and, therefore, the right-hand side of (2) vanishes, as required. This proves that the projection I π z is given by the Formula (2) and, thus, π z is given by the right-hand side of (1). □
Our main use of this theorem is to extend a Hamiltonian system on S to all of P. That is, given a Hamiltonian function H defined on P, we compute the right-hand side of Equation (1) for arbitrary z P to obtain a vector field defined on a neighborhood of the submanifold S (or level surface in this case) that will automatically be equal to X H | S on S itself. Thus, the Dirac formula gives us a way to extend Hamiltonian vector fields on S to the ambient space.

2.2. Feedback Integrators for Mechanical Systems with Holonomic Constraints

Given a Hamiltonian function H on a 2 n -dimensional symplectic manifold P and its Hamiltonian vector field X H on P, consider a Hamiltonian system Σ S on a ( 2 n 2 k ) -dimensional symplectic submanifold S of P whose Hamiltonian function is the restriction H | S of H to S. So, the equations of motion of Σ S can be written as
z ˙ = X H | S ( z ) , z S .
Assume that S is expressed as a level of a function f = ( f 1 , , f 2 k ) : P R 2 k . Thanks to the Dirac Formula (1), the dynamical system (3) extends to P as
x ˙ = X H ( x ) i , j = 1 2 k C i j ( x ) { H , f i } ( x ) X f j ( x ) , x P .
It is easy to verify that S is an invariant manifold of (4). We now wish to integrate the dynamics (4) from a point x 0 in S.
Suppose that we can embed the manifold P in some Euclidean space R m with m 2 n and extend the vector field (4) to R m , not necessarily in the Dirac way. Denote the extended vector field by X and write the corresponding dynamical system as
x ˙ = X ( x ) , x R m .
Suppose that both functions f and H also extend to the ambient Euclidean space R m in such a way that D f · X = 0 and H · X = 0 in R m , and that the manifold S is still a level set of the extension of f in R m . As a result, the functions f and H are first integrals of (5) and the manifold S is an invariant manifold of (5).
Suppose that there are first integrals I = ( I 1 , , I ) of (5) other than H and f, where the function I may include part of a function on R m whose level set defines P in R m . Define a function F : R m R 2 k + 1 + by F = ( f , H , I ) , and a function V on R m by
V ( x ) = 1 2 ( F ( x ) F ( x 0 ) ) K ( F ( x ) F ( x 0 ) )
where K = K T is a ( 2 k + 1 + ) × ( 2 k + 1 + ) constant positive definite symmetric matrix. The minimum value of V is 0 and the function V satisfies
V 1 ( 0 ) = { x R m x S , H ( x ) = H ( x 0 ) , I ( x ) = I ( x 0 ) } .
The set V 1 ( 0 ) is an invariant manifold of (5) since V · X = ( F ( x ) F ( x 0 ) ) T K · D F · X = ( F ( x ) F ( x 0 ) ) T K · 0 = 0 .
Modify the dynamical system (5) as follows:
x ˙ = X ( x ) L ( x ) V ( x ) , x R m ,
where L is a map from R m into the set of m × m positive definite symmetric matrices and V is computed as V ( x ) = D F ( x ) T K ( F ( x ) F ( x 0 ) ) . Since 0 is the minimum value of V, the gradient V vanishes on V 1 ( 0 ) , which implies that the two dynamical systems (5) and (6) coincide on V 1 ( 0 ) . It follows that the set V 1 ( 0 ) is an invariant manifold of (6). Along any trajectory of (6), V ( t ) is an nonincreasing function of t since d V d t = V , ( X L V ) = V , L V 0 ; so, the trajectory converges to V 1 ( 0 ) as t or it remains close to V 1 ( 0 ) if it starts near V 1 ( 0 ) . Under some conditions on V, the set V 1 ( 0 ) becomes a unique attractor of (6) in a neighborhood of itself in R m ; refer to [1] for those conditions. Due to the invariance of V 1 ( 0 ) and the coincidence of (5) and (6) on V 1 ( 0 ) , integrating (5) from x 0 and integrating (6) from x 0 produce the same trajectory. Numerically, however, integrating (6) has the following advantage: if the trajectory deviates from V 1 ( 0 ) at some numerical integration step, then it will get pushed back toward the attractor V 1 ( 0 ) , thus remaining on the manifold S and preserving all the first integrals; refer to [1,23] for a rigorous explanation. Although [1,23] provides some sufficient conditions for V 1 ( 0 ) to be a local attractor of (6), in practice, it is not necessary to verify them. The procedure outlined in this section is good enough for applications, which will be illustrated with the planar and spherical pendulum in the next sections.

2.3. The Spherical Pendulum

We build a feedback integrator for the spherical pendulum [24], which is a Hamiltonian system with a holonomic constraint, and compare its performance with that of such geometric numerical integration methods as the RATTLE method, the Lie–Trotter splitting method, and the Strang splitting method. The phase space of the spherical pendulum is T * S 2 , which is a symplectic submanifold of T * R 3 = R 3 × R 3 . This submanifold is globally defined as the ( 2 , 0 ) -level set of the function f = ( f 1 , f 2 ) : T * R 3 R 2 with > 0 , where f 1 ( q , p ) = | | q | | 2 and f 2 ( q , p ) = q · p for ( q , p ) T * R 3 . We fix the Hamiltonian H ( q , p ) = | | p | | 2 / 2 m + m g q 3 . Restricted to T * S 2 , this gives the Hamiltonian of the spherical pendulum under gravity, with its S 1 symmetry. In order to write down the extended Hamiltonian vector field, we compute the following:
{ H , f 1 } = 2 f 2 / m , { H , f 2 } = p 2 / m + m g q 3 , { f 1 , f 2 } = 2 f 1 , X H = ( p / m , m g e 3 ) , X f 1 = ( 0 , 2 q ) , X f 2 = ( q , p ) ,
and
[ C i j ] = 0 1 / ( 2 f 1 ) 1 / ( 2 f 1 ) 0 ,
where e 3 = ( 0 , 0 , 1 ) . Hence, by the Dirac formula, the spherical pendulum system extends to T * R 3 as
x ˙ = X ( x )
where
X ( x ) = 1 m p f 2 m f 1 q , m g e 3 + f 2 m f 1 p + p 2 m + m g q 3 1 f 1 q
for x = ( q , p ) T * R 3 . The extended system has four first integrals on T * R 3 : the constraint functions f 1 and f 2 , the Hamiltonian H, and the vertical component J ( q , p ) = q 1 p 2 q 2 p 1 of angular momentum.
Choose a point ( q 0 , p 0 ) T * S 2 , and define a function V : T * R 3 R by
V ( q , p ) = 1 2 k 1 | Δ f 1 | 2 + 1 2 k 2 | Δ f 2 | 2 + 1 2 k 3 | Δ H | 2 + 1 2 k 4 | Δ J | 2 ,
where k i > 0 , i = 1 , , 4 are constants; Δ f i = f i ( q , p ) f i ( q 0 , p 0 ) , i = 1 , 2 ; Δ H = H ( q , p ) H ( q 0 , p 0 ) ; and Δ J = J ( q , p ) J ( q 0 , p 0 ) . It is easy to verify that 0 is the minimum value of V and
V 1 ( 0 ) = { ( q , p ) T * R 3 f 1 ( q , p ) = 2 , f 2 ( q , p ) = 0 , H ( q , p ) = H 0 , J ( q , p ) = J 0 } ,
where H 0 = H ( q 0 , p 0 ) and J 0 = J ( q 0 , p 0 ) . The gradient of V is computed as
V = k 1 Δ f 1 f 1 + k 2 Δ f 2 f 2 + k 3 Δ H H + k 4 Δ J J ,
where f 1 = ( 2 q , 0 ) , f 2 = ( p , q ) , H = ( m g e 3 , p / m ) , and J = ( p × e 3 , e 3 × q ) . Then, the feedback integrator corresponding to the function V for the spherical pendulum is given by
x ˙ = X ( x ) V ( x )
or
q ˙ = 1 m p f 2 m f 1 q 2 k 1 Δ f 1 q k 2 Δ f 2 p k 3 m g Δ H e 3 k 4 Δ J ( p × e 3 )
p ˙ = m g e 3 + f 2 m f 1 p + p 2 m + m g q 3 1 f 1 q k 2 Δ f 2 q k 3 Δ H m p k 4 Δ J ( e 3 × q ) .
We now compare the feedback integrator with the RATTLE method, the Lie–Trotter method, and the Strang method on the spherical pendulum system. The RATTLE algorithm is given on p. 246 in [5] and is known to be symplectic and convergent of order two [5]. The Lie–Trotter method and the Strang method are so-called splitting methods and are explained on pp. 253–254 in [5], where the two splitting methods yield first- and second-order numerical integrators, respectively [5]. For the splitting methods, we split the Hamiltonian function H into the kinetic function H [ 1 ] ( q , p ) = | | p | | 2 / 2 m and the potential function H [ 2 ] ( q , p ) = m g q 3 , and we note that the dynamics associated to H [ 1 ] and H [ 2 ] can be integrated analytically. For numerical simulation, choose the parameter values m = g = = 1 for convenience, and the initial points q ( 0 ) = ( 0 , 1 , 0 ) and p ( 0 ) = ( 1 , 0 , 1 ) on T * S 2 . The corresponding initial values of the first integrals f 1 , f 2 , H, and J are f 10 = 1 , f 20 = 0 , H 0 = 2 , and J 0 = 1 , respectively. We fix the time step size h = 10 3 and the time interval [ 0 , 100 ] for integration by all four methods. For the feedback integrator, the usual Euler scheme is used to integrate (7) and (8) with the following gain values: k 1 = 50 , k 2 = 50 , k 3 = 50 , and k 4 = 50 .
The computational results are plotted in Figure 1, Figure 2, Figure 3, Figure 4 and Figure 5. In Figure 1, the trajectories q ( t ) = ( q 1 ( t ) , q 2 ( t ) , q 3 ( t ) ) of the configuration variables generated by the four methods are plotted. The feedback integrator with the Euler, RATTLE, Lie–Trotter, and Strang splitting methods all generate similar trajectories. Figure 2 and Figure 3 show the plots of the deviations | Δ f 1 ( t ) | and | Δ f 2 ( t ) | from the constraint manifold T * S 2 . The result by the RATTLE method is the best, and the trajectory produced by the feedback integrator with Euler remains close to T * S 2 with the step size h = 10 3 taken into account. Likewise, the trajectories by the Lie–Trotter method and the Strang method stay close to T * S 2 , because the flows of the split Hamiltonians H [ 1 ] and H [ 2 ] each preserve T * S 2 , as does their composition.
The feedback integrator with Euler and the RATTLE method perform well in preservation of the values of the Hamiltonian H as do the other two methods, as shown in Figure 4. The vertical component J of angular momentum is well-preserved by all four methods as shown in Figure 5, where it is noticeable that the Lie–Trotter method and the Strang method perform very well in preservation of J.
The computational results imply that the feedback integrator with the first-order Euler scheme performs well on the spherical pendulum system in comparison with the well-known RATTLE method that is of second order, and to the Lie–Trotter method and the Strang method. An advantage of the feedback integrator over the other three methods is that it does not require any projection or solving of algebraic equations to stay on the holonomic constraint manifold. Further, it does not require any special integration algorithms, and simply employs well-known integration algorithms available such as Euler, Runge–Kutta, or Matlab ode45. Moreover, unlike the Lie–Trotter and Strang methods, which require a particular splitting of the dynamics into two parts that are separately integrable, the feedback integrator can be made to work for any constrained dynamical system by modifying the vector field as in (6).

2.4. The Simple Pendulum

One potential drawback of the feedback integration method is that the feedback vector field is more complex due to the presence of the stabilizing forces. For example, the right-hand side of (7) and (8) is a sum of five gradients (one gradient of the Hamiltonian and four gradients of the feedback potential) and this cost compounds for higher-order numerical methods, which typically require several force evaluations per integration step. This added computational cost must be taken into account when comparing the error profile of feedback integrators with that of standard methods, such as RATTLE, which only evaluate the gradient of the Hamiltonian (but possibly multiple times per integration step).
In this section, we show that the increase in complexity is compensated by the approximate conservation properties of the integrator, and in particular, we show that feedback integrators are at least as effective as standard integrators when the computational budget is taken into account. We compare the performance of feedback integrators of different orders with that of the RATTLE method and show that the increase in computational cost is balanced by the fact that feedback integrators require fewer force evaluations overall to achieve a given accuracy.
To assess the global error of the various integrators, we take recourse to a simpler mechanical system, for which exact solutions are known and can be approximated with arbitrary precision: the simple gravitational pendulum. The simple pendulum consists of a mass m that is free to swing at the end of a rigid rod of length under the influence of gravity. The motion of the pendulum takes place entirely in a fixed plane, denoting the angle between the horizontal and the position of the pendulum by θ , is determined by
θ ¨ + g cos θ = 0 ,
where g is the gravitational acceleration. The dynamics of the pendulum as a constrained system can be derived directly via a calculation as in the previous section, or by observing that the spherical pendulum naturally moves in a fixed plane if the initial position, momentum, and direction of gravity are all coplanar. Either way, the extended Hamiltonian vector field is readily seen to be
X ( x ) = 1 m p f 2 m f 1 q , m g e y + f 2 m f 1 p + p 2 m + m g y 1 f 1 q ,
where q = ( x , y ) and p = ( p x , p y ) are the coordinates and momenta of the pendulum, respectively, and e y = ( 0 , 1 ) . The first integrals on T * R 2 are the Hamiltonian H ( q , p ) = p 2 / 2 m + m g y and the constraint functions
f 1 ( q , p ) = q 2 = x 2 + y 2 and f 2 ( q , p ) = q · p = x p x + y p y .
Similar to the spherical pendulum, we consider these three conserved quantities and for given initial values ( q 0 , p 0 ) T * S 2 , we define the function V : T * R 2 R given by
V ( q , p ) = 1 2 k 1 | Δ f 1 | 2 + 1 2 k 2 | Δ f 2 | 2 + 1 2 k 3 | Δ H | 2 ,
where k 1 , k 2 , k 3 are positive constants; Δ f i = f i ( q , p ) f i ( q 0 , p 0 ) ; i = 1 , 2 ; and Δ H = H ( q , p ) H ( q 0 , p 0 ) . The feedback integrator for the simple pendulum then becomes
q ˙ = 1 m p f 2 m f 1 q 2 k 1 Δ f 1 q k 2 Δ f 2 p k 3 m g Δ H e 2
p ˙ = m g e 2 + f 2 m f 1 p + p 2 m + m g y 1 f 1 q k 2 Δ f 2 q k 3 Δ H m p .
While the pendulum can, in principle, be integrated exactly, obtaining the solution as a function of time is not straightforward and requires inverting the elliptic integral of the first kind. To avoid this difficulty, we integrate the Equation (9) for θ using a high-order Runge–Kutta method with tolerance set to 10 13 . For all numerical simulations, we set m = g = = 1 , and we use q ( 0 ) = ( 1 , 0 ) and p ( 0 ) = ( 0 , 0 ) as the initial conditions. For the feedback method, we use three off-the-shelf numerical schemes to integrate (10) and (11): (a) the forward Euler method, (b) the explicit 4th-order Runge–Kutta method (RK4), and (c) the Dormand–Prince method of 8th order (DOP853). Note that the Dormand–Prince method is a variable step-size integrator while all the others use a fixed step size.
Figure 6 (left) shows the trajectory error after one period of the pendulum between the standard RATTLE integrator on the one hand, and the three feedback integrators on the other hand, as a function of the step size. The global error decreases as the step size decreases, in line with the order of the method. This is to be expected, since the feedback approach merely modifies the vector field to be integrated, but does not otherwise alter the underlying numerical method.
The vector field integrated by the feedback methods is more complex; thus, one can ask what the impact is on the total execution time. To investigate this question, we give each method a fixed computational budget and modify the integrator code so that each evaluation of the vector field reduces the computational budget by a certain amount. For evaluations of the feedback vector field (10) and (11) the cost per evaluation is set to 4, since the vector field consists of 4 gradients, whereas for the evaluation of the gradient of the Hamiltonian (needed e.g., by RATTLE), the cost per evaluation is 1. We then compare the performance of the integrators over one period of the pendulum and adjust the step size (for RATTLE and the feedback Euler and RK4 methods) or the tolerance (for the feedback-DOP853 method, which uses a variable step size) so that the computational budget is exhausted over one period.
The result is shown in Figure 6 (right), which shows the global error as a function of the computational budget. Note that the error for the feedback-DOP853 method stabilizes somewhat below 10 12 . This is roughly the point where we encounter the limits in the accuracy of the exact trajectories (which were obtained by numerical integration of (9), where the tolerance was set to 10 13 ).
We see that feedback integrators are able to integrate the dynamics of the underlying system accurately (i.e., with low error) and efficiently (using comparable or lower numbers of force evaluations), compared with specific holonomic integrators. In terms of computational efficiency, higher-order integrators such as the feedback-DOP853 integrator clearly achieve better results than others (lower-order feedback methods and RATTLE). This again demonstrates one of the key benefits of the feedback integrator method, showing that it is possible to use any standard numerical integration scheme to obtain approximate constraint preservation, without loss of accuracy or computational efficacy.

3. Conclusions

We have presented a general framework to extend the feedback integrators [1] to systems with holonomic constraints. Beginning with a symplectic submanifold S in the symplectic manifold P, where S—the holonomic constraint—is the level set of f 1 f 2 k , on which we have a Hamiltonian H | S . We use the Dirac formula to extend the vector field X H | S to a vector field X on P. We then apply the feedback integrator to X using an embedding of P in Euclidean space.
More specifically, in the case that the symplectic manifold S is of the form T * Q , where Q embeds in R n and T * Q embeds in T * R n as the level set of functions f 1 , f 2 k , on which we have an extended Hamiltonian H, we can compute the extension of the vector field X H | T * Q directly from the Dirac formula (1). With this vector field, now defined on a Euclidean space, T * R n , we construct the function V : T * R n R 0 , whose 0 level set contains the dynamic invariants for a given initial point in S, where the set of dynamic invariants includes the functions f 1 , f 2 k . The feedback-integrator-modified vector field, constructed from the extended vector field X by adding the negative gradient of V, is fed into any integrator, for example, first-order Euler. In addition to the dynamic invariants, the integrator automatically respects, from the construction of the modified vector field, the holonomic constraints. The resulting integrator has superior performance to even symplectic integrators, which depend, typically, on implicit solvers. As future work, we plan to examine the problem of optimal control on manifolds and its relationship with Hamiltonian mechanics from the viewpoint of feedback integrators about which some preliminary works have been carried out in [25,26,27]. We also plan to examine the effect of the magnitude of feedback integrator term on the precision of numerical integration and to extend feedback integrators to field theory.

Author Contributions

Conceptualization, D.E.C. and M.P.; methodology, D.E.C., M.P. and J.V.; software, D.E.C., M.P. and J.V.; validation, D.E.C., M.P. and J.V.; formal analysis, D.E.C., M.P. and J.V.; investigation, D.E.C., M.P. and J.V.; writing—original draft preparation, D.E.C., M.P. and J.V.; writing—review and editing, D.E.C., M.P. and J.V.; visualization, D.E.C., M.P. and J.V.; supervision, D.E.C.; project administration, D.E.C.; funding acquisition, D.E.C. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by the IITP grant (No.20210005900012003).

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Not applicable.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Chang, D.E.; Jiménez, F.; Perlmutter, M. Feedback integrators. J. Nonlinear Sci. 2016, 26, 1693–1721. [Google Scholar] [CrossRef]
  2. Chang, D.E.; Perlmutter, M. Feedback integrators for nonholonomic mechanical systems. J. Nonlinear Sci. 2019, 29, 1165–1204. [Google Scholar] [CrossRef]
  3. Kloeden, P.E.; Lorenz, J. Stable attracting sets in dynamical systems and in their one-step discretizations. Siam J. Numer. Anal. 1986, 23, 986–995. [Google Scholar] [CrossRef]
  4. Kloeden, P.E.; Lorenz, J. A note on multistep methods and attracting sets of dynamical systems. Numer. Math. 1989, 56, 667–673. [Google Scholar] [CrossRef]
  5. Hairer, E.; Wanner, G.; Lubich, C. Geometric Numerical Integration: Structure-Preserving Algorithms for Ordinary Differential Equations, 2nd ed.; Springer: Berlin/Heidelberg, Germany, 2006. [Google Scholar]
  6. Hairer, E.; Lubich, C.; Roche, M. The Numerical Solution of Differential-Algebraic Systems by Runge-Kutta Methods; Lecture Notes in Mathematics; Springer: Berlin/Heidelberg, Germany, 1989; Volume 1409. [Google Scholar] [CrossRef]
  7. Kunkel, P.; Mehrmann, V. Differential-Algebraic Equations: Analysis and Numerical Solution; EMS Textbooks in Mathematics, European Mathematical Society (EMS): Zürich, Switzerland, 2006. [Google Scholar] [CrossRef]
  8. Leimkuhler, B.; Reich, S. Simulating Hamiltonian Dynamics; Cambridge Monographs on Applied and Computational Mathematics; Cambridge University Press: Cambridge, UK, 2004; Volume 14. [Google Scholar]
  9. Marsden, J.E.; West, M. Discrete Mechanics and Variational Integrators. Acta Numer. 2001, 10, 357–514. [Google Scholar] [CrossRef]
  10. Man, S.; Gao, Q.; Zhong, W. Variational Integrators in Holonomic Mechanics. Mathematics 2020, 8, 1358. [Google Scholar] [CrossRef]
  11. Andersen, H.C. Rattle: A “Velocity” Version of the Shake Algorithm for Molecular Dynamics Calculations. J. Comput. Phys. 1983, 52, 24–34. [Google Scholar] [CrossRef]
  12. McLachlan, R.I.; Modin, K.; Verdier, O.; Wilkins, M. Symplectic Integrators for Index One Constraints. Siam J. Sci. Comput. 2013, 35, A2150–A2162. [Google Scholar] [CrossRef] [Green Version]
  13. McLachlan, R.I.; Modin, K.; Verdier, O.; Wilkins, M. Geometric Generalisations of SHAKE and RATTLE. Found. Comput. Math. 2014, 14, 339–370. [Google Scholar] [CrossRef]
  14. Barbero-Liñán, M.; de Diego, D.M. Retraction Maps: A Seed of Geometric Integrators. Found. Comput. Math. 2022. [Google Scholar] [CrossRef]
  15. Leok, M.; Zhang, J. Discrete Hamiltonian Variational Integrators. Ima J. Numer. Anal. 2011, 31, 1497–1532. [Google Scholar] [CrossRef]
  16. Colombo, L.; de Diego, D.M. Higher-Order Variational Problems on Lie Groups and Optimal Control Applications. J. Geom. Mech. 2014, 6, 451. [Google Scholar] [CrossRef]
  17. Wenger, T.; Ober-Blöbaum, S.; Leyendecker, S. Construction and Analysis of Higher Order Variational Integrators for Dynamical Systems with Holonomic Constraints. Adv. Comput. Math. 2017, 43, 1163–1195. [Google Scholar] [CrossRef]
  18. Tran, B.; Leok, M. Multisymplectic Hamiltonian Variational Integrators. Int. J. Comput. Math. 2022, 99, 113–157. [Google Scholar] [CrossRef]
  19. Demoures, F.; Gay-Balmaz, F. Multisymplectic Variational Integrators for Fluid Models with Constraints. In Geometric Science of Information; Nielsen, F., Barbaresco, F., Eds.; Springer International Publishing: Cham, Switzerland, 2021; Volume 12829, pp. 283–291. [Google Scholar] [CrossRef]
  20. Bhatt, A. Projected Exponential Runge–Kutta Methods for Preserving Dissipative Properties of Perturbed Constrained Hamiltonian Systems. J. Comput. Appl. Math. 2021, 394, 113556. [Google Scholar] [CrossRef]
  21. Abraham, R.; Marsden, J. Foundations of Mechanics, 2nd ed.; Addison-Wesley Publishing Company: Redwood City, CA, USA, 1980. [Google Scholar]
  22. Marsden, J.; Ratiu, T. Introduction to Mechanics and Symmetry, 2nd ed.; Springer: New York, NY, USA, 1999. [Google Scholar]
  23. Chang, D.E. On controller design for systems on manifolds in Euclidean space. Int. J. Robust Nonlinear Control. 2018, 28, 4981–4998. [Google Scholar] [CrossRef]
  24. Marsden, J. Lectures on Mechanics; Cambridge University Press: Cambridge, UK, 1992. [Google Scholar]
  25. Chang, D.E. A simple proof of the Pontryagin maximum principle on manifolds. Automatica 2011, 47, 630–633. [Google Scholar] [CrossRef]
  26. Rojo, A.; Bloch, A. The Principle of Least Action: History and Physics; Cambridge University Press: Cambridge, UK, 2018. [Google Scholar]
  27. Phogat, K.S.; Chang, D.E. Model predictive regulation on manifolds in Euclidean space. Sensors 2022, 22, 5170. [Google Scholar] [CrossRef] [PubMed]
Figure 1. The trajectories of the position q ( t ) = ( q 1 ( t ) , q 2 ( t ) , q 3 ( t ) ) , 0 t 100 , of the spherical pendulum generated by four different methods with step size h = 10 3 : a feedback integrator with the Euler scheme, the RATTLE method, the Lie–Trotter splitting method, and the Strang splitting method.
Figure 1. The trajectories of the position q ( t ) = ( q 1 ( t ) , q 2 ( t ) , q 3 ( t ) ) , 0 t 100 , of the spherical pendulum generated by four different methods with step size h = 10 3 : a feedback integrator with the Euler scheme, the RATTLE method, the Lie–Trotter splitting method, and the Strang splitting method.
Sensors 22 06487 g001
Figure 2. The trajectories of the deviation, | Δ f 1 ( t ) | = | f 1 ( t ) f 1 ( 0 ) | , 0 t 100 , of the spherical pendulum from the constraint set, q 2 = 1 , generated by four different methods with step size h = 10 3 : a feedback integrator with the Euler scheme, the RATTLE method, the Lie–Trotter splitting method, and the Strang splitting method. A logarithmic scale is used on the y-axis.
Figure 2. The trajectories of the deviation, | Δ f 1 ( t ) | = | f 1 ( t ) f 1 ( 0 ) | , 0 t 100 , of the spherical pendulum from the constraint set, q 2 = 1 , generated by four different methods with step size h = 10 3 : a feedback integrator with the Euler scheme, the RATTLE method, the Lie–Trotter splitting method, and the Strang splitting method. A logarithmic scale is used on the y-axis.
Sensors 22 06487 g002
Figure 3. The trajectories of the deviation, | Δ f 2 ( t ) | = | f 2 ( t ) f 2 ( 0 ) | , 0 t 100 , of the spherical pendulum from the constraint set, q · p = 0 , generated by four different methods with step size h = 10 3 : a feedback integrator with the Euler scheme, the RATTLE method, the Lie–Trotter splitting method, and the Strang splitting method. A logarithmic scale is used on the y-axis.
Figure 3. The trajectories of the deviation, | Δ f 2 ( t ) | = | f 2 ( t ) f 2 ( 0 ) | , 0 t 100 , of the spherical pendulum from the constraint set, q · p = 0 , generated by four different methods with step size h = 10 3 : a feedback integrator with the Euler scheme, the RATTLE method, the Lie–Trotter splitting method, and the Strang splitting method. A logarithmic scale is used on the y-axis.
Sensors 22 06487 g003
Figure 4. The trajectories of the energy error, | Δ H ( t ) | = | H ( t ) H ( 0 ) | , 0 t 100 , of the spherical pendulum generated by four different methods with step size h = 10 3 : a feedback integrator with the Euler scheme, the RATTLE method, the Lie–Trotter splitting method, and the Strang splitting method. A logarithmic scale is used on the y-axis.
Figure 4. The trajectories of the energy error, | Δ H ( t ) | = | H ( t ) H ( 0 ) | , 0 t 100 , of the spherical pendulum generated by four different methods with step size h = 10 3 : a feedback integrator with the Euler scheme, the RATTLE method, the Lie–Trotter splitting method, and the Strang splitting method. A logarithmic scale is used on the y-axis.
Sensors 22 06487 g004
Figure 5. The trajectories of the momentum error, | Δ J ( t ) | = | J ( t ) J ( 0 ) | , 0 t 100 , of the spherical pendulum generated by four different methods with step size h = 10 3 : a feedback integrator with the Euler scheme, the RATTLE method, the Lie–Trotter method, and the Strang method. A logarithmic scale is used on the y-axis.
Figure 5. The trajectories of the momentum error, | Δ J ( t ) | = | J ( t ) J ( 0 ) | , 0 t 100 , of the spherical pendulum generated by four different methods with step size h = 10 3 : a feedback integrator with the Euler scheme, the RATTLE method, the Lie–Trotter method, and the Strang method. A logarithmic scale is used on the y-axis.
Sensors 22 06487 g005
Figure 6. Left: The global trajectory error after one period of the pendulum (approximately 7.4 time units) as a function of the step size (DOP853 is not included as it is a variable step size method). As the step size decreases, the global error decreases at a rate proportional to the error of the method. Right: The global error, but now as a function of the computational budget (number of force evaluations). Larger computational budgets correspond to smaller step sizes and, hence, lower errors, but take into account the fact that the feedback methods involve more force evaluations. Despite the overhead, feedback integrators are able to do at least as well as, or better than, RATTLE.
Figure 6. Left: The global trajectory error after one period of the pendulum (approximately 7.4 time units) as a function of the step size (DOP853 is not included as it is a variable step size method). As the step size decreases, the global error decreases at a rate proportional to the error of the method. Right: The global error, but now as a function of the computational budget (number of force evaluations). Larger computational budgets correspond to smaller step sizes and, hence, lower errors, but take into account the fact that the feedback methods involve more force evaluations. Despite the overhead, feedback integrators are able to do at least as well as, or better than, RATTLE.
Sensors 22 06487 g006
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Chang, D.E.; Perlmutter, M.; Vankerschaver, J. Feedback Integrators for Mechanical Systems with Holonomic Constraints. Sensors 2022, 22, 6487. https://doi.org/10.3390/s22176487

AMA Style

Chang DE, Perlmutter M, Vankerschaver J. Feedback Integrators for Mechanical Systems with Holonomic Constraints. Sensors. 2022; 22(17):6487. https://doi.org/10.3390/s22176487

Chicago/Turabian Style

Chang, Dong Eui, Matthew Perlmutter, and Joris Vankerschaver. 2022. "Feedback Integrators for Mechanical Systems with Holonomic Constraints" Sensors 22, no. 17: 6487. https://doi.org/10.3390/s22176487

APA Style

Chang, D. E., Perlmutter, M., & Vankerschaver, J. (2022). Feedback Integrators for Mechanical Systems with Holonomic Constraints. Sensors, 22(17), 6487. https://doi.org/10.3390/s22176487

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop