Section: Overall Objectives
General scope and motivations
Nonsmooth dynamics concerns the study of the time evolution of systems that are not smooth in the mathematical sense, i.e., systems that are characterized by a lack of differentiability, either of the mappings in their formulations, or of their solutions with respect to time. The class of nonsmooth dynamical systems recovers a large variety of dynamical systems that arise in many applications. The term “nonsmooth”, as the term “nonlinear”, does not precisely define the scope of the systems we are interested in but, and most importantly, they are characterized by the mathematical and numerical properties that they share. To give more insight of what are nonsmooth dynamical systems, we give in the sequel a very brief introduction of their salient features. For more details, we refer to [1], [2] [60], [77], [93], [62], [37].
A flavor of nonsmooth dynamical systems
As a first illustration, let us consider a linear finitedimensional system described by its state $x\left(t\right)\in {I\phantom{\rule{0.166667em}{0ex}}\phantom{\rule{0.166667em}{0ex}}R}^{n}$ over a timeinterval $t\in [0,T]$:
$\dot{x}\left(t\right)=Ax\left(t\right)+a,\phantom{\rule{1.em}{0ex}}A\in {I\phantom{\rule{0.166667em}{0ex}}\phantom{\rule{0.166667em}{0ex}}R}^{n\times n},\phantom{\rule{0.166667em}{0ex}}a\in {I\phantom{\rule{0.166667em}{0ex}}\phantom{\rule{0.166667em}{0ex}}R}^{n},$  (1) 
subjected to a set of $m$ inequality (unilateral) constraints:
$y\left(t\right)=Cx\left(t\right)+c\ge 0,\phantom{\rule{1.em}{0ex}}\phantom{\rule{1.em}{0ex}}C\in {I\phantom{\rule{0.166667em}{0ex}}\phantom{\rule{0.166667em}{0ex}}R}^{m\times n},\phantom{\rule{0.166667em}{0ex}}c\in {I\phantom{\rule{0.166667em}{0ex}}\phantom{\rule{0.166667em}{0ex}}R}^{m}.$  (2) 
If the constraints are physical constraints, a standard modeling approach is to augment the dynamics in (1) by an input vector $\lambda \left(t\right)\in {I\phantom{\rule{0.166667em}{0ex}}\phantom{\rule{0.166667em}{0ex}}R}^{m}$ that plays the role of a Lagrange multiplier vector. The multiplier restricts the trajectory of the system in order to respect the constraints. Furthermore, as in the continuous optimization theory, the multiplier must be signed and must vanish if the constraint is not active. This is usually formulated as a complementarity condition:
which models the onesided effect of the inequality constraints. The notation $y\ge 0$ holds component–wise and $y\perp \lambda $ means ${y}^{T}\lambda =0$. All together we end up with a Linear Complementarity System (LCS) of the form,
$\left\{\begin{array}{c}\dot{x}\left(t\right)=Ax\left(t\right)+a+B\lambda \left(t\right)\hfill \\ y\left(t\right)=Cx\left(t\right)+c\hfill \\ 0\le y\left(t\right)\perp \lambda \left(t\right)\ge 0\hfill \end{array}\right.$  (4) 
where $B\in {I\phantom{\rule{0.166667em}{0ex}}\phantom{\rule{0.166667em}{0ex}}R}^{n\times m}$ is the matrix that models the input generated by the constraints. In a more general way, the constraints may also involve the Lagrange multiplier,
$y\left(t\right)=Cx\left(t\right)+c+D\lambda \left(t\right)\ge 0,\phantom{\rule{1.em}{0ex}}\phantom{\rule{1.em}{0ex}}D\in {I\phantom{\rule{0.166667em}{0ex}}\phantom{\rule{0.166667em}{0ex}}R}^{m\times m},$  (5) 
leading to a general definition of LCS as
$\left\{\begin{array}{c}\dot{x}\left(t\right)=A\phantom{\rule{0.166667em}{0ex}}x\left(t\right)+a+B\phantom{\rule{0.166667em}{0ex}}\lambda \left(t\right)\hfill \\ y\left(t\right)=C\phantom{\rule{0.166667em}{0ex}}x\left(t\right)+c+D\phantom{\rule{0.166667em}{0ex}}\lambda \left(t\right)\hfill \\ 0\le y\left(t\right)\perp \lambda \left(t\right)\ge 0.\hfill \end{array}\right.$  (6) 
The complementarity condition, illustrated in Figure 1 is the archetype of a nonsmooth graph that we extensively use in nonsmooth dynamics. The mapping $y\mapsto \lambda $ is a multivalued (setvalued) mapping, that is nonsmooth at the origin. It has a lot of interesting mathematical properties and reformulations that come mainly from convex analysis and variational inequality theory. Let us introduce the indicator function of ${I\phantom{\rule{0.166667em}{0ex}}\phantom{\rule{0.166667em}{0ex}}R}_{+}$ as
${\Psi}_{{I\phantom{\rule{0.166667em}{0ex}}\phantom{\rule{0.166667em}{0ex}}R}_{+}}\left(x\right)=\left\{\begin{array}{c}0\phantom{\rule{4.pt}{0ex}}\text{if}\phantom{\rule{4.pt}{0ex}}x\ge 0,\hfill \\ +\infty \phantom{\rule{4.pt}{0ex}}\text{if}\phantom{\rule{4.pt}{0ex}}x<0.\hfill \end{array}\right.$  (7) 
This function is convex, proper and can be subdifferentiated [66]. The definition of the subdifferential of a convex function $f:{I\phantom{\rule{0.166667em}{0ex}}\phantom{\rule{0.166667em}{0ex}}R}^{m}\to I\phantom{\rule{0.166667em}{0ex}}\phantom{\rule{0.166667em}{0ex}}R$ is defined as:
$\partial f\left(x\right)=\{{x}^{\u2606}\in {I\phantom{\rule{0.166667em}{0ex}}\phantom{\rule{0.166667em}{0ex}}R}^{m}\mid f\left(z\right)\ge f\left(x\right)+{(zx)}^{\top}{x}^{\u2606},\forall z\}.$  (8) 
A basic result of convex analysis reads as
$0\le y\perp \lambda \ge 0\u27fa\lambda \in \partial {\Psi}_{{I\phantom{\rule{0.166667em}{0ex}}\phantom{\rule{0.166667em}{0ex}}R}_{+}}\left(y\right)$  (9) 
that gives a first functional meaning to the setvalued mapping $y\mapsto \lambda $. Another interpretation of $\partial {\Psi}_{{I\phantom{\rule{0.166667em}{0ex}}\phantom{\rule{0.166667em}{0ex}}R}_{+}}$ is based on the normal cone to a closed and nonempty convex set $C$:
${N}_{C}\left(x\right)=\{v\in {I\phantom{\rule{0.166667em}{0ex}}\phantom{\rule{0.166667em}{0ex}}R}^{m}{v}^{\top}(zx)\le 0\phantom{\rule{0.277778em}{0ex}}\text{for}\phantom{\rule{4.pt}{0ex}}\text{all}\phantom{\rule{0.277778em}{0ex}}z\in K\}.$  (10) 
It is easy to check that $\partial {\Psi}_{{I\phantom{\rule{0.166667em}{0ex}}\phantom{\rule{0.166667em}{0ex}}R}_{+}}={N}_{{I\phantom{\rule{0.166667em}{0ex}}\phantom{\rule{0.166667em}{0ex}}R}_{+}}\left(x\right)$ and it follows that
$0\le y\perp \lambda \ge 0\u27fa\lambda \in {N}_{{I\phantom{\rule{0.166667em}{0ex}}\phantom{\rule{0.166667em}{0ex}}R}_{+}}\left(y\right).$  (11) 
Finally, the definition of the normal cone yields a variational inequality:
The relations (11) and (12) allow one to formulate the complementarity system with $D=0$ as a differential inclusion based on a normal cone (see (15)) or as a differential variational inequality. By extending the definition to other types of convex functions, possibly nonsmooth, and using more general variational inequalities, the same framework applies to the nonsmooth laws depicted in Figure 2 that includes the case of piecewise smooth systems.
The mathematical concept of solutions depends strongly on the nature of the matrix quadruplet $(A,B,C,D)$ in (6). If $D$ is a positive definite matrix (or a $P$matrix), the Linear Complementarity problem
admits a unique solution $\lambda \left(x\right)$ which is a Lipschitz continuous mapping. It follows that the Ordinary Differential Equation (ODE)
is a standard ODE with a Lipschitz righthand side with a ${C}^{1}$ solution for the initial value problem. If $D=0$, the system can be written as a differential inclusion in a normal cone as
$\dot{x}\left(t\right)+Ax\left(t\right)+a\in B{N}_{{I\phantom{\rule{0.166667em}{0ex}}\phantom{\rule{0.166667em}{0ex}}R}_{+}}\left(Cx\left(t\right)\right),$  (15) 
that admits a solution that is absolutely continuous if $CB$ is a definite positive matrix and the initial condition satisfies the constraints. The time derivative $\dot{x}\left(t\right)$ and the multiplier $\lambda \left(t\right)$ may have jumps and are generally considered as functions of bounded variations. If $CB=0$, the order of nonsmoothness increases and the Lagrange multiplier may contain Dirac atoms and must be considered as a measure. Higher–order index, or higher relative degree systems yield solutions in terms of distributions and derivatives of distributions [30].
A lot of variants can be derived from the basic form of linear complementarity systems, by changing the form of the dynamics including nonlinear terms or by changing the complementarity relation by other multivalued maps. In particular the nonnegative orthant may be replaced by any convex closed cone $K\subset {I\phantom{\rule{0.166667em}{0ex}}\phantom{\rule{0.166667em}{0ex}}R}^{m}$ leading to complementarity over cones
where ${K}^{\u2606}$ its dual cone given by
${K}^{\u2606}=\{x\in {I\phantom{\rule{0.166667em}{0ex}}\phantom{\rule{0.166667em}{0ex}}R}^{m}\mid {x}^{\top}y\ge 0\phantom{\rule{4.pt}{0ex}}\text{for}\phantom{\rule{4.pt}{0ex}}\text{all}\phantom{\rule{4.pt}{0ex}}y\in K\}.$  (17) 
In Figure 2, we illustrate some other basic maps that can used for defining the relation between $\lambda $ and $y$. The saturation map, depicted in Figure 2(a) is a single valued continuous function which is an archetype of piecewise smooth map. In Figure 2(b), the relay multifunction is illustrated. If the upper and the lower limits of $\lambda $ are respectively equal to 1 and $1$, we obtain the multivalued sign function defined as
$Sgn\left(y\right)=\left\{\begin{array}{cc}1,\hfill & y>0\hfill \\ [1,1],\hfill & y=0\hfill \\ 1,\hfill & y<0.\hfill \end{array}\right.$  (18) 
Using again convex analysis, the multivalued sign function may be formulated as an inclusion into a normal cone as
More generally, any system of the type,
$\left\{\begin{array}{c}\dot{x}\left(t\right)=Ax\left(t\right)+a+B\lambda \left(t\right)\hfill \\ y\left(t\right)=Cx\left(t\right)+a\hfill \\ \lambda \left(t\right)\in Sgn\left(y\right(t\left)\right),\hfill \end{array}\right.$  (20) 
can reformulated in terms of the following setvalued system
$\left\{\begin{array}{c}\dot{x}\left(t\right)=Ax\left(t\right)+a+B\lambda \left(t\right)\hfill \\ y\left(t\right)=Cx\left(t\right)+c\hfill \\ y\left(t\right)\in {N}_{{[1,1]}^{m}}\left(\lambda \left(t\right)\right).\hfill \end{array}\right.$  (21) 
The system (21) appears in a lot of applications; among them, we can cite the sliding mode control, electrical circuits with relay and Zener diodes [26], or mechanical systems with friction [28].
Though this class of systems seems to be rather specific, it includes as well more general dynamical systems such as piecewise smooth systems and discontinuous ordinary differential equations. Indeed, the system (20) for scalars $y$ and $\lambda $ can be viewed as a discontinuous differential equation:
$\dot{x}\left(t\right)=\left\{\begin{array}{ccc}Ax+a+B\hfill & \phantom{\rule{4.pt}{0ex}}\text{if}\phantom{\rule{4.pt}{0ex}}& Cx+c>0\hfill \\ Ax+aB\hfill & \phantom{\rule{4.pt}{0ex}}\text{if}\phantom{\rule{4.pt}{0ex}}& Cx+c<0.\hfill \end{array}\right.$  (22) 
One of the most wellknown mathematical framework to deal with such systems is the Filippov theory [60] that embed the discontinuous differential equations into a differential inclusion. In the case of a single discontinuity surface given in our example by $S=\{x\mid Cx+c=0\}$, the Filippov differential inclusion based on the convex hull of the vector fields in the neighborhood of $S$ is equivalent to the use of the multivalued sign function in (20). Conversely, as it has been shown in [35], a piecewise smooth system can be formulated as a nonsmooth system based on products of multivalued sign functions.
Nonsmooth Dynamical systems in the large
Generally, the nonsmooth dynamical systems we propose to study mainly concern systems that possess the following features:

A nonsmooth formulation of the constitutive/behavioral laws that define the system. Examples of nonsmooth formulations are piecewise smooth functions, multi–valued functions, inequality constraints, yielding various definitions of dynamical systems such as piecewise smooth systems, discontinuous ordinary differential equations, complementarity systems, projected dynamical systems, evolution or differential variational inequalities and differential inclusions (into normal cones). Fundamental mathematical tools come from convex analysis [86], [67], [66], complementarity theory [55], and variational inequalities theory [59].

A concept of solutions that does not require continuously differentiable functions of time. For instance, absolutely continuous, Lipschitz continuous functions or functions of local bounded variation are the basis for solution concepts. Measures or distributions are also solutions of interest for differential inclusions or evolution variational inequalities.
Nonsmooth systems versus hybrid systems
The nonsmooth dynamical systems we are dealing with, have a nonempty intersection with hybrid systems and cyberphysical systems, as it is briefly discussed in Sect. 3.2.4. Like in hybrid systems, nonsmooth dynamical systems define continuous–time dynamics that can be identified to modes separated by guards, defined by the constraints. However, the strong mathematical structure of nonsmooth dynamical systems allows us to state results on the following points:

Mathematical concept of solutions: wellposedness (existence, and possibly, uniqueness properties, (dis)continuous dependence on initial conditions).

Dynamical systems theoretic properties: existence of invariants (equilibria, limit cycles, periodic solutions,...) and their stability, existence of oscillations, periodic and quasiperiodic solutions and propagation of waves.

Control theoretic properties: passivity, controllability, observability, stabilization, robustness.
These latter properties, that are common for smooth nonlinear dynamical systems, distinguish the nonsmooth dynamical systems from the very general definition of hybrid or cyberphysical systems [39], [65]. Indeed, it is difficult to give a precise mathematical concept of solutions for hybrid systems since the general definition of hybrid automata is usually too loose.
Numerical methods for nonsmooth dynamical systems
To conclude this brief exposition of nonsmooth dynamical systems, let us recall an important fact related to numerical methods. Beyond their intrinsic mathematical interest, and the fact that they model real physical systems, using nonsmooth dynamical systems as a model is interesting, because it exists a large set of robust and efficient numerical techniques to simulate them. Without entering into deeper details, let us give two examples of these techniques:

Numerical time integration methods: convergence, efficiency (order of consistency, stability, symplectic properties). For the nonsmooth dynamical systems described above, there exist event–capturing time–stepping schemes with strong mathematical results. These schemes have the ability to numerically integrate the initial value problem without performing an event location, but by capturing the event within a time step. We call an event, or a transition, every change into the index set of the active constraints in the complementarity formulation or in the normal cone inclusion. Hence these schemes are able to simulate systems with a huge number of transitions or even worth finite accumulation of events (Zeno behavior). Furthermore, the schemes are not suffering from the weaknesses of the standard schemes based on a regularization (smoothing) of the multivalued mapping resulting in stiff ordinary differential equations. For the time–integration of the initial value problem (IVP), or Cauchy problem, a lot of improvements of the standard time–stepping schemes for nonsmooth dynamics (Moreau–Jean timestepping scheme) have been proposed in the last decade, in terms of accuracy and dissipation properties [23], [24], [87], [88], [25], [54], [50], [89], [52]. An important part of these schemes has been developed by members of the BIPOP team and has been implemented in the Siconos software (see Sect. 5.1).

Numerical solution procedure for the time–discretized problem, mainly through wellidentified problems studied in the optimization and mathematical programming community. Another very interesting feature is the fact that the discretized problem that we have to solve at each time–step is generally a wellknown problem in optimization. For instance, for LCSs, we have to solve a linear complementarity problem [55] for which there exist efficient solvers in the literature. Comparing to the brute force algorithm with exponential complexity that consists in enumerating all the possible modes, the algorithms for linear complementarity problem have polynomial complexity when the problem is monotone.
In the Axis 2 of the research program (see Sect. 3.3), we propose to perform new research on the geometric timeintegration schemes of nonsmooth dynamical systems, to develop new integration schemes for Boundary Value Problem (BVP), and to work on specific methods for two timediscretized problems: the Mathematical Program with Equilibrium Constraints (MPEC) for optimal control and Second Order Cone Complementarity Problems (SOCCP) for discrete frictional contact systems.