text:
\section{Introduction}
Sensitivity analysis plays a key role in a wide range of computational engineering problems such as design optimization, optimal control, and implicit time integration methods, by providing derivative information for gradient based algorithms and methods. Sensitivity analysis quantifies the effect of small changes in the system parameters onto the outputs of interest ~\cite{Sandu2015dynamic}. Specifically, in the design of mechanical systems sensitivity analysis reveals the system parameters that affect the given performance criterion the most, thus providing directions for mechanical design improvements. Sensitivity analysis enables gradient-based optimization by providing the derivative of the cost function with respect to design variables.
In adaptive control systems sensitivity analysis allows to assess  stability by accounting for the effects of system disturbances and system parameters inaccuracies.

The most widely used techniques for sensitivity analysis are the direct and the adjoint methods.  These approaches are complementary, as the direct sensitivity provides information on how parametric uncertainties propagate through the system dynamics, while the adjoint method is suitable for inverse modeling, in the sense that it can be used to identify the origin of uncertainty in a given model output ~\cite{sandu2003direct}.

Numerical approximations of sensitivities are often calculated by finite difference methods where the deviation of the state trajectories are evaluated after system parameters disturbances or variations in the initial conditions are added to the system. Because of the simplicity of this method, which does not require any additional inputs other than the provided model, this technique is broadly used. However, the accuracy of the results is severely limited by the perturbation size, and by the roundoff and cancellation errors ~\cite{chang2016design}.

This paper develops a general and unify formulation for {\it direct} sensitivity analysis for hybrid dynamical systems.
In the context of this study, the term hybrid refers to a continuous system that encounters a finite number of events where some of the state variables jump to different values; the dynamics of a hybrid system is piecewise-smooth in time. Sensitivity analysis for hybrid systems has been an active topic of research ~\cite{Barton1998,Barton1999, barton2002modeling, tolsma2002hidden, rozenvasser1967general, saccon2014sensitivity, hiskens2006sensitivity, hiskens2000trajectory, taringoo2010geometry}.

We are especially interested in hybrid mechanical systems where sensitivity analysis involves the time, the position coordinates, the velocity coordinates, and the system parameters. We treat unconstrained systems modeled by ODEs, as well as constrained multibody systems modeled by differential algebraic equations (DAEs) and ODE penalty formulation ~\cite{Sandu2015dynamic,kurdila1993sufficient}. The sensitivity of the time of event, and the jumps in the sensitivities of the state variables at the time of event, are available in the literature. Here a new graphical proof of the jumps in sensitivities at the time of an event is employed, which helps to better understand the conditions for the jump in the sensitivities.  This paper provides a unified methodology for determining the system solutions, their sensitivities, and sensitivities of a cost function for different types of events.  The first type of event is caused by an external impulse (e.g., a contact) leading to a sudden change of velocities. The second type of event is caused by a sudden change of the equations of motion.  The third type of event is caused by  a sudden change in the kinematic constraints.

The paper is organized as follows. A review of the direct sensitivity approach for smooth ODE systems, along with the quadrature variable of the running cost function, is introduced in Section \ref{sec:Direct sensitivity analysis for smooth ODE systems}.  The  extension of this approach to hybrid ODE systems is presented in Section \ref{sec:direct-sensitivity-hybrid}. The sensitivity analysis for smooth constrained rigid multibody dynamic systems is reviewed in Section \ref{sec:multibody-smooth}.  Sensitivity analysis methodology for hybrid constrained systems is developed in Section \ref{sec:multibody-constrained}.
Sensitivity analysis methodology for constrained systems with a sudden change of the of equation of motions is developed in Section \ref{sec:Direct sensitivity analysis with transition functions}. The proposed  sensitivity analysis methodologies in Section \ref{sec:multibody-constrained} are applied to a five-bar mechanism  in Section \ref{sec:numerics}. Conclusions are drawn in Section \ref{sec:conclusions}.

\section{Direct sensitivity analysis for smooth ODE systems}
\label{sec:Direct sensitivity analysis for smooth ODE systems}

We start the discussion with a review of direct sensitivity analysis for dynamical systems governed by smooth ODEs.
\subsection{Smooth ODE system dynamics}

In this study we consider second order systems of ordinary differential equations of the form:
\begin{equation}
\label{eq:EOM}
{\mathsf{M}}\left(t,{q},{\rho}\right) \cdot \ddbq ={\mathsf{F}} \left(t,{q},\dbq,{\rho}\right),
\quad t_0 \leq t \leq  t_F, \quad {q}(t_0)={q}_0({\rho}),\quad\dbq(t_0)=\dbq_0({\rho}),
\end{equation}
or equivalently:
\begin{equation}
\label{eq:EOM-ODE}
\ddbq ={\mathsf{M}}^{-1} \left(t,{q},{\rho}\right)  \cdot {\mathsf{F}}\left(t,{q},\dbq,{\rho}\right)
=: \feom \left(t,{q},\dbq,{\rho}\right),
\end{equation}
that arise from the description of the dynamics of mechanical systems. In \eqref{eq:EOM-ODE} $t \in \mathds{R}$ is time, ${q} \in \mathds{R}^{n}$ is the generalized position vector and $\dbq \in \mathds{R}^{n}$ is the generalized velocity vector, $n$ is the dimension of generalized coordinates, and ${\rho} \in \mathds{R}^{p} $ is the vector of system parameters, where $p$ is the number of parameters. The dot notation ($\dot{\square}$ or $\ddot{\square}$) indicates the total (first or second order) derivative of a function or variable with respect to time. Subscripts indicate partial derivative with respect to a quantity, unless stated otherwise. The mass matrix $\mathsf{M}:\mathds{R} \times \mathds{R}^{n} \times \mathds{R}^{p} \rightarrow \mathds{R}^{n \times n}$ is assumed to be smooth with respect to all its arguments, invertible, and with an inverse $\mathsf{M}^{-1}$ that is also smooth with respect to all arguments. The right hand side function $\mathsf{F}:\mathds{R} \times \mathds{R}^{n} \times \mathds{R}^{n} \times \mathds{R}^{p} \rightarrow \mathds{R}^{n}$ represents external and internal generalized forces and is assumed to be smooth with respect to all its arguments.

The state trajectories are obtained by integrating the equations of motion \eqref{eq:EOM-ODE}, which depend on the system parameters ${\rho}$. Consequently, the state trajectories (the solutions of the equations of motion) depend implicitly on time and on the parameters, ${q}={q}(t, {\rho})$ and $\dbq=\dbq(t, {\rho})$. The state trajectories also depend implicitly on the initial conditions of \eqref{eq:EOM-ODE}. For clarity we denote the velocity state variables by ${v}=\dbq \in \mathds{R}^{n}$.

Sensitivity analysis computes derivatives of the solutions of \eqref{eq:EOM-ODE} with respect to  the system parameters:
\begin{equation}
\label{eq:sensitivity-of-solutions}
{ Q}(t,{\rho}) := D_\brho{q}(t) := \frac{d {q}}{d {\rho}}(t, {\rho}) \in \mathds{R}^{n \times p},
\quad {V}(t,{\rho}) := D_\brho{v}(t) = \frac{d {v}}{d {\rho}}(t, {\rho}) \equiv \frac{d \dbq}{d {\rho}}(t, {\rho})=\dot{{ Q}}(t,{\rho})  \in \mathds{R}^{n \times p}.
\end{equation}

The second order ODE~\eqref{eq:EOM-ODE} can transformed into a first order reduced system as follows. With the velocity state variables ${v}:=\dbq \in \mathds{R}^{n}$ the system \eqref{eq:EOM-ODE} can be written in the form:
\begin{equation}
\label{eq:EOM-ODE-first-order}
\begin{bmatrix}
\mathsf{I} & \mathsf{0} \\
\mathsf{0} & \mathsf{M}(t,{q},{\rho})
\end{bmatrix}
\begin{bmatrix}
 \dbq \\
 {\dot \bv}
\end{bmatrix}  =
\begin{bmatrix}
 {v} \\
 \mathsf{F}(t,{q},{v},{\rho})
\end{bmatrix}
\quad \Leftrightarrow \quad
\begin{bmatrix}
 \dbq \\
 {\dot \bv}
\end{bmatrix}
=
\begin{bmatrix}
 {v} \\
\feom \left(t,{q},{v},{\rho}\right)
\end{bmatrix},
\qquad
\begin{bmatrix}
 {q}(t_0) \\
 {v}(t_0)
\end{bmatrix}
=
\begin{bmatrix}
 {q}_0({\rho}) \\
 {v}_0({\rho})
\end{bmatrix}.
\end{equation}
\begin{definition}[Cost function]
\label{def:CostFunction}
Consider a smooth `trajectory cost function' ${g}:\mathds{R}^{1+3n+p} \rightarrow \mathds{R}$ and a smooth `terminal cost function' $w:\mathds{R} \times \mathds{R}^{1+2n+p} \rightarrow \mathds{R}$.
A general cost function is defined as the sum of the costs along the trajectory plus the cost at the terminal point of the solution:
\begin{eqnarray}
\label{eq:CostFunction}
&&\psi=  \int_{t_0}^{t_F} {\bg\fincost \  dt} + w \bigl(t_F,\bq(t_F,\brho),\bv(t_F,\brho),\brho \bigr).
\end{eqnarray}
\end{definition}
\begin{remark}[Accelerations in the cost function]
\label{rem:Accelerations}
Note that the trajectory cost function \eqref{eq:CostFunction} includes accelerations via ${\dot \bv}$. Accelerations are not independent variables and they can be resolved in terms of positions and velocities:
\[
g\left(\,t,\,{q},\,{v},\,{\dot \bv},\,{\rho}\,\right) = g\bigl(\,t,\,{q},\,{v},\,\feom(t,{q},{v},{\rho}),\,{\rho}\,\bigr) = \tilde{g}\bigl(\,t,\,{q},\,{v},\,{\rho}\,\bigr).
\]
We prefer to keep accelerations as an explicit argument in the cost function \eqref{eq:CostFunction} in order to give additional flexibility in practical applications. However, we will need to resolve the sensitivities of acceleration in terms of other sensitivities in subsequent calculations.
\end{remark}

To further simplify the notation we define the `quadrature' variable $z(t) \in \mathds{R}$ as follows:
\begin{equation}
\label{eqn:quadrature-variable}
\begin{split}
{z}(t,{\rho}) &:=  \int_{t_0}^{t} { \tilde{g}\bigl(t,{q}(t,{\rho}),{v}(t,{\rho}),{\rho}\bigr) \,  dt} \quad \Leftrightarrow \quad \\
{\dot \bz}(t,{\rho}) &=  \bg\fincost = \tilde{g}\bigl(t,{q}(t,{\rho}),{v}(t,{\rho}),{\rho}\bigr),  \quad t_0 \leq t \leq  t_F, \quad z(t_0,{\rho})=0.
\end{split}
\end{equation}
The cost function \eqref{eq:CostFunction} reads:
\begin{equation}
\label{eqn:quadrature-cost-function}
\psi ={z}(t_F) + w \bigl(t_F,\bq(t_F,\brho),\bv(t_F,\brho),\brho \bigr).
\end{equation}
\begin{definition}[The canonical ODE system]
\label{def:canonical-ode}
The canonical  system is obtained by combining the first order ODE dynamics \eqref{eq:EOM-ODE-first-order}  with equation \eqref{eqn:quadrature-variable} for the `quadrature' variable':
\begin{subequations}
\label{eq:canonical-ode}
\begin{eqnarray}
\label{eq:canonical-ode-system}
{x}(t) := \begin{bmatrix}  {{q}}(t) \\ {{v}}(t)  \\ {z}(t)   \end{bmatrix}  \in \mathds{R}^{2 n+1};
\qquad
\dot{{x}}
= \begin{bmatrix} {\dot \bv} \\ \feom\bigl(t,{q}(t,{\rho}),{v}(t,{\rho}),{\rho}\bigr)  \\ \tilde{g}{\left(\,t,\,\bq,\,\bv,\,\brho \, \right)} \end{bmatrix} = F(t,{x},{\rho}), \quad t_0 \le t \le t_F,
\quad
{x}(t_0,{\rho}) = \begin{bmatrix} {q}_0({\rho}) \\ {v}_0({\rho}) \\ 0 \end{bmatrix}.
\end{eqnarray}
The canonical cost function \eqref{eqn:quadrature-cost-function} is purely a terminal cost function:
\begin{equation}
\label{eqn:canonical-cost-function}
\psi ={z}(t_F,{\rho}) + w \bigl(t_F,\bq(t_F,\brho),\bv(t_F,\brho),\brho \bigr) = W\bigl( x(t_F,{\rho}), {\rho} \bigr).
\end{equation}
\end{subequations}
\end{definition}

\subsection{Direct sensitivity approach for smooth ODE systems}
\label{sec:direct-sensitivity-continuous}

\begin{definition}[The sensitivity analysis problem]
\label{def:Sensitivity}
Our goal in this work is to perform a sensitivity analysis of the cost function, i.e., to compute the total derivative of the cost function \eqref{eq:CostFunction} with respect to model parameters ${\rho}$:
\begin{equation}
\label{eq:gradient-cost-function}
D_\brho \psi = \frac{d\, \psi}{d\, {\rho}}  \in \mathds{R}^{1 \times p}.
\end{equation}
Note that the cost function \eqref{eq:CostFunction} depends on the system parameters ${\rho}$ directly (through the direct dependency of $g$ and $w$ on ${\rho}$) as well as indirectly (through the dependency of ${q}$ and ${v}$ on ${\rho}$). The sensitivity \eqref{eq:gradient-cost-function} needs to account for all the direct and the indirect dependencies.
\end{definition}

\begin{remark}[The direct sensitivity analysis approach]
In order to compute the sensitivities \eqref{eq:gradient-cost-function} we take a variational calculus approach \cite{Sandu_2014_sensitivity_ODE_multibody,Sandu2015dynamic,Sandu_2017_vehicle-optimization}.  Infinitesimal changes in parameters
\[
{\rho} \to {\rho} + \delta{\rho} \in \Re^p,
\]
lead to a total change in the cost function as follows:
\[
 \psi \to \psi + \delta \psi, \quad
\delta\psi = \sum_{i=1}^p \frac{d \psi}{d {\rho}_i} \cdot \delta{\rho}_i =
D_\brho \psi \cdot \delta{\rho}.
\]
The direct sensitivity analysis computes each element $(D_\brho \psi)_i = \partial \psi/\partial {\rho}_i$ of the derivative \eqref{eq:gradient-cost-function} by accounting for changes in the cost function that result from changing each individual parameter $\delta{\rho}_i$, $i=1 \dots p$.
\end{remark}

\subsection{Direct sensitivity analysis with respect to system parameters solved analytically}
\begin{definition}[The tangent linear model (TLM)] Consider the `position sensitivity' matrix $Q(t,{\rho})$ and the `velocity sensitivity' matrix  $V(t,{\rho})$ defined in \eqref{eq:sensitivity-of-solutions}:
\begin{subequations}
\label{eq:sensitivity-matrices}
\begin{eqnarray}
\label{eq:sensitivity-matrix-Q}
{ Q}_i(t,{\rho}) &:= \frac{d\, {q}(t,{\rho})}{d\, {\rho}_i} \in \mathds{R}^{n}, ~~ i=1,\dots,p; \quad
{ Q}(t,{\rho}) := \begin{bmatrix} { Q}_1(t,{\rho}) \cdots { Q}_p(t,{\rho}) \end{bmatrix} \in \mathds{R}^{n \times p}, \\
\label{eq:sensitivity-matrix-V}
{V}_i(t,{\rho}) &:= \frac{d\, {v}(t,{\rho})}{d\, {\rho}_i} \in \mathds{R}^{n}, ~~ i=1,\dots,p; \quad
{V}(t,{\rho}) := \begin{bmatrix} {V}_1(t,{\rho}) \cdots {V}_p(t,{\rho}) \end{bmatrix} \in \mathds{R}^{n \times p}.
\end{eqnarray}
\end{subequations}
These sensitivities evolve in time according to the {\it tangent linear model} (TLM) equations \cite{Sandu_2014_sensitivity_ODE_multibody,Sandu2015dynamic,Sandu_2017_vehicle-optimization}, obtained by differentiating the equations of motion \eqref{eq:EOM-ODE} with respect to the parameters:
\begin{subequations}
\label{eq:TLM}
\begin{equation}
\begin{cases}
\dot{{ Q}}_i &\displaystyle = \frac{d\dbq}{d {\rho}_i} = \frac{d{v}}{d {\rho}_i} = {V}_i, \\[8pt]
\dot{{V}}_i &\displaystyle =  \frac{d{\dot \bv}}{d {\rho}_i} = { \fun_\bq^{\scalebox{0.5}{\rm eom}}} \left(t,{q},{v},{\rho}\right) \cdot \frac{d{q}}{d {\rho}_i}
+ { \fun_\bv^{\scalebox{0.5}{\rm eom}}} \left(t,{q},{v},{\rho}\right) \cdot \frac{d{v}}{d {\rho}_i}
+ {\fun_{\brho_i}^{\scalebox{0.5}{\rm eom}}} \left(t,{q},{v},{\rho}\right) \\
&=  { \fun_\bq^{\scalebox{0.5}{\rm eom}}} \left(t,{q},{v},{\rho}\right) \cdot { Q}_i
+ { \fun_\bv^{\scalebox{0.5}{\rm eom}}} \left(t,{q},{v},{\rho}\right) \cdot {V}_i
+ {\fun_{\brho_i}^{\scalebox{0.5}{\rm eom}}} \left(t,{q},{v},{\rho}\right),
\end{cases}
\qquad i=1,\dots,p,
\quad t_0 \le t \le t_F,
\end{equation}
with the initial conditions
\begin{equation}
\label{eq:TLM-IC}
Q_i(t_0,{\rho}) = \frac{d{q}_0}{d {\rho}_i}, \quad
V_i(t_0,{\rho}) = \frac{d{v}_0}{d {\rho}_i}, \quad
i=1,\dots,p.
\end{equation}
\end{subequations}
The expressions ${ \fun_\bq^{\scalebox{0.5}{\rm eom}}}$, ${ \fun_\bv^{\scalebox{0.5}{\rm eom}}}$, and ${\fun_{\brho_i}^{\scalebox{0.5}{\rm eom}}}$ denote the partial derivatives of $\feom$ with respect to the subscripted variables.
\end{definition}

\begin{remark}
The partial derivatives $\partial \feom/\partial \zeta$ are obtained by differentiating  \eqref{eq:EOM-ODE} with respect to $\zeta \in \{ {q},{v},{\rho} \}$:
\begin{equation}
\label{eq:EOM-ODE}
\frac{\partial \feom}{\partial \zeta}
=\frac{\partial (\mathsf{M}^{-1}  \, \mathsf{F})}{\partial \zeta}
= -\mathsf{M}^{-1} \mathsf{M}_\zeta \mathsf{M}^{-1}\,\mathsf{F}
+ \mathsf{M}^{-1}\,\mathsf{F}_\zeta
= \mathsf{M}^{-1}\,\left( \mathsf{F}_\zeta - \mathsf{M}_\zeta\,\feom \right)
= \mathsf{M}^{-1}\,\left( \mathsf{F}_\zeta - \mathsf{M}_\zeta\,{\dot \bv} \right).
\end{equation}
\end{remark}

\begin{definition}[The quadrature sensitivity]
\label{The quadrature sensitivity}
Similarly, let the `quadrature sensitivity' vector $Z(t,{\rho})$ be the Jacobian of the `quadrature' variable $z(t,{\rho})$ ~\eqref{eqn:quadrature-variable} with respect to the parameters ${\rho}$:
\[
Z_i(t,{\rho}) := \frac{\partial z(t,{\rho})}{\partial {\rho}_i}, ~~ i=1,\dots,p; \quad
Z(t,{\rho}) := \nabla_{\rho} z(t,{\rho}) = \begin{bmatrix} Z_1(t,{\rho}) \cdots Z_p(t,{\rho}) \end{bmatrix} \in \mathds{R}^{1 \times p}.
\]
The time evolution equations of the quadrature sensitivities are given by the TLM obtained by differentiating \eqref{eqn:quadrature-variable} with respect to the parameters:
\begin{equation}
\label{eq:TLM-quadrature0}
\begin{split}
{\dot \bZ}_i  =  \frac{d\, g\left(\,t,\,\bq,\,\bv,\,{\dot \bv},\, \brho\,\right)}{d\,{\rho}_i} &= g_{q}\cdot { Q}_{i} +  g_{v}\cdot {V}_{i}  + g_{{\dot \bv}}\cdot\dot{{V}}_{i} + g_{{\rho}_i} \\
&= \big(g_{q} +  g_{\dot \bv}\,{ \fun_\bq^{\scalebox{0.5}{\rm eom}}}\big)\cdot { Q}_{i}
+  \big(g_{v} +  g_{\dot \bv}\,{ \fun_\bv^{\scalebox{0.5}{\rm eom}}}\big)\cdot {V}_{i}+ g_{{\rho}_i} +  g_{\dot \bv}\cdot{\fun_{\brho_i}^{\scalebox{0.5}{\rm eom}}}, \\
& \quad t_0 \le t \le t_F, \quad  { Z}_i(t_0,{\rho}) = 0, \quad i=1,\dots,p.
\end{split}
\end{equation}
\end{definition}
\begin{definition} [Canonical sensitivity ODE]
\label{def:Canonical sensitivity ODE}
The solutions given by Eq.~\eqref{eq:canonical-ode-system}, the TLM given by Eq.~\eqref{eq:TLM}, and the sensitivity quadrature equations \eqref{eq:TLM-quadrature} need to be solved together forward in time, leading to the canonical sensitivity ODE that computes the derivatives of the cost function with respect to the system parameters ${\rho}$ for smooth systems:
\begin{eqnarray}
\label{eq:canonical-ode-sensitivity}
\begin{bmatrix}
\dbq \\ {\dot \bv} \\ {\dot \bz} \\  \bigl[\dot{{ Q}}_i \bigr]_{i=1,\dots,p}\\
 \bigl[\dot{{V}}_i \bigr]_{i=1,\dots,p}\\  \bigl[{\dot \bZ}_i\bigr]_{i=1,\dots,p}
\end{bmatrix}
=
\begin{bmatrix}
{v} \\
\feom\\
 \tilde{{g}}\\
 \bigl[{V}_i \bigr]_{i=1,\dots,p}\\
 \bigl[{ \fun_\bq^{\scalebox{0.5}{\rm eom}}} { Q}_i +
{ \fun_\bv^{\scalebox{0.5}{\rm eom}}} {V}_i +
{\fun_{\brho_i}^{\scalebox{0.5}{\rm eom}}} \bigr]_{i=1,\dots,p} \\
 \bigl[\big(g_{q} +  g_{\dot \bv}\,{ \fun_\bq^{\scalebox{0.5}{\rm eom}}}\big)\cdot { Q}_{i}
+  \big(g_{v} +  g_{\dot \bv}\,{ \fun_\bv^{\scalebox{0.5}{\rm eom}}}\big)\cdot {V}_{i}+ g_{{\rho}_i} +  g_{\dot \bv}\cdot{\fun_{\brho_i}^{\scalebox{0.5}{\rm eom}}} \bigr]_{i=1,\dots,p}
\end{bmatrix},
\end{eqnarray}
where the state vector of the canonical sensitivity ODE is :
\begin{eqnarray}
\label{eq:canonical-sensitivity-state}
{X} = \left[   \, {q}^{\rm T}, \, {v}^{\rm T}, \, {z}, ~
{ Q}_1^{\rm T},  \dots , { Q}_p^{\rm T}, \,
{V}_1^{\rm T},  \dots , {V}_p^{\rm T}, \,
{ Z}_1, \dots, { Z}_p
\right] ^{\rm T}  \in \mathds{R}^{(n+1)(p+1)}.
\end{eqnarray}
\end{definition}

\begin{remark}  [The sensitivities of the cost function]
\label{rem:sensitivities-of-the-cost-function}
Once the quadrature sensitivities \eqref{eq:TLM-quadrature0} have been calculated, the sensitivities of the cost function with respect to each parameter are computed as follows:
\[
\begin{split}
\frac{d\,\psi}{d\, {\rho}_i} &= { Z}_i(t_F) + w_{q}\bigl(\,t_F, \,{q}(t_F,{\rho}),\, {v}(t_F,{\rho}),\, {\rho}\,\bigr)\cdot Q_i(t_F,{\rho})
+ w_{v}\bigl(\,t_F, \,{q}(t_F,{\rho}),\, {v}(t_F,{\rho}),\, {\rho}\,\bigr)\cdot V_i(t_F,{\rho}) \\
& \qquad + w_{{\rho}_i}\bigl(\,t_F, \,{q}(t_F,{\rho}),\, {v}(t_F,{\rho}),\, {\rho}\,\bigr),  \qquad i=1,\dots,p.
\end{split}
\]
\end{remark}

\subsection{Direct sensitivity analysis with respect to system parameters solved with the complex finite difference method}
\label{complex fdf}
An accurate numerical method for sensitivity analysis  of a smooth ODE system with respect to the system parameters ${\rho}$ is the complex finite difference method \cite{Sandu_2014_sensitivity_ODE_multibody,Sandu2015dynamic,Sandu_2017_vehicle-optimization}. Add a small {\it complex} perturbation to one parameter:
\[
\tilde{{\rho}}_j = \begin{cases}
{\rho}_j  & \textnormal{for}~~j \ne \ell, \\
{\rho}_\ell + \mathfrak{i}\,\Delta{\rho}& \textnormal{for}~~j = \ell,
\end{cases}
\quad j=1,\dots,p,
\]
and solve the canonical ODE system \eqref{eq:canonical-ode-system} for this perturbed values of the parameters to obtain:
\[
{q}(t,\tilde{{\rho}}), \quad {v}(t,\tilde{{\rho}}), \quad {z}(t,\tilde{{\rho}}),
\quad \psi(\tilde{{\rho}})={z}(t_F,\tilde{{\rho}}) + w \bigl(\, {q}(t_F,\tilde{{\rho}}),\, {v}(t_F,\tilde{{\rho}}),\,\tilde{{\rho}}\, \bigr).
\]
The sensitivities are approximated numerically by the imaginary parts of the state variables:
\begin{equation}
\label{eq:ODE_sensitivity_complex0}
{ Q}_\ell(t,{\rho})  \approx - \frac{\operatorname{imag}\bigl({{q}(t,\tilde{{\rho}})}\bigr)}{\Vert{\Delta {\rho}}\Vert}, \quad
{V}_\ell(t,{\rho})  \approx - \frac{\operatorname{imag} \bigl({v}(t,\tilde{{\rho}}) \bigr)}{\Vert{\Delta {\rho}}\Vert}, \quad
\frac{d\, \psi}{d\, {\rho}_\ell}  \approx - \frac{\operatorname{imag} \bigl(\psi(\tilde{{\rho}}) \bigr)}{\Vert{\Delta {\rho}}\Vert}.
\end{equation}

We next discuss an approach to sensitivity analysis that accounts for discontinuities in the state variables.

\section{Direct sensitivity analysis for hybrid ODE systems}
\label{sec:direct-sensitivity-hybrid}
\subsection{Hybrid ODE systems}
\begin{definition}[Hybrid dynamics]
\label{def:time-of-invent}
A hybrid mechanical system is  a piecewise-in-time continuous dynamic ODE described by \eqref{eq:EOM-ODE} that exhibits discontinuous dynamic behavior in the generalized velocity state vector at a finite number of time moments (no zeno phenomena \cite{Pace2017}).  Each such moment is a `time of event' ${t_{\rm eve}}$ and corresponds to a triggering event described by the equation:
\begin{equation}
\label{eq:event_function}
r \big( \bq|_{\timp} \big) = 0,
\end{equation}
where $r : \mathds{R}^{n} \to \mathds{R}$ is a smooth `event function'.
\end{definition}

\begin{remark}
In the context of this paper, we assume that there are no grazing phenomena where the system trajectory would make tangential contact with an the event triggering hypersurface \cite{Pace2017}.
\end{remark}

\begin{definition}[Characterization of an event]
\label{def:characterize-invent}
For hybrid systems variables can change values during the event. For this reason we distinguish between the value of a variable right before the event $\boldsymbol{x}|_{\timp}^{-}$, and its value right after the event $\boldsymbol{x}|_{\timp}^{-}$:
\[
\boldsymbol{x}|_{\timp}^{-} := \lim_{\varepsilon>0,~\varepsilon\to 0}\,\boldsymbol{x}({t_{\rm eve}}-\varepsilon), \quad
\boldsymbol{x}|_{\timp}^{+} := \lim_{\varepsilon>0,~\varepsilon\to 0}\,\boldsymbol{x}({t_{\rm eve}}+\varepsilon).
\]
The limits exist since the evolution of the system is smooth in time both before and after the event.

We consider an event happening at time ${t_{\rm eve}}$ that applies a finite energy impulse force to the system. Such an impulse force does not change the generalized position state variables, and therefore:
\begin{equation}
\label{eq:jump-in-position}
{\bq \timpplus} = \bq \timpminus = {q}|_{\timp}.
\end{equation}
However, the finite energy event can abruptly change the generalized velocity state vector $\dbq$ from its value $\bv \timpminus$ right before the event to a new value ${\bv \timpplus} $ right after the event.  The `jump function' at the time of event ${t_{\rm eve}}$ characterizes the change in the generalized velocity during the event:
\begin{equation}
\label{eq:jump-in-velocity}
{\bv \timpplus}=h \Big({{t_{\rm eve}}},\, {q}|_{\timp},\,\bv \timpminus,\,{{\rho}} \Big)  \quad \Leftrightarrow \quad  \dbq |_{\timp}^{+}=h \Big({{t_{\rm eve}}},\, {q}|_{\timp},\,\dbq  |_{\timp}^{-},\,{{\rho}} \Big) .
\end{equation}
\end{definition}

\begin{remark}[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_\ell : \mathds{R}^n \to \mathds{R}$, $\ell=1,\dots,e$. The detection of the next event, which can be one of the possible $e$ options, is described by:
\[
r_1 \big(  \bq|_{\timp} \big) \cdot r_2 \big( \bq|_{\timp} \big)~ \dots ~ r_e \big(  \bq|_{\timp}\big)= 0,
\]
and if event $\ell$ takes place then $r_\ell = 0$ and the corresponding jump is:
\[
{\bv \timpplus} = h_\ell \Big(\bq|_{\timp},\,\bv \timpminus\Big).
\]
\end{remark}

\begin{remark}[Numerical implementation of invents]
Numerical solutions of hybrid systems use an event detection mechanism. The event function \eqref{eq:event_function} is implemented in the numerical time solver such that the integrator is stopped at the solution of \eqref{eq:event_function}. The jump function \eqref{eq:jump-in-velocity} is implemented as a callback function that is executed after the event is detected. The numerical integration resumes with new initial conditions after the jump.
\end{remark}

\begin{definition}[Twin perturbed systems]
\label{def:perturbed-twins}
Consider two versions of the system \eqref{eq:EOM-ODE} with identical dynamics and initial conditions, but with different parameters values ${\rho}_1$ and ${\rho}_2$, respectively. Without loss of generality in this proof we consider the scalar parameter case $p=1$; the general equation \eqref{eq:jump-position-sensitivity} can be proven element by element by considering sensitivities with respect to individual parameters. The two parameters are infinitesimally small perturbations $\delta{\rho}$ of the reference parameter value ${\rho}$:
\[
{\rho}_1 = {\rho} - \frac{\delta{\rho}}{2}; \qquad  {\rho}_2 = {\rho} + \frac{\delta{\rho}}{2}.
\]
We denote by ${q}_1(t) = {q}(t,{\rho}_1)$, ${v}_1(t) = {v}(t,{\rho}_1)$, and ${z}_1(t) = {z}(t,{\rho}_1)$ the position and velocity states, and the quadrature variable of the first system, respectively. We denote by ${q}_2(t) = {q}(t,{\rho}_2)$, ${v}_2(t) = {v}(t,{\rho}_2)$ , and ${z}_2(t) = {z}(t,{\rho}_2)$ the position and velocity states, and the quadrature variable of the second system, respectively.

Assume that the sign of the perturbation $\delta{\rho}$ is such that $\tau_{2} > {t_{\rm eve}} > \tau_{1}$, and denote $\delta\tau = \tau_{2} - \tau_{1}$. Since $\delta{\rho}$ is infinitesimally small, so is $\delta\tau$. The trajectories of the positions ${q}(t)$, $ {q}_1(t)$, and ${q}_2(t)$, as well as the trajectories of the velocities ${v}(t)$, ${v}_1(t)$, and $ {v}_2(t)$, are schematically illustrated in Fig.~\ref{fig:im1} and Fig.~\ref{fig:im2}, respectively. As shown in Fig.~\ref{fig:im1} the first system meets the event described by the function \eqref{eq:event_function} at the time of event ${t_{\rm eve}}({\rho}_1) = \tau_{1}$, when its position state is ${\bq_1}|_{\tau_1}$. The second system meets the event at time ${t_{\rm eve}}({\rho}_2) =\tau_{2}$, when its position state is ${\bq_2}|_{\tau_2}$. Note that in the limit of vanishing $\delta{\rho}$ we have:
\begin{equation}
\label{eqn:twin-system-limits}
\delta{\rho} \to 0 \quad \Rightarrow \quad
{\bv_1}|_{\tau_1}^{-} \to \bv \timpminus, ~~
{\bv_1}|_{\tau_2} \to {\bv \timpplus}, ~~
{\bv_1}|_{\tau_1}^{+} \to {\bv \timpplus}; \quad
{\bv_2}|_{\tau_1}\to \bv \timpminus, ~~
{\bv_2}|_{\tau_2}^{-} \to \bv \timpminus, ~~
{\bv_2}|_{\tau_2}^{+} \to {\bv \timpplus}.
\end{equation}
\end{definition}

Let ${ Q}|_{\timp}^{+},{ Q}|_{\timp}^{-} \in \mathds{R}^{n{\times}p}$  be the sensitivities of the generalized positions \eqref{eq:sensitivity-matrix-Q} right before and right after the event, respectively. Let ${V}|_{\timp}^{+},{V}|_{\timp}^{-} \in \mathds{R}^{n{\times}p}$  be the sensitivities of the generalized velocities \eqref{eq:sensitivity-matrix-V} right before and right after the event, respectively. Our methodology to find these sensitivities is to first evaluate the states ${q}_1(t),{v}_1(t)$ and ${q}_2(t),{v}_2(t)$ of each of the twin perturbed systems at both $\tau_1$ and $\tau_2$. The sensitivities are obtained from their definition by taking differences of states of the two systems, dividing them by the perturbation in parameters, and taking the limits, for example:
\[
\begin{split}
Q_i|_{\timp}^{-} = \lim_{\delta{\rho}_i \to 0}\, \frac{{\bq_2}|_{\tau_1} - {\bq_1}|_{\tau_1}}{ \delta{\rho}_i},\quad
V_i|_{\timp}^{+} = \lim_{\delta{\rho}_i \to 0}\,  \frac{{\bv_2}|_{\tau_2}^{+} - {\bv_1}|_{\tau_2}}{ \delta{\rho}_i},
\quad i=1,\dots,p.
\end{split}
\]
\begin{figure} [H]
	\begin{center}
	\includegraphics[width=80mm]{Figures/position2.pdf}
	\end{center}
	\caption{Schematic visualization of the jump in the sensitivity of the position.}
	\label{fig:im1}
\end{figure}
\begin{figure} [H]
	\begin{center}
	\includegraphics[width=80mm]{Figures/velocity2.pdf}
	\end{center}
	\caption{Schematic visualization of the jump in the sensitivity of the velocity. }
	\label{fig:im2}
\end{figure}

\subsection{The sensitivity of the time of event with respect to the system parameters }

\begin{theorem}[Sensitivity of the time of event \cite{Donze07systematicsimulation, backer1966jump, hiskens2000trajectory,Barton1999,Barton1998,rozenvasser1967general}]
\label{thm:sensitivity-time-of-invent}
Let $r(\cdot) \in \mathds{R}$ be the scalar event function defined by Eq.~\eqref{eq:event_function}, and $dr/d{q} \in \mathds{R}^{1 \times n }$ be its Jacobian. The sensitivity of the time of event with respect to the system parameters is:
\begin{eqnarray}
\label{eq:sensitivity-time-of-invent}
\frac{d\timp}{d\brho}=
- \, \dfrac{\displaystyle \frac{d r}{d{q}}\left(\bq|_{\timp}\right)  \cdot { Q}|_{\timp}^{-}}{\displaystyle \frac{d r}{d{q}}\left(\bq|_{\timp}\right)  \cdot  \bv \timpminus }
\in \mathds{R}^{1 \times p}.
\end{eqnarray}
\end{theorem}
\begin{pf}
The time at which the event function becomes zero is indirectly dependent on the system parameters ${\rho}$. We evaluate the derivative of equation \eqref{eq:event_function} with respect to the system parameters:
\begin{equation}
\label{eq:event_function1}
0 = r \ \bigl( {q}({t_{\rm eve}},{\rho}) \bigr) \quad \Rightarrow \quad
0=\frac{dr}{d\brho}=\frac{d r}{d{q}}\, \left( \frac{d{q}}{d{\rho}} + \dbq \, \frac{d\timp}{d\brho} \right).
\end{equation}
Rearrange the terms in Eq.~\eqref{eq:event_function1} to obtain Eq.~\eqref{eq:sensitivity-time-of-invent}.
\qed
\end{pf}

\subsection{The jump in the sensitivity of the position state vector due to the event}
This section provides the jumps in the sensitivities of the position state vector ${q}(t)$ at the time of event  \cite{Donze07systematicsimulation, backer1966jump, hiskens2000trajectory,Barton1999,Barton1998,rozenvasser1967general}. Due to the nonzero inertia the position state variable is continuous in time \eqref{eq:jump-in-position}. However, its sensitivity can be discontinuous at the time of event, as established next.

\begin{theorem}[Jump in position sensitivity \cite{Donze07systematicsimulation, backer1966jump, hiskens2000trajectory, barton2002modeling}]
\label{thm:jump-in-position-sensitivity}
Let ${\bv \timpplus},\bv \timpminus \in \mathds{R}^{n}$ be the generalized velocity state vectors  after and before the invent, respectively; the corresponding velocity jump function was introduced in Eq.~\eqref{eq:jump-in-velocity}. Let ${ Q}|_{\timp}^{+}$ and ${ Q}|_{\timp}^{-} \in \mathds{R}^{n{\times}p}$  be the sensitivities of the generalized position state vectors after and before the invent, respectively. The jump equation of the sensitivities of the generalized position state vector is:
\begin{eqnarray}
\label{eq:jump-position-sensitivity}
{ Q}|_{\timp}^{+} = { Q}|_{\timp}^{-} - \bigg( {\bv \timpplus} - \bv \timpminus \bigg) \cdot \frac{d\timp}{d\brho}.
\end{eqnarray}
\end{theorem}
\begin{pf}
Consider the twin perturbed systems from Definition \ref{def:perturbed-twins}. The evolution of positions is illustrated in Fig.~\ref{fig:im1}, where the two different dashed line trajectories represent the position variables of the two perturbed systems.
The jump in the velocity state variables occurs at time $\tau_1$ only for the first system. The position variables at time $\tau_2$ for both systems are:
\begin{eqnarray}
\label{eq:jump-in-position-sensitivity-1}
\begin{split}
{\bq_1}|_{\tau_2} &= {\bq_1}|_{\tau_1} + h\big( {\bv_1}|_{\tau_1}^{-} \big) \, \delta{\tau}, \\
{\bq_2}|_{\tau_2}  &= {\bq_2}|_{\tau_1} + {\bv_2}|_{\tau_1}\, \delta{\tau}.
\end{split}
\end{eqnarray}
Subtract the two equations and scale by the perturbation in the parameters:
\begin{equation}
\label{eq:jump-in-position-sensitivity-2}
\frac{ {\bq_2}|_{\tau_2} - {\bq_1}|_{\tau_2} }{\delta{\rho}}=
-\left( {\bv_1}|_{\tau_1}^{+} - {\bv_2}|_{\tau_1} \right)\,\frac{\delta{\tau}}{\delta{\rho}}  +
\frac{{\bq_2}|_{\tau_1} - {\bq_1}|_{\tau_1}}{\delta{\rho}}.
\end{equation}
Using \eqref{eqn:twin-system-limits} and taking the limit $\delta{\rho} \to 0$ in \eqref{eq:jump-in-position-sensitivity-2} we obtain \eqref{eq:jump-position-sensitivity}. The trajectory state differences are illustrated by the vertical lines in Fig.~\ref{fig:im1}.
\end{pf}

\subsection{The jump in the sensitivity of the velocity state vector due to the event}
This section provides the jumps in the sensitivities of the velocity state vector ${v}(t)$ at the time of event  \cite{hiskens2000trajectory,Barton1999,Barton1998,rozenvasser1967general} corresponding to the jump function \eqref{eq:jump-in-velocity}.

\begin{theorem}[Jump in velocity sensitivity.]  \cite{hiskens2000trajectory,Barton1999,Barton1998,rozenvasser1967general}
\label{thm:jump-in-velocity-sensitivity}
Let ${\bV \afterimp},{\bV \beforeimp}  \in \mathds{R}^{n{\times}p}$ be the sensitivities of the generalized position state vectors after and before the invent, respectively. Let ${\bv \timpplus}$ and $\bv \timpminus \in \mathds{R}^{n}$ be the velocity state vectors after and before the invent affected by the jump function Eq.~\eqref{eq:jump-in-velocity} , respectively.
Let $\ddbq|_{\timp}^{+}$ and $\ddbq|_{\timp}^{-} \in \mathds{R}^{n}$ be the generalized acceleration state vectors after and before the invent, respectively.  The jump equation of the sensitivities of the generalized velocity state vector is:
\begin{eqnarray}
\label{eq:jump-in-velocity-sensitivity}
{\bV \afterimp} &=& h_{q}|_{\timp}^{-} \cdot {\bQ \beforeimp} +h_{v}|_{\timp}^{-} \cdot  {\bV \beforeimp}
+\left( h_{q}|_{\timp}^{-}  \cdot \bv \timpminus -\ddbq |_{\timp}^{+} + h_{v}|_{\timp}^{-}  \cdot \ddbq |_{\timp}^{-}  {+ h_t}|_{\timp}^{-} \right) \cdot \frac{d\timp}{d\brho} {+ h_{\rho}|_{\timp}^{-}},
\end{eqnarray}
where the Jacobians of the jump function are:
\[
\begin{split}
h_t|_{\timp}^{-} &:= \frac{\partial h}{\partial t}\big({t_{\rm eve}}, \bq|_{\timp}, \bv \timpminus ,{\rho}\big) \in \Re^{f}, \qquad
h_{q}|_{\timp}^{-} := \frac{\partial h}{\partial{q}}\big({t_{\rm eve}}, \bq|_{\timp}, \bv \timpminus ,{\rho}\big)\in \Re^{f \times n}, \\
h_{v}|_{\timp}^{-} &:= \frac{\partial h}{\partial{v}}\big({t_{\rm eve}},\bq|_{\timp}, \bv \timpminus,{\rho} \big)\in \Re^{f \times f}, \qquad
h_{\rho}|_{\timp}^{-} := \frac{\partial h}{\partial{\rho}}\big({t_{\rm eve}},\bq|_{\timp}, \bv \timpminus,{\rho} \big)\in \Re^{f \times p}.
\end{split}
\]
\end{theorem}

\begin{pf}
We consider again the twin perturbed systems from Definition \ref{def:perturbed-twins}. The jumps in velocities are illustrated in Fig.~\ref{fig:im2}. The velocities for each system are determined as follows:
\begin{equation}
\label{eq:jump sensitivity velocity}
\begin{split}
{\bv_1}|_{\tau_2} &={\bv_1}|_{\tau_1}^{+}  + \feom \Big(\,\tau_1, \,{\bq_1}|_{\tau_1}, \,{\bv_1}|_{\tau_1}^{+}, \,{\rho}_1\, \Big) \, \delta{\tau},\\
 &=h\big( \, \tau_1,  \, {\bq_1}|_{\tau_1},  \,{\bv_1}|_{\tau_1}^{-}, \,{\rho}_1\, \big) + \feom \Big(\,\tau_1, \,{\bq_1}|_{\tau_1}, \,h\Big(\, \tau_1, \,{\bq_1}|_{\tau_1}, \,{\bv_1}|_{\tau_1}^{-}, \,{\rho}_1\, \Big), \,{\rho}_1\, \Big) \, \delta{\tau},\\
{\bv_2}|_{\tau_2}^{+} &= h \ \left(\, \tau_2,  \,{\bq_2}|_{\tau_2}, \,{\bv_2}|_{\tau_2}^{-},  \,{\rho}_2\,  \,\right) \\
 &= h \ \Big(\, \tau_2,  \,{\bq_2}|_{\tau_1}+{\bv_2}|_{\tau_1} \ \delta{\tau},\; {\bv_2}|_{\tau_1}  + \feom \Big(\tau_1, \,{\bq_2}|_{\tau_1}, \,{\bv_2}|_{\tau_1}, \,{\rho}_2 \Big) \, \delta{\tau}, \,{\rho}_2\,  \, \Big) \\
&\approx  h\big(\, \tau_2,  \,{\bq_2}|_{\tau_1},{\bv_2}|_{\tau_1}, \,{\rho}_2\,  \,  \big) \; + \;
\frac{d h}{d {q}}\big({\bq_2}|_{\tau_1},{\bv_2}|_{\tau_1} \big) \cdot {\bv_2}|_{\tau_1} \ \delta{\tau} \\
&+\frac{d h}{d {v}}\big({\bq_2}|_{\tau_1},{\bv_2}|_{\tau_1} \big)  \cdot
 \feom \Big(\,\tau_1, \,{\bq_2}|_{\tau_1}, \,{\bv_2}|_{\tau_1}, \,{\rho}_2\, \Big) \, \delta{\tau},
\end{split}
\end{equation}
where ${\feom}$ is the instantaneous acceleration of the system from Eq.~\eqref{eq:EOM-ODE}. The last relation represents a linearization (first order Taylor expansion) that is infinitely accurate since $\delta\tau$ is infinitesimally small.
The scaled difference between the velocity state vectors at the time of event is :
\begin{eqnarray*}
\frac{{\bv_2}|_{\tau_2}^{+} - {\bv_1}|_{\tau_2}}{ \delta{\rho}} &\approx&
\frac{h\big(\,\tau_2, \,{\bq_2}|_{\tau_1},{\bv_2}|_{\tau_1}, \,{\rho}_2\,\big) - h\big(\,\tau_1, \,{\bq_1}|_{\tau_1},{\bv_1}|_{\tau_1}^{-}, \,{\rho}_1\,\big) }{\delta{\rho}} \\
&& -\feom \Big(\, \tau_1,\, {\bq_1}|_{\tau_1},\, h \big({\bq_1}|_{\tau_1}, {\bv_1}|_{\tau_1}^{-}\big),\, {\rho}_1\, \Big)\, \frac{\delta{\tau}}{\delta{\rho}} \\
&& +\frac{d h}{d {q}}\big({\bq_2}|_{\tau_1},{\bv_2}|_{\tau_1} \big) \cdot {\bv_2}|_{\tau_1}
\cdot \frac{\delta{\tau}}{\delta{\rho}}
+\frac{d h}{d {v}}\big({\bq_2}|_{\tau_1},{\bv_2}|_{\tau_1}\big) \cdot
 \feom \Big(\,\tau_1, \,{\bq_2}|_{\tau_1}, \,{\bv_2}|_{\tau_1}, \,{\rho}_2\, \Big) \cdot \frac{\delta{\tau}}{\delta{\rho}} \\
&\approx&
\frac{d h}{d t}\big({\bq_1}|_{\tau_1},{\bv_1}|_{\tau_1}^{-} \big) \cdot \frac{\tau_2 - \tau_1}{ \delta{\rho}} +
\frac{d h}{d {q}}\big({\bq_1}|_{\tau_1},{\bv_1}|_{\tau_1}^{-} \big) \cdot \frac{{\bq_2}|_{\tau_1} - {\bq_1}|_{\tau_1}}{ \delta{\rho}} +
\frac{d h}{d {v}}\big({\bq_1}|_{\tau_1},{\bv_1}|_{\tau_1}^{-} \big) \cdot \frac{{\bv_2}|_{\tau_1} - {\bv_1}|_{\tau_1}^{-}}{ \delta{\rho}} \\
&& + \frac{d h}{d {\rho}}\big({\bq_1}|_{\tau_1},{\bv_1}|_{\tau_1}^{-} \big)
-\feom \Big(\tau_1,\,{\bq_1}|_{\tau_1},~ h \big({\bq_1}|_{\tau_1},{\bv_1}|_{\tau_1}^{-}),\,{\rho}_1 \Big)\, \frac{\delta{\tau}}{\delta{\rho}} \\
&& +\frac{d h}{d {q}}\big({\bq_2}|_{\tau_1},{\bv_2}|_{\tau_1} \big) \cdot {\bv_2}|_{\tau_1}
\cdot \frac{\delta{\tau}}{\delta{\rho}}
 +\frac{d h}{d {v}}\big({\bq_2}|_{\tau_1},{\bv_2}|_{\tau_1} \big) \cdot
 \feom \Big(\tau_1, \,{\bq_2}|_{\tau_1}, \,{\bv_2}|_{\tau_1}, \,{\rho}_2 \Big) \cdot  \frac{\delta{\tau}}{\delta{\rho}}.
\end{eqnarray*}
Taking the limit $\delta{\rho} \to 0$ and using Eq. \eqref{eqn:twin-system-limits} yields:
\[
\begin{split}
&\frac{d h}{d {q}}\big({\bq_1}|_{\tau_1},{\bv_1}|_{\tau_1}^{-} \big) \to
h_{q}|_{\timp}^{-}, \quad
\frac{d h}{d {q}}\big({\bq_2}|_{\tau_1},{\bv_2}|_{\tau_1} \big) \to
h_{q}|_{\timp}^{-}, \\
&\frac{d h}{d {v}}\big({\bq_1}|_{\tau_1},{\bv_1}|_{\tau_1}^{-} \big) \to
h_{v}|_{\timp}^{-}, \quad
\frac{d h}{d {v}}\big({\bq_2}|_{\tau_1},{\bv_2}|_{\tau_1} \big) \to
h_{v}|_{\timp}^{-},\\
&\frac{d h}{d t}\big({\bq_1}|_{\tau_1},{\bv_1}|_{\tau_1}^{-} \big) \to
h_t|_{\timp}^{-}, \quad
\frac{d h}{d {\rho}}\big({\bq_1}|_{\tau_1},{\bv_1}|_{\tau_1}^{-} \big) \to
h_{\rho}|_{\timp}^{-},\\
&\feom \Big(\tau_1, \,{\bq_2}|_{\tau_1}, \,{\bv_2}|_{\tau_1}, \,{\rho}_2 \Big)
\to \feom \Big({t_{\rm eve}}, \,\bq \timpminus, \,\bv \timpminus, \,{\rho} \Big)
= \ddbq |_{\timp}^{-},\\
&\feom \Big(\tau_1, \,{\bq_1}|_{\tau_1},{\bv_1}|_{\tau_1}^{+}, \,{\rho}_1 \Big)
 \to \feom \Big({t_{\rm eve}}, \,{\bq \timpplus}, \,{\bv \timpplus}, \,{\rho} \Big)
= \ddbq |_{\timp}^{+}.
\end{split}
\]
which leads to Eq.~\eqref{eq:jump-in-velocity-sensitivity}.

For simplicity we denote the derivatives of the jump function  with respect to $\zeta \in \{ t, {q},{v},{\rho} \}$ by:
\[
\begin{split}
&\frac{d h}{d \zeta}\big(\tau_1, \,{\bq_1}|_{\tau_1}, \,{\bv_1}|_{\tau_1}^{-}, \,{\rho}_1\, \big) =
\frac{d h}{d \zeta}\big({\bq_1}|_{\tau_1},{\bv_1}|_{\tau_1}^{-} \big) \\
&\frac{d h}{d {q}}\big(\tau_1, \,{\bq_2}|_{\tau_1}, \,{\bv_2}|_{\tau_1}, \,{\rho}_2\, \big) =
\frac{d h}{d \zeta}\big({\bq_2}|_{\tau_1},{\bv_2}|_{\tau_1} \big) \\
\end{split}
\]
\qed
\end{pf}

\begin{theorem}[Events that only change the acceleration]
\label{thm:ODE-change-in-acceleration}
We now consider an event where the system undergoes a sudden change of the equation of motions \eqref{eq:EOM} at ${t_{\rm eve}}$:
\[
\ddbq |_{\timp}^{-} = \feom^- \big({t_{\rm eve}},\bq|_{\timp},\bv|_{\timp}, {\rho}\big) =:  \feom^-|_{\timp}
\quad \stackrel{\rm event}{\longrightarrow} \quad
\ddbq |_{\timp}^{+} = \feom^+ \big({t_{\rm eve}},\bq|_{\timp},\bv|_{\timp}, {\rho}\big)=:  \feom^+|_{\timp}.
\]
There is no abrupt jump in the system velocity, ${\bv \timpplus}=\bv \timpminus$, and therefore the jump function \eqref{eq:jump-in-velocity} is identity. Let ${\bV \afterimp},{\bV \beforeimp}  \in \mathds{R}^{n{\times}p}$ be the sensitivities of the generalized position state vectors right after and before the event, respectively.  Let $\ddbq|_{\timp}^{+}$ and $\ddbq|_{\timp}^{-} \in \mathds{R}^{n}$ be the generalized acceleration state vectors right after and before the event, respectively.
The jump equation of the sensitivities of the generalized velocity state vector is:
\begin{eqnarray}
\label{eq:jump-in-velocity-sensitivity-changeacc}
{\bV \afterimp} = {\bV \beforeimp} -\left( \ddbq |_{\timp}^{+} - \ddbq |_{\timp}^{-} \right) \cdot \frac{d\timp}{d\brho}
= {\bV \beforeimp} -\left( \feom^+|_{\timp} - \feom^-|_{\timp} \right) \cdot \frac{d\timp}{d\brho}.
\end{eqnarray}
\end{theorem}

\begin{pf}
For the type of events under consideration we have that:
\[
\frac{d h}{d {q}}= \mathsf{0}, \quad \frac{d h}{d {v}}=\mathsf{I}.
\]
Using this in Eq. \eqref{eq:jump-in-velocity-sensitivity} leads to Eq.~ \eqref{eq:jump-in-velocity-sensitivity-changeacc}.

\qed
\end{pf}

\subsection{The jump in the sensitivity of the cost functional due to the event}
We now consider the sensitivity of the quadrature variable ${z}(t)$. Due to the integral form of Eq. \eqref{eqn:quadrature-variable} defining ${z}$,  the quadrature variable is continuous in time:
\[
{z \timpplus}= z \timpminus = z|_{\timp}.
\]
However, its sensitivity can be discontinuous at the event time, as established next.

\begin{theorem}[Jump in quadrature sensitivity.]
\label{thm:jump-in-quadrature-sensitivity}
Let ${\bZ \afterimp}$ and ${\bZ \beforeimp}$, with  ${ Z} \in \mathds{R}^{p}$, be the sensitivities of the quadrature variable ${z}(t)$ (Definition \ref{The quadrature sensitivity}) right after and right before the event, respectively.
Let
\[
\bg \timpplus := \tilde{g} \big(\timp,\bq|_{\timp},{\bv \timpplus}, \brho\big),
\quad
\bg \timpminus := \tilde{g} \big(\timp,\bq|_{\timp},\bv \timpminus, \brho\big),
\]
be the running cost function evaluated right after and right before the event, respectively.
The sensitivity of the cost functional changes during the event as follows:
\begin{eqnarray}
\label{eq:jump-in-quadrature-sensitivity}
{\bZ \afterimp}= {\bZ \beforeimp} - \Big( \bg \timpplus - \bg \timpminus \Big) \cdot \frac{d\timp}{d\brho}.
\end{eqnarray}
\end{theorem}
\begin{pf}
Consider again the twin perturbed systems from Definition \ref{def:perturbed-twins}, and evaluate the associated quadrature variables \eqref{eqn:quadrature-variable} at the event:
\begin{equation}
\begin{split}
{\bz_1}|_{\tau_2} &=
{\bz_1}|_{\tau_1}
+ \int_{\tau_1}^{\tau_2}  \tilde{g}\bigl(t,{q}_1(t),{v}_1(t),{\rho}_1\bigr) \, dt
= {\bz_1}|_{\tau_1}
+ \tilde{g}\bigl(\,\tau_1,\,{q}_1|_{\tau_1},\,{v}_1|_{\tau_1}^+,\,{\rho}_1\,\bigr) \, \delta{\tau}, \\
{\bz_2}|_{\tau_2}  &=
{\bz_2}|_{\tau_1}
+ \int_{\tau_1}^{\tau_2}  \tilde{g}\bigl(t,{q}_2(t),{v}_2(t),{\rho}_2\bigr)\, dt
= {\bz_2}|_{\tau_1}
+ \tilde{g}\bigl(\,\tau_2,\,{q}_2|_{\tau_2},\,{v}_2|_{\tau_2}^-,\,{\rho}_2\,\bigr) \, \delta{\tau}.
\end{split}
\end{equation}
Subtract the two equations and scale by the parameter perturbation to obtain:
\[
\begin{split}
\frac{ {\bz_2}|_{\tau_2} - {\bz_1}|_{\tau_2} }{\delta{\rho}}  &=
\frac{ {\bz_2}|_{\tau_1}- {\bz_1}|_{\tau_1} }{\delta{\rho}}
 + \left(
\tilde{g}\bigl(\,\tau_2,\,{q}_2|_{\tau_2},\,{v}_2|_{\tau_2}^-,\,{\rho}_2\,\bigr)
-
\tilde{g}\bigl(\,\tau_1,\,{q}_1|_{\tau_1},\,{v}_1|_{\tau_1}^+,\,{\rho}_1\,\bigr)
\right)\, \frac{\delta\tau}{\delta{\rho}}.
\end{split}
\]
Taking the limit $\delta{\rho} \to 0$ leads to Eq.~\eqref{eq:jump-in-quadrature-sensitivity}.
\qed
\end{pf}
\begin{remark}
When there are multiple events along the trajectory jumps in sensitivity ~\eqref{eq:jump-in-quadrature-sensitivity} will happen for each one.  The jump of the quadrature variable ${ Z}$ is governed by the values of the cost function $g$ before and after the event.
\end{remark}

\section{Direct sensitivity analysis for constrained multibody systems with smooth trajectories}
\label{sec:multibody-smooth}

This section reviews the direct sensitivity analysis for constrained systems governed by differential algebraic equations (DAEs). The presentation follows the authors' earlier work
\cite{Sandu_2013_sensitivity_ODE_multibody,Sandu_2014_sensitivity_ODE_multibody,zhu2014mbsvt,Zhu_2014_PhD,Sandu2015dynamic,Sandu_2017_vehicle-optimization}.

\subsection{Representation of constrained multibody systems}
Constrained multibody systems must satisfy the following kinematic constraints:
\begin{subequations}
\label{eq:ConstraintsEq}
\begin{eqnarray}
\label{eq:ConstraintsEq-position}
\mathsf{0}  &=& { \Phi}, \\
\label{eq:ConstraintsEq-velocity}
\mathsf{0} &=& {\dot \bPhi} ={\bPhi}_{\bq} \, \dbq + { \Phi}_t \quad \Rightarrow \quad  {\bPhi}_{\bq} {v} = -{ \Phi}_t, \\
\label{eq:ConstraintsEq-acceleration}
\mathsf{0} &=& {\ddot \bPhi} ={\bPhi}_{\bq} \, \ddbq + {\bPhi}_{\bq , \, \bq} \, ( \dbq, \dbq ) + \bPhi_{t, \, \bq } \, \dbq + { \Phi}_{t, \, t} \quad \Rightarrow \quad  {\bPhi}_{\bq} \, {\dot \bv} = -  ( {\bPhi}_{\bq}  \, {v})\, {v}  - \bPhi_{t, \, \bq } \, {v} - { \Phi}_{t, \, t} := \mathsf{C}.
\end{eqnarray}
\end{subequations}
Here \eqref{eq:ConstraintsEq-position} is a holonomic position constraint equation ${ \Phi}(t,{q},{\rho})=\mathsf{0}$, where ${ \Phi} : \mathds{R}^{1+n+p} \to \mathds{R}^{m}$ is a smooth `position constraint' function. The velocity \eqref{eq:ConstraintsEq-velocity} and the acceleration \eqref{eq:ConstraintsEq-acceleration} kinematic constraints are found by differentiating the position constraint with respect to time.

There are two main approaches to solve such systems, the DAE approach through direct inclusion of the algebraic constraints in the dynamics, and the ODE approach through either following locally the independent coordinates (Maggi) or through a penalty formulation.

\subsection{Direct sensitivity analysis for smooth systems in the index-3 differential-algebraic formulation}

\begin{definition}[Constrained multibody dynamics: the index-3 DAE formulation]
A constrained rigid multibody dynamics system is described by the following index-3 differential-algebraic equations (DAEs) \cite{Sandu_2014_sensitivity_ODE_multibody}:
\begin{equation}
\label{eq:EOM-DAE-index3}
\begin{cases}
\dbq &= {v}, \\
{\mathsf{M}}\left(t,{q},{\rho}\right) \cdot {\dot \bv} &={\mathsf{F}} \left(t,{q},{v},{\rho}\right) + {\bPhi}_{\bq}^{\rm T}\left(t,{q},{\rho}\right)\cdot\lambda, \\
{ \Phi}\left(t,{q},{\rho}\right) &= \mathsf{0},
\end{cases}
\quad t_0 \leq t \leq  t_F,  \quad {q}(t_0)={q}_0({\rho}),\quad{v}(t_0)={v}_0({\rho}).
\end{equation}
Unlike the ODE formulation \eqref{eq:EOM-ODE} the position vector of the system \eqref{eq:EOM-DAE-index3} is constrained by the equation  \eqref{eq:ConstraintsEq-position}. The joint forces ${ \Phi}_{q}^{\rm T}\,\lambda$ ensure that the system solution obeys the constraints at all points along the trajectory, and $\lambda \in \mathds{R}^{m}$ are Lagrange multipliers  associated with the position constraint \eqref{eq:ConstraintsEq-position}.
\end{definition}

Sensitivities of the position and velocity state variables are defined in \eqref{eq:sensitivity-of-solutions}. In addition, we need to consider the sensitivity of the Lagrange multipliers with respect to system parameters:
\begin{equation}
\label{eq:sensitivity-of-multipliers}
\Lambda(t,{\rho}) := D_\brho\lambda(t) := \frac{d \lambda}{d {\rho}}(t, {\rho}) \in \mathds{R}^{m \times p}.
\end{equation}

\begin{definition}[TLM of the index-3 DAE formulation]
Sensitivities of solutions \eqref{eq:sensitivity-of-solutions} and multipliers \eqref{eq:sensitivity-of-multipliers} of the system \eqref{eq:EOM-DAE-index3} with respect to parameters evolve according to the tangent linear model derived in \cite{Sandu_2013_sensitivity_ODE_multibody,Sandu_2014_sensitivity_ODE_multibody,zhu2014mbsvt,Zhu_2014_PhD,Sandu2015dynamic,Sandu_2017_vehicle-optimization}:
\begin{equation*}
\begin{cases}
&\dot{{ Q}} = {V}, \\
&\mathsf{M}\cdot \dot{{V}} =
\mathsf{F}_{v}\cdot {V}
-\left(\mathsf{M}_{{q}}\, {\dot{{v}}} + { \Phi}_{{q},{q}}^{\rm T}\, \lambda -\mathsf{F}_{q} \right) \cdot { Q}
- {\bPhi}_{\bq}^{\rm T}\cdot \Lambda + \mathsf{F}_{{\rho}}-\mathsf{M}_{{\rho}}\, {\dot{{v}}}
- { \Phi}_{{q},{\rho}}^{\rm T} \lambda, \\
&{\bPhi}_{\bq} \cdot { Q} = -{ \Phi}_{{\rho}},
\end{cases}
\end{equation*}
with initial conditions given by Eq. \eqref{eq:TLM-IC}.
\end{definition}

\subsection{Direct sensitivity analysis for smooth systems in the index-1 differential-algebraic formulation}

\begin{definition}[Constrained multibody dynamics: the index-1 DAE formulation]
The index-1 formulation of the equations of motion is obtained by replacing the position constraint \eqref{eq:ConstraintsEq-position} in \eqref{eq:EOM-DAE-index3} with the acceleration constraint \eqref{eq:ConstraintsEq-acceleration}:
\begin{equation}
\label{eq:EOM-DAE-index1}
\begin{bmatrix}
\mathsf{I} & \mathsf{0} & \mathsf{0} \\
\mathsf{0} & {\mathsf{M}}\left(t,{q},{\rho}\right) & {\bPhi}_{\bq}^{\rm T}\left(t,{q},{\rho}\right) \\
\mathsf{0} & {\bPhi}_{\bq}\left(t,{q},{\rho}\right) & \mathsf{0}
\end{bmatrix}
\cdot
\begin{bmatrix}
\dbq \\ {\dot \bv} \\ \lambda
\end{bmatrix}
=
\begin{bmatrix}
{v} \\
{\mathsf{F}} \left(t,{q},{v},{\rho}\right)  \\
\mathsf{C}\left(t,{q},{v},{\rho}\right)
\end{bmatrix},
\quad t_0 \leq t \leq  t_F,  \quad {q}(t_0)={q}_0({\rho}),\quad{v}(t_0)={v}_0({\rho}),
\end{equation}
or equivalently,
\begin{equation}
\label{eq:DAE-index1}
\dbq = {v}, \qquad
\begin{bmatrix}
{\dot \bv} \\ \lambda
\end{bmatrix}
=
\begin{bmatrix}
{\mathsf{M}} & {\bPhi}_{\bq}^{\rm T} \\
{\bPhi}_{\bq} & \mathsf{0}
\end{bmatrix}^{-1}
\cdot
\begin{bmatrix}
{\mathsf{F}}   \\
\mathsf{C}
\end{bmatrix}=
\begin{bmatrix}
{ \fun^{\scalebox{0.4}{\rm DAE-{\scalebox{1.5}{${\dot \bv}$}}}}}   \\
{ \fun^{\scalebox{0.4}{\rm DAE-{\scalebox{1.5}{$\lambda$}}}}}
\end{bmatrix}=
\fdae {\left(\,t,\,\bq,\,\bv,\,\brho \, \right)}.
\end{equation}
The algebraic equation has the form ${ \fun^{\scalebox{0.4}{\rm DAE-{\scalebox{1.5}{$\lambda$}}}}} -\lambda=0$.
\end{definition}

\begin{definition}[TLM of the index-1 DAE formulation]
Sensitivities of solutions \eqref{eq:sensitivity-of-solutions} and multipliers \eqref{eq:sensitivity-of-multipliers} of the system  \eqref{eq:EOM-DAE-index1} with respect to parameters evolve according to the tangent linear model derived in \cite{Sandu_2013_sensitivity_ODE_multibody,Sandu_2014_sensitivity_ODE_multibody,zhu2014mbsvt,Zhu_2014_PhD,Sandu2015dynamic,Sandu_2017_vehicle-optimization}:
\begin{equation*}
\label{eq:TLM-DAE-index1}
\begin{bmatrix}
\mathsf{I} & \mathsf{0} & \mathsf{0} \\
\mathsf{0} & {\mathsf{M}} & {\bPhi}_{\bq}^{\rm T} \\
\mathsf{0} & {\bPhi}_{\bq} & \mathsf{0}
\end{bmatrix}
\cdot
\begin{bmatrix}
\dot{{ Q}} \\ \dot{{V}} \\ \Lambda
\end{bmatrix}
=
\begin{bmatrix}
{V} \\
\mathsf{F}_{v}\cdot  {V}
-\left(\mathsf{M}_{{q}} \, {\dot \bv} + {\bPhi}_{\bq , \, \bq}^{\rm T}\, \lambda  -\mathsf{F}_{q} \right) \cdot { Q}
+ \mathsf{F}_{{\rho}}-\mathsf{M}_{{\rho}} \,\dot{{v}}
- { \Phi}_{{q}, \, {\rho}}^{\rm T}\, \lambda\\
\mathsf{C}_{v} \cdot {V}
-\left({\bPhi}_{\bq , \, \bq} \, {\dot \bv}  -\mathsf{C}_{{q}} \right) \cdot { Q} +
\mathsf{C}_{{\rho}} - \bPhi_{\bq, \, \brho} \, {\dot \bv}
\end{bmatrix},
\end{equation*}
with initial conditions given by Eq. \eqref{eq:TLM-IC}.
\end{definition}

\begin{definition}[Cost function]
\label{def:CostFunction-DAE}
Following Definition \ref{def:CostFunction-DAE}, consider a smooth scalar ``trajectory cost function'' ${g}$ and a smooth scalar ``terminal cost function'' $w$. A general cost function has the form:
\begin{eqnarray*}
\label{eq:CostFunction-DAE}
&&\psi=  \int_{t_0}^{t_F} {\bg\fincostdae \  dt} + w \bigl(t_F,\bq(t_F,\brho),\bv(t_F,\brho),\brho \bigr).
\end{eqnarray*}
\end{definition}
Note that the trajectory cost function \eqref{eq:CostFunction} depends on both accelerations ${\dot \bv}$ and on the Lagrange multipliers $\lambda$. These are not independent variables and they can be resolved in terms of positions and velocities using \eqref{eq:EOM-DAE-index1}, to obtain an equivalent regular trajectory cost function:
\[
{g}\bigl(\,t,\,{q},\,{v},\,{\dot \bv}(t,{q},{v},{\rho}),\,\lambda(t,{q},{v},{\rho}),\,{\rho}\,\bigr) = \tilde{g}\bigl(\,t,\,{q},\,{v},\,{\rho}\,\bigr).
\]
We keep accelerations and Lagrange multipliers (constraint forces) as explicit parameters in the cost function \eqref{eq:CostFunction-DAE} in order to give additional flexibility in practical applications. In addition, we define the `quadrature' variable $z(t) \in \mathds{R}$ as follows:
\begin{subequations}
\label{eqn:DAE-quadrature-variable}
\begin{eqnarray}
\label{eqn:DAE-quadrature-variable-int}
{z}(t,{\rho}) &:=&  \int_{t_0}^{t} { \tilde{g}\bigl(t,{q}(t,{\rho}),{v}(t,{\rho}),{\rho}\bigr) \,  dt} \quad \Leftrightarrow \quad \\
\label{eqn:DAE-quadrature-variable-dif}
{\dot \bz}(t,{\rho}) &=&  {g}\bigl(\,t,\,{q},\,{v},\,{\dot \bv},\,\lambda,\,{\rho}\,\bigr)  = \tilde{g}\bigl(t,{q}(t,{\rho}),{v}(t,{\rho}), {\rho}\bigr),  \quad t_0 \leq t \leq  t_F, \quad z(t_0,{\rho})=0.
\end{eqnarray}
\end{subequations}

\begin{definition}[The DAE quadrature sensitivity]
Similarly, let the `quadrature sensitivity' vector $Z(t,{\rho})$ be the Jacobian of the `quadrature' variable $z(t,{\rho})$ Eq. ~\eqref{eqn:DAE-quadrature-variable-int} with respect to the parameters ${\rho}$:
\[
Z_i(t,{\rho}) := \frac{\partial z(t,{\rho})}{\partial {\rho}_i}, ~~ i=1,\dots,p; \quad
Z(t,{\rho}) := \nabla_{\rho} z(t,{\rho}) = \begin{bmatrix} Z_1(t,{\rho}) \cdots Z_p(t,{\rho}) \end{bmatrix} \in \mathds{R}^{1 \times p}.
\]
The time evolution equations of the quadrature sensitivities are given by the TLM obtained by differentiating
\eqref{eqn:DAE-quadrature-variable-dif} with respect to the parameters:
\begin{equation}
\label{eq:TLM-quadrature}
\begin{split}
{\dot \bZ}_i & =  \frac{d\, {g}\bigl(\,t,\,{q},\,{v},\,{\dot \bv},\,\lambda,\,{\rho}\,\bigr)}{d\,{\rho}_i} \\
&= g_{q}\cdot { Q}_{i} +  g_{v}\cdot {V}_{i}  + g_{{\dot \bv}}\cdot \frac{d { \fun^{\scalebox{0.4}{\rm DAE-{\scalebox{1.5}{${\dot \bv}$}}}}}}{d {\rho}}
+g_{\lambda}\cdot \frac{d { \fun^{\scalebox{0.4}{\rm DAE-{\scalebox{1.5}{$\lambda$}}}}}}{d {\rho}} + g_{{\rho}_i} \\
&= \big(g_{q} +  g_{{\dot \bv}} \cdot { \fun_\bq^{\scalebox{0.4}{\rm DAE-{\scalebox{1.5}{${\dot \bv}$}}}}} +g_{\lambda} \cdot { \fun_\bq^{\scalebox{0.4}{\rm DAE-{\scalebox{1.5}{$\lambda$}}}}} \big)\cdot { Q}_{i}
+  \big(g_{v} +  g_{{\dot \bv}} \cdot { \fun_\bv^{\scalebox{0.4}{\rm DAE-{\scalebox{1.5}{${\dot \bv}$}}}}} +g_{\lambda} \cdot { \fun_\bv^{\scalebox{0.4}{\rm DAE-{\scalebox{1.5}{$\lambda$}}}}} \big)\cdot {V}_{i}
+   g_{{\dot \bv}}\cdot{ \fun_\bq^{\scalebox{0.4}{\rm DAE-{\scalebox{1.5}{$\lambda$}}}}} +g_{\lambda}\cdot{ \fun_\bv^{\scalebox{0.4}{\rm DAE-{\scalebox{1.5}{$\lambda$}}}}}+g_{{\rho}_i} , \\
& i=1,\dots,p, \quad t_0 \le t \le t_F, \quad  { Z}_i(t_0,{\rho}) = 0.
\end{split}
\end{equation}
\end{definition}

\begin{definition} [Canonical index-1 sensitivity DAE]
\label{def:canonical-dae-sensitivity}
The canonical DAE system for the solution given by \eqref{eq:DAE-index1}, the DAE TLM given by \eqref{eq:TLM-DAE-index1}, and the sensitivity quadrature equations given by \eqref{eq:TLM-quadrature} need to be solved together forward in time, leading to the canonical sensitivity DAE that computes the derivatives of the cost function with respect to the system parameters ${\rho}$ for smooth systems:
\begin{eqnarray}
\label{eq:canonical-DAE-sensitivity}
\begin{bmatrix}
\dbq \\ {\dot \bv} \\ \lambda \\ {\dot \bz} \\
\big[ \dot{{ Q}}_i  \big]_{i=1,\dots,p} \\
\big[ {\dot{\bV}}_i  \big]_{i=1,\dots,p} \\
\big[ \Lambda_i  \big]_{i=1,\dots,p} \\
\big[ {\dot \bZ}_i  \big]_{i=1,\dots,p}
\end{bmatrix}=
\begin{bmatrix}
{v} \\
{ \fun^{\scalebox{0.4}{\rm DAE-{\scalebox{1.5}{${\dot \bv}$}}}}}\\
{ \fun^{\scalebox{0.4}{\rm DAE-{\scalebox{1.5}{$\lambda$}}}}}\\
 \tilde{{g}} \\
\left[  {V}_i  \right]_{i=1,\dots,p}  \\
\left[ { \fun_\bq^{\scalebox{0.4}{\rm DAE-{\scalebox{1.5}{${\dot \bv}$}}}}} { Q}_i +{ \fun_\bv^{\scalebox{0.4}{\rm DAE-{\scalebox{1.5}{${\dot \bv}$}}}}} {V}_i +{ \fun_{\brho_i}^{\scalebox{0.4}{\rm DAE-{\scalebox{1.5}{${\dot \bv}$}}}}}  \right]_{i=1,\dots,p} \\
\left[ { \fun_\bq^{\scalebox{0.4}{\rm DAE-{\scalebox{1.5}{$\lambda$}}}}} { Q}_i +{ \fun_\bv^{\scalebox{0.4}{\rm DAE-{\scalebox{1.5}{$\lambda$}}}}} {V}_i +{ \fun_{\brho_i}^{\scalebox{0.4}{\rm DAE-{\scalebox{1.5}{$\lambda$}}}}}  \right]_{i=1,\dots,p} \\
\left[  \big(g_{q} +  g_{{\dot \bv}} \, { \fun_\bq^{\scalebox{0.4}{\rm DAE-{\scalebox{1.5}{${\dot \bv}$}}}}} +g_{\lambda} \, { \fun_\bq^{\scalebox{0.4}{\rm DAE-{\scalebox{1.5}{$\lambda$}}}}} \big)\cdot { Q}_{i}
+  \big(g_{v} +  g_{{\dot \bv}} \, { \fun_\bv^{\scalebox{0.4}{\rm DAE-{\scalebox{1.5}{${\dot \bv}$}}}}} +g_{\lambda} \, { \fun_\bv^{\scalebox{0.4}{\rm DAE-{\scalebox{1.5}{$\lambda$}}}}} \big)\cdot {V}_{i}
+   g_{{\dot \bv}}\,{ \fun_\bq^{\scalebox{0.4}{\rm DAE-{\scalebox{1.5}{$\lambda$}}}}} +g_{\lambda}\,{ \fun_\bv^{\scalebox{0.4}{\rm DAE-{\scalebox{1.5}{$\lambda$}}}}} +g_{{\rho}_i}  \right]_{i=1,\dots,p}
\end{bmatrix},
\end{eqnarray}
where the state vector of the canonical index-1 sensitivity DAE \eqref{eq:canonical-DAE-sensitivity} is :
\begin{eqnarray}
\label{eq:canonical-DAEsensitivity-state}
{X} = \left[   \, {q}^{\rm T}, \, {v}^{\rm T}, \lambda^{\rm T}, \, {z}, ~
{ Q}_1^{\rm T},  \dots , { Q}_p^{\rm T}, \,
{V}_1^{\rm T},  \dots , {V}_p^{\rm T}, \,
\Lambda_1^{\rm T},  \dots , \Lambda_p^{\rm T}, \,
{ Z}_1, \dots, { Z}_p
\right] ^{\rm T}  \in \mathds{R}^{(n+1)(p+1)},
\end{eqnarray}
and the derivatives of the DAE function are:
\[
{ \fun_\bq^{\scalebox{0.4}{\rm DAE}}}=
\begin{bmatrix}
{\mathsf{M}} & {\bPhi}_{\bq}^{\rm T} \\
 {\bPhi}_{\bq} & \mathsf{0}
\end{bmatrix}^{-1}
\begin{bmatrix}
\mathsf{F}_{q} -\mathsf{M}_{{q}} \, {\dot \bv} - {\bPhi}_{\bq , \, \bq} ^{\rm T}\, \lambda
\\
\mathsf{C}_{{q}}- {\bPhi}_{\bq , \, \bq} \dot{{v}}
\end{bmatrix}
, \;
{ \fun_\bv^{\scalebox{0.4}{\rm DAE}}}= \begin{bmatrix}
{\mathsf{M}} & {\bPhi}_{\bq}^{\rm T} \\
 {\bPhi}_{\bq} & \mathsf{0}
\end{bmatrix}^{-1}
\begin{bmatrix}
\mathsf{F}_{v} \\
\mathsf{C}_{v}
\end{bmatrix}
,\;
{ \fun_\brho^{\scalebox{0.4}{\rm DAE}}}=\begin{bmatrix}
{\mathsf{M}} & {\bPhi}_{\bq}^{\rm T} \\
 {\bPhi}_{\bq} & \mathsf{0}
\end{bmatrix}^{-1}
\begin{bmatrix}
\mathsf{F}_{{\rho}}-\mathsf{M}_{{\rho}} \, {\dot \bv}
- \bPhi_{\bq, \, \brho}^{T}\, \lambda\\
\mathsf{C}_{{\rho}} - \bPhi_{\bq, \, \brho} \, {\dot \bv}
\end{bmatrix}.
\]
\end{definition}

\subsection{Direct sensitivity analysis for smooth systems in the penalty ODE formulation}
\begin{definition}[Constrained multibody dynamics: the penalty ODE formulation]
Define the extended mass matrix $\overline{\mathsf{M}}:\mathds{R} \times \mathds{R}^{n} \times \mathds{R}^{n} \times \mathds{R}^{p} \rightarrow \mathds{R}^{n \times n}$ as:
\begin{subequations}
\label{eq:EOM-ODE-penalty}
\begin{equation}
\label{eq:EOM-ODE-penalty-mass}
\overline{\mathsf{M}}\left(t,{q},{v},{\rho}\right) := \mathsf{M}\left(t,{q},{v},{\rho}\right) +{\bPhi}_{\bq}^{\rm T}\left(t,{q},{v},{\rho}\right)\cdot \alpha\cdot {\bPhi}_{\bq}\left(t,{q},{v},{\rho}\right),
\end{equation}
where $\alpha \in \Re^{m \times m}$ is the penalty factor of the ODE penalty formulation.
Define the extended right hand side function $\overline{\mathsf{F}}:\mathds{R} \times \mathds{R}^{n} \times \mathds{R}^{n} \times \mathds{R}^{p} \rightarrow \mathds{R}^{n}$ as:
\begin{equation}
\label{eq:EOM-ODE-penalty-force}
\overline{\mathsf{F}}\left(t,{q},{v},{\rho}\right) := {\mathsf{F}}\left(t,{q},{v},{\rho}\right)
-{\bPhi}_{\bq}^{\rm T}\cdot \alpha\cdot \left({\dPhi}_{\bq}\, {v} +
{\dot \bPhi}_{t}+2 \, \xi\, \omega\, {\dot \bPhi}+\omega^2 {{ \Phi}}\right),
\end{equation}
where $\xi\in \Re$ and $\omega\in \Re$ are the natural frequency and damping ratio coefficients of the formulation, respectively, and  ${\dot \bPhi}$ is the total time derivative of the kinematic constraints. The algebraic position constraints \eqref{eq:ConstraintsEq-position} are removed and an auxiliary spring-damper force is added in \eqref{eq:EOM-ODE-penalty-force} to prevent the system from deviating away from the constraints.

In the penalty formulation the EOM of a constrained rigid multibody system is the second order ODE:
\begin{equation}
\begin{cases}
\dbq &= {v}, \\
{\dot \bv} &= {\feom} {\left(\,t,\,\bq,\,\bv,\,\brho \, \right)} = \overline{\mathsf{M}}^{-1}{\left(\,t,\,\bq,\,\bv,\,\brho \, \right)} \cdot
\overline{\mathsf{F}}{\left(\,t,\,\bq,\,\bv,\,\brho \, \right)}.
\end{cases}
\end{equation}
The Lagrange multipliers associated to the constraint forces are estimated as follows:
\begin{equation}
\label{eq:lambda}
{\lambda}^{*}=\alpha\, \left(\, \ddot{{ \Phi}}+
2 \, \xi\, \omega\, {\dot \bPhi}+\omega^2 \,{ \Phi} \, \right)\,.
\end{equation}
\end{subequations}
\end{definition}
The cost function \eqref{eq:CostFunction-DAE} is formulated using the Lagrange multiplier estimates \eqref{eq:lambda}, i.e., using the trajectory cost function ${g}\left(\,t,\,{q},\,{v},\,{\dot \bv},\,\lambda^{*},\,{\rho}\,\right)$.
Sensitivities \eqref{eq:sensitivity-of-solutions} of the position and velocity state variables of the system \eqref{eq:EOM-ODE-penalty} with respect to parameters evolve according to the tangent linear model derived in \cite{Sandu_2013_sensitivity_ODE_multibody,Sandu_2014_sensitivity_ODE_multibody,zhu2014mbsvt,Zhu_2014_PhD,Sandu2015dynamic,Sandu_2017_vehicle-optimization}:
\begin{equation}
\label{eq:TLM-ODE-penalty}
\begin{cases}
\dot{{ Q}} &= {V}, \\
\overline{\mathsf{M}}\cdot \dot{{V}} &=
\left( \overline{\mathsf{F}}_{q}  - \overline{\mathsf{M}}_{q} \, {\dot \bv}\right)\cdot { Q}
+\overline{\mathsf{F}}_{v}\cdot  {V}
+ \overline{\mathsf{F}}_{\rho} -\overline{\mathsf{M}}_{{\rho}} \, {\dot \bv},
\end{cases}
\qquad t_0 \le t \le t_F,
\end{equation}
with initial conditions given by Eq.~\eqref{eq:TLM-IC}. The derivatives
$\overline{\mathsf{F}}_{q}, \overline{\mathsf{F}}_{v}, \overline{\mathsf{F}}_{\rho}, \overline{\mathsf{M}}_{{\rho}},$ and $ \overline{\mathsf{M}}_{q}$ are given in  \ref{sec:AppendixA}.

\begin{definition}[Canonical ODE sensitivity system]
The canonical sensitivity ODE that computes the derivatives of the cost function with respect to the system parameters ρ for the smooth ODE penalty system Eq.~\eqref{eq:EOM-ODE-penalty} is the same than the ODE canonical system presented in \eqref{eq:canonical-ode-sensitivity} and extended to the cost function \eqref{eq:CostFunction-DAE} formulated using the Lagrange multiplier estimates.
\begin{eqnarray}
\label{eq:canonical-ode-penalty-sensitivity}
\begin{bmatrix}
\dbq \\ {\dot \bv} \\ {\dot \bz} \\ \big[ \dot{{ Q}}_i \big]_{i=1,\dots,p}\\
\big[ \dot{{V}}_i \big]_{i=1,\dots,p}\\
\big[ {\dot \bZ}_i \big]_{i=1,\dots,p}
\end{bmatrix} =
\begin{bmatrix}
{v} \\
\feom\\
 \tilde{{g}}\\
\big[  {V}_i \big]_{i=1,\dots,p}\\
\big[ { \fun_\bq^{\scalebox{0.5}{\rm eom}}} \cdot { Q}_i +
{ \fun_\bv^{\scalebox{0.5}{\rm eom}}} \cdot {V}_i +
{\fun_{\brho_i}^{\scalebox{0.5}{\rm eom}}} \big]_{i=1,\dots,p} \\
\big[ \big(g_{q} +  g_{\dot \bv} \cdot { \fun_\bq^{\scalebox{0.5}{\rm eom}}} + g_{{\lambda}^{*}} \cdot {\lambda}_{{q}}^{*} \big)\cdot { Q}_{i}
+  \big(g_{v} +  g_{\dot \bv} \cdot { \fun_\bv^{\scalebox{0.5}{\rm eom}}} +g_{{\lambda}^{*}} \cdot {\lambda}_{{v}}^{*} \big)\cdot {V}_{i}+
\big(g_{{\rho}_i} +  g_{\dot \bv}\cdot{\fun_{\brho_i}^{\scalebox{0.5}{\rm eom}}} +  g_{\dot \bv}\cdot{\lambda}_{{\rho}}^{*}\big) \big]_{i=1,\dots,p}
\end{bmatrix},
\end{eqnarray}
where
\[
{ \fun_\bq^{\scalebox{0.5}{\rm eom}}}=\overline{\mathsf{M}}^{\rm -1}\left(\overline{\mathsf{F}}_{q}-
\overline{\mathsf{M}}_{{q}} {\dot \bv}\right), \quad
{ \fun_\bv^{\scalebox{0.5}{\rm eom}}}=\overline{\mathsf{M}}^{\rm -1}\overline{\mathsf{F}}_{v},\quad
{\fun_{\brho_i}^{\scalebox{0.5}{\rm eom}}}=\overline{\mathsf{M}}^{\rm -1}\left(\overline{\mathsf{F}}_{\rho} - \overline{\mathsf{M}}_{\rho} {\dot \bv} \right),
\]
and with the initial conditions given by Eq. \eqref{eq:TLM-IC}. The derivatives ${\lambda}_{{q}}^{*}, {\lambda}_{{v}}^{*}$ and  $ {\lambda}_{{\rho}}^{*}$ are given in \ref{sec:AppendixA}.
\end{definition}
\begin{remark} The sensitivity of the estimated Lagrange multipliers
	%
\begin{equation}
\label{eq:sensitivity-of-estimated-multipliers}
\Lambda^*(t,{\rho}) := D_\brho\lambda^*(t) := \frac{d \lambda^*}{d {\rho}}(t, {\rho}) \in \mathds{R}^{m \times p}
\end{equation}
is calculated as:
		\begin{equation}
		\Lambda_i^{*}  =  {\lambda}_{{q}}^{*}  \; { Q}_i +{\lambda}_{{v}}^{*}\; {V}_i +{\lambda}_{{\rho}_i}^{*}, \quad i=1,\dots,p.
		\end{equation}
	\end{remark}

\section{Direct sensitivity analysis for hybrid constrained multibody systems}
\label{sec:multibody-constrained}

We now discuss constrained multibody systems when the dynamics is piecewise continuous in time.

\subsection{Coordinates partitioning for hybrid multibody systems}

The direct sensitivity analysis for a constrained rigid hybrid multibody dynamic system requires to find the jump conditions at the time of event. For this we need to distinguish between dependent and independent state variables and their sensitivities.

Assume that the Jacobian of the position constraint \eqref{eq:ConstraintsEq-position} has full row rank at a given configuration, $\operatorname{rank}({ \Phi}_{q})=m$. One can rearrange the columns and split the Jacobian in two submatrices:
\begin{equation}
\label{eqn:split-constraint-jacobian}
{ \Phi}_{q} \cdot \mathsf{P}^{\rm T} = \big[ { \Phi}_{{q}_{\rm dep}} ~~ { \Phi}_{{q}_{\rm dof}} \big], \quad
 { \Phi}_{{q}_{\rm dep}} \in \mathds{R}^{m \times m}, \quad
 { \Phi}_{{q}_{\rm dof}} \in \mathds{R}^{m \times f}, \quad f = n-m,
\end{equation}
such that the first block ${ \Phi}_{{q}_{\rm dep}}$ is nonsingular. Here $\mathsf{P} \in \mathds{R}^{n \times n}$ is a permutation matrix, obtained by permuting rows of identity matrix; the multiplication ${ \Phi}_{q} \cdot \mathsf{P}$ performs a permutation of the columns of ${ \Phi}_{q}$.

By the implicit function theorem one can partition locally the position state variables into independent coordinates ${q}_{\rm dof}  \in \mathds{R}^{f}$ (the local `degrees of freedom' of the system) and dependent coordinates ${q}_{\rm dep} \in \mathds{R}^{m}$, and solve for the dependent ones in terms of the degrees of freedom:
\[
{ \Phi}(t,{q}) = \mathsf{0} \quad \textnormal{and} \quad { \Phi}_{{q}_{\rm dep}}(t,{q})~ \textnormal{nonsingular}
\quad \Rightarrow \quad
{q}_{\rm dep} = \zeta\bigl( t,{q}_{\rm dof}  \bigr).
\]
This induces a  corresponding local partitioning of the state variables into independent components ${q}_{\rm dof},{v}_{\rm dof}  \in \mathds{R}^{f}$ and dependent components ${q}_{\rm dep},{v}_{\rm dep} \in \mathds{R}^{m}$:
\begin{equation}
\label{eq:coordinate-partitioning}
\mathsf{P} \cdot {q} =  \begin{bmatrix} \mathsf{P}_{\rm dep} \\ \mathsf{P}_{\rm dof}\end{bmatrix}  \cdot {q} = \begin{bmatrix} {q}_{\rm dep} \\ {q}_{\rm dof} \end{bmatrix}, \qquad
\mathsf{P} \cdot {v} =  \begin{bmatrix} \mathsf{P}_{\rm dep}  \\ \mathsf{P}_{\rm dof}  \end{bmatrix} \cdot {v} = \begin{bmatrix} {v}_{\rm dep} \\ {v}_{\rm dof} \end{bmatrix},
\end{equation}
where $\mathsf{P}_{\rm dep} \in \mathds{R}^{m \times n}$ and $\mathsf{P}_{\rm dof} \in \mathds{R}^{f \times n}$ consist the first $m$ and the last $f$ rows of $\mathsf{P}$, respectively.
Let:
\begin{equation}
\label{eq:R-matrix}
\mathsf{R} :=  -{ \Phi}_{{q}_{\rm dep}}^{-1}\,{ \Phi}_{{q}_{\rm dof}}  \in \mathds{R}^{m \times f}.
\end{equation}
The velocity constraint equation \eqref{eq:ConstraintsEq-velocity} becomes:
\begin{equation}
\label{eq:ConstraintsEq-velocity-dof}
 { \Phi}_{{q}_{\rm dep}}\, {v}_{\rm dep} + { \Phi}_{{q}_{\rm dof}}\, {v}_{\rm dof} = -{ \Phi}_t
\quad \Rightarrow \quad
 {v}_{\rm dep} = -{ \Phi}_{{q}_{\rm dep}}^{-1}\, \bigl( { \Phi}_{{q}_{\rm dof}}\, {v}_{\rm dof} +  { \Phi}_t \bigr)
 = \mathsf{R}\,{v}_{\rm dof}  - { \Phi}_{{q}_{\rm dep}}^{-1}\,{ \Phi}_t.
\end{equation}
Similarly, the acceleration constraint equation \eqref{eq:ConstraintsEq-acceleration} becomes:
\begin{equation}
\label{eq:ConstraintsEq-acceleration-dof}
{ \Phi}_{{q}_{\rm dep}}\, {\dot \bv}_{\rm dep} + { \Phi}_{{q}_{\rm dof}}\, {\dot \bv}_{\rm dof} = \mathsf{C}
\quad \Rightarrow \quad
{\dot \bv}_{\rm dep} = -{ \Phi}_{{q}_{\rm dep}}^{-1}\, \bigl( { \Phi}_{{q}_{\rm dof}}\, {\dot \bv}_{\rm dof} - \mathsf{C} \bigr)
= \mathsf{R}\,{\dot \bv}_{\rm dof} + { \Phi}_{{q}_{\rm dep}}^{-1}\,\mathsf{C}.
\end{equation}
From \eqref{eqn:split-constraint-jacobian}, \eqref{eq:coordinate-partitioning}, and \eqref{eq:R-matrix} we have that:
\[
{ \Phi}_{{q}_{\rm dof}}  = { \Phi}_{q} \cdot \mathsf{P}_{\rm dof}^{\rm T} \in \mathds{R}^{m \times f}, \quad
{ \Phi}_{{q}_{\rm dep}}  = { \Phi}_{q} \cdot \mathsf{P}_{\rm dep}^{\rm T} \in \mathds{R}^{m \times m}, \quad
\mathsf{R} = -\bigl(  { \Phi}_{q} \cdot \mathsf{P}_{\rm dep}^{\rm T}  \bigr)^{-1}\cdot { \Phi}_{q} \cdot \mathsf{P}_{\rm dof}^{\rm T}.
\]

\subsection{Representation of constrained hybrid  multibody systems}

The hybrid dynamics of a constrained mechanical system  refers to the smooth system defined in Section ~\ref{sec:multibody-smooth} subjected to a finite number of events, as discussed in Definition \ref{def:time-of-invent}. Each event \eqref{eq:event_function} happening at the `time of event' ${t_{\rm eve}}$ introduces a kink in the trajectory of the mechanical system. At each event the velocity state vector of an {\it unconstrained} system undergoes a jump \eqref{eq:jump-in-velocity} that can be arbitrary, i.e., can be described by any smooth function $h(\cdot)$. In case of a {\it constrained} system we need a more comprehensive understanding of the event.

\begin{definition}[Characterization of an event for constrained multibody systems]
\label{def:characterize-impulsive jump event}
During an event at time ${t_{\rm eve}}$ a constrained mechanical system undergoes a sudden change in state characterized as follows:
\begin{itemize}
\item The constraints may change at the time of event (e.g., when a walking humanoid robot changes its supporting foot at each step). Consequently, the position constraint function \eqref{eq:ConstraintsEq-position}  changes from ${ \Phi}^- : \mathds{R}^{1+n+p} \to \mathds{R}^{m-}$ before invent to ${ \Phi}^+ : \mathds{R}^{1+n+p} \to \mathds{R}^{m+}$ after invent:
\begin{equation*}
{ \Phi}^-(t,{q},{\rho}) \stackrel{\rm event}{\longrightarrow} { \Phi}^+(t,{q},{\rho}).
\end{equation*}
The two constraint functions are different, and in particular the number of constraints can differ, $m^+ \neq m^-$.
\item The Jacobians of the position constraints  before and after the invent have full row ranks at the invent configuration $\bq|_{\timp}$:
\[
\operatorname{rank}\big({ \Phi}_{q}^-(t,{q},{\rho})\big)=m^{-}, \quad
\operatorname{rank}\big({ \Phi}_{q}^+(t,{q},{\rho})=m^{+}.
\]
\item Since the constraints can be different after and before the event, the partitions of variables into independent and dependent can also differ. We denote by $\square_\textnormal{dof-},\square_\textnormal{dep-}$ the independent and dependent components before the event, and by $\square_\textnormal{dof+},\square_\textnormal{dep+}$ the independent and dependent components after the event:
\begin{equation}
\label{eq:velocity-partitioning}
\mathsf{P}^{-} \cdot {v} = \begin{bmatrix} {v}_\textnormal{dep-} \\ {v}_\textnormal{dof-} \end{bmatrix} \in \Re^n, \quad {v}_\textnormal{dof-} \in \Re^{f-},
\qquad
\mathsf{P}^{+} \cdot {v} =   \begin{bmatrix} {v}_\textnormal{dep+} \\ {v}_\textnormal{dof+} \end{bmatrix} \in \mathds{R}^n,
\quad{v}_\textnormal{dof+}\in \mathds{R}^{f+}.
\end{equation}
Here $\mathsf{P}^{-}$ and $\mathsf{P}^{+}$ are the permutation matrices that select the dependent and independent coordinates before and after the event, respectively. The dimensions of the velocity degrees of freedom vectors are  $f^{-} = n-m^{-}$ and $f^{+} = n-m^{+}$ before and after the event, respectively.

\item The generalized position state variables remain the same \eqref{eq:jump-in-position}, i.e., ${q}|_{\timp}^{+} = {q}|_{\timp}^{-} = \bq|_{\timp}$. Consequently, the state at the time of event $\bq|_{\timp}$ need to satisfy both constraint functions:
\begin{equation}
{ \Phi}^-|_{\timp}^{-}:={ \Phi}^-\left({t_{\rm eve}},\bq|_{\timp}, {\rho} \right) = 0,
\qquad
{ \Phi}^+|_{\timp}^{+}:={ \Phi}^+\left({t_{\rm eve}},\bq|_{\timp}, {\rho} \right) = 0.
\end{equation}
At the invent the system moves from one constraint manifold to another, and $\bq|_{\timp}$ is on the intersection of the two manifolds.

\item The jump in velocity from right before the invent to right after the invent is defined in terms of the independent components, i.e., in terms of the velocity degrees of freedom:
\begin{equation}
\label{eq:jump-in-velocity-constrained-dof}
{v}_\textnormal{dof+}|_{\timp}^{+} =h \Big({{t_{\rm eve}}} , \bq|_{\timp}, {v}_\textnormal{dof-}|_{\timp}^{-}, {{\rho}} \Big), \qquad
h : \mathds{R}^{1+n+f^-+p} \to  \mathds{R}^{f^+}.
\end{equation}
The jump function \eqref{eq:jump-in-velocity-constrained-dof} is assumed to be smooth. Note that its formulation is not unique, since it depends on the selections of the degrees of freedom that are not unique.

\item The velocity state vectors satisfy the velocity kinematic constraints \eqref{eq:ConstraintsEq-velocity-dof}. Consequently, the jumps in velocity \eqref{eq:jump-in-velocity} cannot be arbitrary for the dependent components. The dependent components of velocity are obtained from solving the velocity constraints \eqref{eq:ConstraintsEq-velocity-dof}:
\begin{equation}
\label{eq:jump-in-velocity-constrained-dep}
\begin{split}
 {v}_{\rm dep+}|_{\timp}^{+} &= -\left({ \Phi}^+_{{q}_{\rm dep+}}|_{\timp}^{+}\right)^{-1}\cdot \left( { \Phi}^+_{{q}_{\rm dof+}}|_{\timp}^{+}\, {v}_{\rm dof+}|_{\timp}^{+}  +  { \Phi}^+_t|_{\timp}^{+} \right) \\
 &= \mathsf{R}^+|_{\timp}^{+}\,  {v}_{\rm dof+}|_{\timp}^{+} -\left({ \Phi}^+_{{q}_{\rm dep+}}|_{\timp}^{+}\right)^{-1}\cdot   { \Phi}^+_t|_{\timp}^{+}.
 \end{split}
\end{equation}
Here $\mathsf{R}^\pm$ are the matrices corresponding to the constraints ${ \Phi}^\pm$.
\end{itemize}
\end{definition}

\begin{remark}[Collision events]
\label{rem:jump-in-collision}
The proposed formalism \eqref{eq:jump-in-velocity-constrained-dof}--\eqref{eq:jump-in-velocity-constrained-dep} covers the case of elastic contact/collision/impact without change in the set of constraint equations,  ${ \Phi}^+\equiv{ \Phi}^-$. The impulsive (external) contact forces act to change the independent components of the velocity state \eqref{eq:jump-in-velocity-constrained-dof}.
\end{remark}

\begin{remark}[Hybrid DAE jump formulation]
\label{rem:jump-in-robotics}
The proposed formalism \eqref{eq:jump-in-velocity-constrained-dof}--\eqref{eq:jump-in-velocity-constrained-dep} also covers the case where the event consists solely of a change of constraints ${ \Phi}^+\ne{ \Phi}^-$, without any external force to modify the independent velocities. This type of event appears mainly in the humanoid robotics field where general and relative coordinates are used and inelastic collisions are considered. A popular approach in robotics  is to solve for the DAE involving impulsive forces in the constraints at the time of event  \cite{Kolathaya2016, Posa2013}:
\begin{subequations}
\label{eq:DAE-impulse}
\begin{equation}
\label{eq:DAE-impulse-1}
\begin{bmatrix}
{\mathsf{M}|_{\timp}} & ({\bPhi}_{\bq}^+)^{\rm T}|_{\timp} \\
{\bPhi}_{\bq}^+|_{\timp} & \mathsf{0}
\end{bmatrix}
\cdot
\begin{bmatrix}
{\bv \timpplus} \\ \delta\lambda
\end{bmatrix}
=
\begin{bmatrix}
 \mathsf{M}|_{\timp} \cdot \bv \timpminus \\
-{ \Phi}_t^+|_{\timp}
\end{bmatrix},
\end{equation}
or equivalently,
\begin{equation}
\label{eq:DAE-impulse-2}
\begin{bmatrix}
{\bv \timpplus} \\ \delta\lambda
\end{bmatrix}
=
\begin{bmatrix}
{\mathsf{M}|_{\timp}} & ({\bPhi}_{\bq}^+)^{\rm T}|_{\timp} \\
{\bPhi}_{\bq}^+|_{\timp} & \mathsf{0}
\end{bmatrix}^{-1}
\cdot
\begin{bmatrix}
 \mathsf{M}|_{\timp} \cdot \bv \timpminus \\
-{ \Phi}_t^+|_{\timp}
\end{bmatrix}=
\begin{bmatrix}
\fdaeimpv {\left(\,{t_{\rm eve}},\,\bq|_{\timp} ,\,\bv \timpminus,\,{\rho} \, \right)} \\
\fdaeimpl {\left(\,{t_{\rm eve}},\,\bq|_{\timp} ,\,\bv \timpminus,\,{\rho} \, \right)}
\end{bmatrix}.
\end{equation}
\end{subequations}
Here ${\bv \timpplus}$ contains both the independent and dependent coordinates. We see that the second equation in \eqref{eq:DAE-impulse-1} automatically imposes the velocity constraint \eqref{eq:ConstraintsEq-velocity}.

Our formalism covers this approach by defining the jump function given by Eq. \eqref{eq:jump-in-velocity-constrained-dof} as:
\[
{v}_{\rm dof+}|_{\timp}^{+} = \mathsf{P}_{\rm dof+}\,\fdaeimpv {\left(\,{t_{\rm eve}},\,\bq|_{\timp} ,\,\bv \timpminus,\,{\rho} \, \right)}
=: h\left(\,{t_{\rm eve}},\,\bq|_{\timp} ,\,{v}_{\rm dof-}|_{\timp}^{-},\,{\rho} \, \right).
\]
\end{remark}

\subsection{The jump in the sensitivity of the position state vector}
The jump conditions at the time of event in the sensitivity state vector for a constrained rigid multibody involve finding the sudden change in values of the sensitivity with respect to the system parameters ${\rho}$ of the position and the dependent and independent velocity state variables due the impulsive jump of the independent velocity state variables.

\begin{remark}[Partitioning of sensitivity matrices]
The partitioning of state variables into dependent and independent \eqref{eq:velocity-partitioning} induces a similar partitioning of the sensitivity matrices \eqref{eq:sensitivity-of-solutions}:
\begin{equation}
\label{eq:sensitivity-partitioning}
\mathsf{P} \cdot { Q} = \begin{bmatrix} { Q}_\textnormal{dep} \\ { Q}_\textnormal{dof} \end{bmatrix} \in \mathds{R}^{n \times p}, \quad { Q}_\textnormal{dof} \in \mathds{R}^{f \times p}, \qquad
\mathsf{P} \cdot {V} = \begin{bmatrix} {V}_\textnormal{dep} \\ {V}_\textnormal{dof} \end{bmatrix} \in \mathds{R}^{n \times p},
\quad{V}_\textnormal{dof} \in \mathds{R}^{f \times p}.
\end{equation}
Differentiation of the position constraint equation \eqref{eq:ConstraintsEq-position} with respect to the system parameters ${\rho}$ gives:
\[
\mathsf{0}  = \frac{d { \Phi}(t,{q}(t,{\rho}),{\rho})}{d {\rho}} = { \Phi}_{q} \cdot { Q} +  { \Phi}_{\rho} =
{ \Phi}_{{q}_\textnormal{dep}} \cdot { Q}_\textnormal{dep} + { \Phi}_{{q}_\textnormal{dof}} \cdot { Q}_\textnormal{dof} +  { \Phi}_{\rho},
\]
and therefore:
\begin{equation}
\label{eq:dependent-Q}
{ Q}_\textnormal{dep}= -{ \Phi}_{{q}_\textnormal{dep}}^{-1}\, \left( { \Phi}_{{q}_\textnormal{dof}} \cdot { Q}_\textnormal{dof} +{ \Phi}_{\rho} \right) = \mathsf{R} \cdot { Q}_\textnormal{dof} -
 { \Phi}_{{q}_\textnormal{dep}}^{-1}\, { \Phi}_{\rho}.
\end{equation}
Similarly, differentiation of the velocity constraint equation \eqref{eq:ConstraintsEq-velocity} with respect to the system parameters gives:
\begin{eqnarray*}
\mathsf{0} &=& \frac{d}{d {\rho}} \left( {\bPhi}_{\bq}(t,{q}(t,{\rho}),{\rho}) \, {v}(t,{\rho}) + { \Phi}_t(t,{q}(t,{\rho}),{\rho}) \right) \\
&=& { \Phi}_{{q}} \cdot {V} + ({ \Phi}_{{q},{q}}\,{v} + { \Phi}_{{q},t}) \cdot { Q} + { \Phi}_{{\rho},{q}}\,{v} + { \Phi}_{{\rho},t} \\
&=& { \Phi}_{{q}_\textnormal{dep}} \cdot {V}_\textnormal{dep} + { \Phi}_{{q}_\textnormal{dof}} \cdot {V}_\textnormal{dof}
 + \big({ \Phi}_{{q},{q}}\,{v} +  { \Phi}_{{q},t} \big) \cdot { Q}
+  { \Phi}_{{\rho},{q}} \, {v}  + { \Phi}_{{\rho},t},
 \end{eqnarray*}
and therefore:
\begin{equation}
\label{eq:dependent-V}
{V}_\textnormal{dep} = \mathsf{R} \cdot {V}_\textnormal{dof}
- { \Phi}_{{q}_\textnormal{dep}}^{-1}
 \left(\big({ \Phi}_{{q},{q}}\,{v} +  { \Phi}_{{q},t} \big) \cdot { Q}
 +  { \Phi}_{{\rho},{q}} \, {v}  + { \Phi}_{{\rho},t}
 \right).
\end{equation}

\end{remark}

\begin{remark}[Sensitivity of the time of event for constrained systems]
The time of event depends only on the position state and on the event function \eqref{eq:event_function}. Consequently, the sensitivity of the time of event for constrained systems is the same as for unconstrained systems, and is given by Eqn.  \eqref{eq:sensitivity-time-of-invent} in Theorem \ref{thm:sensitivity-time-of-invent}.
\end{remark}

\begin{theorem}[Jump in position sensitivity for constrained system]
\label{thm:jump-in-position-sensitivity-constrained}
Let ${ Q}|_{\timp}^{+}$ and ${ Q}|_{\timp}^{-} \in \mathds{R}^{n{\times}p}$  be the sensitivities of the generalized position state vectors right after and right before the event, respectively.  The independent components of the sensitivity of the generalized positions right after the event are:
\begin{subequations}
\label{eq:jump-position-sensitivity-constrained}
\begin{equation}
\label{eq:jump-position-sensitivity-constrained-dof}
{ Q}_\textnormal{dof+}|_{\timp}^{+} = { Q}_\textnormal{dof+}|_{\timp}^{-}  - \bigg( {v}_\textnormal{dof+}|_{\timp}^{+} - {v}_\textnormal{dof+}|_{\timp}^{-} \bigg) \cdot \frac{d\timp}{d\brho}.
\end{equation}
The dependent components of the sensitivity of the generalized positions right after the event are given by equation \eqref{eq:dependent-Q}, using the after-event constraints:
\begin{equation}
\label{eq:jump-position-sensitivity-constrained-dep}
{ Q}_\textnormal{dep+}|_{\timp}^{+} =  \mathsf{R}^+|_{\timp}^{+} \cdot { Q}_\textnormal{dof+}|_{\timp}^{+} - \left. \left({ \Phi}^+_{{q}_\textnormal{dep+}}|_{\timp}^{+}\right)^{-1}\, { \Phi}^+_{\rho} \right|_{\timp}^{+}.
\end{equation}
\end{subequations}
\end{theorem}
\begin{pf}
The proof of the jump in the independent coordinates \eqref{eq:jump-position-sensitivity-constrained-dof} follows closely the proof of Theorem \ref{thm:jump-in-position-sensitivity}. The equation for dependent coordinates \eqref{eq:jump-position-sensitivity-constrained-dep} follows from the linearized position constraint equation \eqref{eq:dependent-Q}.
\qed
\end{pf}

\begin{remark}
From \eqref{eqn:split-constraint-jacobian} and \eqref{eq:coordinate-partitioning} we can rewrite \eqref{eq:jump-position-sensitivity-constrained-dof} as:
\begin{equation}
\mathsf{P}^+_\textnormal{dof+} \cdot \bigg({ Q}|_{\timp}^{+} - { Q}|_{\timp}^{-}  \bigg)
= - \mathsf{P}^+_\textnormal{dof+} \cdot \bigg( {v}|_{\timp}^{+} - \bv \timpminus\bigg) \cdot \frac{d\timp}{d\brho}.
\end{equation}
\end{remark}

\subsection{The jump in the sensitivity of the velocity state vector}
\begin{theorem}[Jump in velocity sensitivity for a constrained system]
\label{thm:jump-in-velocity-sensitivity-constrained}
Let ${V}|_{\timp}^{+}, {V}|_{\timp}^{-} \in \mathds{R}^{n{\times}p}$ be the sensitivity matrices of the generalized velocity state vectors after and before the invent, respectively. The independent coordinates of the velocity sensitivities right after the event are given by:
 %
\begin{subequations}
\label{eq:jump-in-velocity-sensitivity-constrained}
\begin{equation}
\label{eq:jump-in-velocity-sensitivity-constrained-dof}
\begin{split}
{V}_\textnormal{dof+}|_{\timp}^{+}&= h_{{q}}|_{\timp}^{-} \cdot {\bQ \beforeimp}
+h_{{v}_\textnormal{dof-}}|_{\timp}^{-} \cdot  {V}_\textnormal{\rm dof-}|_{\timp}^{-}
+\Bigl( h_{{q}}|_{\timp}^{-}\cdot \bv \timpminus-\ddbq_\textnormal{dof+}|_{\timp}^{+} +h_{{v}_\textnormal{dof-}}|_{\timp}^{-}\cdot \ddbq_\textnormal{dof-}|_{\timp}^{-} {+ h_t|_{\timp}^{-}} \Bigr) \cdot \frac{d\timp}{d\brho}
{+ h_{\rho}|_{\timp}^{-}},
\end{split}
\end{equation}
where the Jacobians of the jump function are:
\[
\begin{split}
h_{{q}}|_{\timp}^{-} &:= \frac{\partial h}{\partial{q}}\big( {q}|_{\timp},{v}_\textnormal{dof-}|_{\timp}^{-}, {\rho} \big) \in \Re^{f^+ \times n},
\qquad
h_{{v}_\textnormal{dof-}}|_{\timp}^{-} := \frac{\partial h}{\partial{v}_\textnormal{dof-}}\big( {q}|_{\timp},{v}_\textnormal{dof-}|_{\timp}^{-}, {\rho} \big) \in \Re^{f^+ \times f^-}.
\\
h_t|_{\timp}^{-} &:= \frac{\partial h}{\partial t}\big( {q}|_{\timp},{v}_\textnormal{dof-}|_{\timp}^{-}, {\rho} \big) \in \Re^{f}, \qquad \;\;\;\;
h_{\rho}|_{\timp}^{-} := \frac{\partial h}{\partial{\rho}}\big( {q}|_{\timp},{v}_\textnormal{dof-}|_{\timp}^{-}, {\rho} \big)\in \Re^{f \times p}.
\end{split}
\]
The dependent components  of the velocity sensitivities right after the event are calculated via \eqref{eq:dependent-V}, using the after-event constraints:
\begin{equation}
\label{eq:jump-in-velocity-sensitivity-constrained-dep}
{V}_\textnormal{dep+}|_{\timp}^{+} =
\left. - \big( { \Phi}^+_{{q}_\textnormal{dep+}} |_{\timp}^{+} \big)^{-1}
 \left(
{ \Phi}^+_{{q}_\textnormal{dof+}} \cdot {V}_\textnormal{dof+}
 + \big({{\bPhi}^+_{\bq , \, \bq}} \,{v} +  {\bPhi^+_{t, \, \bq }} \big) \cdot { Q}
 +  {\bPhi^+_{\bq, \, \brho}} \, {v}  + {\bPhi^+_{t,\, \brho }}
 \right)\right|_{\timp}^{+}.
\end{equation}
\end{subequations}
\end{theorem}
\begin{pf}
The proof of the jump in the independent coordinates \eqref{eq:jump-in-velocity-sensitivity-constrained-dof} follows closely the proof of Theorem \ref{thm:jump-in-velocity-sensitivity}. The equation for dependent coordinates \eqref{eq:jump-in-velocity-sensitivity-constrained-dep} follows from the linearized velocity constraint equation \eqref{eq:dependent-V}.
\qed
\end{pf}

\subsection{The jump in the sensitivity of the velocity state vector using the hybrid DAE jump formulation}

Consider the case of a sudden change in constraints discussed in Remark \ref{rem:jump-in-robotics}.
The jump in the velocity sensitivity for constrained system due to impulsive forces is determined as follows:
\[
\begin{split}
\begin{bmatrix}
{\mathsf{M}|_{\timp}^{+}} & {{\bPhi}^+_{\bq}}^{\rm T}|_{\timp}^{+} \\
{{\bPhi}^+_{\bq}}|_{\timp}^{+} & \mathsf{0}
\end{bmatrix}
\cdot
\begin{bmatrix}
{\bV \afterimp} \\ \delta\Lambda
\end{bmatrix}
=&
-\begin{bmatrix}
 {\mathsf{M}|_{\timp}^{+}}_{{q}}|_{\timp}^{+} \cdot ({\bv \timpplus} - \bv \timpminus) + {{\bPhi}^+_{\bq , \, \bq}}^{\rm T}|_{\timp}^{+} \cdot \delta\lambda \\
 {{\bPhi}^+_{\bq , \, \bq}}|_{\timp}^{+} \cdot {\bv \timpplus}
\end{bmatrix} \cdot {\bQ \afterimp} +
\begin{bmatrix}
\mathsf{M}|_{\timp}^{+} \\
\mathsf{0}
\end{bmatrix} \cdot {\bV \beforeimp} \\&
-\begin{bmatrix}
 {\mathsf{M}}_{{\rho}}|_{\timp}^{+} \cdot {\bv \timpplus} +  {\bPhi^+_{\bq, \, \brho}} ^{\rm T}|_{\timp}^{+} \cdot \delta\lambda \\
 {\bPhi^+_{\bq, \, \brho}}|_{\timp}^{+} \cdot {\bv \timpplus} + {\bPhi^+_{t, \, \bq }}|_{\timp}^{+} \cdot \bv \timpminus
 + {\bPhi^+_{t, \,\bv }}|_{\timp}^{+} \cdot \bv \timpminus  + {\bPhi^+_{t,\, \brho }}|_{\timp}^{+} \cdot \bv \timpminus
\end{bmatrix},
\end{split}
\]
or equivalently:
\[
\begin{split}
\begin{bmatrix}
{\bV \afterimp} \\ \delta\Lambda
\end{bmatrix}
=&
-
\begin{bmatrix}
{\mathsf{M}} & {{\bPhi}^+_{\bq}}^{\rm T} \\
{{\bPhi}^+_{\bq}} & \mathsf{0}
\end{bmatrix}^{-1}
\begin{bmatrix}
 {\mathsf{M}}_{{q}} \cdot ({\bv \timpplus} - \bv \timpminus) + {{\bPhi}^+_{\bq , \, \bq}}^{\rm T} \cdot \delta\lambda \\
 {{\bPhi}^+_{\bq , \, \bq}} \cdot {\bv \timpplus}
\end{bmatrix}\cdot {\bQ \afterimp} +
\begin{bmatrix}
{\mathsf{M}} & {{\bPhi}^+_{\bq}}^{\rm T} \\
{{\bPhi}^+_{\bq}} & \mathsf{0}
\end{bmatrix}^{-1}
\begin{bmatrix}
\mathsf{M} \\
\mathsf{0}
\end{bmatrix} \cdot {\bV \beforeimp}
\\& -
\begin{bmatrix}
{\mathsf{M}} & {{\bPhi}^+_{\bq}}^{\rm T} \\
{{\bPhi}^+_{\bq}} & \mathsf{0}
\end{bmatrix}^{-1}
\begin{bmatrix}
 {\mathsf{M}}_{{\rho}} \cdot {\bv \timpplus} +  {\bPhi^+_{\bq, \, \brho}}^{\rm T} \cdot \delta\lambda \\
 {\bPhi^+_{\bq, \, \brho}} \cdot {\bv \timpplus}  + {\bPhi^+_{t,\, \brho }} \cdot \bv \timpminus
\end{bmatrix}
-
\begin{bmatrix}
{\mathsf{M}} & {{\bPhi}^+_{\bq}}^{\rm T} \\
{{\bPhi}^+_{\bq}} & \mathsf{0}
\end{bmatrix}^{-1}
\begin{bmatrix}
 \mathsf{0} \\
  {\bPhi^+_{t, \, \bq }} \cdot \bv \timpminus
 + {\bPhi^+_{t, \,\bv }} \cdot \bv \timpminus
\end{bmatrix},
\end{split}
\]
which simplifies to:
\[
\begin{bmatrix}
{\bV \afterimp} \\ \delta\Lambda
\end{bmatrix}
=
{ f}_{q}^{\scalebox{0.6}{\rm DAE-imp}} \cdot {\bQ \afterimp}
+{ f}_{\bv \timpminus}^{\scalebox{0.6}{\rm DAE-imp}} \cdot {\bV \beforeimp}
+ { f}_{{\rho}}^{\scalebox{0.6}{\rm DAE-imp}}
+ { f}_{t}^{\scalebox{0.6}{\rm DAE-imp}}
\]

\subsection{The jump in the sensitivity of the Lagrange multipliers}

\begin{remark}
\label{rem:jump in DAE Lagrange multipliers}
When the DAE formalism is selected to model the smooth dynamics of a constrained mechanical system, the jump in the sensitivity of the Lagrange multipliers \eqref{eq:sensitivity-of-multipliers} from $\Lambda|_{\timp}^{-} \to \Lambda|_{\timp}^{+}$  at the time of event is:
\begin{equation}
\Lambda_i|_{\timp}^{+}  =\Lambda_i|_{\timp}^{-}  +{ \fun_\bq^{\scalebox{0.4}{\rm DAE-{\scalebox{1.5}{$\lambda$}}}}} |_{\timp}^{+} \; { Q}_i|_{\timp}^{+} +{ \fun_\bv^{\scalebox{0.4}{\rm DAE-{\scalebox{1.5}{$\lambda$}}}}}|_{\timp}^{+} \; {V}_i|_{\timp}^{+} +{ \fun_{\brho_i}^{\scalebox{0.4}{\rm DAE-{\scalebox{1.5}{$\lambda$}}}}}|_{\timp}^{+}, \quad i=1,\dots,p
\end{equation}
 \end{remark}

\begin{remark}
\label{rem:jump in estimated Lagrange multipliers}
When the ODE penalty formalism is selected to model the smooth dynamics of a constrained mechanical system, the jump in the sensitivity of the estimated Lagrange multipliers \eqref{eq:sensitivity-of-estimated-multipliers} from $\Lambda^{*}|_{\timp}^{-} \to \Lambda^{*}|_{\timp}^{+}$  at the time of event is:
\begin{equation}
\Lambda_i^{*}|_{\timp}^{+}  =\Lambda_i^{*}|_{\timp}^{-}  +{\lambda}_{{q}}^{*} |_{\timp}^{+} \; { Q}_i|_{\timp}^{+} +{\lambda}_{{v}}^{*}|_{\timp}^{+} \; {V}_i|_{\timp}^{+} +{\lambda}_{{\rho}_i}^{*}|_{\timp}^{+}, \quad i=1,\dots,p
\end{equation}
\end{remark}

\subsection{The sensitivity of the cost function for hybrid systems}

\begin{remark}
\label{rem:jump-in-the-cost-function}
The formalism that computes the sensitivities of the cost function with respect to parameters for hybrid systems does not change from the formalism presented for smooth systems illustrated in Remark \ref{rem:sensitivities-of-the-cost-function}. Indeed, the sensitivities of the cost function sum all the sensitivities of the trajectories and the quadrature variables evaluated at the final time. Any jump in the sensitivities of the trajectories and  quadrature variables were anteriorly computed. The jump in the sensitivities of the quadrature variables are given by \eqref{eq:jump-in-quadrature-sensitivity}.
\end{remark}

\section{Direct sensitivity analysis for constrained mechanical systems with transition functions}
\label{sec:Direct sensitivity analysis with transition functions}

The transition function refers to a sudden change of the governing function or vector field. In this section we discuss about direct sensitivity analysis for constrained mechanical with jump-discontinuity in the acceleration caused by a sudden change of the equation of motions at the time of event.
\begin{definition}[Change of EOM at the time of event]
Unlike previous methodology, the ODE penalty formulation Eq.~\eqref{eq:EOM-ODE-penalty} incorporates the kinematic constraints (position, velocity, and acceleration constraint equations) into the equation of motions and stabilize them over time. Therefore, any change in the set of kinematic constraints involves a change in the equation of motions, thus, a change in the acceleration vector {blue}{(or right-hand side function)}. Because  the ODE penalty formulation is a control based constraint stabilization method, the position constraint is not satisfied exactly right after the sudden change in the set of kinematic constraints:
\begin{equation}
{ \Phi}^-|_{\timp}^{-} = 0 , \quad { \Phi}^+|_{\timp}^{+} \neq 0.
\end{equation}
This differs from  \eqref{eq:DAE-impulse} as there are no instantaneous kinematic jump in the velocity state variable.
\end{definition}

\begin{theorem}[Jump in the velocity sensitivity for constrained system due to the change of equation of motions]
\label{thm:jump-in-velocity-sensitivity-constrained-EOMchanged}
Let ${V}|_{\timp}^{+}, {V}|_{\timp}^{-} \in \mathds{R}^{n{\times}p}$ be the sensitivity matrices of the generalized velocity state vectors after and before the invent, respectively. Let the event characterized as a change of equation of motions due to the change of constraints including in the equation of motions. The sensitivities of the independent  velocities right after the event are given by
\begin{equation}
{V}_\textnormal{dof-}|_{\timp}^{+} ={V}_\textnormal{dof-}|_{\timp}^{-} + \left(  \ddbq_\textnormal{dof-}|_{\timp}^{-} - \ddbq_\textnormal{dof-}|_{\timp}^{+} \right) \cdot \frac{d\timp}{d\brho}.
\end{equation}
\end{theorem}
\begin{proof}
The proof of the jump in the independent coordinates (58a) follows closely the proof of Theorem \ref{thm:ODE-change-in-acceleration}.
\end{proof}
\begin{remark}
The sensitivities of the dependent  velocities right after the event are given by \eqref{eq:jump-in-velocity-sensitivity-constrained-dep}. As well, the sensitivities of the position right after the event for the independent and dependent variables are given by \eqref{eq:jump-position-sensitivity-constrained-dof} and  \eqref{eq:jump-position-sensitivity-constrained-dep}, respectively.
\end{remark}

\begin{remark}
The sudden change in the equation of motions at the event time is caused by a sudden change of forces acting on the system, such as constraint forces, friction forces, or a change of masses. The proposed formalism to calculate the jump conditions for systems{blue}{ with discontinuous right-hand sides} remains valid for any type of change of the equation of motions.
\end{remark}

\begin{remark} The proposed formalism in calculating the jump conditions for systems with jerk discontinuity incorporates  Remarks  \ref{rem:jump in estimated Lagrange multipliers} and \ref{rem:jump-in-the-cost-function}.
\end{remark}

\section{Case study: sensitivity analysis of a five-bar mechanism}
\label{sec:numerics}
The five-bar mechanism with two degrees of freedom, shown in Fig.~\ref{Fig_fivebar}, is used as a case study to illustrate the sensitivity analysis approach for hybrid constrained multibody systems developed herein.

The mechanism has five revolute joints located at points A, 1, 2, 3, and B; the masses of the bars are $m_1=1 \; kg$, $m_2=1.5 \; kg$, $m_3=1.5 \; kg$, $m_4=1 \; kg$; the polar moments of inertia are assumed to be ideal, with uniform distribution of mass; the two springs have stiffness coefficients of $k_1=k_2=100 \; N/m$ and natural lengths of $L_{01}=2.2360 \; m$ and $L_{02}=2.0615 \; m$.

The state vector includes the natural coordinates of the point 1, 2, and 3 of the mechanism with ${q}= \big[  {q}^{\rm T}_{1} ~ {q}^{\rm T}_{2} ~ {q}^{\rm T}_{3} \big]^{\rm T}$ where  ${q}_{2}= \big[ x_{2} ~ y_{2}  \big]^{\rm T}$ represents the independent coordinates, and ${q}_{1}= \big[  x_{1} ~ y_{1}  \big]^{\rm T}$ with  ${q}_{3}= \big[ x_{3} ~ y_{3} \big]^{\rm T}$ represent the dependent coordinates. This implies that the coordinates of point 2 are the DOF of the system on which the  coordinates of points 1 and 3 are dependent. The dependent coordinates are solved using the constraints of the system. The constraints are defined according to the fixed lengths between each set of points, as follows:
\begin{align}
	\label{eq:Constraints 5 bar}
    {{ \Phi}} &= \begin{bmatrix}
           \|{q_A-q_1}\|^2  -L_{A1}^2\\
			\|{q_2-q_1}\|^2  -L_{21}^2\\
           \|{q_3-q_2}\|^2  -L_{32}^2\\
			\|{q_B-q_3}\|^2  -L_{B3}^2\\
         \end{bmatrix}=0,
\end{align}
with the lengths $L_{A1}=L_{B3}=1.4142 \; m$ ; $L_{A1}=L_{B3}=1.8027 \; m $ and the ground points ${q}_{A}= \big[ \begin{array}{c c} -0.5 & 0 \end{array} \big]^{\rm T}$; ${q}_{B}= \big[ \begin{array}{c c} 0.5 & 0 \end{array} \big]^{\rm T}$. Finally, the point 2 of the five-bar mechanism hits the ground at -2.35 m along the vertical y axis. The event (touching the ground) is simulated through the event function $r(\cdot)$ described in Eq.~\eqref{eq:event_function}. Once the event is detected, the vertical velocity of point 2 jumps to its opposite value, while its horizontal velocity remains the same.

The system of equation \eqref{eq:MBD canonical ODE system} is simulated with a time span of five seconds. Fig.~\ref{fig2} shows the residuals of the constraint equations. The position constraints are satisfied within an error of $10^{-6}$, while the velocity constraints within and error of $10^{-5}$, which is satisfactory.
\begin{figure} [H]
\centering
	\begin{subfigure}[t]{.49\textwidth}
	\centering
	\raisebox{-4.3mm}{\includegraphics[width=0.9\textwidth]{Figures/five-bar.pdf}}
	\caption{Diagram of the five-bar mechanism}
	\label{Fig_fivebar}
	\end{subfigure}
	%
	\begin{subfigure}[t]{.49\textwidth}
	\centering
	\includegraphics[width=\textwidth]{Figures/Constraints.pdf}
	\captionsetup{margin=1cm}
	\caption{The position and the velocity constraint residuals for the five-bar mechanism}
	\label{fig2}
	\end{subfigure}
\caption{Structure of the five-bar mechanism}
\label{sfig1}
\end{figure}
The trajectories of the position and velocity of point 2 of the five-bar mechanism along the vertical y axis are shown in Fig.~\ref{fig3} and Fig.~\ref{fig4}, respectively. These results show that  point 2's vertical position bounces at -2.35m, as expected, and its vertical velocity jumps at each time of event.
\begin{figure} [H]
\centering
	\begin{subfigure}{.49\textwidth}
	\centering
	\includegraphics[width=\textwidth]{Figures/PositionState.eps}
	\captionsetup{margin=1cm}
	\caption{The vertical position of the bottom point $ y_2$ of the five-bar mechanism}
	\label{fig3}
	\end{subfigure}
	%
	\begin{subfigure}{.49\textwidth}
	\centering
	\includegraphics[width=\textwidth]{Figures/VelocityState.eps}
	\captionsetup{margin=1cm}
	\caption{The vertical velocity of the bottom point $\dot y_2$ of the five-bar mechanism}
	\label{fig4}
	\end{subfigure}
\caption{Trajectories of the bottom point the five-bar mechanism}
\label{sfig2}
\end{figure}

The trajectories of the sensitivity of the position and velocity of point 2 of the five-bar mechanism along the vertical y axis are shown in Fig.~\ref{fig5} and Fig.~\ref{fig6}, respectively. The analytical sensitivity is represented by the continuous line, while the central finite difference sensitivity is represented by the dashed line. There is an excellent correlation between the numerical and the analytical sensitivities, with a difference between the two trajectories of less than $0.1\%$. Note that the numerical sensitivity of the velocity of point 2 along the vertical axis tends to be really large in magnitude, $1/\varepsilon$ at each time of event. This is shown by the vertical dashed lines and it is due to the fact that the difference between the trajectories $v({\rho}+\varepsilon)$ and $v({\rho}-\varepsilon)$ increases considerably during $\Delta t$, as shown in Fig.~\ref{fig:im2}. This result shows that the novel analytical sensitivity method presented in this paper is considerably more robust than the numerical method, as it correctly calculates the sensitivity jumps and accurately determines the sensitivity trajectories after each invent without any delta-like jumps.
\begin{figure} [H]
\centering
	\begin{subfigure}{.49\textwidth}
	\centering
	\includegraphics[width=\textwidth]{Figures/SensitivityPositionState.eps}
	\captionsetup{margin=1.1cm}
	\caption{The sensitivity of the position of the bottom point $y_2$ of the five-bar mechanism}
	\label{fig5}
	\end{subfigure}
	%
	\begin{subfigure}{.49\textwidth}
	\centering
	\includegraphics[width=\textwidth]{Figures/SensitivityVelocityState.eps}
	\captionsetup{margin=1.1cm}
	\caption{The sensitivity of the velocity of the bottom point $\dot y_2$ of the five-bar mechanism}
	\label{fig6}
	\end{subfigure}
\caption{Sensitivity analysis of the bottom point the five-bar mechanism}
\label{sfig3}
\end{figure}

The trajectories of the quadrature variables  $\bz = \int_{t_0}^{t_F} \dot y_2 \ {\rm dt}$  and $\bz = \int_{t_0}^{t_F} \ddot y_2 \ {\rm dt}$ are shown in Fig.~\ref{fig7} and Fig.~\ref{fig8}, respectively. Note that $\bz = \int_{t_0}^{t_F} \dot y_2 \ {\rm dt}$  matches the trajectory of the position of point 2 along the vertical axis in Fig.~\ref{fig3}, while $\bz = \int_{t_0}^{t_F} \ddot y_2 \ {\rm dt}$ does not completely match the trajectory of the velocity of point 2 in Fig.~\ref{fig4}. This is due to the fact that the velocity variable is affected by the impulse function at the time of event, while the quadrature variable is not. The trajectories after each invent differ by a constant since the quadrature variable evaluates the integral of the acceleration of point 2.
\begin{figure} [H]
\centering
	\begin{subfigure}{.49\textwidth}
	\centering
	\includegraphics[width=\textwidth]{Figures/Quadrature2.eps}
	\captionsetup{margin=1cm}
	\caption{The quadrature variable $\bz = \int_{t_0}^{t_F} \dot y_2 \ {\rm dt}$ }
	\label{fig7}
	\end{subfigure}
	%
	\begin{subfigure}{.49\textwidth}
	\centering
	\includegraphics[width=\textwidth]{Figures/Quadrature.eps}
	\captionsetup{margin=1cm}
	\caption{The quadrature variable $\bz = \int_{t_0}^{t_F} \ddot y_2 \ {\rm dt}$ }
	\label{fig8}
	\end{subfigure}
\caption{The quadrature variables of the five-bar mechanism.}
\label{sfig4}
\end{figure}

The trajectories of the sensitivities of the quadrature variables  $\bz = \int_{t_0}^{t_F} \dot y_2 \ {\rm dt}$ and $\bz = \int_{t_0}^{t_F} \ddot y_2 \ {\rm dt}$ are presented in  Fig.~\ref{fig9} and Fig.~\ref{fig10}, respectively. The sensitivities are with respect to the system parameters ${\rho}= \big[L_{01} ~~ L_{02} \big]$. The results of the analytical and numerical sensitivity analysis of the five-bar mechanism  highlight the quasi-perfect correlation between the numerical and analytical sensitivities with a difference between the two trajectories of less than $0.1\%$. Note that a similar observation with the one previously made is valid here: the sensitivity of $\bz = \int_{t_0}^{t_F} \dot y_2 \ {\rm dt}$ shown in Fig.~\ref{fig9} matches the trajectory of the sensitivity of the position illustrated in Fig.~\ref{fig5}, while the sensitivity of $\bz = \int_{t_0}^{t_F} \ddot y_2 \ {\rm dt}$ shown in Fig.~\ref{fig10} does not completely match the trajectory of the sensitivity of the velocity from Fig.~\ref{fig6} because of the impulse function.
\begin{figure} [H]
\centering
	\begin{subfigure}{.49\textwidth}
	\centering
	\includegraphics[width=\textwidth]{Figures/QuadratureSensitivityZ.eps}
	\captionsetup{margin=1cm}
	\caption{The sensitivity of $\bz = \int_{t_0}^{t_F} \dot y_2 \ {\rm dt}$ with ${ Z}=\frac{d{z}}{d{\rho}}$.}
	\label{fig9}
	\end{subfigure}
	%
	\begin{subfigure}{.49\textwidth}
	\centering
	\includegraphics[width=\textwidth]{Figures/QuadratureSensitivityZ2.eps}
	\captionsetup{margin=1cm}
	\caption{The sensitivity of $\bz = \int_{t_0}^{t_F} \ddot y_2 \ {\rm dt}$ with ${ Z}=\frac{d{z}}{d{\rho}}$}
	\label{fig10}
	\end{subfigure}
\caption{Sensitivity analysis of the quadrature variables of the five-bar mechanism}
\label{sfig5}
\end{figure}

\section{Conclusions}
\label{sec:conclusions}

An important ingredient in the analysis and design of multibody dynamical systems is the ability to compute sensitivities of model states with respect to system parameters. Direct and adjoint sensitivity analyses for mechanical systems with smooth trajectories have been studied in the literature \cite{Sandu_2013_sensitivity_ODE_multibody,Sandu_2014_sensitivity_ODE_multibody,zhu2014mbsvt,Zhu_2014_PhD,Sandu2015dynamic,Sandu_2017_vehicle-optimization}. This work focuses on the case where trajectories are non-smooth, and specifically, on hybrid systems where a finite number of events can instantly change the velocity state variables.

We develop herein a complete mathematical framework needed to compute the sensitivities of model states, and of general cost functionals, with respect to model parameters for hybrid dynamical systems. The paper first considers multibody system is described by ordinary differential equations, as is the case when the dynamics is modeled using recursive coordinate systems. Next, the paper multibody systems described by differential algebraic equations, as is the case when the constraints imposed by the joints holding the system together are explicitly modeled. At the time moments when events happen, and there are points of discontinuity in the forward trajectory, the evolution of sensitivity variables is also discontinuous. Jump conditions for sensitivity variables are established for both the ordinary and the differential algebraic cases. These jump conditions specify the values of the sensitivities right after the event, given their value right before the event, and the characteristics of the impact.

A five-bar mechanism with non-smooth contacts is used as a case study. The analytical sensitivities obtained by the proposed methodology are validated against numerical sensitivities computed by real and complex finite differences.

Current  efforts  focus on  applying the  hybrid  system sensitivity analysis methodology to robotics, where sensitivities of the performance of a robotic system with respect to changes in the system configuration, under non-smooth impact conditions, are a topic of great interest. Future work will consider extending the current framework to perform adjoint sensitivity analysis of hybrid mechanical systems.

\clearpage

\section*{Acknowledgments}
\label{sec:acknowledgments}

This work was supported in part by awards NSF DMS--1419003, NSF CCF--1613905, AFOSR DDDAS 15RT1037,
by the Computational Science Laboratory, and by the Advanced Vehicle Systems Laboratory at Virginia Tech.

\section{Introduction}

\file{elsarticle.cls} is a thoroughly re-written document class
for formatting \LaTeX{} submissions to Elsevier journals.
The class uses the environments and commands defined in \LaTeX{} kernel
without any change in the signature so that clashes with other
contributed \LaTeX{} packages such as \file{hyperref.sty},
\file{preview-latex.sty}, etc., will be minimal.
\file{elsarticle.cls} is primarily built upon the default
\file{article.cls}.  This class depends on the following packages
for its proper functioning:

\begin{enumerate}
\item \file{pifont.sty} for openstar in the title footnotes;
\item \file{natbib.sty} for citation processing;
\item \file{geometry.sty} for margin settings;
\item \file{fleqn.clo} for left aligned equations;
\item \file{graphicx.sty} for graphics inclusion;
\item \file{txfonts.sty} optional font package, if the document is to
  be formatted with Times and compatible math fonts;
\item \file{hyperref.sty} optional packages if hyperlinking is
  required in the document.

\end{enumerate}

All the above packages are part of any standard \LaTeX{} installation.
Therefore, the users need not be bothered about downloading any
extra packages.  Furthermore, users are free to make use of \textsc{ams}
math packages such as \file{amsmath.sty}, \file{amsthm.sty},
\file{amssymb.sty}, \file{amsfonts.sty}, etc., if they want to.  All
these packages work in tandem with \file{elsarticle.cls} without
any problems.

\section{Major Differences}

Following are the major differences between \file{elsarticle.cls}
and its predecessor package, \file{elsart.cls}:

\begin{enumerate}[\textbullet]
\item \file{elsarticle.cls} is built upon \file{article.cls}
while \file{elsart.cls} is not. \file{elsart.cls} redefines
many of the commands in the \LaTeX{} classes/kernel, which can
possibly cause surprising clashes with other contributed
\LaTeX{} packages;

\item provides preprint document formatting by default, and
optionally formats the document as per the final
style of models $1+$, $3+$ and $5+$ of Elsevier journals;

\item some easier ways for formatting \verb+list+ and
\verb+theorem+ environments are provided while people can still
use \file{amsthm.sty} package;

\item \file{natbib.sty} is the main citation processing package
  which can comprehensively handle all kinds of citations and
works perfectly with \file{hyperref.sty} in combination with
\file{hypernat.sty};

\item long title pages are processed correctly in preprint and
  final formats.

\end{enumerate}

\section{Installation}

The package is available at author resources page at Elsevier
(\url{http://www.elsevier.com/locate/latex}).
It can also be found in any of the nodes of the Comprehensive
\TeX{} Archive Network (\textsc{ctan}), one of the primary nodes
being
\url{http://www.ctan.org/tex-archive/macros/latex/contrib/elsevier/}.
Please download the \file{elsarticle.dtx} which is a composite
class with documentation and \file{elsarticle.ins} which is the
\LaTeX{} installer file. When we compile the
\file{elsarticle.ins} with \LaTeX{} it provides the class file,
\file{elsarticle.cls} by
stripping off all the documentation from the \verb+*.dtx+ file.
The class may be moved or copied to a place, usually,
\verb+$TEXMF/tex/latex/elsevier/+,
or a folder which will be read
by \LaTeX{} during document compilation.  The \TeX{} file
database needs updation after moving/copying class file.  Usually,
we use commands like \verb+mktexlsr+ or \verb+texhash+ depending
upon the distribution and operating system.

\section{Usage}\label{sec:usage}
The class should be loaded with the command:

\begin{vquote}
 \documentclass[<options>]{elsarticle}
\end{vquote}

\noindent where the \verb+options+ can be the following:

\begin{description}

\item [{\tt\color{verbcolor} preprint}]  default option which format the
  document for submission to Elsevier journals.

\item [{\tt\color{verbcolor} review}]  similar to the \verb+preprint+ option, but
  increases the baselineskip to facilitate easier review process.

\item [{\tt\color{verbcolor} 1p}]  formats the article to the look and feel of the final
  format of model 1+ journals. This is always single column style.

\item [{\tt\color{verbcolor} 3p}] formats the article to the look and feel of the final
  format of model 3+ journals. If the journal is a two column
model, use \verb+twocolumn+ option in combination.

\item [{\tt\color{verbcolor} 5p}] formats for model 5+ journals. This is always
  of two column style.

\item [{\tt\color{verbcolor} authoryear}] author-year citation style of
  \file{natbib.sty}. If you want to add extra options of
\file{natbib.sty}, you may use the options as comma delimited
strings as arguments to \verb+\biboptions+ command. An example
would be:
\end{description}
\begin{vquote}
 \biboptions{longnamesfirst,angle,semicolon}
\end{vquote}

\begin{description}
\item [{\tt\color{verbcolor} number}] numbered citation style. Extra options
  can be loaded with\linebreak \verb+\biboptions+ command.

\item [{\tt\color{verbcolor} sort\&compress}] sorts and compresses the
numbered citations. For example, citation [1,2,3] will become [1--3].

\item [{\tt\color{verbcolor} longtitle}] if front matter is unusually long, use
  this option to split the title page across pages with the correct
placement of title and author footnotes in the first page.

\item [{\tt\color{verbcolor} times}] loads \file{txfonts.sty}, if available in
  the system to use Times and compatible math fonts.

\item[] All options of \file{article.cls} can be used with this
  document class.

\item[] The default options loaded are \verb+a4paper+, \verb+10pt+,
  \verb+oneside+, \verb+onecolumn+ and \verb+preprint+.

\end{description}

\section{Frontmatter}

There are two types of frontmatter coding:
\begin{enumerate}[(1)]
\item each author is
connected to an affiliation with a footnote marker; hence all
authors are grouped together and affiliations follow;
\item authors of same affiliations are grouped together and the
relevant affiliation follows this group. An example coding of the first
type is provided below.
\end{enumerate}

\begin{vquote}
 \title{This is a specimen title\tnoteref{t1,t2}}
 \tnotetext[t1]{This document is a collaborative effort.}
 \tnotetext[t2]{The second title footnote which is a longer
    longer than the first one and with an intention to fill
    in up more than one line while formatting.}
\end{vquote}

\begin{vquote}
 \author[rvt]{C.V.~Radhakrishnan\corref{cor1}\fnref{fn1}}
 \ead{cvr@river-valley.com}

 \author[rvt,focal]{K.~Bazargan\fnref{fn2}}
 \ead{kaveh@river-valley.com}

 \author[els]{S.~Pepping\corref{cor2}\fnref{fn1,fn3}}
 \ead[url]{http://www.elsevier.com}
\end{vquote}

\begin{vquote}
 \cortext[cor1]{Corresponding author}
 \cortext[cor2]{Principal corresponding author}
 \fntext[fn1]{This is the specimen author footnote.}
 \fntext[fn2]{Another author footnote, but a little more
             longer.}
 \fntext[fn3]{Yet another author footnote. Indeed, you can have
    any number of author footnotes.}

 \address[rvt]{River Valley Technologies, SJP Building,
    Cotton Hills, Trivandrum, Kerala, India 695014}
 \address[focal]{River Valley Technologies, 9, Browns Court,
    Kennford, Exeter, United Kingdom}
 \address[els]{Central Application Management,
    Elsevier, Radarweg 29, 1043 NX\\
    Amsterdam, Netherlands}

\end{vquote}

The output of the above TeX source is given in Clips~\ref{clip1} and
\ref{clip2}. The header portion or title area is given in Clip~\ref{clip1} and
the footer area is given in Clip~\ref{clip2}.

\vspace*{6pt}
\deforange{blue!70}
\src{Header of the title page.}
\includeclip{1}{132 571 481 690}{els1.pdf}
\deforange{orange}

\deforange{blue!70}
\src{Footer of the title page.}
\includeclip{1}{122 129 481 237}{els1.pdf}
\deforange{orange}
\pagebreak

Most of the commands such as \verb+\title+, \verb+\author+,
\verb+\address+ are self explanatory.  Various components are
linked to each other by a label--reference mechanism; for
instance, title footnote is linked to the title with a footnote
mark generated by referring to the \verb+\label+ string of
the \verb=\tnotetext=.  We have used similar commands
such as \verb=\tnoteref= (to link title note to title);
\verb=\corref= (to link corresponding author text to
corresponding author); \verb=\fnref= (to link footnote text to
the relevant author names).  \TeX{} needs two compilations to
resolve the footnote marks in the preamble part.
Given below are the syntax of various note marks and note texts.

\begin{vquote}
  \tnoteref{<label(s)>}
  \corref{<label(s)>}
  \fnref{<label(s)>}
  \tnotetext[<label>]{<title note text>}
  \cortext[<label>]{<corresponding author note text>}
  \fntext[<label>]{<author footnote text>}
\end{vquote}

\noindent where \verb=<label(s)>= can be either one or more comma
delimited label strings. The optional arguments to the
\verb=\author= command holds the ref label(s) of the address(es)
to which the author is affiliated while each \verb=\address=
command can have an optional argument of a label. In the same
manner, \verb=\tnotetext=, \verb=\fntext=, \verb=\cortext= will
have optional arguments as their respective labels and note text
as their mandatory argument.

The following example code provides the markup of the second type
of author-affiliation.

\begin{vquote}
\author{C.V.~Radhakrishnan\corref{cor1}\fnref{fn1}}
 \ead{cvr@river-valley.com}
 \address{River Valley Technologies, SJP Building,
   Cotton Hills, Trivandrum, Kerala, India 695014}
\end{vquote}

\begin{vquote}
\author{K.~Bazargan\fnref{fn2}}
 \ead{kaveh@river-valley.com}
 \address{River Valley Technologies, 9, Browns Court, Kennford,
   Exeter, UK.}
\end{vquote}

\begin{vquote}
\author{S.~Pepping\fnref{fn1,fn3}}
 \ead[url]{http://www.elsevier.com}
 \address{Central Application Management,
   Elsevier, Radarweg 43, 1043 NX Amsterdam, Netherlands}
\end{vquote}

\begin{vquote}
\cortext[cor1]{Corresponding author}
\fntext[fn1]{This is the first author footnote.}
\fntext[fn2]{Another author footnote, this is a very long
  footnote and it should be a really long footnote. But this
  footnote is not yet sufficiently long enough to make two lines
  of footnote text.}
\fntext[fn3]{Yet another author footnote.}
\end{vquote}

The output of the above TeX source is given in Clip~\ref{clip3}.

\vspace*{12pt}
\deforange{blue!70}
\src{Header of the title page..}
\includeclip{1}{132 491 481 690}{els2.pdf}
\deforange{orange}

The frontmatter part has further environments such as abstracts and
keywords.  These can be marked up in the following
manner:

\begin{vquote}
 \begin{abstract}
  In this work we demonstrate the formation of a new type of
  polariton on the interface between a ....
 \end{abstract}
\end{vquote}

\begin{vquote}
 \begin{keyword}
  quadruple exiton \sep polariton \sep WGM

  \PACS 71.35.-y \sep 71.35.Lk \sep 71.36.+c
 \end{keyword}
\end{vquote}

\noindent Each keyword shall be separated by a \verb+\sep+ command.
\textsc{pacs} and \textsc{msc} classifications shall be provided in
the keyword environment with the commands \verb+\PACS+ and
\verb+\MSC+ respectively.  \verb+\MSC+ accepts an optional
argument to accommodate future revisions.
eg., \verb=\MSC[2008]=. The default is 2000.\looseness=-1

\section{Floats}
{Figures} may be included using the command, \verb+\includegraphics+ in
combination with or without its several options to further control
graphic. \verb+\includegraphics+ is provided by \file{graphic[s,x].sty}
which is part of any standard \LaTeX{} distribution.
\file{graphicx.sty} is loaded by default. \LaTeX{} accepts figures in
the postscript format while pdf\LaTeX{} accepts \file{*.pdf},
\file{*.mps} (metapost), \file{*.jpg} and \file{*.png} formats.
pdf\LaTeX{} does not accept graphic files in the postscript format.

The \verb+table+ environment is handy for marking up tabular
material. If users want to use \file{multirow.sty},
\file{array.sty}, etc., to fine control/enhance the tables, they
are welcome to load any package of their choice and
\file{elsarticle.cls} will work in combination with all loaded
packages.

\section[Theorem and ...]{Theorem and theorem like environments}

\file{elsarticle.cls} provides a few shortcuts to format theorems and
theorem-like environments with ease. In all commands the options that
are used with the \verb+\newtheorem+ command will work exactly in the same
manner. \file{elsarticle.cls} provides three commands to format theorem or
theorem-like environments:

\begin{vquote}
 \newtheorem{thm}{Theorem}
 \newtheorem{lem}[thm]{Lemma}
 \newdefinition{rmk}{Remark}
 \newproof{pf}{Proof}
 \newproof{pot}{Proof of Theorem \ref{thm2}}
\end{vquote}

The \verb+\newtheorem+ command formats a
theorem in \LaTeX's default style with italicized font, bold font
for theorem heading and theorem number at the right hand side of the
theorem heading.  It also optionally accepts an argument which
will be printed as an extra heading in parentheses.

\begin{vquote}
  \begin{thm}
   For system (8), consensus can be achieved with $\|T_{\omega z}$
   ...
     \begin{eqnarray}\label{10}
     ....
     \end{eqnarray}
  \end{thm}
\end{vquote}

Clip~\ref{clip4} will show you how some text enclosed between the
above code looks like:

\vspace*{6pt}
\deforange{blue!70}
\src{{\ttfamily\color{verbcolor}\expandafter\@gobble\string\\ newtheorem}}
\includeclip{2}{1 1 453 120}{jfigs.pdf}
\deforange{orange}

The \verb+\newdefinition+ command is the same in
all respects as its\linebreak \verb+\newtheorem+ counterpart except that
the font shape is roman instead of italic.  Both
\verb+\newdefinition+ and \verb+\newtheorem+ commands
automatically define counters for the environments defined.

\vspace*{12pt}
\deforange{blue!70}
\src{{\ttfamily\color{verbcolor}\expandafter\@gobble\string\\ newdefinition}}
\includeclip{1}{1 1 453 105}{jfigs.pdf}
\deforange{orange}

The \verb+\newproof+ command defines proof environments with
upright font shape.  No counters are defined.

\vspace*{6pt}
\deforange{blue!70}
\src{{\ttfamily\color{verbcolor}\expandafter\@gobble\string\\ newproof}}
\includeclip{3}{1 1 453 65}{jfigs.pdf}
\deforange{orange}

Users can also make use of \verb+amsthm.sty+ which will override
all the default definitions described above.

\section[Enumerated ...]{Enumerated and Itemized Lists}
\file{elsarticle.cls} provides an extended list processing macros
which makes the usage a bit more user friendly than the default
\LaTeX{} list macros.   With an optional argument to the
\verb+\begin{enumerate}+ command, you can change the list counter
type and its attributes.

\begin{vquote}
 \begin{enumerate}[1.]
 \item The enumerate environment starts with an optional
   argument `1.', so that the item counter will be suffixed
   by a period.
 \item You can use `a)' for alphabetical counter and '(i)' for
   roman counter.
  \begin{enumerate}[a)]
    \item Another level of list with alphabetical counter.
    \item One more item before we start another.
    \begin{enumerate}[(i)]
     \item This item has roman numeral counter.
     \item Another one before we close the third level.
    \end{enumerate}
    \item Third item in second level.
  \end{enumerate}
 \item All list items conclude with this step.
 \end{enumerate}
\end{vquote}

\vspace*{12pt}
\deforange{blue!70}
\src{List -- Enumerate}
\includeclip{4}{1 1 453 185}{jfigs.pdf}
\deforange{orange}

Further, the enhanced list environment allows one to prefix a
string like `step' to all the item numbers.  Take a look at the
example below:

\begin{vquote}
 \begin{enumerate}[Step 1.]
  \item This is the first step of the example list.
  \item Obviously this is the second step.
  \item The final step to wind up this example.
 \end{enumerate}
\end{vquote}

\deforange{blue!70}
\src{List -- enhanced}
\includeclip{5}{1 1 313 83}{jfigs.pdf}
\deforange{orange}

\vspace*{-18pt}

\section{Cross-references}
In electronic publications, articles may be internally
hyperlinked. Hyperlinks are generated from proper
cross-references in the article.  For example, the words
\textcolor{black!80}{Fig.~1} will never be more than simple text,
whereas the proper cross-reference \verb+\ref{tiger}+ may be
turned into a hyperlink to the figure itself:
\textcolor{blue}{Fig.~1}.  In the same way,
the words \textcolor{blue}{Ref.~[1]} will fail to turn into a
hyperlink; the proper cross-reference is \verb+\cite{Knuth96}+.
Cross-referencing is possible in \LaTeX{} for sections,
subsections, formulae, figures, tables, and literature
references.

\section[Mathematical ...]{Mathematical symbols and formulae}

Many physical/mathematical sciences authors require more
mathematical symbols than the few that are provided in standard
\LaTeX. A useful package for additional symbols is the
\file{amssymb} package, developed by the American Mathematical
Society. This package includes such oft-used symbols as
$\lesssim$ (\verb+\lesssim+), $\gtrsim$ (\verb+\gtrsim+)  or
$\hbar$ (\verb+\hbar+). Note that your \TeX{}
system should have the \file{msam} and \file{msbm} fonts installed. If
you need only a few symbols, such as $\Box$ (\verb+\Box+), you might try the
package \file{latexsym}.

Another point which would require authors' attention is the
breaking up of long equations.  When you use
\file{elsarticle.cls} for formatting your submissions in the
\verb+preprint+ mode, the document is formatted in single column
style with a text width of 384pt or 5.3in.  When this document is
formatted for final print and if the journal happens to be a double column
journal, the text width will be reduced to 224pt at for 3+
double column and 5+ journals respectively. All the nifty
fine-tuning in equation breaking done by the author goes to waste in
such cases.  Therefore, authors are requested to check this
problem by typesetting their submissions in final format as well
just to see if their equations are broken at appropriate places,
by changing appropriate options in the document class loading
command, which is explained in section~\ref{sec:usage},
\nameref{sec:usage}. This allows authors to fix any equation breaking
problem before submission for publication.
\file{elsarticle.cls} supports formatting the author submission
in different types of final format.  This is further discussed in
section \ref{sec:final}, \nameref{sec:final}.

\section{Bibliography}

Three bibliographic style files (\verb+*.bst+) are provided ---
\file{elsarticle-num.bst}, \file{elsarticle-num-names.bst} and
\file{elsarticle-harv.bst} --- the first one for the numbered scheme, the
second for the numbered with new options of \file{natbib.sty} and the
last one for the author year scheme.

In \LaTeX{} literature, references are listed in the
\verb+thebibliography+ environment.  Each reference is a
\verb+\bibitem+ and each \verb+\bibitem+ is identified by a label,
by which it can be cited in the text:

\verb+\bibitem[Elson et al.(1996)]{ESG96}+ is cited as
\verb+\citet{ESG96}+.

\noindent In connection with cross-referencing and
possible future hyperlinking it is not a good idea to collect
more that one literature item in one \verb+\bibitem+.  The
so-called Harvard or author-year style of referencing is enabled
by the \LaTeX{} package \file{natbib}. With this package the
literature can be cited as follows:

\begin{enumerate}[\textbullet]
\item Parenthetical: \verb+\citep{WB96}+ produces (Wettig \& Brown, 1996).
\item Textual: \verb+\citet{ESG96}+ produces Elson et al. (1996).
\item An affix and part of a reference:
\verb+\citep[e.g.][Ch. 2]{Gea97}+ produces (e.g. Governato et
al., 1997, Ch. 2).
\end{enumerate}

In the numbered scheme of citation, \verb+\cite{<label>}+ is used,
since \verb+\citep+ or \verb+\citet+ has no relevance in the numbered
scheme.  \file{natbib} package is loaded by \file{elsarticle} with
\verb+numbers+ as default option.  You can change this to author-year
or harvard scheme by adding option \verb+authoryear+ in the class
loading command.  If you want to use more options of the \file{natbib}
package, you can do so with the \verb+\biboptions+ command, which is
described in the section \ref{sec:usage}, \nameref{sec:usage}.  For
details of various options of the \file{natbib} package, please take a
look at the \file{natbib} documentation, which is part of any standard
\LaTeX{} installation.

\subsection*{Displayed equations and double column journals}

Many Elsevier journals print their text in two columns. Since
the preprint layout uses a larger line width than such columns,
the formulae are too wide for the line width in print. Here is an
example of an equation  (see equation 6) which is perfect in a
single column preprint format:

\bigskip
\setlength\Sep{6pt}
\src{See equation (6)}
\deforange{blue!70}
\includeclip{4}{134 391 483 584}{els1.pdf}
\deforange{orange}

\noindent When this document is typeset for publication in a
model 3+ journal with double columns, the equation will overlap
the second column text matter if the equation is not broken at
the appropriate location.

\vspace*{6pt}
\deforange{blue!70}
\src{See equation (6) overprints into second column}
\includeclip{3}{61 531 532 734}{els-3pd.pdf}
\deforange{orange}

\pagebreak

\noindent The typesetter will try to break the equation which
need not necessarily be to the liking of the author or as it
happens, typesetter's break point may be semantically incorrect.
Therefore, authors may check their submissions for the incidence
of such long equations and break the equations at the correct
places so that the final typeset copy will be as they wish.

\section{Final print}\label{sec:final}

The authors can format their submission to the page size and margins
of their preferred journal.  \file{elsarticle} provides four
class options for the same. But it does not mean that using these
options you can emulate the exact page layout of the final print copy.
\lmrgn=3em
\begin{description}
\item [\texttt{1p}:] $1+$ journals with a text area of
384pt $\times$ 562pt or 13.5cm $\times$ 19.75cm or 5.3in $\times$
7.78in, single column style only.

\item [\texttt{3p}:] $3+$ journals with a text area of 468pt
$\times$ 622pt or 16.45cm $\times$ 21.9cm or 6.5in $\times$
8.6in, single column style.

\item [\texttt{twocolumn}:] should be used along with 3p option if the
journal is $3+$ with the same text area as above, but double column
style.

\item [\texttt{5p}:] $5+$ with text area of 522pt $\times$
682pt or 18.35cm $\times$ 24cm or 7.22in $\times$ 9.45in,
double column style only.
\end{description}

Following pages have the clippings of different parts of
the title page of different journal models typeset in final
format.

Model $1+$ and $3+$  will have the same look and
feel in the typeset copy when presented in this document. That is
also the case with the double column $3+$ and $5+$ journal article
pages. The only difference will be wider text width of
higher models.  Therefore we will look at the
different portions of a typical single column journal page and
that of a double column article in the final format.

\vspace*{2pc}

\begin{center}
\hypertarget{bsc}{}
\hyperlink{sc}{
{\bf [Specimen single column article -- Click here]}
}

\vspace*{2pc}

\hypertarget{bsc}{}
\hyperlink{dc}{
{\bf [Specimen double column article -- Click here]}
}
\end{center}

\newpage
\vspace*{-2pc}
\src{}\hypertarget{sc}{}
\deforange{blue!70}
\hyperlink{bsc}{\includeclip{1}{121 81 497 670}{els1.pdf}}
\deforange{orange}

\newpage

\src{}\hypertarget{dc}{}
\deforange{blue!70}
\hyperlink{bsc}{\includeclip{1}{55 93 535 738}{els-3pd.pdf}}
\deforange{orange}

\end{document}

\section{}
\label{}

\section{}
\label{}

redpajama_set_name: RedPajamaArXiv