Simulation of Nonholonomic Mechanical Systems Using Algorithmic Differentiation

Simulation of Nonholonomic Mechanical Systems Using Algorithmic Differentiation Matthias Franke ∗ Tobias Zaiczek ∗ Klaus R¨ obenack ∗∗ ∗ Frauhofer Institute Integrated Circuits, Design Automation Division, Zeunerstraße 38, 01069 Dresden, Germany (e-mail: {Matthias.Franke, Tobias.Zaiczek}@eas.iis.fraunhofer.de). ∗∗ Institute of Control Theory, Electrical Engineering Department, Technische Universit¨ at Dresden, Germany, (e-mail: [email protected])

Abstract: The purpose of this contribution is to introduce a new approach for simulating mechanical systems with nonholonomic constraints. The technique of algorithmic differentiation is employed together with Hamel ’s equations of motion permitting to incorporate the nonholonomic constraints directly instead of generating a system of differential-algebraic equations. To apply the proposed simulation approach, the user only needs to state the kinetic and potential energies, the nonconservative forces, and the nonholonomic constraints of the system to be simulated. The procedure is demonstrated by an illustrating example of a nonholonomic moving robot. Keywords: Mechanical system, nonholonomic constraint, Lagrange multiplier, algorithmic differentiation, quasi-velocities, Hamel’s equations 1. INTRODUCTION We consider mechanical systems with n degrees of freedom and p nonholonomic constraints. The system trajectories evolve on the n-dimensional smooth configuration manifold Q with local coordinates q = (q 1 , . . . , q n ). The Lagrangian L, which is defined as kinetic minus potential energy, maps from the tangent bundle T Q, with local coordinates (q 1 , . . . , q n , q˙1 , . . . , q˙n ), to the real numbers, see e.g. Abraham and Marsden (1978) for details. The nonholonomic constraints are represented by a system of one-forms ω i ∈ T ∗ Q, i = 1, . . . , p, in coordinates ω i = Pn i k ω (q)dq , which restrict the solutions to a subbundle k=1 k of the tangent bundle T Q. In coordinates this restriction reads as n X ω i cq˙ = ω i (q) ˙ = ωki (q)q˙k = 0, i = 1, . . . , p, k=1

i.e., the constraints are linear with respect to the velocities q˙ = (q˙1 , . . . , q˙n ). Constraints of that type are called nonholonomic if they are not integrable: Pp There are no mappings hi : Q → R such that dhi = j=1 µij ω j , where d denotes the exterior derivative on Q. The coefficients µij play the role of an integrating factor, i.e., there are also no integrable linear combinations 1 of the constraints. One can check if such an integrable linear combination exists by a dual version of Frobenius’ theorem, see Lee (2002) for example. Thus in the case of nonholonomic constraints the condition i

1

p

dω ∧ ω ∧ · · · ∧ ω = 0,

i = 1, .., p

is violated. 1

Only regular matrices (µij )i=1,...,p j=1,...,p are taken into account.

One way of dealing with this type of system is utilizing Lagrange multipliers to represent the reaction forces of the constraints, see e.g. Goldstein (1980), Neimark and Fufaev (1972). This leads to the well-known system of equations p X ˙ ∂L(q, q) ˙ d ∂L(q, q) − = f + λi ωki (q), k dt ∂ q˙k ∂q k i=1 n X

ωki (q)q˙k = 0,

k = 1, . . . , n,

i = 1, . . . , p.

(1)

k=1

Here fk is the nonconservative generalized force that corresponds to the coordinate q k . By introducing the variables v k = q˙k , this system of equations can be transformed to 2n first order differential equations together with p algebraic constraints. The numeric solution of such a system of differential algebraic equations (DAE, index 2 in general) is more involved than the solution of a system of ordinary differential equations (ODE-system), see Griepentrog and M¨arz (1986); Hairer et al. (1989); Brenan et al. (1996); R¨obenack et al. (2001); Kunkel and Mehrmann (2006). Moreover, the dimension of the problem has increased since the p Lagrange multipliers appear in the equations. Even though it is possible to eliminate the Lagrange multipliers by an appropriate projection to the admissible subbundle of T Q, we shall pursue a more direct approach due to Hamel (1904) here. After presenting Hamel ’s equations in Section 2, we give a short introduction on algorithmic differentiation. It is also shown in Section 3 how algorithmic differentiation is used to evaluate the equations of motion. Finally in Section 4, an illustrating example including a detailed derivation of the equations of motion by Hamel’s approach is presented.

2. HAMEL’S EQUATIONS AND QUASI-VELOCITIES Hamel’s equations are a generalization of the Lagrange equations of motion, see Neimark and Fufaev (1972). Since they also cover the special case of Euler’s equations of rigid body motion, they originally had been called LagrangeEuler equations by Hamel (1904). The crucial advantage of Hamel ’s approach is that so-called quasi-velocities n X ω i (q) ˙ = ωki (q)q˙k , i = 1, . . . , n (2) k=1

can be used, which are not necessarily a differential 2 of a corresponding coordinate on Q. This change of velocity coordinates is required to be invertible 3 : n X q˙k (ω) = $ik (q)ω i . (3) i=1

Moreover, it is reasonable to choose the first p velocities according to the nonholonomic constraints (1). We can formulate the Langrangian depending on the new coor˜ 1 , . . . , q n , ω 1 , . . . , ω n ) = L(q, ˜ ω) = L(q, q(ω)). dinates: L(q ˙ Note that the linear terms of ω i , i = 1, .., p must not be neglected in the kinetic energy even if the nonholonomic constraints (1) hold, see Hamel (1949) for a detailed discussion on that issue. The equations of motion can subsequently be written as n ˜ ω) ∂ q˙j (ω) ˜ ω) X ∂ L(q, d ∂ L(q, − + k dt ∂ω k ∂q j | ∂ω {z } j=1 j $k (q)

+

n X

n X

j=1 l=p+1

(4)

˜ ω) j ∂ L(q, γl,k (q)ω l = f˜k , ∂ω j

k = p + 1, . . . , n, where n n X X ∂ q˙i (ω) f = $ki (q)fi f˜k = i k ∂ω i=1 i=1 j and γl,k (q) is given by  n  X ∂ωνj (q) ∂ωσj (q) j γl,k (q) = − $lσ (q)$kν (q). σ ν ∂q ∂q ν,σ=1

These equations together with n X q˙i = $ki (q)ω k ,

i = 1, . . . , n,

time, the system contains plenty of partial derivatives of first and second order, which can be calculated efficiently using algorithmic differentiation. 3. ALGORITHMIC DIFFERENTIATION In order to simulate the equations of motion, we have to compute first and second order derivatives. Numerical differentiation by divided differences is not well-suited for this task due to truncation and cancellation errors. Therefore, the derivatives occurring in these equations are usually computed symbolically with computer algebra packages or libraries. We suggest the use of an alternative differentiation technique known as automatic or algorithmic differentiation, see Griewank and Walther (2008). Assume that the function under consideration is a sequence of elementary functions and operations. Derivatives of this function can be calculated by applying elementary differentiation rules to this sequence. In algorithmic differentiation, all intermediate values are floating point numbers instead of symbolic expressions. Thus, compared to symbolic calculations, this approach requires less memory for the exact calculation of the partial derivatives. Consider a smooth map F : Ω → Rm defined on an open subset Ω ⊆ Rn . In the so-called forward mode of algorithmic differentiation, each elementary statement in the program or algorithm describing F is not only evaluated with respect to the function value, but also with respect to the derivative value using elementary differentiation rules together with the chain rule. For a given point x ∈ Ω, where the function is to be evaluated, and a direction described by the vector v ∈ Rn , the forward mode yields the function value z ∈ Rm and the directional derivative w ∈ Rm simultaneously: z = F (x) 0

w = F (x)v . (5)

(6)

k=p+1

describe the motion of the mechanical system. By solving the equations (4) for ω˙ k , k = p + 1, .., n, we get an explicit system of 2n − p first order differential equations. Hence, the dimension has been decreased by 2p compared to the Lagrange multiplier approach. Instead of calculating Hamel’s equations (4),(6) in advance, we want to do this automatically during the simulation. After resolving the total derivative with respect to 2 From a geometric point of view, it might be more adequate to employ the first jet bundle instead of the tangent bundle: quasivelocities do in this terms not correspond to a prolongation of a section on T × Q with time manifold T , see for example Saunders (1989). 3 As is commonly done in the physics literature, we denote one-forms and coordinates by the same symbol.

(7a) (7b)

Algorithmic differentiation can be implemented by several techniques such as source code transformation, operator overloading or expression templates, see Corliss et al. (2002). The implementation of the forward mode can easily be sketched using operator overloading in C++. Assume that the function to be differentiated is given as C or C++ code using the floating point type double for real-valued variables. We have to replace all occurances of double by a new class, e.g. ddouble, containing not only the function value but also the value of the directional derivative for each intermediate variable: class ddouble { public: double val; // function value double der; // derivative value } Then, one has to overload all relevant floating point operations and functions such as +, −, ∗, / and sin, cos, pow, . . . such that apart from the function values the derivative values are calculated according to the elementary dif-

ferentiation rules. This approach is very similar to the implementation of complex variables in C++. The name forward mode is derived from the fact that the program flows for the function and the derivative are evaluated in the same direction. In the reverse mode, the derivative calculation is carried out in reverse order compared to the function evaluation. The reverse mode can be understood as a generalization of the back-propagation algorithm known from neuronal networks. For a given covector z¯ ∈ (Rm )∗ , the reverse mode yields the weighted derivative x ¯ ∈ (Rn )∗ as x ¯ = z¯F 0 (x) , (8) where (Rm )∗ denotes the dual space of the vector space Rm . If the elements of Rm are written as column vectors, the elements of (Rm )∗ can be described by row vectors. Forward and reverse mode result in a first order derivative in terms of a vector w or a covector x ¯, respectively. To compute the full Jacobian matrix we either calculate (7a) for all n unit vectors of Rn or all m unit covectors of (Rm )∗ . Therefore, the forward mode is used if n ≤ m, whereas the reverse mode is used for n > m. If the map F is a functional, i.e., m = 1, we set z¯ = 1 and obtain the full gradient F 0 (x) evaluated at x ∈ Ω in a single pass of the reverse mode. For this reason, the reverse mode is often used in optimization to calculate derivative values of a cost functional. Moreover, forward and reverse mode can be combined to obtain second order derivatives. We apply the reverse mode not only to the result of the pure function evaluation (7a), but to the whole result (7) of the forward mode. In addition to (8) we get y¯ = z¯F 00 (x)v , (9) which is the second order derivative of F at x. In particular, an explicit calculation of the derivative tensor F 00 , which is a tensor of contravariant degree 1 and covariant degree 2, is not required. The abovementioned combination of forward and reverse mode yields the contracted expression y¯ ∈ (Rn )∗ immediately. In case of m = 1, F 00 (x) becomes the Hessian matrix. Recall that a first order derivative F 0 (x) is formally a linear map. Similarly, a second order derivative F 00 (x) can be interpreted as a bilinear map, i.e., a map which is linear in each of its two arguments. Representing these arguments by two vectors u, v ∈ Rn , the evaluation of the bilinear map F 00 (x) can be written as  i=1,...,m n 2 i X ∂ F (x) j k  y = F 00 (x) (u, v) =  u v . (10) ∂xj ∂xk j,k=1

In contrast to (9), where the tensor F 00 (x) is contracted by a vector v ∈ Rn and a covector z¯ ∈ (Rm )∗ , we now contract F 00 (x) by two vectors u, v ∈ Rn resulting in a vector y ∈ Rm . The bidirectional derivative (10) can be computed by several generalizations of the forward mode, see Griewank et al. (2000). We use the differentiation package ADOL-C, see Griewank et al. (1996), which utilizes operator overloading in C++. The code to be differentiated has to be prepared as follows: First, all relevant variables, which are usually implemented

Fig. 1. Workflow for derivative calculation with ADOL-C using the floating point type double, have to be reassigned as variables of type adouble. Second, we have to classify the input and output variables, see Walther et al. (2003) for details. From this, ADOL-C generates a computational trace of the map F , which is stored in a global serial data structure called tape and can be accessed via an integervalued handle tag. Then, various types of derivatives can be calculated using drivers. The C++ drivers forward and reverse provide an interface for the forward and reverse mode. Moreover, there are several drivers of classical derivatives, which are often used in optimization. For example, the gradient g and the Hessian matrix H of a function accessible via the handle tag can be computed using the following drivers: gradient(tag,n,x[n],h[n]) hessian(tag,n,x[n],H[n][n]) The matrix vector products (7b) and (8) corresponding to directional and weighted derivatives, respectively, are calculated with jac_vec(tag,m,n,x[n],v[n],w[m]) vec_jac(tag,m,n,x[n],zbar[m],xbar[n]) The second order derivative (9) can be evaluated with lagra_hess_vec(tag,m,n,x[n],zbar[m],ybar[n]) These ADOL-C drivers can be called easily from MATLAB using appropriate wrapper functions, see R¨obenack et al. (2011). To explain the usage of ADOL-C we consider a simple function F : R2 → R defined by p z = F (x1 , x2 ) = x1 · sin(x2 ) , which should be evaluated at the point x1 = 4, x2 = 5. The required code modifications to calculate derivatives with ADOL-C are illustrated in Fig. 1. Algorithmic differentiation has been employed in Palis and Palis (2010) for the simulation of mechanical system based on Hamilton’s equations of motion, where the authors used their own implementation of the reverse mode in MATLAB. The simulation of holonomic systems using the Euler-Lagrange equations has been treated in R¨obenack

et al. (2011). We want to extend these results to the simulation of nonholonomic systems using Hamel’s Equations (4) and (6). We assume that we have the computational traces of the ˜ in quasi-velocity coordinates and both the Lagrangian L transformation (2) and its inverse (3). The full gradient ! ˜ ∂L ˜ ∂L , ∂q ∂ω of the Lagrangian can be computed in a single pass of the reverse mode. Since the transformation (2) and its inverse (3) are linear with respect to the velocities and the quasivelocities, respectively, the forward mode can be used for an efficient computation of the matrix entries ∂ω j (q) ˙ ∂ q˙σ (ω) i = ω = $lσ (q) . (q) and ν ∂ q˙ν ∂ω l j In order to calculate the quantities γl,k , cf. Equation (5), we introduce  n  2 j X ∂ ω (q) ˙ j $lσ (q)$kν (q), (11) βl,k = ∂ q˙ν ∂q σ ν,σ=1

k, l = p + 1, . . . , n, j = 1, . . . , n, j j j j such that γl,k = βl,k − βk,l . The skew symmetry γl,k = j −γk,l thus becomes obvious and we only need to comj pute βl,k , k 6= l. The sum (11) represents a bidirectional second order derivative of the transformation (2) having the form (10). The vectors u and v are replaced by ($lσ (q))σ=1,...n and ($kν (q))ν=1,...,n , respectively, with l and k fixed. In ADOL-C, these second order derivatives can be calculated using the function tensor eval, which evaluates derivative tensors in certain directions, see Griewank et al. (1996, 2000). Alternatively, these directional derivatives can be obtained with the driver function forward using higher order univariate Taylor series, see R¨obenack (2008).

Fig. 2. Mobile robot rolling on the horizontal plane. The wheels (light blue) are assumed to be ideal disks. The center of gravity of the main body (dark blue) is labeled by CG. Its distance from the wheels’ axle is l. this body are Jx , Jy , and Jz , while the wheels are assumed to be ideal disks with principal moments of inertia 21 mw r2 and 14 mw r2 , respectively. The distance between the center of gravity of the main body and the wheels’ common axis is denoted by l. In order to manipulate the robot, the torques τ1 and τ2 can be applied between the body and the right and left wheel, respectively. Since slipping is not permitted, the following nonholonomic constraints hold:

4. EXAMPLE We consider a single-axle mobile robot with two independently driven wheels rolling on a horizontal plane under the influence of the gravitational field as illustrated in Figure 2. The dimension of the configuration manifold is n = 5 and the number of constraints is p = 2. This example is an adapted version of the system considered by Franke et al. (2009). The body of the robot as well as each of the wheels are assumed to be rigid. The position of the center of the axle can be expressed in terms of the coordinates x and y. These are the Cartesian coordinates of the horizontal plane. The wheels’ angles with respect to the main body are described by ϕ1 and ϕ2 . Furthermore, the angles ϑ and α denote the rotation of the robot’s body around the vertical axis and, respectively, the wheels’ axis. The constants r and b denote the wheels’ radius and the length of the axle, respectively. Moreover, the y-axis of the global coordinate system is defined to be parallel to the robot’s axle for ϑ = 0. The gravitational field acts in vertical direction with gravity constant g. If α = 0 and ϑ = 0 the principal axes of inertia of the main body are supposed to be parallel to the axes of the global coordinate system. The corresponding principal moments of inertia of

ω 1 (q) ˙ = x˙ sin ϑ − y˙ cos ϑ = 0, r ω 2 (q) ˙ = rα˙ + (ϕ˙ 1 + ϕ˙ 2 ) − (x˙ cos ϑ + y˙ sin ϑ) = 0. 2 To get an invertible transformation, the quasi-velocities ω 1 , ω 2 are accompanied by ω 3 (q) ˙ = α, ˙ ω 4 (q) ˙ = x˙ cos ϑ + y˙ sin ϑ, and r ω 5 (q) ˙ = (ϕ˙ 1 − ϕ˙ 2 ), b where ϑ = r(ϕ1 − ϕ2 )/b. To give an illustrative meaning for these variables, ω1 can be interpreted as the velocity of the axle’s center point in the direction of the wheels’ axis and ω 4 as its (nonholonomic) path velocity. The velocity coordinate ω 2 gives a measure for the slip of the wheels with respect to the axle’s center point. Moreover, ω 3 and ω 5 represent the angular velocities of the main body around the wheels’ axis and the vertical line, respectively. The Lagrangian in terms of the quasi-velocities ω i is then given by

˜ ω) = m0 + 2mw (ω 1 )2 + mw (ω 2 )2 + L(q, 2 2 Jy + m0 l2 3 2 m0 + 3mw 4 2 + (ω ) + (ω ) + 2 2 J5 + (ω 5 )2 + m0 l(ω 3 ω 4 cos α − ω 1 ω 5 sin α)+ 2 − mw ω 2 ω 4 − m0 lg cos α (12) where m0 , mw are the masses of the body and one wheel, respectively. The moment of inertia J5 is defined by 3 1 J5 = (Jx + m0 l2 ) sin2 α + Jz cos2 α + mw b2 + mw r2 . 4 2 Since the quadratic terms of the variables ω 1 and ω 2 vanish in Hamel’s equations, the first two expressions in (12) can be omitted for the derivation of the equations of motion. Finally, the nonholonomic system can be simulated by evaluating the partial derivatives of the Lagrangian in every integration step. In order to illustrate Hamel’s approach, the equations of motion of the mobile robot are stated in the following paragraphs. For the given change of velocity variables all j 1 1 coefficients γl,k vanish except γ2,3 and γ3,2 , which take the values −1 and +1, respectively. The generalized forces corresponding to the new velocities are given by f˜3 = −(τ1 + τ2 ), f˜4 = (τ1 + τ2 )/r, and f˜5 = b(τ1 − τ2 )/(2r). Hence, the equations of motion can be written as f˜3 = (Jy + m0 l2 )ω˙ 3 + m0 l cos αω˙ 4 + − 12 (Jx − Jz + m0 l2 ) sin(2α)(ω 5 )2 − m0 lg sin α, f˜4 = m0 l cos αω˙ 3 + (m0 + 3mw )ω˙ 4 +   + m0 l sin α (ω 5 )2 − (ω 3 )2 , and
f˜5 = J5 ω˙ 5 + Jx − Jz + m0 l2 sin(2α)ω 3 ω 5 + − m0 l sin αω 4 ω 5 , This system can be easily solved for (ω˙ 3 , ω˙ 4 , ω˙ 5 ) yielding a system of ordinary differential equations. These equations are completed by a system of kinematic differential equations   ϕ1 − ϕ2 ω4 x˙ = cos r b   ϕ1 − ϕ2 y˙ = sin r ω4 b α˙ = ω 3 1 b ϕ˙ 1 = −ω 3 + ω 4 + ω 5 r 2r 1 4 b 3 ϕ˙ 2 = −ω + ω − ω 5 r 2r which is the inverse of the velocity transformation taking the nonholonomic constraints into account, compare Equation (6). Thus, we have generated a system of 2n − p = 8 ordinary differential equations describing the motion of the robot.

5. SUMMARY AND OUTLOOK In this paper, the approach of Hamel (1904) has been employed to derive the equations of motion of nonholonomic mechanical systems. Compared to the Lagrange multiplier approach, this leads to a system of ordinary differential equations instead of differential algebraic equations. Moreover, the dimension of this system is lower. Subsequently, the technique of algorithmic differentiation has been used to generate a numerical simulation model of the nonholonomic system. Using this technique, it is only necessary to state the total energy of the system, the nonconservative forces, and the nonholonomic constraints to simulate the system. In this respect, our contribution extends the results of Palis and Palis (2010) and R¨obenack et al. (2011). During the last years, algorithmic differentiation was also used for the direct computation of nonlinear control laws based on Lie derivatives, see R¨obenack (2005a,b). Further research could be devoted to the computation of nonlinear controllers for mechanical systems with nonholonomic constraints. ACKNOWLEDGEMENTS The authors gratefully acknowledge support from Dr. Jan Winkler, TU Dresden, who implemented the interface between ADOL-C and MATLAB. The first author is indebted to Prof. Joachim Rudolph, Saarland University, for bringing his attention to nonholonomic systems and for some illuminating discussions. Furthermore, the third author would like to thank his colleague, Prof. Thomas Sattel, TU Ilmenau, for an interesting discussion on nonholonomic systems during the 56th International Scientific Colloquium, which contributed to the paper. REFERENCES Abraham, R. and Marsden, J. (1978). Foundations of Mechanics. Addison-Wesley Pub. Co. Brenan, K.E., Campbell, S.L., and Petzold, L.R. (1996). Numerical Solution of Initial-Value Problems in Differential-Algebraic Equations. SIAM, Philadelphia, 2nd edition. Corliss, G., Faure, C., Griewank, A., Hasco¨et, L., and Naumann, U. (eds.) (2002). Automatic Differentiation: From Simulation to Optimization. Springer-Verlag, New York. Franke, M., Rudolph, J., and Woittennek, F. (2009). Motion planning and feedback control of a planar robotic unicycle model. In Proc. Int. Conf. on Methods and Models in Automation and Robotics, volume 14. Goldstein, H. (1980). Classical Mechanics. AddisonWesley series in physics. Addison-Wesley Pub. Co. Griepentrog, E. and M¨arz, R. (1986). DifferentialAlgebraic Equations and Their Numerical Treatment, volume 88 of Teubner-Texte zur Mathematik. Teubner Verlagsgesellschaft, Leipzig. Griewank, A., Juedes, D., and Utke, J. (1996). ADOL-C: A package for automatic differentiation of algorithms written in C/C++. ACM Trans. Math. Software, 22, 131–167. Griewank, A., Utke, J., and Walther, A. (2000). Evaluating higher derivative tensors by forward propagation of univariate Taylor series. Mathematics of Computation, 69, 1117–1130.

Griewank, A. and Walther, A. (2008). Evaluating Derivatives: Principles and Techniques of Algorithmic Differentiation. SIAM, 2nd edition. Hairer, E., Lubich, C., and Roche, M. (1989). The Numerical Solution of Differential-Algebraic Systems by Runge-Kutta Methods, volume 1409 of Lecture Notes in Mathematics. Springer-Verlag, Berlin. ¨ Hamel, G. (1904). Uber die virtuellen Verschiebungen in der Mechanik. Mathematische Annalen, 59, 416–434. Hamel, G. (1949). Theoretische Mechanik. Springer. Kunkel, P. and Mehrmann, V. (2006). DifferentialAlgebraic Equations. EMS Publishing House, Z¨ urich. Lee, J.M. (2002). Introduction to smooth manifolds. Springer. Neimark, J. and Fufaev, N. (1972). Dynamics of Nonholonomic Systems. Translations of mathematical monographs. American Mathematical Society. Palis, S. and Palis, F. (2010). Mechanical system simulation via automatic differentiation. In Proc. of Advanced Problems of Mechanics. R¨obenack, K. (2005a). Computation of Lie derivatives of tensor fields required for nonlinear controller and observer design employing automatic differentiation. Proc. in Applied Mathematics and Mechanics, 5(1), 181– 184. R¨obenack, K. (2005b). Regler- und Beobachterentwurf f¨ ur nichtlineare Systeme mit Hilfe des Automatischen Differenzierens. Shaker Verlag, Aachen. R¨obenack, K. (2008). Computation of multiple Lie derivatives by algorithmic differentiation. J. of Computational and Applied Mathematics, 213(2), 454–464. R¨obenack, K., Bausch, T., and Uhlig, A. (2001). Struktureller Index von Deskriptorvariablen. In K. Panreck and F. D¨orrscheidt (eds.), Simulationstechnik, 15. Symposium in Paderborn, ASIM 2001, Frontiers in Simulation, 569–574. SCS – The Society for Modeling and Simulation International. R¨obenack, K., Winkler, J., and Knoll, C. (2011). Direct simulation of mechanical control systems using algorithmic differentiation. In Proc. 56th International Scientific Colloquium. Ilmenau. Saunders, D.J. (1989). The Geometry of Jet Bundles. Cambridge University Press. Walther, A., Griewank, A., and Vogel, O. (2003). ADOLC: Automatic differentiation using operator overloading in C++. Proc. in Applied Mathematics and Mechanics, 2(1), 41–44.