Antoine Olivier
email: olivier@ljll.math.upmc.fr

Thomas Haberkorn
email: thomas.haberkorn@univ-orleans.fr

Emmanuel Trélat
email: emmanuel.trelat@upmc.fr

Eric Bourgeois
email: eric.bourgeois@cnes.fr

David-Alexis Handschuh
email: david-alexis.handschuh@cnes.fr.

Éric Bourgeois

Redundancy implies robustness for bang-bang strategies

published or not. The documents may come   L'archive ouverte pluridisciplinaire

Introduction 1.Overview of the method

To introduce the subject, we explain our approach on the control problem consisting of steering the finite-dimensional nonlinear control system ẋ(t) = f (t, x(t), u(t)),

from a given x(0) = x 0 to the target point x(t f ) = x f , with a scalar control u that can only switch between two values, say 0 and 1. The general method, as well as all assumptions, will be written in details in a further section. Let E(x 0 , t f , u) = x(t f ) be the end-point mapping, where x(•) is the solution of (1) starting at x(0) = x 0 and associated with the control u. One aims at finding a bang-bang control u, defined on [0, t f ] for some final time t f > 0, such that E(x 0 , t f , u) = x f . Many problems impose to implement only bang-bang controls, i.e., controls saturating the constraints but not taking any intermediate value. These are problems where only external actions of the kind on/off can be applied to the system.

Of course, such bang-bang controls can usually be designed by using optimal control theory (see [START_REF] Lee | Foundations of optimal control theory[END_REF][START_REF] Pontryagin | The mathematical theory of optimal processes[END_REF][START_REF] Trélat | Contrôle optimal[END_REF]). For instance, solving a minimal time control problem, or a minimal L 1 norm as in [START_REF] Caponigro | Sparse stabilization and optimal control of the Cucker-Smale model[END_REF], is in general a good way to design bang-bang control strategies. However, due to their optimality status, such controls often suffer from a lack of robustness with respect to uncertainties, model 1 errors, deviations from the target. Moreover, when the Pontryagin maximum principle yields bangbang controls, such controls have in general a minimal number of switchings: in dimension 3 for instance, it is proved in [START_REF] Krener | The structure of small-time reachable sets in low dimensions[END_REF][START_REF] Kupka | Geometric theory of extremals in optimal control problems. I. The fold and Maxwell case[END_REF][START_REF] Schättler | On the local structure of time-optimal bang-bang trajectories in R 3[END_REF] (see also [START_REF] Bonnard | Singular trajectories and their role in control theory[END_REF][START_REF] Bonnard | Optimal control of the atmospheric arc of a space shuttle and numerical simulations with multiple-shooting method[END_REF][START_REF] Trélat | Optimal control and applications to aerospace: Some results and challenges[END_REF] for more details on this issue) that, locally, minimal time trajectories of single-input control-affine systems have generically two switchings. Taking into account the free final time, this makes three degrees of freedom, which is the minimal number to generically make the trajectory reach a target point in R 3 , i.e., to solve three (nonlinear) equations.

In these conditions, a natural idea is to add redundancy to such bang-bang strategies, by enforcing the control to switch more times than necessary. These additional switching times are introduced by needle-like variations, as in the classical proof of the Pontryagin maximum principle (see [START_REF] Lee | Foundations of optimal control theory[END_REF][START_REF] Pontryagin | The mathematical theory of optimal processes[END_REF]).

We recall that a needle-like variation π 1 = (t 1 , δt 1 , u 1 ) of a given control u is the perturbation u π1 of the control u given by

u π1 (t) = u 1 if t ∈ [t 1 , t 1 + δt 1 ], u(t) otherwise, (2)
where t 1 ∈ [0, t f ] is the time at which the spike variation is introduced, δt 1 is a real number of small absolute value that stands for the duration of the variation, and u 1 ∈ [0, 1] is some arbitrary element of the set of values of controls. When δt 1 < 0, one replaces the interval [t 1 , t 1 + δt 1 ] with [t 1 + δt 1 , t 1 ] in [START_REF] Pontryagin | The mathematical theory of optimal processes[END_REF]. It is well known that, if |δt 1 | is small enough, the control u π1 is admissible (that is, the associated trajectory solution of ( 1) is well-defined on [0, t f ]) and generates a trajectory x π1 (•), which can be viewed as a perturbation of the nominal trajectory x(•) associated with the control u, and which steers the control system to the final point

E(x 0 , t f , u π1 ) = E(x 0 , t f , u) + |δt 1 | v π1 (t f ) + o(δt 1 ), (3)
where the so-called variation vector v π1 (•) is the solution of some Cauchy problem related to a linearized system along x(•) (see [START_REF] Lee | Foundations of optimal control theory[END_REF][START_REF] Pontryagin | The mathematical theory of optimal processes[END_REF][START_REF] Silva | Smooth regularization of bang-bang optimal control problems[END_REF] and Proposition 1). Recall that the first Pontryagin cone K(t f ) is the smallest closed convex cone containing all variation vectors v π1 (t f ); it serves as a local convex estimate of the set of reachable points at time t f (with initial point x 0 ). Assume that the nominal control u, which steers the system from x 0 to the target point x f , is bang-bang and switches N times between the extreme values 0 and 1 over the time interval [0, t f ]. We denote by T = (t 1 , . . . , t N ) the vector consisting of its switching times 0 < t 1 < • • • < t N < t f . Then the control u can equivalently be represented by the vector T , provided one makes precise the value of u(t) for t ∈ (0, t 1 ). One can also add new switching times: for instance if u(t) = 0 for t ∈ (0, t 1 ), given any s 1 ∈ (0, t 1 ), the needle-like variation π 1 = (s 1 , δs 1 , 1) (with |δs 1 | small enough) is a bang-bang control having two new switching times at s 1 and s 1 + δs [START_REF] Lee | Foundations of optimal control theory[END_REF] .

In what follows, we designate a bang-bang control either by u or by the set T = (t 1 , . . . , t N ) of its switching times. This is with a slight abuse because we should also specify the value of u along the first bang arc. But we will be more precise, rigorous and general in a further section. The endpoint mapping is then reduced to the switching times, and one has E(x 0 , t f , T ) = x f . A variation δT = (δt 1 , . . . , δt N ) of the switching times generates N variation vectors (v 1 (t f ), . . . , v N (t f )), and the corresponding bang-bang trajectory reaches at time t f the point (see Figure 1, where two variations vectors are displayed, for two switching times t 1 and t 2 )

E (x 0 , t f , T + δT ) = x f + δt 1 • v 1 (t f ) + • • • + δt N • v N (t f ) + o( δT ).
Therefore the end-point mapping E is differentiable with respect to T , and Notice that compared to (3), the absolute values disappear. We will prove this result in details further in the paper. In particular, the range of this differential is the first Pontryagin cone K(t f ) (see also [START_REF] Silva | Smooth regularization of bang-bang optimal control problems[END_REF]). Obviously, the more switching times (i.e., degrees of freedom), the more accurate the approximation of the reachable set.

∂E ∂T (x 0 , t f , T ) • δT = δt 1 • v 1 (t f ) + • • • + δt N • v N (t f ). (4) • x 0 t 1 t 2 t 1 + δt 1 t 2 + δt 2 • • x f • δt 1 • v 1 (t f ) δt 2 • v 2 (t f )
We now add redundant switching times (s 1 , . . . , s ) for some ∈ IN in order to generate more degrees of freedom to solve the control problem E (x 0 , t f , (t 1 , . . . , t N , s 1 , . . . , s )) = x f . We order the times in the increasing order and we still denote by T the vector of all switching times.

Redundancy creates robustness. We will see further that these redundant switching times contribute to make the trajectory robust to external disturbances or model uncertainties, we will develop a method to tune the switching times in order to absorb these perturbations and steer the system to the desired target

x f ∈ R n .
Here, in this still informal introduction, we show how to use the additional switching times to make the system reach targets x f + δx f in a neighborhood of x f . The idea is to solve the nonlinear system of equations

E(x 0 , t f , T + δT ) = x f + δx f .
Using (4), we propose to solve, at the first order,

∂E ∂T (x 0 , t f , T ) • δT = δx f , (5)
which makes n equations with N + degrees of freedom. We assume that N + is (possibly much) larger than n and that the matrix in [START_REF] Krener | The structure of small-time reachable sets in low dimensions[END_REF] is surjective. Then one can solve (5) by using the Moore-Penrose pseudo-inverse ∂E ∂T † of ∂E ∂T (see [START_REF] Golub | Matrix computations[END_REF], or see [START_REF] Beutler | The operator theory of the pseudo-inverse. I. Bounded operators[END_REF][START_REF] Beutler | The operator theory of the pseudo-inverse. II. Unbounded operators with arbitrary range[END_REF] for a theory in infinite dimension), which yields the solution of minimal Euclidean norm

δT = ∂E ∂T † • δx f ,
and we have

δT 2 δx f 2 σ min , (6)
where σ min is the smallest positive singular value of ∂E ∂T . This estimate gives a natural measure for robustness, that we will generalize.

The two main contributions of this paper are:

• the idea of adding redundant switching times in order to make a nominal bang-bang control more robust, while keeping it as being bang-bang;

• the design of a practical tracking algorithm, consisting of solving an overdetermined nonlinear system by least-squares, thus identifying a robustness criterion that we optimize.

They are developed in a rigorous and general context in the core of the paper.

State of the art on robust control design

There is an immense literature on robust control theory, with many existing methods in order to efficiently control a system subjected to uncertainties and disturbances. Whereas there are many papers on H 2 and H ∞ methods, except a few contributions in specific contexts, we are not aware of any general theory allowing one to tackle perturbations by using only bang-bang controls. This is the focus of this paper.

Let us however shortly report on robustness methods when one is not bound to design bangbang controls. In [START_REF] Koh | A smooth path tracking algorithm for wheeled mobile robots with dynamic constraints[END_REF], a path-tracking algorithm with bang-bang controls is studied, for a double integrator and a wheeled robot. The technique relies heavily on the expression of the equations and does not apply to more general systems. In [START_REF] Singh | Robust time-optimal control -Frequency domain approach[END_REF], the authors build a robust minimal time control for spacecraft's attitude maneuvers by canceling the poles of some transfer function. A remarkable fact is that the robustified control presents more switchings than the minimal time control. In this case, the robustness is evaluated as the maximum amplitude on a Bode diagram (see also [START_REF] Qiang | Robust time-optimal control of uncertain flexible spacecraft[END_REF] and [START_REF] Bong | Robust time-optimal control of uncertain structural dynamic systems[END_REF] for similar works). In [START_REF] You | Robust, near time-optimal control of nonlinear second order systems with model uncertainty[END_REF], the authors observe that bang-bang controls are intrinsically not robust, and use pieces of singular trajectories (hence, not bang-bang) to overcome this issue.

In the H 2 and H ∞ theories, control systems are often written in the frequency domain using the Laplace transform. For a transfer matrix G(s), the two classical measures for performance are (see [START_REF] Doyle | State-space solutions to standard h 2 and h ∞ control problems[END_REF][START_REF] Zhou | Robust and Optimal Control[END_REF]) the H 2 norm and the H ∞ norm respectively:

G 2 = 1 2π +∞ -∞ Trace(G(jω)G(jω) * )dω 1/2 and G ∞ = sup ω∈R σ(G(jω)),
where σ(G) is the largest singular value of G.

In the linear quadratic theory, the question of optimal tracking has been widely addressed: given a reference trajectory ξ(•), we track it with a solution of some control system ẋ(t) = f (x(t), u(t)), minimizing a cost of the form

t f 0 x(t) -ξ(t) 2 W + u(t) 2 U dt + x(t f ) -ξ(t f ) 2 Q ,
with weighted norms (see [START_REF] Brian | Linear optimal control[END_REF][START_REF] Kwakernaak | Linear optimal control systems[END_REF][START_REF] Trélat | Contrôle optimal[END_REF]). The first term in the integral measures how close one is to the reference trajectory, the second one measures a L 2 norm of the control (energy), and the third one accounts for the distance at final time between the reference trajectory ξ(•) and x(•).

Then, the control can be expressed as a feedback function of the error x(t)ξ(t), involving the solution of some Riccati equation. In [START_REF] Novel | Control Theory for Engineers: A Primer[END_REF][START_REF] Hassan | Nonlinear systems[END_REF], the authors investigate the question of stabilizing around a slowly time-varying trajectory. They also introduce uncertainties on the model and study the sensitivity of the system to those uncertainties. In the case of the existence of a delay on the input, a feedback law is proposed. In [START_REF] Lin | Robust Control Design: An Optimal Control Approach[END_REF][START_REF] Tan | An optimal control approach to robust tracking of linear systems[END_REF], uncertainties p are introduced in a linear system ẋ(t) = A(p)x(t) + Bu(t), and a tracking algorithm is suggested, under matching conditions on the uncertainties or not (see also [START_REF] Abdallah | Survey of robust control for rigid robots[END_REF] for a survey on robust control for rigid robots).

In the late 1970's, H ∞ control theory developed. The control system is often described by a plant G and a controller K. Then, the dependency of the error z (to be minimized) on the input v can be written as z = F (G, K)v. The H ∞ control problem consists of finding the best controller K such that the H ∞ norm of the matrix F (G, K) is minimized:

F (G, K) ∞ = sup ω∈R σ(F (G, K)(jω)).
It can be interpreted as the maximum gain from the input v to the output z. This criterion was introduced in order to deal with uncertainties on the model (on the plant G). In [START_REF] Zames | Feedback and optimal sensitivity: Model reference transformations, multiplicative seminorms, and approximate inverses[END_REF], the author introduced the notion and highlighted the connection with robustness. In [START_REF] Doyle | State-space solutions to standard h 2 and h ∞ control problems[END_REF], a link is shown between the existence of such a controller and conditions on the solutions of two Riccati equations. Following a notion introduced in [START_REF] Gahinet | A convex parametrization of h ∞ suboptimal controllers[END_REF], the linear matrix inequality (LMI) approach was introduced in [START_REF] Gahinet | A linear matrix inequality approach to h ∞ control[END_REF], and used in [START_REF] Apkarian | A Spectral Quadratic-SDP Method with Applications to Fixed-Order H 2 and H ∞ Synthesis[END_REF][START_REF] Apkarian | Nonsmooth h ∞ synthesis[END_REF] to solve the H ∞ synthesis. The Riccati equations are replaced with Riccati inequalities, whose set of solutions parameterizes the H ∞ controllers (see also [START_REF] Boyd | Linear Matrix Inequalities in System and Control Theory[END_REF] for the use of LMIs in control theory). The papers [START_REF] Doyle | Multivariable feedback design: Concepts for a classical/modern synthesis[END_REF][START_REF] Mcfarlane | A loop-shaping design procedure using h infinity synthesis[END_REF][START_REF] Xie | Robust h infin; control for linear systems with norm-bounded time-varying uncertainty[END_REF] present design procedures in this context to elaborate the feedback controller K. In [START_REF] Ge | Robust H ∞ state feedback control for linear systems with state delay and parameter uncertainty[END_REF], the theory is extended to systems with parameters uncertainties and state delays, as well as in [START_REF] Xu | Robust stochastic stabilization and control of uncertain neutral stochastic time-delay systems[END_REF], with stochastic uncertainty.

In many optimal control problems, the application of the Pontryagin maximum principle leads to bang-bang control strategies, and the classical H 2 and H ∞ theories were not designed for such a purpose. But the optimal trajectories are in general not robust. Adding needle-like variations is therefore a way to improve robustness, and is the main motivation of this paper. Of course, the method applies to any bang-bang control strategy, not necessarily optimal.

The approach that we suggest in this paper combines an off-line treatment of the control strategies, with a feedback algorithm based on the structure of the control. We emphasize here that this algorithm preserves the bang-bang structure of the control. It consists of applying a nominal control strategy (that needs to be computed a priori ), and adjusting it in real time, allowing one to track a nominal trajectory. The off-line method takes a solution of the control problem and makes it more robust by adding additional switching times (i.e., redundancy), which can be seen as additional degrees of freedom. Note that our analysis is done in the state space, without needing to consider the frequency domain. A key ingredient to the method is the use of needle-like variations.

Structure of the paper

The paper is organized as follows. In Section 2, we develop an algorithm to steer a perturbed system to the desired final point. The method is similar to the one presented in Section 1.1, except that we need to consider a backward problem. Indeed, the final point is fixed, and perturbations appear all along the trajectory. Besides, our measure for robustness comes out naturally in view of [START_REF] Kupka | Geometric theory of extremals in optimal control problems. I. The fold and Maxwell case[END_REF]. Having identified the robustness criterion, we show in Section 3 how to add redundant switching times, leading one to solve a finite-dimensional nonlinear optimization problem. In Section 4, we provide some numerical illustrations on the attitude control problem of a 3-dimensional rigid body.

Tracking algorithm

Setting. In this paper, we consider the control system

ẋ(t) = f (t, x(t), u(t)), (7)
where

f is a smooth function R × R n × R m → R n , the state x(•) ∈ R n , the control u(•) ∈ L ∞ ([0, t f ]; Ω), and Ω is the subset of R m : [a 1 , b 1 ] × • • • × [a m , b m ].
We make two additional hypothesis: the controls we consider are "bang-bang", with a finite number of switching times:

(H 1 ) ∀i ∈ 1, m , u i (t) ∈ {a i , b i }, a.e. (H 2 ) ∀i ∈ 1, m , u i does not chatter.
A control is chattering when it switches infinitely many times over a compact time interval (see [START_REF] Zhu | Minimum time control of the rocket attitude reorientation associated with orbit dynamics[END_REF][START_REF] Fuller | Study of an Optimum Non-linear Control System[END_REF]). Therefore, our method does not apply to those controls. However, when the solution of an optimal control problem chatters, provided that it is possible, one could consider a sub-optimal solution, with only a finite number of switching times.

In the context of optimal control, we will denote the cost under the form

C(u) = t f 0 f 0 (t, x(t), u(t)) dt. (8)
We recalled in the introduction the (classical) definitions of the end-point mapping, of a needlelike variation (2) and the expansion of the end-point mapping subject to a needle-like variation (3).

Reduced end-point mapping

In this subsection, we give the definition of the reduced end-point mapping, and show a differentiability property.

Let us consider a bang-bang control u(•), and its associated trajectory x(•). For the sake of simplicity, we make the additional assumption that for every switching time t j , one and only one component of the control commutes. Therefore, provided we specify the initial value of each component, the control u is entirely characterized by the switching times of its components and can be represented by a vector:

((u 10 , . . . , u m0 ), (t 1 , i 1 ) , . . . , (t N , i N ) , t f ) ∈ Ω × R 2N +1 , where u i0 ∈ {a i , b i } is the initial value for the control u i (•) (i ∈ 1, m ), N
is the total number of switching times, t f is the final time, and i j is the component of the control that switches at time t j . As this representation entirely characterizes the control, we will use indistinctly the notation u and ((u 10 , . . . , u m0 ), (t 1 , i 1 ) , . . . , (t N , i N ) , t f ) to speak about the control whose components switch at the times t j . In the literature, ((t 1 , i 1 ) , . . . , (t N , i N )) is often called a switching sequence.

Remark 1. Had we wanted to allow simultaneous switching of multiple components, we would need to consider controls represented by: ((u 10 , . . . , u m0 ), (t 1 , I 1 ) , . . . , (t N , I N ) , t f ) , where I j ⊂ 1, m represents the set of components that switch at time t j .

Definition 1 (Reduced end-point mapping). We define the reduced end-point mapping by

E(x 0 , (u 10 , . . . , u m0 ), (t 1 , i 1 ) , . . . , (t N , i N ) , t f ) = x u (x 0 , t f ),
where u is the control represented by ((u 10 , . . . , u m0 ), (t 1 , i 1 ) , . . . , (t N , i N ) , t f ), and x u (x 0 , t f ) is the associated state at time t f , starting at x 0 .

Note that in [START_REF] Maurer | Optimization methods for the verification of second order sufficient conditions for bangbang controls[END_REF][START_REF] Maurer | Second Order Sufficient Conditions for Time-Optimal Bang-Bang Control[END_REF], the authors also reduce a bang-bang control to its switching points, in order to formulate an optimization problem in finite-dimension.

In the following, when writing this reduced end-point mapping, we may consider that the initial point x 0 is fixed, as well as the way the components of the control switch (i.e., we consider that the N-tuple (i 1 , . . . , i N ) is fixed), the initial values u i0 and the final time t f . In this context, we may forget them in the notations, and denote the reduced end-point mapping by

E(t 1 , . . . , t N ) = x u (t f ).
A remarkable fact is that the reduced end-point-mapping is differentiable. Compared to the expansion (3) with respect to a needle-like variation, the sign of δt does not matter. For the sake of completeness, we give the proof in appendix.

Proposition 1. The reduced end-point mapping is differentiable, and

dE(t 1 , . . . , t N ) = v 1 (t f ) • • • v N (t f ) ∈ M n,N (R),
where v j (•) (j ∈ 1, N ) is the solution of the Cauchy problem, defined for t t j :

vj (t) = ∂f ∂x (t, x(t), u(t))v j (t) v j (t j ) = f (t j , x(t j ), (. . . , a ij , . . .)) -f (t j , x(t j ), u(t + j )) if u ij switches from a ij to b ij . f (t j , x(t j ), (. . . , b ij , . . .)) -f (t j , x(t j ), u(t + j )) if u ij switches from b ij to a ij .
The notation (. . . , a ij , . . .) (resp. (. . . , b ij , . . .)) is used to show a difference with u(t + j ) (resp. u(t - j )) on the i j -th component only.

Remark 2. In the special (and important in practice) case of a control-affine system

ẋ(t) = f 0 (x(t)) + m k=1 u k (t)f k (x(t)),
the initial condition on v j can be written much more easily:

v j (t j ) = (a ij -b ij )f ij (x(t j )) if u ij switches from a ij to b ij . (b ij -a ij )f ij (x(t j )) if u ij switches from b ij to a ij .

Absorbing perturbations

As explained in the introduction, we present in this paper a closed-loop method to actually steer the system towards a point x f , with bang-bang controls, even in the presence of perturbations. First, for the sake of simplicity, we will explain how to control the system to some point x f +δx f . We will see that this idea can be adapted for our purpose of controlling a perturbed trajectory, by simply reversing the time.

Perturbations on the final point. We briefly generalize the problem introduced in the introduction. Let u = (u 10 , . . . , u m0 ),

t 1 , i 1 , . . . , t N , i N , t f ∈ Ω × R 2N +1
be a control such that x u (t f ) = x f . That is, using the definition of Subsection 2.1, we have that

E(x 0 , (u 10 , . . . , u m0 ), t 1 , i 1 , . . . , t N , i N , t f ) = x f .
Or, considering that the final time t f , the initial point x 0 , the components (i 1 , . . . , i N ) and the initial values (u 10 , . . . , u m0 ) are fixed,

E(t 1 , . . . , t N ) = x f .
Let δx f be some perturbation of the final point x f . We look for a vector δT = (δt 1 , . . . , δt N ) so that the system reaches the target point x f + δx f :

E(t 1 + δt 1 , . . . , t N + δt N ) = x f + δx f .
As we have shown in Proposition 1 the differentiability of the reduced end-point mapping, we can write

E(t 1 + δt 1 , . . . , t N + δt N ) = E(t 1 , . . . , t N ) + dE(t 1 , . . . , t N ) • δT + o( δT ).
At order one, the solution is given by the solution of the linear equation

dE(t 1 , . . . , t N ) • δT = δx f .
It is natural to target the final point x f +δx f while shifting the switching times as little as possible.

That is, we look for the solution of minimal euclidian norm of the previous equation, which is given by δT = dE(t 1 , . . . , t N ) † • δx f . Therefore, we have shown how to compute, at order one, the correction to apply to control the system to some point x f + δx f : it boils down to solving a least-squares problem. Let us keep in mind that our definitive goal is to control systems that are perturbed all along their trajectory, to a fixed final point x f . In other words, from a perturbed point x(t) + δx(t) at some time t ∈ [0, t f ), we want to absorb the perturbation δx(t) and still reach the final point x f . Even if this is a slightly different setting, we show that we can apply the same idea if we look at a backward problem.

Absorbing a perturbation at time t. Let (x(•), u(•)) be a nominal solution of the control system [START_REF] Schättler | On the local structure of time-optimal bang-bang trajectories in R 3[END_REF]. We assume that when applying in practice the control u = T , because of model uncertainties and perturbations, we observe a perturbed trajectory x per (t) = x(t) + δx(t).

Let t ∈ [0, t f ]. Starting from the perturbed point x(t) + δx(t), which stands as a new initial point, we want to reach the final point x f in time t ft. Hence, we look for a control u + δu such that

E(x(t) + δx(t), u + δu, t f -t) = x f .
Assume for a moment that the perturbation of the control δu is small in L ∞ norm. Then, at least formally, one can write

E(x(t), u, t f -t) + ∂E ∂x 0 (x(t), u, t f -t) • δx(t) + ∂E ∂u (x(t), u, t f -t) • δu + o( δx(t) + δu ) = x f .
Therefore, at order one, we look for a solution of the (linear) equation

∂E ∂x 0 (x(t), u, t f -t) • δx(t) + ∂E ∂u (x(t), u, t f -t) • δu = 0. (9)
However, we do not want, in this paper, to apply small perturbations in the L ∞ norm, as they would not result in bang-bang controls (However, this is similar to what is done while performing a Ricatti procedure to stabilize a system or track a reference trajectory). Nevertheless, reducing the end-point mapping to the switching times enables us to preserve the bang-bang structure: in the formalism previously introduced, we need to solve the nonlinear system of equations

E(x(t) + δx(t), T + δT , t f -t) = x f .
The equation [START_REF] Bonnard | Optimal control of the atmospheric arc of a space shuttle and numerical simulations with multiple-shooting method[END_REF] becomes

∂E ∂T (x(t), T , t f -t) • δT = - ∂E ∂x 0 (x(t), T , t f -t) • δx(t), (10)
where the expression ∂E/∂T is given by Proposition 1.

A backward problem. Solving this equation requires the computation of the partial differential ∂E/∂x 0 at the initial point x(t). We will see now that it can be overcome by introducing a backward problem. Of course, the two formulations are equivalent.

Definition 2 (Backward end-point mapping). Let u = (t 1 , . . . , t N ) be a bang-bang control, and t ∈ [0, t f ]. We define the backward end-point mapping by

Ẽ(t, t 1 , . . . , t N ) = x(t f -t),
where x(•) is the solution to the Cauchy problem

ẋ(t) = -f (t f -t, x(t), u(t f -t)), x(0) = x f .
Note that for the nominal trajectory (x(•), u(•)), we have that

Ẽ(t, t 1 , . . . , t N ) = x(t).
Indeed, we have in this case that x(t) = x(t ft): if we integrate the nominal system backward, starting from the point x f during a time period t ft, we end up at point x(t).

Remark 3. Let t ∈ [0, t f ], and j be the smallest index such that t j > t (with the convention that j = N + 1 if t > t N ). Then, note that t 1 , . . . , t j-1 do not play any role in the computation of Ẽ(t, t 1 , . . . , t N ). The differential of Ẽ can be computed with the Proposition 1. It is a matrix of size n × (Nj + 1).

In this context, the problem of adjusting the system back towards x f writes: at time t, find (t j , . . . , t N ) such that Ẽ(t, t 1 , . . . , t N ) = x per (t). [START_REF] Silva | Smooth regularization of bang-bang optimal control problems[END_REF] We see that reversing the time, we place ourselves in the setting previously described of aiming at a perturbed final point. Therefore, we have the following proposition. Proposition 2. At order one in δx, the solution of minimal norm of the problem (11) is given by

T + δT , with δT = d Ẽ(t, T ) † • δx(t), (12)
where d Ẽ(t, T ) † denotes the pseudo-inverse of d Ẽ(t, T ). Moreover, we have the estimate

δT 2 1 σ min (t) δx(t) 2 , ( 13
)
where σ min (t) is the smallest positive singular value of d Ẽ(t, T ).

Proof. The scheme of the proof has already been exposed previously in the paper. However, we write it extensively here. Let δT = T -T . The problem writes Ẽ(t, T + δT ) = x per (t).

According to Proposition 1, the backward end-point mapping is differentiable (and we also know how to compute its derivative), so

Ẽ(t, T + δT ) = Ẽ(t, T ) + d Ẽ(t, T ) • δT + o( δT ) = x(t) + d Ẽ(t, T ) • δT + o( δT ).
So, at order one, the problem writes

d Ẽ(t, T ) • δT = δx(t). (14)
It is well known (see [START_REF] Allaire | Algèbre linéaire numérique[END_REF] for instance), that the solution of minimal norm of this equation is

δT = d Ẽ(t, T ) † • δx(t). Besides, let σ max (t) > • • • > σ min (t) > 0 denote the positive singular values of d Ẽ(t, T ). We have that d Ẽ(t, T ) † 2 = 1/σ min (t) ( • 2
for a matrix denotes the induced norm corresponding to the euclidean norm), so that

δT 2 = d Ẽ(t, T ) † • δx(t) 2 d Ẽ(t, T ) † 2 • δx(t) 2 δx(t) 2 σ min ,
which concludes the proof.

Remark 4. We have the relation that, for all vector of switching times T

E( Ẽ(t, T ), T , t f -t) = x f .
Differentiating this equality with respect to T , we have that, for all δT

∂E ∂x 0 ( Ẽ(t, T ), T , t f -t) • d Ẽ(t, T ) • δT + ∂E ∂T ( Ẽ(t, T ), T , t f -t) • δT = 0.
Replacing the second term by its value in [START_REF] Trélat | Optimal control and applications to aerospace: Some results and challenges[END_REF], it follows that

∂E ∂x 0 ( Ẽ(t, T ), T , t f -t) • d Ẽ(t, T ) • δT = ∂E ∂x 0 ( Ẽ(t, T ), T , t f -t) • δx(t).
It is easy to show that ∂E/∂x 0 can be expressed as the resolvent of a linearized system. Therefore, the matrix ∂E/∂x 0 is invertible, and the equations ( 10) and ( 14) are equivalent. But solving [START_REF] Beutler | The operator theory of the pseudo-inverse. II. Unbounded operators with arbitrary range[END_REF] only requires to compute the derivative of Ẽ. This is what we do in the following.

Remark 5. Note that it might not always be possible to find a solution to the equation d Ẽ(t, T ) • δT = δx(t). This may happen for instance if t > t N -n+1 , i.e., we do not have enough degrees of freedom left to absorb the perturbation δx(t) ∈ R n . However, we can still give a meaning to the equation d Ẽ • δT = δx(t). We look for a solution of the least-square problem:

min δT ∈R N d Ẽ(t, t 1 , . . . , t N ) • δT -δx(t) 2 2
, for which δT = d Ẽ(t, t 1 , . . . , t N ) † • δx(t) is still the solution of minimal norm (see [START_REF] Allaire | Algèbre linéaire numérique[END_REF]). We see here emerging the idea that the number of switching times (i.e., degree of freedom) left at time t, is going to be an important factor to track the system back towards the final point x f . Numerical algorithm. At time t, Equation ( 12) provides us with a formula to adjust the control so that the perturbed trajectory eventually reaches x f . But it certainly does not enable us to face perturbations that would happen after time t. In order to absorb perturbations all along the trajectory, we suggest the following algorithm: Let T be an initial control. Given an integer s and a subdivision 0 < τ 1 < • • • < τ s < t f of the interval [0, t f ], we adjust the control at each τ i for all i ∈ 1, s . That is, for each i ∈ 1, s , we measure the drift δx(τ i ) = x per (τ i )x ref (τ i ), and compute the differential of the backward end-point mapping d Ẽ(τ i , t 1 , . . . , t N ). We deduce from [START_REF] Golub | Matrix computations[END_REF] that the correction to apply is then δT = d Ẽ(τ i , t 1 , . . . , t N ) † • δx(τ i ). We then update the control by considering the new vector of switching times T + δT . Remark 6. When computing the correction T + δT , it may happen that the new switching times are not ordered, i.e., there exists some integer j ∈ 1, N -1 such that t j+1 < t j . In this case, we consider that the correction is not physically acceptable, and we reject it. (Note that in some cases, we may want to continue the integration of the system even if two switching times are not ordered. In that case, we can always use the last admissible control, where all the switching times are ordered.) Remark 7. The computation of the differential d Ẽ(t, t 1 , . . . , t N ) is done via the integration of a system of ordinary differential equations, which can be done efficiently and quickly using numerical integrators. However, the size of the system (as well as the time required to compute the pseudoinverse) directly depends on the number of switching times N and on the state dimension n.

Promoting robustness

Intuitively, we want to say that a control is robust whenever the correction δT required to absorb the perturbation δx(t) is small. Since we have shown the estimate δT 2 δx(t) 2 /σ min (t), a robust trajectory is then one for which the values of 1/σ min (t) remain small along the trajectory. Definition 3. We define the following cost, that we will use to characterize the robustness of a trajectory

C r (t 1 , . . . , t N ) = t N 0 1 σ min (t) 2 dt. ( 15
)
Remark 8. In the previous definition, the upper bound in the integral is t N , because for t > t N , the backward end-point mapping derivative d Ẽ(t, t 1 , . . . , t N ) is not defined, and neither is σ min (t).

For some reason, we may only want to have robustness up until some time t < t N . Then the previous definition would become t 0 1/σ min (t) 2 dt. In this section, we show how the switching times of a trajectory can be chosen to build one that is more robust. We also suggest a new way to design a trajectory, by adding redundant switching times, that give us more degrees of freedom. Note also that we will start from a solution of an optimal control problem, because it is of high importance in practice, but the method generally applies when starting from any control, as long as it satisfies the hypothesis (H 1 ) and (H 2 ). Starting from an initial control such that E(t 1 , . . . , t N ) = x f , we look for redundant switching times (s 1 , . . . , s l ) such that E(t 1 , . . . , t N , s 1 , . . . , s l ) = x f , while minimizing the cost (15) that accounts for robustness:

C r (t 1 , . . . , t N , s 1 , . . . , s l ).

An auxiliary optimization problem

Let us consider a bang-bang trajectory (satisfying the hypothesis (H 1 ) and (H 2 )) of the control system [START_REF] Schättler | On the local structure of time-optimal bang-bang trajectories in R 3[END_REF], optimal for the cost [START_REF] Bonnard | Singular trajectories and their role in control theory[END_REF]. That is, u = ((u 10 , . . . , u m0 ), t

Let us emphasize the fact that reducing the control to its switching times enables us to reduce a problem in infinite dimension

min u∈L ∞ ([0,t f ];Ω) C(u) s.t. E(u) = x f
to a finite number of non-linear problems under non-linear constraints in finite dimension, provided we set N , as we left aside chattering trajectories.

In order to make the control more robust we suggest to solve the following problem. We fix the components of the control (i 1 , . . . , i N ), and we introduce the cost that accounts for the robustness of a trajectory:

min (t1,...,t N ) λ 1 C(t 1 , . . . , t N ) + λ 2 C r (t 1 , . . . , t N ), s.t. E(t 1 , . . . , t N ) = x f
where λ 1 and λ 2 are two parameters, chosen to give more or less importance to the different costs. For instance, if λ 1 λ 2 , the solution is close to the initial one (t 1 , . . . , t N ).

Redundancy creates robustness

Let us consider a control u = ((u 10 , . . . , u m0 ), (t 1 , i 1 ) , . . . , (t N , i N ) , t f ). In order to reduce the optimization space, we will consider in the following subsection that the initial control values (u 10 , . . . , u m0 ), the components (i 1 , . . . , i N ) and the final time t f are fixed, so we will forget them in the notations.

We propose here to go further in order to improve the robustness of the corresponding trajectory. We do so by adding needles to some components of the control. By needle, we mean a short impulse on one of the control. Let us denote by l the number of needles we are willing to add. It means that we look for additional switching times [(s 1 , s 2 ), . . . , (s 2l-1 , s 2l )] and components of the control (j 1 , . . . , j l ), so that for all i ∈ 1, l , (s 2i-1 , s 2i ) are switching times for the j i -th components of the control (see Figure 2). It aims at giving us more degrees of freedom while trying to absorb perturbations δx by moving the switching times (T , S) = (t 1 , . . . , t N , (s 1 , s 2 ), . . . , (s 2l-1 , s 2l )). Thus, we are solving the optimization problem min (j1,...,j l ) min (T ,S)

λ 1 C(T , S) + λ 2 C r (T , S). s.t. E(T , S) = x f ( 17
)
Remark 9. If the original bang-bang control strategy ū does not come from an optimization process, that is there is no cost C associated with it, we can still consider problem [START_REF] Qiang | Robust time-optimal control of uncertain flexible spacecraft[END_REF] but with λ 1 = 0.

Let us denote by T the solution of problem ( 16), and by (T , S ) the solution of problem [START_REF] Qiang | Robust time-optimal control of uncertain flexible spacecraft[END_REF]. Then, we have that C(T ) C(T , S ).

It means that the solution (T , S ) is sub-optimal with respect to the initial cost C. However, this sub-optimality comes with a gain in terms of robustness. Besides, the loss of optimality (and therefore gain in robustness) can be controlled by the choice of the coefficients λ 1 and λ 2 . This problem is a mixed problem, with integer variables (the components (j 1 , . . . , j l )), and continuous variables (the switching times (t 1 , . . . , t N , (s 1 , s 2 ), . . . , (s 2l-1 , s 2l ))). However, if the components are fixed, we only have to solve a non-linear problem subject to non-linear constraints in finite dimension min (T ,S)

λ 1 C(T , S) + λ 2 C r (T , S). s.t. E(T , S) = x f ( 18
)
We used an interior-point algorithm to solve [START_REF] Bong | Robust time-optimal control of uncertain structural dynamic systems[END_REF]. In [START_REF] Xu | Optimal control of switched systems: new results and open problems[END_REF][START_REF] Zhu | Optimal control of hybrid switched systems: A brief survey[END_REF], gradient-based algorithms are shown to be effective to solve such problems, when the sequence of indices (j 1 , . . . , j l ) is fixed. Therefore a "naïve" way to proceed, if m denotes the number of components of the control, is to solve m l optimization problems, which is extremely costly if m or l is big. A compromise has to be found between the potential benefit in robustness and the computational cost. Such a compromise will however depend on the particular problem at hand, so we do not elaborate too much on this issue and give an example in Section 4. Let us cite [START_REF] Chyba | Design and implementation of time efficient trajectories for autonomous underwater vehicles[END_REF][START_REF] Chyba | Increasing underwater vehicle autonomy by reducing energy consumption[END_REF], where the authors parametrize an optimal control problem (for the time-minimal and L 1 problem) with the switching times of the controls. They simplify its complex structure by fixing the number of switching times, and wonder how many switching times are required to obtain a cost close to the optimal one : the result is striking as 2 or 3 may be enough. However, they know from an a priori study the value of the optimal L 1 or time-minimal cost, and therefore can stop adding switching times when reaching a given percentage of this optimal value of the criterion. In our problem, we do not know what is the optimal value of the criterion we identified to quantify the robustness of a trajectory. It becomes necessary to find another way to decide how many needles to add. One could consider tackling directly Problem 17, a combinatorial optimization problem (which is a class of problem known to be hard to solve). Recent years have seen the development of advanced numerical procedures to deal with the combinatorial nature of those problem at a reasonable computational cost. We give more details on this issue at the end of this section.

Remark 10. Let us make here a remark on the ordering of the switching times. In the vector (T , S) are stored the switching times t i and s i that represent the control u. Those swicthing times are not necessarily ordered during or after the optimization process, so let T = (τ 1 , . . . , τ N +2l ) denote the ordered equivalent to (T , S). So far, we have made the implicit assumption that when we perform the numerical integration of the system, the switching times are ordered: τ i+1τ i 0 for all i ∈ 0, N + 2l -1 . We recall that our goal is to absorb perturbations δx. As explained in Subsection 2.2, we compute at order one the correction to apply δT = dE(T) † • δx. At this point, we could have that T + δT does not satisfy this ordering property. Then, we consider that T + δT is not admissible, and an estimate like [START_REF] Beutler | The operator theory of the pseudo-inverse. I. Bounded operators[END_REF] would not hold.

In the following, in order to guarantee that we do not have an interchanging of the switching times (at least for small perturbations), we add an additional constraint whilst elaborating the robustified trajectory (u(•), x(•)) at (17):

τ i+1 -τ i η for all i ∈ 0, N + 2l -1 , (19)
for some η > 0, where T = (τ 1 , . . . , τ N +2l ) denotes the re-ordering of the vector (T , S). In that way, we ensure that two consecutive switching times (T and S combined) are at least distant of η. Thus, if δx is small enough, the elements of the vector T + dE(T) † • δx remain in ascending order. Besides, such a constraint is often highly justified in practice, for instance if a physical system has to spend some minimum time η before it switches to another mode. For example, in Section 4, the attitude control of a rigid body is studied. In real life, because of robustness issues and mechanical constraints, nozzles on a space launcher have indeed a minimum activation time.

Remark 11. Let t f denote the final time. If η is the minimal time between two switchings in [START_REF] You | Robust, near time-optimal control of nonlinear second order systems with model uncertainty[END_REF], then the total number of switchings N + 2l has an upper bound of t f /η .

The elaboration of a robust trajectory in ( 17) can be seen as an optimal control problem of switched-mode dynamical system. A recent survey on switched systems can be found in [START_REF] Zhu | Optimal control of hybrid switched systems: A brief survey[END_REF]. This theory deals with control systems where the dynamics can only take a finite number of modes. To determine the command law, one has to determine the switching times, as well as the different modes of the system. If the modes are fixed (in our case, it means that the components (i 1 , . . . , i N , j 1 , . . . , j l ) are fixed), it is often called a timing-optimization problem ; if not, a scheduling optimization problem. In [START_REF] Piccoli | Necessary conditions for hybrid optimization[END_REF][START_REF] Sussmann | Set-valued differentials and the hybrid maximum principle[END_REF], necessary conditions are derived, for trajectories of hybrid systems considering a fixed sequence of modes of finite length (in our setting, it corresponds to the Problem ( 18)). In [START_REF] Ali | Optimal control of switched dynamical systems under dwell time constraints[END_REF][START_REF] Wardi | Optimal control of switched-mode dynamical systems[END_REF], the authors develop numerical algorithms to solve both the timing and the scheduling problems. Their techniques rely heavily on gradient-like methods. However, the latter problem is much more complex because of its discrete nature: indeed the procedure needs to account for both continuous and discrete control variables, and can therefore be seen as a combinatorial optimization problem. Note that the paper [START_REF] Ali | Optimal control of switched dynamical systems under dwell time constraints[END_REF] deals with dwell time constraints. It consists in imposing a threshold η between two consecutive switching times which is the constraint we introduced at [START_REF] You | Robust, near time-optimal control of nonlinear second order systems with model uncertainty[END_REF]. Let us also mention other techniques to solve scheduling optimization problems, like zoning algorithms [START_REF] Shaikh | Optimality zone algorithms for hybrid systems computation and control: From exponential to linear complexity[END_REF], or relaxation methods, where discrete variables are temporarly relaxed into continuous variables [START_REF] Sorin | Optimal control of switching systems[END_REF].

Numerical results

In order to illustrate the results of Sections 2 and 3, we consider the problem of the attitude control of a rigid body. Let ω = (ω 1 , ω 2 , ω 3 ) be the angular velocity of the body with respect to a frame fixed on the body. Introducing the inertia matrix I, the Euler's equation for a rigid body, subjected to torques (b 1 , . . . , b m ), writes:

I ω = Iω ∧ ω + m k=1 b k .
In the case when the axes of the body frame are the axes of inertia of the body, the matrix I is diagonal: I = diag(I 1 , I 2 , I 3 ). The controlled Euler's equations can then be reduced to

ω(t) = f (ω(t), u(t)) ,
where for 1 k m, u k (t) ∈ {0, 1} almost everywhere, and the function f describing the dynamics writes:

f (ω 1 , ω 2 , ω 3 , u 1 , u 2 , u 3 , u 4 ) =    α 1 ω 2 ω 3 + m k=1 b k 1 u k α 2 ω 1 ω 3 + m k=1 b k 2 u k , α 3 ω 1 ω 2 + m k=1 b k 3 u k (20)
with α 1 = (I 2 -I 3 )/I 1 , α 2 = (I 3 -I 1 )/I 2 and α 3 = (I 2 -I 1 )/I 3 . This is with a slight abuse in the notations, because we still denote by b k the normalized vector (b k

1 /I 1 , b k 2 /I 2 , b k 3 /I 3
). The controllability of such a system has been studied in [START_REF] Bonnard | Mécanique céleste et contrôle des véhicules spatiaux[END_REF]. Let us mention here the papers [START_REF] Krstic | Inverse optimal stabilization of a rigid spacecraft[END_REF][START_REF] Outbib | Stabilizability of the angular velocity of a rigid body revisited[END_REF][START_REF] Windeknecht | Optimal stabilization of rigid body attitude[END_REF], that implement, in the special case of the stabilization of a rigid spacecraft, methods to stabilize the spacecraft towards the point (0, 0, 0), but once again, the controls used are not bang-bang. Note that ( 20) is a control-affine system, and therefore, Remark 2 applies.

In the following, we consider the numerical values

α 1 = 1, α 2 = -1, α 3 = 1, b 1 = [2, 1, 0.3], b 2 = [-2, -1, -0.3], b 3 = [0, 0, 1] and b 4 = [0, 0, -1],
and initial and final conditions x 0 = (0, 0, 0) and x f = (0.4, -0.3, 0.4).

We start by building an optimal trajectory for the L 1 cost t f 0 4 j=1 |u j (t)|dt + t f (the presence of t f ensures us not to obtain a trajectory with infinite final time). The resolution of such a problem with a L 1 cost can be numerically challenging. Numerical methods in optimal control are often categorized in two categories: direct methods and indirect methods. Whereas direct methods consist in a total discretization of the state and control spaces, indirect methods exploit Pontryagin maximum principle. (see [START_REF] Trélat | Optimal control and applications to aerospace: Some results and challenges[END_REF] for a survey on numerical methods in optimal control). The aim of the following subsection is to explain briefly the principle of a continuation method.

Computing the nominal trajectory

The nominal trajectory, optimal for the L 1 cost, is computed with a continuation procedure. The idea of such a procedure is to solve first an "easier" problem, and deform it step by step to solve the targeted problem. We introduce the continuation parameter λ ∈ [0, 1], and we consider the optimal control problem (P λ ) of steering the system (20) from x 0 to x f , by minimizing the cost

λ t f 0 4 i=1 |u j (t)| 2 dt + (1 -λ) t f 0 4 i=1 |u j (t)| dt + t f .
When λ = 0, we recognize our problem. For some λ ∈ [0, 1], solving problem (P λ ) is done by finding the zeros of a shooting function that results from the application of Pontryagin maximum principle. Solving a shooting problem is done with Newton like methods. Such methods are highly sensitive to their initialization, that can be very difficult, especially in the case of the minimization of the L 1 norm t f 0 |u(t)|dt. The continuation procedure is introduced to overcome this difficulty. For λ = 1, the cost is stricly convex in the controls, and writes

t f 0 4 i=1 |u j (t)| 2 dt + t f ,
for which the initialization of the induced shooting method is much easier. Therefore, we solve a sequence of optimal control problems, for values of λ decreasing from 1 to 0. The result of the shooting problem for some λ ∈]0, 1] serves as the initialization of another problem with λ < λ.

Robustifying the nominal trajectory

From this L 1 -minimal trajectory, represented on Figure 3, with three switching times that we denote (t 1 , t 2 , t 3 ) we build a new trajectory by solving the problem (17) with 3 needles (i.e., l = 3), λ 1 = λ 2 = 1, and taking η = 0.05 in Equation [START_REF] You | Robust, near time-optimal control of nonlinear second order systems with model uncertainty[END_REF]. As explained in Remark 5, we see that it is worthwile to have the additional switching times available as long as possible. That is, we force the additional switchings to occur after t 3 . Keeping in mind Equation [START_REF] You | Robust, near time-optimal control of nonlinear second order systems with model uncertainty[END_REF], this constraint can be written:

t i+1 -t i η (∀i ∈ 1, 3 ), s 1 -t 3 η, s i+1 -s i η (∀i ∈ 1, 6 ).
We find that the optimal triplet is (j 1 , j 2 , j 3 ) = (1, 4, 2), for which we have C = 0.77 and C r = 2.22.

We found this optimal triplet by exploring the 4 3 = 64 possibilities. We then used the heuristic that this solution would make a good choice to start looking for the solution with 4 needles (as it would have been to costly to examine the 4 4 = 256 possibilities). However we could not make the cost dicrease significantly (the best cost we found was C r = 2.07). This heuristic is very similar to what is used in Branch and Bound methods. Besides, as an element of comparison, the optimal couple when adding only two needles is (j 1 , j 2 ) = (1, 4), for which C r = 4.25, and the optimal solution when adding only on needle is j 1 = 2, for which C r = 30.28. Thus, we notice a substantial improvement when increasing the number of needles from 1 to 2 and from 2 to 3, whereas it seems less profitable to add a fourth one. We therefore stopped at 3 needles. The controls are displayed on Figure 3, and the components 1, 2 and 4, on which needles have been added, are represented in red.

In order to represent perturbations, we consider that the principal moments of inertia can vary, causing the coefficients α 1 , α 2 and α 3 to vary. Thus we consider the perturbed dynamics

f per (t, ω 1 , ω 2 , ω 3 , u 1 , u 2 , u 3 , u 4 ) =    α per,ε 1 (t)ω 2 ω 3 + m k=1 b k 1 u k α per,ε 2 (t)ω 1 ω 3 + m k=1 b k 2 u k , α per,ε 3 (t)ω 1 ω 2 + m k=1 b k 3 u k (21)
so that ε models the size of the perturbation. More precisely, we take α per,ε i (t) = α i + εh i (t), where h i (•) is some periodic function satisfying h i ∞ 1 (note that the exact expression of h i is not relevant here, as it is supposed to model any perturbation of the α i ). We denote by x per the solution of the Cauchy problem ẋ(t) = f per (t, x(t), u(t)), x(0) = x 0 .

We denote by x cor the corrected trajectory computed with our algorithm. We show, on Figure 5, the three trajectories, for ε = 0.78 and a cost C r = 2.22. We can see the perturbed trajectory x per drifting away from the reference trajectory x ref and away from the final point x f , whereas the corrected trajectory x cor eventually reaches a point very close to x f . Actually, for the trajectories represented on Figure 5, we have that x cor (t f )x f / x f = 5.5 × 10 -3 , whereas x per (t f )x f / x f = 1.3 × 10 -1 . Our algorithm has indeed been able to adjust the perturbed trajectory back towards x f .

One may wonder how this method behaves with respect to the choice of ε. As explained in Remark 6, we stop if two switching times are interchanged, that is, if δT is too big, as the initial vector of switching times satisfies a gap property [START_REF] You | Robust, near time-optimal control of nonlinear second order systems with model uncertainty[END_REF]. Actually, this is not strictly true, as we could have a "big" correction that does not change the ascending order of the switching times, for instance if we shift all the switching times in the same direction. However, we experimentally notice that the cost C r has an impact on the size of the perturbation we are able to absorb.

We build several trajectories, for which we apply our algorithm for increasing values of ε, until the algorithm fails as explained in Remark 6, for some ε max . We plot on Figure 4 ε max with respect to the cost C r (that is, for a given cost C r , ε max is the smallest value for which there is an interchanging of switching times). Even if the curve is not decreasing (for the reason explained above), we can see that having a low cost C r enables us to absorb bigger perturbations.

On Figure 6, we show the relative error x(t f )x f / x f for the perturbed x per and corrected x cor trajectories, for several values of ε. As we apply order one corrections, we see that our method shows better results for small values of ε, but also gives very satisfactory results for larger values of ε.

Conclusion

Starting with the expansion of the end-point mapping with respect to a needle like variation, we have shown in this paper how redundant switching times can be added in order to make a control more robust, for general control systems of the form ẋ(t) = f (t, x(t), u(t)). Those additional switching times can be seen as extra degrees of freedom meant to help us absorb perturbations. A potential application is to start from a bang-bang solution of an optimal control problem, that is usually not robust, and make it more robust. Then the gain in robustness compensates for the loss in optimality.

In the presence of a perturbation δx, the correction to apply to the switching times is the solution of an equation dE • δT = δx. It is natural to try to solve this equation while shifting the switching times as little as possible. The least-squares problem formulation is then the appropriate setting to find the solution of minimal (euclidian) norm of the previous equation, and it is given by δT = dE † • δx, for which we have the norm estimation δT 2 δx 2 /σ min . This enabled us to identify the measure for robustness: 1 σ min (t) 2 dt.

The numerical example studied in Section 4 is academic, and was used to legitimize the theoretical ideas explained previously. In a future work, we aim at applying the method to the complete (and more complex) attitude control system of a three-dimensional rigid body, for which we wish to control the angular velocity, as well as the orientation with respect to a fixed reference frame. To the three velocity variables will be added three angles to parametrize the orientation of the body. Thus, a challenge will come from the dimension of the state space (6), as well as the potentially bigger number of needle-like variations required to robustify a trajectory.

A Proof of proposition 1

In order to prove the differentiability of the end-point mapping, we start with the differentiability with respect to one component. The proof relies heavily on the expansion (3), that we recall first. )) if u i1 switches from a i1 to b i1 . f (t 1 , x(t 1 ), (. . . , b i1 , . . .))f (t 1 , x(t 1 ), u(t - 1 )) if u i1 switches from b i1 to a i1 .

Proof. Assume that at time t 1 the control u i1 switches from a i1 to b i1 , and that δt 1 > 0. Let us define the needle-like variation π = (t 1 , δt 1 , a i1 ) for the i 1 -th component of the control. Then, the control u π is represented by the vector (t 1 + δt 1 , . . . , t N , t f ) (figure 7): adding the needle-like variation π to the i 1 -th component, with value a i1 and length δt 1 is equivalent to shifting the opening time to t 1 + δt 1 . Thus, we have that u(t + 1 ) i1 = b i1 and u π (t + 1 ) i1 = a i1 . Hence, we obtain that, according to lemma 1

x π (t f ) = x(t f ) + δt 1 • v 1 (t f ) + o(δt 1 ), ( 22
)
where v 1 (•) is the solution of the Cauchy problem: v1 (t) = ∂f ∂x (t, x(t), u(t))v 1 (t), v 1 (t 1 ) = f (t 1 , x(t 1 ), u π (t + 1 ))f (t 1 , x(t 1 ), u(t + 1 )) = f (t 1 , x(t 1 ), (. . . , a i1 , . . .))f (t 1 , x(t 1 ), (. . . , b i1 , . . .)). (Between u π (t + 1 ) and u(t + 1 ), only the i 1 -th component differs.) If δt 1 < 0, define the variation π = (t 1 , δt 1 , 1) for the i 1 -th component of the control. Then again, the control u π is represented by the vector (t 1 + δt 1 , . . . , t N , t f ) (figure 7). Thus, we have that u(t - 1 ) = a i1 and u π (t - 1 ) i1 = 1. Thanks to lemma 1, we obtain that

x π (t f ) = x(t f ) -δt 1 • w 1 (t f ) + o(δt 1 ), (23)
where w 1 (•) is the solution of the Cauchy problem: ẇ1 (t) = ∂f ∂x (t, x(t), u(t))w 1 (t), w 1 (t 1 ) = f (t 1 , x(t 1 ), u π (t - 1 ))f (t 1 , x(t 1 ), u(t - 1 )) = f (t 1 , x(t 1 ), (. . . , b i1 , . . .))f (t 1 , x(t 1 ), (. . . , a i1 , . . .)) = -v 1 (t 1 ). Thus, by uniqueness we have w 1 = -v 1 , and from ( 22) and ( 23), we obtain:

x π (t f ) = x(t f ) + δt 1 • v 1 (t f ) + o(δt 1 ).
We can proceed the exact same way if at t 1 , the control u i1 switches from b i1 to a i1

The general result at proposition 1 follows by an immediate iteration.

Figure 1 :

 1 Figure 1: Changing the switching times induces a displacement at the final time.

  (a) Initial control. (b) Adding four switching times i.e., 2 needles.

Figure 2 :

 2 Figure 2: Principle of adding needles.

4 ( 4 (

 44 a) Controls for the minimal L 1 trajectory. C = 0b) Controls with three needles. C = 0.77, Cr = 2.22

Figure 3 :

 3 Figure 3: Improving the robustness of a trajectory adding needles.

Figure 4 :

 4 Figure 4: Size of the maximal perturbation absorbed with respect to the robustness of a trajectory

Lemma 1 .

 1 Let t 1 ∈ [0, t f [, and let u π1 (•) be a needle-like variation of u(•), with π 1 = (t 1 , δt 1 , u 1 ). Then x π1 (t f ) = x(t f ) + |δt 1 | v π1 (t f ) + o(δt 1 )

Figure 5 :

 5 Figure 5: Reference, perturbed and corrected trajectories for ε = 0.78, C r = 2.22.

  δt 1 < 0.

Figure 7 :

 7 Figure 7: Shifting an opening time is equivalent to add a needle.

  N ) min (t1,...,t N ) C(t 1 , . . . , t N ). s.t. E(t 1 , . . . , t N ) = x f

[START_REF] Lee | Foundations of optimal control theory[END_REF]

, i 1 , . . . , t N , i N , t f ) is an optimal solution of the optimization problem min (i1,...,i

  Figure 6: Tracking results for several values of ε.where v π1 (•) is the solution of a Cauchy problem on [t 1 , t f ] Proposition 3. We denote by u the control (t 1 , . . . , t N , t f ) and x(•) the associated trajectory of the control system. Let δt 1 ∈ R be small enough. ThenE(t 1 + δt 1 , t 2 , . . . , t N , t f ) = E(t 1 , . . . , t N , t f ) + δt 1 • v 1 (t f ) + o(δt 1 ),where v 1 (•) is the solution of the Cauchy problem on [t 1 , t f ]:

	10 -1
	f
	f / x
	-x
	f )
	x(t
				x per x cor
				ε
			vπ1 (t) =	∂f ∂x	(t, x(t), u(t))v π1 (t),
	v1 (t) =	∂f ∂x	(t, x(t), u(t))v

v π1 (t ) = f (t 1 , x(t 1 ), u 1 )f (t 1 , x(t 1 ), 1 )). i (t), v 1 (t 1 ) = f (t 1 , x(t 1

), (. . . , a i1 , . . .))f (t 1 , x(t 1 ), u(t + 1

Acknowledgements

This study has been performed in the frame of the CNES Launchers Research & Technology program.


Guillaume Chèze
email: guillaume.cheze@math.univ-toulouse.fr

Luca Amodei
email: luca.amodei@math.univ-toulouse.fr

Guillaume Ch

How to cut a cake with a gram matrix

In this article we study the problem of fair division. In particular we study a notion introduced by J. Barbanel that generalizes super envy-free fair division. We give a new proof of his result. Our approach allows us to give an explicit bound for this kind of fair division. Furthermore, we also give a theoretical answer to an open problem posed by Barbanel in 1996. Roughly speaking, this question is: how can we decide if there exists a fair division satisfying some inequalities constraints? Furthermore, when all the measures are given with piecewise constant density functions then we show how to construct effectively such a fair division.

Introduction

In the following X will be a measurable set. This set represents an heterogeneous good, e.g. a cake, that we want to divide between n players. A division of the cake is a partition X = ⊔ n i=1 X i , where each X i is a measurable subset of X. After the division X i is given to the i-th player. A natural and old problem is: how to get a fair division? This problem appears when we study division of land, time or another divisible resource between different agents with different points of view. These problems appear in the economics, mathematics, political science, artificial intelligence and computer science literature, see [Mou03, BCE + 16]. In order to study this problem, to each player is associated a non-atomic probability measure µ i . Thus, in particular µ i (X) = 1, and µ i (A ⊔ B) = µ i (A) + µ i (B), where A, and B are disjoint measurable sets. These measures represent the preference of each player. Severall notions of fair divisions exist:

• Proportional division: ∀i, µ i (X i ) ≥ 1/n.

• Exact division: ∀i, ∀j, µ i (X j ) = 1/n.

• Equitable division: ∀i, ∀j, µ i (X i ) = µ j (X j ).

• Envy-free division: ∀i, ∀j, µ i (X i ) ≥ µ i (X j ). All these fair divisions are possible, see e.g [Ste48, DS61, BT96, RW98, Chè17, SHS, Wel85]. Some fair divisions are possible under some conditions. Barbanel has shown in [START_REF] Barbanel | Super envy-free cake division and independence of measures[END_REF] that a super envy-free division is possible if and only if the measures µ i are linearly independent. We recall the definition of a super envy-free fair division:

• Super envy-free division: ∀i, ∀j = i, µ i (X i ) > 1/n > µ i (X j ). Actually, if the measures are linearly independent then there exists a real δ > 0 such that: µ i (X i ) ≥ 1/n + δ, and µ i (X j ) ≤ 1/n -δ/(n -1).

Date: July 6, 2017.

We can define an even more demanding fair division. For example, we can imagine that the first player would like to get a partition such that: µ 1 (X 1 ) = 1/n + 3δ, µ 1 (X 3 ) = 1/n + 2δ, µ 1 (X 4 ) = µ 1 (X 5 ) = 1/n + δ, µ 1 (X 2 ) = µ 1 (X 6 ) = 1/n -6δ.

This means that the third, the forth and the fifth player are friends with the first player, but the second and sixth player are not friends with this player. Furthermore, the first player prefers the third to the forth and the fifth player. We can also imagine that the other players have also preferences between the other players. These kinds of conditions have also been studied by Barbanel in [Bar96a,[START_REF] Barbanel | The geometry of efficient fair division[END_REF]. This leads to a notion of fair division that we call hyper envy-free.

Definition 1. Consider a matrix K = (k ij ) ∈ M n (R), such that for all i = 1, . . . , n, n j=1 k ij = 0 and a point p = (p 1 , . . . , p n ) such that n j=1 p j = 1, with p i ≥ 0. We say that a partition X = ⊔ n i=1 X i is hyper envy-free relatively to K and p when there exists a real number δ > 0 such that

µ i (X j ) = p j + k ij δ.
For example, for a super envy-free division we have

k ij = -1/(n -1), k ii = 1 and p = (1/n, . . . , 1/n).
Barbanel has given a criterion for the existence of an hyper envy-free division. Unfortunately, the proof of this result does not give an explicit bound on δ and a natural question is: How big δ can be?

In this paper we give a new proof of Barbanel's result. Furthermore, our strategy generalizes the approach for computing super envy-free fair division given by Webb in [START_REF] Webb | An algorithm for super envy-free cake division[END_REF]. This allows us to give a bound on δ in terms of the measures.

Another question related to fair division and asked by Barbanel is the following, see [START_REF] Barbanel | On the possibilities for partitioning a cake[END_REF]:

Suppose that p = (p 1 , . . . , p n ) is a point such that p 1 + • • • + p n = 1 with p i positive and r ij are n 2 relations in {<, =, >}. How can we decide if there exists a partition X = ⊔ n i=1 X i such that µ i (X j ) r ij p j ?
This problem is also a generalization of the super envy-free fair division problem. Indeed, super envy-free fair division corresponds to the situation where p is (1/n, . . . , 1/n), r ii is the relation ">" and r ij is "<". We give a theoretical answer to Barbanel's question in the last section of this article.

The organization of our article is the following. In the next section we present our toolbox. We recall the Dvoretzky, Wald, Wolfowitz's theorem that will be the main ingredient of our proof. In Section 2, we prove the existence of an hyper envy-free division under some linear conditions on the measures µ i . This gives a new proof of Barbanel's theorem. A direct consequence of our construction gives a bound on δ. At last, in Section 3 we solve the open question asked by Barbanel. Furthermore, when all the measures are given by piecewise constant density functions we show how to construct effectively an hyper envy-free fair division. This gives a method to construct a partition such that µ i (X j ) r ij p j . 1. Our toolbox Definition 2. A matrix M = (µ i (X j )) is said to be a sharing matrix when X = ⊔ n i=1 X i is a partition of X. A sharing matrix is a row stochastic matrix, this means that each coefficients are nonnegative and the sum of the coefficients of each row is equal to 1. In the following, we will use classical results about this kind of matrices. We recall below without proofs these results.

Lemma 3. We denote by e the vector (1, . . . , 1) T . If S = (s ij ) and T = (t ij ) are matrices, and s, t ∈ R are such that S e = s e and T e = t e, then:

(1) N = S T is a matrix such that N e = st e.

(2) If S is invertible, then s = 0 and S -1 is such that S -1 e = 1 s e. In the following we define a new measure. This measure will be usefull to write each µ i (X j ) in term of the same measure.

Definition 4. We denote by µ the measure

µ = µ 1 + • • • + µ n . The Radon-Nikodym derivative dµ i /dµ is denoted by f i .
The Radon-Nikodym derivative f i exists because µ i is absolutely continuous relatively to µ. By definition we have for any measurable subset A ⊂ X:

µ i (A) = A f i (x)dµ.
Definition 5. The Dvoretzky, Wald, Wolfowitz set (DWW set) is the set of all matrices X η j (x)dµ i where η 1 , . . . , η n , are positive functions such that η

1 (x) + • • • + η n (x) = 1.
Thanks to the previous notations a matrix in the DWW set can be written

X η j (x)f i (x)dµ .
The following deep result will be one of the main ingredient of our proof.

Theorem 6 (Dvoretzky, Wald, Wolfowitz [START_REF] Dvoretzky | Relations among certain ranges of vector measures[END_REF]). The DWW set is the set of sharing matrices. This theorem is classical when we study fair division problems. As a corollary we get that the set of sharing matrices is a convex set. Barbanel's proof uses only this corollary.

In his book Barbanel gives the following lemma, see [Bar05, Lemma 9.10].

Lemma 7. The functions f i satisfy:

f 1 (x) + • • • + f n (x) = 1, 0 ≤ f i (x) ≤ 1, ∀x ∈ X.
We can then construct a sharing matrix thanks to the density functions f i . The following Gram matrix will be the basis of our construction.

Corollary 8. The Gram matrix

G(f 1 , . . . , f n ) = X f j (x)f i (x)dµ , is a sharing matrix.
We remark that the matrix G(f 1 , . . . , f n ) is well defined since by Lemma 7,

f i ∈ L ∞ (X).
Thanks to the Dvoretsky, Wald, Wolfowitz's theorem we know that G(f 1 , . . . , f n ) is a sharing matrix. Thus, there exists a partition X = ⊔ n i=1 X i such that

X f j (x)f i (x)dµ = µ i (X j ).
Remark 9. The matrix G(f 1 , . . . , f n ) gives a proportional fair division. Indeed, we apply Cauchy-Schwartz's inequality in L 2 (X) where

< f, g >= X f (x)g(x)dµ,
and we get < 1, f i >≤ f i 1 . However

1 2 = X dµ = X dµ 1 + • • • + dµ n = n and < 1, f i >= X f i dµ = µ i (X) = 1.
This gives

f i 2 ≥ 1/n. Thus all diagonal elements of G(f 1 , . . . , f n ) satisfy f i 2 = µ i (X i ) ≥ 1/n.
This gives a proportional fair division.

This matrix has also other usefull properties. For example, this matrix is symmetric. Now, we recall a classical property of a Gram matrix.

Lemma 10. Let x denote the vector x = (x 1 , . . . , x n ) T , and ker(G) the kernel of a matrix G. We have the equivalence

x ∈ ker G(f 1 , . . . , f n ) ⇐⇒ n i=1 x i f i = 0 in L 2 (X) ⇐⇒ n i=1 x i µ i = 0.
Proof. We have:

G(f 1 , . . . , f n ) x = < f 1 , n j=1 x j f j >, . . . , < f n , n j=1 x j f j > T . Thus x ∈ ker G(f 1 , . . . , f n ) means < f i , n j=1
x j f j >= 0, ∀i = 1, . . . , n. Combining these equalities we get n j=1 x j f j 2 = 0 which implies n j=1 x j f j = 0 in L 2 (X). The converse implication is straightforward. Furthermore, we have the following equivalence:

n i=1 x i f i = 0 in L 2 (X) ⇐⇒ for all measurable set A, A n i=1 x i f i (x)dµ = 0 ⇐⇒ for all measurable set A, n i=1 x i A f i (x)dµ = 0 ⇐⇒ for all measurable set A, n i=1 x i µ i (A) = 0.
This gives the desired result.

In terms of the measures µ i this property implies that G(f 1 , . . . , f n ) is nonsingular if and only if the measures µ i are independent.

Lemma 11. We denote by G(f 1 , . . . , f n ) + the Moore-Penrose pseudo-inverse of G(f 1 , . . . , f n ) and by e the vector e = (1, . . . , 1) T . We have

G(f 1 , . . . , f n ) + e = e.
Proof. For a matrix G we have the following classical property, see e.g. [BIG03, Ex. 25, p.115]:

(⋆)

G + = lim ρ→0 ρ>0 G T (GG T + ρI) -1 .
As G(f 1 , . . . , f n ) is symmetric and row stochastic, then G(f 1 , . . . , f n ) T is also row stochastic. Thus, by Lemma 3, the product

G(f 1 , . . . , f n ) G(f 1 , . . . , f n ) T is a row stochastic matrix. Then, G(f 1 , . . . , f n ) G(f 1 , . . . , f n ) T + ρ I e = (1 + ρ) e.
Since ρ > 0, the matrix

G(f 1 , . . . , f n ) G(f 1 , . . . , f n ) T + ρ I is nonsingular. Indeed, G(f 1 , . . . , f n ) G(f 1 , . . . , f n )
T is a semi-definite positive matrix and thus the matrix

G(f 1 , . . . , f n ) G(f 1 , . . . , f n ) T + ρ I is definite positive. Using again Lemma 3, we deduce that G(f 1 , . . . , f n ) G(f 1 , . . . , f n ) T + ρ I -1 e = 1 1+ρ e. Since G(f 1 , . . . , f n ) T is row stochastic, then property (⋆) gives the desired conclusion.
The next lemma shows how to get new sharing matrices. This idea was already present in Webb's algorithm for computing a super envy-free fair division, see [START_REF] Webb | An algorithm for super envy-free cake division[END_REF].

Lemma 12. Let M ∈ M n (R) be a sharing matrix and S ∈ M n (R) be a row stochastic matrix. Then M S is a sharing matrix.

Proof. Suppose that the sharing matrix is given by M = µ i (X j ) and denote by S = (s ij ) the stochastic matrix. We set η j = n k=1 s kj χ X k , where χ X k is the indicator function associated to X k . We remark that η j (x) is positive and

η 1 (x) + • • • + η n (x) = 1 because X = ⊔ n i=1 X i
is a partition and S is a row stochastic matrix. Furthermore, we have

X η j (x)f i (x)dµ = n k=1 s kj X χ X k (x)f i (x)dµ = n k=1 µ i (X k )s kj .
This gives:

η j (x)f i (x)dµ ij = M S.
Then by the Dvoretzky, Wald, Wolfowitz's theorem M S is a sharing matrix.

Existence of hyper envy-free division and an explicit bound

2.1. Sufficient conditions. Definition 13. A matrix K = (k ij ) is said to be proper relatively to the measures µ 1 , . . . , µ n when the two following conditions hold:

(1) for all i = 1, . . . , n,

n j=0 k ij = 0, (2) if n i=0 λ i µ i = 0, where λ i ∈ R, then n i=0 λ i k ij = 0, for all j = 1, . . . , n.
With this definition we can state our extension of Barbanel's theorem.

Theorem 14. If we denote by P the matrix where each row is (p 1 , . . . , p n ), and n i=1 p i = 1, with p i > 0, and K is a proper matrix, then for

0 ≤ δ ≤ B := min i (p i ) max ij | G(f 1 , . . . , f n ) + K ij |
we have:

(1) G(f 1 , . . . , f n ) + (P + δK) is a row stochastic matrix.

(2) G(f 1 , . . . , f n ) G(f 1 , . . . , f n ) + (P + δK) = (P + δK).

(3) P +δK is a sharing matrix, i.e. there exists an hyper envy-free fair division relatively to K and p, with δ ≤ B.

The last item is an effective statement of Barbanel's theorem. Our approach allows us to give a bound B on δ. The case δ = 0 is classical and follows form the convexity of the set of sharing matrices.

Proof. By Lemma 11 we have G(f 1 , . . . , f n ) + e = e. Using Lemma 3 and the properties of P and K we get

G(f 1 , . . . , f n ) + (P + δK) e = e.
In order to prove the first item we now have to show that the coefficients of G(f 1 , . . . , f n ) + (P + δK) are non-negative for δ sufficiently small. We set

G(f 1 , . . . , f n ) + = (g + ij ) and K = (k ij ). The ij-th coefficient of G(f 1 , . . . , f n ) + (P + δK) is of the following form: n l=1 g + il p j + δ n l=1 g + il k lj = p j + δ n l=1 g + il k lj ,
since G(f 1 , . . . , f n ) + is row stochastic. Now, we denote by c ij the coefficients of G(f 1 , . . . , f n ) + K. We have:

p j + δc ij ≥ min(p i ) + δc ij .
Thus, if δ satisfies min(p i ) + δc ij ≥ 0, we obtain the desired result. This condition is trivially satisfied when c ij ≥ 0.

In the general case, we have: min

(p i ) + δc ij ≥ 0 ⇐⇒ min(p i ) |c ij | ≥ δ.
Thus, if δ satisfies the bound given in the theorem, then all the coefficients of G(f 1 , . . . , f n ) + (P + δK) are non-negative.

In order to simplify the notation, we denote by G the Gram matrix G(f 1 , . . . , f n ). In order to show the second item of the theorem, we use the well known equality satisfied by the pseudo-inverse G + : GG + = Π Im(G) , where Π Im(G) is the orthogonal projection onto the image of G (that we denoted by Im(G)). We are going to show the inclusion Im(P + δK) ⊂ Im(G) and this will imply the desired equality.

Let x = (x 1 , . . . , x n ) T a vector in ker(G). From Lemma 10 we have n i=1 x i µ i = 0 which implies x T K = 0 since K is a proper matrix. The matrix G is symmetric, therefore x T G = 0 and x T Ge = x T e = 0 which implies x T P = 0. This gives x T (P + δK) = 0 and we conclude that x ∈ Im(P + δK) ⊥ (the orthogonal complement of Im(P + δK)). We have obtained the inclusion ker(G) ⊂ Im(P + δK) ⊥ . Thus ker(G) ⊥ ⊃ Im(P + δK). Furthermore, we have ker(G) ⊥ = Im(G T ) = Im(G), since G is symmetric. This gives Im(G) ⊃ Im(P + δK), and thus the second item is proved.

The third item is a direct consequence of the previous items and Lemma 12.

The previous theorem shows that if the measures are linearly independent then when can get an hyper envy-free fair division for all matrix K = (k ij ) such that n j=1 k ij = 0, see Definition 13. Intuitively, if the measures are "nearly" linearly dependent then δ will be small. The following gives a precise statement for this intuition.

As we suppose that the measures are linearly independent, G = G(f 1 , . . . , f n ) is inversible, see Lemma 10. We denote by gij the coefficients of G -1 . Cramer's rule gives the following equality

gij = det(G ji ) det(G) ,
where G ji is the matrix where we have substituted in the matrix G the i-th column by the column vector e j = (0, . . . , 0, 1, 0, . . . , 0) T . Thus, if gij = 0 then we have

0 = det(G) - 1 gij det(G ji ),
and by linearity of the determinant relatively to the i-th column we deduce:

0 = det G - 1 gij E ji ,
where E ji is the matrix with just one 1 in the i-th column and j-th row and 0 elsewhere. We deduce that G -1/g ij E ij is a singular matrix. We denote by Σ n the set of all n by n singular matrices.

It follows 1/g ij E ji ≥ d(G, Σ n ), where . is a norm in the space of square matrices and d(G, Σ n ) is the distance relatively to this norm between the matrix G and the set Σ n . We suppose that E ji = 1. This assumption is not restrictive since all the classical norms satisfy this property. Thus 1

|g ij | ≥ d(G, Σ n ).
We deduce that the coefficients c ij of G(f 1 , . . . , f n ) -1 K satisfy:

|c ij | ≤ n max ij |k ij | max ij gij ≤ n max ij |k ij | d(G, Σ n ) . Thus min i (p i ) |c ij | ≥ min i (p i ) d(G, Σ n ) n max ij |k ij | .
This gives the following corollary:

Corollary 15. Let . be a norm in the space of n×n matrices such that E ij = 1 where E ij is the matrix with one 1 in the i-th row and j-th column and 0 elsewhere. We denote by d the associated distance. Let Σ n be the set of singular n × n matrices.

If the measures are linearly independent, then for all p and all K = (k ij ) such that n j=1 k ij = 0, for all i = 1, . . . , n, there exists an hyper envy-free fair division with δ satisfying the inequality:

δ ≤ min i (p i ) d(G, Σ n ) n max ij |k ij | .
This bound is smaller than the previous one but it gives a natural relation between the size of δ, p, K, n and the linear independency of the measures. The coefficient d(G, Σ n ) describes the linear independency of the measures. Indeed, we use a Gram matrix because the density functions belong to an infinite dimensional vector space and thus we cannot express the determinant of these functions in a finite basis of the ambient space. The Gram matrix encodes all the linear relations between the measures. It is natural to control how the measures are "linearly independent" in term of the distance between this matrix and the set of singular matrices Σ n .

Our corollary can be used for all classical norms in the space of n × n matrices. The choice of a norm and then of a metric corresponds to different points of view. As an example, if we consider the norm . ∞ we could measure the quality of the corresponding fair division in term of the distance d ∞ (M, I) = M -I ∞ , where I is the identity matrix. Indeed, I corresponds to the case where each player thinks that he or she has the maximum of his or her utility function.

If M = µ i (X j ) is a sharing matrix, then d ∞ (M, I) = max i n j=1 |µ i (X j ) -δ ij |.
Then 1 -d ∞ (M, I) can be seen as a kind of Rawlsian social welfare function. Thus the choice of a metric is related to a social welfare function.

Necessary conditions.

In order to make more natural the previous conditions on the matrix K, we give the following necessary condition: Lemma 16. Suppose that P + δK is a sharing matrix, where P is the matrix where each row is (1/n, . . . , 1/n) and δ > 0. Then K = (k ij ) is a proper matrix.

Proof. As P is a stochastic matrix, and P + δK is a sharing matrix, we deduce that n j=1 k ij = 0. Now suppose that we have the following relation between the measure µ i : n i=0 λ i µ i = 0. This gives the following relation:

(♯) n i=0 λ i (1/n + δk ij ) = 0,
because there exists measurable sets X j such that µ i (X j ) = 1/n+δk ij since P +δK is a sharing matrix. Furthermore, as n i=0 λ i µ i = 0, we get that n i=1 λ i µ i (X) = 0. Thus n i=1 λ i = 0. As δ > 0, the equation (♯) gives the desired result:

i=0 λ i k ij = 0.

Effective methods

In the following, we are going to study how to test if a fair division satisfying some inequalities constraints exists. Then we will show how to get this kind of fair division when the density functions are piecewise constant.

How to test the existence of a fair division given by a relation matrix.

In this section, we use the previous theorem in order to study another fair division problem. Suppose that p = (p 1 , . . . , p n ) is a point such that p 1 + • • • + p n = 1 with p i positive and r ij are n 2 relations in {<, =, >}. Barbanel has asked the question: How can we decide if there exists a partition X = ⊔ n i=1 X i such that µ i (X j ) r ij p j ?

If there exists a proper matrix K = (k ij ) such that k ij r ij 0 then by Theorem 14 there exists a partition X = ⊔ n i=1 X i such that µ i (X j ) r ij p j . Conversely, if X = ⊔ n i=1 X i is a partition such that µ i (X j ) r ij p j then the matrix K = (k ij ) where k ij = µ i (X j ) -p j is a proper matrix such that k ij r ij 0. Thus as already remarked by Barbanel, the previous question can be restated as: How can we decide if there exists a proper matrix K satisfying k ij r ij 0? Suppose that we know a basis {λ 1 , . . . , λ l } for the relations between the measures µ i , where λ α = (λ α1 , . . . , λ αn ). Then, we just have to solve the following system:

(S)      ∀i, n j=1 k ij = 0, ∀j, ∀α, n i=1 λ αi k ij = 0, ∀i, ∀j, k ij r ij 0.
We thus have to solve a linear system of equalities and inequalities. This can be done for example by using the Fourier-Motzkin elimination method, see e.g. [START_REF] Kuhn | Solvability and consistency for linear equations and inequalities[END_REF]. Now the question is : How can we find a basis {λ 1 , . . . , λ l } of relations between the µ i ?

If the measure µ i are given thanks to a density function ϕ i (x) relatively to a given measure, e.g. the Lebesgue measure, then we can answer the previous question. Indeed, we can set

f i (x) = ϕ i (x) n i=1 ϕ i (x)
,

for x such that n i=1 ϕ i (x) = 0 and f i (x) = 0 otherwise. Indeed, the density of µ = µ 1 + • • • + µ n is ϕ 1 (x) + • • • + ϕ n (x).
Thus, if ϕ i are given explicitly then we can deduce an explicit formula for the f i . This gives

X f i f j dµ = X ϕ i (x) ϕ j (x) n i=1 ϕ i (x) dx.
If these expressions are computable then we have an explicit form for G(f 1 , . . . , f n ). Thus in this situation the basis {λ 1 , . . . , λ l } can be computed easily. Indeed, this basis is a basis of ker G(f 1 , . . . , f n ) thanks to Lemma 10. In conclusion, the system S can be solved with linear algebra only.

Remark 17. As mentioned before, the previous strategy gives a theoretical method for solving Barbanel's problem. Our solution do not give an algorithm in the Robertson-Webb model, see [START_REF] Robertson | Cake-cutting algorithms -be fair if you can[END_REF], [BCE + 16, Chapter 13]. However, if X = [0, 1] and the ϕ i are polynomials or simple functions i.e. linear combination of indicator functions, then G(f 1 , . . . , f n ) can be constructed effectively and our method can be used.

Example 18. Suppose that X = [0, 1] and that the three measures are given with their densities ϕ 1 (x) = 10 χ [0,1/10] (x), ϕ 2 (x) = 10 9 χ [1/10,1] (x), ϕ 3 (x) = χ [0,1] (x). Suppose also that p = (1/3, 1/3, 1/3) and that R 1 = (r

(1) ij ) R 1 = (r (2) ij ) are the following matrix R 1 =   > = < > > < < < >   , R 2 =   > = < < > > < > >   .
As we can compute easily the expression

X ϕ i (x)ϕ j (x) n i=1 ϕ i (x)
dx. The Gram matrix associated to this situation is:

G(f 1 , f 2 , f 3 ) =          
          .
We can remark that this matrix is a sharing matrix corresponding to a proportinal fair division as proved in Remark 9, but it is not an envy-free fair division.

A basis for the kernel of this matrix is given by {(1, 9, -10)}. Thus the only relation between the measure is: µ 1 + 9µ 2 = 10µ 3 . Now, we have to solve the following system when s = 1 or s = 2:

(S s )      ∀i, 3 j=1 k ij = 0, ∀j, 3 i=1 λ i k ij = 0, ∀i, ∀j, k ij r (s) ij 0,
where (λ 1 , λ 2 , λ 3 ) = (1; 9; -10). When s = 1 the system have no solution. Indeed, we must have k 11 > 0, k 21 > 0, k 11 + 9k 21 = 10k 31 , and k 31 < 0 which gives a contradiction. When s = 2 we can compute solutions. For example the following matrix K = (k ij ) is a solution of (S 2 ):

K =         1 0 -1 - 1 3 1 9 2 9 - 2 10 1 10 1 10         .
3.2. The case of piecewise constant density functions. In this last section, we consider the case where the density functions are piecewise constant. This means ϕ i (x) = j c ij χ Iij (x), where χ Iij is the indicator function of the interval

I ij = [x ij , x i(j+1) ], x ij ∈ [0, 1].
Suppose that the measures µ i are given by these density functions and that we want to construct an hyper envy-free fair division relatively to a given matrix K = (k ij ) and p = (p 1 , . . . , p n ). We denote by I the set of intervals [x kl , x mn ] such that there exists no x ij satisfying x kl < x ij < x mn . This means that on each interval I ∈ I the density functions are constant.

Now, we explain how to construct an hyper envy-free fair division: Since all the densities ϕ i are constant on I we have for all subintervals J of I:

µ i (J)ℓ(I) = ℓ(J)µ i (I),
where ℓ(I) is the length of I for the Lebesgue measure.

Thus if µ i (I) = 0 we have:

µ i (J) µ i (I) = ℓ(J) ℓ(I) .
Now, consider α 1,I , . . . , α n,I such that α i,I ≥ 0 and i α i,I = 1. We can divide each I ∈ I in order to have a partition of I into n subintervals I 1 , . . . , I n satisfying ℓ(I j ) = α j,I ℓ(I). The previous relation implies (♭) µ i (I j ) = α j,I µ i (I) for all i, j.

Thus we set X j = ⊔ I∈I I j , and we get:

X = ⊔ n j=1 X j and µ i (X j ) = I∈I α j,I µ i (I).
Now, in order to find an hyper envy-free fair division relatively to K and p we have to solve:

(S ′ )          ∀I ∈ I, n i=1 α i,I = 1, ∀i ∈ {1, . . . , n}, ∀I ∈ I, α i,I ≥ 0, ∀i, j ∈ {1, . . . , n}, I∈I α j,I µ i (I) = p j + k ij δ, δ > 0.
Example 19. With the density functions and the matrix K = (k ij ) of Example 18, and p = (1/3, 1/3, 1/3), we have I = {[0, 1/10], [1/10, 1]} and the system (S ′ ) gives:

                   3 i=1 α i,[0,1/10] = 1 ∀i ∈ {1, . . . , 3}, α i,[0,1/10] ≥ 0, 3 i=1 α i,[1/10,1] = 1 ∀i ∈ {1, . . . , 3}, α i,[1/10,1] ≥ 0, ∀i, j ∈ {1, . . . , 3}, I∈I α j,I µ i (I) = 1/3 + k ij δ, δ > 0.
We get the following solution: We set a := α 1,[0,1/10] and a must satisfies the condition a ∈]1/3, 2/3[. Then α 2,[0,1/10] = 1/3, α 3,[0,1/10] = 2/3 -a, α 1,[1/10,1] = 4/9 -1/3a, α 2,[1/10,1] = 8/27+1/9a, α 3,[1/10,1] = 7/27+2/9a, δ = a -1/3. Thus, with α = 1/2 we get:

a := α 1,[0,1/10] = 1/2, α 2,[0,1/10] = 1/3, α 3,[0,1/10] = 1/6, α 1,[1/10,1] = 5/18, α 2,[1/10,1] = 19/54, α 3,[1/10,1] = 10/27, δ = 1/6.
We divide [0, 1/10] and [1/10, 1] in order to respect the condition (♭). For the interval I = [0, 1/10] we get: I 1 = [0, 1/20], I 2 = [1/20, 1/20 + 1/3 × 1/10] = [1/20, 1/12] and I 3 = [1/12, 1/10]. With the same method we divide [1/10, 1] and then we get:

X 1 = 0, 1 20 ⊔ 1 10 , 7 20 , X 2 = 1 20 , 1 12 ⊔ 7 20 , 2 3 , X 3 = 1 12 , 1 10 ⊔ 2 3 , 1 .
By construction, this partition gives an hyper envy-free fair division relatively to K and p = (1/3, 1/3, 1/3) with δ = 1/6. Indeed, we can check that the sharing matrix associated to the previous partition is Therefore, this partition gives also a fair division respecting the inequalities given in the matrix R 2 .

          1
In conclusion, when the density functions are piecewise constant then we have a method to compute an hyper envy-free fair division. Thus we can : decide if a fair division given by a relation matrix exists, and then construct a partition giving this fair division.

Remark 20. The pseudo-inverse of the matrix G(f 1 , f 2 , f 3 ) obtained in Example 18 and Example 19 is Thus Theorem 14 gives the existence of an hyper envy-free fair division relatively to K and p with δ < 455/1536 ≈ 0.296 without computing a solution of (S ′ ). In the above example, the solution of (S ′ ) gives δ = a -1/3 with a < 2/3. Thus in this particular situation an hyper envy-free fair division exists with δ < 1/3. This shows that our bound B is not optimal in this case, but the order of magnitude of B is not too small.

G(f 1 , f 2 , f 3 ) + =      

Conclusion

The key ingredient of this article is the Gram matrix G(f 1 , . . . , f n ). Thanks to this matrix we can construct the sharing matrix P + δK, see Theorem 14, and we can also compute the relations between the measures, see Lemma 10. This matrix seems very useful and we have already noticed that this matrix corresponds to a proportional and symmetric (i.e. µ i (X j ) = µ j (X i )) fair division, see Remark 9. Thus a natural question appears: Is G(f 1 , . . . , f n ) the optimal fair division for a certain criterion and what is the meaning of this criterion in terms of fair division?

.

  The matrixG(f 1 , f 2 , f 3 ) + K is G(f 1 , f 2 , f 3 ) + K = f 1 , . . . , f n ) + K ij | = 512/455.As mini (p i ) = 1/3, we get B = 455/1536.

Acknowledgments:

The authors thank L.-M. Pardo for fruitful discussions during the preparation of this article.