Amandine Aftalion
email: amandine.aftalion@uvsq.fr.

J Fr Éd

Éric Bonnans
email: frederic.bonnans@inria.fr.

OPTIMIZATION OF RUNNING STRATEGIES BASED ON ANAEROBIC ENERGY AND VARIATIONS OF VELOCITY. REVISED VERSION 08 MAY 2014

Keywords: Running race, anaerobic energy, energy recreation, optimal control, singular arc, state constraint, optimality conditions

The aim of this work is to present a model relying on a system of ordinary differential equations describing the evolution of the anaerobic energy and the velocity of a runner. We perform numerical simulations and rigorous analysis to deduce an optimal running strategy. Our model relies on the fundamental principle of dynamics, energy conservation, a hydraulic analogy, and control conditions. We take into account the resistive force, the propulsive force and the variations in the volume of oxygen used per unit time. Our main result is to show that varying one's velocity rather than running at a constant velocity allows to run longer. Using optimal control theory, we present proofs on the structure of the optimal race and relate the problem to a relaxed formulation, where the propulsive force represents a probability distribution rather than a function of time. Our mathematical analysis leads us to introduce a bound on the variations of the propulsive force to obtain a more realistic model which displays oscillations of the velocity. We also present numerical simulations of our system which qualitatively reproduce quite well physiological measurements on real runners. We show how, by optimizing over a period, we recover the oscillations of speed.

1. Introduction. In sports training, the use of various measurements devices to track speed and calories and to identify the way bodies react to efforts has greatly spread. The issue of the optimal strategy of a runner, given the distance or time to run, is still a major challenge. The motivation of this paper is to provide a reliable system of differential equations describing the evolution of the anaerobic energy and the velocity of a runner. We are going to use the optimal control theory and algorithms.

A pioneering mathematical work is that of Keller [START_REF] Joseph | Optimal velocity in a race[END_REF] relying on Newton's law of motion and energy conservation. For sufficiently long races (> 291m), his analysis leads to an optimal race in three parts 1. initial acceleration with a propulsive force at its maximal value, 2. constant speed during the major part of the race, 3. final small part with constant energy equal to zero. Though this analysis reproduces quite well the record times for distances up to 10km, it has some weaknesses. First, physiological measurements [START_REF] Billat | Differential modeling of anaerobic and aerobic metabolism in the 800-m and 1,500-m run[END_REF][START_REF] Hanon | Pacing strategy and vo2 kinetics during a 1500-m race[END_REF] show that runners do not keep a constant speed but vary their speed by an order of 10%. Next, Keller assumes that the runner keeps a constant value of V O2 max , the maximal oxygen uptake, whereas this value gradually increases to its maximal value, and then drops at the end of the race (Figure 1 of [START_REF] Billat | Differential modeling of anaerobic and aerobic metabolism in the 800-m and 1,500-m run[END_REF] and Figure 1 of [START_REF] Hanon | Pacing strategy and vo2 kinetics during a 1500-m race[END_REF]). Some authors [START_REF] Behncke | A mathematical model for the force and energetics in competitive running[END_REF][START_REF] Mathis | The effect of fatigue on running strategies[END_REF][START_REF] Woodside | The optimal strategy for running a race (a mathematical model for world records[END_REF] have tried to improve Keller's model, but still relying on the same strategy and mathematical arguments, leading to an almost constant speed. Other references concerning the optimality of a run include [START_REF] Morton | A 3-parameter critical power model[END_REF][START_REF] Morton | A new modelling approach demonstrating the inability to make up for lost time in endurance running events[END_REF][START_REF] Ward-Smith | A mathematical theory of running, based on the first law of thermodynamics, and its application to the performance of world-class athletes[END_REF].

In this paper, we will rely on Keller's equations [START_REF] Joseph | Optimal velocity in a race[END_REF] and we will try to improve them using a hydraulic analogy and physiological indications described in [START_REF] Morton | The critical power and related whole-body bioenergetic models[END_REF]. Nevertheless, the formula of [START_REF] Morton | The critical power and related whole-body bioenergetic models[END_REF] yield averaged values while we want to make instantaneous energy balance taking into account optimal control theory. We aim at fully accounting for measurements of [START_REF] Billat | Differential modeling of anaerobic and aerobic metabolism in the 800-m and 1,500-m run[END_REF][START_REF] Hanon | Pacing strategy and vo2 kinetics during a 1500-m race[END_REF]. We will also provide a proof that Keller's optimal race has exactly the three parts: indeed, Keller makes up a race with three optimal pieces that he computes, without proving the optimality.

The paper is organized as follows: in section 2, we present our new models together with numerical simulations of optimal velocity and anaerobic energy. In section 3, we describe our mathematical analysis and proofs while a conclusion is derived in section 4.

2. Numerical presentation of the models. In this section, we will present our ideas for an improved modeling of races and our numerical findings based on the Bocop toolbox for solving optimal control problems [START_REF] Bonnans | Bocop v1.0.3: A collection of examples[END_REF]. This software combines a user friendly interface, general Runge-Kutta discretization schemes described in [START_REF] Hager | Runge-Kutta methods in optimal control and the transformed adjoint system[END_REF][START_REF] Bonnans | Computation of order conditions for symplectic partitioned Runge-Kutta schemes with application to optimal control[END_REF], and the numerical resolution of the discretized problem using the nonlinear programming problems solver IPOPT [START_REF] Wächter | On the implementation of an interior-point filter line-search algorithm for large-scale nonlinear programming[END_REF].

Our first equation is the equation of motion, as in Keller's paper:

dv dt + v τ = f (t) (2.1)
where t is the time, v(t) is the instantaneous velocity, f (t) is the propulsive force and v/τ is a resistive force per unit mass. The resistive force can be modified to include another power of v. Note that we could take into account a changing altitude, by adding to the right hand side a term of the form -g sin α(d(t)), where α(d(t)) is the slope at distance d(t). We can relate sin α(d) to A(d), the altitude of the center of mass of the runner at distance d, by sin α(d(t)) = A (d(t))

1 + (A (d(t))) 2 .

For most races, one can assume that sin α(d(t)) ∼ A (d(t)) and here, we assume for simplicity that A is constant along the race. Next, we have to establish an equation governing the energy. In fact, human energy can be split into aerobic energy called e ae , which is the energy provided by oxygen consumption, and anaerobic energy e an , which is provided by glycogen and lactate. A very good review on different types of modeling can be found in [START_REF] Morton | The critical power and related whole-body bioenergetic models[END_REF] and a more general reference is [START_REF] Billat | Physiologie et méthodologie de l'entraînement: De la théorie à la pratique[END_REF]. In his paper [START_REF] Joseph | Optimal velocity in a race[END_REF], Keller claims that his energy equation only deals with aerobic energy: he speaks of oxygen supplies. In fact, as we will show below, using a hydraulic analogy, we believe that it well describes the accumulated oxygen deficit: e 0 an -e an (t), where e 0 an is the value at t = 0 of e an (t) the anaerobic energy.

Constraints have to be imposed; the force is controlled by the runner but it cannot exceed a maximal value:

0 ≤ f (t) ≤ f max . (2.
2)

The aim is to minimize the time T , given the distance d = T 0 v(t) dt, with the initial conditions: v(0) = 0, e an (0) = e 0 an under the constraint e an (t) ≥ 0.

(2.

3)

The rest of this section will be devoted to finding a good equation for e an (t). 2.1. First model: how Keller's model describes the accumulated oxygen deficit. Let us use a hydraulic analogy to account for Keller's equations and justify our improvements. This hydraulic analogy is described in [START_REF] Morton | A 3-parameter critical power model[END_REF][START_REF] Morton | The critical power and related whole-body bioenergetic models[END_REF] in order to develop a three parameter critical power model: the equations in [START_REF] Morton | A 3-parameter critical power model[END_REF][START_REF] Morton | The critical power and related whole-body bioenergetic models[END_REF] are on averaged values of the energy and the power, while we will use instantaneous values. We assume that the anaerobic energy has finite capacity modeled by a container of height 1 and surface A p . When it starts depleting by a height h, then what is called in physiology the accumulated oxygen deficit is e 0

an -e an and we relate e 0 an -e an = A p h,

where e 0 an = A p is the initial supply of anaerobic energy, therefore, A p has in fact the dimension of an energy and h is a nondimensionalized height.

We assume that the aerobic energy is of infinite capacity and flows at a maximal rate of σ through R 1 (a connecting tube of fixed diameter illustrated in Figure 2.1). Here, O is the infinite aerobic container, P is the finite capacity anaerobic container, h is the height of depletion of the anaerobic container. An important point is the height at which the aerobic container is connected to the anaerobic one. If we assume in this first model that it is connected at height 1 (at the top of the anaerobic container, and not ϕ for the moment as on Figure 2.1), then it means that the aerobic energy always flows at rate σ. Note that σ is the energetic equivalent per unit time of V O2, the volume of oxygen used by unit of time. This equivalent can be determined thanks to the Respiratory Exchange Ratio and depends on the intensity of effort. Nevertheless, a reasonable average value is that 1l of oxygen produces 20kJ [START_REF] Peronnet | Table of nonprotein respiratory quotient: an update[END_REF]. The available flow at the bottom of the anaerobic container, through T is

W = σ + A p dh dt = σ - de an dt .
Since the energy is used at a rate f v, where v is the velocity and f is the propulsive force, we have that W is equal to the available work capacity hence to f v. This allows us to find the equation governing the evolution of the anaerobic energy

de an dt = σ -f v. (2.5)
We point out that this is exactly the energy equation studied by Keller, except that we have explained that it models the accumulated oxygen deficit, while Keller describes it as the aerobic energy. Some improvements are needed for this model to better account for the physiology:

• change the height where the aerobic container is connected. This implies that the flow of aerobic energy is not constant at the beginning of the race. • take into account that when the energy supply is low, then a physiological control mechanism implies that the flow of energy drops significantly. In other words, V O2 does not keep its maximal value at the end of the race. Before improving the model, we describe our numerical simulations of (2.1)-(2.2)-(2.3)-(2.5) using bocop. We plot, in figure 2.2, the velocity v, the force f , the accumulated oxygen deficit (AOD) e 0

an -e an . We have added σ, though it is constant, just to be consistent with the next figures. We take σ = 49m 2 s -3 , f max = 9.6m s -2 , e 0 an = 1400m 2 s -3 , τ = 0.67s and d = 1500m. The optimal time is 248.21s, and we have 2000 discretization steps, i.e. the time step is close to 0.12s. We display in figure 3.1 a detailed view of what happens at the end of the race for the AOD and the force. We point out that Keller [START_REF] Joseph | Optimal velocity in a race[END_REF] and Woodside [START_REF] Woodside | The optimal strategy for running a race (a mathematical model for world records[END_REF] choose σ = 41.56m 2 s -3 , f max = 12.2m s -2 , e 0 an = 2409m 2 s -3 and τ = 0.892s. They choose these values by minimizing the sum of the squares of the relative errors compared to world records for a series of races. If we take for σ either the value of Keller, or the experimental value of [START_REF] Hanon | Pacing strategy and vo2 kinetics during a 1500-m race[END_REF] corresponding to a V O2 of 66ml mn -1 kg -1 , we do not match as well the velocity curve of [START_REF] Hanon | Pacing strategy and vo2 kinetics during a 1500-m race[END_REF] for our final model. So we have adapted our values of σ, f max , e an and τ to minimize the error on the velocity curve of Figure 2.3 compared to [START_REF] Hanon | Pacing strategy and vo2 kinetics during a 1500-m race[END_REF].

We observe that the race splits into three parts • The race starts with a strong acceleration, the velocity increases quickly and the force is at its maximal value, • for the major part of the race, the force is at an intermediate constant value, the velocity is constant with value close to 6.06 m/s. We will see that this corresponds to what is called, in the optimal control theory, a singular arc. • during the last part corresponding to the last two seconds, the velocity decreases, the force sharply decreases, the energy reaches 0 and then stays at the zero level (AOD is constant equal to e 0 an ), and the force slightly increases again.

We insist on the fact that this is the first simulation not based on the hypothesis that there are three parts of the above type. Also, we can optimize either on the time to run or the distance to run, where all previous simulations had to fix the time to run and optimize on the distance. Even if this is mathematically equivalent in terms of optimization, fixing the distance requires an extra parameter in the simulations. The next models introduce improvements.

2.2. Second model: improving the initial phase to reach V O2 max . The experimental results of [START_REF] Billat | Differential modeling of anaerobic and aerobic metabolism in the 800-m and 1,500-m run[END_REF][START_REF] Hanon | Pacing strategy and vo2 kinetics during a 1500-m race[END_REF] show that the rate of oxygen uptake V O2 is not constant throughout the race but rises steadily from an initial value of about 10ml min -1 kg -1 to its maximum value 66ml min -1 kg -1 over the first 20 to 40 seconds of the race. To model this effect, we now assume that the aerobic container is connected to the anaerobic container at a height ϕ ∈ (0, 1). This implies that there is an initial phase of the race where the flow from the aerobic container to the anaerobic one is no longer σ, but is proportional to the difference of fluid heights in the containers. We assume additionally that there is a residual value at σ r so that,

σ(h) = σ r + (σ -σ r ) h 1-ϕ when h < 1 -ϕ σ when h ≥ 1 -ϕ. (2.6)
This is illustrated in Figure 2.1. We still have the same balance on total work capacity namely

W = f v = σ(h) + A p dh dt . (2.7)
So this and (2.4) lead to the following equations for e an :

de an dt = (σ r + λ(σ -σ r )(e 0 an -e an )) -f v when λ(e 0 an -e an ) < 1, (2.8)
where 1/λ = A p (1 -ϕ). To match the results of [START_REF] Billat | Differential modeling of anaerobic and aerobic metabolism in the 800-m and 1,500-m run[END_REF][START_REF] Hanon | Pacing strategy and vo2 kinetics during a 1500-m race[END_REF], we expect that λ(e 0 an -e an ) reaches 1 in about 20 to 40 seconds, so that we choose ϕ = 0.7, and σ r = 10W kg -1 .

In the second phase, when λ(e 0 an -e an ) has reached 1, we are back to equation (2.5), that is

de an dt = σ -f v when λ(e 0
an -e an ) > 1.

(2.9)

2.3. Third model: drop in V O2 at the end of the race. We want to keep the same initial phases as in the previous model, but take into account that there are limitations when the energy supply is small. The experimental results of [START_REF] Hanon | Pacing strategy and vo2 kinetics during a 1500-m race[END_REF] show that the rate of oxygen uptake V O2 falls slightly at the end of the race, dropping from 66ml min -1 kg -1 to about 60ml min -1 kg -1 over the last 200m to 250m. Morton [START_REF] Morton | The critical power and related whole-body bioenergetic models[END_REF] models this by supposing that the work capacity is proportional to the anaerobic energy e an when e an is small. We prefer to assume that there is a drop in the rate of oxygen uptake V O2 (and therefore the aerobic power σ) when e an is too small. So we add a last phase to the run: when e an /e 0 an < e crit , then σ is replaced by σ e an e 0 an e crit

+ σ f (1 - e an e 0 an e crit ).
We choose for instance e crit = 0.15 and σ f = 0.75σ to fit the drop in [START_REF] Hanon | Pacing strategy and vo2 kinetics during a 1500-m race[END_REF]. We add the final stage:

de an dt = σ e an e 0 an e crit + σ f (1 - e an e 0 an e crit ) -f v when e an e 0
an < e crit .

(2.10)

The coupling of the 3 equations (2.8), (2.9) when ean e 0

an > e crit on the one hand, and (2.10) on the other hand, leads to a better running profile. This model encompasses physiological observations that σ or V O2 is not constant and provides a velocity profile with initial acceleration, then deceleration and acceleration again before the final sprint. This is closer to real races than the constant velocity. It takes care of fatigue with a much better physiological description than [START_REF] Behncke | A mathematical model for the force and energetics in competitive running[END_REF][START_REF] Mathis | The effect of fatigue on running strategies[END_REF]. This model can be summarized as follows an -e an ) ≥ 1

de an dt = σ(e an ) -f v (2.
(σ r + λ(σ -σ r )(e 0
an -e an ) if λ(e 0 an -e an ) < 1

(2.12) together with (2.1)-(2.2)-(2.3). Numerically, in order to avoid piecewise linear functions, we make a regularization in order to have a C ∞ function.

Let us now describe our numerical findings. The results are displayed in Figure 2.3. Since σ is not constant, the singular arc has no longer a constant velocity. Let us be more specific:

• The very first part of the race is still at maximal force with a strong acceleration, • then the force smoothly decreases to its minimal value at the middle of the race, and so does the velocity, • then the velocity and force smoothly increase again • the last part of the race is at maximal force, corresponding to the final sprint, with a final zero energy arc. The final time is 246.961.

2.4. Fourth model: energy recreation when slowing down. On sufficiently long races, it has been observed that slowing down recreates anaerobic energy. This is the motivation of intermittent training. Morton and Billat [START_REF] Morton | The critical power model for intermittent exercise[END_REF] find evidence of energy recreation in intermittent exercises. This can be understood using the hydraulic analogy: if the tap T is closed, then vessel A n refills by virtue of the flow through R 1 . Therefore, on a regular run, the mathematical model has to encompass a new term recreating energy when slowing down. The experimental results for the velocity profiles in [START_REF] Billat | Differential modeling of anaerobic and aerobic metabolism in the 800-m and 1,500-m run[END_REF][START_REF] Hanon | Pacing strategy and vo2 kinetics during a 1500-m race[END_REF] show that the velocity is not constant during the bulk of the race, as predicted by the previous models, but oscillates about a constant value. The recreation of anaerobic energy when slowing down can produce such oscillations, but only if the total Hamiltonian is non-convex. This can be achieved by replacing (2.11) by

de an dt = σ(e an ) + η(a) -f v (2.13)
where a = dv dt is the acceleration, if the energy recreation term η(a) is a power of a strictly larger than 1 when a is negative. We choose η(a) = ca 2 -where a -is the We will see that this provides a simulation of velocity which is consistent with real runs.

Because of the term η(a), the hamiltonian gets non convex, so that by Pontryagin's maximum principle (see our analysis in section 3.3.2) the optimal solution oscillates between the maximal and minimal value of the force (i.e. f max and 0). This is in fact to be understood in a relaxed sense, as a probability of taking the maximal and minimal values of the force. However, a runner cannot change his propulsive force instantaneously for several reasons: the information to vary its propulsive force takes some time to reach the brain, and the dynamics of the bone and muscles take some time. We choose to take this into account by bounding the derivative of f :

| df dt | ≤ C. (2.15)
The simulations of (2. We see that the force, having a bounded derivative, does not oscillates between its maximal value and 0, but between its maximal value, and some lower value, the derivative of the force reaching its bounds. Consequently, the velocity oscillates and so does the energy which gets recreated. These oscillation reproduce qualitatively the measurements of [START_REF] Billat | Differential modeling of anaerobic and aerobic metabolism in the 800-m and 1,500-m run[END_REF][START_REF] Hanon | Pacing strategy and vo2 kinetics during a 1500-m race[END_REF]. 2.5. A periodic pattern. The previous experiments show a behavior of the optimal control which looks, in the time interval [START_REF] Osmolovskii | Applications to regular and bang-bang control[END_REF]240], i.e., except for the initial and final part of the trajectory, close to a periodic one. We have approximately e an (20) = 950 and e an (240) = 100, and so the average decrease per unit time is e d = 3.86. We observe that over this time interval the speed varies between 5.8 and 6.3 m/s, and the force varies between 8.2 and 9.6.

This leads us to consider the problem of maximizing the average speed over a period T (the period itself being an optimization parameter): the periodicity conditions apply to the speed and force, and the energy is such that e(0) = e(T ) + T e d . In other words, we wish to solve the following optimal control problem:

       Min - 1 T T 0 v(t)dt; v(0) = v(T ); e(0) = e(T ) + T e d ; f (0) = f (T ). v(t) = f (t) -v(t) τ ; ė(t) = σ + η(a(t)) -f (t)v(t); 0 ≤ f (t) ≤ f M ; | ḟ (t)| ≤ 1, for a.e. t ∈ (0, T ). (2.16)
We can fix the initial energy to 0. Now we can compare figure 2.5, where we made a zoom on the solution computed before over the time interval [113,119], with the solution of the periodic problem, displayed in figure 2.6, with period 3.5. We observe a good agreement between those two figures, which indicates that computing over a period may give a good approximation of the optimal trajectory.

3. Mathematical analysis. We have to study optimal control problems with a scalar state constraint and a scalar control, which, in some of the models enters linearly into the state equation. We mention among others the related theoretical studies by Bonnans and Hermant [START_REF] Bonnans | Second-order analysis for optimal control problems with pure state constraints and mixed control-state constraints[END_REF] about state constrained problems, Maurer [START_REF] Maurer | On optimal control problems with bounded state variables and control appearing linearly[END_REF], who considers problems with bounded state variables and linear control, Felgenhauer [START_REF] Felgenhauer | Structural stability investigation of bang-singular-bang optimal controls[END_REF] about the stability of singular arcs, and the two recent books by Osmolovskii and Maurer [START_REF] Osmolovskii | Applications to regular and bang-bang control[END_REF] and Schättler and Ledzewicz [START_REF] Schättler | Geometric optimal control[END_REF].

Statement of the model. we consider the following state equation

ḋ = v; v = f -φ(v); ė = σ(e) -f v, (3.1)
where the drag function φ satisfies

φ is a C 2 function; φ(0) = 0, φ positive, vφ (v) nondecreasing. ( 3
For v ≥ 1, we have that φ (v) ≥ φ (1)/v, so that φ(v) ≥ φ(1) + φ (1) log v. Conse- quently, φ(v) ↑ +∞ when v ↑ +∞. (3.3)
We assume for the moment (in section 3.3 we will discuss a more general recreation model) that the recreation function σ(e) satisfies σ(e) is C 2 and nonnegative.

We will often mention Keller's model which corresponds to the case where φ(v) = v/τ and σ(e) is a positive constant.

(3.5)

As before, the initial condition is

d(0) = 0; v(0) = 0; e(0) = e 0 > 0, (3.6)
and the constraints are

0 ≤ f (t) ≤ f M ; e(t) ≥ 0; t ∈ (0, T ); d(T ) = D. (3.7)
The optimal control problem is to minimize the final time: Note that we could as well take the final constraint as -d(T ) ≤ -D. Writing an inequality in this way yields the sign of the Lagrange multiplier. We establish in Appendix A that the optimal solutions (for a given distance D) with corresponding time T are also solutions of the problem of maximizing the distance over a time interval T . So, in the following, we will rather use this second formulation. Some useful tools of the optimal control theory are introduced and discussed in Appendix D.

Main results

. If 0 ≤ a < b ≤ T is such that e(t) = 0 for t ∈ [a, b], but e(•)
does not vanish over an interval in which [a, b] is strictly included, then we say that (a, b) is an arc with zero energy. Similarly, if f (t) = f M a.e. over (a, b) but not over an open interval strictly containing (a, b), we say that (a, b) is an arc with maximal force. We define in a similar way an arc with zero force. We say that a (resp. b) is the entry (resp. exit) time of the arc. A singular arc is one over which the bound constraints are not active (the associated multiplier is a.e. equal to 0).

If the distance is small enough, then the optimal strategy consists in setting the force to its maximal value. Let D M > 0 be the supremum (assumed to be finite) of the distance for which this property holds. Using standard arguments based on minimizing sequences and weak topology (based on the fact that the control enters linearly into the state equation), the following can be easily proved:

Lemma 3.1. The above problem has at least one optimal solution. Theorem 3.2. Assume that D > D M . Then: (i) An optimal trajectory starts with a maximal force arc, and is such that e(T ) = 0. (ii) If σ is a positive constant, an optimal trajectory has the following structure: a maximal force arc, followed or not by a singular arc, and a zero energy arc.

This will be a consequence of Theorem 3.16 and Remark 3.3. Since, by proposition A.1, the optimal solutions are solutions of problems of maximizing the achieved distance in a given time, we consider in the sequel problems with a given final time T such that the strategy of constant maximal force is not feasible.

Recreation when decreasing the speed.

3.3.1. Framework. We next consider a variant of the previous model where the dynamics of the energy is a sum of functions of the energy and the acceleration a(t) = f (t) -φ(v(t)). We can write the dynamics are as follow:

v = f -φ(v), ė = σ(e) + η(a) -f v, (3.9)
with initial conditions v(0) = 0; e(0) = e 0 > 0, (3.10)

The constraints are

0 ≤ f (t) ≤ f M ; e(t) ≥ 0; t ∈ (0, T ). ( 3.11)
The recreation optimal control problem is as follows: This is a C 1 and nondecreasing function of v, with value 0 at 0. The following function appears in the analysis of singular arcs :

Min - T 0 v (
Q(v) := vφ (v) + φ(v) -R (v)(1 -φ(v)/f M ) + φ (v)R(v)/f M . = (vφ(v)) + (φ(v)R(v)) /f M -R (v). (3.15)
Let v M denote the supremum of speeds that can be reached by taking f (t) = f M for all t, defined in (B.1). In some of the results, we need to assume that

Q is an increasing function of v over [0, v M ]. (3.16)
In Keller's model, we have that R(v) = 0 and Q(v) = 2v/τ , so that the previous hypothesis holds. Remark 3.3. More generally, if η has the following structure η = cη, with c ≥ 0, η(s) > 0 for all s > 0, (3.17)

since by (3.2) we have that vφ (v) + φ(v) is increasing, we have that (3.16) holds if c > 0 is small enough.

A solution of (3.12) called an optimal control solution does not necessarily exist:

Theorem 3.4. Assume that (3.16) holds, that (3.17) holds with c > 0, and that σ(•) is a positive constant. Then the optimal control problem (3.12) has no solution.

This will be proved at the end of the section, as a consequence of the analysis of the relaxed problem that we now perform. The theorem motivates the introduction of a relaxed problem.

Relaxed problem.

In the relaxed formulation for which we refer to [START_REF] Ekeland | Convex analysis and variational problems[END_REF], we replace the control f (t) with a probability distribution Ξ(t, f ) with values in [0, f M ]. Denoting by IE Ξ(t) the expectation associated with this probability measure, and by Ξ(t) the expectation of f at time t, the state equation becomes

v(t) = Ξ(t) -φ(v(t)), ė(t) = σ(e(t)) -Ξ(t)v(t) + IE Ξ(t) η(f -φ(v(t))).
(3.18)

The relaxed optimal control problem is The Hamiltonian function is the same as for the non relaxed version, and its expression is

H[p](f, v, e) := -v + p v (f -φ(v)) + p e (σ(e) -f v + η(f -φ(v)). (3.20)
The costate equation is therefore, omitting time arguments: The variables p v and p e are called costate variables. By standard arguments based on minimizing sequences, we obtain that Lemma 3.5. For c ≥ 0 small enough, the relaxed optimal control problem (3.19) has at least one solution.

-ṗv = -1 -p v φ (v) -p e Ξ -p e φ (v)IE Ξ(t) η (f -φ(v)
By lemma D.2, any feasible point of the relaxed optimal control problem (3.19) is qualified.

Lemma 3.6. For c ≥ 0 small enough, the relaxed optimal control problem (3.19) has the same value as the (non relaxed) problem (3.12), and therefore, any solution of (3.12) is solution of (3.19).

Proof. Let Ξ be a feasible point of the relaxed problem. For c ≥ 0 small enough, by lemma D.2, this point is qualified, so that for any ε > 0, there exists a feasible point Ξ of the relaxed problem such that Ξ -Ξ ≤ ε in the norm of L ∞ (0, T, M ([0, T ])), with associated state (v , e ) so that the state constraint e(t) ≥ 0 holds with strict inequality for all t ∈ [0, T ]. It follows [START_REF] Ekeland | Convex analysis and variational problems[END_REF] that Ξ can be approximated by a classical control f whose state (v f , e f ) satisfies in the uniform norm v f -v + e f -e ≤ ε f , for arbitrary ε f > 0. When ε f ↓ 0 we have that f is feasible for the (unrelaxed) optimal control problem, and the associated cost converges to the one associated with Ξ . The conclusion follows.

We have then that Pontryagin's principle holds in qualified form, i.e., with each optimal trajectory (f, v, e) is associated at least one multiplier (p, µ) such that the relaxed control minimizes the Hamiltonian, in the sense that Ξ(t, •) has support in argmin{H[p(t)](f, v(t), e(t)); f ∈ [0, f M ]}, for a.a. t. (3.25) Lemma 3.7. Let t 1 ∈ (0, T ] be such that e(t 1 ) = 0. Then: (i) we have that

v(t 1 ) ≥ σ(0)/f M . (3.26)
(ii) An arc of maximal force cannot start at time t 1 , or include time t 1 .

Proof. (i) If (3.26) does not hold, then σ(0) -v(t 1 )f M > 0, and therefore ė(t) ≥ σ(t) -v(t)f M > 0 for t close enough to t 1 , contradicting the fact that e(t) reaches its minimum at time t 1 . (ii) Let t 1 contradict point (ii). By point (i), ė(t 1 +) = σ(0) -f M v(t) ≤ 0, but since the energy must be nonnegative, this is an equality, and so ė(t + ) = 0. Therefore,

ë(t + 1 ) = σ (0) ė(t + 1 ) -f M v(t + 1 ) = -f M v(t + 1 ) < 0, (3.27)
implying that the energy cannot be positive after time t 1 . This gives the desired contradiction. Lemma 3.8. We have that (i) T ∈ supp(dµ), so that e(T ) = 0, and

(ii) p e (t) < 0, t ∈ [0, T ). (iii) p v (t) < 0, t ∈ [0, T ). (3.28)
Proof. (i) If T ∈ supp(dµ), let t e := max supp(dµ). Then t e ∈ (0, T ) and e(t e ) = 0: indeed, if t e = T , then min t e(t) > 0 and, since by assumption the strategy of maximal force is not feasible, one can easily construct by slightly increasing the force a feasible trajectory that reaches a smaller cost. We analyze what happens over (t e , T ). Since p e has derivative ṗe = -p e σ (e(t)) and is continuous with zero value at time T , it vanishes, and so, ṗv = 1 + p v φ (v). Since p v (T ) = 0, this implies that p v has negative values. As p e vanishes and p v is negative, the Hamiltonian is equal to -v + p v (f -φ(v)) and has a unique minimum at f M . It follows that (t e , T ) is included in a maximal force arc, which, since e(t e ) = 0, is in contradiction with lemma 3.7(ii). Point (i) follows. (ii) If p e (t a ) ≥ 0 for some t a ∈ (0, T ), then, by the costate equation, p e should vanish on (t a , T ] and so (t a , T ) would not belong to the support of dµ, in contradiction with (i). This proves (ii). (iii) Let on the contrary t c ∈ [0, T ) be such that p v (t c ) ≥ 0. Then the Hamiltonian is a sum of nondecreasing functions of f and has a unique minimum point at 0, and since p e (t c ) < 0, the same holds for t close enough to t c . Therefore t c belongs to a zero force arc, along which

ṗv = 1 + p v φ (v) + p e η (-φ(v))φ (v) (3.29)
remains positive (remember that p e (t) < 0 and that η is nonincreasing), and so this arc cannot end before time T . But then we cannot meet the final condition p v (T ) = 0. The conclusion follows.

Since p e (t) < 0 over (0, T ), we deduce by (3.25) that when η(a) > 0 for all a < 0, H is a concave function of f with minima over [0, f M ] at {0, f M } for a.a. t, and so, by (3.25):

Corollary 3.9. If η(a) > 0 for all a < 0, an optimal control for the relaxed control problem (3.19) has support over {0, f M } for a.a. t.

Reformulation of the relaxed problem.

In view of the previous corollary, when η(a) > 0 for all a < 0, we may as well restrict the study of the relaxed optimal control problem (3.19) to the case when the relaxed control has values in {0, f M }. Such a relaxed control can be parameterized by its expectation f (t) at any time t: the probability to take the value 0 is 1 -f (t)/f M . Remembering the definition of the recovery obtained with a zero force in (3.14), the dynamics can now be written in the form, since η(f M -φ(v)) = 0 along the trajectory: (3.31) Remark 3.10. When η identically vanishes, the above problem still makes sense and coincides with the formulation of the original model of section 3.1. So we will be able to apply the results of this section to Keller's problem.

   v(t) = f (t) -φ(v(t)), ė(t) = σ(e(t)) -f (t)v(t) + 1 - f (t) f M R(v(t)), (3.30
The Hamiltonian is

H R = -v + p v (f -φ(v)) + p e (σ(e) -f v) + p e 1 - f f M R(v). (3.32)
The costate equations are

   -ṗv = -1 -p v φ (v) -p e f + p e 1 - f f M R (v),
-dp e = p e σ (e)dt -dµ.

(3.33)

Of course we recover the expressions obtained in section 3.3.2 in the particular case of relaxed controls with support in {0, f M }, and therefore all lemmas of this section are still valid. By construction, the Hamiltonian is an affine function of the control f , and the switching function is

Ψ R = H R f = p v -p e (v + R(v)/f M ). (3.34)
Lemma 3.11. Any optimal trajectory starts with a maximal force arc. Proof. Since v(0) = 0, we have that Ψ R (0) = p v (0) is negative by lemma 3.8(iii). The conclusion follows from the fact that the optimal control minimizes the Hamiltonian, see e.g. [START_REF] Bonnans | Second-order analysis for optimal control problems with pure state constraints and mixed control-state constraints[END_REF].

The following hypothesis implies that the energy is nonzero along a zero force arc.

Either σ(0) > 0, or η is nonzero over R -.

(3.35) Lemma 3.12. Let (3.35) hold. Then µ has no jump over [0, T ]. Proof. a) Let t ∈ [0, T ) be such that [µ(t)] > 0. Necessarily e(t) = 0, and so t ∈ (0, T ), and Ψ R (t -) ≤ 0 (since otherwise t belongs to a zero force arc and then by (3.35) we cannot have e(t) = 0). We have that [p e (t)] = [µ(t)] > 0, and so [Ψ R (t)] = -(v(t) + R(v)/f M )[p e (t)] < 0, implying Ψ R (t + ) < 0. Therefore, for some ε > 0, (t, t + ε) is included in a maximal force arc, in contradiction with lemma 3.7(ii). b) If [µ(T )] = 0, since lim p v (t) = 0 and lim p e (t) = -[µ(T )] when t → T , and v(T ) > 0, we get lim t↑T Ψ R (t) = [µ(T )](v(T ) + R(T )/f M ) > 0, meaning that the trajectory ends with a zero force arc, but then by (3.35), the energy cannot vanish at the final time, contradicting lemma 3.8(i).

By lemma 3.12, the switching function is continuous. When the state constraint is not active, its derivative satisfies

ΨR = 1 + p v φ (v) + p e f -p e 1 - f f M R (v) +p e σ (e) (v + R(v)/f M ) -p e (1 + R (v)/f M ) (f -φ(v)).
(3.36)

The coefficient of f vanishes as expected and we find that

ΨR = 1 + p v φ (v) + p e φ(v) + p e σ (e)(v + R(v)/f M ) -p e R (v)(1 -φ(v)/f M ). (3.37)
Subtracting Ψ R φ (v) in order to cancel the coefficient of p v , we obtain that

ΨR -Ψ R φ (v) = 1 + p e σ (e) (v + R(v)/f M ) + p e Q(v), (3.38)
where Q was defined in (3.15). Lemma 3.13. For given 0 ≤ t 1 < t 2 ≤ T , let (t 1 , t 2 ) be included in a singular arc over which σ (e(t)) is equal to 0. Then, over (t 1 , t 2 ), Q(v) is constant and, if the function Q is not constant on any interval, v is constant. 3.4. Bounding variations of the force. It seems desirable to avoid discontinuities of the force that occur with the previous model, and for that we introduce bounds on ḟ . The force becomes then a state and the new control ḟ is denoted by g. So the state equation is (note that we have taken here

σ = 0) v = f -φ(v); ė = σ(e) + η(a) -f v; ḟ = g, (3.40)
with constraints

0 ≤ f ≤ f M ; e ≥ 0; g m ≤ g ≤ g M . (3.41)
We minimize as before -T 0 v(t)dt. The Hamiltonian is (3.43)

H = -v + p v (f -φ(v)) + p e (η(f -φ(v)) -f v) + p f g. ( 3
The state constraint e ≥ 0 is of second order, and we may expect a jump of the measure µ at time T . The final condition for the costate are therefore

p v (T ) = 0; p e (T ) = 0; p f (T ) = 0. (3.44)
We may expect and will assume that the above two inequalities are strict. By the analysis of the previous section we may expect that the optimal trajectory is such that g is bang-bang (i.e., always on its bounds), except if a state constraint is active (the state constraints now include bound constraints on the force), as is confirmed in our numerical experiments. The bounds on the derivative of the force modify the last part of the race compared to Keller's (see figure 3.1): the drop in force cannot be so sharp.

Conclusion.

We have established a system of ordinary differential equations governing the evolution of the velocity v, the anaerobic energy e an , and the propulsive force f . This is based on the equation of motion (relating the acceleration a = dv/dt to the propulsive force and the resistive force) and a balance of energy. Several constraints have to be taken into account: the propulsive force is positive and less than a maximal value, its derivative has to be bounded, the anaerobic energy is positive. Keller [START_REF] Joseph | A theory of competitive running[END_REF][START_REF] Joseph | Optimal velocity in a race[END_REF] used in his model the evolution of the aerobic energy, which is not satisfactory. Here, using a hydraulic analogy initiated by Morton [START_REF] Morton | A 3-parameter critical power model[END_REF][START_REF] Morton | The critical power and related whole-body bioenergetic models[END_REF], we manage to write an equation for the instantaneous accumulated oxygen deficit instead. In our model, in difference with respect to Keller's, we introduce variations in σ, modeling the oxygen uptake, V O2: indeed, one of the roles of the anaerobic energy is to compensate the deficit in oxygen uptake, V O2, which has not reached its maximal value at the beginning of the race. Conversely, when the anaerobic energy gets too low, the oxygen uptake V O2 cannot be maintained to its maximal value. We make two further extensions: we introduce a physiological observation that energy is recreated when the acceleration is negative, that is when decreasing the speed, and the fact that the derivative of the propulsive force has to be bounded.

Our model could be used, in the simulations, given the velocity profile of a runner, to compute the evolution of his anaerobic energy. This is an important challenge for sportsmen to determine instantaneous anaerobic energy consumption.

In this paper, we use our system for the optimization of strategy in a race: given a distance, we want to find the optimal velocity leading to the shortest run. Our main results are illustrated in Figures 2.3, 2.4. Without the energy recreation term, we find that the race starts with a strong acceleration to achieve a peak velocity, then the race is run at an average velocity, first with a decreasing pace, and then an increasing pace, before the final sprint. Our numerical simulations on the final model where energy recreation is involved provide oscillations of the velocity that qualitatively reproduce the physiological measurements of [START_REF] Billat | Differential modeling of anaerobic and aerobic metabolism in the 800-m and 1,500-m run[END_REF][START_REF] Hanon | Pacing strategy and vo2 kinetics during a 1500-m race[END_REF].

Using optimal control theory, we manage to get rigorous proofs of most of our observations. We prove in particular that in the case of Keller, the race is made up of exactly three parts: run at maximal propulsive force, run at constant speed (corresponding to a singular arc), run at zero energy. It cannot be made of any other arcs. For this purpose, we relate the problem to a relaxed formulation, where the propulsive force represents a probability distribution rather than a function of time. We also find that the concavity of the Hamiltonian results in speed oscillations and we show how, by reducing the problem on optimizing over a period, we recover the latter.

Appendix A. Abstract distance and time functions. In this section, we establish, in a general setting, the relation between the distance and time functions, defined in (3.1), (3.6), (3.7). Let the following be the control and state spaces : Let f ∈ U M be a feasible control, i.e., which is such that f ∈ U M and e[f ] ∈ K. The constraints are said to be qualified at f (see [START_REF] Robinson | First order conditions for general nonlinear optimization[END_REF] or [7, section 2. In other word, the variation δf of the control is compatible with the control constraints, and the linearized state δe allows to reach the interior of the set of feasible states. Remember that e(0) > 0. Lemma D.1. If c is small enough, the optimal control problem (3.12) is qualified. Proof. a) We first obtain the result when c = 0. If e(t) is always positive the qualification holds with δf = 0. Otherwise, let t a be the smaller time at which the energy vanishes with δf = -f . Obviously f + δf ∈ U M , and since δf is a.e. nonpositive, so is δv. Next, since c = 0, we have that δe(t) = σ (e(t))δe(t) -δf (t)v(t) -f (t)δv(t), t ∈ (0, T ); δe(0) = 0, (D. Taking ε ∈ (0, t a -t f ), it follows that δe(t) > 0 over [t a , T ], and so e(t) + δe(t) is positive over [0, T ], and hence uniformly positive as was to be shown. b) We now show that for c > 0 small enough the qualification can be obtained, again by taking δf = -f . Given a sequence c k of positive number converging to 0 and f k ∈ U M , we may extract a subsequence such that f k converges to f in L ∞ weak * , and since c k ↓ 0, the associated states (v k , e k ) uniformly converge to the associated state (v, e). Since δf k converges to δf in L ∞ weak * , we deduce that (δv k , δe k ) uniformly converge to (δv, δe). By point (a), e k + δe k is (uniformly) positive over [0, T ], as was to be shown. We now consider the relaxed formulation of section 3.3.2. Lemma D.2. If c is small enough, the optimal control problem (3.19) is qualified.

Proof. The proof is essentially the same, up to technical details (the main point is that for the variation of the control we still take the opposite of the control), and is left to the reader.

from 50 m to 275 km). Math. Comput. Modelling, 15(10): 1-12, 1991.

Fig. 2 . 1 .

 21 Fig. 2.1. Scheme of the container modeling.

Fig. 2 . 2 .

 22 Fig. 2.2. Simulations of (2.1)-(2.2)-(2.3)-(2.5). Plot of the velocity v, the accumulated oxygen deficit (AOD) e 0 an -ean, the propulsive force f , and σ vs time.

≥

  e crit and λ(e 0

Fig. 2 . 3 .

 23 Fig. 2.3. Variable σ: Plot of the velocity v, the accumulated oxygen deficit (AOD) e 0an -ean, the propulsive force f , and σ vs time.

Fig. 2 . 4 .

 24 Fig. 2.4. Recreation when slowing down: speed, energy, force, σ.

  Min T ; s.t. (3.1) and (3.6)-(3.7).(3.8)

  t)dt; s.t. (3.9)-(3.11) hold. (3.12) We will assume that η is a convex and C 1 function, that vanishes over R + . (3.13) This implies that η in nonincreasing. A typical example is η(a) = c|a -| β , with c ≥ 0, β ≥ 1, and a -:= min(a, 0). Let us denote the maximal recovery function (obtained with a zero force) by R(v) := η(-φ(v)). (3.14)

  dt; s.t. (3.18) and (3.10)-(3.11) hold. (3.19)

  ), -dp e = p e σ (e)dt -dµ, (3.21) with final conditions p v (T ) = p e (T ) = 0. (3.22) Here dµ, identified to the bounded variation function µ on [0, T ], is a Borel measure (that can be interpreted as a Lagrange multiplier) associated to the state constraint -e ≤ 0. It satisfies the properties of nonnegativity and complementarity, i.e., dµ ≥ 0; T 0 e(t)dµ(t) = 0. (3.23) Note that equivalent relations are, denoting by supp(•) the support of a measure: dµ ≥ 0; supp(dµ) ⊂ {t ∈ [0, T ]; e(t) = 0}. (3.24)

  dt; s.t. (3.30) and (3.10) hold, and 0 ≤ f ≤ f M a.e., e ≥ 0 on [0, T ].

Fig. 3 . 1 .

 31 Fig. 3.1. Race problem with Keller's model: zoom on end of race: AOD and force.

  .42) The costate equation -ṗ = H y are now    -ṗv = -1 -p v φ (v) -p e (η (a)φ (v) + f ), -dp e = σ (e)p e dt -dµ, -ṗf = p v + p e (η (a) -v).

U

  T := L ∞ (0, T ); Y T := L ∞ (0, T ; R n ). (A.1) Given δ : R n → R, F : R × R n → R n , and K T ⊂ U T × Y T ,we consider the "abstract" problems of minimal time (P D ) Min T ; ẏ(t) = F (f (t), y(t)), for a.a. t ∈ (0, T ), y(0) = y 0 ; (f, y) ∈ K T ; δ(y(T )) = D, and of maximal distance ( PT ) Max δ(y(T )); ẏ(t) = F (f (t), y(t)), t ∈ (0, T ), y(0) = y 0 ; (f, y) ∈ K T .

3 . 4 ]

 34 ) if there exists δf ∈ U such thatf + δf ∈ U M ; e[f ] + δe[δf ] ∈ int(K).(D.4)

  Zoom of the case with recreation when slowing down (2000 time steps)

				force f(t)							velocity v(t)
	9.8
	9.6						6.20
	9.4						6.15
	9.2						6.10
	9.0						6.05
	8.8
							6.00
	8.6
							5.95
	8.4
							5.90
	8.2
	8.0						5.85
	7.8						5.80
	113	114	115	116	117	118	119	113	114	115	116	117	118	119
		Fig. 2.5.

  .2)

					force f(t)									velocity v(t)
	10.0								6.4
	9.8
									6.3
	9.6
	9.4								6.2
	9.2								6.1
	9.0
									6.0
	8.8
	8.6								5.9
	8.4
									5.8
	8.2
	8.0								5.7
	2	3	4	5	6	7	8	9	10	2	3	4	5	6	7	8	9	10
							Fig. 2.6. Optimization over a period

  [START_REF] Bonnans | Computation of order conditions for symplectic partitioned Runge-Kutta schemes with application to optimal control[END_REF] which implies δe ≥ 0. Let t b be the essential supremum of times for with f (t) = 0 for a.a. t ∈ (0, t b ). Clearly, t b < t a , and for any ε > 0, there exists α > 0 such that v(t) ≥ α and -δv(t) ≥ α, for all t ∈ [t b + ε, T ].

(D.6) Let C := max{|σ (e(t))|; t ∈ [0, T ]}. Then for t in (t f + ε, T ]: δe(t) ≥ σ (e(t))δe(t) -f (t)δv(t) ≥ -Cδe(t) + 2αf (t).

(D.7)

Acknowledgments. The first author would like to thank Christophe Clanet for pointing out Keller's contributions, Véronique Billat, Christine Hanon for discussions on the physiological models and the bibliography, Daphne Giorgi, Pierre Martinon for their help on bocop and Nadège Arnaud who made figure 2.1. The second author acknowledges support from the EU 7th Framework Programme (FP7-PEOPLE-2010-ITN), under GA number 264735-SADCO and thanks Bastien Berret for his comments. Finally, both authors thank the two referees for their useful remarks.

Proof. Along a singular arc, dµ vanishes, so that p e is constant. By (3.38), so is also Q(v) = -1/p e . The conclusion follows.

Remark 3.14. If η is analytic over R -and φ is analytic over R + , then Q is analytic over R + , and therefore is either constant over R, or not constant on any interval of R + .

We say that t ∈ (0, T ) is a critical time if Ψ R (t) = 0, and we say that t± is energy free if e(t) > 0 for t in ( t, t ± ε), for ε > 0 small enough (but it may happen that e( t± ) = 0. As for (3.38), we have that, for such times the existence of left or right derivatives: ΨR t± = 1 + p e (t ± )σ (e(t)) (v(t) + R(v(t))/f M ) + p e (t ± )Q(v(t)), (3.39) where p e (t ± ) is the right or left limit of p e at time t, and Q(•) was defined in (3.15). Lemma 3.15. Let σ be a positive constant, and (3.16) hold. Then, along an optimal trajectory: (i) There is no zero force arc, and hence, Ψ R (t) ≤ 0, for all t ∈ [0, T ]. (ii) The only maximal force arc is the one starting at time 0.

Proof. (i) Let (t a , t b ) be a zero force arc, over which necessarily Ψ R is nonnegative. By lemma 3.11, t a > 0, and so Ψ R (t a ) = 0 and ΨR (t a+ ) ≥ 0. Since σ is constant and positive, over the arc, e(t) > 0, p e is constant, the speed is decreasing, and so by (3.39), we have that ΨR (t b-) > ΨR (t a+ ) ≥ 0 meaning that the zero force arc cannot end before time T , contradicting the final condition e(T ) = 0. (ii) On a maximal force arc (t a , t b ) with t a > 0, since the speed increases, (3.39) implies ΨR (t b-) < ΨR (t a+ ) ≤ 0, and since Ψ R t ≤ 0 along the maximal force arc, it follows that Ψ R (t b ) < 0, meaning that the maximal force arc ends at time T . But then [p e (T )] = [µ(T )] > 0, in contradiction with lemma 3.12.

Theorem 3.16. Let σ be a positive constant, and (3.16) hold. Then an optimal trajectory has the following structure: maximal force arc, followed or not by a singular arc, and a zero energy arc.

Proof. The existence of a maximal force arc starting at time 0 is established in lemma 3.11. Let t a ∈ (0, T ) be its exit point (t a = T is not possible since T > T M ), and let t b ∈ (0, T ) be the first time at which the energy vanishes (that t b < T follows from lemmas 3.8(i) and 3.12). If t a < t b , over (t a , t b ), by lemma 3.15, Ψ R is equal to zero and hence, (t a , t b ) is a singular arc. Finally let us show that the energy is zero on (t b , T ). Otherwise there would exist t c , t d with t b ≤ t c < t d ≤ T such that e(t c ) = e(t d ) = 0, and e(t) > 0, for all t ∈ (t c , t d ). By lemma 3.15, (t c , t d ) is a singular arc, over which ė = σ -f v is constant, which gives a contradiction since the energy varies along this arc. The result follows.

Proof. [Proof of theorem 3.4] By lemma 3.6, any solution of the classical problem is solution of the relaxed one. By theorem 3.16, the trajectory must finish with an arc of zero energy, over which 0

). The r.h.s. is a strictly decreasing function of f (t). We deduce that f (t) is a continuous function of v(t), and hence, of time over this arc. On the other hand, since p e < 0 a.e., the Hamiltonian is a concave function of f which is not affine on [0, f M ], and so attains its minima at either 0 or f M . Therefore f (t) is constant and equal to either 0 or f M over the zero energy arc. For f (t) = 0 we have that ė(t) is positive. That f (t) = f M is not possible since we know that the trajectory has only one maximal force arc. We have obtained the desired contradiction.

Remark 3.17. We plot in figure 3.1 a zoom on the end of the race. While we have proved that for the continuous problem, there is a switching time from the singular arc, with constant speed, to the zero energy arc, we observe in the discretized problem a progressive transition between 246 and 248 seconds.

In our examples, y = (h, v, e) and δ(y) = y 1 is the distance.

Proposition A.1. Denote by T (D) and D(T ) the optimal values of the above problems. Assume that these functions are finitely valued, increasing and continuous over R + with value 0 at 0. Then (a) T (D) is the inverse function of D(T ), and (b) any optimal solution of ( PT ) (resp. (P D )) is a solution of (P D(T ) ) (resp. ( PT (D) )).

Proof. (a1) Given ε > 0, let ( f , ŷ) be as above and such that T ≤ T (D) + ε. Then

The first inequality is due to the fact that the trajectory ( f , ŷ) is feasible for ( P T ), and the second one holds since D is nondecreasing. Passing to the limit when ε ↓ 0 and using the continuity of D, we deduce that

(a2) Given ε > 0, let ( f , ȳ) be as above and such that D(T ) -ε ≤ D. Then

The first inequality holds since T is nondecreasing, and the second one is due to the fact that the trajectory ( f , ȳ) is feasible for (P D ). Passing to the limit when ε ↓ 0 and using the continuity of T , we deduce that Appendix B. Strategy of maximal force. The strategy of maximal force is the one for which the force always has its maximal value. Then speed is an increasing function of time, with positive derivative, and asymptotic value

We first discuss the existence of a critical distance D M at which the energy vanishes, under the following hypothesis, implying that, when reaching the maximal speed with a maximal force, the energy decreases :

With the strategy of maximal force, if (B.3) holds, the energy cannot remain nonnegative for all time t ≥ 0. (ii) If the energy vanishes at time t M , then the maximal force strategy does not respect the contraint of nonnegative energy over [0, t] for any t > t M .

Proof. (i) By (B.3), there exists ε M > 0 such that sup e≥0 σ(e)) + ε < f M v M . For large enough time, ė(t) ≤ -ε M so that e(t) → -∞; point (i) follows. (ii) If the conclusion does not hold, then e(t) attains a local minimum over (0, τ ) at time t M , and so we have e(t M ) = 0 and ė(t M ) = σ(0) -f M v(t M ) = 0. Since the speed has a positive derivative, it follows that

and therefore in any case e(t) < 0 for t > t M , close to t M , which gives the desired contradiction.

Appendix C. Link with the Bellman function.

Consider the following extension of the relaxed problem (3.19), where the problem is still to compute the maximal distance one can run in time T , but with general intial conditions v 0 = v 0 , e 0 = e 0 > 0. The value function V (T, v 0 , e 0 ) (equal to the opposite of the optimal distance) is called in this setting the Bellman value. The analysis of section 3 is based on the (negative) sign of the costate variables. It is known that, under some hypotheses, the costate is equal to the the gradient of the Bellman value. So, it is of interest to check in a direct way (without using the optimal control theory) the following result.

Lemma C.1. The function V (T, v 0 , e 0 ) : R 3 + → R is a decreasing function of each of its three arguments.

Proof. That D is a increasing function of T is easy to prove. When changing the initial energy from e 0 to ê0 > e 0 , given an optimal control and state ( f , v, ē), we have that f is feasible for the new problem. Indeed, the state has the same speed v and a new energy ê that must satisfy ê(t) ≥ ē(t) for all t ≥ 0, since if equality holds at some time t e > 0, then ê(t) = ē(t) for all t > t e . It follows that D is an increasing function of e 0 . Finally, let us change the initial speed to v0 > v 0 . If, for a zero force strategy, the corresponding speed v is always greater than v over (0, T ), the conclusion holds. Otherwise, let t a ∈ (0, T ) be such that when applying the zero force over (0, t a ), we have that v(t a ) = v(t a ). Define the strategy f to have value 0 over (0, t a ), and to be equal to f otherwise.

Clearly, the distance at time t a is greater than the corresponding one for the original strategy, and the energy denoted by ê satisfies ê(t a ) > ē(t a ) (equality is not possible since it would mean that f (t) = 0 = f (t) for all t ∈ (0, t a ), but then v(t a ) > v(t a )). Since we know that D is an increasing function of energy, the conclusion follows.

Appendix D. Qualification. We consider the model with energy recreation of section 3.3. We assume that the functions σ, Φ and η are of class C


Elie Awwad
email: elie.awwad@telecom-paristech.fr

Yves Jaouën

Rekaya-Ben Ghaya

Othman

Erwan Pincemin

Polarization-Time Coded OFDM for PDL Mitigation in Long-Haul Optical Transmission Systems

We experimentally demonstrate the potential of Polarization-Time codes in mitigating PDL in long-haul transmissions. The Silver code exhibits the best performance (2dB Q-penalty reduction at a PDL of 6dB). Moreover, for inline PDL, it improves the mean Q-factor by 0.6dB while significantly narrowing the variance of the corresponding distribution.

Introduction

Polarization-multiplexed (PolMux) coherent optical systems combined with DSP techniques has established itself as a key solution for current and future ultra-high bit rate optical transmissions 1 . Among the linear impairments affecting these systems, Polarization Dependent Loss (PDL) was the remaining unsolved performance-degrading effect until very recently, when advanced coding techniques have been suggested 2 . Space-Time (ST) codes have been designed for wireless MIMO communications in order to exploit all the degrees of freedom of these systems. The Golden and Silver codes are the best performing codes designed for MIMO systems with 2 emitting and 2 receiving antennas in a random fading environment. They have been adapted to the optical PolMux channel and are referred to as Polarization-Time (PT) codes. PT codes have proven to be efficient in mitigating PDL [START_REF] Mumtaz | Proc. ICC'10[END_REF][START_REF] Mumtaz | Proc. SPPcom'11[END_REF] . PT coding consists in sending different linear combinations of symbols S 1 ,S 2, S 3, S 4 (i.e., QPSK symbols) on 2 polarizations (pol 1 , pol 2 ) during 2 symbol times (T 1 , T 2 ). The corresponding transmitted codeword matrix X is given by 2 :

  2 1 4 3 2 1 4 4 3 2 1 2 4 3 2 1 3 4 3 2 1 1 pol pol ) S , S , S , S ( f ) S , S , S , S ( f ) S , S , S , S ( f ) S , S , S , S ( f X        T1 T2
The Golden and Silver codes are full-rate, redundancy-free codes, i.e. the transmitted bit rate is the same as for a PolMux system with independent QPSK symbols sent on each polarization, whereas the Alamouti code is a half-rate code. .

Experimental setup

The transmission scheme, shown in Fig. 1, is composed of a PolMux OFDM transmitter, a recirculating loop and a dual-polarization coherent receiver. A total effective bit rate of 25Gbps is achieved, leading to 100Gbps if a multi-band OFDM approach 6 is considered. The OFDM signals transmitted on two orthogonal polarization states consist in 256 subcarriers including 194 data subcarriers and 10 pilot subcarriers for common phase error estimation. In order to compare the performance of the PT codes, we alternatively modulate the data subcarriers with QPSK (No-coding), Silver-, Golden-and Alamouti-coded symbols. To insure the same spectral efficiency for all codes, Alamouti-coded subcarriers use 16-QAM symbols whereas the Golden and Silver codewords are built using QPSK symbols. An 18-sample cyclic prefix is appended to each OFDM symbol to absorb inter-symbolinterference. A 22-symbols training sequence is inserted each 100 OFDM symbols for time synchronization, frequency offset compensation and channel estimation. PT coding correlates the polarization-multiplexed symbols. Hence, four DACs are necessary to generate two complex baseband OFDM signals. Having a unique AWG (Tektronix 7122B) with two outputs, we will consider OFDM symbols satisfying the Hermitian symmetry property. This constraint does not affect our analysis knowing that the complex conjugate of the subcarriers will be treated independently of the original subcarriers at the receiver side. The signals have a ~9GHz bandwidth and the raw bit rate is 36Gbps taking into account the different transmission overheads. Up-conversion to the optical domain is done by two single-drive Mach-Zehnder modulators (MZM). Then the two optical signals are combined by a Polarization Beam Combiner and multiplexed among seven 50GHz-spaced wavelengths. The path difference between the two tributaries of the PolMux OFDM signal is accurately compensated. The transmission line consists in a recirculating loop which contains 2 spans of 100 km of SMF, a PDL element of 2dB, a polarization scrambler and an optical band-pass filter. After additional ASE noise loading at the receiver, the desired wavelength (λ = 1542.9 nm) is filtered and the signal is detected with a dual-polarization coherent receiver. An External cavity laser ECL (Δν =100 kHz) is used as an optical carrier at the transmitter and as a local oscillator at the receiver. A real-time Tektronix oscilloscope samples the received signal at 50 GSa/s. 2.10 6 samples (~10 5 bits/code) are recorded for each channel realization and offline-processing 1,7 Fig. 2: BER curves versus OSNR0.1nm with PDL = 0,3, 6dB for the different schemes: no coding, Silver, Golden and Alamouti codes. ending with ML decoding is carried to measure the BER and the corresponding Q-factor in dB ( and ).

Comparison of OSNR degradation with PDL

Before transmitting the PT-coded OFDM signals in the recirculating loop, we conduct preliminary measurements of the BER evolution as a function of the OSNR for 3 different PDL values at the transmitter: 0, 3 and 6 dB. An attenuator is inserted in one branch of the PolMux transmitter to emulate the worst OSNR degradation of one polarization tributary

5

. The loop is replaced by a polarization scrambler and the signal is then combined with an ASE noise source. The measured BER curves are plotted in Fig. 2. We note firstly that the Silver code performs better than the Golden code and secondly that the performance of the Alamouti code is independent of the PDL value. The observed coding gains have been theoretically explained in 8 . In Fig. 3, for a PDL ranging from 0 to 6dB, we show the Q-factor penalties obtained with QPSK, Silver and Golden codes with regard to a PDL-free case at an OSNR of 13dB. We see that the Silver code performs the best reducing the Q-penalty by 2dB for a PDL = 6dB.

Lumped PDL

After these preliminary tests, we will transmit the PT-coded OFDM in the recirculating loop of Fig. 1 and replace the in-line PDL element by an attenuator having the same insertion loss. PDL is emulated at the transmitter side as in the previous section. BER measurements are carried out after a 1000km-transmission for a span input power of the modulated channel ranging from -15 to 5 dBm. For each input power, up to 50 realizations are considered depending on the target BER. The obtained bellshaped curves are plotted in Fig. 4 for the QPSK and Silver-coded schemes at a PDL of 0, 3 and 6dB. We firstly note that all curves have their optimum operating point at -3 dBm. However in presence of PDL, the uncoded scheme severely suffers (Q = 10.5dB at -3 dBm for PDL = 6dB) Fig. 3: Q-penalty versus PDL at OSNR0.1nm = 13dB with regard to a PDL-free scenario.

while the Silver code completely compensates the penalty at PDL = 3dB and nearly all penalty at PDL = 6dB (Q = 12.3dB at -3dBm). Moreover, an important observation is that the PT-coded OFDM does not induce any extra penalties in the non-linear propagation regime.

Fig. 4:

Comparison of BER evolution versus launched input power after 1000km for the Silver-coded and QPSK schemes (PDL at the transmitter = 0, 3, 6dB).

Distributed in-line PDL

PDL being a distributed effect induced by in-line components along the optical transmission link, we are interested in testing the capability of codes in mitigating in-line PDL. Therefore, we remove the emulated PDL at the transmitter and insert a PDL element of 2dB into the loop. Then, we extract the Q-factor distribution for the different coding schemes by measuring 2000 Q-factor values at the optimum operating point corresponding to an input power of -3dBm. The measured OSNR at the receiver is 12dB after additional ASE noise loading in order to evaluate the complete Q-factor distribution. The obtained Q-factor distributions are shown in Fig. 5. After 5 loops, we measure PDL using the estimated 2x2 channel matrices and obtain a Maxwellian probability distribution function of PDL with a 4.2dB mean (the expected theoretical mean being 4dB [START_REF] Mecozzi | [END_REF] ). The inset in Fig. 5 shows the measured probability distribution and the fitted Maxwellian distribution. We notice that the Q-factor distribution for the QPSK-coded subcarriers has a mean of 11.2dB and a standard deviation (std) of 0.55dB whereas the mean Q-factor of the Silver-and Golden-coded subcarriers is 11.8dB (the same observed value for the PDL-free case when OSNR = 12dB) and the Q-factor distributions are narrower. Note, as well, that the Silver code gives a distribution slightly narrower than the one obtained with the Golden code: std of 0.35dB and 0.37dB respectively. Regarding the Alamouti code, we observe a mean Q-factor of 8.3dB due to the use of 16-QAM symbols instead of QPSK to guarantee the same spectral efficiency. However, its Q-factor distribution is the narrowest of all with a std of 0.24dB. This observation can be explained by the powerful orthogonal structure of the Alamouti codeword matrix that makes its performance independent of the amount of accumulated PDL in the link 8 . The observed small variance of its Q-factor distribution can be explained by the increased sensitivity of 16-QAM modulation to phase noise and non-linear effects.

Conclusions

We have experimentally demonstrated, for the first time, the potential of PT codes to mitigate PDL in long-haul optical links, as well as their resistance to non-linear effects. An enhanced performance evaluated through Q-factor measurements was observed when PT codes are used: a narrower Q-factor distribution is obtained and the mean Q-factor is improved.

Fig. 1 :

 1 Fig. 1: Experimental setup. (ECL: External Cavity Laser, AWG: Arbitrary Waveform Generator, MUX: Multiplexer, AO: Acousto-Optical Modulator, PS: Polarization Scrambler, OBPF: Optical Band-Pass Filter, ASE: Accumulated Spontaneous Emission source, LO: Local Oscillator, OSA: Optical Spectrum Analyzer, OSC: Tektronix 50GS/s Oscilloscope)

Acknowledgements

This work has been partially supported by the Celtic-Plus SASER-SIEGFRIED project.


Elie Awwad
email: awwad@telecom-paristech.fr

Rekaya-Ben Ghaya

Yves Othman

Jaouën
email: jaouen@telecom-paristech.fr

Design Criterion of Polarization-Time Codes for Optical Fiber Channels

Keywords: Polarization Dependent Loss, Error Probability, Polarization-Time Coding, Optical Fiber Communication, MIMO

Coherent detection with Polarization Multiplexing (PolMux) is the most promising technique for future optical fiber transmission systems. However, the optical channel suffers from non-unitary impairments known as Polarization Dependent Loss (PDL). Space-Time coding, originally designed for wireless Rayleigh fading channels, was proven to be capable of mitigating PDL. Coding gains of ST codes were evaluated through simulations and experiments that showed differences in their performance on the optical channel and on the wireless channel. In this paper, we derive an upper bound of the pairwise error probability of an optical channel considering the PDL effect. This upper bound explains the performance of ST codes used to mitigate PDL and yields the design criterion that a code should satisfy in order to completely mitigate PDL.

I. INTRODUCTION

From long-haul links to metropolitan and access networks, optical fibers are today's indispensable carriers of data traffic over the world. They are able to transport light over long distances with a low attenuation. Due to this property and to the THz bandwidth in the near infrared region of the electromagnetic spectrum, binary modulation schemes such as On-Off Keying carrying 1 bit/symbol were first used to send data on fiber-optic systems. However, the rapidly growing Internet traffic urged the search for techniques that increase the capacity of the optical fiber channel. For this purpose, coherent detection started replacing conventional direct detection schemes. The main advantages of coherent detection are the enhanced receiver sensitivity and the introduction of multilevel modulation formats that achieve high spectral efficiency. Moreover, research in Digital Signal Processing (DSP) algorithms and development of high-speed electronics offered the possibility of compensating for fiber impairments in the electrical domain at the receiver rather than optically [START_REF] Taylor | Coherent detection method using DSP for demodulation of signal and subsequent equalization of propagation impairments[END_REF]- [START_REF] Savory | Digital filters for coherent optical receivers[END_REF], allowing the implementation of forward error correction codes and advanced equalization techniques in a DSP-based emitter and receiver.

Apart from the amplitude and the phase of an electromagnetic field, the polarization state is an available degree of freedom that is not yet used. With Polarization Multiplexing (PolMux) [START_REF] Jansen | 16x52.5-Gb/s, 50-GHz spaced, POLMUX-CO-OFDM transmission over 4,160 km of SSMF enabled by MIMO processing[END_REF], [START_REF] Jansen | 10x121.9-Gb/s PDM-OFDM Transmission with 2-b/s/Hz Spectral Efficiency over 1,000 km of SSMF[END_REF], we can double the spectral efficiency of an optical fiber link by sending independent information on two orthogonal polarization states. However, during propagation through the optical channel, the transmitted signal suffers several impairments. We focus on the linear impairments that include dispersive effects: deterministic Chromatic Dispersion (CD) and random Polarization Mode Dispersion (PMD), as well as fiber attenuation and random fading effects known as Polarization Dependent Loss (PDL). To mitigate the linear impairments of the optical channel and recover the transmitted data, a polarization diversity coherent receiver along with DSP algorithms are needed [START_REF] Savory | Digital filters for coherent optical receivers[END_REF]- [START_REF] Jansen | 10x121.9-Gb/s PDM-OFDM Transmission with 2-b/s/Hz Spectral Efficiency over 1,000 km of SSMF[END_REF]. While dispersive effects have been thoroughly studied and completely compensated, mitigation of PDL was still not considered and system margins were left in order to absorb the induced penalties.

A PolMux scheme can be described as a 2 × 2 Multi-Input-Multi-Output (MIMO) system where 2 polarization states are used at the emission and the reception [START_REF] Shieh | Theoretical and experimental study on PMD-supported transmission using polarization diversity in coherent optical OFDM systems[END_REF]. Instead of the spatial dimension (antennas) in a wireless scheme, we use here the polarization dimension. Hence, Space-Time (ST) coding can be applied to mitigate PDL by taking advantage of the diversity of the MIMO configuration. The motivation behind the use of ST codes, or more appropriately Polarization-Time (PT) codes in this context, is the random nature of PDL. Indeed, some optical elements attenuate unequally the polarization states of the transmitted signal. This along with random polarization rotations results in breaking the orthogonality of these states. Therefore, it is more likely to recover the original data symbols by sending linear combinations of these symbols at different time slots and on two different polarizations.

The performance of ST codes such as the Golden code [START_REF] Belfiore | The golden code: a 2x2 fullrate space-time code with nonvanishing determinants[END_REF], Silver code [START_REF] Tirkkonen | Improved MIMO performance with nonorthogonal space-time block codes[END_REF] and Alamouti code [START_REF] Alamouti | A simple transmit diversity technique for wireless communications[END_REF] was evaluated on an optical channel with PDL through simulations [START_REF] Mumtaz | Space-Time Codes for Optical Fiber Communication with Polarization Multiplexing[END_REF] and experimental demonstrations [START_REF] Mumtaz | Experimental Demonstration of PDL Mitigation using Polarization-Time Coding in PDM-OFDM Systems[END_REF]. The Silver code performed better than the Golden code, unlike the case of the wireless Rayleigh fading channel. Besides, it was found that the performance of Alamouti code is independent of the amount of PDL in the link [START_REF] Mumtaz | Alamouti Code against PDL in Polarization Multiplexed Systems[END_REF]. In this work, we aim to explain these observations through error probability calculus. We also check whether the Silver code is the optimal code for an optical channel with PDL or better coding schemes can be found.

Accordingly, the paper is organized as follows: in Section II, we recall the main linear impairments of an optical fiber link, then we present in Section III a mathematical model of an optical channel with PDL that will be used in Section IV to compute an upper bound of the error probability expression. Afterwards, in Section V, we derive the design criterion of codes that minimize the error probability of a transmission on an optical fiber channel with PDL and we investigate in Section VI the performance obtained with the existing codes. Finally, we conclude and set forth the perspectives of our work.

II. OPTICAL FIBER EFFECTS

A. Fiber attenuation & dispersive effects

The loss coefficient of Standard Single Mode Fibers (SSMF) is equal to 0.2dB/km. For long-haul links spanning thousands of kilometers, optical amplifiers are periodically inserted to raise the signal power to its initially injected value at the transmitter. As for chromatic dispersion (CD) and polarization mode dispersion (PMD), both time-domain and frequencydomain equalizations exist [START_REF] Savory | Digital filters for coherent optical receivers[END_REF]. However, a special multicarrier format, Orthogonal Frequency Division Multiplexing (OFDM) has proven to be capable of eliminating all interference caused by dispersive effects using a cyclic prefix [START_REF] Shieh | Theoretical and experimental study on PMD-supported transmission using polarization diversity in coherent optical OFDM systems[END_REF].

B. Polarization Dependent Loss

While the dispersion compensation was thoroughly studied, PDL effects were not being considered until recently. In fact, PDL is mainly introduced by inline components such as amplifiers, isolators, couplers and not the fiber itself. The individual PDL value of each component is kept as low as possible. However, many of these components are found in long-haul optical links leading to a large accumulated PDL. Unlike unitary dispersive effects, PDL causes Optical Signalto-Noise Ratio (OSNR) distortions. The information-carrying polarizations experience unequal loss and crosstalk [START_REF] Duthel | Impact of Polarization Dependent Loss on Coherent POLMUX-NRZ-DQPSK[END_REF].

III. OPTICAL CHANNEL MODEL

In order to proceed with an error probability calculus, we will first mathematically define the studied optical channel.

A. Optical channel with PDL

An overview of previous works on PDL, especially by A. Mecozzi, M. Shtaif [START_REF] Mecozzi | The statistics of polarization-dependent loss in optical communication systems[END_REF] and N. Gisin [START_REF] Gisin | Statistics of Polarization Dependent Loss[END_REF], shows that PDL must be described statistically. However, we consider in our work a simple representation, commonly found in literature [START_REF] Garcia | Feedforward carrier recovery for polarization demultiplexed signals with unequal signal to noise ratios[END_REF], [START_REF] Vassilieva | Interplay between PDL and nonlinear effects in coherent polarization multiplexed systems[END_REF], that takes into account the OSNR inequality between two polarization states as well as their crosstalk. This representation had been used in [START_REF] Mumtaz | Space-Time Codes for Optical Fiber Communication with Polarization Multiplexing[END_REF] to simulate the performance of PT coding in mitigating PDL and an experimental demonstration [START_REF] Mumtaz | Experimental Demonstration of PDL Mitigation using Polarization-Time Coding in PDM-OFDM Systems[END_REF] validated the obtained numerical results. PDL can be described by the following transmission matrix:

H P DL = R α  √ 1 - 0 0 √ 1 +  R -1 α ℎ R α =  cos() -sin() sin() cos()  (1)
The diagonal matrix gives the imbalanced attenuation values of the least and most attenuated polarization states. The rotation matrix R α describes a mismatch between the polarization states of the incident signal and the axes of the PDL component. is defined through Γ dB = 10 10 1+γ 1-γ where Γ dB ≥ 0 is the PDL coefficient in dB and consists of the ratio between the highest and the lowest gains of the link. Hence, 0 ≤ ≤ 1 with = 0 being a PDL-free case and = 1 corresponding to a perfect polarizer. We consider in this work constant values of Γ dB and the angles of the rotation matrices are drawn from a uniform distribution over [0 : 2].

B. PolMux OFDM channel model

In a PolMux OFDM transmission, an OFDM signal is sent on each polarization. At the transmitter, the information symbols are assigned to different subcarriers in the frequency domain then the signal is converted to time domain by an inverse Fast Fourier Transform (iFFT). In order to absorb the interference induced by CD and PMD, a suitable cyclic prefix (CP) is added to the beginning of each OFDM symbol. At the receiver, the CP is removed and an FFT operation converts the signal back into the frequency domain. Despite the redundancy introduced by the CP, this prefix eliminates all Inter-Symbol-Interference and allows to correct any distortion caused by a linear dispersive channel with a single-tap frequency equalizer. Each subcarrier sees a flat channel and the received symbols are given by [START_REF] Shieh | Theoretical and experimental study on PMD-supported transmission using polarization diversity in coherent optical OFDM systems[END_REF]:

Y k,i = exp( i ) exp( D ( k ))H k X k,i + N k,i (2)
where X k,i is the symbol on the th subcarrier of the th OFDM symbol. Y k,i is the received symbol. i is the laser phase noise, D ( k ) is the phase shift owing to chromatic dispersion common to both polarizations. H k is the 2×2 Jones matrix of the link, resulting from a concatenation of PMD and PDL components. N k,i represents the additive noise. We consider next that the CP is larger than the delay spread induced by CD and PMD. Hence, it absorbs all resulting interference and the remaining effect is the non-unitary PDL that we shall regard in the error probability calculus.

IV. DERIVATION OF AN UPPER BOUND OF THE ERROR

PROBABILITY

To get an insight on the behavior of PT codes on an optical link with PDL, we consider a PolMux OFDM transmission and use H P DL defined in (1) as the channel matrix. Therefore, each subcarrier can be seen as a slow fading and frequency non-selective channel:

Y 2×T = H P DL X 2×T + N 2×T (3)
is the temporal length of the code. X is the PT codeword matrix. Y is the received matrix. When = 1, X and Y are vectors and X contains two QAM symbols (No Coding). When > 1, PT coding is used ( = 2 for the codes investigated in this paper). The noise N is modeled as additive white with independent and identically distributed circular complex Gaussian components (0, 2 2 ).

According to the channel model in (3), a maximum likelihood (ML) decoder estimates X ′ with X ′′ as:

X ′′ = X ′ ∈ ‖Y -H P DL X ′ ‖ 2 (4)
where is the set of all possible codewords and ‖ • ‖ is the Euclidian norm. We suppose that we have perfect channel state information, i.e. H P DL is known. The error probability is defined as:

P error = {X ′′ ̸ = X ′ } =  X ′ ∈ {X ′ } {X ′′ ̸ = X ′ |X ′ } (5)
For equiprobable codewords, the error probability can be upper-bounded by [START_REF] Proakis | Digital Communications[END_REF]:

error ≤ 1 card()  X ′ ,X ′′ ∈C, X ′′ ̸ =X ′ (X ′ → X ′′ ) (6)
where card() is the cardinality of and (X ′ → X ′′ ) is the pairwise error probability obtained by averaging the conditional pairwise error probability defined as:

(X ′ → X ′′ |H P DL ) = {‖Y -H P DL X ′′ ‖ 2 ≤ ‖Y -H P DL X ′ ‖ 2 |X ′ , H P DL } (7)
When computing the pairwise error probability, we suppose that X ′ and X ′′ are the only possible codewords in the codeword space. Hence, we can write the pairwise error probability in terms of the Gaussian tail function :

(X ′ → X ′′ |H P DL ) = ( ‖H P DL (X ′ -X ′′ )‖ 2 ) (8) ℎ
() = 1 √ 2  +∞ x exp -t 2 /2 (9)
Using Chernoff's bound () ≤ 1 2 exp -x 2 2 ∀, and averaging over H P DL , we get:

(X ′ → X ′′ ) ≤ E H [exp(- ‖H PDL (X ′ -X ′′ )‖ 2 8 2 )] (10)
where E H [] is the averaging operation over all possible channel realizations. Let X = X ′ -X ′′ be the difference of two codewords and

D =  √ 1 - 0 0 √ 1 +
 , we have:

‖H P DL (X ′ -X ′′ )‖ 2 = ‖H P DL X‖ 2 = ‖R α DR -1 α X‖ 2 = ‖DR -1 α X‖ 2 (11)
Now, let = -. We can write

R β X = U =  ⃗ 1 ⃗ 2  , ⃗ 1,2
being row vectors, we get:

‖DU‖ 2 = ‖⃗ 1 ‖ 2 + ‖⃗ 2 ‖ 2 + (‖⃗ 2 ‖ 2 -‖⃗ 1 ‖ 2 ) (12)
However, ‖⃗ 1 ‖ 2 + ‖⃗ 2 ‖ 2 = ‖U‖ 2 = ‖X‖ 2 because a rotation matrix is unitary. Besides, we develop ‖⃗ 2 ‖ 2 -‖⃗ 1 ‖ 2 and obtain:

‖⃗ 2 ‖ 2 -‖⃗ 1 ‖ 2 = () = (‖⃗ 2 ‖ 2 -‖⃗ 1 ‖ 2 ) cos(2) + 2(⟨⃗ 1 , ⃗ 2 ⟩) sin(2) (13) with X =  ⃗ 1 ⃗ 2
, ⃗ 1,2 being row vectors. Let:

= ‖⃗ 2 ‖ 2 -‖⃗ 1 ‖ 2 (14) = 2(⟨⃗ 1 , ⃗ 2 ⟩) (15)
Substituting ( 12) to ( 15) in ( 10), we get:

(X ′ → X ′′ ) = E H [ (X ′ → X ′′ |H P DL )] ≤ exp(- ‖X‖ 2 8 2 )E H [exp(- 8 2 ( cos(2) + sin(2)))] (16)
The last step consists of averaging over H P DL . Considering constant values of Γ dB , we are left with one random variable that varies uniformly on the interval [0 : 2]. Using the following trigonometric identity:

(2) + (2) =  2 + 2 cos(2 + ) ℎ =  -( b a ) > 0 -( b a ) < 0  (17)
We obtain:

(X ′ → X ′′ ) ≤ exp(- ‖X‖ 2 8 2 ) 1 2  2π+φ φ exp(- 8 2  2 + 2 ()) (18)
The function inside the integral is periodic of period 2 and has an axis of symmetry = in the interval [0 : 2]. Using these two properties, (18) becomes:

(X ′ → X ′′ ) ≤ exp(- ‖X‖ 2 8 2 ) 1  π 0 exp(- 8 2  2 + 2 ()) (19)
Hence,

(X ′ → X ′′ ) ≤ exp(- ‖X‖ 2 8 2 ) 0 ( 8 2  2 + 2 ) (20)
where 0 () is the 0 th order modified Bessel function of the first kind. We can approximate the error probability expression in (20) for high SNR values by using a first order approximation of 0 () when → ∞:

0 () ≈ exp() √ 2 (21)
Then, we get:

exp(- ‖X‖ 2 8 2 ) 0 ( 8 2  2 + 2 ) ≈ exp(- ‖X‖ 2 - √ 2 + 2 8 2 ) (22)
V. DESIGN CRITERION In Fig. 1, we see that 0 () is monotonously increasing for ≥ 0 and hence has its minimum for = 0. This corresponds to and in (20) equal to zero and as a result, the obtained error probability expression will be independent of PDL. Consequently, we define the following design criterion: Proposition 1. A Polarization-Time code completely mitigates PDL if and only if all codeword differences satisfy:

1) = ‖⃗ 2 ‖ 2 -‖⃗ 1 ‖ 2 = 0 and 2) = 2(⟨⃗ 1 , ⃗ 2 ⟩) = 0 .
When the design criterion is met, we recover the performance for two parallel additive white Gaussian channels which is the best performance we can achieve:

error,AW GN ≈ exp -‖X‖ 2 8σ 2 (23)
We clearly see that the resulting criteria are completely different from the rank and the minimum determinant criteria for a Rayleigh fading channel [START_REF] Tarokh | Space-time codes for high data rate wireless communication: performance criteria in the presence of channel estimation errors, mobility, and multiple paths[END_REF], defining respectively a diversity gain and a coding gain. If we compare the approximation of the error probability expression at high SNR (22) to the one obtained in the case of a 2 × 2 MIMO Rayleigh fading channel, we notice different behaviors: the error probability of the Rayleigh fading channel decays as SNR -2r [START_REF] Tarokh | Space-time codes for high data rate wireless communication: performance criteria in the presence of channel estimation errors, mobility, and multiple paths[END_REF] where is the rank of the matrix X, whereas the error probability of the PDL channel decays exponentially as a function of the SNR. Hence, Space-Time codes bring no diversity gain to the optical channel with PDL. Furthermore, if one of the design criterion is not satisfied, (22) points out that the minimum value of ‖X‖ 2 -√ 2 + 2 should be maximized in order to minimize the error probability.

VI. CODE PERFORMANCE ANALYSIS

In this section, we analyze the performance of both coded and uncoded schemes using the derived upper bound of the pairwise error probability expression in (22). All the reported performance in the next subsections match with the performance presented in [START_REF] Mumtaz | Space-Time Codes for Optical Fiber Communication with Polarization Multiplexing[END_REF].

In Table I we use 4-QAM information symbols to build the codewords of the Silver and the Golden codes (full-rate codes with 2 symbols/cu), and 16-QAM symbols for the Alamouti code because of the redundancy in its structure (1 symbol/cu).

We consider an average symbol energy S = 1 for all constellations.

A. Alamouti code

In Fig. 2, we report the performance of Alamouti code for different PDL values (Γ = 3, 6, 10) obtained through Monte Carlo simulations. We notice that the performance of the Alamouti code is independent of the amount of PDL in the link. The orthogonality of the Alamouti codeword matrix induces = = 0 for all possible codeword differences. This explains why it performs the same independently of PDL. However, its performance is affected by the use of 16-QAM symbols giving a squared minimal distance ‖X‖ 2 = 0.8.

B. Silver code

The Silver code is not optimal in mitigating PDL since it does not satisfy the derived design criterion. Unlike the Alamouti code, the Silver code has only some codeword differences having = = 0. In Fig. 3, we plot the performance of the Silver code for different PDL values. We see that the code mitigates almost all PDL when the PDL coefficient is equal to 3 and 6dB. A slight difference in performance is observed for these two PDL values. Whereas for a PDL of 10dB, the code is not able to completely palliate PDL. Looking at

C. Golden code and uncoded scheme

From Fig. 4, we can see that the Silver code outperforms the Golden code for Γ = 6, and both reduce the penalty that PDL causes to the uncoded scheme. Again, this result can be explained by looking at Table I. Indeed, the minimum value of ‖X‖ 2 -√ 2 + 2 is the greatest for the Silver code followed by the Golden code and then the uncoded scheme.

In conclusion, we were able to explain, in terms of error probability bounds, the performance of the Alamouti, Silver and Golden codes when applied on an optical channel with PDL. The Golden and the Silver codes were designed to satisfy the rank and the minimum determinant criteria for a wireless channel that are no more relevant for the optical channel.

VII. CONCLUSION AND PERSPECTIVES

We have established an upper bound of the error probability expression for an optical channel with PDL that allows to predict and compare the performance of any Polarization-Time code used to mitigate PDL by computing a certain distance that depends on the codewords and the PDL value. This expression also yields the design criterion required to construct codes that guarantee a PDL-independent error probability matching the performance for an additive white Gaussian channel. In our future work, we look forward to elaborate the error probability calculus by including the statistics of the PDL coefficient Γ dB that is, under some conditions, Maxwellian distributed [START_REF] Mecozzi | The statistics of polarization-dependent loss in optical communication systems[END_REF] and design specific codes for the optical channel.

Fig. 1 .

 1 Fig. 1. Plot of the 0 th order modified Bessel function of the first kind.

10 Fig. 2 .

 102 Fig. 2. Bit Error Rate as a function of SNR for the Alamouti code, obtained through Monte Carlo simulations.

Γ dB = 6 Fig. 4 .

 64 Fig. 4. Performance comparison of 4-QAM (No coding), Silver, Golden and Alamouti coding schemes when Γ dB = 6.

  , we report the minimum values of ‖X‖ 2 -√ 2 + 2 for each investigated ST code at four different PDL values. To fill the table, we set a spectral efficiency of

		Γ = 0dB	Γ = 3dB	Γ = 6dB	Γ = 10dB
	Silver Code	2	2	2	1.23
	Golden Code	2	1.7	1.46	1.07
	No Coding	2	1.34	0.8	0.38
	Alamouti Code	0.8	0.8	0.8	0.8

  Table I, we find that the minimum value of ‖X‖ 2 -√ 2 + 2 for the Silver code is the same at a PDL of 3 and 6dB, and is equal to 2. Whereas at a PDL of 10dB, this minimum falls to 1.23 which explains the obtained performance. Fig. 3. Bit Error Rate as a function of SNR for the Silver code, obtained through Monte Carlo simulations.

		10 0			Bit Error Rate, Silver Code
		10 -1
		10 -2
		10 -3
	BER	10 -4
		10 -5		Silver Code, Γ dB = 0
				Silver Code, Γ dB	= 3
		10 -6		Silver Code, Γ dB = 6
		10 -7		Silver Code, Γ dB = 10
		10 -8
		0	2	4	6	8	10	12	14	16
					E s /2σ 2

bits per channel use (/cu) for all coding schemes. Hence, TABLE I MINIMUM VALUE OF ‖X‖ 2γ √ a 2 + b 2 FOR DIFFERENT CODING SCHEMES AT DIFFERENT PDL VALUES