Xi-Ming M Sun

Xue-Fang F Wang

Frédéric Mazenc
email: frederic.mazenc@l2s.centralesupelec.fr

Delay-hybrid-dependent stability for systems with large delays

Keywords: Stability, Large delays, Switching, Lyapunov-Krasovskii functional

This paper investigates a problem of stability analysis for a class of nonlinear systems with a time-varying delay taking both large and small values in an alternating manner. This problem is a delay-dependent stability analysis one but it shares some features with delay-independent ones. Hence traditional stability analysis techniques for systems with delay do not apply. By building upon our previous works, we first introduce the concept of delay-hybriddependent stability, which grasps the features of the delays described above. To ease the problem we investigate, we represent the studied system as a system with a switched delay. Then by using switching techniques and Lyapunov-Krasovskii functionals (LKFs), we provide a new stability criterion. Next, in the linear context, we show how the needed LKFs can be constructed by developing an original design of LKFs.

I. INTRODUCTION

Time-delays are unavoidable in many practical control systems such as, networked control systems (NCSs) [START_REF] Heemels | Networked control systems with communication constraints: Tradeoffs between transmission intervals, delays and performance[END_REF], temperature control systems [START_REF] Sun | Absolutely exponential stability and temperature control for gas chromatograph system under dwell time switching techniques[END_REF], multi-agent systems [START_REF] Tang | Leader-following consensus of a class of stochastic delayed multi-agent systems with partial mixed impulses[END_REF], [START_REF] Wu | Stability analysis of stochastic delayed systems with an application to multi-agent systems[END_REF], sampled control systems [START_REF] Mazenc | Reduction model approach for linear systems with sampled delayed inputs[END_REF], biological systems [START_REF] Müller | Estimation of delay times in biological systems[END_REF], pneumatic systems [START_REF] Hong | Stabilizing network control for pneumatic systems with time-delays[END_REF]. The presence of time-delays often leads to poor performances and even instability [START_REF] Gu | Stability of time-delay systems[END_REF]. Consequently, the stability analysis of time-delay systems has attracted considerable attention in last two decades (see, for example, [START_REF] Park | A delay-dependent stability criterion for systems with uncertain time-invariant delays[END_REF]- [START_REF] Karimi | Local capacity H ∞ control for production networks of autonomous work systems with time-varying delays[END_REF] and references cited therein) in both the delay-independent and delay-dependent stability context. It is well-known that delay-dependent stability conditions are in general less conservative than delayindependent ones, especially for systems with small delays. Let us recall that the stability of a system is called delay-independent if this system is stable no matter what the size of its delay is and it is called delay-dependent if the stability property holds for delays belonging to specific intervals only. Thus, the maximum allowable size of the delay is a fundamental parameter for assessing the conservatism of a delay-dependent stability criterion. It is worth observing that many works have been done to reduce the conservatism of delay-dependent stability results for time-delay systems (see, for example, [START_REF] Tang | Leader-following consensus of a class of stochastic delayed multi-agent systems with partial mixed impulses[END_REF], [START_REF] Wu | Stability analysis of stochastic delayed systems with an application to multi-agent systems[END_REF], [START_REF] Medvedeva | Synthesis of razumikhin and lyapunov-krasovskii approaches to stability analysis of time-delay systems[END_REF]- [START_REF] Karimi | Stability analysis and control synthesis of neutral systems with time-varying delays and nonlinear uncertainties[END_REF]).

However, with the development of NCSs [START_REF] Gao | Stabilization of networked control systems with a new delay characterization[END_REF], [START_REF] Tang | Tracking control of networked multi-agent systems under new characterizations of impulses and its applications in robotic systems[END_REF], uncertain kinematics systems [START_REF] Xiao | Tracking control of robotic manipulators with uncertain kinematics and dynamics[END_REF], [START_REF] Wang | Adaptive neural control of nonlinear systems with unknown control directions and input dead-zone[END_REF], fault tolerant control [START_REF] Xiao | Velocity-free fault tolerant and uncertainty attenuation control for a class of nonlinear systems[END_REF], [START_REF] Blanke | Diagnosis and fault-tolerant control[END_REF], and controller/actuator failure systems [START_REF] Xiao | Attitude stabilization of spacecraft under actuator saturation and partial loss of control effectiveness[END_REF], [START_REF] Sun | Stability of systems with controller failure and time-varying delay[END_REF], a new class of time-delay systems, whose stability cannot be established by the traditional stability criteria due to the presence of time-varying delays, has to be considered. For instance, in NCSs, due to the effect of the consecutive packet dropouts and the application of zero-order hold, in some local time periods, there will appear large delays (see, for example [START_REF] Sun | Stability analysis for networked control systems based on average dwell time method[END_REF]) whose upper bound may be larger than the maximum allowable constant delay that is deduced from classical delay-dependent stability criteria [START_REF] Gu | Stability of time-delay systems[END_REF], [START_REF] Xu | Improved delay-dependent stability criteria for time-delay systems[END_REF], [START_REF] Pang | Data-based predictive control for networked nonlinear systems with network-induced delay and packet dropout[END_REF]. In the research of controller/actuator failure problems, the controller or actuator failure may lead to locally large delay periods [START_REF] Sun | Stability of systems with controller failure and time-varying delay[END_REF]. In addition, with the development of the Internet of things, the large delays may appear during the process of the data transmission to the terminal (see, for example, [START_REF] Hale | Cloud-based centralized/decentralized multi-agent optimization with communication delays[END_REF]). Hence, it is of both theoretical and practical significance to investigate this kind of system model and make an analysis of such systems. In fact, some results of such systems has been successfully applied to multi-agent systems (see, for example, [START_REF] Wu | Stability analysis of stochastic delayed systems with an application to multi-agent systems[END_REF]). Indeed, the delay of these systems presents a very specific property: roughly speaking, in some time intervals, it is small and in others it is large. The stability of these systems is not surprising: the intuition suggests that if a system is stable when its delay is small, then if its delay is small in sufficiently long and frequent time intervals, the systems may be stable, no matter how large is the delay in some other intervals. This motivates the present paper where a new stability analysis for systems of this type is proposed.

For a system with a delay that is sometimes large and sometimes small, we call the two different delay periods large delay periods (LDPs) and small delay periods (SDPs) and we call the system itself system with mixed delay periods.

The key ideas of our new result are the following. First, we observe that the change from one delay period to the other can be described as a switching behavior. By representing the delay as a switched delay, we represent a system with mixed delay periods as a switched system and call the stability of such a system delay-hybrid-dependent stability. This crucial new idea will enable us to study the stability of these systems by resorting to the switched system theory, for which several efficient tools are available (for more information on switched systems, see for instance [START_REF] Liberzon | Switching in systems and control[END_REF]- [START_REF] Caliskan | Dwell-time computation for stability of switched systems with time delays[END_REF] and references cited therein). It turns out that the switched system we obtain is composed of a stable and an unstable subsystems. The fundamental idea we will use to study this system consists in obtaining for the solutions a decay rate estimate for the stable subsystem, an increase rate estimate of the solutions for the unstable one and next to use these estimates to deduce the stability of the system under assumptions pertaining to the switching sequences (basically a restriction on the frequency and the length of LDPs is imposed).

The present paper can be seen as a review, comprehensive extension and simplification of our previous works. Let us explain how. We have done some related works [START_REF] Sun | Delay-dependent stability for discrete systems with large delay sequence based on switching techniques[END_REF]- [START_REF] Sun | Input-to-state stability for nonlinear systems with large delay periods based on switching techniques[END_REF] to deal with the presence of large delays. In [START_REF] Sun | Delay-dependent stability for discrete systems with large delay sequence based on switching techniques[END_REF], the large delay sequence concept was firstly introduced for analyzing the stability of discrete-time systems with time-delays and the normal delay system was first turned to a class of switched delay system. Then, the concepts of the length rate and frequency of the large delay period were extended to linear continuous system in [START_REF] Sun | Stability analysis for systems with large delay period: A switching method[END_REF]. Furthermore, the concept of the large delay period was extended to nonlinear delay systems to establish an Input-to-State Stability (ISS) property in [START_REF] Sun | Input-to-state stability for nonlinear systems with large delay periods based on switching techniques[END_REF]. Let us observe that the methods used in these contributions to construct Lyapunov functionals and to analyze the stability of systems with large delay are rather complex in papers mentioned above [START_REF] Sun | Delay-dependent stability for discrete systems with large delay sequence based on switching techniques[END_REF]- [START_REF] Sun | Input-to-state stability for nonlinear systems with large delay periods based on switching techniques[END_REF]. In addition, the given stability conditions are all based on LMI forms which cannot incarnate the inherent features of Lyapunov functions and the computation on LMIs may be a little bit complex. To obtain simpler results, in the present paper, for the first time, we propose the delay-hybrid-dependent stability notion, which is different from the classical delay-dependent and delay-independent stability notions. In certain applications, the new stability concept is very helpful, such as networked controlled systems, the standard delay-dependent and delay-independent analysis are too restrictive. This is because the delay signal might have large peaks, while the frequency of these large peaks may be very low. In such scenarios, the usage of the delayhybrid-dependent stability notion might be helpful in reducing the conservatism in the results. This new stability concept can be viewed as the summarization of our previous stability definitions since in our previous works, we in fact gave some examples related to large delays and some analysis results for the large delayed system but did not propose the delay-hybrid-dependent stability criterion clearly. Moreover, with a view to further studies, pertaining for example to controller designs, we give a constructed Lyapunov-Krasovskii functional technique in the paper to conveniently check the delay-hybrid-dependent stability, which is never used in the analysis of large delayed case before although some of techniques may be used to analyze the small delay case. Besides, for a class of linear systems ẋ(t) = Ax(t) + Bx(td σ (t) (t)), the designed Lyapunov-Krasovskii functional technique proposed in the present paper does not assume that the matrix A is Hurwitz, which is more efficient in practice compared with our previous works. So such a technique gives a relaxed condition compared with our previous works given in LMI forms, and may be more easier to utilize or extend the delay-hybrid-dependent stability. Finally, we illustrate our theoretical results with two examples.

The remainder of this paper is organized as follows. Section II introduces the problem and some definitions. Section III gives the main results. Section IV demonstrates the effectiveness of the proposed method through two examples. Finally, in Section V we draw some concluding remarks.

A. Notation

The argument of the functions will be omitted whenever no confusion can arise from the context. For a vector x ∈ R n , the symbol | • | denotes the Euclidean norm of x. P > 0 is used to denote a symmetric positive definite matrix P ∈ R n×n and λ min (P) denotes the minimum eigenvalue of P. L 2 ([-τ, 0]; R n ) is used to denote the space of the square integrable functions φ :

[-τ, 0] → R n . For any function φ ∈ C 0 ([-τ, 0], R n ), we let |φ | = sup θ ∈[-τ,0] |φ (θ )|.
For a given fixed time-delay τ, we define x t (θ ) = x(t + θ )(θ ∈ [-τ, 0]). The space of functions W [-τ, 0] denotes the Banach Space of absolutely continuous functions φ :

[-τ, 0] → R n with φ ∈ L 2 ([-τ, 0]; R n ), equipped with the norm |φ | W = max θ ∈[-τ,0] |φ (θ )| + 0 -τ | φ (s)| 2 ds 1 2 .
For a given fixed time-delay τ, we define x t (θ ) = x(t + θ )(θ ∈ [-τ, 0]). We end this section with a definition of derivative for functionals, which will be instrumental in the forthcoming Lyapunov analyzes.

Let

V (t, φ , ϕ) : R ≥0 × E in → R + with E in = W [-τ, 0] × L 2 ([-τ, 0], R n
) be a continuous functional. Define (see, e.g., [START_REF] Fridman | On input-to-state stability of systems with time-delay: A matrix inequalities approach[END_REF])

V (t, φ , φ ) = lim sup h→0 + V (t + h, x t+h (t, φ ), ẋt+h (t, φ )) -V (t, φ , φ ) h , (1)
where x t (t 0 , φ ), for t ≥ t 0 , is a solution with the initial condition x t 0 = φ ∈ W .

II. PROBLEM FORMULATION

A. The studied systems

We consider the nonlinear system:

ẋ(t) = f (x(t), x(t -d(t))) (2)
where x ∈ R n is the state vector, f :

R n × R n → R n is of class C 1 with f (0, 0) = 0, d(t)
(t i ) such that 0 ≤ t 0 < t 1 < t 2 < • • • and for all t ∈ [t 2 j ,t 2 j+1 ), d(t) ≤ h. (3)
Moreover, there is ν > 0 such that for all i ∈ N, t i+1t i ≥ ν.

Remark 1:

The case where it is assumed that d(t) > h when t ∈ [t 2 j+1 ,t 2 j+2 ) is entirely similar to the one of Assumption 1 and so we will not consider it.

We define for later use two sets:

E s = +∞ i=0 [t 2i ,t 2i+1 ) and E l = +∞ i=0 [t 2i+1 ,t 2i+2 ) (4)
and we recall the following definitions: Definition 1: [START_REF] Sun | Stability analysis for systems with large delay period: A switching method[END_REF] Under Assumption 1, the time interval [t 2 j ,t 2 j+1 ) for any j ∈ N is called small delay period and the time interval [t 2 j+1 ,t 2 j+2 ) for any j ∈ N is called large delay period. We introduce two definitions:

Definition 2: [38] For any T 2 > T 1 ≥ 0, T ld p (T 1 , T 2 ) denotes the total time length of LDPs in time interval [T 1 , T 2 ) and T sd p (T 1 , T 2 ) denotes the total time length of SDPs in [T 1 , T 2 ).
Definition 3: For any T 2 > T 1 ≥ 0, N s (T 1 , T 2 ) denotes the number of t i belonging to (T 1 , T 2 ). Definition 4: We say that system (2) is a system with mixed delay periods when the delay is defined over the interval [t 0 , +∞) and operates in an alternating manner between SDPs and LDPs. The stability of such a system shares features with the delay-dependent and delay-independent stability and this stability property is called delayhybrid-dependent stability and we say that the system is delay-hybrid-dependently stable. Next based on the above description and definitions, hybrid representation of the studied systems is given.

B. Hybrid representation of the system

When LDPs and SDPs occur in alternative way in system (2), switching techniques can be adopted to analyze its delay-hybrid-dependent stability properties. To do this, we rewrite system (2) as the following switched system:

ẋ(t) = f (x(t), x(t -d σ (t) (t))) x(θ ) = φ (θ ), θ ∈ [-d b , 0] (5)
where φ denotes a continuously differentiable vector-valued initial function defined over [-d b , 0], σ : [0, +∞) → M = {1, 2} is a function called switching signal, defined as follows: σ (t) = 1 when t ∈ E s and σ (t) = 2 when t ∈ E l and d(t) := d 1 (t) when t ∈ E s and d(t) := d 2 (t) when t ∈ E l . Thus d(t) := d σ (t) (t) for all t ≥ 0, which implies that system (5) is identical to system (2).

In next section, stability criteria with some restrictions on the time length and the number of occurrence of large delays of nonlinear and linear systems are given.

III. DELAY-HYBRID-DEPENDENT STABILITY RESULTS

A. General stability result

We give in this section delay-hybrid-dependent stability conditions for the general nonlinear system (5): we are ready to state and prove a stability result which, basically, shows that the exponential stability of system ( 5) is guaranteed under certain restrictions on the length and the frequency of the LDPs. Before giving the main theorem, we first present some conditions:

C 1 : The inequalities γ 1 |φ (0)| 2 ≤ V j (t, φ , φ ) ≤ γ 2 |φ | 2 W , j = 1, 2 (6)
and

V 1 (t, φ , φ ) ≤ µV 2 (t, φ , φ ), V 2 (t, φ , φ ) ≤ µV 1 (t, φ , φ ) (7)
hold for all t ≥ 0 and

φ ∈ W [-τ, 0]. C 2 :
The derivatives of V 1 and V 2 along the trajectories of (5) satisfy

V1 (t, φ , φ ) ≤ -α 1 V 1 (t, φ , φ ), ∀t ∈ E s , V2 (t, φ , φ ) ≤ α 2 V 2 (t, φ , φ ), ∀t ∈ E l . (8)
C 3 : Let R(t 0 ,t) = -α 1 T sd p (t 0 ,t) + α 2 T ld p (t 0 ,t). (9)
There are α * > 0, β * > 0, η i > 0, i = 1, 2 such that for all t ≥ t 0 , the inequalities

R(t 0 ,t) ≤ -α * (t -t 0 ) + η 1 , N s (t 0 ,t) ≤ β * (t -t 0 ) + η 2 (10)
are satisfied.

C 4 : The inequality ln(µ)β * < α * (11)
is satisfied.

Theorem 1: Let system (5) satisfy Assumption 1. Assume that there are two functionals V 1 and V 2 satisfying the properties described in (1) and positive real numbers γ i , α i , i = 1, 2 and µ ≥ 1 such that conditions (6)- [START_REF] Aleksandrov | Delay-independent stability conditions for some classes of nonlinear systems[END_REF] hold. Then for all t ≥ t 0 ,

|x(t)| 2 ≤ e (ln(µ)β * -α * )(t-t 0 ) γ * |x t 0 | 2 W with γ * = e η 1 +ln(µ)(η 2 +1) γ 2 γ 1 . Proof: See Appendix A. Remark 2:
The fact that the inequalities in [START_REF] Hong | Stabilizing network control for pneumatic systems with time-delays[END_REF] need to be satisfied with the same constant µ is not restrictive since if two Lyapunov functionals V 1 and V 2 are such that

V 1 (t, φ , φ ) ≤ µ 1 V 2 (t, φ , φ ), V 2 (t, φ , φ ) ≤ µ 2 V 1 (t, φ , φ ) (12)
then, the functionals

V 1 = µ 2 µ 1 V 1 , V 2 = V 2
satisfy the inequalities [START_REF] Hong | Stabilizing network control for pneumatic systems with time-delays[END_REF] with µ = √ µ 1 µ 2 . Moreover, ( 7) is just used to restrain the boundedness of V 1 , V 2 and make sure it will not go to infinity at switching points. So it is enough to adopt the same µ on V 1 , V 2 . Although we use different µ in [START_REF] Hong | Stabilizing network control for pneumatic systems with time-delays[END_REF], in order to obtain exponential stability by iterative method, we need to take the maximum value among different µ. Besides, for positive definite functions V 1 , V 2 , the existence of the ratio bound µ in ( 7) is a sufficient condition for the existence of the comparison functions γ 1 , γ 2 in [START_REF] Müller | Estimation of delay times in biological systems[END_REF].

For clearity of this idea, by constructing Lyapunov functionals, in next subsection, the stability analysis of the linear system with mixed delay periods is shown.

B. Illustration

In this section, we show how Theorem 1 can be applied to linear systems with mixed delay periods by constructing LKFs satisfying conditions C 1 and C 2 . For the sake of simplicity, we present a result for time invariant systems that is based on LKFs which can be constructed via rather simple calculations. More refined, but more lengthy, constructions leading to less conservative results with respect to the size of the delays in small periods can be proposed.

Consider the system ẋ(t) = Ax(t)

+ A d x(t -d σ (t) (t)), (13)
where d 1 , d 2 and σ are the same as those in Section II-B and the matrices A and A d are constant and such that there are a constant c a > 0 and a symmetric positive definite matrix P such that

PH + H P ≤ -2c a P , (14)
where H = A + A d and I is the identity matrix. To simplify the calculations, we assume that d 2 is of class C 2 and 2 . Moreover, we assume that

d < 1 and | d2 (t)| ≤ d with d = 2|PA| + 1 4 |A| 2 + |PA d | + |PA d |+ 1 4 |A d | 2 1-d (1 -d )
h ≤ min 1, 1 4 √ c b , c a λ min (P) 16c b |H| 2 , 1 8|A d | 2 , 1 3c a ( 15
)
with

c b = |A T d PA d | c a .
Remark 3: A crucial aspect of the family we consider is that we do not assume that the matrix A is Hurwitz, which is frequently encountered in practice and implies that in general this family is not delay independent stable. Now, we show how we can construct Lyapunov functionals ensuring that the conditions C 1 and C 2 of Theorem 1 are satisfied by system [START_REF] Karimi | H ∞ control of parameter-dependent state-delayed systems using polynomial parameter-dependent quadratic functions[END_REF].

First step. Let us study system (13) when t ∈ E s . Let

ν 1 (x) = x Px. ( 16
) Since ẋ(t) = Hx(t) -A d t t-d 1 (t) ẋ( )d , (17)
the derivative of ν 1 along (17) satisfies

ν1 (t) ≤ -2c a x(t) Px(t) -2x(t) PA d t t-d 1 (t) ẋ( )d ≤ -c a x(t) Px(t) + 1 c a |A T d PA d | t t-d 1 (t) ẋ( )d 2 ≤ -c a x(t) Px(t) + c b h t t-h | ẋ( )| 2 d .
Let us introduce the operator

ν 2 (x t , ẋt ) = ν 1 (x(t)) + 2c b h t t-h t m | ẋ( )| 2 d dm. Then ν2 (t) ≤ -c a ν 1 (x(t)) -c b h t t-h | ẋ( )| 2 d + 2c b h 2 Hx(t) -A d t t-d σ (t) (t) ẋ( )d 2 ≤ -c a ν 1 (x(t)) -c b h t t-h | ẋ( )| 2 d + 4c b h 2 |H| 2 |x(t)| 2 + |A d | 2 h t t-h | ẋ( )| 2 d .
From ( 16), it follows that

ν2 (t) ≤ -c a ν 1 (x(t)) -c b h t t-h | ẋ( )| 2 d + 4c b h 2 |H| 2 λ min (P) ν 1 (x(t)) + h|A d | 2 t t-h | ẋ( )| 2 d . Since (15) implies that 16c b h 2 |H| 2 λ min (P) ≤ c a and 8h 2 |A d | 2 ≤ 1 the inequality ν2 (t) ≤ - 3c a 4 ν 1 (x(t)) - c b 2 h t t-h | ẋ( )| 2 d is satisfied. Noticing that 1 2c b h 2 (ν 2 (x t , ẋt ) -ν 1 (x(t))) ≤ t t-h | ẋ( )| 2 d , we deduce that ν2 (t) ≤ - 3c a 4 ν 1 (x(t)) - 1 4h 2c b h t t-h t m | ẋ( )| 2 d dm. Therefore ν2 (t) ≤ -3c c ν 2 (x t , ẋt ), (18)
with c c = c a 4 because h ≤ 1 3c a .

Second step. Now, let us study system (13) when t ∈ E l . Immediate calculations give

ν2 (t) ≤ 2x(t) P[Ax(t) + A d x(t -d 2 (t))] + 2c b h 2 | ẋ(t)| 2 ≤ 2x(t) P[Ax(t) + A d x(t -d 2 (t))] + 2c b h 2 |Ax(t) + A d x(t -d 2 (t))| 2 ≤ (2|PA| + 4c b h 2 |A| 2 + |PA d |)|x(t)| 2 + (|PA d | + 4c b h 2 |A d | 2 )|x(t -d 2 (t))| 2 . Since h ≤ 1 4 √ c b , the inequality ν2 (t) ≤ c d ν 1 (x(t)) + c e ν 1 (x(t -d 2 (t))) with c d = 2|PA|+ 1 4 |A| 2 +|PA d | λ min (P)
and c e =

|PA d |+ 1 4 |A d | 2
λ min (P) , is satisfied. Let

ν 3 (t, x t , ẋt ) = ν 2 (x t , ẋt ) + c e 1 -ḋ2 (t) t t-d 2 (t) ν 1 (x(m))dm. Then ν3 (t) ≤ c d ν 1 (x(t)) + c e 1-ḋ2 (t) ν 1 (x(t)) + c e d2 (t) (1-ḋ2 (t)) 2 t t-d 2 (t) ν 1 (x(m))dm. ≤ c f ν 1 (x(t)) + c e d (1-d ) 2 t t-d 2 (t) ν 1 (x(m))dm ≤ c f ν 3 (t, x t , ẋt ), with c f = c d + c e 1-d .
Third step. It is very difficult to find or it does not exist two constants µ 1 > 0 and µ 2 > 0 such that ν 2 and ν 3 satisfy inequalities of the type ( 12), therefore they cannot be directly used to construction LKFs ensuring that the conditions C 1 and C 2 are satisfied. This leads us to construct a new LKFs in the SDPs. Let us consider again the case where t ∈ E s . Let

ν 4 (t, x t , ẋt ) = ν 2 (x t , ẋt ) + c c t t-d 2 (t) ν 1 (x(m))dm + c c 1 + d b t t-d 2 (t) t s ν 1 (x(m))dmds.
Then, from [START_REF] Wei | A novel memory filtering design for semi-markovian jump time-delay systems[END_REF] and the fact that d < 1, it follows that

ν4 (t) ≤ -2c c ν 2 (x t , ẋt ) -c c (1 -ḋ2 (t))ν 1 (x(t -d 2 (t)) + c c 1+d b d b ν 1 (x(t)) -c c 1+d b (1 -ḋ2 (t)) t t-d 2 (t) ν 1 (x(m))dm ≤ -c c 2+d b 1+d b ν 2 (x t , ẋt ) + 1-d 2+d b t t-d 2 (t) ν 1 (x(m))dm . Using ν 4 (t, x t , ẋt ) ≤ ν 2 (x t , ẋt ) + c c 1+2d b 1+d b t t-d 2 (t) ν 1 (x(m))dm, (19)
we obtain: ν4 (t) ≤ -c g ν 4 (t, x t , ẋt ),

with c g = c c 2+d b 1+d b min 1, (1+d b )(1-d ) c c (1+2d b )(2+d b ) . Now, observe that ν 3 (t, x t , ẋt ) ≤ ν 2 (x t , ẋt ) + c e 1-d t t-d 2 (t) ν 1 (x(m))dm, ν 2 (x t , ẋt ) + c e t t-d 2 (t) ν 1 (x(m))dm ≤ ν 3 (t, x t , ẋt ), ν 2 (x t , ẋt ) + c c t t-d 2 (t) ν 1 (x(m))dm ≤ ν 4 (t, x t , ẋt ).
These inequalities and [START_REF] Zhao | Adaptive fuzzy hierarchical sliding-mode control for a class of mimo nonlinear time-delay systems with input saturation[END_REF] imply that ν 3 (t, x t , ẋt ) ≤ µ 1 ν 4 (t, x t , ẋt ) and ν 4 (t, x t , ẋt ) ≤ µ 2 ν 3 (t, x t , ẋt )

with µ 1 = max 1, c e c c (1-d )
and

µ 2 = max 1, c c (1+2d b ) c e (1+d b ) . Final Step. The conditions C 1 and C 2 of Theorem 1 are satisfied with V 1 (t, φ , ϕ) = µ 1 µ 2 ν 4 (t, φ , ϕ) , V 2 (t, φ , ϕ) = ν 3 (t, φ , ϕ)
and α 1 = c g and α 2 = c f . Indeed, for all (φ , ϕ) ∈ E in and t ≥ 0,

V 1 (t, φ , ϕ) ≤ µ 1 µ 2 µ 2 V 2 (t, φ , ϕ) = µV 2 (t, φ , ϕ) V 2 (t, φ , ϕ) ≤ µ 1 ν 3 (t, φ , ϕ) = µV 1 (t, φ , ϕ) with µ = √ µ 1 µ 2 .
Remark 4: From this illustration, one can know that our method is more effective in practice from four aspects: 1) It has less computational complexity.

2) The forms of Lyapunov functions are very simple and the stability criteria in present paper are easily checked.

3) It does not assume that the matrix A is Hurwitz for linear system. 4) It has robustness which will be shown in the following example 2.

IV. EXAMPLES

A. Example 1

In this section, we use an example to show the effectiveness of our results. We consider the system ẋ = -x(td σ (t) (t)) [START_REF] Wang | Approximation-based adaptive fuzzy tracking control for a class of nonstrict-feedback stochastic nonlinear time-delay systems[END_REF] where x ∈ R and d σ (t) (t) is a periodic function of period 20 such that d 1 (t) = 1 9 and d 2 (t) = d b > 1 9 for all t ≥ 0 and t 0 = 0, t 1 = p, p ∈ (0, 20). We adopt the notation of the previous sections and choose

ν 1 (x) = x 2 . Then A = 0, A d = -1, H = -1, d = 0, c a =

B. Example 2

In this subsection, we use the same example in [START_REF] Sun | Stability analysis for systems with large delay period: A switching method[END_REF] to compare with our previous results. Consider system (13) with A = -2 0 0 -0.9

,

A d = -1 0 -1 -1 .
Here we give α 1 = 0.038, α 1 shows that the considered system is exponentially stable without considering large delay case. Fig. 2 shows that system is not stable when only consider large delay case. Fig. 3 shows that, although large delays occur, the system is still exponentially stable under a certain restriction on the time length and the number of occurrence of large delays. In Fig. 3, switching signal σ (t) is random one of some switching signals satisfying inequalities [START_REF] Xu | Improved delay-dependent stability criteria for time-delay systems[END_REF]. From simulation, it can be seen that T ld p ≤ 42.2 and N s = 18. It means that the LDP can be permitted to occur nine times in each 100s (second) and the permitted length of LDP can reach to 42.2. If choose the same parameters α 1 = 0.14, α 2 = 0.3 with [START_REF] Sun | Stability analysis for systems with large delay period: A switching method[END_REF], we can obtain the same results. However in [START_REF] Sun | Stability analysis for systems with large delay period: A switching method[END_REF], it only can be permitted to occur once in each 100s and the permitted length of LDP just only reach to 13.1818. It is clearly shown the method proposed in this paper is more effective. In order to show the robustness of our method, we consider the following system

ẋ(t) = Ax(t) + A d x(t -d(t)) + Dω(t)
where A and A d are the same as the above, D is constant matrix, ω(t) is the disturbance input.

In the simulation, we set ω(t) = 0.5e -0.5t and D = [1 1] T . The simulation results are given in Fig. 4, Fig. 5 and Fig. 6. It is clear to see from figures 4, 5 and 6 that the method proposed in the present paper has certain robustness.

V. CONCLUSION

This paper has studied the stability problem of systems with a delay of a type frequently encountered in engineering problems: systems with mixed delay periods, i.e. time-varying delays with large values in small intervals and small values in long intervals. We have represented systems with this type of delay as switched systems. For them, we proposed a general stability result based on Lyapunov functionals and gave a new stability concept which is called as delay-hybrid-dependent stability. In certain applications, this new stability notion might be helpful in reducing the conservatism in the results. Numerical results have been used to illustrate the feasibility and robustness of the proposed method.

There are some remaining issues, which can be the subjects of future works, e.g., how to extend this result the stochastic delay problem (see for example, [START_REF] Tang | Leader-following consensus of a class of stochastic delayed multi-agent systems with partial mixed impulses[END_REF], [START_REF] Wu | Stability analysis of stochastic delayed systems with an application to multi-agent systems[END_REF]), how to design observer and how to propose approach for systems of neutral type in the presence of switched time-varying delays. These topics motivate our future research.

APPENDIX

First, let us observe that the properties C 1 and C 2 ensure that system (5) is forward complete. Now, based on the Lyapunov functionals V 1 (t, φ , φ ) and V 2 (t, φ , φ ) defined in [START_REF] Müller | Estimation of delay times in biological systems[END_REF], we propose a piecewise Lyapunov functional for system [START_REF] Mazenc | Reduction model approach for linear systems with sampled delayed inputs[END_REF]. Let V (t) = V σ (t) (t, φ , φ ).

According to (8), we have

V (t) ≤ e -α
for all s ∈ N. Consider t 0 ≤ t with t ∈ [t 2 j+1 ,t 2 j+2 ) and t 0 ∈ [0,t 2 j+1 ). In view of ( 7) and ( 21), we deduce that, V (t) ≤ e α 2 (t-t 2 j+1 ) V 2 (t 2 j+1 ) ≤ e α 2 (t-t 2 j+1 ) µV 1 (t - 2 j+1 ) ≤ e α 2 (t-t 2 j+1 ) µe -α 1 (t 2 j+1 -t 2 j ) V 1 (t 2 j ) ≤ • • • ≤ µ N s (t 0 ,t) e -α 1 T sd p (t 0 ,t) e α 2 T ld p (t 0 ,t) V q (t 0 ), where q = 1 or 2. Next, considering the case where t ∈ [t 2 j ,t 2 j+1 ), we deduce easily that for all t ≥ t 0 ≥ 0, V (t) ≤ µ N s (t 0 ,t)+1 e R(t 0 ,t) V 1 (t 0 ), where R is the function defined in [START_REF] Park | A delay-dependent stability criterion for systems with uncertain time-invariant delays[END_REF]. Next, using [START_REF] Xu | Improved delay-dependent stability criteria for time-delay systems[END_REF], we deduce that V (t) ≤ e ln(µ)β * (t-t 0 )+ln(µ)(η 2 +1) e -α * (t-t 0 )+η 1 V 1 (t 0 ).

We deduce from C 1 that γ 1 |x(t)| 2 ≤ V (t) ≤ e (ln(µ)β * -α * )(t-t 0 ) e η 1 +ln(µ)(η 2 +1) γ 2 |x t 0 | 2 W . Clearly, the same inequalities are satisfied for t 0 ≤ t such that there is k ∈ N such that t 0 ∈ [t k ,t k+1 ) and t ∈ [t k ,t k+1 ). This allows us to conclude.

  denotes nonnegative time-varying delay. Along the paper, we let system (2) satisfy the following assumption: Assumption 1: The function d(t) is piecewise-C 1 and there are constants d ≥ 0 and d b > 0 such that | ḋ(t)| ≤ d almost everywhere and d(t) ≤ d b for all t ≥ 0. There is a constant h ∈ [0, d b ) and a sequence

√ 5 , α 1 = c g and α 2 = 9 4 .

 514 1, c b = 2, c c = 1 4 , c d = 1 and c e = 5 4 , c f = 9 4 , c g = 2+d b 4(1+d b ) min 1, 4 1+d b(1+2d b )(2+d b ) , µ = We have, for all t 0 ≥ 0, t ≥ t 0 , R(t 0 ,t) = -c g T sd p (t 0 ,t) + 9 4 T ld p (t 0 ,t) ≤ -α * (tt 0 ) + η 1 with α * = 2+d b 4(1+d b ) min 1, 4 1+d b (1+2d b )(2+d b ) p -9 4 (20p), η 1 = d b and N s (t 0 ,t) ≤ β * (tt 0 ) + 2with β * = 1 20 and η 2 = 2. We deduce from (11) that if p is such thatln( √ 5) 1 20 < 2+d b 4(1+d b ) min 1, 4 1+d b (1+2d b )(2+d b ) p -94 (20p) then the origin of (20) is GES. This condition rewrites as lnNotice that when d b = 1 2 , the system ẋ = -x(td 2 (t)) is unstable and system (20) is GES when the inequality 3 satisfied, which is satisfied if 17 ≤ p.

= 33. 8 ,

 8 d b = 10. Choose α * = 0.1, β * = 0.035, c a = 1.899, η 1 = 1450, η 2 = 16, then we obtain P = 4.1802 -2.2091 * 2.4372 , c b = 1.356, c c = 0.475, c d = 25.97, c e = 4.487, c f = 33.8, c g = 0.038, µ 1 = 11.8, µ 2 = 1, µ = 3.435. From (15) we have h ≤ min{1, 0.214, 0.117, 0.21, 0.175}. The small delays and large delays are set as d 1 (t) = 0.1|sin(t)| and d 2 (t) = 0.5|sin(t)| + 9.5, respectively. The simulation results are shown in Fig. 1, Fig. 2 and Fig. 3. Fig.

Fig. 1 .Fig. 2 .Fig. 3 .

 123 Fig. 1. State trajectories under small delay: d 1 (t) = 0.1|sint|.

Fig. 4 .

 4 Fig. 4. State trajectories under small delay: d 1 (t) = 0.1|sint|.

Fig. 5 .

 5 Fig. 5. State trajectories under large delay: d 2 (t) = 9.5 + 0.5|sint|.

Fig. 6 .

 6 Fig. 6. State trajectories under switching signal σ (t).

  1 (t-t 2s ) V (t 2s ), t ∈ [t 2s ,t 2s+1 ), e α 2 (t-t 2s+1 ) V (t 2s+1 ), t ∈ [t 2s+1 ,t 2s+2 ),

This work were supported by the National Natural Science Foundation of China under Grant 61325014 and the PHC Bosphore 2016. Title: Switched nonlinear and time-varying systems with arbitrarily long delay.


Jean-David Benamou
email: jean-david.benamou@inria.fr

AND FRANC Thomas O Gallou Ët
email: thomas.gallouet@inria.fr

¸ois-Xavier Vialard
email: fxvialard@normalesup.org

SECOND ORDER MODELS FOR OPTIMAL TRANSPORT AND CUBIC SPLINES ON THE WASSERSTEIN SPACE

On the space of probability densities, we extend the Wasserstein geodesics to the case of higher-order interpolation such as cubic spline interpolation. After presenting the natural extension of cubic splines to the Wasserstein space, we propose a simpler approach based on the relaxation of the variational problem on the path space. We explore two different numerical approaches, one based on multi-marginal optimal transport and entropic regularization and the other based on semi-discrete optimal transport.

Introduction

We propose a variational method to generalize cubic splines on the space of densities using multimarginal optimal transport. In short, the proposed method consists in minimizing, on the space of measures on the path space, under marginal constraints, the norm squared of the acceleration. In this setting, we show that two numerical approaches, classical in optimal transportation can be applied. One is based on entropic regularization and the Sinkhorn Algorithm, the other relies on the Semi-Discrete formulation of Optimal Transportation and the computation of Laguerre cells, a classical problem in computationnal geometry. We showcase our methodology on 1D and 2D data.

In the past few years, higher-order interpolations methods have been investigated for applications in computer vision or medical imaging, for time-sequence interpolation or regression. The most usual setting is when data are modeled as shapes, which can be understood as objects embedded in the Euclidean space with no preferred parametrization: space of unparametrized curves or surfaces, or images are some of the most important examples. These examples are infinite dimensional but the finite dimensional case of a Riemannian manifold was interesting for camera motion interpolation as first introduced in [START_REF] Noakes | Cubic splines on curved spaces[END_REF] and further developed in [START_REF] Camarinha | Splines of class C k on non-euclidean spaces[END_REF][START_REF] Crouch | The dynamic interpolation problem: On Riemannian manifold, Lie groups and symmetric spaces[END_REF]. Motivated by different applications, the problem of interpolation between two shapes is usually treated via the use of a Riemannian metric on the space of shapes and computing a geodesic between the two shapes. From a mathematical point of view, shape spaces are often infinite dimensional and thus, non-trivial analytical questions arise such as existence of minimizing geodesics or global well-posedness of the initial value problem associated with geodesics. A finite dimensional approximation is still possible such as in [START_REF] Vialard | Shape Splines and Stochastic Shape Evolutions: A Second Order Point of View[END_REF], in which spline interpolation is proposed for a diffeomorphic group action on a finite dimensional manifold. It has been extended for invariant higher-order lagrangians in [START_REF] Gay-Balmaz | Invariant Higher-Order Variational Problems[END_REF][START_REF] Gay-Balmaz | Invariant Higher-Order Variational Problems II[END_REF] on a group, still finite dimensional. A numerical implementation of the variational and shooting splines has been developed in [START_REF] Singh | Splines for diffeomorphisms[END_REF] with applications to medical imaging. The question of existence of an extremum is not addressed in these publications. An attempt is given in [START_REF] Tahraoui | Riemannian cubics on the group of diffeomorphisms and the Fisher-Rao metric[END_REF] where the exact relaxation of the problem is computed in the case of the group of diffeomorphisms of the unit interval. In a similar direction, in [START_REF] Heeren | Variational time discretization of Riemannian splines[END_REF], the authors discuss the convergence of the discretization of cubic splines in some particular infinite dimensional Riemannian context on the space of shapes.

As a shape space, we are interested in this article in probability measures endowed with the Wasserstein metric. Since the Wasserstein metric shares some similarities with a Riemannian metric on this space of probability densities, it is natural to study further higher-order models in this context. Our motivation is to answer the following practical question of the extension of cubic splines to the Wasserstein space and their numerical computation.

We present in Section 2 the notion of cubic splines on a Riemannian manifold and detail its variational formulation in Hamiltonian coordinates. We then discuss independently in Section 3 a geometric approach to the Wasserstein space that will be useful for the introduction of our proposed method detailed in Section 4. Finally in Sections 5 we present the numerical entropic relaxation method and an alternative numerical method based on semi-discrete optimal transport. The reader not interested in geometric interpretation can skip directly to Section 4.

To the best of our knowledge, this question has not been yet addressed in the literature on optimal transport until very recently in two independant and simultaneous preprints : [START_REF] Tryphon | Measure-valued spline curves: An optimal transport viewpoint[END_REF] and [START_REF] Vialard | Second order models for optimal transport and cubic splines on the wasserstein space[END_REF] (this paper). Both work share the same idea of relaxing the cubic spline formulation in the space of measure using multi-marginal optimal transport. Our paper however explores a larger hierarchy of models and several numerical methods.

Cubic splines on Riemannian manifolds

In this section, we present Riemannian cubics, which are the extension of variational splines to a Riemannian manifold (M, g) where g is the Riemannian metric. Variational cubic splines on a Riemannian manifold are the minimizers of the acceleration; that is, denoting D Dt the covariant derivative, minimization on the set of curves x : [0, T ] → M of the functional (2.1)

E(x) = 1 0 g(x) D Dt ẋ, D Dt ẋ dt ,
subject to constraints on the path such as constraints on the tangent space, (x(t i ), ẋ(t i )) are prescribed for a collection of times t i ∈ [0, 1], or constraints on the positions such as x(t i ) = x i . Under mild conditions on the constraints, if M is complete, minimizers exist, for instance in the case of constraints on the tangent space mentioned above. A pathological case where minimizers might not exist is when the initial speed is not prescribed. Consider for instance the two dimensional torus, where lines of irrational slopes are dense, it is possible to show that for any collection of points which do not lie on a line, the infimum of E is 0 while it is never reached, see [START_REF] Heeren | Variational time discretization of Riemannian splines[END_REF]. The Euler-Lagrange equation associated to the functional E is (2.2)

D 3 Dt 3 ẋ -R ẋ, D Dt ẋ ẋ = 0 ,
where R is the curvature tensor of the Riemannian manifold M . Note that this equation is similar to a Jacobi field equation. We now formulate the variational problem in coordinates. In a coordinate chart around a point x(t) ∈ M , the geodesic equations are given by

(2.3) D Dt ẋ = ẍ + Γ(x)( ẋ, ẋ) = 0 ,
where Γ is a short notation for the Christoffel symbols associated with the Levi-Civita connection.

It is a second-order differential equation which is conveniently written as a first-order differential equation, via the Hamiltonian formulation. Again in local coordinates on T * M the cotangent bundle of M , the geodesic equation can be written as

(2.4) ṗ + ∂ x H = 0 ẋ -∂ p H = 0 ,
where H(x, p) = 1 2 g(x) -1 (p, p). Note that, the ODE (2.3) can be obtained from the Hamiltonian system using ẋ = g(x) -1 p. From these two equivalent formulations (2.3) and (2.4), it can be shown that g -1 (x)( ṗ + ∂ x H) = D Dt ẋ. Therefore, it proves that the variational spline problem can be rewritten in Hamiltonian coordinates as follows

inf u 1 0 g(x) -1 (a, a) dt , under the constraint ẋ -g(x) -1 p = 0 ṗ + ∂ x H(x, p) = a ,
with initial conditions x(0) = x 0 and p(0) = p 0 . It is natural to ask whether such variational problems carry over in infinite dimensional situations such as the Wasserstein space, which will be discussed in the rest of the paper.

A formal application of spline interpolation to the Wasserstein space

It is well known that the Hamiltonian formulation of geodesics on the Wasserstein space, define over a riemannian manifold M , are (3.1) ρ + ∇ • (ρ∇φ) = 0 φ + 1 2 |∇φ| 2 = 0 , where ρ : M → R ≥0 and φ : M → R implicitly time dependant are respectively a probability density and a function. Note that these equations are valid when working with smooth densities. The Hamiltonian is the following,

(3.2) H(ρ, φ) = 1 2 M |∇φ| 2 ρ dµ 0 ,
where µ 0 is a reference measure on M .

Remark 1. Taking the gradient of the equation governing φ, and denoting v = ∇φ, we get Burger's equation:

(3.3) v + (v, ∇)v = 0 ,
where in coordinates, the operator (v, ∇) is defined as (v, ∇)w . = n i=1 v i ∇w i where v, w are vector fields and n is the dimension of the M . In Lagrangian coordinates, this equation implies that (3.4) φ = 0 , where ϕ(t) : M → M is the Lagrangian flow associated with v ( φ = v • ϕ), which is well-defined under sufficient regularity conditions.

Remark 2. For the Wasserstein case, the operator is given by g(ρ) -1 φ = -∇ • [ρ∇φ] so that the (formal) computation of the covariant derivative D Dt ρ on the Wasserstein space is:

(3.5) D Dt ρ = -∇ • [ρ (v + (v, ∇)v)] ,
where v = ∇φ is the horizontal lift associated with ρ, that is ρ + ∇ • (ρ∇φ) = 0. This result is proven rigorously in [START_REF] Lott | Some geometric calculations on Wasserstein space[END_REF].

From a control viewpoint, we aim at minimizing 1 2 1 0 H(ρ, a) dt for the control system:

(3.6) ρ + ∇ • (ρ∇φ) = 0 φ + 1 2 |∇φ| 2 = a ,
where a is a time dependent function defined on M . Alternatively, in terms of the variables (ρ, φ), this amounts to minimize (3.7)

1 0 M |∇[ φ + 1 2 |∇φ| 2 ]| 2 ρ dµ 0 dt ,
under the continuity equation constraint ρ + ∇ • (ρ∇φ) = 0. It is a nonconvex optimization problem in the couple (ρ, φ). The key issue here is that the variational problem itself is a priori not well-posed since our formulation is valid in a smooth setting and to make it rigorous on the space of measures, the tight relaxation of this problem is needed. However, we do not address this issue in our work and in the next section we turn our attention to a simple relaxation of the problem which is probably not tight.

A hierarchy of relaxed models

4.1. Context. We recall the classical optimal transport setting. We have the following well known equivalence [START_REF] Otto | The geometry of dissipative evolution equations: The porous medium equation[END_REF][START_REF] Villani | Optimal transport: old and new[END_REF] (4.1)

W 2 2 (ρ 0 , ρ 1 ) = inf ϕ 1 0 M | φ| 2 dµ 0 dt = inf ρ,v 1 0 M |v| 2 dρ dt = inf ρ 1 0 inf v M |v| 2 dρ dt = inf ρ,∇φ 1 0 M |∇φ| 2 dρ dt Under constraints that [ϕ(t)] * µ 0 = ρ(t) for t = 0, 1 ([ϕ(t)] * µ 0 is the image measure of µ 0 : M f (y) d[ϕ(t)] * µ 0 (y) = f (T (x)) dµ(x) for every measurable function f : M → R ) and the continuity equation ρ + ∇ • (ρv) = ρ + ∇ • (ρ∇φ) = 0
with fixed initial and final conditions ρ(0) = ρ 0 and ρ(1) = ρ 1 .

Moreover, geodesics in the space of densities for the Wasserstein metric are given by [ϕ(t)] * µ 0 = ρ(t) and the associated displacement maps satisfy v • ϕ = φ.

The last equality in (4.1) exactly says that the infimum inf v(t) M |v(t)| 2 dρ(t) among all v(t) satisfying the continuity equation at each time t is achieved when v(t) is a gradient. This property is a consequence of a Riemannian submersion and ∇φ is called the horizontal lift of ρ. It is this last formulation that formally gives a Riemannian structure on the space of probability measures. See the remark 1 below for more details on the geometrical structure.

For higher-order variational problems, e.g. the minimization of the acceleration, the reduction in the last inequality does not holds true in general, even if the Riemannian submersion structure is present as shown in [START_REF] Gay-Balmaz | Invariant Higher-Order Variational Problems II[END_REF]. It means in the case of acceleration that, a priori, with the same constraint as for (4.1) :

(4.2) inf ϕ 1 0 M | φ| 2 dµ 0 dt = inf ρ,v 1 0 M | v + (v, ∇)v| 2 dρ dt = inf ρ,∇φ 1 0 M | φ + (∇φ, ∇)∇φ| 2 dρ dt,
where we have used that φ

= v • ϕ + (v • ϕ, ∇)v • ϕ.
Remark 1. From a geometrical point of view, (4.1) says the Wasserstein space can be seen, at least formally, as a homogeneous space as described in [15, Appendix 5] and originally in [START_REF] Otto | The geometry of dissipative evolution equations: The porous medium equation[END_REF]. Consider the group of (smooth) diffeomorphisms of M a closed manifold, Diff(M ), and the space of (smooth) probability densities Dens(M ). The space of densities is endowed with a Diff(M ) action defined by the pushforward, that is to a given ϕ ∈ Diff(M ) and ρ ∈ Dens(M ), the pushforward of ρ by ϕ is Jac(ϕ -1 )ρ • ϕ -1 . By Moser's lemma, this action is transitive, thus making the space of densities as a homogeneous space. More importantly, there exists a compatible Riemannian structure between Diff(M ) and Dens(M ). Once having chosen a reference density µ 0 , the L 2 (M, µ 0 ) metric on the diffeomorphism group descends to the Wasserstein L 2 metric on the space of densities, or in other words, the pushforward action ϕ → ϕ * µ 0 is a Riemannian submersion. An important property of Riemannian submersion is that geodesics on Dens(M ) are in correspondence with geodesics on the group, given by horizontal lift. This property is actually contained in Brenier's polar factorization theorem, which shows that the horizontal lift is the gradient of a convex function.

The Monge formulation.

In Section 3 we used the formal Riemannian structure on the set of probability measure to define an intrinsic notion of splines, (3.7) is indeed the RHS of inequality (4.2). In this section we propose a simpler alternative definition of Wasserstein splines based on the LHS of inequality (4.2).

Definition 1 (Monge formulation). Let 0 = t 0 < . . . < t n = 1, n ≥ 2 and ρ 1 , . . . , ρ n be n probability measures on M . Minimize, among time dependent maps ϕ(t) : M → M , (4.3)

1 0 M | φ| 2 dµ 0 dt , under the marginal constraints ϕ(t i ) * µ 0 = ρ i . This minimizing problem is denoted by (M S).
It is a Monge formulation of the variational problem, similar to standard optimal transport. On a Riemannian manifold M , the notation φ stands for D Dt φ. By the change of variable with the map ϕ, the problem can be written in Eulerian coordinates, that is using the vector field associated with the Lagrangian map ϕ, ∂ t ϕ = v • ϕ, one aims at minimizing for (ρ, u)

(4.4) 1 0 M |u| 2 ρ dµ 0 dt under the constraints (4.5) ρ + div(ρv) = 0 v + (v, ∇)v = u ,
with the marginals constraints ρ(t i ) = ρ i .

Remark 2. Remark that formally when v = ∇φ, this new model reduces to the formulation (3.7). Therefore, it justifies the fact that Problem (4.3) is a relaxation of (3.7). However, as already mentioned, this relaxation is probably not tight. Another formal geometric argument in the direction of proving that the two formulations are different is that the Wasserstein space has nonnegative curvature if the underlying space M has nonnegative curvature, but the space of maps in the Euclidean space is flat. Therefore, the two Euler-Lagrange equations (2.2) lead to a different evolution equations: for instance, if M is the Euclidean space then the Euler-Lagrange equation for the second model is simply .... ϕ = 0, which is a priori different from the splines Euler-Lagrange equation in the Wasserstein case.

4.3.

The Kantorovich relaxation. Since, as is well-known in standard optimal transport, the Monge formulation is not well-posed for general given margins ρ 1 , . . . , ρ n , we propose instead to solve yet another relaxation of the problem on the space of curves which takes the form: Definition 2 (Kantorovich relaxation). Let 0 = t 1 < . . . < t n = 1, n ≥ 3 and ρ 1 , . . . , ρ n be n probability measures on M .

Minimize on the space of probability measures on the path space

H 2 ([0, 1], M ) denoted by H in short, (4.6) min µ H |ẍ| 2 dµ(x) ,
which is a linear functional of dµ. The curves of densities is given by its marginals in time

(4.7) t → ρ(t)µ 0 := [e t ] * (µ) , e t is the evaluation function at time t : if γ ∈ H 2 ([0, 1], M ) ⊂ C 0 ([0, 1], M ) then e t (γ) = γ(t, .) ∈ M .
The notation [e t ] * µ is the image measure by the map e t defined by duality :

M f (y) d[e t ] * µ(y) = H f (e t (x)) dµ(x) for every measurable function f : M → R. Note that x is a path on [0, 1] × M while y is a point on M .
With these notations, the marginal constraint at given time t i are

(4.8) [e t ] * (µ) = ρ i µ 0 .
By standard arguments, the Kantorovich relaxation admits minimizers under general hypothesis on the manifold M , which we do not detail here. It is straightforward to check that existence of minimizers holds when M = R d .

As expected, the Kantorovich formulation is the relaxation of the Monge formulation in Definition 1.

Theorem 1. Let M = R d , 0 = t 1 < . . . < t n = 1, n ≥ 3 and ρ 1 , . . . , ρ n ∈ be n probability measures on R d with compact support and ρ 1 being atomless. Then, under the constraints (4.8), the infimums of the variational problem (4.3) and (4.6) coincide, moreover, the infimum is attained for the latter.

Proof. See the proof of a more general result in Appendix A.

First we remark that we can reformulate both the Monge and Kantorovich problems on the set of cubic splines. It is the purpose of the following lemmas and corollaries, whose proofs are straightforward.

Definition 3 (Cubic interpolant). Let (x 1 , . . . , x n ) ∈ R d be n given points and (t 1 < . . . < t n ) be n timepoints. There exists a unique cubic spline minimizing the acceleration of the curve x(t) such that x(t i ) = x i . This unique curve is called cubic interpolant and is denoted by c x1,...,xn , depending implicitly on the timepoints.

Lemma 2. When the supports of the measures ρ i are compact on R d , the support of every minimizing µ in Definition 2 is included in the set the cubic interpolants c x1,...,xn for (x 1 , . . . , x n ) ∈ Supp(ρ 1 ) × . . . × Supp(ρ n ).

Proof. The constraints are the marginal constraints [e ti ] * (µ) = ρ i for i ≥ 3 which implies that set of paths charged by an optimal measures satisfies x(t i ) ∈ Supp(ρ i ). In particular, any path in this set can be replaced by its minimal spline energy, the cubic interpolant c x1,...,xn .

Corollary 3. As a consequence, the set of paths charged by an optimal plan are uniformly C 2 and for every smooth function η : R d → R with compact support, the map t → µ(t), η is C 2 .

Proof. The set of cubic interpolants is compact since the map (x 1 , . . . , x n ) → c x1,...,xn is continuous from R dn to the space of C 2 fonctions (solution of an invertible linear system) and Supp(ρ i ) are compact. Therefore, the set of maps are uniformly C 1 . The last point follows directly.

Corollary 4. The Kantorovich problem in Definition 2 on R d reduces to a multimarginal optimal transport problem, as follows, let c(x 1 , . . . , x n ) be the continuous cost of the cubic interpolant at times t 1 , . . . , t n , the minimization of (4.6) reduces to the minimization of (4.9)

M n c(x 1 , . . . , x n ) dπ(x 1 , . . . , x n ) (K)
on the space of probability measures π ∈ P(M n ) and under the marginal constraints (p i ) * (π) = ρ i where p i is the projection of the i th factor.

Proof. Direct consequence of Lemma 2.

Similarly

Corollary 5. The Monge problem in Definition 1 on R d reduces to a Monge multimarginal optimal transport problem, as follows, let c(x 1 , . . . , x n ) be the continuous cost of the cubic interpolant at times t 1 , . . . , t n , the minimization of (4.3) reduces to the minimization of

(4.10) M c (x, ϕ(t 1 , x), . . . , ϕ(t n , x)) dµ 0 (x),
on the space of path ϕ ∈ C 2 ([0, 1], M ) (or even cubic splines) and under the marginal constraints

(ϕ(t i )) * µ 0 = ρ i .
The dual formulation of the minimization problem (K) is also well known [16, Theorem 2.1]

Definition 4 (Kantorovich dual problem (KP )). Let Q = φ i ∈ L 1 (ρ i µ 0 ) , i = 1..n be the space of integrable n-uplet. Maximize on Q (4.11) n i=1 M φ i ρ i µ 0 , under the constraint n i=1 φ i (x i ) ≤ c(x 1 , ..., x n ).
And the following duality results holds true: Proposition 6. There exists a n-uplet (φ i ) i=1..n ∈ Q optimal for (KP ). Moreover (K)=(KP ) and for any π optimal in (4.9) there holds n 1 φ i (x i ) = c(x 1 , ..., x n ), π almost everywhere. A natural question is whether the solution of the Kantorovich problem (K) is admissible in the Monge formulation (M S) (Definition 1). With the formulation reduced above the spline, given by (4.9) and (4.10), one can try to apply existing theory to answer to this question, see [START_REF] Kim | A general condition for monge solutions in the multi-marginal optimal transport problem[END_REF][START_REF] Pass | Multi-marginal optimal transport: Theory and applicationsâˆ[END_REF] and references therein for precise criterion. However our cost does not satisfy any of those known criterion. In fact, we have the following result which proves that the relaxation to plans are necessary even in the context of Theorem 1. Proposition 7. (Counter Example) Given the three-marginals problems of minimizing the acceleration, there exist data (ρ 0 , ρ 1 , ρ 2 ) such that ρ 0 is atomless and such that the solution of (K) is not a (measurable) Monge map.

Proof. Consider ρ 0 (x) = 1 [-1,1] and the Dirac masses a = δ 1 and b = δ -1 and the maps T a , T b that respectively pushforward ρ 0 onto a and b. These maps are uniquely determined and affine. Consider now

ρ 2 = 1 2 (T a ) * ρ 0 + 1 2 (T b ) * ρ 0 = a 2 + b 2 . Then, introducing (T 1/2 ) = 1 2 (Id +T ), we consider ρ 1 = 1 2 (T 1/2 a ) * ρ 0 + 1 2 (T 1/2 b ) * ρ 0 , note that it is equal to ρ 0 since the maps T 1/2
a,b are affine. By construction, the minimization of the acceleration for (ρ 0 , ρ 1 , ρ 2 ) is null since it is a mixture of plans supported by straight lines. If there existed an optimal Monge solution it is necessarily supported by only one map denoted by T and since the cost is null, the map at time 1/2 is necessarily T 1/2 defined above. The preimage of 1 (resp. -1) by T is a measurable set A (resp. B). Then, necessarily, ρ 1 = (T 1/2 ) * χ A + (T 1/2 ) * χ B , and in fact, T |A = T a and T |B = T b (since the image of the map is known). Therefore, we have

ρ 1 = 2χ A • (T 1/2 a ) -1 + 2χ B • (T 1/2 b ) -1 which is not equal to the uniform Lebesgue measure on [-1, 1]. Remark 3.
It is an open question to prove or disprove a similar result when the final density ρ 2 is atomless. The counterexample explained above strongly uses the fact that the final density is a sum of Dirac masses and it might not be robust when replacing the final density by a uniform density on a small interval.

4.4.

The corresponding interpolation problem on the tangent space. The relaxed problem on the space of curves can be used to define variational interpolation problem on the phase space, or more precisely on the tangent space T M . Since the space H 2 ([0, T ], M ) is contained in C 1 ([0, T ], M ), one can formulate the optimal transport problem on phase space (identified with the tangent space) for the acceleration cost.

Definition 5 (Optimal transport on phase space). Let ρ0 , ρ1 be two probability measures on T M . Minimize on the space of probability measures on H, where

j t : H 2 ([0, T ], M ) → T M is defined by j t (x) = (x(t), ẋ(t)). time R • -1 • 1 • t = 2 • t = 1 • t = 0 Figure 1.
The inital density at time 0 is described with a mixture of two densities colored in red and blue which are evolving indepently along straight lines in time. The blue density is mapped onto -1 and the red density is mapped onto 1. The acceleration cost is null and the proof of Proposition 7 shows that it is not possible to reproduce the density at time 1/2 by a map.

Proposition 8 (Optimal interpolation on phase space). The support of every optimal solution is contained in the set of cubic splines interpolating between (x, v) ∈ Supp( ρ0 ) and (y, w) ∈ Supp( ρ1 ). Moreover if M = R d and if ρ0 has density with respect to the Lebesgue measure, then the unique solution to Problem (4.12) is characterized by a map ϕ : T M → T M .

Remark that the optimal solution in the last part of Proposition 8 provides an interpolation on the phase space using [j t ] * (µ).

Proof. The proof of the first part is similar to Lemma 2 and the second part follows by application of Brenier's theorem since the total cost of the cubic splines between (x, v) and (y, w) can be explicitly computed as (4.14) c ph ((x, v), (y, w)) = 12|x -y| 2 + 4(|v| 2 + |w| 2 + v, w + 3 v + w, x -y ) and satisfies the twisted condition, so [30, Theorem 10.28] applies.

Note that this problem is very different from using the Wasserstein distance on P(T M ) where the tangent space T M is endowed with the direct product metric. Indeed, the cost c ph does not vanish on the diagonal (x, v) = (y, v) contrarily to the quadratic cost on T M .

Interestingly, let us remark that the multimarginal problem can be recast as the minimization problem on Π ∈ P(T M × . . . × T M n times

), denoting Π ti,ti+1 the pushforward on T M × T M at times

(t i , t i+1 ), (4.15) min π n-1 i=1 Π ti,ti+1 , c ph ((x i , v i ), (x i+1 , v i+1 ))
under the constraints that [e ti ] * (Π i,i+1 ) = ρ i . From the numerical point of view, this rewriting might be useful since the cost used on the multimarginal problem is now separable in time. This relaxation to the tangent space is used in the semidiscrete algorithm in Section 5.3.1. Obviously, up to the minimization on the variables v i , we retrieve the minimization problem (K) since one has a cost c which is defined on M n

c ph ((x i , v i ), (x i+1 , v i+1 ))
where the index i runs over the marginals.

Numerical Study

We have discussed several variational relaxation of the classical definition of splines, applied to the Wasserstein space of densities. At least two different numerical techniques from Optimal Transportation can be used in this setting. We apply the Entropic regularisation and Sinkhorn (briefly recalled in appendix B first to a simple Hermite interpolation problem (section 5.1) and then in section to the multimarginal problem (4.9). In section 5.3, we use the semi-discrete Optimal Transportation approach in the spirit of [START_REF] Mérigot | Minimal geodesics along volume-preserving maps, through semidiscrete optimal transport[END_REF] directly to problem (4.6) without the time discretisation in (4.9).

Hermite interpolation.

In this section, we are interested in the problem of interpolation on the phase space described in the previous. The marginals [e t ] * (µ) are densities defined on the tangent space T M . If we only specify the marginals at time 0 and 1 as empirical measures: [e 0 ] * (µ) = k i=1 α i δ xi δ vi and [e 1 ] * (µ) = k j=1 β j δ yj δ wj , as explained in Section 4.4, we can simplify the Kantorovich using the exact L 2 norm of the acceleration of the spline between (x v ) and (y, w), whose cost is given in Formula (4.14). Again, let us underline that this cost is not a Riemannian cost on the tangent space of R d since if v = w and x, y are close, the cost is dominated by the term 4(|v| 2 + |w| 2 + v, w ) which need not be zero. Then, the Kantorovich problem reduces to the minimization of (5.1)

k,l i,j=1 π i,j c((x i , v i ), (y j , w j )) , under the constraints (5.2) k i=1 π i,j = β j l j=1 π i,j = α i .
It is straightforward to apply entropic regularization/Sinkhorn in this case which amounts to add, for a positive parameter ε, ε i,j π i,j log(π i,j ) to the previous linear functional and to numerically solve the corresponding variational problem with the Sinkhorn algorithm [START_REF] Sinkhorn | Diagonal equivalence to matrices with prescribed row and column sums[END_REF][START_REF] Cuturi | Sinkhorn distances: Lightspeed computation of optimal transport[END_REF] (See also appendix B where Sinkhorn algorithm is detailed in the more general multimarginal case). It is interesting to note that the choice of ε is more delicate than in the standard case of a quadratic distance cost.

In Figure 2, we present the convergence rate of this method with respect to two different values of ε and the most likely deterministic plan given the optimal plan π ε . Note that this entropic regularization method scales with the number of points as N 2 and is valid in every dimension.

5.2.

MultiMarginal formulation. This is the direct discretization of (4.6) which avoids working in phase space with the cost (4.16) thus enabling fast computations in 2D. In what follows, the time cylinder [0, 1] × M is discretized in time as i=0,N M i , the product space of N + 1 copies of M at each of the N + 1 time steps. We will use a regular time step discretization τ i = i dτ where dτ = 1 N . Using a classic finite difference approach, the time discretization of (4.6) is (5.3) min

µ dτ i=0,N Mi c dτ (x 1 , ..., x N ) dµ dτ (x 1 , .., x N ) ,
where µ dτ now spans the space of probability measures on i=0,N M i representing the space of piecewise linear curves passing through x 0 , x 1 , ..., x N at times τ 0 , ..., τ N .

A straighforward computation gives (5.4) c dτ (x 1 , ..., x N ) := i=1,N -1 For all times, marginals (4.7) are computed as :

x i+1 + x i-1 -2 x i 2 dτ 3
(5.5)

τ j → i =j Mi dµ dτ (x 1 , ..x N )
In order to simplify the presentation we will assume that the marginal constraints (4.8) are set at times t 1 , ..t n which coincide with times steps of the discretization (of course n < N , meaning the number of constraint is not the same as the number of time steps).

In short, there exist (j 1 , ..j n ) ∈ [0, N ] such that (t 1 , .., t n ) = (τ j1 , ..., τ jn ).

The constraint (4.8) becomes for all k = 1, ..n

(5.6) i =j k Mi dµ dτ (x 1 , ..x N ) = ρ j k (x j k )
where ρ j k is the prescribed density to interpolate at time τ j k = t k .

The time discretized problem is the multimarginal problem (5.3 -5.6).

The simplest space discretization strategy is to use a regular cartesian grid. In dimension 2 and for M = [0, 1] 2 and at time t i , the grid will be denoted

x αi,βi = (α i h, β i h) for (α i , β i ) ∈ [0, N x ] and h = 1
Nx , a = {α i } and b = {β i } will be the vectors of indices. The time and space discretization of the problem then becomes

(5.7) min T a,b C a,b T a,b
Where T is the N × N x × N x tensor of grid values µ dτ (x α1,β1 , .., x α N ,β N ) and

(5.8) C a,b = c dτ (x α1,β1 , .., x α N ,β N )
The marginals (5.5) at all times τ j are given by (5.9)

a\{αj }, b\{βj } T a,b
The constraints (5.6) therfore becomes for all k

(5.10)

a\{αj k }, b\{βj k } T a,b = ρ j k (x αj k ,βj k ) a \ {α j k } denotes the set of indices a minus α j k .
The Entropic regularized problem is Numerical Simulations. 1D case: We present, figures 3 and 4, a 1D test case to highlight some of the qualitative properties of the cubic splines interpolation on the space of densities. We consider four interpolation time points and the corresponding data are mixture of Gaussians of different standard deviations. We use a discretization of 140 points on the interval [0, 1] with 16 time steps. The doted line represent the reconstructed density curve in time. This experiment shows that the mass can concentrate or diffuse in some situation.

Another important point here is that the entropic regularization parameter has an important impact on this concentration/diffusion effects: we show the simulations for ε = 0.002 and ε = 8.10 -5 . In the simulation with a large ε, the concentration effect is not present and it is due to the diffusion on the path space. 2D case: We present a 2D test case which computes a Wasserstein spline in the sense of (5.7) interpolating four Gaussian identical densities at time 1, 5, 13, and 17, see figure 5. We use a time step dτ = 1 and 17 N = 17 time steps. The space discretization is N x = 50. The entropic regularization parameter is = 0.002, note that the stability of the method depends on this parameter. It also generates artificial diffusion as it becomes more costly top concentrate the available mass on fewer Euclidean splines between the points of the support of the four Gaussians. We can compute the interpolating densities at intermediate times using (5.9) but is more interesting to represent in figure 6 the contour line of the third quartile, i.e. the highest values of the densities representing 1/4 of the total mass. Comparing with figure 7, it seems clear that the Entropy diffusion spreading pollutes the solution of the original problem (without entropic regularization).

We compare this solution with the classical Quadratic cost Optimal Transport interpolation, i.e. with the speed instead of the acceleration in the cost. More precisely taking :

(5.12)

c dτ (x 1 , ..., x N ) := i=0,N -1 x i+1 -x i 2 dτ
As expected the mass follows respectively the linear interpolation or the Euclidean spline interpolation of the center of the Gaussians which are represented as thick red lines in figure 5.

Finally we show the convergence of the Sinkhorn iterate for both simulations in figure 6. The convergence is much slower for the speed case but we did not optimize the implementation which does not need tensors and instead just used a degraded version of the acceleration code. This may be the reason for this strange difference. 5.3. Semi-Discrete approach. We propose another numerical scheme based on the semi-discrete approach introduced by Mérigot in [START_REF] Mérigot | A multiscale approach to optimal transport[END_REF] in dimension 2 and developed by Levy [START_REF] Lévy | A numerical algorithm for l2 semi-discrete optimal transport in 3d[END_REF] in dimension 3.

Here we approximate the optimal plan π in the formulation (4.9) by a sum of N tensor product of diracs masses. That is Remark 4. Since there is a unique corresponds between n points X 1 j , . . . , X n j and the spline c X 1 j ,...,X n j passing through these points at time (t 1 , . . . .t n ) the measure π N can also be seen as N direct masses defined over the set of splines:

π N = N j=1 n i=1 1 N δ X i j = N j=1 1 N δ (X 1 j ,...,X n j ) .

Initial time data and targets

π N = N j=1 1 N δ c X 1 j ,...,X n j .
We then have to relax the constraint

(p i ) * (π) = ρ i since (p i ) * (π N ) = N j=1
1 N δ X i j cannot be absolutely continuous. It leads to the following variational problem. Definition 6 (Semi-discrete variational problem). Let > 0, 0 = t 1 < . . . < t n = 1, n ≥ 3 and (ρ i ) i=1...n be n absolutely continuous measures. Recall that c(Y 1 , . . . , Y n ) is the cost of the cubic spline passing through the points (Y 1 , . . . , Y n ) at time (t 1 , . . . .t n ). Let Then the semi-discrete variational problem, (SDV), is given by (5.13) (SDV ) = min where W 2 is the classical Wasserstein distance given by the quadratic cost.

Q N =    N j=1 1 N δ (X 1 j ,...,X n j ) (X j ) j=1,...,N ∈ M n    .
Q N 1 N N j=1 c(X 1 j , . . . , X n j ) + n i=1 1 2 2 W 2 2   N j=1 1 N δ X i j , ρ i   ,
The main drawback of this method is that, as illustrated in the numerical simulations below, the problem (SDV ) is not convex. 5.3.1. Implementation. In order to solve numerically the minimization problem (SDV ) we use the reformulation of the spline cost in the phase space, that is in R d , with t i+1 -t i = δ i :

(5.14) c (Y 1 , . . . , Y n ) = min (V1,..Vn)∈(R d ) n n-1 i=1 1 δ 3 i c ph [(Y i , δ i V i ) , (Y i+1 , δ i V i+1 )]
where

(5.15) c ph [(x, v), (y, w)] = 12|x -y| 2 + 4(|v| 2 + |w| 2 + v, w + 3 v + w, x -y ).
The advantage of the formulation (5.14) is that the cost is separable in the phase space and the gradient with respect to speeds and positions is easy to compute. We thus implement a gradient descent in the phase space using the lbfgs function in python. We compute the gradient by automatic differentiation. The Wasserstein terms in the minimization problem (5.13) depends only on the positions and are computed thanks to Mérigot Library [1] in dimension 2. To do simulations in dimension 3 one has to use Lévy Library [2]. The density constraints ρ i are given trough linear functions on a triangulation. Remark 5. Other problems can be addressed using similar optimization problem as in Definition 6. For instance the quadratic cost in (5.13) leads to Wasserstein interpolation. We can also interpolate with curves as smooth as we want, using for instance the L 2 norm of the derivative of order m of the curve or even other classical interpolating curves. 5.3.2. Numerical simulations. We propose three numerical simulations, one to compare the qualitative results with respect to the multi marginal approach and especially Figure 5. A second one in order to illustrate the non-convexity issue and a third one for applications in images. The rotation case: Figure 7. In this case we compute Wasserstein splines passing through four gaussians with variance 15 and center of masses respectively (0, 2), (10, 0), [START_REF] Galichon | Matching with Trade-offs: Revealed Preferences over Competiting Characteristics[END_REF][START_REF] Camarinha | Splines of class C k on non-euclidean spaces[END_REF], (0, 4) with constraint parameter = 10 -3 . The number of points is 2000. In this case the result is a global minimizer and is not sensible to initialization. The lack of convexity is not an issue. Compare to . Left: sample of each density constraints ρ i , i=1,2,3,4. Right: Some trajectory of diracs masses randomly chosen, marginals at the constrained time 0, 1, 2, 3 and marginals at time 0.5, 1.2, 1.5, 1.7, 2.5. Second Line : the same configuration as in figure 5.

Figure 5, this approach gives a better a approximation of the intermediate densities especially with less diffusion. The crossing case: Figure 8,9. Here we compute Wasserstein splines starting from a mixture of two gaussians with centrer (0, -1), (0, 1) and variance 15 then passing through a gaussian with center (0, 0) and variance 15 and finishing at a translation of the initial mixture. The number of points is 2000, will value 1 or 1000.

We expect the global minimizer to be straight lines crossing around the middle constraint and with a low cost. Numerically depending on the initial conditions, we can recover different local minimizers, the local minimum which is reached is extremely correlated with the initial coupling. In Figure 8 we observe that changing but keeping a similar initial coupling, all points are given by a quantization of the middle density with a random enumeration and 0 initial speed, yields to a similar local minimum.

Finding a good initial coupling is the hard part in order reach the global maximum. One solution is to initialize with points close to each other and a very large . Then one as to add some noise in the gradient and decreases slowly . Unfortunately we didn't find a systematic approach for this random multi-scale method and one as to fit the parameters case by case. In Figure 9 the global minimizer is achieved by first computing the spline with a relaxed constraint, i.e. large , only for the final time ( in pratice = [1000, 1000, 1]. Then we use this result, which has the good initial coupling, as and initial condition and set = 1000 for all the constraints. We also compare this results with the interpolation with a different initial condition and the Wasserstein geodesics. In all these simulations we clearly observe that particles can cross along the dynamic appart from the optimal transport inthis situation.

Note that this spline approach is related to the problem of finding minimal geodesics along volume preserving maps done by Mérigot and Mirebeau [START_REF] Mérigot | Minimal geodesics along volume preserving maps, through semidiscrete optimal transport[END_REF] : in their work the constraints ρ i are the Lebesgue measure, the cost is changed by the quadratic cost between two points and they have a coupling constraint. Therefore their minimization problem is also non convex but the coupling is given as a constraint so the non convexity issue didn't rise as clearly as in this spline problem. Image interpolation: pour l'instant c'est pas presentable, ca passe vraiment au milieu. Je vais relancer dans la semaine mais je propose de faire une version sans.

Remark 6 (Extrapolation). The minimization of the acceleration can be used to provide time extrapolation of Wasserstein geodesic in a natural way: particles follow straight lines. This can be implemented in a 3-marginal problem with the acceleration cost c(x 1 , x 2 , x

3 ) = 1 λ 2 |x 3 -2x 2 + x 1 | 2 + 1 λ |x 2 -x 1 | 2
under marginal constraints at time 1 and 2. Note that, in the spline model, the formulation we proposed does not prevent particles from crossing each other. They are completely independent. Therefore, the particles following simply geodesic lines and after a shock, the evolution is not geodesic in the Wasserstein sense (since shocks do not occur but at initial and final times). The implementation of time extrapolation using entropic regularization is straightforward. Figures 10 and11 show some experiments on [0, 1] discretized with 100 points and ε = 0.015. The translation experiment recovers what is expected however the effect of the diffusion can be seen with a twice larger ε. We also show two other simulations, one is a splitting simulation and the last one is a merging of two "bumps" into a single one. The extrapolation shows an other bimodal distribution which is explained by particle crossings. Note that this extrapolation scheme may proven useful in the development of higher-order schemes for the JKO algorithm.

Perspectives

In this paper, we presented natural approaches to define cubic splines on the space of probability measures. We have presented a Monge formulation and its Kantorovich relaxation on the path space as well as their corresponding reduction on minimal cubic spline interpolation. We leave for future work theoretical questions such as the study of conditions under which the existence of a Monge map as a minimizer occurs, as well as the relaxation of cubic spline in the Wasserstein metric. Our main contributions focus on the numerical feasibility of the minimization of the acceleration on the path space with marginal constraints. We have developed the entropic regularization scheme for the

(B.3) U k j k = ρ j k (x αj k ,βj k )( a\{αj k }, b\{βj k } K a,b Π k ∈{1,..n}\k U k j k ) -1
Sinkhorn algorithm simply amounts to perform a Gauss-Seidel type iterative resolution of the system (B.3) and therefore consists in computing the sums on the right-hand side and then perform the (grid) point wise division. B.2. Implementation. In dimension 2, each unknown U k has dimension N 2

x , the cost of one full Gauss Seidel cycle, i.e. on Sinkhorn iteration on all unknowns, will therefore be n × N 2

x × the cost to compute the tensor matrix products in the denominator of (B.3). Remember that n is the number of time steps with constraints and N the total number of time steps. The given tensor Kernel K a,b is a priori a large N × N x × N x tensor with indices a, b = α 1 , ..α N , β 1 , .., β N . It can however advantageously be tensorized both along dimensions and also margins. First, using (5.4-5.8) we see that the Kernel is the product of smaller tensors K a,b = Π i=1,N -1 K 0 i-1,i,i+1 , with K 0 i-1,i,i+1 := e -1 dτ 3 x α i+1 ,β i+1 +x α i-1 ,β i-1 -2 x α i ,β i

2

.

Moreover as we chose to work on a cartesian grid at all time steps, K 0 tensorize again into K 0 i-1,i,i+1 = K α i-1,i,i+1 K β i-1,i,i+1 with K α i-1,i,i+1 := e -h 2 dτ 3 αi+1+αi-1-2αi 2

Finally our large kernel K a,b can be represented a the product of 2 (N -2) identical tensors of size N x × N x × N x . Assuming a cubic cost n 3 for the multiplication of two (n × n) matrix, we see oru algorithm is of order O(N N 4 x ) in dimension 2.

  (x) , which is a linear functional of µ under the marginal constraints [j 0 ] * (µ) = ρ0 , and [j 1 ] * (µ) = ρ1 , (4.13)

( 4 .

 4 [START_REF] Kim | A general condition for monge solutions in the multi-marginal optimal transport problem[END_REF]) c(x 0 , . . . , x n ) = min v0,...,vn n-1 i=1

Figure 2 .

 2 Figure 2. Convergence (left) and Hermite interpolation problem between Two empirical measure in phase space (right). We represent the most likely splines in the position space.

  {C a,b T a,b + T a,b log(T a,b )} and easier to solve. See Appendix B for a description of Sinkhorn algorithm.

Figure 3 .

 3 Figure 3. Four interpolation timepoints, 1, 6, 11, 16 and representation of the four density configurations, as well as 6 intermediate times. The doted line represent the reconstructed density curve in time. This experiment underlines that the spline curve has more smoothness in time and can present some concentration or diffusion effects depending on the data which would not be present for the usual Wassertein geodesic. The entropic regularization parameter is ε = 8.10 -5 .

Figure 4 .

 4 Figure 4. The same experiment with a larger entropic regularization parameter ε = 0.002. As expected, we observe less concentration of mass.

Figure 5 .

 5 Figure 5. Spline interpolation of Four Gaussians with 17 times steps. Left : the data and the linear and classic cubic spline interpolation of the of Gaussian center point. Right : the level curve of the third quartile of the density every 2 time steps, in solid line for our Spline Wasserstein interpolation and in dashed line for the classic quadratic cost (speed) interpolation.

Figure 6 .

 6 Figure 6. Convergence, i.e. Infinity norm of the difference of the Dual unknown between to Sinkhorn iteration. This is computed every 10 iterations. Left :for the acceleration cost, right : for the speed cost .

Figure 7 .

 7 Figure 7. Spline interpolation for gaussians with 2000 Dirac masses for each measure, = 1O -3. Left: sample of each density constraints ρ i , i=1,2,3,4. Right: Some trajectory of diracs masses randomly chosen, marginals at the constrained time 0, 1, 2, 3 and marginals at time 0.5, 1.2, 1.5, 1.7, 2.5. Second Line : the same configuration as in figure5.

Figure 8 .

 8 Figure 8. Spline interpolation for a mixture of gaussians with 2000 Dirac masses. Same initial coupling for both figure. Left: = 1. Right: = 1000.

Figure 9 .

 9 Figure 9. Spline interpolation for a mixture of gaussians with 2000 Dirac masses for each measure. = 1000. Top Left: Initialization with a good coupling, total cost = 302. Top Right: Initialization with a quantization of the middle density and no speed, total cost = 804 (local minima). Bottom: Interpolation with the Wasserstein geodesic. = 1000, cost = 930.

Figure 10 .K

 10 Figure 10. Extrapolation of a translation with two different ε = 0.015 and ε = 0.03

approach to the unbalanced case in the spirit of [START_REF] Chizat | An Interpolating Distance between Optimal Transport and Fisher-Rao[END_REF] seems challenging due to the this dimensionality constraint and could be achieved within the semi-discrete setting.

In the Lagrangian setting, i.e. semi-discrete method, the extrapolation of a Wasserstein geodesic between ρ 0 and ρ 1 is obtained using three positions with the following formulation : let

where d is the distance on M and c(X 1 j , X 2 j , X 3 j ) the cost of the cubic spline. In particular this formulation forces the curve to be a Wasserstein geodesic between ρ 1 and ρ 2 , using the quadratic cost, and let free the final marginal. The implementation is completely similar as in Section 5.3 and the trajectory of each dirac masses is a straight line.

Appendix A. Proof of Theorem 1

The proof is a rewriting of the proof of [START_REF] Santambrogio | Optimal transport for applied mathematicians[END_REF]Theorem 1.33] when the initial and final spaces do not have the same dimension. In particular we prove that transport plans concentrated on a graph of a map T : R d → R p are dense into transport plans in R d × R p and deduce, taking p = (n -1)d, that for any continuous cost the multimarginal Kantorovich problem is the relaxation of the multimarginal Monge problem. Theorem 9. Let M = R d and c : M n → R be a continuous cost fonction. Let (ρ i ) i∈1,...,n be n probability measures on M . We define the Monge Problem (M c ) as

, where p i is the projection of the i th factor. Then, if all (ρ i ) i∈1,...,n have compact support and ρ 1 is atomless there holds (M c ) = (K c ).

In order to prove Theorem 9 we first remark that [25, Corrollary 1.29 and Theorem 1.32 ] have their multimarginal counterpart. Lemma 10. Let µ ∈ P(R d ) be atomless measure and ν ∈ P(R p ), then there exists a transport map T : R d → R p such that T * µ = ν.

Proof of Lemma 10. Let σ d : R d → R (resp σ p : R p → R) be an injective Borel map with Borel inverse (see [START_REF] Santambrogio | Optimal transport for applied mathematicians[END_REF]Lemma 1.28] for instance for a very simple proof of existence in this case). Since µ is atomless (σ d ) * µ is also atomless. Let t : R → R be the optimal transport map from (σ d ) * µ to (σ p ) * ν for the quadratic cost.

Theorem 11. With the notation of Theorem 9, if the support of all ρ i are included in a compact domain then the set of plans Π T induced by a transport is dense, for the weak topology, in the set of plans Π whenever ρ 1 is atomless.

Remark 7. Theorem 11 is in fact very general, one can consider M N be only Polish spaces for instance. Then there exists invertible Borel maps from M (resp N) to [0, 1]. This is enough to obtain Lemma 10. Then one just need to consider a uniformly small partition of Ω to prove the density Theorem 11.

Proof of Theorem 11. Again the proof is based on [START_REF] Santambrogio | Optimal transport for applied mathematicians[END_REF]Theorem 1.32]. In particular the strategy of the proof is to approach a transport plan by transport maps defined on small sets on which the measure is preserved.

We consider a compact domain

For any m set a partition of Ω p (resp Ω q ) into (disjoint) sets K i,m (resp L j,m ) with diameter smaller than 1/2m. Then C i,j,m = K i,m × L j,m is a partition of Ω into sets with diameter smaller than 1/m. Let π i,m be the restriction of π on K i,m × Ω p and µ i,m = (p R d ) * (π i,m ) and ν i,m = (p R d ) * (π i,m ). Since µ is atomless µ i,m = µ |Ki,m is also atomless and thanks to Lemma 10 there exists t i,m such that (t i,m ) * µ i,m = ν i,m . By definition (A.1)

where t m is define on Ω by t |Ki,m = t i,m . In particular (t m ) * (µ) = ν. Equation (A.1) and the definition of the partition sets C i,j,m implies that (Id, t m ) * (µ) weakly converges toward π as m + ∞ (they give same masses to any set of the partition). See [Theorem 1.31]santambrogio2015optimal for instance. To finish the proof let us remark that we can set p = d(n -1) then µ = ρ 1 is atomless and t m : R d → R d(n-1) defines (t 2,n , ..., t n,m ).

Proof of Theorem 9 . The continuity of the cost c and the density Theorem 11 implies that (K c ) ≤ (M c ). Since the converse is always true we have (M c ) = (K c ).

Remark 8. Theorem 1 is a consequence of Theorem A since both the Monge and the Kantorovich (Definition 1 and 2) problems reduces on M n with the spline cost which is continuous (see Corollary 4 and 5.

Appendix B. Entropic Regularisation and Sinkhorn

B.1. Entropic regularization and Sinkhorn algorithm. The linear programming problems (5.7-5.10) is extremely costly to solve numerically and a natural strategy, which has received a lot of attention recently following the pionneering works of [START_REF] Galichon | Matching with Trade-offs: Revealed Preferences over Competiting Characteristics[END_REF] and [START_REF] Cuturi | Sinkhorn distances: Lightspeed computation of optimal transport[END_REF] is to approximate these problems by strictly convex ones by adding an entropic penalization. It has been used with good results on a number of multi-marginal optimal transport problems [START_REF] Benamou | Iterative Bregman projections for regularized transportation problems[END_REF] [4] [START_REF] Benamou | Generalized incompressible flows, multi-marginal transport and Sinkhorn algorithm[END_REF]. Here is a rapid and simplified description, see the references above for more details.