The 2nd Joint International Conference on Multibody System Dynamics May 29-June 1, 2012, Stuttgart, Germany

Computation of independent sensitivities using Maggi’s formulation María D. Gutiérrez-López*, Alfonso Callejo*, Javier García de Jalón*# *Computational Mechanics Group INSIA, Universidad Politécnica de Madrid Ctra. Valencia km 7, 28031, Madrid, Spain {md.gutierrez, a.callejo}@upm.es

#

Department of Applied Mathematics ETSII, Universidad Politécnica de Madrid José Gutiérrez Abascal 2, 28006, Madrid, Spain [email protected]

ABSTRACT The utility of sensitivity analyses in mechanical engineering for various purposes is very well-known. In the field of multibody dynamics, sensitivity methods tend to be very complex and numerical differentiation is often used instead. In this work, a method based on the differentiation of a recursive Maggi’s formulation is presented. Efficiency and accuracy are assessed by implementing and comparing three different types of differentiation techniques (numerical, direct and automatic) with two numerical examples (a 3D double-pendulum and a 3D four-bar mechanism). Results show how the integration of independent sensitivities through Maggi’s scheme is accurate and stable. 1. INTRODUCTION In many multibody applications such as mechanism design, system optimization, parameter identification and optimal control, it is very useful to know how the model inputs affect the model outputs. Sensitivity analysis quantifies these dependencies by computing the variation of the response with respect to the design parameters. This calculation requires the derivatives of the positions, velocities, accelerations, and Lagrange multipliers, with respect to the design parameters. Existing methods for the computation of these sensitivities are often complicated, inefficient and not valid for large systems [1]. Several methods are usually considered for the computation of sensitivities: numerical differentiation, direct differentiation [7], [11], [12], [14], adjoint variable method [2], [3], [5], [6], [13] and automatic differentiation [4]. This work presents an original algorithm for the computation of sensitivities through direct differentiation (DD) and automatic differentiation (AD). In many of the approaches described in the literature (for instance, [14], [15]), the Lagrange’s equations of the first kind are differentiated with respect to the design parameters to obtain the sensitivities. However, this approach leads to a set of differential-algebraic sensitivity equations of high dimension, which has made the direct differentiation method not very popular in the past [14]. In order to reduce the computational effort associated with direct differentiation methods, in this article the motion differential equations are written using a minimal set of joint coordinates, following a recursive double-step Maggi’s approach [9]. The formulation is efficient, and, at the same time, the underlying recursivity is easily obtained through the so-called path matrix and the use of natural coordinates for the initial geometry. The theoretical approach is completed with the analysis of two academic examples: a 3D double pendulum and a 3D four-bar mechanism. The results obtained by direct differentiation are verified using numerical differentiation. Finally, a comparison with an automatic differentiation method is carried out. 2. DYNAMIC ANALYSIS THROUGH MAGGI’S EQUATIONS A robust formulation of the motion differential equations is crucial for a good sensitivity analysis. In this section, an efficient version of Maggi’s formulation is described. Starting from the descriptor form of the motion differential equations together with the kinematic constraints, the index-3 system of differential algebraic equations (DAEs) is transformed into a system of state-space ordinary differential equations (ODEs) by means of the coordinate partitioning method [16]-[19]. Then, the generalized velocities and accelerations are integrated over time by avoiding the solution of the finite displacement problem. Terms such as generalized positions, velocities and accelerations, and inertia terms, are computed recursively.

More details on this formulation can be found at [9]. Overall, in spite of the numerical drift associated with the non-enforcement of position constraints and the need to re-compute the basis of the Jacobian null-space every time step, the formulation shows great efficiency and stability. 2.1. Recursive Index-3 motion differential equations In order to apply recursion techniques, the system is considered as a tree-configured multibody system. In the case of closed-chain systems, the joints and rods closing the kinematic chains and are conveniently enforced through constraint equations. Firstly, Cartesian coordinates are used to define the velocity and acceleration of bodies:  s i  Zi    ,  ωi 

  i   si  Z i  ω

(1)

where s i and si are respectively the velocity and acceleration of the point attached to body i that instantaneously coincides with the origin of the inertial reference frame. In this way, all bodies share the same reference point coordinates, which brings important advantages. The recursive expression of the Cartesian velocities and accelerations of body i in terms of those of body i  1 are: Zi  Zi 1  bi zi ,

Z i  Z i 1  bi  zi  d i

(2)

Note the lack of transformation matrices in the previous equations. Vector bi represents the velocity of the point of body i that coincides with the origin of the global reference frame when zi  1 and z j  0, j  i . Vectors bi and di depend on the type of joint between i and i  1 . Only revolute and prismatic joints (and combinations thereof) are considered. Vectors ZT  {Z1T , ZT2 ,..., ZTn } and  T  {Z 1T , Z  T2 ,..., Z Tn } group the system velocities and accelerations, being n the number of moving bodZ ies. The virtual power of the inertia and external forces of the open-chain system can be expressed as: n

 Zi T  Mi Z i  Qi   0

(3)

i 1

M  diag(M1 , M 2 ,..., M n ),  mi I 3 Mi =   mi g i

mi g i  , J i  mi g i g i 

QT  {Q1T , QT2 ,..., QTn }

(4)

 iω  i gi Fi  mi ω   Qi         m   g F ω J ω g ω ω g  i i i i i i i i i i 

(5)

where M i   66 and Qi   61 are, respectively, the inertia matrix and the vector of external and velocity-dependent inertia forces acting on body i , mi is the mass, J i  33 is the inertia tensor, gi is the position of the center of gravity (COG), ωi is the angular velocity and Fi is the applied force. The upper bar indicates reference to the origin of the global reference frame. The upper tilde transforms the vector  . into the associated skew-symmetric matrix, such that, for generic vectors α and β , α  β  αβ A velocity transformation R   6nn between Cartesian and relative velocities z is now introduced. Bodies (and their corresponding input joints) are numbered from the leaves to the root of the spanning tree. The j-th column of matrix R are the Cartesian velocities of all the bodies that are upwards of joint j when a unit relative velocity is introduced in j keeping null the remaining relative velocities; due to the fact that the origin of the global reference frame is the reference point for all bodies, these Cartesian velocities happen to be bi for all bodies, according to Equation (2). Z  Rz  T diag  b1 , b 2 ,..., b n  z  TR d z

(6)

  TR d   d z Z z  TR

(7)

The connectivity of the mechanism has been defined through an upper triangular path matrix T   nn . Element Tij is 1 if body i is upwards of joint j and 0 otherwise. Substituting (6) and (7) into Equation (3), multiplying by RTd TT and knowing that the relative open-chain virtual velocities are independent, the dependent virtual velocities Zi are eliminated and a new set of differential equations is obtained:  d z ) RTd M  R d  z  RTd (Q  M  R

(8)

where M   TT MT and Q  TT Q are the accumulated inertia matrix and force vector respectively. Finally, the closure-of-the-loop constraint equations Φ   m1 have to be added; three types are considered: revolute joint, prismatic joint and rod (slender body with two spherical joints and a negligible moment of inertia around the direction of the axis). For the sake of brevity, only the revolute joint constraint is detailed. Constraints are defined in terms of natural coordinates (i.e., Cartesian coordinates of points and unit vectors) attached to the bodies, and then analytically transformed (when necessary) into relative coordinate expressions. The five independent constraints of a revolute joint between bodies j and k, in terms of the joint position r j ,k and the parallel unit vectors u j ,k , can be extracted from:  r j  rk  Φ   u j  uk 

(9)

The Jacobian matrix of the constraint equations can be obtained by using the chain rule of differentiation: Φ z  Φr

j

r j r j rk rk  Φr  Φr  Φr k j k z z z z

(10)

Closure-of-the-loop constraints and rod constraints are grouped in a vector of system constraint equations Φ   m1 , whose Jacobian matrix is Φ z   mn . Further details can be found at [9].The double differentiation of Φ(z, t ) with respect to time yields the velocity-level and acceleration-level constraints:  z z  Φ t Φ(z, t )  0, Φ z z  Φt , Φ z  z  Φ

(11)

Equations (8) and (11) constitute a set of recursive index-3 DAEs defining the motion of the multibody system. They can be solved efficiently using Maggi’s approach, as explained in the following sub-section. 2.2. Coordinate partitioning and Maggi’s approach

The numerical solution of Equations (8) and (11) can be facilitated by eliminating the Lagrange multipliers using a basis of the Jacobian matrix null space, which transforms the system of DAEs into a system of ODEs. Using the generalized coordinate partitioning method [16]-[19] based on Gaussian elimination with full pivoting, the relative coordinates z are partitioned into dependent (z d   m1 ) and independent (z i   f 1 ) sets, being f  n  m the number of degrees of freedom; that way, the dependent velocities can be computed from the independent ones as:  Φ dz

 z d  Φiz   i   0  z d  (Φ dz ) 1 Φiz z i  z 

(12)

where matrix Φ dz is assumed to be invertible. Using this expression, we can numerically obtain a transformation matrix R z   n f relating dependent and independent relative velocities (and accelerations by differentiating with respect to time): 1  z d     Φ dz  Φiz  i  z  R z z i z   i      z   If

(13)

 z z i  z  R z  zi  R

(14)

Introducing (14) into (8) and pre-multiplying by RTz , the final set of recursive motion differential equations (with one ODE for each independent acceleration) is obtained:  d z  R d R  z z i ) RTz RTd M  R d R z  z i  RTz RTd Q  RTz RTd M  (R

(15)

For the sake of convenience, some terms are grouped and the dependencies are made explicit: ˆ (q, q , t )  Pˆ (q, q , t ) ˆ (q) M z i (t )  Q

(16)

Equation (16) can be integrated over time using an appropriate time-integration scheme. Maggi’s formulation implies the configuration of the state vector yT  {zT , z iT } instead of the somewhat more typical yT  {z iT , z iT } . The latter requires the solution of a nonlinear system of equations (the finite displace-

ment problem) in order to compute the dependent positions z from the independent ones z i every time step, whereas that process is avoided by the former. On the other hand, the former suffers from a numerical drift coming from the fact that the position-level constraints are not enforced. 3. DIRECT DIFFERENTIATION OF DYNAMIC EQUATIONS

In this section, the equations of motion (16) coming from Maggi’s formulation are differentiated to obtain an expression that allows computing sensitivities z p , z p and zp . Vector p contains the parameters of the system, such as inertias, masses, spring stiffnesses, damping coefficients, etc. A single generic parameter p is considered for the purpose of theoretical development. Differentiating Equation (16) with respect to p and applying the chain rule, the following expression is obtained: ˆ p Q ˆ qq p  Q ˆ q q p  Pˆ p  Pˆ q q p  Pˆ q q p  M ˆ ip  Q ˆ p z i  (Mz ˆ i )q q p Mz

(17)

where subscripts denote partial differentiation and the notation of a bar under a term (_) indicates that it is treated as constant in the partial differentiation. Equation (17) is a set of ODEs that can (and must [14]) be integrated together with the set of ODEs in Equation (16) for the independent sensitivities zip , as outlined by [8]. For a better presentation of the method, Equation (17) is going to be simplified by introducing total derivatives with respect to design parameter p: ˆ ˆ ˆ ˆ ip  dQ  dP  dM  Mz zi dp dp dp

(18)

where the new terms introduced in Equation (18) have the following expressions: ˆ dQ ˆ p Q ˆ qq p  Q ˆ q q p , Q dp

ˆ dM ˆ p  ˆ i )q q p ,  zi  M z i  (Mz dp

dPˆ ˆ  Pp  Pˆ q q p  Pˆ q q p dp

(19)

In the sequel, a procedure for obtaining the total derivatives of the above equation is presented. In order to compute those derivatives, the sensitivities of Cartesian positions, velocities and accelerations must be calculated first. For brevity in the exposition, the influence of rods is not detailed. 3.1. Sensitivities of Cartesian positions For open-chain systems (and for closed-loop systems, after removing certain joints and/or rods), the sensitivities of Cartesian positions with respect to the design parameters can be computed recursively from the relative coordinates and their sensitivities. The sensitivity with respect to parameter p of the input point rk of body i is obtained from the sensitivity of the output point rj of body i–1, depending on the type of joint between both bodies: rkp  r jp (revolute joint)

(20)

rkp  r jp  u jp zi  u j zip (prismatic joint)

(21)

The sensitivity with respect to the design parameter p of any point rs of body i, can be obtained by using the following equation: rsp  rkp 

dA i r  Ai rsp dp s

(22)

where rs is a vector with the coordinates of point s expressed in a local reference frame attached to the moving body and Ai is the matrix needed to transform vector rs to the global reference frame. The sensitivities of rs with respect to the design parameters are easy to define. They are generally null, except when the design parameter p is one of the coordinates of point s in the local reference frame. The expression required to compute the derivative of Ai with respect to parameter p depends on the type of joint i. The derivatives of the unit vectors can be computed in a similar way by using the following expression:

uip 

dA i u  A i uip dp i

(23)

3.2. Sensitivities of Cartesian velocities and accelerations The recursive expressions for the computation of the sensitivities of Cartesian velocities and accelerations with respect to design parameter p are the following: dZi dZi 1 dbi   z  bi z ip , dp dp dp i

dZ i dZ i 1 dbi dd    z  bi  z ip  i dp dp dp i dp

(24)

where the derivatives of bi and di depend on the type of joint i between the contiguous bodies. In the case of closed-loop systems, the sensitivities of the dependent relative velocities must be firstly obtained for the computation of the sensitivities of the Cartesian velocities. For open-chain systems all relative velocities are independent and this step can be skipped. Taking derivatives of Equation (13) with respect to p, one can obtain the expression that relates the sensitivities of the relative velocities to the independent relative velocities and their sensitivities:  z dp  dR z i z p   i   z  R z z ip dp  z p 

(25)

where the total derivative of the transformation matrix R z with respect to p can be computed as:  dR z   Φ dz   dp 



1

dΦ dz  Φdz dp



1

Φiz   Φ dz



1

03

dΦiz  dp   

(26)

The derivative of the Jacobian matrix of constraint equations with respect to the design parameter is obtained from the sensitivities of the Cartesian positions. For sake of brevity, the expressions needed for its computation are not included here. Taking into account the sensitivities of the dependent relative velocities calculated according to Equation (25), the sensitivities of the Cartesian velocities can be computed using Equation (24). Considering Equations (7), (13) and (14), Cartesian accelerations can be written as follows:   TR d R z   d Rz  Rd R  z  z i Z zi  T  R

(27)

 d Rz  Rd R  z  appearing in Equation (15) can be computed Therefore, the derivative of the product T  R  as the sensitivity of the Cartesian accelerations Z when the relative independent accelerations zi are null.

For closed-loop systems, the sensitivity of the dependent relative accelerations corresponding to null  d Rz  Rd R  z  as explained before. Dependent zi must be computed to obtain the derivative of T  R relative accelerations are related to relative velocities by the following equation when the independent relative accelerations are null:  dz Φ zd  z d   Φ

d  iz   z   0 Φ  z i  

(28)

Taking derivatives of the previous equation with respect to the design parameter p, the following expression is obtained:  z dp   Φ zd



 zd  dΦ    dp

1 

 iz   z d  dΦ d    Φ dp   z i   z

d d   iz   z p   dΦ z  Φ zd   z i dp  p 

(29)

This equation allows computing the sensitivities of the dependent relative accelerations when the independent accelerations are null.

3.3. Derivatives of the accumulated inertia matrix and force vectors After having computed the sensitivities of the Cartesian positions and the Cartesian velocities, the total derivatives of the inertia matrix M i of body i and of the vector of external forces Qi with respect to a parameter p can be computed by differentiating the matrices of Equation (5).  d z  R d R  z z i ) appearing in Equation (15) will be denoted as P  from now on. SimiThe product M  (R larly, vector Pi is defined as:  d z i  R d R  z z i ) Pi  M i T(R

(30)

The derivative of vector Pi can be computed as: i dPi dM i  dZ Zi  M i  dp dp dp

(31)

 i are the Cartesian accelerations when the independent relative accelerations are null. where Z The total derivative of the accumulated inertia matrix M i with respect to a parameter p can be easily computed by adding the derivatives of the inertia matrices of all elements upwards of joint i. The derivatives of the accumulated external forces Qi and the accumulated inertia forces Pi are computed in a similar way. 3.4. Miscellaneous derivatives Once the total derivatives of the accumulated inertia matrix and external forces with respect to parameter p have been computed, the derivatives of products RTd M  R d and RTd Q Σ can be obtained. The total derivative of product RTd M  R d with respect to p is computed by using the path matrix T to assemble a set of derivatives. Provided that body j is under body i, the element located in row i and column j of that matrix (and vice versa) is computed as: d  bTi M i b j dp

  dbTi M b dp

i

j

 bTi

db j dM i b j  bTi M i dp dp

(32)

The rest of elements of that matrix are null. Similarly, element i of the derivative of RTd Q Σ with respect to parameter p is obtained through the following expression: d  bTi Qi  dp



dbTi  dQi Qi  bTi dp dp

(33)

The element i of the total derivative of product RTd P  with respect to parameter p is computed as: d  bTi Pi  dp



dbTi  dPi Pi  bTi dp dp

(34)

Finally, the total derivatives of products RTz  RTd M  R d  R z , RTz  RTd Q Σ  and RTz  RTd P   can be obtained through the following equations: ˆ d  RTd M  R d  dM dRTz dR  T T  R M R R  R R z  RTz  RTd M  R d  z  d d  z z dp dp dp dp

(35)

ˆ dRTz d  RTd Q Σ  dQ  RTd Q Σ   RTz  dp dp dp

(36)

d  RTd P   dPˆ dRTz T  T  R P  R  d  z dp dp dp

(37)

These derivatives are used to form the set of ODEs given by Equation (18), from which the sensitivities of the independent accelerations can be computed. y x3 x4

x4 x3 x4 1

x3 x4 ln x3

x1 x2

x3

x4 sin x2

 x1 x22

x21 x1

cos x2 x2

Figure 1. Computational graph of (x1/x2)^sin(x2). 4. AUTOMATIC DIFFERENTIATION OF DYNAMIC EQUATIONS Automatic differentiation (AD) is an alternative computational-mathematical technique for the differentiation of computer functions. It is based on the decomposition of computer functions in elementary arithmetic operations (addition, subtraction, product, division) and calls to library functions (sine, cosine, exponential, etc.), and on the systematic application of the chain rule of differentiation [10]. With it, any mathematical code, no matter how complex, can be differentiated automatically. In this section, AD is used to differentiate Equation (16) and obtain sensitivities z p , z p and zp .

4.1. Computational graphs AD is based on the computational graph representation of the function under differentiation. Computational graphs are a type of directed acyclic graphs (DAGs). They represent the sequence of operations performed by the computer in the evaluation of a function. The result of each elementary operation is stored in a separate vertex, and the flow of variables is displayed as directed edges between vertices. Consider, as an example, the following scalar function: x sin x2 y  f ( x1 , x2 )   1   x2 

(38)

This function can be rewritten so that the result of each arithmetic operation or library function is stored in a new intermediate variable x j , with index j greater than the one of the variables it depends on: x3 

x1 , x4  sin x2 x2 y  x3 x4

(39) (40)

where x1 and x2 are the independent variables, x3 and x4 are intermediate variables and y  x5 is the dependent variable. Figure 1 shows the computational graph of the previous sample function. Bottom vertices are associated with independent variables, the top vertex with the dependent variable and middle vertices with intermediate variables. The edges linking the vertices represent the direction of the data flow. The value on each edge is the partial derivative of the target vertex with respect to the source edge.

4.2. Modes of differentiation In the example, the derivative of the dependent variable y with respect to x1 can be computed by applying the chain rule of differentiation systematically from the bottom vertices to the top vertex, or vice versa. Accordingly, there are two modes of AD: forward and reverse. Only the former is explained here. In the forward mode, the chain rule is evaluated from the independent to the dependent variables, and the com-

putational load is proportional to the number of independent variables. The derivative of a vertex (i.e., of the associated variable) is the addition of the contributions of the edges that arrive to that vertex: x1 

dx1 1 dx1

x3 

dx3 x3 dx1 x3 dx2    ( x21 ) x1  ( x1 x22 ) x2 dx1 x1 dx1 x2 dx1

(41)

x2 

dx2 0 dx1

x4 

dx4 x4 dx2   (cos x2 ) x2 dx1 x2 dx1

(42)

y 

dy y dx3 y dx4    ( x4 x3 x4 1 ) x3  ( x3 x4 ln x3 ) x4 dx1 x3 dx1 x4 dx1

(43)

Each edge contributes with the total derivative of the source vertex multiplied by the partial derivative associated with the edge. Thus, the derivatives can be propagated together with the evaluation of the intermediate variables, as Equation (41) shows. Operator d (.) d (.) denotes total derivatives, and operator  (.)  (.) denotes partial derivatives.

4.3. Implementation There are two main ways of implementing AD: using operator overloading tools and source transformation tools. Both have the same theoretical foundation and AD modes. Source transformation methods generate a new code with the explicit expressions of the derivatives. This new code has to be called by the program where the derivatives are needed and conveniently compiled. On the other hand, operator overloading tools transform each variable of the function into a structure whose member variables are the value of the variable and the value of the derivatives with respect to the independent variables. Then, arithmetic operations and library functions are overloaded so that they handle those structures and compute both the result of the operation and the corresponding derivative, following the rules of differentiation. The original program remains almost unchanged, because all those substitutions of variables, operators and functions are carried out at compile time. This way, the implementation of AD is very versatile and valid for nearly all functions. However, the execution is rather slow because of the indirect costs of overloading operators and functions (the cost of indirect addressing and the runtime overhead). In the present work, the C-language operator overloading tool ADOL-C has been implemented [17]. The efficiency results may be considered as a limit inferior of the efficiency that can be obtained through source transformation tools like TAPENADE or ADIC. The function to be differentiated with respect to the vector of parameters p, as detailed in Sections 2 and 3, is:  zT y  f (y , p), y   iT  z

  z T  f ,     ˆ (z, z , p, t )  Pˆ (z, z , p, t ))  ˆ (z, p) 1 (Q  M 

(44)

Regarding the mode of differentiation, as in general real-life models the number of inputs of f (the number of parameters) is a lot smaller than the number of outputs (2n), the forward mode is used. The first time function f is executed, the computational graph (which is nothing but the sequence of basic operations and calls to library functions) is recorded by the AD tool. Later on, this registry can be used to execute the function again or to compute the derivatives of the function with respect to the inputs; in this case, those derivatives are equivalent to the Jacobian matrix of f with respect to p, referred to as g:  zTp  y p  g (y p , p), y p   iT   z p 

(45)

Motion differential equations (44) and sensitivity equations (45) can therefore be integrated jointly using the same time integration scheme and without having to differentiate the equations of motion by hand.

5. NUMERICAL EXAMPLES In this section, two examples are used to validate the presented method: a 3D double-pendulum and a 3D four-bar mechanism. The sensitivities of those mechanisms are computed by using direct (DD), numerical

(ND) and automatic differentiation (AD), in order to assess efficiency and accuracy. As AD is the only one implemented in C language (the rest have been programmed in MATLAB), the CPU times have to be compared with caution. In the ND approach the size of the perturbation is 10-6 and forward differences are used. A 4th-order Runge-Kutta integrator with a time-step of 10 ms has been used for both examples. Gravity acts in the –Z-direction. All simulations have been run on an Intel® Core™ i7 at 2.93 GHz with 6 GB and Windows 7 Professional. A 5-second simulation has been run for each analysis. Z

Z

(a)

(b) Y c, k

Joint 1 Bar 1

Joint 4 Y

Joint 1

z4

c, k Ix , m

Ix , m

Bar 1

Joint 2

Bar 3

z2

Bar 2 U

Bar 2

S

Joint 2

Joint 3

Figure 2. Double-pendulum (a) and four-bar mechanism (b). 5.1. Double-pendulum The double-pendulum has two bodies, two parallel revolute joints and two degrees of freedom. Figure 2(a) shows the initial configuration. Bar 1 has an initial angular velocity of –1 rad/s. Both bars have a length of 1 m and a mass of 1 kg. A spring-damper set is attached to joint 1, and has the following properties: k=1 Nm, c=1 Nsm. Position, velocity and acceleration sensitivities with respect to the following parameters are computed by using DD, ND and AD:    

Mass of bar 1 (m). Inertia moment of bar 1 (Ix). Damping of joint 1 (c). Stiffness of joint 1 (k).

Figure 3 shows the sensitivity of z2 with respect to m. A good agreement between the results obtained by ND, DD and AD is achieved. Figure 4 shows the relative error as error(AD)  (DD AD) DD 1 and error(ND)  (DD ND) DD 1 in the case of z2 m . The error of ND, although of relatively low importance, is the highest of the three errors. -4

0.8

1.2

ND DD AD

0.6

x 10

ND error AD error

1

0.4 0.8

relative error (m)

z 2m

0.2 0 -0.2 -0.4

0.4

0.2

-0.6 -0.8

0.6

0

1

2

3

4

5

t [s]

Figure 3. Sensitivity of z2 with respect to m.

0

0

1

2

3

4

t [s]

Figure 4. Relative error of z2m.

5

Table 1 shows CPU times of the dynamic analysis (first row), the sensitivity analysis (second row) and both (third row) with each of the three methods. Last row contains the weight of the sensitivity computation relative to the total time. Four analysis, each of them adding a parameter to the computation, have been run. Comparing ND and DD (both implemented in MATLAB), the efficiency of DD is higher than the efficiency of ND when the number of parameters is higher. However, when only one or two parameters are analysed, ND times are shorter.

Table 1. CPU times of the double-pendulum. CPU time [s] Dynamics Sensitivity Total Sens./Total

ND 3.182 2.871 6.050 47%

m DD 3.182 3.900 7.080 55%

AD 0.004 0.133 0.137 97%

ND 2.527 4.868 7.400 66%

m, Ix DD 2.527 5.335 7.860 68%

AD 0.004 0.142 0.146 97%

ND 3.573 8.642 12.22 71%

m, Ix, c DD AD 3.573 0.004 7.004 0.152 10.58 0.156 66% 97%

ND 2.605 11.71 14.32 82%

m, Ix, c, k DD 2.605 10.09 12.70 79%

AD 0.004 0.254 0.258 98%

5.2. Four-bar mechanism The four-bar mechanism is a closed-loop system with three moving bodies, four joints (revolute, universal, spherical and revolute) and one degree of freedom. Figure 2(b) shows the configuration at the initial position. Bar 1 has an initial angular velocity of –1 rad/s. The system is converted into an open-chain system by removing the spherical joint between bars two and three. All bars have the same properties (length, mass and inertia moments) as the bars of the previous example. A spring-damper set is attached to joint 1, with the same properties as in the double-pendulum. Again, the following parameters are computed using ND, DD and AD:    

Mass of bar 1 (m). Inertia moment of bar 1 (Ix). Damping of joint 1 (c). Stiffness of joint 1 (k).

Figure 5 shows the sensitivity of z4 with respect to m. A good agreement between the results obtained by ND, DD and AD is achieved. Figure 6 shows the relative error as error(AD)  (DD AD) DD 1 and error(ND)  (DD ND) DD 1 . The error of ND, although of relatively low importance, is the highest of the three errors. Table 2 shows the CPU times the same way the previous example did.

Table 2. CPU times of the four-bar mechanism. CPU time [s] Dynamics Sensitivity Total Sens./Total

ND 5.632 5.678 11.31 50%

m DD 5.632 8.361 13.99 60%

AD 0.010 0.233 0.243 96%

ND 6.021 11.72 17.74 66%

m, Ix DD 6.021 11.00 17.02 65%

AD 0.010 0.244 0.254 96%

m, Ix, c ND DD 6.162 6.162 17.08 15.02 23.24 21.19 73% 71%

AD 0.010 0.380 0.390 98%

ND 6.133 20.86 26.99 77%

m, Ix, c, k DD 6.133 16.67 22.80 73%

AD 0.010 0.264 0.274 96%

6. CONCLUSIONS The coordinate partitioning method has been used to formulate sensitivity equations. That way, only independent sensitivities are integrated, and Maggi’s strategy is used for both the integration of the motion differential equations and the sensitivity equations. Three methods of differentiation (numerical, direct and automatic) and two numerical examples (a 3D double-pendulum and a 3D four-bar mechanism) have been implemented to explore the effect of the differentiation method on accuracy and efficiency. Direct and automatic differentiation yield accurate sensitivities, whereas numerical differentiation, as expected, entails an error difficult to estimate and control. Direct differentiation is more efficient than numerical differentiation when the number of parameters is large, and the development time with automatic differentiation is far shorter than with direct differentiation. Finally, the effect of the number of parameters on the computation time is smallest with automatic differentiation. In short, a step has been taken towards a more efficient, robust and simple computation of design sensitivities in the field of multibody dynamics.

-5

6

0.2 ND DD AD

0.15

x 10

ND error AD error

5

0.1

4 relative error (m)

z 4m

0.05 0 -0.05 -0.1

2

1

-0.15 -0.2

3

0

1

2

3

4

5

t [s]

Figure 5. Sensitivity of z4 with respect to m.

0

0

1

2

3

4

5

t [s]

Figure 6. Relative error of z4m.

ACKNOWLEDGEMENTS The authors acknowledge the support of the Ministry of Science and Innovation of Spain under Research Project TRA2009-14513-C02-01 (OPTIVIRTEST), and of the Education Department of the Government of Navarre, Spain.

REFERENCES [1] Banerjee, J. M.; McPhee, J.: Graph-Theoretic Sensitivity Analysis of Multibody Systems. In J.-C. Samin and P. Fisette (Eds.) Multibody Dynamics 2011 ECCOMAS Thematic Conference, 2011. [2] Bestle, D.; Eberhard, P.: Analyzing and Optimizing Multibody Systems. Mechanics of Structures & Machines, Vol 20 (1), pp. 67-92, 1992. [3] Bestle, D.; Seybold, J.: Sensitivity Analysis of Constrained Multibody Systems. Archive of Applied Mechanics, Vol. 62, pp. 181-190, 1992. [4] Bischof, C. H.: On the automatic differentiation of computer programs and an application to multibody systems. Proc. IUTAM Symposium on Optimization of Mechanical Systems, Dordrecht : Kwer, pp. 41-48, 1996. [5] Cao, Y.; Li, S.; Petzold, L.: Adjoint sensitivity analysis for differential-algebraic equations : algorithm and software. Journal of Computational and Applied Mathematics, Vol. 149, pp. 171-191, 2002. [6] Cao, Y.; Li, S.; Petzold, L.; Serban, R.: Adjoint sensitivity analysis for differential-algebraic equations : the adjoint DAE system and its numerical solution. SIAM Journal on Scientific Computing. Vol. 24, No. 3, pp. 1076-1089, 2003. [7] Chang, C.O.; Nikravesh P.E.: Optimal Design of Mechanical Systems with Constraint Violation Stabilization Method. Journal of Mechanisms, Transmissions, and Automation in Design, Vol. 107, No. 4, pp. 493–498, 1985. [8] García de Jalón, J.; Bayo, E.: Kinematic and Dynamic Simulation of Multi-Body Systems – The Real-Time Challenge. Springer-Verlag, New York, p. 365, 1994. [9] García de Jalón, J.; Callejo, A.; Hidalgo, A.F.: Efficient solution of Maggi’s equations. To be published in Journal of Computational and Nonlinear Dynamics, DOI: 10.1115/1.4005238, 2011. [10] Griewank, A.: Evaluating derivatives - Principles and Techniques of Algorithmic Differentiation. SIAM, Frontiers in Applied Mathematics, 2000.

[11] Haug, E. J.: Design sensitivity analysis of dynamic systems. C.A. Mota-Soares, Ed., ComputerAided Design : Structural and Mechanical Systems. Springer-Verlang, Berlin, 1987. [12] Krischnaswani, P.; Whage, R. A.; Haug, E. J.: Design Sensitivity Analysis of Constrained Dynamic Systems by Direct Differentiation. Technical Report No. 83-5, Cemter for Computer-Aided Design, The University of Iowa, Iowa City, Iowa, 1983. [13] Maly, T.; Petzold, L. R.: Numerical methods and software for sensitivity analysis of differentialalgebraic systems. Journal of Applied Numerical Mathematics, Vol. 20, Issue 1-2, 1996. [14] Serban, R.; Freeman, J.S.: Direct differentiation methods for the design sensitivity of multibody dynamic systems. In ASME Design Engineering Technical Conferences and Computers in Engineering Conference, 96-DETC/DAC-1087, pp. 1–7, 1996. [15] Serban, R.; Feeman, J.S.: Identification and Identifiability of Unknown Parameters in Multibody Dynamic Systems. Multibody System Dynamics, Vol. 5, pp. 335–350, 2001. [16] Serna, M.A., Avilés, R., García de Jalón, J.: Dynamic analysis of plane mechanisms with lower pairs in basic coordinates. Mech. Mach. Theory, Vol. 17, pp. 397–403, 1982. [17] Walther, A., Griewank, A.: ADOL-C: A Package for the Automatic Differentiation of Algorithms Written in C/C++. Version 2.0.0, 2008. [18] Wang, X.; Haug, E.J.; Pan, W.; Implicit Numerical Integration for Design Sensitivity Analysis of Rigid Multibody Systems. Mechanics Based Design of Structures and Machines, Vol. 33, pp. 1–30, 2005. [19] Wehage, R.A., Haug, E.J.: Generalized coordinate partitioning for dimension reduction in analysis of constrained dynamic systems. ASME J. Mech. Des., Vol. 104, pp. 247–255, 1982.