In this chapter, we give some hints on the automatic computation of the number of index sets, the number of derivatives in the {\tt Interaction} and the levelMin and LevelMax.


\section{Why is the  relative degree not relevant ?}
In this section, we give a very brief overview of the notion of relative degree.

\subsection{First order Linear complementary systems}

 A  Linear Complementarity System (LCS) is defined by
\begin{equation}
  \label{eq:LCS-bis}
  \begin{cases}
    \dot x = A x +B \lambda \\
     y = C x + D \lambda\\
    0 \leq  y \perp \lambda \geq 0 \\
  \end{cases}
\end{equation} 
 \begin{definition}[Relative degree in the SISO case]
      Let us consider a linear system in state representation given by the quadruplet $(A,B,C,D) \in \RR^{n\times n}\times\RR^{n \times m}\times \RR^{m\times n}\times\RR^{m\times m} $:
      \begin{equation}
        \label{eq:LS}
        \begin{cases}
          \dot x = A x +B \lambda \\
          y = C x + D \lambda
        \end{cases}
      \end{equation}
      \begin{itemize}
      \item In the Single Input/ Single Output (SISO) case ($m=1$), the relative
        degree is defined by the first non zero Markov parameters :
        \begin{equation}
          \label{eq:Markov-Parameter}
          D, CB, CAB, CA^2B, \ldots, CA^{r-1}B, \ldots
        \end{equation}
      \item In the multiple input/multiple output (MIMO) case ($m>1$), an \textit{uniform} relative degree is defined as follows.
        If $D$ is non singular, the relative degree is equal to $0$. Otherwise, it is assumed to be the first positive integer $r$ such that 
        \begin{equation}
          \label{eq:mimo-r}
          CA^{i}B =0, \quad i=0\ldots q-2
        \end{equation}
        while
        \begin{equation}
          \label{eq:mimo-r2}
          CA^{r-1}B \text{ is non singular}.
        \end{equation}
      \end{itemize}
    \end{definition}
    The Markov parameters arise naturally when we derive with respect
    to time the output $y$,
    \begin{eqnarray*}
      \label{eq:y-derive}
      y &=& C x + D \lambda \\
      \dot y &=& CA x + CB \lambda, \text{ if } D= 0  \\
      \ddot y &=& CA^2 x + CAB \lambda, \text{ if }  D=0, CB=0\\
      &\ldots& \\
      y^{(r)} &=& CA^{r} x + CA^{r-1}B \lambda, \text{ if } D=0, CB=0, CA^{r-2}B=0, r=1\ldots r-2 \\
      &\ldots&
    \end{eqnarray*}
    and the first non zero Markov parameter allows us to define the
    output $y$ directly in terms of the input $\lambda$.

In continuous time, the nature of solutions depends strongly on the relative degree. When we want to perform the time--integration of such systems, we need also to reduce the relative degree or to known it to correctly operate.

We can observe that the relative degree $0$ is well defined only by the relation ($D$ nonsingular) and by the nonsmooth law. Indeed, let us imagine that the nonsmooth law is defined by $0\leq\dot y \perp \lambda \geq 0 $. We can easily see that the relative degree is reduced.

In the MIMO, the computation of non uniform relative degree is hard task. This is also the case for nonlinear systems.



\subsection{Second order Lagrangian systems}


Let us consider a second order linear and time-invariant Lagrangian dynamical system (see \S~\ref{Sec:LagrangianLineatTIDS})
\begin{equation}
  \label{eq:rd1}
  \begin{cases}
    M \dot v + C v + K q = F_{Ext}(t) + p \\
    \dot q = v
  \end{cases}
\end{equation}
together with a Lagrangian linear relation
\begin{equation}
  y= Cq + e + D \lambda + Fz,
  \label{eq:rd2}
\end{equation}
\begin{equation}
  p = C^t \lambda
\label{eq:rd3}  
\end{equation}
and a simple nonsmooth law,
\begin{equation}
  0\leq y \perp \lambda \geq 0
\label{eq:rd4}  
\end{equation}

If $D>0$, the relative degree is uniformly zero and the system can be solved without deriving the output~(\ref{eq:rd2}). Indeed, we known that the solution of the LCP
\begin{equation}
  0\leq Cq + e + D \lambda + Fz, \perp \lambda \geq 0
\label{eq:rd5}  
\end{equation}
is unique and Lipschitz with respect to $q$. It can be denoted as $\lambda(q) = \mbox{SOL}(D,Cq + e +Fz)$. Therefore, the differential equation~(\ref{eq:rd1}) reduces to a standard ODE with a Lipschitz RHS
 \begin{equation}
  \label{eq:rd6}
  \begin{cases}
    M \dot v + C v + K q = F_{Ext}(t) + C^t \lambda(q)  \\
    \dot q = v
  \end{cases}
\end{equation}

In the case that we deal with unilateral contact, we usually have $D=0$ and the relative degree of the system is $2$. In this case, the output has to be differentiated as
\begin{equation}
  \label{eq:rd7}
   \dot y= C \dot q,
\end{equation}
and an impact law has to added, for instance the newton's impact law
\begin{equation}
  \label{eq:rd8}
  \text{ if } y=0, \text{ when } \dot y^+= -e y^-
\end{equation}
In the same vein, the equations of motion (\ref{eq:rd1}) is not sufficient since the velocity may encounter jumps. The dynamics is usually replaced by a measure differential equation of the form
\begin{equation}
  \label{eq:rd10}
  \begin{cases}
    M dv + C v^+(t) dt + K q(t) dt = F_{Ext}(t)dt + di \\
    \dot q = v
  \end{cases}
\end{equation}
where $di$ is the measure that can be related to $p$ thanks to
\begin{equation}
  \label{eq:rd11}
  di = p dt + \sigma \delta _{t^*}
\end{equation}
is only one jump is expected at ${t^*}$.


\subsection{Conclusion for the implementation}
From the continuous time mathematical analysis, the relative degree is very important to know if we have to compute the derivatives of the output $y^{(n)}$ and to consider various levels for the input $p , \sigma, ....$

However in the numerical practice,  the time --discretization makes an assumption on the relative degree and treats the nonsmooth law at different levels. The resulting time discretized system posseses more or less variables.

Consider for instance  (\ref{eq:rd1}) in the case of the Moreau scheme
\begin{subnumcases}{\label{eq:MoreauTS}}
  M(v_{k+1}-v_k)  + h  (K q_{k+\theta}+ C v_{k+\theta}) = p_{k+1} = G(q_{k+1}) \lambda_{k+1},\quad\,\\[1mm] 
  q_{k+1} = q_{k} + h v_{k+\theta}, \quad \\[1mm]
  \dot y_{k+1} = G^\top(q_{k+1})\, v_{k+1} \\[1mm]
  \begin{array}{l}
    \text{if }\quad\bar y^\alpha_{k+1} \leq 0 \text{ then } 0 \leq \dot y^\alpha_{k+1} + e  \dot y^\alpha_{k} \perp \lambda^\alpha_{k+1}  \geq 0, \\[1mm]
    \text{otherwise  } \lambda^\alpha_{k+1}  =0.\label{eq:MoreauTSd}
  \end{array}, \alpha \in \mathcal I
\end{subnumcases} 
and the Schatzman--Paoli  scheme
\begin{subnumcases}{}
  M(q_{k+1}-2q_{k}+q_{k-1})  + h^2 (K q_{k+\theta}+ C v_{k+\theta})  =  p_{k+1},\quad\,\\ \notag\\ 
  v_{k+1}=\Frac{q_{k+1}-q_{k-1}}{2h}, \\ \notag \\
  y_{k+1} = h\left(\Frac{q_{k+1}+e q_{k-1}}{1+e}\right) \\
  p_{k+1}= G\left(\Frac{q_{k+1}+e q_{k-1}}{1+e}\right) \lambda_{k+1} \\
  0 \leq y_{k+1}  \perp\lambda_{k+1} \geq 0 .
\end{subnumcases}

We can see easily that the number of derivatives (or levels) that we store for $y$ and $\lambda$ is independent on the relative degree but is chosen by the {\tt OneStepIntegrator} with respect to the type of systems.

\section{How to define and compute the various levels and the number of indexSets }

\subsection{$y$ related variables}

The size of the vector {\tt y} in the {\tt Interaction} depends on
\begin{itemize}
\item the {\tt  OneStepIntegrator} type.
  \begin{itemize}
  \item see the difference between the Moreau and Schatzman Paoli
    scheme,
  \item plan the time--discontinuous Galerkin scheme
  \item plan the Higher Order Moreau sweeping process (HOSP)
  \end{itemize}
\item the {\tt  Simulation} type.
  \begin{itemize}
  \item In {\tt Timestepping} or {\tt Event-driven} we do not need the same number of stored $y$
  \end{itemize}

\item the {\tt NonSmoothLaw} type.
  \begin{itemize}
  \item If we consider some cases with or without friction in {\tt
      Timestepping} or {\tt Event-driven}, we need to adapt the number
    of stored $y$
  \end{itemize}

\end{itemize}

Since the various levels of  {\tt y} are used to build the index sets we will need from $0$ to a computed size that depends on the previous criteria. Only a part will be used in the {\tt OneStepNSProblem}.

\subsection{$\lambda$ related variables}

The size of the vector {\tt lambda} in the {\tt Interaction} depends on the same criteria than in the previous section.  Only, the number of lambda is not the same as {\tt y} since a multiplier {\tt lambda[i]} is not necessarily related to {\tt y[i]}

\section{Rules for implementation}

We can define new members in {\tt Interaction}:
\begin{itemize}
\item {\tt \_lowerlevelForOutput}, this value is to $0$ by default
\item {\tt \_upperlevelForOutput},  this value must be computed at initialization with respect to the previous criteria
\item {\tt \_lowerlevelForInput},  this value must be computed at initialization with respect to the previous criteria
\item {\tt \_upperlevelForInput},  this value must be computed at initialization with respect to the previous criteria
\end{itemize}




This level are computed in {\tt Simulation::ComputeLevelsForInputAndOutput}. A visitor is used for the {\tt OneStepIntegrator}. Furthermore, four global levels are computed 
\begin{itemize}
\item {\_levelMinForOutput} this value is the minimum level for the output {\tt Interaction::\_lowerlevelForOutput}  for all the interactions
\item {\_levelMaxForOutput} this value is the maximum level for the output {\tt Interaction::\_upperlevelForOutput}  for all the interactions
\item {\_levelMinForInput} this value is the minimum level for the output {\tt Interaction::\_lowerlevelForInput}  for all the interactions
\item {\_levelMaxForInput} this value is the maximum level for the output {\tt Interaction::\_upperlevelForInput}  for all the interactions
\end{itemize}




\begin{itemize}
\item the values {\tt y[i]} must be initialized from {\tt \_lowerlevelForOutput} to {\tt \_upperlevelForOutput}.
\item the values {\tt lamdba[i]} must be initialized from {\tt \_lowerlevelForInput} to  {\tt \_upperlevelForInput}.
\item the values {\tt y[i]} in {\tt Interaction} must be used in priority to store the i-th derivative of $y$. When it is needed, higher index $i$ can be used for other triggering variables. For instance, for an Event--Driven scheme with a Lagrangian systems with friction, sliding velocity must be stored.
\item the values of {\tt lamdba[i]} must stored the various multiplier for the nonsmooth law. We affect the same index $i$ as for the level of {\tt y[i]} present in the corresponding nonsmooth law.
\item The number of {\tt IndexSets} should follows {\tt \_levelMaxForY}.
\end{itemize}



For the dynamical systems :
\begin{itemize}
\item The number of levels for {\tt \_r} and {\tt \_p} in the DS should follow {\tt \_lowerlevelForInput} and {\tt \_upperlevelForOutput} of the associated interactions. This is done in {\tt Interaction::initialize}.
\item A new variable should be added in the LagrangianDS to store the multiplier at the position level ({\tt \_tau} ?) to avoid the use of {\tt \_p[0]}. Indeed, we will continue to assume that {\tt \_p} is the input in the equation of motion. For {\tt lambda} we can  use {\tt lambda[0]} 
\end{itemize}

TODO LIST AND QUESTIONS
\begin{itemize}
\item What about the case of multiples interactions on a DS with various {\tt \_lowerlevelForInput} and {\tt \_upperlevelForOutput} ? Normally, all the levels should be correctly initialized thanks to the proposed implementation (r2821)
\item {\tt DynamicalSystem::\_r} should be a VectorOfVectors
\item {\tt DynamicalSystem::\_r} is split in {\tt LagrangianDS}. a first part is {\tt LagrangianDS::\_p}. The other is not implemented !! {\tt LagrangianDS::\_tau} ?
\end{itemize}


%%% Local Variables: 
%%% mode: latex
%%% TeX-master: "DevNotes"
%%% End: