Adjoint Sensitivity Analysis of Hybrid Multibody Dynamical Systems
CComputational Science Laboratory TechnicalReport CSL-TR-18 - February 21, 2018
Sebastien Corner, Corina Sandu, Adrian Sandu “ Adjoint Sensitivity Analysis of HybridMultibody Dynamical Systems ” Computational Science LaboratoryComputer Science DepartmentVirginia Polytechnic Institute and State UniversityBlacksburg, VA 24060Phone: (540)-231-2193Fax: (540)-231-6075Email: [email protected]
Web: http://csl.cs.vt.edu
Compute the Future a r X i v : . [ m a t h . O C ] F e b djoint Sensitivity Analysis of Hybrid Multibody Dynamical Systems Sebastien Corner a,b , Adrian Sandu b , Corina Sandu a a Department of Mechanical Engineering, Virginia Tech, Blacksburg, VA 24061 b Computational Science Laboratory, Department of Computer Science, Virginia Tech, Blacksburg, VA 24061
Abstract
Sensitivity analysis of multibody systems computes the derivatives of general cost functions that depend on the sys-tem solution with respect to parameters or initial conditions. This work develops adjoint sensitivity analysis forhybrid multibody dynamic systems. Hybrid systems are characterized by trajectories that are piecewise continuousin time, with finitely-many discontinuities being caused by events such as elastic / inelastic impacts or sudden changesin 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 ad-joint 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-smoothcontacts. Keywords:
Direct and adjoint sensitivity analysis, hybrid dynamics, jump conditions, constrained multibody systems
1. Introduction
Hybrid dynamical systems are characterized piecewise-in-time smooth trajectories, with discontinuities appearingat a finite number of time moments as a consequence of external events. The discontinuities are characterized by ajump in the generalized velocity variables, e.g., caused by an impact or / and an abrupt change on the right hand sideof the equation of motion.Sensitivity analysis aims to quantify the e ff ect of small changes in the system parameters (inputs) on a generalcost function (outputs) [1]. Sensitivity analysis is essential in solving computational engineering problems such asdesign and control optimization, implicit time integration methods, and deep learning. Finite di ff erence methods thatapproximate the sensitivities by the di ff erence between perturbed and nominal trajectories are often inaccurate [2].Two complementary approaches to sensitivity analysis are widely used, the direct and adjoint methods. While theyprovide the same derivatives, their approach and computational costs are di ff erent. Direct sensitivity propagates smallperturbations forward through the system dynamics, while the adjoint method performs an inverse modeling thatattempts 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 beenstudied in the literature [4–13]. The jump conditions of the direct sensitivities for hybrid ODE systems were first Email addresses: [email protected] (Sebastien Corner), [email protected] (Adrian Sandu), [email protected] (Corina Sandu)
Preprint submitted to Nonlinear Analysis: Hybrid Systems February 21, 2018 resented by Becker [13] in 1966 and a year latter by Rozenvasser [8]. Thirty years later, Gal´an el al. [5] presentedsu ffi cient conditions for the existence and uniqueness of these jump equations. Jump conditions involve the sensitivityof the time of event, and the jumps in the sensitivities of the state variables at the time of event. Within the sameperiod, Hiskens applied this theory to power switching systems [11]. The jump conditions of the adjoint sensitivitiesfor hybrid ODE systems with discontinuities in the right-hand side and with switching manifold parameters werepresented by Stewart [14] and Taringoo [15], respectively. Recently, Zhang et al. [16] derived the jump conditions foradjoints of di ff erential-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 multibodydynamic systems and general cost functions. The framework includes both unconstrained and constrained mechanicalsystems. The direct sensitivity analysis was developed in [17], where a new graphical proof of the jump conditions fordirect sensitivity variables was given. Jump conditions for constrained mechanical systems with change of mechanismand 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 systemswith non-smooth trajectories. The approach taken here is as follows. An event detection mechanism (e.g., embeddedin the numerical integration method) finds the time of the next event. At this time moment the trajectories of thegeneralized position variables are continuous but not di ff erentiable, while the trajectories of the generalized velocitiesare discontinuous due to either impulsive external forces or to abrupt changes of the right-hand side. The jumpconditions that map the direct sensitives from right before the event to right after the event can be formalized with thehelp of a jump sensitivity matrix. The jumps in the adjoint sensitivity variables are obtained via the transpose of thisjump sensitivity matrix.The paper is organized as follows. A review of the direct and adjoint sensitivity analyses for smooth dynamicalsystems with extended general cost functions is given in Section 2. The framework for direct and adjoint sensitivityanalyses for hybrid dynamical systems is discussed in Section 3. The methodology is applied to compute sensitivitiesof a five-bar mechanism with non-smooth contacts in Section 4. Conclusions are drawn in Section 5.
2. 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 sensitivityanalysis for dynamical systems governed by smooth second order systems of ordinary di ff erential equations (2ndorder ODEs). More details of this method can be found here (cite the paper). This method is extended to multiple costfunctions that contain argument function. 2 .1. Smooth ODE system dynamics and extended cost functions We consider an unconstrained mechanical system governed by the second order ordinary di ff erential equation(ODE): M ( t , q , ρ ) · ¨ q = F ( t , q , ˙ q , ρ ) , t ≤ t ≤ t F , q ( t ) = q ( ρ ) , ˙ q ( t ) = ˙ q ( ρ ) , ⇔ ¨ q = M − ( 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 independentparameters 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: ψ ( ρ ) = (cid:90) t F t ˜ g j (cid:0) τ, q , v , ρ (cid:1) d τ + ˜ w (cid:0) t F , q t F , v t F , ρ (cid:1) , q t F : = q ( t F , ρ ) , v t F : = v ( t F , ρ ) . (2)The function ˜ g : R + n + p → R n c is a vector of ‘trajectory cost functions’, and ˜ w : R + n + p → R n c is a vectorof ‘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 via ˙ v . Accelerations are not independentvariables and can be resolved in terms of positions and velocities ˙ v : = f eom ( t , q , ˙ q , ρ ) ∈ R n . The cost functions can alsoinclude arguments ˜ u (cid:0) t , q , v , ρ (cid:1) = u ( t , q , v , ˙ v , ρ ) that depend on the solution and on the acceleration. Our notationencompasses these cases by defining:˜ g (cid:0) t , q , v , ρ (cid:1) = g (cid:0) t , q , v , ˙ v , ρ, u ( t , q , v , ˙ v , ρ ) (cid:1) , ˜ w (cid:0) t F , q t F , v t F , ρ (cid:1) = w (cid:0) t F , q t F , v t F , ˙ v t F , ρ, u ( t F , q t F , v t F , ˙ v t F , ρ ) (cid:1) . (3)All functions are considered to be smooth. Definition 1 (Sensitivity analysis problem).
The sensitivity analysis problem is to compute the derivatives of themodel outputs (2) with respect to model parameters: d ψ d ρ : = (cid:34) d ψ d ρ · · · d ψ d ρ p (cid:35) ∈ R n c × p . (4) Definition 2 (The canonical ODE system).
To simplify the representation of the system we define the vector of3quadrature’ variables z ∈ R n c as follows: z ( t , ρ ) : = (cid:90) tt ˜ g (cid:0) τ, q , v , ρ (cid:1) d τ ⇔ ˙ z ( t , ρ ) = ˜ g (cid:0) t , q , v , ρ (cid:1) , t ≤ t ≤ t F , z ( t , ρ ) = , (5)which leads the vector of cost function (2) at final time to become: ψ = z ( t F ) + ˜ w (cid:0) t F , q t F , v t F , ρ (cid:1) . (6)Next, we add dummy evolution equations ρ (cid:48) = x ( t ) : = (cid:20) q ( t ) T v ( t ) T ρ ( t ) T z ( t ) T (cid:21) T ∈ R (2 n + p + n c ) × together with the ‘canonical ODE system’ that describes its evolution:˙ x = vf eom (cid:0) t , q , v , ρ (cid:1) p × ˜ g ( t , q , v , ρ ) : = F ( t , x ) ∈ R n + p + n c , t ≤ t ≤ t F , x ( t ) : = q ( ρ ) v ( ρ ) ρ n c × . (7) 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: Q i ( t , ρ ) : = d q ( t , ρ ) d ρ i ∈ R n , i = , . . . , p ; Q ( t , ρ ) : = (cid:20) Q ( t , ρ ) · · · Q p ( t , ρ ) (cid:21) ∈ R n × p , (8a) V i ( t , ρ ) : = d v ( t , ρ ) d ρ i ∈ R n , i = , . . . , p ; V ( t , ρ ) : = (cid:20) V ( t , ρ ) · · · V p ( t , ρ ) (cid:21) ∈ R n × p , (8b) Γ i ( t , ρ ) : = d ρ ( t , ρ ) d ρ i ∈ R p , i = , . . . , p ; Γ ( t , ρ ) : = (cid:20) Γ · · · Γ p (cid:21) = I p × p , (8c) Z i ( t , ρ ) : = ∂ z ( t , ρ ) ∂ρ i ∈ R n c , i = , . . . , p ; Z ( t , ρ ) : = (cid:20) Z ( t , ρ ) · · · Z p ( t , ρ ) (cid:21) ∈ R n c × p . (8d)The direct sensitivity for ODE systems, referred to as the Tangent Linear Model (TLM), computes the sensitivitymatrix X = (cid:104) Q T , V T , Γ , Z T (cid:105) T ∈ R (2 n + p + n c ) × p . obtained by di ff erentiating the canonical ODE system (7) with respect4o the parameters:˙ X = ˙ Q ˙ V ˙ Γ ˙ Z = Vf eom q Q + f eom v V + f eom ρ p × p ˜ g q Q + ˜ g v V + ˜ g ρ , t ≤ t ≤ t F , X ( t ) : = dq ( ρ ) d ρ dv ( ρ ) d ρ I p × p n c × p ∈ R (2 n + p + n c ) × p . (9)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 denotethe partial derivatives of f eom with respect to the subscripted variables. The detailed calculation of these expressionsand the remaining partial derivatives is explained in Appendix A. Once the sensitivities (9) have been calculated, thesensitivities of the cost functions (4) with respect to parameters are computed as follows: d ψ d ρ = Z ( t F ) + (cid:104) ˜ w q · Q + ˜ w v · V + ˜ w ρ (cid:105) t F ∈ R n c × p . (10)We note that the TLM system (9) can be written in matrix form as follows: ˙ Q ˙ V ˙ Γ ˙ Z = n × n I n × n n × p n × n c f eom q f eom v f eom ρ n × n c p × n p × n p × p p × n c ˜ g q ˜ g v ˜ g ρ n c × n c · QV Γ Z , t ≤ t ≤ t F . (11) In this section we provide the system of equations that governs the adjoint sensitivity analysis for smooth ODEsystems.
Definition 3 (Adjoint sensitivity analysis).
Apply the chain rule di ff erentiation to the total sensitivity of the costfunction (4): d ψ d ρ = d ψ d x ( t , ρ ) · d x ( t , ρ ) d ρ = λ T ( t , ρ ) · X ( t , ρ ) , (12)5here λ = ( d ψ/ d x ) T = (cid:104) λ Q T , λ V T , λ Γ T , λ Z T (cid:105) T is defined as: λ Qj ( t , ρ ) : = (cid:16) d ψ j d q ( t ,ρ ) (cid:17) T ∈ R n × , j = , . . . , n c ; λ Q ( t , ρ ) : = (cid:20) λ Q ( t , ρ ) · · · λ Qn c ( t , ρ ) (cid:21) ∈ R n × n c , (13a) λ Vj ( t , ρ ) : = (cid:16) d ψ j d v ( t ,ρ ) (cid:17) T ∈ R n × , j = , . . . , n c ; λ V ( t , ρ ) : = (cid:20) λ V ( t , ρ ) · · · λ Vn c ( t , ρ ) (cid:21) ∈ R n × n c , (13b) λ Γ j ( t , ρ ) : = (cid:16) d ψ j d ρ (cid:17) T ∈ R p × , j = , . . . , n c ; λ Γ ( t , ρ ) : = (cid:20) λ Γ ( t , ρ ) · · · λ Γ n c ( t , ρ ) (cid:21) ∈ R p × n c , (13c) λ Zj ( t , ρ ) : = (cid:16) d ψ j d z ( t ,ρ ) (cid:17) T ∈ R n c × , j = , . . . , n c ; λ Z ( t , ρ ) : = (cid:20) λ Z · · · λ Zn c (cid:21) = I n c × n c . (13d)Note that, from (5)–(6) ψ = z ( t , ρ ) + (cid:90) t F τ ˜ g (cid:0) τ, q , v , ρ (cid:1) d τ + ˜ w (cid:0) t F , q t F (cid:1) , 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 : d ψ d ρ = λ Q ( t , ρ ) T · Q ( t , ρ ) + λ V ( t , ρ ) T · V ( t , ρ ) + λ Γ ( t , ρ ) T + λ Z ( t , ρ ) T · Z ( t , ρ ) . (14)Evaluating (14) at t = t F leads to the direct sensitivity approach: d ψ d ρ = λ Q ( t F , ρ ) T · Q ( t F , ρ ) + λ V ( t F , ρ ) T · V ( t F , ρ ) + λ Γ ( t F , ρ ) T · Γ ( t F , ρ ) + λ Z ( t F , ρ ) T · Z ( t F , ρ ) . (15)By comparing this equation with (10) one obtains the values of the adjoint variables at the final time t F : λ Q ( t F , ρ ) = ˜ w T q (cid:12)(cid:12)(cid:12) t F , λ V ( t F , ρ ) = ˜ w T v (cid:12)(cid:12)(cid:12) t F , λ Γ ( t F , ρ ) = ˜ w T ρ (cid:12)(cid:12)(cid:12) t F , λ Z ( t F , ρ ) = I n c × n c . (16)The equation (14) evaluated at t = t F leads to the direct sensitivity approach: d ψ d ρ = ˜ w q (cid:12)(cid:12)(cid:12) t F · Q ( t F , ρ ) + ˜ w v (cid:12)(cid:12)(cid:12) t F · V ( t F , ρ ) + ˜ w ρ (cid:12)(cid:12)(cid:12) t F · I p × p + I n c × n c · Z ( t F , ρ ) = ˜ w q (cid:12)(cid:12)(cid:12) t F · Q ( t F , ρ ) + ˜ w v (cid:12)(cid:12)(cid:12) t F · V ( t F , ρ ) + ˜ w ρ (cid:12)(cid:12)(cid:12) t F + Z ( t F , ρ ) . (17)6valuating (14) at t = t leads to the adjoint sensitivity approach: d ψ d ρ = λ Q ( t , ρ ) T · Q ( t , ρ ) + λ V ( t , ρ ) T · V ( t , ρ ) + λ Γ ( t , ρ ) T · I p × p + λ Z ( t , ρ ) T · Z ( t , ρ ) = λ Q ( t , ρ ) T · dq ( ρ ) d ρ + λ V ( t , ρ ) T · dv ( ρ ) d ρ + λ Γ ( t , ρ ) T · I p × p + I n c × n c · n v × n c = λ Q ( t , ρ ) T · dq ( ρ ) d ρ + λ V ( t , ρ ) T · dv ( ρ ) d ρ + λ Γ ( t , ρ ) T . (18)Note that the adjoint variables are initialized at t = t F . However, their values at t = t are the ones needed forcomputing the desired sensitivities. Definition 4 (The canonical adjoint sensitivity for ODE systems).
The evolution of adjoint variables for ODE sys-tems is governed by the following continuous adjoint model: ˙ λ Q ˙ λ V ˙ λ Γ ˙ λ Z = − n × n f eom q T n × p ˜ g T q I n × n f eom v T n × p ˜ g T v p × n f eom ρ T p × p ˜ g T ρ n c × n n c × n n c × p n c × n c · λ Q λ V λ Γ λ Z , t F ≥ t ≥ t , λ ( t F , ρ ) : = ˜ w T q ( t F , ρ )˜ w T v ( t F , ρ )˜ w T ρ ( t F , ρ )I n c × n c ∈ R (2 n + p + n c ) × n c . (19)The adjoint sensitivities (19) are solved backward in time. In this study, we consider hybrid ODE systems characterized by piecewise-in-time smooth dynamics describedby (1), and that exhibit discontinuous dynamic behavior (jump or non-smoothness) in the generalized velocity statevector at a finite number of time moments (no zeno phenomena [18]). Each such moment is corresponds to an eventtriggered by the event equation: r (cid:0) q | t eve (cid:1) = , (20)where t eve is the ‘time of event’ and r : R n → R is a smooth ‘event function’. Note that grazing phenomena are notconsidered 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 ε> , ε → x ( t eve − ε ) , and its value rightafter the event by x | + t eve : = lim ε> , ε → x ( t eve + ε ) . The limits exist since the evolution of the system is smooth intime 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 . Thisis a consequence of the event changing the energy of the system by a finite amount.7
Also due to the finite energy change during the event, the quadrature variable is continuous in time, z | + t eve = z | − t eve = z | t eve . • An event that applies a finite energy impulse force to the system can abruptly change the generalized velocitystate vector ˙ q , from its value v | − t eve right before the event to a new value v | + t eve right after the event. The change invelocity is characterized by the ‘jump function’: v | + t eve = h (cid:16) t eve , q | t eve , v | − t eve , ρ (cid:17) ⇔ ˙ q | + t eve = h (cid:16) t eve , q | t eve , ˙ q | − t eve , ρ (cid:17) . (21) • An event where the system undergoes a sudden change of the equation of motions (1) at t eve is characterized bythe equations:¨ q | − t eve = f eom − (cid:0) t eve , q | t eve , v | t eve , ρ (cid:1) = : f eom − | t eve event −→ ¨ q | + t eve = f eom + (cid:0) t eve , q | t eve , v | t eve , ρ (cid:1) = : f eom + | t eve . (22) Remark 1 (Multiple events).
In many cases, the change can be triggered by one of multiple events. Each individualevent is described by the event function r (cid:96) : R n → R , (cid:96) = , . . . , e . The detection of the next event (20), which canbe one of the possible e options, is described by Π ei = r i (cid:0) q | t eve (cid:1) =
0, and if event (cid:96) takes place, then r (cid:96) = v | + t eve = h (cid:96) (cid:16) q | t eve , v | − t eve (cid:17) , or the corresponding change in the equations of motion(22) is ¨ q | + t eve = : f eom + (cid:96) | t eve . 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 theevent, respectively. Let Z | + t eve and Z | − t eve , with Z ∈ R p , be the sensitivities of the quadrature variable z ( t ) after and beforethe 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: dt eve d ρ = − drdq (cid:0) q | t eve (cid:1) · Q | − t eve drdq (cid:0) q | t eve (cid:1) · v | − t eve ∈ R × p . (23a)where dr / dq ∈ R × n is the Jacobian of the event function.8 The jump equation of the sensitivities of the generalized position state vector is: Q | + t eve = Q | − t eve − (cid:18) v | + t eve − v | − t eve (cid:19) · dt eve d ρ . (23b) • The jump equation of the sensitivities of the generalized velocity state vector is: V | + t eve = h q | − t eve · Q | − t eve + h v | − t eve · V | − t eve + (cid:16) h q | − t eve · v | − t eve − ¨ q | + t eve + h v | − t eve · ¨ q | − t eve + h t | − t eve (cid:17) · dt eve d ρ + h ρ | − t eve , (23c)where the Jacobians of the jump function are: h t | − t eve : = ∂ h ∂ t (cid:0) t eve , q | t eve , v | − t eve , ρ (cid:1) ∈ R f × , h q | − t eve : = ∂ h ∂ q (cid:0) t eve , q | t eve , v | − t eve , ρ (cid:1) ∈ R f × n , h v | − t eve : = ∂ h ∂ v (cid:0) t eve , q | t eve , v | − t eve , ρ (cid:1) ∈ R f × f , h ρ | − t eve : = ∂ h ∂ρ (cid:0) t eve , q | t eve , v | − t eve , ρ (cid:1) ∈ R f × p . (23d) • The sensitivity of the cost function changes during the event is : Z | + t eve = Z | − t eve − (cid:16) g | + t eve − g | − t eve (cid:17) · dt eve d ρ . (23e)where g | + t eve : = ˜ g (cid:0) t eve , q | t eve , v | + t eve , ρ (cid:1) , g | − t eve : = ˜ g (cid:0) t eve , q | t eve , v | − t eve , ρ (cid:1) , (23f)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 writtencompactly in matrix form as X | + t eve = S · X | − t eve , where S is the generalized sensitivity jump matrix: Q | + t eve V | + t eve Γ | + t eve Z | + t eve = (cid:16) Q | + t eve (cid:17) Q | − t eve n × n n × p n × n c (cid:16) V | + t eve (cid:17) Q | − t eve h v | − t eve h ρ | − t eve n × n c p × n p × n I p × p p × n c (cid:16) Z | + t eve (cid:17) Q | − t eve n c × n n c × p I n c × n c (cid:124) (cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32) (cid:123)(cid:122) (cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32) (cid:125) S eve · Q | − t eve V | − t eve Γ | − t eve Z | − t eve . (24a)9rom (23a) we have that: (cid:32) dt eve d ρ (cid:33) Q | − t eve = − drdq (cid:0) q | t eve (cid:1) drdq (cid:0) q | t eve (cid:1) · v | − t eve ∈ R × n . (24b)The Jacobians (cid:16) Q | + t eve (cid:17) Q | − t eve and (cid:16) V | + t eve (cid:17) Q | − t eve are: (cid:16) Q | + t eve (cid:17) Q | − t eve = I − (cid:18) v | + t eve − v | − t eve (cid:19) · (cid:32) dt eve d ρ (cid:33) Q | − teve (24c) = I + (cid:18) v | + t eve − v | − t eve (cid:19) · drdq (cid:0) q | t eve (cid:1) drdq (cid:0) q | t eve (cid:1) · v | − t eve ∈ R n × n , and (cid:16) V | + t eve (cid:17) Q | − t eve = h q | − t eve + (cid:16) h q | − t eve · v | − t eve − ¨ q | + t eve + h v | − t eve · ¨ q | − t eve + h t | − t eve (cid:17) · (cid:32) dt eve d ρ (cid:33) Q | − t eve (24d) = h q | − t eve − (cid:16) h q | − t eve · v | − t eve − ¨ q | + t eve + h v | − t eve · ¨ q | − t eve + h t | − t eve (cid:17) · drdq (cid:0) q | t eve (cid:1) drdq (cid:0) q | t eve (cid:1) · v | − t eve ∈ R n × n , (24e)respectively. The Jacobian (cid:16) Z | + t eve (cid:17) Q | − t eve is: (cid:16) Z | + t eve (cid:17) Q | − t eve = − (cid:16) g | + t eve − g | − t eve (cid:17) · (cid:32) dt eve d ρ (cid:33) Q | − t eve = (cid:16) g | + t eve − g | − t eve (cid:17) · drdq (cid:0) q | t eve (cid:1) drdq (cid:0) q | t eve (cid:1) · v | − t eve ∈ R n c × n . (24f) Theorem 1 (Adjoint sensitivity jump matrix) . Let λ Q | − t eve ∈ R n × n c , λ V | − t eve ∈ R n × n c , λ Γ | − t eve ∈ R p × n c and λ Z | − t eve ∈ R n c × n c be the adjoint sensitivities before the time of event respectively, and λ | − t eve = (cid:104) λ Q | − t eve T λ V | − t eve T λ Γ | − t eve T λ Z | − t eve T (cid:105) T . Let λ Q | + t eve ∈ R n × n c , λ V | + t eve ∈ R n × n c , λ Γ | + t eve ∈ R p × n c and λ Z | + t eve ∈ R n c × n c be the adjoint sensitivities after the time of eventrespectively, and λ | + t eve = (cid:104) λ Q | + t eve T λ V | + t eve T λ Γ | + t eve T λ Z | + t eve T (cid:105) T .The adjoint sensitivity jump equations at the time of an event are: λ | − t eve = S T eve · λ | + t eve ∈ R (2 × n + p + n c ) × n c (25)10 here 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 thesensitivity state matrix with the adjoint sensitive state matrix is constant at any time, λ | + t eve T · X | + t eve = λ | − t eve T · X | − t eve . Using(24a), the previous relationship is equivalent to λ | + t eve T · S eve · X | − t eve = λ | − t eve T · X | − t eve . Since this holds for any matrix X | − t eve it follows that λ | + t eve T · S eve = λ | − t eve T , which is equivalent to (25). Remark 2.
From (24) and (25) the adjoint sensitivity jump equations for ODE systems without constraints are: λ Q | − t eve = (cid:16) Q | + t eve (cid:17) T Q | − t eve · λ Q | + t eve + (cid:16) V | + t eve (cid:17) T Q | − t eve · λ V | + t eve + (cid:16) Z | + t eve (cid:17) T Q | − t eve · λ Z | + t eve (26a) λ V | − t eve = ( h v | − t eve ) T · λ V | + t eve (26b) λ Γ | − t eve = ( h ρ | − t eve ) T · λ V | + t eve + λ Γ | + t eve (26c) λ Z | − t eve = λ Z | + t eve (26d)
3. Sensitivity analysis for constrained multibody dynamical systems and extended cost functions
We consider constrained multibody systems that satisfy the following kinematic constraints: = Φ , (27a) = ˙ Φ = Φ q ˙ q + Φ t ⇒ Φ q v = − Φ t , (27b) = ¨ Φ = Φ q ¨ q + Φ q , q ( ˙ q , ˙ q ) + Φ t , q ˙ q + Φ t , t ⇒ Φ q ˙ v = − ( Φ q v ) v − Φ t , q v − Φ t , t : = C . (27c)Here (27a) is a holonomic position constraint equation Φ ( t , q , ρ ) = , where Φ : R + n + p → R m is a smooth ‘positionconstraint’ function. The velocity (27b) and the acceleration (27c) kinematic constraints are found by di ff erentiatingthe position constraint with respect to time. Remark 3.
Formalisms for constrained multibody systems may involve Lagrangian coe ffi cients µ : R + n + p → R m that provide the necessary forces to satisfy the kinematic constraints [17]. Our notation encompasses the case wherethe cost function penalizes the accelerations ˙ v and the joint forces via the Lagrangian coe ffi cients µ :˜ g (cid:0) t , q , v , ρ (cid:1) = g (cid:0) t , q , v , ˙ v ( t , q , v , ρ ) , ρ, µ ( t , q , v , ρ ) (cid:1) . (28)11t is shown in 3 that the terminal cost function ˜ w cannot directly depend on the acceleration ˙ v or on the Lagrangecoe ffi cients µ , and therefore the derivatives are ˜ w ˙ v = w µ =
0. In a di ff erent notation, such result is alsoshown in [19]. Using equation (10) we see that the final condition for the adjoint of the algebraic variables µ is zero, λ Λ t F = ˜ w µ (cid:12)(cid:12)(cid:12) t F = Define the extended mass matrix M : R × R n × R n × R p → R n × n and the extended right hand side function F : R × R n × R n × R p → R n as: M ( t , q , v , ρ ) : = M ( t , q , v , ρ ) + Φ T q ( t , q , v , ρ ) · α · Φ q ( t , q , v , ρ ) , (29) F ( t , q , v , ρ ) : = F ( t , q , v , ρ ) − Φ T q · α · (cid:16) ˙ Φ q v + ˙ Φ t + ξ ω ˙ Φ + ω Φ (cid:17) , (30)where α ∈ R m × m is the penalty factor of the ODE penalty formulation, ξ ∈ R and ω ∈ R are the natural frequency anddamping ratio coe ffi cients of the formulation, respectively. The functions Φ , ˙ Φ , ¨ Φ : R × R n × R n × R p → R m are theposition, velocity and acceleration kinematic constraints, respectively. The penalty formulation of a constrained rigidmultibody system is written as a first order ODE: ˙ q = v , ˙ v = f eom ( t , q , v , ρ ) = M − ( t , q , v , ρ ) · F ( t , q , v , ρ ) . (31)The Lagrange multipliers associated to the constraint forces are estimated as µ ∗ = α (cid:16) ¨ Φ + ξ ω ˙ Φ + ω Φ (cid:17) . Thesensitivities of the state variables of the system with respect to parameters evolve according to the tangent linear modelderived in [1, 20–23]. Since the penalty formulation (31) evolves as an ODE, we can compute the direct sensitivitiesusing (11) with f eom q = F q − M q ˙ v , f eom v = F v , f eom ρ = F ρ − M ρ ˙ v , as shown in Appendix A. The derivatives M ρ , M q , F q , F v , and F ρ are given in [17]. Similarly, one can compute the adjoint sensitivities of the penalty formulation (31)using (19). 12 .3. Direct and adjoint sensitivity analysis for smooth systems in the index-1 di ff erential-algebraic formulation Definition 6 (Constrained multibody dynamics: the index-1 DAE formulation).
The index-1 formulation of theequations of motion is obtained by replacing the position constraint (27a) with the acceleration constraint (27c):
I 0 00 M ( t , q , ρ ) Φ T q ( t , q , ρ ) Φ q ( t , q , ρ ) · ˙ q ˙ v µ = v F ( t , q , v , ρ ) C ( t , q , v , ρ ) , t ≤ t ≤ t F , q ( t ) = q ( ρ ) , v ( t ) = v ( ρ ) . (32)The algebraic equation has the form f DAE- µ − µ = Definition 7 (Tangent linear index-1 DAE).
Sensitivities of solutions (8) and multipliers: Λ i ( t , ρ ) : = d µ ( t , ρ ) d ρ i ∈ R m , i = , . . . , p ; (33)of the system (32) with respect to parameters evolve according to the tangent linear model derived in [1, 20–23]: ˙ Q ˙ V ˙ ΓΛ ˙ Z = Vf DAE- ˙ v q Q + f DAE- ˙ v v V + f DAE- ˙ v ρ p × p f DAE- µ q Q + f DAE- µ v V + f DAE- µ ρ (cid:0) g q + g ˙ v f DAE- ˙ v q + g µ f DAE- µ q (cid:1) · Q + (cid:0) g v + g ˙ v f DAE- ˙ v v + g µ f DAE- µ v (cid:1) · V + (cid:0) g ρ + g ˙ v f DAE- ˙ v ρ + g µ f DAE- µ ρ (cid:1) . (34)It is shown in Appendix A that equation (34) can be written in matrix form as follows: ˙ Q ˙ V ˙ ΓΛ ˙ Z = n × n I n × n n × p n × m n × n c f DAE- ˙ v q f DAE- ˙ v v f DAE- ˙ v ρ n × m n × n c p × n p × n p × p p × m p × n c f DAE- µ q f DAE- µ v f DAE- µ ρ m × m m × n c ˜ g q ˜ g v ˜ g ρ n c × m n c × n c · QV ΓΛ Z , (35)with initial conditions given by Eq. (9). Using Appendix A, the derivatives of the DAE function are: f DAE q = M Φ T q Φ q − F q − M q ˙ v − Φ T q , q µ C q − Φ q , q ˙ v , f DAE v = M Φ T q Φ q − F v C v , f DAE ρ = M Φ T q Φ q − F ρ − M ρ ˙ v − Φ Tq , ρ µ C ρ − Φ q , ρ ˙ v . efinition 8 (Continuous adjoint index-1 DAE system). The continuous adjoint di ff erential equation correspond-ing to the index-1 DAE tangent linear model (35) is: ˙ λ Q ˙ λ V ˙ λ Γ λ Λ ˙ λ Z = − n × n f DAE- ˙ v q T n × p f DAE- µ q T ˜ g T q I n × n f DAE- ˙ v v T n × p f DAE- µ v T ˜ g T v p × n f DAE- ˙ v ρ T p × p f DAE- µ ρ T ˜ g T ρ m × n m × n m × p m × m m × n c n c × n n c × n n c × p n c × m n c × n c · λ Q λ V λ Γ λ Λ λ Z , t F ≥ t ≥ t , λ ( t F , ρ ) : = ˜ w T q ( t F , ρ )˜ w T v ( t F , ρ )˜ w T ρ ( t F , ρ )0 m × n c I n c × n c ∈ R (2 n + p + m + n c ) × n c . (36)Noting from Remark 3 that the algebraic equation in (36) reads: λ Λ ( t ) = , t F ≥ t ≥ t , the index-1 adjoint DAE (36) can be reduced to the following adjoint ODE: ˙ λ Q ˙ λ V ˙ λ Γ ˙ λ Z = − n × n f DAE- ˙ v q T n × p ˜ g T q I n × n f DAE- ˙ v v T n × p ˜ g T v p × n f DAE- ˙ v ρ T p × p ˜ g T ρ n c × n n c × n n c × p n c × n c · λ Q λ V λ Γ λ Z , t F ≥ t ≥ t , λ ( t F , ρ ) : = ˜ w T q ( t F , ρ )˜ w T v ( t F , ρ )˜ w T ρ ( t F , ρ )I n c × n c ∈ R (2 n + p + n c ) × n c . (37) We now discuss constrained dynamical systems when the dynamics is piecewise smooth in time. Performing asensitivity analysis for a constrained rigid hybrid multibody dynamic system requires finding the jump conditions atthe time of event. These jump equations are explained in our previous work [17]. We summarize below the jumpequations 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 bothconstraint functions Φ − | − t eve : = Φ − (cid:0) t eve , q | t eve , ρ (cid:1) = , and Φ + | + t eve : = Φ + (cid:0) t eve , q | t eve , ρ (cid:1) = . • The velocity state variables jump from their values right before the event to right after the event according tothe jump equation: v dof + | + t eve = h (cid:16) t eve , q | t eve , v dof- | − t eve , ρ (cid:17) , h : R + n + f − + p → R f + . (38)The jump function (38) is assumed to be smooth and defined in terms of the velocity degrees of freedom (theindependent components). 14 The jumps in velocity cannot be arbitrary for the dependent components. They are dependent of the degree offreedom and are obtained from solving the velocity constraints leading to: v dep + | + t eve = − (cid:16) Φ + q dep + | + t eve (cid:17) − · (cid:16) Φ + q dof + | + t eve v dof + | + t eve + Φ + t | + t eve (cid:17) = R + | + t eve v dof + | + t eve − (cid:16) Φ + q dep + | + t eve (cid:17) − · Φ + t | + t eve . (39)Where R ± corresponds to the null space of the constraints if the constraints are scleronomic (non explicitly timedependent).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 velocitystate. The impulsive (external) contact forces act to change the DOF components without changing the set ofconstraint equations, Φ + ≡ Φ − . • The case where the event consists solely of an inelastic collisions and a change of constraints Φ + (cid:44) Φ − , withoutany external force modifying the independent velocities. The impulsive (internal) constraints forces at the timeof event are solved by using a popular approach in robotics [24]: M | t eve ( Φ + q ) T | t eve Φ + q | t eve · v | + t eve δµ = M | t eve · v | − t eve − Φ + t | t eve , (40a)or, equivalently, v | + t eve δµ = M | t eve ( Φ + q ) T | t eve Φ + q | t eve − · M | t eve · v | − t eve − Φ + t | t eve = f DAE-imp- v (cid:16) t eve , q | t eve , v | − t eve , ρ (cid:17) f DAE-imp- µ (cid:16) t eve , q | t eve , v | − t eve , ρ (cid:17) . (40b)The second equation (40b) imposes the velocity constraint on both independent and dependent coordinates,which is covered by our formalism as: v dof + | + t eve = P dof + f DAE-imp- v (cid:16) t eve , q | t eve , v | − t eve , ρ (cid:17) = : h (cid:16) t eve , q | t eve , v dof − | − t eve , ρ (cid:17) , (41)where P = P dep P dof is a permutation matrix that partitions the state variables into dependent and independentvariables.Finally, the jump conditions at the time of event in the sensitivity state matrix are:15 The independent components of the sensitivity of the generalized positions right after the event: Q dof + | + t eve = Q dof + | − t eve − (cid:18) v dof + | + t eve − v dof + | − t eve (cid:19) · dt eve d ρ . (42a)which are equivalent to: P + dof + · (cid:18) Q | + t eve − Q | − t eve (cid:19) = − P + dof + · (cid:18) v | + t eve − v | − t eve (cid:19) · dt eve d ρ . (42b) • The dependent components of the sensitivity of the generalized positions right after the event: Q dep + | + t eve = R + | + t eve · Q dof + | + t eve − (cid:16) Φ + q dep + | + t eve (cid:17) − Φ + ρ (cid:12)(cid:12)(cid:12)(cid:12) + t eve . (42c) • The independent coordinates of the velocity sensitivities right after the event, V dof + | + t eve = h q | − t eve · Q | − t eve + h v dof- | − t eve · V dof- | − t eve (43) + (cid:16) h q | − t eve · v | − t eve − ¨ q dof + | + t eve + h v dof- | − t eve · ¨ q dof- | − t eve + h t | − t eve (cid:17) · dt eve d ρ + h ρ | − t eve , where the Jacobians of the jump function are: h q | − t eve : = ∂ h ∂ q (cid:0) q | t eve , v dof- | − t eve , ρ (cid:1) ∈ R f + × n , h v dof- | − t eve : = ∂ h ∂ v dof- (cid:0) q | t eve , v dof- | − t eve , ρ (cid:1) ∈ R f + × f − . h t | − t eve : = ∂ h ∂ t (cid:0) q | t eve , v dof- | − t eve , ρ (cid:1) ∈ R f , h ρ | − t eve : = ∂ h ∂ρ (cid:0) q | t eve , v dof- | − t eve , ρ (cid:1) ∈ R f × p . • The dependent components of the velocity sensitivities right after the event, V dep + | + t eve = − (cid:0) Φ + q dep + | + t eve (cid:1) − (cid:16) Φ + q dof + · V dof + + (cid:0) Φ + q , q v + Φ + t , q (cid:1) · Q + Φ + q , ρ v + Φ + t , ρ (cid:17)(cid:12)(cid:12)(cid:12)(cid:12) + t eve . (44) Definition 9 (The generalized sensitivity jump matrix for elastic impact).
The jump equations (38)–(44) for con-strained systems can be written compactly in matrix form as a jump of the state sensitivity matrix X at the time of the16vent, X | + t eve = S eve · X | − t eve , where S eve represents the generalized jump sensitivity matrix: Q dep + | + t eve Q dof + | + t eve V dep + | + t eve V dof + | + t eve Γ | + t eve Z | + t eve = (cid:16) Q dep + | + t eve (cid:17) Q | − t eve ( n − f ) × ( n − f ) ( n − f ) × f D ( n − f ) × n c (cid:16) Q dof + | + t eve (cid:17) Q | − t eve f × ( n − f ) f × f f × p f × n c (cid:16) V dep + | + t eve (cid:17) Q | − t eve ( n − f ) × ( n − f ) (cid:16) V dep + | + t eve (cid:17) V dof + | − t eve K ( n − f ) × n c (cid:16) V dof + | + t eve (cid:17) Q | − t eve f × ( n − f ) h v | − t eve h ρ | − t eve f × n c p × n p × ( n − f ) p × f I p × p p × n c (cid:16) Z | + t eve (cid:17) Q | − t eve n c × ( n − f ) n c × f n c × p I n c × n c (cid:124) (cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32) (cid:123)(cid:122) (cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32) (cid:125) S eve · Q dep + | − t eve Q dof + | − t eve V dep + | − t eve V dof + | − t eve Γ | − t eve Z | − t eve . (45)The Jacobians of the jump equations with respect to the sensitivity state before the time of event are: (cid:16) Q dof + | + t eve (cid:17) Q | − t eve = P + dof + I n × n − (cid:18) v | + t eve − v | − t eve (cid:19) · (cid:32) dt eve d ρ (cid:33) Q | − teve ( P − ) T ∈ R f × n , (46)with (cid:32) dt eve d ρ (cid:33) Q | − t eve = − drdq (cid:0) q | t eve (cid:1) drdq (cid:0) q | t eve (cid:1) · v | − t eve ∈ R × n . (47)It follows that: (cid:16) Q dep + | + t eve (cid:17) Q | − t eve = R + | + t eve · (cid:16) Q dof + | + t eve (cid:17) Q | − t eve ∈ R ( n − f ) × n , (48)and (cid:16) V dof + | + t eve (cid:17) Q | − t eve = h q | − t eve + (cid:16) h q | − t eve · v | − t eve − ¨ q dof + | + t eve + h v dof- | − t eve · ¨ q dof- | − t eve + h t | − t eve (cid:17) · (cid:32) dt eve d ρ (cid:33) Q | − t eve · ( P − ) T ∈ R f × n . (49)Rewriting (44) as: (cid:16) V dep + | + t eve (cid:17) = (cid:16) R + · V dof + + R + · Q + C (cid:17)(cid:12)(cid:12)(cid:12)(cid:12) + t eve ∈ R ( n − f ) × p (50)17r, equivalently, as: (cid:16) V dep + | + t eve (cid:17) = R + · V dof + + R + · ( P − ) T · Q dep + | − t eve Q dof + | − t eve + C (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) + t eve , R + | + t eve = − (cid:0) Φ + q dep + (cid:1) − (cid:16) Φ + q , q v + Φ + t , q (cid:17)(cid:12)(cid:12)(cid:12)(cid:12) + t eve , C = − (cid:0) Φ + q dep + (cid:1) − (cid:16) Φ + q , ρ v + Φ + t , ρ (cid:17)(cid:12)(cid:12)(cid:12)(cid:12) + t eve ∈ R ( n − f ) × p , (51)we find the following expressions for the Jacobians: (cid:16) V dep + | + t eve (cid:17) Q | − t eve = R + | + t eve · (cid:16) V dof + | + t eve (cid:17) Q | − t eve + R + | + t eve · ( P − ) T R ( n − f ) × n , (cid:16) V dep + | + t eve (cid:17) V dof + | − t eve = R + | + t eve · h v | − t eve . (52)The expressions for D and K in (45) are: D = − (cid:16) Φ + q dep + (cid:17) − Φ + ρ (cid:12)(cid:12)(cid:12)(cid:12) + t eve ∈ R ( n − f ) × p , K = C + R + · h ρ | − t eve + R + · ( P − ) T · D f × p (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) + t eve ∈ R ( n − f ) × p . (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 veloc-ity sensitivity for constrained systems due to impulsive forces, presented in [17], is determined as follows: V | + t eve δ Λ = − M Φ + q T Φ + q − M q · ( v | + t eve − v | − t eve ) + Φ + q , q T · δλ Φ + q , q · v | + t eve · Q | + t eve + M Φ + q T Φ + q − M0 · V | − t eve − M Φ + q T Φ + q − M ρ · v | + t eve + Φ + q , ρ T · δλ Φ + q , ρ · v | + t eve + Φ + t , ρ · v | − t eve − M Φ + q T Φ + q − Φ + t , q · v | − t eve + Φ + t , v · v | − t eve , (54)which simplifies to: V | + t eve δ Λ = f DAE-imp q · Q | + t eve + f DAE-imp v | − t eve · V | − t eve + f DAE-imp ρ + f DAE-imp t . (55)Thus, the jump the velocity state variables at the time of event is V | + t eve = f DAE-imp- v q | + t eve Q | + t eve + f DAE-imp- v v | + t eve V | + t eve + f DAE-imp- v ρ | + t eve , (56)18nd the jump in the sensitivity of the Lagrange multipliers from Λ | − t eve → Λ | + t eve is: Λ | + t eve = Λ | − t eve + f DAE-imp- µ q | + t eve Q | + t eve + f DAE-imp- µ v | + t eve V | + t eve + f DAE-imp- µ ρ | + t eve . (57)The corresponding sensitivity jump matrix (45) is: Q dep + | + t eve Q dof + | + t eve V | + t eve Λ | + t eve Γ | + t eve Z | + t eve = (cid:16) Q dep + | + t eve (cid:17) Q | − t eve ( n − f ) × n ( n − f ) × m D ( n − f ) × n c (cid:16) Q dof + | + t eve (cid:17) Q | − t eve f × n f × m f × p f × n c f DAE-imp- v q | + t eve f DAE-imp- v v | + t eve n × m f DAE-imp- v ρ | + t eve n × nc f DAE-imp- µ q | + t eve f DAE-imp- µ v | + t eve m × m f DAE-imp- µ ρ | + t eve m × nc p × n p × n p × m I p × p p × n c (cid:16) Z | + t eve (cid:17) Q | − t eve n c × n n c × m n c × p I n c × n c (cid:124) (cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32) (cid:123)(cid:122) (cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32) (cid:125) S eve · Q dep + | − t eve Q dof + | − t eve V | − t eve Λ | − t eve Γ | − t eve Z | − t eve . (58) Definition 11 (Jump in adjoint sensitivity for constrained systems with elastic impact).
The transpose of the di-rect sensitivity jump matrix S Teve (25) associated with an elastic impact (45) is: S Teve = (cid:16) Q dep + | + t eve (cid:17) T Q | − t eve (cid:16) Q dof + | + t eve (cid:17) T Q | − t eve (cid:16) V dep + | + t eve (cid:17) T Q | − t eve (cid:16) V dof + | + t eve (cid:17) T Q | − t eve n × p (cid:16) Z | + t eve (cid:17) T Q | − t eve ( n − f ) × ( n − f ) ( n − f ) × f ( n − f ) × ( n − f ) ( n − f ) × f ( n − f ) × p ( − n − f ) × n c f × ( n − f ) f × f (cid:16) V dep + | + t eve (cid:17) T V dof + | − t eve (cid:104) h v | − t eve (cid:105) T f × p f × p D T p × f K T (cid:104) h ρ | − t eve (cid:105) T I p × p p × n c n c × ( n − f ) n c × f n c × ( n − f ) n c × f n c × p I n c × n c (59)From the adjoint sensitivity equation (25) the jumps in adjoint variables for ODE systems with constraints undergoingan elastic impact are: λ Q | − t eve = (cid:20)(cid:16) Q dep + | + t eve (cid:17) T Q | − t eve (cid:16) Q dof + | + t eve (cid:17) T Q | − t eve (cid:21) · λ Q | + t eve + (cid:20)(cid:16) V dep + | + t eve (cid:17) T Q | − t eve (cid:16) V dof + | + t eve (cid:17) T Q | − t eve (cid:21) · λ V | + t eve + (cid:16) Z | + t eve (cid:17) T Q | − t eve · λ Z | + t eve λ V | − t eve = ( n − f ) × f ( n − f ) × f (cid:16) V dep + | + t eve (cid:17) T V dof + | − t eve (cid:104) h v | − t eve (cid:105) T · λ V | + t eve λ Γ | − t eve = (cid:20) D T p × f (cid:21) · λ Q | + t eve + (cid:20) K T h ρ | − t eve T (cid:21) · λ V | + t eve + λ Γ | + t eve λ Z | − t eve = λ Z | + t eve (60)19 efinition 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 Teve (25) associated with (58) is: S Teve = (cid:16) Q dep + | + t eve (cid:17) T Q | − t eve (cid:16) Q dof + | + t eve (cid:17) T Q | − t eve (cid:104) f DAE-imp- v q | + t eve (cid:105) T (cid:104) f DAE-imp- µ q | + t eve (cid:105) T n × p (cid:16) Z | + t eve (cid:17) T Q | − t eve n × ( n − f ) n × f (cid:104) f DAE-imp- v v | + t eve (cid:105) T (cid:104) f DAE-imp- µ v | + t eve (cid:105) T n × ( p ) n × ( n c ) m × ( n − f ) m × f m × n m × m m × ( p ) m × ( n c ) D T p × f (cid:104) f DAE-imp- v ρ | + t eve (cid:105) T (cid:104) f DAE-imp- µ ρ | + t eve (cid:105) T I p × p p × n c n c × ( n − f ) n c × f n c × n n c × m n c × ( p ) I n c × n c (61)Since the adjoints of the algebraic Lagrange variables are zero (Remark 3), the adjoint sensitivity equations (25)provide the jump equations for adjoint variables at the time of event: λ Q | − t eve = (cid:20)(cid:16) Q dep + | + t eve (cid:17) T Q | − t eve (cid:16) Q dof + | + t eve (cid:17) Q | − t eve (cid:21) · λ Q | + t eve + (cid:104) f DAE-imp- v q | + t eve (cid:105) T · λ V | + t eve + (cid:16) Z | + t eve (cid:17) T Q | − t eve · λ Z | + t eve ,λ V | − t eve = (cid:104) f DAE-imp- v v | + t eve (cid:105) T · λ V | + t eve ,λ Γ | − t eve = (cid:20) D T p × f (cid:21) · λ Q | + t eve + (cid:104) f DAE-imp- v ρ | + t eve (cid:105) T · λ V | + t eve + λ Γ | + t eve ,λ Z | − t eve = λ Z | + t eve . (62) Remark 4 (Sensitivities of the cost function).
Once the evolution of the sensitivities of the direct or adjoint sensi-tivities are computed, the sensitivity of the cost function with respect to parameters d ψ/ d ρ is obtained from equations(15) and (18). Note that the evolution of the direct and adjoint sensitivities involve is piecewise continuous in time,with jumps occurring at each event.
4. Case study: sensitivity analysis of a five-bar mechanism
The five-bar mechanism, presented in Fig. 1a, is used as a case study to validate the adjoint sensitivity methodin computing the sensitivity of cost functions with respect to parameters for hybrid constrained dynamical systems.The mechanism has two degrees of freedom, five revolute joints located at points A, 1, 2, 3, and B; the masses ofeach bars are m = kg , m = . kg , m = . kg , m = kg ; the polar moments of inertia for each bars haveuniformly distributed mass; the two springs have sti ff ness coe ffi cients of k = k = N / m and natural lengths of L = . m and L = . m . The state vector q = (cid:2) q T1 q T2 q T3 (cid:3) T includes the natural coordinates of the point1, 2, and 3 of the mechanism. The coordinates q = (cid:2) x y (cid:3) T are independent and defines the DOF of the system,while the coordinates q = (cid:2) x y (cid:3) T with q = (cid:2) x y (cid:3) T are dependent. The constraint equations, used to solve for the20ependent coordinates, are defined according to the fixed lengths between each set of points, as follows: Φ = (cid:107) q A − q (cid:107) − L A (cid:107) q − q (cid:107) − L (cid:107) q − q (cid:107) − L (cid:107) q B − q (cid:107) − L B = , (63)with the lengths L A = L B = . m ; L A = L B = . m and the ground points q A = (cid:2) − . (cid:3) T ; q B = (cid:2) . (cid:3) T . The study focuses on point 2 of the five-bar mechanism. This points hits the ground at -2.35 malong the vertical y axis which is detected by an the event function r ( · ) described in Eq. (20). At the time of event,the vertical velocity of point 2 jumps to its opposite value, while its horizontal velocity remains the same. The ODEforward system, the direct and adjoint sensitivity are simulated with a time span of five seconds. The residuals of theconstraint equations, presented in Fig. 1b, shows that the position and the velocity constraints are satisfied within anerror of 10 − and 10 − , respectively, which is satisfactory. (a) Diagram of the five-bar mechanism. Time, t [s] -2024681012 P o s i t i o n C o n s t r a i n t s , Φ ,[ m ] × -6 -2024681012 V e l o c i t y C o n s t r a i n t s ,˙ Φ ,[ m / s ] × -6 PositionVelocity (b) The position and the velocity constraintresiduals for the five-bar mechanism.
Figure 1: Structure of the five-bar mechanism.
The trajectories of the vertical position and velocity of point 2 of the five-bar mechanism are presented in Fig. 2aand Fig. 2b, respectively. As expected, the point 2’s vertical position bounces at -2.35m, and its vertical velocityjumps at each time of event with v | + t eve = v | − t eve . 21 Time, t [s] -2.5-2.4-2.3-2.2-2.1-2 P o s i t i o n o f p o i n t , y , [ m ] (a) The position of point 2. Time, t [s] -2-1.5-1-0.500.511.52 V e l o c i t y o f p o i n t , ˙ y , [ m / s ] (b) The velocity of point 2. Figure 2: The state variables of point 2 of the five-bar mechanism.
The trajectory of the quadrature variable z ( t ) = (cid:82) tt ˙ y ( τ ) d τ of the five-bar mechanism and its sensitivity are shownin Fig. 3a and Fig. 3b, respectively. The same analysis is presented in Fig. 4a and Fig. 4b for the quadrature variable z ( t ) = (cid:82) tt ¨ y d τ . The direct sensitivity is represented by the continuous line, while the central finite di ff erence sensitivityis represented by the dashed line. Both solutions were solved forward in time. The adjoint sensitivity is presentedas well, and was solved backwards in time. As presented in our previous paper, the direct di ff erentiation methodto compute the sensitivity of the cost function with discontinuities in the velocity state variables of the mechanism isvalidated. 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 di ff erence method. One main conclusionof our previous paper was to state that our proposed direct sensitivity method in computing the sensitivity of the costfunction with discontinuities in the velocity state variables was more robust than the numerical method. Indeed, thedirect method accurately determines the jump in the sensitivities and their trajectories. This after each event, withoutany delta-like jumps in magnitude 1 /ε that occurs in the numerical method at each time of event. This validated directsensitivity method is now compared to the proposed adjoint method in computing the sensitivity of the cost functionwith discontinuities in the velocity state variables. The results presented in Fig. 3a and Fig. 4a show that the adjointand direct method exactly converge to the same sensitivity cost number with a di ff erence of less than 0.01 %. Thisconvergence in both methods validates the adjoint sensitivity method in computing the sensitivity of the cost functionwith discontinuities in the trajectories. 22 Time, t [s] -30-20-100102030 D i r ec t s e n s i t i v i t y m e t h o d , Z ( t ) -30-20-100102030 A d j o i n t s e n s i t i v i t y m e t h o d , λ Γ ( t ) Proposed formulationsCentral Finite dif.
X: 5Y: 12.31X: 0Y: 12.34 (a) Direct and adjoint sensitivities.
Time, t [s] -0.3-0.25-0.2-0.15-0.1-0.050 T h e q u a d r a t u r e v a r i a b l e , z ( t ) (b) The quadrature variable z ( t ). Figure 3: Sensitivity analysis of the five-bar mechanism with z ( t ) = (cid:82) tt ˙ y ( τ ) d τ . Time, t [s] -300-200-1000100200 D i r ec t s e n s i t i v i t y m e t h o d , Z ( t ) -300-200-1000100200 A d j o i n t s e n s i t i v i t y m e t h o d , λ Γ ( t ) Proposed formulationsCentral Finite dif.
X: 5Y: -319.7X: 0Y: -320.6 (a) Direct and adjoint sensitivities.
Time, t [s] -14-12-10-8-6-4-202 T h e q u a d r a t u r e v a r i a b l e , z ( t ) (b) The quadrature variable z ( t ). Figure 4: Sensitivity analysis of the five-bar mechanism with z ( t ) = (cid:82) tt ¨ y d τ . Note that z ( t ) = (cid:82) tt ¨ y ( τ ) d τ does not completely match the trajectory of the velocity of point 2 in Fig. 2b. Indeed,the point2’s velocity jumps at the time of event, while the quadrature variable does not. The quadrature variableevaluates the integral of the acceleration of point 2 only.The same analysis is provided with the quadrature variable z ( t ) = (cid:82) tt ¨ y ( τ ) + ˙ y ( τ ) d τ in Fig. 5a and Fig. 5b. Theadjoint and direct method converge to the same sensitivity cost number with a di ff erence of less than 0.01 %.23 Time, t [s] -100001000200030004000 D i r ec t s e n s i t i v i t y m e t h o d , Z ( t ) -100001000200030004000 A d j o i n t s e n s i t i v i t y m e t h o d , λ Γ ( t ) Proposed formulationsCentral Finite dif.
X: 5Y: 1703X: 0Y: 1721 (a) Direct and adjoint sensitivities.
Time, t [s] T h e q u a d r a t u r e v a r i a b l e , z ( t ) (b) The quadrature variable z ( t ). Figure 5: Sensitivity analysis of the five-bar mechanism with z ( t ) = (cid:82) tt ¨ y ( τ ) + ˙ y ( τ ) d τ .
5. Conclusions
Gradient based algorithms are widely used in computational engineering problems such as design and controloptimization, implicit time integration methods, and deep learning. Sensitivity analysis plays a key role in this typeof algorithms as it provides the necessary derivative information. In the context of dynamical systems governed byordinary or di ff erential algebraic equations, sensitivity analysis computes the derivatives of general cost functions thatdepend 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 theliterature [1, 20–23]. Our earlier work has extended the direct sensitivity analysis to hybrid multibody dynamicsystems 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 dynamicsystems modeled by ordinary di ff erential equations and by index-1 di ff erential algebraic equations. A very generalformulation of the cost functions is used. For the hybrid systems considered herein discontinuities in the forward tra-jectories appear at time moments triggered by an event. Jump conditions for adjoint sensitivity variables are providedfor mechanical systems with and without constraints. These jump conditions handle the change in the sensitivitiescaused 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 afive-bar mechanism with non-smooth contacts. The direct and adjoint sensitivities computed by the proposed math-ematical framework are validated against numerical sensitivities calculated by real finite di ff erences. The results ofthis study show that all the alternative analyses provide the same sensitivities of the general cost function with respect24o model parameters, within an error of 0.01%.Future work will extend the mathematical framework to calculate adjoint sensitivities of hybrid mechanical sys-tems with respect to actuation functions. These sensitivities will allow to solve optimal control problems for hybridsystems. A. Calculation of partial derivatives used in sensitivity analysesRemark 5.
The expressions f eom q , f eom v , and f eom ρ i denote the partial derivatives of f eom with respect to the subscriptedvariables. The partial derivatives ∂ f eom /∂ζ are obtained by di ff erentiating f eom with respect to ζ ∈ { q , v , ρ } : ∂ f eom ∂ζ = ∂ ( M − F ) ∂ζ = − M − M ζ M − F + M − F ζ = M − (cid:16) F ζ − M ζ f eom (cid:17) = M − (cid:16) F ζ − M ζ ˙ v (cid:17) . (A.1) Remark 6.
The expressions ˜ g q , ˜ g v , and ˜ g ρ i denote the partial derivatives of ˜ g with respect to the subscripted variables.The partial derivatives ∂ ˜ g /∂ζ are obtained by di ff erentiating (1) with respect to ζ ∈ { q , v , ρ } :˜ g ζ = g ζ + g ˙ v f eom ζ + g ˜ u ˜ u ζ (A.2) = g ζ + g ˙ v f eom ζ + g u u ζ + g u u ˙ v f eom ζ , which leads to: (cid:2) ˜ g q Q i + ˜ g v V i + ˜ g ρ i (cid:3) i = ,..., p = (cid:2)(cid:0) g q + g ˙ v f eom q + g u u q + g u u ˙ v f eom q (cid:1) · Q i + (cid:0) g v + g ˙ v f eom v + g u u v + g u u ˙ v f eom v (cid:1) · V i (A.3) + g ρ i + g ˙ v · f eom ρ i + g u u ρ + g u u ˙ v f eom ρ i (cid:3) i = ,..., p . Remark 7.
Similarly, the expressions ˜ w q , ˜ w v , and ˜ w ρ i denote the partial derivatives of ˜ w with respect to the subscriptedvariables. The partial derivatives ∂ ˜ w /∂ζ are obtained by di ff erentiating w with respect to ζ ∈ { q , v , ρ } :˜ w ζ = w ζ + w ˙ v f eom ζ + w ˜ u ˜ u ζ (A.4) = w ζ + w ˙ v f eom ζ + w u u ζ + w u w ˙ v f eom ζ . B. Adjoint of the algebraic Lagrangian coe ffi cient Methods to compute the adjoint of an index-1 DAE available in the literature [19, 25, 26] use the followingapproach. Define the Lagrangian using the multipliers µ Q , µ V , µ Γ that correspond to the constraints posed by the25ndex-1 DAE equations (32): µ Q µ V µ Γ T · I 0 00 M ( t , q , ρ ) Φ T q ( t , q , ρ ) Φ q ( t , q , ρ ) · ˙ q ˙ v µ − v F ( t , q , v , ρ ) C ( t , q , v , ρ ) . (B.1)We rearrange equation (B.1) as follows: µ Q µ V µ Γ · I 0 00 M ( t , q , ρ ) Φ T q ( t , q , ρ ) Φ q ( t , q , ρ ) · ˙ q ˙ v µ − I 0 00 M ( t , q , ρ ) Φ T q ( t , q , ρ ) Φ q ( t , q , ρ ) − · v F ( t , q , v , ρ ) C ( t , q , v , ρ ) (B.2) = µ Q µ V µ Γ T · I 0 00 M ( t , q , ρ ) Φ T q ( t , q , ρ ) Φ q ( t , q , ρ ) · ˙ q ˙ v µ − Vf DAE- ˙ v f DAE- µ (B.3) = λ Q λ V λ Γ T · ˙ q ˙ v µ − Vf DAE- ˙ v f DAE- µ . (B.4)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: λ Q λ V λ Γ = I 0 00 M ( t , q , ρ ) Φ T q ( t , q , ρ ) Φ q ( t , q , ρ ) · µ Q µ V µ Γ . The adjoint DAE equations and boundary conditions in the “ µ formulation” [19, 25, 26] can be derived from theequations and boundary conditions in the “ λ formulation” discussed in this paper, and vice-versa.26 omenclatureDimensions n The number of generalized coordinates p The number of parameters n c The number of cost functions m The number of equations of constraints
Dynamics f eom The function solving the Equation of Motion ∈ R × R n × R n × R p → R n q , ˙ q ∈ R n The generalized position and velocity state vector¨ q ∈ R n The generalized acceleration state vector z ∈ R n c The vector of quadrature variables x ∈ R (2 n + p + n c ) The state vector of the canonical ODE t eve ∈ R The time of event ρ ∈ R p The vector of system parameters F The generalized force vector ∈ R × R n × R n × R p → R n M The generalized smooth and invertible Mass matrix ∈ R × R n × R p → R n × n Φ The equations of constraints ∈ R × R n × R p → R m General ˙ (cid:3) or ¨ (cid:3) The total (first or second order) derivative of a function or variable with respect to time (cid:3) ζ,φ
Double subscripts indicates a three-dimensional Jacobian with respect to a quantity ζ and φ ,unless stated otherwise (cid:3) ζ 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 vector ˙ q with respect to the vector of system parameters ρ X ∈ R (2 n + p + n c ) × p The sensitivity matrix of the x state vector with respect to the vector of system parameters ρ dt eve / d ρ ∈ R × p The sensitivity of the time of event t eve with respect to the vector of system parameters ρλ ∈ R (2 n + 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 cknowledgments This project has been partially funded by the European Union Horizon 2020 Framework Program, Marie SkłodowskaCurie actions, under grant agreement no. 645736, Project EVE, Innovative Engineering of Ground Vehicles with in-tegrated Active Chassis Systems. It was also supported in part by awards NSF DMS–1419003, NSF CCF–1613905,NSF ACI–1709727, AFOSR DDDAS 15RT1037, by the Computational Science Laboratory.
ReferencesReferences [1] Y. Zhu, D. Dopico, C. Sandu, A. Sandu, Dynamic response optimization of complex multibody systems in a penalty formulation using adjointsensitivity, ASME. J. Comput. Nonlinear Dynam. 10 (3) (2015) 031009. doi:10.1115/1.4029601 .[2] K.-H. Chang, Chapter 18 - Structural Design Sensitivity Analysis, Academic Press, Boston, 2015. doi:http://dx.doi.org/10.1016/B978-0-12-382038-9.00018-1 .[3] A. Sandu, D. N. Daescu, G. R. Carmichael, Direct and adjoint sensitivity analysis of chemical kinetic systems with kpp: Part i—theory andsoftware tools, Atmospheric Environment 37 (36) (2003) 5083 – 5096. doi:http://dx.doi.org/10.1016/j.atmosenv.2003.08.019 .[4] P. I. Barton, R. J. Allgor, W. F. Feehery, S. Gal´an, Dynamic optimization in a discontinuous world, Industrial & Engineering ChemistryResearch 37 (3) (1998) 966–981. arXiv:http://dx.doi.org/10.1021/ie970738y , doi:10.1021/ie970738y .URL http://dx.doi.org/10.1021/ie970738y [5] S. Gal´an, W. F. Feehery, P. I. Barton, Parametric sensitivity functions for hybrid discrete / continuous systems, Applied Numerical Mathematics31 (1) (1999) 17 – 47. doi:http://dx.doi.org/10.1016/S0168-9274(98)00125-1 .URL [6] P. I. Barton, C. K. Lee, Modeling, simulation, sensitivity analysis, and optimization of hybrid systems, ACM Trans. Model. Comput. Simul.12 (4) (2002) 256–289. doi:10.1145/643120.643122 .[7] J. E. Tolsma, P. I. Barton, Hidden discontinuities and parametric sensitivity calculations, SIAM Journal on Scientific Computing 23 (6) (2002)1861–1874. doi:10.1137/S106482750037281X .[8] E. Rozenvasser, General sensitivity equations of discontinuous systems, Automatika i telemekhanika 3 (1967) 52–56.[9] A. Saccon, N. van de Wouw, H. Nijmeijer, Sensitivity analysis of hybrid systems with state jumps with application to trajectory tracking,53rd IEEE Conference on Decision and Control (2014) 3065–3070 doi:10.1109/CDC.2014.7039861 .[10] I. A. Hiskens, J. Alseddiqui, Sensitivity, approximation, and uncertainty in power system dynamic simulation, IEEE Transactions on PowerSystems 21 (4) (2006) 1808–1820. doi:10.1109/TPWRS.2006.882460 .[11] I. A. Hiskens, M. A. Pai, Trajectory sensitivity analysis of hybrid systems, IEEE Transactions on Circuits and Systems I: Fundamental Theoryand Applications 47 (2) (2000) 204–220. doi:10.1109/81.828574 .[12] F. Taringoo, P. Caines, On the geometry of switching manifolds for autonomous hybrid systems, IFAC Proceedings Volumes 43 (12) (2010)35 – 40. doi:http://dx.doi.org/10.3182/20100830-3-DE-4013.00008 .[13] W. Backer, Jump conditions for sensitivity coe ffi cients, Sensitivity methods in control theory (Symp. Dubrovnik 1964; L. Radanovi´c, ed.)(1964) pp. 168–175.
14] D. E. Stewart, M. Anitescu, Optimal control of systems with discontinuous di ff erential equations, Numerische Mathematik 114 (4) (2010)653–695. doi:10.1007/s00211-009-0262-2 .URL http://dx.doi.org/10.1007/s00211-009-0262-2 [15] F. Taringoo, P. E. Caines, The sensitivity of hybrid systems optimal cost functions with respect to switching manifold parameters, in: R. Ma-jumdar, P. Tabuada (Eds.), Hybrid Systems: Computation and Control, Springer Berlin Heidelberg, Berlin, Heidelberg, 2009, pp. 475–479.[16] H. Zhang, S. Abhyankar, E. Constantinescu, M. Anitescu, Discrete adjoint sensitivity analysis of hybrid dynamical systems with switching,IEEE Transactions on Circuits and Systems I: Regular Papers 64 (2017) 13. doi:10.1109/TCSI.2017.2651683 .[17] S. Corner, C. Sandu, A. Sandu, Modeling and sensitivity analysis methodology for hybrid dynamical systems, arXiv preprintarXiv:1710.04292.[18] A. M. Pace, S. A. Burden, Piecewise - di ff erentiable trajectory outcomes in mechanical systems subject to unilateral constraints, in: Proceed-ings of the 20th International Conference on Hybrid Systems: Computation and Control, HSCC ’17, ACM, New York, NY, USA, 2017, pp.243–252. doi:10.1145/3049797.3049807 .URL http://doi.acm.org/10.1145/3049797.3049807 [19] D. Dopico, A. Sandu, C. Sandu, Y. Zhu, Sensitivity analysis of multibody dynamic systems modeled by odes and daes., In: Terze Z. (eds)Multibody Dynamics. Computational Methods in Applied Sciences, vol 35. doi:978-3-319-07260-9 .URL https://doi.org/10.1007/978-3-319-07260-9_1 [20] D. Dopico, A. Sandu, Y. Zhu, C. Sandu, Direct and adjoint sensitivity analysis of ordinary di ff erential equation multibody formulations,Journal of Computational and Nonlinear Dynamics 10 (1) (2014) 011012 (7 pages). doi:10.1115/1.4026492 .URL http://dx.doi.org/10.1115/1.4026492 [21] Y. Zhu, D. Dopico, A. Sandu, C. Sandu. Mbsvt. a library for the simulation and optimization of multibody systems [online] (2014) [citedJanuary 2015].[22] Y. Zhu, Sensitivity analysis and optimization of multibody systems, Ph.D. thesis, Virginia Tech (2014).[23] Y. Zhu, D. Dopico, C. Sandu, A. Sandu, Optimization of vehicle dynamics based on multibody models using adjoint sensitivity analysis,Journal of Computational and Nonlinear Dynamics submitted.[24] S. Kolathaya, A. D. Ames, Parameter to state stability of control lyapunov functions for hybrid system models of robots, Nonlinear Analysis:Hybrid Systems (2016) – doi:http://dx.doi.org/10.1016/j.nahs.2016.09.003 .[25] P. Ballard, The dynamics of discrete mechanical systems with perfect unilateral constraints, Archive for Rational Mechanics and Analysis154 (3) (2000) 199–274. doi:10.1007/s002050000105 .URL https://doi.org/10.1007/s002050000105 [26] A. S. Scha ff er, On the adjoint formulation of design sensitivity analysis of multibody dynamics, Ph.D. thesis, The University of Iowa (2005).er, On the adjoint formulation of design sensitivity analysis of multibody dynamics, Ph.D. thesis, The University of Iowa (2005).