Enhancement of the Adjoint Method by Error Control of Accelerations for Parameter Identification in Multibody Dynamics

The present paper shows the embedding of the adjoint method in multibody dynamics and its broad applicability for examples for both, parameter identification and optimal control. Especially, in case of parameter identifications in engineering multibody applications, a theoretical enhancement of the proposed adjoint method by an error control of accelerations is inevitable in order to meet the circumstances of experimental studies using acceleration sensors in general.


Introduction
In the explanations of Giles and Pierce [1], the adjoint method is seen as a special case of linear duality, in which the dual problem has to solve only a single linear system. Obviously, huge benefits can be achieved by solving the dual formulation. The adjoint method utilizes this powerful aspect of duality to dramatically improve the efficiency of the computation.
The adjoint method is probably the most efficient way to solve a wide range of optimization problems in engineering sciences. In the sensitivity analysis, the adjoint method has been utilized by various authors, e.g., [12,8,7], and in dynamic programming methods the adjoint approach has been used extensively for the computation of gradients as well [10].
The adjoint method is applied in, e.g., aerodynamic shape optimization by Jameson [13], in which the gradient of an objective function is determined indirectly by solving an adjoint equation which has coefficients determined by the solution of the multibody dynamics equations. This directly corresponds to the gradient technique for trajectory optimization pioneered by Bryson and Ho [9]. The fast computation of the gradients makes optimization computationally feasible even for complex three-dimensional multibody systems.
There are two strategies for this purpose, the equations of motion of the multibody system and adjoint equations may either be separately discretized from their representations as differential-algebraic equations, or, alternatively, the equations of motion of the multibody system may be discretized first, while the discrete adjoint equations are derived directly from the discrete multibody system equations [9].
Previous work on the adjoint method in multibody dynamics can be found in [14], where the solution of inverse dynamics and trajectory optimization problems for multibody systems is reflected by an indirect approach combining optimal control theory with control and adjoint equations with transversality conditions. The design of the indirect method for solving optimal control problems for multibody systems presented in [18] seems to be similar to the idea to solve adjoint equations.
The work by Schaffer [15] proposes a numerical algorithm, the piecewise adjoint method, which formulates the coordinate partitioning underlying ordinary differential equations as a boundary value problem, which is solved by multiple shooting methods.
The group around Petzold, Cao, Li and Serban [17,16] describe forward and adjoint methods for sensitivity analysis for differential-algebraic equations and partial differential equations and emphasizes that the results of sensitivity analysis have a wide range of applications, including model development, optimization, parameter estimation, model simplification, data assimilation, optimal control, uncertainty analysis and experimental design [17].
Eberhard [7] uses the adjoint method for sensitivity analysis in multibody systems interpreted as a continuous, hybrid form of automatic differentiation.
In case of an orientation parametrization of a body without using angles, as e.g. described in the absolute nodal coordinate formulation (ANCF) [22], a gradient-based optimization approach using adjoint equations has been proposed by Held and Seifried [19]. The framework of the ANCF combined with the adjoint method is as well used in the sensitivity analysis by Pi et al. [21] within a first order approach, while in [20] a second order adjoint sensitivity analysis of multi-48 Enhancement of the Adjoint Method by Error Control of Accelerations for Parameter Identification in Multibody Dynamics body systems has been presented within the classical floating frame of reference formulation for multibody systems. For a comparison of the computational costs between the adjoint method and the direct differentiation method see Table 4 in [20]. Despite of the dramatically decrease of computational costs and the resulting potential of the adjoint method in multibody dynamics, the adjoint method is rarely applied. The structure of the equations of motion in a multibody system is usually quite complicated, and the effort to obtain the set of adjoint equations seems tremendously high.
The main goal of the present paper is to show how the adjoint method can be embedded efficiently to a multibody system described by a system of differential-algebraic equations of index 3 for optimal control problems or parameter identification applications. The body description is, in contrast to the approach in [19], in which the flexible multibody system is displayed in minimal coordinates, defined by the position of the center of mass and the four redundant Euler parameters for the orientation parametrization. This choice includes an internal constraint for each body which has to be considered also in the constraint Jacobian of the whole system and of course also influences the optimal control equations. The present paper shows the potential of the adjoint method for solving classical optimization problems in multibody dynamics and enhances the governed equations by an error control of accelerations for parameter identification in order to meet the standards of experimental studies using acceleration sensors in general.

Adjoint gradient computation
The basic idea of the adjoint method, see e.g. [17,16,7], is summarized in brief here. Suppose we have the semi-explicit differential algebraic systeṁ where x(t) ∈ R Nx and λ(t) ∈ R N λ are the vectors of state variables and algebraic variables of a dynamical system such as a multibody system in a redundant formulation, e.g., using Euler parameters for the parametrization of the orientation. Herein, u ∈ R Nu may either describe a vector of (constant) parameters or a vector of control signals. Our goal is to find u such that a functional of the form is minimized. The focus is to determine a set of controls or parameters u in a way that the measured signalss i (t), i = 1, . . . , N s are best approximated by the system outputs s i (x).
In this case one could define the root mean square (RMS) error as our functional to be minimized Various methods are available to compute the argument for which a function or a functional attains a minimum, e.g., the method of the steepest descent, the conjugate gradient method, the Gauss-Newton method or quasi Newton methods. In any cases the gradient of J(u) has to be determined and for this purpose several strategies can be pursued again.
On the one hand, if u is a vector of N u parameters, the sensitivity equations for x u = ∂x/∂u are usually considered, see e.g. [11] or [23]. When choosing this sensitivity approach, the computational effort is equal to solving N u linear sets of equations with the same dimension as Eq. (1). On the other hand, if u(t) is a vector of control signals, the problem could be transformed into a finite dimensional one by discretizing the control signals. The adjoint method is an efficient alternative to compute the gradient of J(u) in both cases.
First, it has to recognized that J(u) does not change if we add Eq. (1) to the integrand with arbitrary functions y(t) and µ(t). The Hamiltonian Let us now consider a forward solution x(t) and λ(t) of the system equations in Eq. (1) for a set of parameters or control variables u. A variation of u about δu will result in variations of x(t) and λ(t) about δx(t) and δλ(t) and, moreover, in a variation of the functional J about δJ. Up to a first order approximation δJ is given by where H u , H x and H λ denote the row vectors with the partial derivatives of H with respect to the components of u, x and λ, e. g., Integrating by parts of the last term of Eq. (7) yields since δx(0) = 0 as the initial state is assumed to be prescribed. It is the key idea of the adjoint method that the computation of δx and δλ can be circumvented if the adjoint variables y(t) and µ(t) are chosen such thaṫ y T = −H x , H λ = 0 and y(T ) = 0 (9) or including an other known fixed value for the boundary condition y(T ). This set of semi-explicit differential algebraic equations is called Then, with x(t), λ(t), y(t) and µ(t) from Eqs. (1) and (9), the variation of J according to Eq. (8) is readily given by In case u is a control signal, the largest possible increase of δJ is achieved, if δu(t) is chosen in the direction of H T u . Hence, H T u may be considered as the gradient of the functional J(u).
It has to be mentioned here that in case u is a control signal, the integral in Eq. (10) is computed as a sum including δu i computed for the i-th time step and ∆t is chosen as the time discretization of the control signal itself. In case u is a vector of parameters, the gradient of the function J(u) is given by ∇J = ∫ T 0 H T u dt. It has to be emphasized here, that only two systems of DAEs must be integrated for computing the direction of the gradient. The main difficulty with multibody dynamics results from the complexity of the original system (1). Hence, many authors focused on two-dimensional examples or rather general aspects, see e.g. [11], [6]. A former paper uses this key idea of the adjoint method and shows how these equations can be implemented in a MBS-code under the assumption that the rigid body rotations are described by Euler parameters, since a redundant formulation gives relatively simple adjoint equations, see [5]. The present paper extends the theory presented in [5] by enhancing the method by an error of accelerations in the cost functional.

The adjoint method enhanced by an error of accelerations for parameter identification
This section discusses the embedding of the adjoint gradient computation to multibody dynamics.
In case the measures included in the system are restricted to position and velocity level, the embedding of the adjoint method to multibody dynamics has been presented already in Nachbagauer et al. [5]. In case the error of accelerations has to be accounted for in the optimization functional, the adjoint method in multibody dynamics presented in [5] must be enhanced. Herein, the enhancement of the adjoint method by error control of accelerations for parameter identification is presented.
Firstly, the adjoint equations for the equations of motion of a multibody system are derived and, secondly, the structure of the boundary conditions is discussed.

Adjoint equations
A mechanical system consisting of rigid and flexible bodies, forces and constraints acting between these bodies can be described by equations of motion in the following form: Here, q denotes a vector of redundant generalized coordinates. They are subject to the holonomic constraints C(q) = 0, which enter the equations of motion via the Jacobian C q and the vector of Lagrange multipliers λ representing the constraint forces in the system. Moreover, the system may incorporate a control u(t) or a vector of parameters u. For simplicity, we suppose that u appears only in the vector f , which is the fact, if u is an actuating force or a stiffness or damping parameter. Using the additional variables v =q, the equations of motion can be reformulated as a first order system of equations aṡ The state vector x defined in the previous section now consists of q and v. The goal of the optimization is to find the parameter (or control) u which minimizes a cost functional of the form (15) in which S(q, v)| T introduces an additional end point term, usually denoted as scrap function. For example, S could describe an end point RMS-error of a sensor variable similar to Eq. (3). The cost functional J is augmented by the system equations as follows The variables p(t), w(t) and µ(t) may be chosen arbitrarily here. The variation of the cost functional J is derived by Integration by parts for the terms including δq and δv is computed by with δq(0) = δv(0) = 0 since initial conditions for the adjoint variables are prescribed as q(0) = q 0 and v(0) = v 0 . Using Eqs. (19) and extracting the terms multiplied with δq, δv, δλ and δu, the variation of the functional J given by Eq. (17) can be rewritten as To eliminate the terms involving δq, δv the adjoint variables p, w and µ may now be defined by equating the respective expressions in square brackets to zero. After transposing these expressions one obtains the following system of adjoint equations: where the abbreviations and simplifications due to the symmetry of the mass matrix M = M T have been used. Note that in case a root mean square error has to be minimized, as similar to Eq. (3),  S(q, v,v) as well. Concluding this derivation, in case Eqs. (21)-(25) are satisfied the variation of the cost functional in Eq. (20) reduces to which directly relates the independent variation δu to the variation of the objective function.

Consistent boundary conditions for the adjoint system
It has to be noticed that the boundary condition (25) for the adjoint variable w is generally incompatible with the constraint equation (23) at time t = T . Only in case S v = 0, i. e. when the scrap function does not depend on v, all equations are satisfied by setting p(T ) = −S T q and w(T ) = 0. However, we run into trouble if S v ̸ = 0. To circumvent this issue an approach similar to the idea of Gear-Gupta-Leimkuhler [3] is applied by considering the constraint equations C(q) = 0 of the original system Eq. (11) at velocity level and its variation with respect to v and q results in This relation is considered at time t = T and is multiplied by an arbitrary vector of numbers ξ leading to Since the expression on the left side is always zero one could add it to Eq. (20) without modifying the variation of our objective functional. Hence, two additional boundary terms in Eq. (20) arise transforming Eqs. (24) and (25) into Using the still undetermined variable ξ, the boundary conditions p(T ), w(T ) can be computed in such a way that also the constraint equation (23) is satisfied at time t = T . For that purpose, the following system of equations has to be solved for p(T ), w(T ) and ξ at time t = T : This means that one may first solve at time t = T : yielding ξ and w(T ), and compute then p(T ) at time t = T from Once w(T ) and p(T ) has been determined in the correct way, the differential-algebraic set of adjoint equations For the adjoint equations of a multibody system the matrices M, C q , A and B from Eqs. (26) and (27) are required along a forward simulation of the equations of motion as well. The matrices M, C q also appear in the equations of motion and are therefore available, while the determination of the Jacobian matrices A and B requires additional computational effort.
In [5], more details on the different computation strategies for the Jacobians have been presented. Beside exact computation of the Jacobian matrices by implementing explicit formulas in the MBS-code, one could also compute the derivatives numerically by substituting the derivatives by finite difference quotients or investigate the technique of automatic differentiation. Despite the explicit formulas seems expensive and susceptible to programming errors, it becomes more attractive in case a highly redundant formulation of the equations of motion is used. The use of Euler parameters for the parametrization of rotational motions simplifies the system matrices such that programming explicit formulas for A and B becomes the most efficient strategy. It should also be noted, that the Jacobian matrices (Mv) q , f q , f v and (C T q λ) q may be required already for the simulation of the multibody system if an implicit integration scheme such as, e. g., the HHT-algorithm [2,4] is applied. Hence, an efficient computation of these matrices is crucial also for solving the equations of motion in forward direction.

Conclusions
The proposed method shows the embedding of the adjoint method in multibody system dynamics and its potential for parameter identification and optimal control on account of its linear structure. The enhanced adjoint method includes the error of accelerations in the optimization functional and therefore delivers an elegant and efficient way to incorporate inverse dynamics for parameter identification in flexible multibody system applications arising from modern engineering problems including acceleration sensors.