{'fraction_non_alphanumeric': 0.0390474147178717, 'fraction_numerical': 0.01036973467782307, 'mean_word_length': 4.6215075376884425, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 12, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Adjoint Sensitivity Analysis of Hybrid Multibody Dynamical Systems

Sensitivity analysis of multibody systems computes the derivatives of general cost functions that depend on the system solution with respect to parameters or initial conditions. This work develops adjoint sensitivity analysis for hybrid multibody dynamic systems. Hybrid systems are characterized by trajectories that are piecewise continuous in time, with finitely-many discontinuities being caused by events such as elastic/inelastic impacts or sudden changes in constraints. The corresponding direct and adjoint sensitivity variables are also discontinuous at the time of events. The framework discussed herein uses a jump sensitivity matrix to relate the jump conditions for the direct and adjoint sensitivities before and after the time event, and provides analytical jump equations for the adjoint variables. The theoretical framework for sensitivities for hybrid systems is validated on a five-bar mechanism with non-smooth contacts.


Introduction
Hybrid dynamical systems are characterized piecewise-in-time smooth trajectories, with discontinuities appearing at a finite number of time moments as a consequence of external events. The discontinuities are characterized by a jump in the generalized velocity variables, e.g., caused by an impact or/and an abrupt change on the right hand side of the equation of motion.
Sensitivity analysis aims to quantify the effect of small changes in the system parameters (inputs) on a general cost function (outputs) [1]. Sensitivity analysis is essential in solving computational engineering problems such as design and control optimization, implicit time integration methods, and deep learning. Finite difference methods that approximate the sensitivities by the difference between perturbed and nominal trajectories are often inaccurate [2].
Two complementary approaches to sensitivity analysis are widely used, the direct and adjoint methods. While they provide the same derivatives, their approach and computational costs are different. Direct sensitivity propagates small perturbations forward through the system dynamics, while the adjoint method performs an inverse modeling that attempts to identify the origin of uncertainty in the model output [3].
The sensitivity analysis with respect to system parameters and initial conditions for hybrid systems has been studied in the literature [4][5][6][7][8][9][10][11][12][13]. The jump conditions of the direct sensitivities for hybrid ODE systems were first presented by Becker [13] in 1966 and a year latter by Rozenvasser [8]. Thirty years later, Galán el al. [5] presented sufficient conditions for the existence and uniqueness of these jump equations. Jump conditions involve the sensitivity of the time of event, and the jumps in the sensitivities of the state variables at the time of event. Within the same period, Hiskens applied this theory to power switching systems [11]. The jump conditions of the adjoint sensitivities for hybrid ODE systems with discontinuities in the right-hand side and with switching manifold parameters were presented by Stewart [14] and Taringoo [15], respectively. Recently, Zhang et al. [16] derived the jump conditions for adjoints of differential-algebraic systems and applied them to large-scale power systems with switching dynamics.
This paper provides a unified mathematical framework for the direct and adjoint sensitivity analysis for multibody dynamic systems and general cost functions. The framework includes both unconstrained and constrained mechanical systems. The direct sensitivity analysis was developed in [17], where a new graphical proof of the jump conditions for direct sensitivity variables was given. Jump conditions for constrained mechanical systems with change of mechanism and dealing with impulsive forces at the time of event was also presented.
This paper extends the mathematical framework to perform adjoint sensitivity analysis for mechanical systems with non-smooth trajectories. The approach taken here is as follows. An event detection mechanism (e.g., embedded in the numerical integration method) finds the time of the next event. At this time moment the trajectories of the generalized position variables are continuous but not differentiable, while the trajectories of the generalized velocities are discontinuous due to either impulsive external forces or to abrupt changes of the right-hand side. The jump conditions that map the direct sensitives from right before the event to right after the event can be formalized with the help of a jump sensitivity matrix. The jumps in the adjoint sensitivity variables are obtained via the transpose of this jump sensitivity matrix.
The paper is organized as follows. A review of the direct and adjoint sensitivity analyses for smooth dynamical systems with extended general cost functions is given in Section 2. The framework for direct and adjoint sensitivity analyses for hybrid dynamical systems is discussed in Section 3. The methodology is applied to compute sensitivities of a five-bar mechanism with non-smooth contacts in Section 4. Conclusions are drawn in Section 5.

Sensitivity analysis for unconstrained mechanical systems and extended cost functions
This section provides a summary of a previously method developed by the authors to implement direct sensitivity analysis for dynamical systems governed by smooth second order systems of ordinary differential equations (2nd order ODEs). More details of this method can be found here (cite the paper). This method is extended to multiple cost functions that contain argument function.

Smooth ODE system dynamics and extended cost functions
We consider an unconstrained mechanical system governed by the second order ordinary differential equation M (t, q, ρ) ·q = F (t, q,q, ρ) , t 0 ≤ t ≤ t F , q(t 0 ) = q 0 (ρ),q(t 0 ) =q 0 (ρ), ⇔q = M −1 (t, q, ρ) · F (t, q,q, ρ) =: f eom (t, q,q, ρ) , (1) where q ∈ R n are the generalized positions, v :=q ∈ R n the generalized velocities, and ρ ∈ R p the time independent parameters of the system. The state trajectories depend implicitly on time and on the parameters, q = q(t, ρ) and v =q(t, ρ). We consider a general system output of the form: The functiong : R 1+2n+p → R n c is a vector of 'trajectory cost functions', andw : R 1+2n+p → R n c is a vector of 'terminal cost functions', and the system output ψ ∈ R n c is a vector of n c 'outputs', i.e., scalar cost functions.
Both the trajectory and terminal cost functions can include accelerations viav. Accelerations are not independent variables and can be resolved in terms of positions and velocitiesv := f eom (t, q,q, ρ) ∈ R n . The cost functions can also include argumentsũ t, q, v, ρ = u ( t, q, v,v, ρ) that depend on the solution and on the acceleration. Our notation encompasses these cases by defining: All functions are considered to be smooth.
Definition 1 (Sensitivity analysis problem). The sensitivity analysis problem is to compute the derivatives of the model outputs (2) with respect to model parameters: Definition 2 (The canonical ODE system). To simplify the representation of the system we define the vector of 'quadrature' variables z ∈ R n c as follows: which leads the vector of cost function (2) at final time to become: Next, we add dummy evolution equations ρ = 0 for the time independent parameters.
Finally, we append the parameters and the quadrature variables to system state to obtain the following extended 'canonical' state vector: together with the 'canonical ODE system' that describes its evolution:

Direct sensitivity analysis for smooth ODE systems and extended cost function
Define the 'position sensitivity' matrix Q(t, ρ), the 'velocity sensitivity' matrix V(t, ρ), the 'quadrature sensitivity' matrix Z(t, ρ), and an identity matrix Γ as the formal sensitivity of the parameters, as: The direct sensitivity for ODE systems, referred to as the Tangent Linear Model (TLM), computes the sensitivity matrix X = Q T , V T , Γ, Z T T ∈ R (2n+p+n c )×p . obtained by differentiating the canonical ODE system (7) with respect to the parameters: The direct sensitivity for ODE systems needs to be solved forward in time. The expressions f eom q , f eom v , and f eom ρ i denote the partial derivatives of f eom with respect to the subscripted variables. The detailed calculation of these expressions and the remaining partial derivatives is explained in Appendix A. Once the sensitivities (9) have been calculated, the sensitivities of the cost functions (4) with respect to parameters are computed as follows: We note that the TLM system (9) can be written in matrix form as follows:

Adjoint sensitivity analysis for smooth ODE systems and extended cost function
In this section we provide the system of equations that governs the adjoint sensitivity analysis for smooth ODE systems.
Definition 3 (Adjoint sensitivity analysis). Apply the chain rule differentiation to the total sensitivity of the cost function (4): where Note that, from (5)-(6) which leads to the relation dψ/dz(t, ρ) = I n c ×n c for any time t.
From (12) we have that for any time t: Evaluating (14) at t = t F leads to the direct sensitivity approach: By comparing this equation with (10) one obtains the values of the adjoint variables at the final time t F : The equation (14) evaluated at t = t F leads to the direct sensitivity approach: Evaluating (14) at t = t 0 leads to the adjoint sensitivity approach: Note that the adjoint variables are initialized at t = t F . However, their values at t = t 0 are the ones needed for computing the desired sensitivities.
Definition 4 (The canonical adjoint sensitivity for ODE systems). The evolution of adjoint variables for ODE systems is governed by the following continuous adjoint model: The adjoint sensitivities (19) are solved backward in time.

Hybrid ODE system dynamics
In this study, we consider hybrid ODE systems characterized by piecewise-in-time smooth dynamics described by (1), and that exhibit discontinuous dynamic behavior (jump or non-smoothness) in the generalized velocity state vector at a finite number of time moments (no zeno phenomena [18]). Each such moment is corresponds to an event triggered by the event equation: where t eve is the 'time of event' and r : R n → R is a smooth 'event function'. Note that grazing phenomena are not considered the in this study. The following quantities are used to characterize an event: • The value of a variable right before the event is denoted by x| − t eve := lim ε>0, ε→0 x(t eve − ε), and its value right after the event by x| + t eve := lim ε>0, ε→0 x(t eve + ε). The limits exist since the evolution of the system is smooth in time both before and after the event.
• The generalized position state variables remain the same after the event as before it, q| + t eve = q| − t eve = q| t eve . This is a consequence of the event changing the energy of the system by a finite amount.
• Also due to the finite energy change during the event, the quadrature variable is continuous in time, • An event that applies a finite energy impulse force to the system can abruptly change the generalized velocity state vectorq, from its value v| − t eve right before the event to a new value v| + t eve right after the event. The change in velocity is characterized by the 'jump function': • An event where the system undergoes a sudden change of the equation of motions (1) at t eve is characterized by the equations: Remark 1 (Multiple events). In many cases, the change can be triggered by one of multiple events. Each individual event is described by the event function r : R n → R, = 1, . . . , e. The detection of the next event (20), which can be one of the possible e options, is described by Π e i=1 r i q| t eve = 0, and if event takes place, then r = 0 and the corresponding jump in velocity (21) is v| + t eve = h q| t eve , v| − t eve , or the corresponding change in the equations of motion (22) isq| + t eve =: f eom + | t eve .

Direct sensitivity analysis for hybrid ODE systems
Let Q| + t eve and Q| − t eve ∈ R n×p be the sensitivities of the generalized position state matrix after and before the event, respectively. Let V| + t eve , V| − t eve ∈ R n×p be the sensitivities of the generalized velocity state matrix after and before the event, respectively. Let Z| + t eve and Z| − t eve , with Z ∈ R p , be the sensitivities of the quadrature variable z(t) after and before the event, respectively. It is shown in [17] that, at the time of the event, we have: • The sensitivity of the time of event with respect to the system parameters is: where dr/dq ∈ R 1×n is the Jacobian of the event function.
• The jump equation of the sensitivities of the generalized position state vector is: • The jump equation of the sensitivities of the generalized velocity state vector is: where the Jacobians of the jump function are: • The sensitivity of the cost function changes during the event is : is the running cost function evaluated right after and right before the event, respectively.
Definition 5 (The generalized jump sensitivity matrix). The direct sensitivity jump equations (23) can be written compactly in matrix form as X| + t eve = S · X| − t eve , where S is the generalized sensitivity jump matrix: From (23a) we have that: The Jacobians Q| + where S T eve is the transpose of the generalized sensitivity jump matrix (24a).
Proof. We start the proof from the following statement provided in [14] that mentions that the dot product of the sensitivity state matrix with the adjoint sensitive state matrix is constant at any time, λ| + It is shown in 3 that the terminal cost functionw cannot directly depend on the accelerationv or on the Lagrange coefficients µ, and therefore the derivatives arewv = 0 andw µ = 0. In a different notation, such result is also shown in [19]. Using equation (10) we see that the final condition for the adjoint of the algebraic variables µ is zero,

Direct and adjoint sensitivity analysis for smooth systems in the penalty ODE formulation
Define the extended mass matrix M : R × R n × R n × R p → R n×n and the extended right hand side function where α ∈ R m×m is the penalty factor of the ODE penalty formulation, ξ ∈ R and ω ∈ R are the natural frequency and damping ratio coefficients of the formulation, respectively. The functions Φ,Φ,Φ : R × R n × R n × R p → R m are the position, velocity and acceleration kinematic constraints, respectively. The penalty formulation of a constrained rigid multibody system is written as a first order ODE: The Lagrange multipliers associated to the constraint forces are estimated as µ * = α Φ + 2 ξ ωΦ + ω 2 Φ . The sensitivities of the state variables of the system with respect to parameters evolve according to the tangent linear model derived in [1,[20][21][22][23]. Since the penalty formulation (31) evolves as an ODE, we can compute the direct sensitivities F q , F v , and F ρ are given in [17]. Similarly, one can compute the adjoint sensitivities of the penalty formulation (31) using (19).

Direct and adjoint sensitivity analysis for smooth systems in the index-1 differential-algebraic formulation
Definition 6 (Constrained multibody dynamics: the index-1 DAE formulation). The index-1 formulation of the equations of motion is obtained by replacing the position constraint (27a) with the acceleration constraint (27c): The algebraic equation has the form f DAE-µ − µ = 0.
Definition 7 (Tangent linear index-1 DAE). Sensitivities of solutions (8) and multipliers: of the system (32) with respect to parameters evolve according to the tangent linear model derived in [1,[20][21][22][23]: It is shown in Appendix A that equation (34) can be written in matrix form as follows: with initial conditions given by Eq. (9). Using Appendix A, the derivatives of the DAE function are: Noting from Remark 3 that the algebraic equation in (36) reads: the index-1 adjoint DAE (36) can be reduced to the following adjoint ODE:

Direct sensitivity analysis for hybrid constrained dynamical systems
We now discuss constrained dynamical systems when the dynamics is piecewise smooth in time. Performing a sensitivity analysis for a constrained rigid hybrid multibody dynamic system requires finding the jump conditions at the time of event. These jump equations are explained in our previous work [17]. We summarize below the jump equations at the time of event: • The generalized position state variables remain the same , i.e., q| + t eve = q| − t eve = q| t eve and need to satisfy both constraint functions Φ − | − t eve := Φ − t eve , q| t eve , ρ = 0, and Φ + | + t eve := Φ + t eve , q| t eve , ρ = 0.
• The velocity state variables jump from their values right before the event to right after the event according to the jump equation: The jump function (38) is assumed to be smooth and defined in terms of the velocity degrees of freedom (the independent components).
• The jumps in velocity cannot be arbitrary for the dependent components. They are dependent of the degree of freedom and are obtained from solving the velocity constraints leading to: Where R ± corresponds to the null space of the constraints if the constraints are scleronomic (non explicitly time dependent).
There are two types of velocity jumps that our formalism covers (38)-(39): • The case where the event consists of an elastic contact/collision/impact on the DOF components of the velocity state. The impulsive (external) contact forces act to change the DOF components without changing the set of constraint equations, Φ + ≡ Φ − .
• The case where the event consists solely of an inelastic collisions and a change of constraints Φ + Φ − , without any external force modifying the independent velocities. The impulsive (internal) constraints forces at the time of event are solved by using a popular approach in robotics [24]: or, equivalently, The second equation (40b) imposes the velocity constraint on both independent and dependent coordinates, which is covered by our formalism as: is a permutation matrix that partitions the state variables into dependent and independent variables.
Finally, the jump conditions at the time of event in the sensitivity state matrix are: • The independent components of the sensitivity of the generalized positions right after the event: which are equivalent to: • The dependent components of the sensitivity of the generalized positions right after the event: (42c) • The independent coordinates of the velocity sensitivities right after the event, • The dependent components of the velocity sensitivities right after the event, Definition 9 (The generalized sensitivity jump matrix for elastic impact). The jump equations (38)-(44) for constrained systems can be written compactly in matrix form as a jump of the state sensitivity matrix X at the time of the event, X| + t eve = S eve · X| − t eve , where S eve represents the generalized jump sensitivity matrix: The Jacobians of the jump equations with respect to the sensitivity state before the time of event are: with It follows that: and Rewriting (44) as: or, equivalently, as: we find the following expressions for the Jacobians: The expressions for D and K in (45) are: (53) Definition 10 (The generalized sensitivity jump matrix for inelastic impact with a sudden change of constraints).
Consider the event consisting of an inelastic collision and a sudden change of constraints (40). The jump in the velocity sensitivity for constrained systems due to impulsive forces, presented in [17], is determined as follows: which simplifies to: Thus, the jump the velocity state variables at the time of event is and the jump in the sensitivity of the Lagrange multipliers from Λ| − t eve → Λ| + t eve is: The corresponding sensitivity jump matrix (45) is:

Adjoint sensitivity analysis for hybrid constrained dynamical systems
Definition 11 (Jump in adjoint sensitivity for constrained systems with elastic impact). The transpose of the direct sensitivity jump matrix S T eve (25) associated with an elastic impact (45) is: From the adjoint sensitivity equation (25) the jumps in adjoint variables for ODE systems with constraints undergoing an elastic impact are:

Definition 12 (Jump in adjoint sensitivity for constrained systems with inelastic impact and a sudden change of constraints).
The transpose of the direct sensitivity jump matrix S T eve (25) associated with (58) is: dependent coordinates, are defined according to the fixed lengths between each set of points, as follows:  The trajectories of the vertical position and velocity of point 2 of the five-bar mechanism are presented in Fig. 2a and Fig. 2b, respectively. As expected, the point 2's vertical position bounces at -2.35m, and its vertical velocity jumps at each time of event with v 2 | + t eve = v 2 | − t eve .  The trajectory of the quadrature variable z(t) = t t 0ẏ 2 (τ) dτ of the five-bar mechanism and its sensitivity are shown in Fig. 3a and Fig. 3b, respectively. The same analysis is presented in Fig. 4a and Fig. 4b for the quadrature variable 2 dτ. The direct sensitivity is represented by the continuous line, while the central finite difference sensitivity is represented by the dashed line. Both solutions were solved forward in time. The adjoint sensitivity is presented as well, and was solved backwards in time. As presented in our previous paper, the direct differentiation method to compute the sensitivity of the cost function with discontinuities in the velocity state variables of the mechanism is validated. The validation comes to the fact that the trajectories of the sensitivity of the quadrature variable Z( t, q, v, ρ ) exactly matches the trajectory of numerical sensitivity computed with a finite difference method. One main conclusion of our previous paper was to state that our proposed direct sensitivity method in computing the sensitivity of the cost function with discontinuities in the velocity state variables was more robust than the numerical method. Indeed, the direct method accurately determines the jump in the sensitivities and their trajectories. This after each event, without any delta-like jumps in magnitude 1/ε that occurs in the numerical method at each time of event. This validated direct sensitivity method is now compared to the proposed adjoint method in computing the sensitivity of the cost function with discontinuities in the velocity state variables. The results presented in Fig. 3a and Fig. 4a show that the adjoint and direct method exactly converge to the same sensitivity cost number with a difference of less than 0.01 %. This convergence in both methods validates the adjoint sensitivity method in computing the sensitivity of the cost function with discontinuities in the trajectories.    The same analysis is provided with the quadrature variable z(t) = t t 0ÿ 2 (τ) 2 +ẏ 2 (τ) 2 dτ in Fig. 5a and Fig. 5b. The adjoint and direct method converge to the same sensitivity cost number with a difference of less than 0.01 %.

Conclusions
Gradient based algorithms are widely used in computational engineering problems such as design and control optimization, implicit time integration methods, and deep learning. Sensitivity analysis plays a key role in this type of algorithms as it provides the necessary derivative information. In the context of dynamical systems governed by ordinary or differential algebraic equations, sensitivity analysis computes the derivatives of general cost functions that depend on the system solution with respect to parameters or initial conditions. Direct and adjoint sensitivity analyses for continuous multibody dynamic systems have been discussed in the literature [1,[20][21][22][23]. Our earlier work has extended the direct sensitivity analysis to hybrid multibody dynamic systems systems that are subject to events such as impacts or sudden changes in constraints [17]. This paper extends the mathematical framework to compute adjoint sensitivities for hybrid multibody dynamic systems modeled by ordinary differential equations and by index-1 differential algebraic equations. A very general formulation of the cost functions is used. For the hybrid systems considered herein discontinuities in the forward trajectories appear at time moments triggered by an event. Jump conditions for adjoint sensitivity variables are provided for mechanical systems with and without constraints. These jump conditions handle the change in the sensitivities caused by the non-smoothness of the forward trajectories at a finite number of events.
We validate the mathematical framework for adjoint sensitivities for hybrid dynamical systems on the study of a five-bar mechanism with non-smooth contacts. The direct and adjoint sensitivities computed by the proposed mathematical framework are validated against numerical sensitivities calculated by real finite differences. The results of this study show that all the alternative analyses provide the same sensitivities of the general cost function with respect to model parameters, within an error of 0.01%. Future work will extend the mathematical framework to calculate adjoint sensitivities of hybrid mechanical systems with respect to actuation functions. These sensitivities will allow to solve optimal control problems for hybrid systems.
A. Calculation of partial derivatives used in sensitivity analyses Remark 5. The expressions f eom q , f eom v , and f eom ρ i denote the partial derivatives of f eom with respect to the subscripted variables. The partial derivatives ∂ f eom /∂ζ are obtained by differentiating f eom with respect to ζ ∈ {q, v, ρ}: Remark 6. The expressionsg q ,g v , andg ρ i denote the partial derivatives ofg with respect to the subscripted variables.
The partial derivatives ∂g/∂ζ are obtained by differentiating (1) with respect to ζ ∈ {q, v, ρ}: which leads to: ..,p = g q + gv f eom q + g u u q + g u uv f eom Remark 7. Similarly, the expressionsw q ,w v , andw ρ i denote the partial derivatives ofw with respect to the subscripted variables. The partial derivatives ∂w/∂ζ are obtained by differentiating w with respect to ζ ∈ {q, v, ρ}: = w ζ + wv f eom ζ + w u u ζ + w u wv f eom ζ .

B. Adjoint of the algebraic Lagrangian coefficient
Methods to compute the adjoint of an index-1 DAE available in the literature [19,25,26] use the following approach. Define the Lagrangian using the multipliers µ Q , µ V , µ Γ that correspond to the constraints posed by the index-1 DAE equations (32): We rearrange equation (B.1) as follows: The adjoint variables λ Q , λ V , λ Γ defined in this paper, and the adjoint variables µ Q , µ V , µ Γ used in the literature (B.1), are related by the following matrix multiplication: The adjoint DAE equations and boundary conditions in the "µ formulation" [19,25,26] can be derived from the equations and boundary conditions in the "λ formulation" discussed in this paper, and vice-versa. Double subscripts indicates a three-dimensional Jacobian with respect to a quantity ζ and φ, unless stated otherwise ζ Subscript indicates partial derivative with respect to a quantity ζ, unless stated otherwise

Sensitivity Analysis
Q ∈ R n×p The sensitivity matrix of the state vector q with respect to the vector of system parameters ρ V ∈ R n×p The sensitivity matrix of the state vectorq with respect to the vector of system parameters ρ X ∈ R (2n+p+n c )×p The sensitivity matrix of the x state vector with respect to the vector of system parameters ρ dt eve /dρ ∈ R 1×p The sensitivity of the time of event t eve with respect to the vector of system parameters ρ λ∈ R (2n+p+n c )×n c The adjoint sensitivity matrix of X λ Q ∈ R n×n c The adjoint sensitivity matrix of Q λ V ∈ R n×n c The adjoint sensitivity matrix of V ψ ∈ R n c The vector of cost functions g ∈ R n c The vector of trajectory cost functions w ∈ R n c The vector of terminal cost functions Z ∈ R n c ×p The sensitivity matrix of the vector of quadrature variables z