G Carlier
email: carlier@ceremade.dauphine.fr

A Lachapelle
email: lachapelle@ceremade.dauphine.fr

A numerical approach for a class of risk-sharing problems

Keywords: risk-sharing, comonotonicity, sup-convolution, calculus of variations, numerical approximation

This paper deals with risk-sharing problems between many agents, each of whom having a strictly concave law invariant utility. In the special case where every agent's utility is given by a concave integral functional of the quantile of her individual endowment, we fully characterize the optimal risk-sharing rules. When there are many agents, these rules cannot be computed analytically. We therefore give a simple convergent algorithm and illustrate it on several examples.

Introduction

Risk-sharing problems have their roots in the seminal works of Arrow [START_REF] Arrow | Uncertainty and the welfare of medical care[END_REF], [START_REF] Arrow | Optimal insurance and generalized deductibles[END_REF] and Borch [START_REF] Borch | Equilibrium in a reinsurance market[END_REF] in insurance and have received a lot of attention since. Starting from the case of two-agents having preferences given by expected utilities, the theory has been developed in recent years in particular to incorporate more general law invariant preferences such as rank-dependent utilities or monetary risk measures in the financial literature (see Dana [START_REF] Dana | A representation result for concave Schur concave functions[END_REF], Jouini, Schachermayer and Touzi, [START_REF] Jouini | Optimal risk sharing for law invariant monetary utility functions[END_REF], Carlier and Dana [START_REF] Carlier | Law invariant concave utility functions and optimization problems with monotonicity and comonotonicity constraints[END_REF], [START_REF] Carlier | Two-Persons Efficient Risk-Sharing and Equilibria for Concave Law-Invariant Utilities[END_REF], [START_REF] Carlier | Optimal demand for contingent claims when agents have law invariant utilities[END_REF] as well as the book by Föllmer and Schied [START_REF] Föllmer | Stochastic finance. An introduction in discrete time[END_REF] and the references therein).

The fact that comonotonicity is a key property in efficient risk-sharing problems is well-known and has been particularly emphasized by Landsberger and Meilijson in [START_REF] Landsberger | Comonotone allocations, Bickel Lehmann dispersion and the Arrow-Pratt measure of risk aversion[END_REF] (see [START_REF] Carlier | Pareto efficiency for the concave order and multivariate comonotonicity[END_REF] for an extension to multivariate risks). For suitable law invariant (or quantile-based) and concave utilities, the abstract risk-sharing problem may be brought down to the maximization of some concave functional over the set of comonotone allocations. When utilities are integral functionals of the quantile (which covers the rank-dependent utility case and more generally the so-called rank-linear utility case), this enables one to rewrite the problem as a tractable variational problem subject to a comonotonicity constraint. This reformulation is a key point in the papers mentioned above to analyze and actually solve many instances of the twoagents case.

The theoretical analysis of the many agents-case is not significantly harder than in the two-agents case. However, the analysis of optimality conditions is more involved since many more cases may arise in the many agents case. Typically, optimal risk-sharing rules exhibit ranges of aggregated risk for which a subset of agents is fully insured by the others. In the two-agents case, optimal solutions are combinations of three regimes: agent 1 insures agent 2, agent 2 insures agent 1 or the solution is interior and thus given by some first-order condition. If there are more agents, there are many more possibilities for the comonotonicity constraint to be binding. We thus claim that the difficulty of the many agents case is in fact a matter of combinatorics and that the search for an efficient computational scheme is therefore natural.

The paper is organized as follows. In section 2, we reformulate a class of risk-sharing problems as tractable variational problems subject to a comonotonicity constraint. More precisely, some notations and preliminaries are given in paragraph 2.1, the two-agents case is then addressed in paragraph 2.2, we finally show a stability by sup-convolution result in paragraph 2.3 which enables us to reformulate the risk-sharing problem as a comonotonicity constrained variational problem. Optimality conditions for such problems are established in section 3. Finally, a simple and easy to implement algorithm is described in section 4 in which convergence is established and various numerical simulations are presented.

Reformulation of a class of risk-sharing problems

Let (Ω, F, P) be a nonatomic probability space i.e. a probability space such that there is no A ∈ F such that P(A) > 0 and P(B) ∈ {0, P(A)} for every B ∈ F with B ⊂ A. The nonatomicity of (Ω, F, P) is well-known to be equivalent to the existence of a uniformly distributed random variable on (Ω, F, P). Given a random variable X on (Ω, F, P), the law of X is denoted L(X). If X and Y are two random variables X on (Ω, F, P), we shall denote by X ∼ Y the fact that L(X) = L(Y ).

In the sequel, we will only consider essentially bounded random variables -or risks-on (Ω, F, P). We will denote for short by L ∞ + the nonnegative cone of L ∞ (Ω, F, P).

Let C be the class of utility functions V : L ∞ + → R that :

• are strictly concave,

• are monotone, i.e. V (X) ≥ V (Y ) whenever (X, Y ) ∈ L ∞ + × L ∞ + and X ≥ Y a.s., • law invariant, i.e. V (X) = V (Y ) whenever (X, Y ) ∈ L ∞ + × L ∞ + and L(X) = L(Y ),
• satisfy the following Fatou property:

V (X) ≥ lim sup n V (X n ) whenever (X n ) n is bounded in L ∞
+ and X n converges a.s. to X.

Given an aggregate risk X 0 ∈ L ∞ + and d + 1 agents, i = 1, ..., d + 1, each of whom having a utility V i in the class C, the optimal risk-sharing of X 0 among those agents is determined by solving the sup-convolution problem:

d+1 i=1 V i (X 0 ) := sup d+1 i=1 V i (X i ), X i ∈ L ∞ + , d+1 i=1 X i = X 0 . (2.1)
The sup-convolution d+1 i=1

V i is said to be exact if the previous supremum is attained (and then it is attained at a unique point by strict concavity). Our aim is to reformulate the previous risk-sharing problem in a more tractable way using the notion of comonotone allocations. Before we do so, we shall need some preliminaries.

Preliminaries

Let X be a bounded random variable on (Ω, F, P) and let F X (t) = P(X ≤ t), t ∈ R denote its distribution function. The generalized inverse of F X (or quantile function of X) is defined by:

F -1 X (t) = inf{z ∈ R : F X (z) > t}, for all t ∈ (0, 1). (2.2)
The random variable X and the quantile function F -1 X are then equimeasurable in the sense that for every continuous function ϕ one has:

E(ϕ(X)) = 1 0 ϕ(F -1 X (t))dt.
Let us now recall the well known Hardy-Littlewood inequality (see [START_REF] Hardy | Inequalities[END_REF], [START_REF] Lorentz | An inequality for rearrangements[END_REF], [START_REF] Cambanis | Stout Inequalities for Ek(X, Y ) when the marginals are fixed[END_REF]).

Proposition 2.1. Let X and Y be in L ∞ + then one has

E(XY ) ≤ 1 0 F -1 X (t)F -1 Y (t)dt = sup Z∼Y E(XZ)
and for every concave u : R → R one has

E(u(X -Y )) ≤ 1 0 u(F -1 X (t) -F -1 Y (t))dt.
Definition 2.2. Let X and Y be in L ∞ (Ω, F, P), then X dominates Y in the sense of second order stochastic dominance (which we denote by X Y ) if E(u(X)) ≥ E(u(Y )) for every concave increasing function u : R → R.

Various characterizations of stochastic dominance are well-known (see [START_REF] Rothschild | Increasing Risk, I. A Definition[END_REF]) among which:

X Y ⇔ t 0 F -1 X (t)dt ≥ t 0 F -1 Y (t)dt, ∀t ∈ [0, 1],
which is equivalent to

1 0 g(t)F -1 X (t)dt ≤ 1 0 g(t)F -1 Y (t)dt, ∀g bounded nondecreasing. (2.3)
The next statement is originally due to Ryff [START_REF] Ryff | Extreme points of some convex subsets of L 1[END_REF] in the case (Ω, F, P) is [0, 1] endowed with its Borel field and the Lebesgue measure, we believe it is well-known but give a short proof for the sake of completeness:

Lemma 2.3. Let X and Y be in L ∞ + , then X Y if and only if there is a sequence Z n of the form Z n = Nn i=1 λ n i Y n i with λ n i ≥ 0, Nn i=1 λ n i = 1 and each Y n i ∼ Y , such that Z n converges a.s. to X.
Proof. Assume first that Z n is of the form mentioned above and Z n converges a.s. to X, then for each concave u, one has E(u(Z n )) ≥ E(u(Y )) and one concludes with Lebesgue's dominated convergence theorem (the Z n 's being uniformly bounded).

Conversely, assume X Y and let us prove that X ∈ K where K is the closed convex hull in

L 1 of the set {Z ∈ L 1 , Z ∼ Y } (clearly K ⊂ L ∞ + since Y ∈ L ∞ + ). If X /
∈ K it follows from Hahn's-Banach theorem that there exists P ∈ L ∞ and ε > 0 such that

E(P X) ≥ sup Z∼Y E(P Z) + ε.
With Hardy-Littlewood's inequality, this yields

1 0 F -1 P F -1 X ≥ 1 0 F -1 P F -1 Y + ε
which contradicts X Y because of inequality (2.3) recalled above. We thus deduce that X is the L 1 (and thus also a.s. taking a subsequence if necessary) limit of a sequence Z n of the form mentioned in the statement of the lemma.

As a consequence of the previous lemma, we have the following compatibility result which is originally due to Dana [START_REF] Dana | A representation result for concave Schur concave functions[END_REF]:

Lemma 2.4. If a utility function V is in the class C then it is compatible with second order stochastic dominance that is V (X) ≥ V (Y ) whenever (X, Y ) ∈ L ∞ + × L ∞ + and X Y .
Proof. Assume X Y , using lemma 2.3, X is an a.s. limit of a sequence Z n as in lemma 2.3. Since V is concave and law invariant, V (Z n ) ≥ V (Y ) and by the Fatou property, we get V (X) ≥ V (Y ).

A notion that will play an important role in the sequel is that of comonotonicity that we now recall, Definition 2.5. Let X 1 and X 2 be two real-valued random variables on (Ω, F, P), then the pair

(X 1 , X 2 ) is called comonotone if (X 1 (ω ′ ) -X 1 (ω))(X 2 (ω ′ ) -X 2 (ω)) ≥ 0 for P ⊗ P-a.e. (ω, ω ′ ) ∈ Ω 2 .
A family of random variable (X 1 , ..., X d+1 ) on (Ω, F, P) is said to be comonotone if (X i , X j ) is comonotone for every (i, j) ∈ {1, ..., d + 1} 2 .

Roughly speaking (X 1 , ..., X d+1 ) is comonotone if all the X i 's evolve in the same direction (or coevolve) and we recall that it is well-known that it is equivalent to the fact that each X i can be written as a nondecreasing function of the sum i X i (see for instance [START_REF] Denneberg | Non-additive Measures and Integral[END_REF]). Given X ∈ L ∞ + , the set of allocations (X 1 , ..., X d+1 ) ∈ L ∞ + that sum to X and are comonotone can then be parametrized by X i := x i (X) where the x i 's are nondecreasing function that sum to the identity map (on the range of X but this can be extended to the whole line if necessary). Obviously, each map x i is 1-Lipschitz and therefore, by Ascoli's theorem, this parametrization of comonotone allocations with a fixed sum is compact.

Comonotonicity is well-known to be a key property in risk-sharing problems as soon as agents have preferences that are compatible with second order stochastic dominance. Indeed, an important result of Landsberger and Meilijson [START_REF] Landsberger | Comonotone allocations, Bickel Lehmann dispersion and the Arrow-Pratt measure of risk aversion[END_REF] states that any allocation is dominated (for second order stochastic dominance) by a comonotone one.

The two-agents case

In this paragraph, we study in details the case of two agents (d = 1) by a rearrangement argument taken from [START_REF] Carlier | Pareto efficiency for the concave order and multivariate comonotonicity[END_REF].

Proposition 2.6. Let V 1 and V 2 be two utilities in the class C and X 0 ∈ L ∞ + , there exists a unique (X 1 , X 2 ) ∈ L ∞ + × L ∞ + such that X 1 + X 2 = X 0 and V 1 V 2 (X 0 ) = V 1 (X 1 ) + V 2 (X 2 ).
Moreover (X 1 , X 2 ) is comonotone hence there exist two nondecreasing functions x 1 and

x 2 : [F -1 X 0 (0), F -1 X 0 (1)] → R such that x 1 (t) + x 2 (t) = t for all t ∈ [F -1 X 0 (0), F -1 X 0 (1)] and (X 1 , X 2 ) = (x 1 (X 0 ), x 2 (X 0 )).
Proof. Let us first suppose that X 0 has no atom (i.e. F X 0 is continuous, or equivalently

F -1 X 0 increasing). Let X 1 and X 2 be in L ∞ + such that X 1 + X 2 = X 0 , it follows from Jensen's inequality that E(X i | X 0 ) X i so that in the sup-convolution, it is enough to maximize over pairs (X 1 , X 2 ) that are measurable functions of X 0 . Let X 1 = f (X 0 ) and X 2 := X 0 -f (X 0 ) be such a pair. Since X 0 is nonatomic there exists a nondecreasing map (or monotone rearrangement) f such that f (X 0 ) ∼ f (X 0 ) with 0 ≤ f (X 0 ) ≤ X 0 (see [19]). Then set Y 1 := f (X 0 ) and Y 2 := X 0 -f (X 0 ). Since Y 1 ∼ X 1 , V 1 (Y 1 ) := V 1 (X 1 ), we now claim that V 2 (Y 2 ) ≥ V 2 (X 2
) and to prove it, thanks to lemma 2.4, it is enough to prove that Y 2 X 2 . Let u be a concave function, we have

E(u(Y 2 )) = E(u(X 0 -f (X 0 )) = 1 0 u(F -1 X 0 (t) -f (F -1 X 0 (t)))dt but since f (F -1 X 0 ) = F -1 X 1 , by proposition 2.1, we get E(u(X 2 )) = E(u(X 0 -X 1 )) ≤ 1 0 u(F -1 X 0 (t) -F -1 X 1 (t))dt = E(u(Y 2 )).
All this proves that given an admissible pair (X 1 , X 2 ) one may find a better one of the form (f (X 0 ), X 0 -f (X 0 )) where f is a non decreasing function on

[F -1 X 0 , F -1 X 0 ] → R such that 0 ≤ f (t) ≤ t for every t ∈ [F -1 X 0 (0), F -1 X 0 (1)] → R.
In particular, one can find a sequence of such maps f n such that

lim n V 1 (f n (X 0 )) + V 2 (X 0 -f n (X 0 )) = V 1 V 2 (X 0 ).
Since the functions f n are all nondecreasing and bounded by X 0 ∞ , it follows from Helly's theorem that there is some (not relabeled) sequence that converges pointwise to some function f , hence

(f n (X 0 ), X 0 -f n (X 0 )) is bounded in L ∞ and converges a.s. to (f (X 0 ), X 0 -f (X 0 )). The Fatou property guar- antees then that V 1 (f (X 0 )) + V 2 (X 0 -f (X 0 )) = V 1 V 2 (X 0 ).
This proves existence of a maximizer in the sup-convolution problem (that is the sup-convolution is exact at X 0 ), uniqueness follows from the strict concavity of V 1 and V 2 . It remains to show that (f (X 0 ), X 0 -f (X 0 )) is comonotone i.e. X 0 -f (X 0 ) is nondecreasing in X 0 . To prove this last claim, we apply the same trick as before : let g be nondecreasing such that 0 ≤ g(t) ≤ t and g(X

0 ) ∼ X 0 -f (X 0 ), then X 0 -g(X 0 ) f (X 0 ) so that V 1 (X 0 -g(X 0 )) + V 2 (g(X 0 )) ≥ V 1 (f (X 0 )) + V 2 (X 0 -f (X 0 ))
and therefore (X 0 -g(X 0 ), g(X 0 )) is optimal as well and by uniqueness this yields X 0 -f (X 0 ) = g(X 0 ). The supremum is thus uniquely attained at some comonotone pair.

Let us now treat the case where X 0 is arbitrary in

L ∞ + . A theorem of Ryff ([19]) enables us to write X 0 = F -1 X 0 (U 0 ) with U 0 uniformly distributed. For n ∈ N * then define X n 0 := F -1 X 0 (U 0 ) + n -1 U 0 , since X n
0 is nonatomic, we deduce from the previous step that there exists, for every n ≥ 0, a pair of nondecreasing 1-Lipschitz functions x n 1 and x n 2 summing to the identity such that

V 1 V 2 (X n 0 ) = V 1 (x n 1 (X n 0 )) + V 2 (x n 2 (X n 0 ))
. By Ascoli's Theorem, we may extract some converging (and not relabeled) subsequence from (x n 1 , x n 2 ) and denote (x 1 , x 2 ) its limit. Now, using the fact that V 1 and V 2 are monotone and that X n 0 ≥ X 0 , we have, for every

X ∈ L ∞ + such that X ≤ X 0 , V 1 (X 0 -X) + V 2 (X) ≤ V 1 (X n 0 -X) + V 2 (X) ≤ V 1 (x n 1 (X n 0 )) + V 2 (x n 2 (X n 0 )). Since (x n 1 (X n 0 ), x n 2 (X n 0 )
) is bounded and converges a.s. (and in fact even in L ∞ ) to (X 1 , X 2 ) = (x 1 (X 0 ), x 2 (X 0 )), we deduce from the Fatou property that the limsup of the right-hand side is less than

V 1 (x 1 (X 0 )) + V 2 (x 2 (X 0 )).
Taking the supremum with respect to X we conclude that the comonotone pair (x 1 (X 0 ), x 2 (X 0 )) is the unique maximizer in the sup-convolution problem for X 0 .

Stability by sup-convolution and existence for d + 1 agents

In this section, our aim is to generalize the existence and characterization result of proposition 2.6 to the case of d + 1 agents. To do so, let us prove the following stability by sup-convolution result.

Proposition 2.7. Let V 1 and V 2 be two utilities in the class C, then V 1 V 2 also belongs to C. Proof. Let X 0 , Y 0 ∈ L ∞ + × L ∞ + , then by proposition 2.6, V 1 V 2 (X 0 ) and V 1 V 2 (Y 0
) are exact that is, there exist nondecreasing functions x 1 , x 2 , and y 1 , y 2 summing to the identity and such that

V 1 V 2 (X 0 ) = V 1 (x 1 (X 0 )) + V 2 (x 2 (X 0 )) and V 1 V 2 (Y 0 ) = V 1 (y 1 (Y 0 )) + V 2 (y 2 (Y 0 )). Let us prove step by step that V 1 V 2 belongs to C.
• Strict concavity: let λ ∈ (0, 1), and assume X 0 = Y 0 , we thus have

(x 1 (X 0 ), x 2 (X 0 )) = (y 1 (Y 0 ), y 2 (Y 0 )), the strict concavity of V 1 and V 2 then yields V 1 V 2 (λX 0 + (1 -λ)Y 0 ) V 1 (λx 1 (X 0 ) + (1 -λ)y 1 (Y 0 )) +V 2 (λx 2 (X 0 ) + (1 -λ)y 2 (Y 0 )) > λ[V 1 (x 1 (X 0 )) + V 2 (x 2 (X 0 ))] +(1 -λ)[V 1 (y 1 (Y 0 )) + V 2 (y 2 (Y 0 ))] so that V 1 V 2 is strictly concave. • Monotonicity: assume X 0 Y 0 and let X ∈ L ∞ + be such X X 0 , since V 1 and V 2 are monotone, we have V 1 (X 0 -X) + V 2 (X) V 1 (Y 0 -X) + V 2 (X) V 1 V 2 (Y 0 )
taking the supremum with respect to X, we deduce that V 1 V 2 is monotone.

• Law invariance property: assume that X 0 ∼ Y 0 . We then have by

law invariance of V 1 and V 2 V 1 V 2 (X 0 ) = V 1 (x 1 (Y 0 )) + V 2 (x 2 (Y 0 )) V 1 V 2 (Y 0 ), reversing the role of X 0 and Y 0 then yields V 1 V 2 (X 0 ) = V 1 V 2 (Y 0 ).
• Fatou property: let (X n ) n be a bounded sequence in L ∞ + that converges a.s. to a limit X. Then for every n ≥ 1, there exist nondecreasing and 1-Lipschitz nonnegative functions

x n 1 , x n 2 such that V 1 V 2 (X n ) = V 1 (x n 1 (X n )) + V 2 (x n 2 (X n )),
and, by Ascoli's Theorem, up to a not relabeled subsequence, x n 1 and x n 2 converge uniformly to some x 1 and x 2 so that x n i (X n ) converges a.s. to x i (X). We thus have lim sup

n V 1 V 2 (X n ) lim sup n V 1 (x n 1 (X n )) + lim sup n V 2 (x n 2 (X n )) ≤ V 1 (x 1 (X)) + V 2 (x 2 (X)) ≤ V 1 V 2 (X)
from which we deduce the Fatou property for

V 1 V 2 .
Inductively, one immediately deduces from the previous proposition that if V 1 , ...., V d+1 all belong to the class C then so does

d+1 i=1 V i . Theorem 2.8. Let X 0 ∈ L ∞ + and V 1 , ..., V d+1 be in the class C. Then there exists a unique X = (X 1 , ..., X d+1 ) ∈ (L ∞ + ) d+1 such that d+1 i=1 X i = X 0 and d+1 i=1 V i (X 0 ) = d+1 i=1 V i (X i ).
Moreover X is comonotone hence there exist d+1 nondecreasing functions

x i : [F -1 X 0 (0), F -1 X 0 (1)] → R such that d+1 i=1 x i (t) = t for all t ∈ [F -1 X 0 (0), F -1 X 0 (1)
] and X = (x 1 (X 0 ), ..., x d+1 (X 0 )).

Proof. For the sake of simplicity, let us prove the Theorem for d + 1 = 3. Introducing the notation W 2 = V 2 V 3 , the maximization problem reads as:

sup{V 1 (X 1 ) + W 2 (Y 2 ); X 1 ∈ L ∞ + , Y 2 ∈ L ∞ + , X 1 + Y 2 = X 0 }
. By proposition 2.7, we know that W 2 ∈ C. Proposition 2.6 then says that there exists a unique solution (

X 1 = x 1 (X 0 ), Y 2 = y 2 (X 0 )) ∈ L ∞ + × L ∞ + such that x 1 ,y 2 are nondecreasing, x 1 (X 0 ) + y 2 (X 0 ) = X 0 , and V 1 W 2 (X 0 ) = V 1 (x 1 (X 0 )) + W 2 (y 2 (X 0 )).
We then solve the sub-problem

sup{V 2 (X 2 ) + V 3 (X 3 ); X 2 ∈ L ∞ + , X 3 ∈ L ∞ + , X 2 + X 3 = y 2 (X 0 )}.
By the same arguments, there exists a unique comonotone solution (

X 2 = z 2 (y 2 (X 0 )), X 3 = z 3 (y 2 (X 0 )) ∈ L ∞ + × L ∞
+ satisfying the same properties as in proposition 2.6. Finally, it is easy to see that ( V i (X 0 ), with the desired properties.

X 1 = x 1 (X 0 ), X 2 = x 2 (X 0 ), X 3 = x 3 (X 0 )),
Inductively, one can prove the theorem for every d. Remark 2.9. One may wonder if it is really necessary to proceed inductively to prove the previous result. On the one hand, we believe that proving stability of the class C by supremal convolution (that is by aggregation) as we did in proposition 2.7 has its own interest. On the other hand, the rearrangement trick we used in the two-agents case does not generalize to more agents. Remark 2.10. In the previous statements, we have used strict concavity of the elements of C for uniqueness purpose only. Strict concavity is in fact quite restrictive since it rules out the case of monetary risk measures. We now claim that strict concavity, though convenient because it implies uniqueness, is in fact not necessary to obtain the existence of at least one optimal comonotone solution in the supremal convolution problem. To see this, it is enough to proceed by approximation and use the compactness of comonotone allocations. By the way, in the next section, we will treat numerically the case of concave and non strictly concave utilities by using such strictly concave perturbations.

Reformulation

It follows from the previous results that the risk-sharing problem (2.1) between d + 1 agents, each of whom has a utility in C, and for aggregate risk X 0 , can be simply reformulated as

sup x∈A J(x), with J(x) := d i=1 V i (x i (X 0 )) + V d+1 X 0 - d i=1 x i (X 0 ) (2.4)
where

A := {x ∈ W 1,∞ ([a, b], R d ) : ẋ ∈ ∆ a.e., x(a) ∈ S}, a := F -1 X 0 (0), b := F -1 X 0 (1)
, and ∆ and S are the two simplices:

∆ := {u ∈ R d + , d i=1 u i ≤ 1}, S := {x ∈ R d + , d i=1 x i ≤ a}.
This reformulation is more tractable than the initial one since maximization is now performed over a more concrete set of functions over an interval, namely the convex and compact (for the uniform norm) set A. As we shall see in the next paragraph, the fact that constraints are given by simplices will be very convenient to express necessary conditions. Of course this reformulation is not very helpful if we keep the previous level of generality on utility functions. This is why we shall now further specify the utilities and restrict them to a subclass for which they have a simple expression in terms of quantiles. From now on, we shall restrict utilities to belong to the class of Rank-Linear Utilities (RLU for short), such utilities are of the form

V L (X) := 1 0 L(t, F -1 X (t))dt. (2.5)
This class of utilities is already quite large since it contains expected utilities as well as the Rank-Dependent Utilities defined by Choquet expectations.

Utilities in this class are obviously law invariant and satisfy the Fatou property as soon as L is continuous (say). RLU have been studied in details in [START_REF] Carlier | Two-Persons Efficient Risk-Sharing and Equilibria for Concave Law-Invariant Utilities[END_REF], where it is proved in particular that for V L defined by (2.5) (with a C 2 function L to simplify), the following statements are equivalent:

1. V L is compatible with second order stochastic dominance and monotone,

2. ∂ x L ≥ 0, ∂ xx L ≤ 0 and ∂ tx L ≤ 0 on [0, 1] × R 3. V L is concave, monotone and σ(L ∞ (Ω), L 1 (Ω)) upper semi-continuous.
Now for a utility of the form (2.5), if x is a nondecreasing function and X 0 is nonatomic with an increasing distribution function (say), one has

V L (X 0 ) = 1 0 L(t, x(F -1 X 0 (t))dt = b a L(F X 0 (s), x(s))dµ 0 (s),
where µ 0 is the law of X 0 . In this framework, (2.4) takes the form of a simple variational problem with

J(x) := b a d i=1 L i (F X 0 (s), x i (s)) + L d+1 (F X 0 (s), s - d i=1
x i (s)) dµ 0 (s).

(2.6) The next paragraphs are precisely devoted to the theoretical and numerical study of such problems (under appropriate regularity and concavity assumptions).

Optimality conditions

In this section, our aim is to give necessary and sufficient conditions for the problem sup where we assume that

F ∈ C([a, b] × R d , R) is such that F (t, .
) is strictly concave and differentiable for every t ∈ [a, b] and ∇ x F is continuous in its two arguments. Under these assumptions, problem (3.1) admits a unique solution x that is characterized as follows: Theorem 3.1. Problem (3.1) admits a unique solution x that is characterized by the following:

• x ∈ A,
• for a.e. t, ẋi (t) = 0 for every i ∈ {1, ..., d} with i / ∈ I(t) where

I(t) := {j ∈ {1, ..., d} : p j (t) = min i=1,...,d p i (t)}
and p is the adjoint variable:

p(t) = - b t ∇ x F (s, x(s))ds. (3.2)
• x i (a) = 0, ∀i ∈ {1, ...d} with i / ∈ I(a).

Proof. Existence and uniqueness follow from standard arguments. Let us denote by x the solution of (3.1), by concavity of J and convexity of A, x is characterized by the variational inequalities

J ′ (x) • (y -x) = b a ∇ x F (t, x(t)) • (y(t) -x(t))dt ≤ 0, ∀y ∈ A. (3.3)
Now for y ∈ A, defining p by (3.2) and integrating by parts yields

J ′ (x) • y = - b a p • ẏ -p(a) • y(a) so that (3.3) becomes b a p • ẋ + p(a) • x(a) = inf y∈A b a p • ẏ + p(a) • y(a) , (3.4)
now the rightmost member of the previous identity is easy to compute by pointwise minimization over the simplex and is achieved exactly at those y's in A such that for a.e. t ∈ [a, b], ẏi = 0 as soon as i / ∈ I(t) and y i (a) = 0 for i / ∈ I(a).

In the algorithm of next section, we shall take full advantage of the fact that (3.4) amounts to linear programming over the simplex and is therefore explicit. In particular, for x ∈ A and p the associated adjoint variable:

p(t) = - b t ∇ x F (s, x(s))ds, (3.5)
this characterizes the set of solutions of the linear programming problem

sup y∈A J ′ (x) • y. (3.6)
In particular, a solution of (3.6) is given explicitly by:

y(t) = y(a) + t a u(s)ds, t ∈ [a, b], (3.7)
with

y i (a) = a #Ip(a) if i ∈ I p (a) 0 otherwise, (3.8)
• start from x 0 ∈ A and define p 0 as the corresponding adjoint:

p 0 (t) = - b t ∇ x F (s, x 0 (s))ds • given (x k , p k ) define y k by y k (t) = y k (a) + t a u k (s)ds, t ∈ [a, b],
with

y k i (a) = a #I p k (t) if i ∈ I p k (a)
0 otherwise and

u i (t) = 1 #I p k (t) if i ∈ I p k (t) 0 otherwise. ( 4.1)
Then set

x k+1 := x k + ρ k (y k -x k )
where

ρ k = argmax ρ∈[0,1] J(x k + ρ(y k -x k ))
which can be computed by simple dichotomy method.

Finally let p k+1 be the adjoint state associated to x k+1 :

p k+1 (t) = - b t ∇ x F (s, x k+1 (s))ds.
The convergence of this simple ascent algorithm is then given by: Theorem 4.1. The sequence x k generated by the previous algorithm converges to x, the solution of (3.1) as k → ∞.

Proof. Since A is compact (for the uniform topology) and the sequence (x k ) k has values in A, it is enough to show that x is the unique cluster point of (x k ) k .

Let x be such a cluster point, and let x kn be a subsequence converging to x, taking a subsequence if necessary we may also assume that x kn+1 converges to some z. Since G has a closed graph, z ∈ G(x) i.e. z has of the form z = x+ρ(y-x) with y maximizing J ′ (x)•y over A and ρ maximizing J(x+ρ(y-x)) over [0, 1]. By construction J(x k ) is nondecreasing, consequently one has J(x) = J(z). Now, if x = x, then by lemma 3.2, one has x / ∈ G(x) so that ρ > 0 and J ′ (x)•(y-x) > 0 which implies that J(z) = J(x+ρ(y-x)) > J(x), giving the desired contradiction. This proves that x is the only cluster point of (x k ) k and thus that the whole sequence (x k ) k converges to x. Remark 4.2. As already mentioned, it may seem very restrictive to assume strict concavity since it excludes monetary utility functions that are relevant in finance. In fact, as we already explained, strict concavity is by no means essential for existence of an optimal comonotone risk-sharing rule. It is however a crucial assumption for convergence of our algorithm. To deal with problems of the form (3.1) with a concave but not strictly concave J it is enough to consider the perturbed problem

inf x∈A J ε (x) with J ε (x) := J(x) + εΦ(x) (4.2)
where ε > 0 is a small perturbation parameter and Φ is a nice (by this, we mean a smooth integral functional) strictly concave perturbation. Then, our algorithm enables us to compute the solution x ε of the perturbed problem (4.2). Now, it is easy to see that under mild assumptions on the perturbation, x ε converges as ε → 0 + to the unique maximizer of Φ among the maximizers of J on A. This perturbation scheme thus naturally selects a particular optimal risk-sharing rule that can be approximated by the algorithm described above.

Numerical experiments

In this final section, we implement the algorithm introduced in section 4.1 and illustrate it on various examples. The numerical results that we obtain confirm the theoretical convergence established in theorem 4.1 and show that our algorithm converges really fast on every example we tested.

More precisely, we first study three particular cases for a uniform distribution of the total risk X 0 on [0, 1]. In a fourth example, we compare the solutions for different distributions of X 0 . The last test aims at showing that our numerical procedure is robust for approaching solutions of risk sharing problems with concave but not strictly concave utilities. But before doing so, we would like to mention that the algorithm gives very good results, in particular very quick convergence. Indeed, it is shown on Figure 1-(a) how fast the convergence of the value is. We can also see on Figure 1-(b) that the difference between two successive iterations converges to zero with about the same speed. In all our tests, the algorithm numerically converges in very few iterations (between 5 and 15).

In the following examples, we always take a = 0 and b = 1 (which is in fact just a matter of normalization).

L i (t, x) = (1 -t) α i (1 + x) β i ,
for i = 1, ..., d+1 which corresponds to RDU with power distortion and power utility indices. It is easy to see that for α i ≥ 0 and 0 < β i < 1, such functions lead to utilities that belong to the class C. In this first numerical experiment, we look at the case d + 1 = 5 and we choose α := (2.8, 1.6, 2.2, 1.05, 1.3) and β := (0.58, 0.52, 0.59, 0.53, 0.57). Figure 2 shows the optimal risk-sharing for these data and represents the cumulative sharing in the sense that the size of the risk supported by an agent is the difference between her labeled line (the lines are labeled on the right hand side of the picture) and the nearest sub-line. Looking at Figure 2, one can for instance notice that Agent 4 fully insures the others for any risk with value smaller than 0.34. Then, for risk values smaller than about 0.7, Agents 3 and 4 insure the others and finally, for large values of the risk, every agent bears a part of the risk. Second case: logarithmic RLU. In this example, we make a different choice for the functions L i . We take two logarithmic functions h 1 (t, x) := log(t 1/2 + x + 1) and h 2 (t, x) := log(t 2/3 + x + 1). Note that this RLU example does not correspond to RDU agents. We look at the case of 8 agents, five of them being of type L i = h 1 (i = 3, 4, 5, 6, 7) and three of them of type L i = h 2 (i = 1, 2, 8). The graph of the optimal risk-sharing is represented in Figure 3. The continuous lines correspond to agents of type h 1 and the discontinuous ones to those of type h 2 . We can see that for risk values between 0 and 0.3, the three agents of type h 2 fully insure the others. Numerically, agents of the same type support the same proportion of risk which agrees with the theory.

Third case: mixing RLU and RDU agents. We now turn to a case involving more agents than in the previous ones: d+1 = 20. The functions L i are chosen between L i (t, x) = (1 -t) α i (1 + x) β i , for some α i and β i satisfying the same assumptions as before, and L i (t, x) = log(t γ i + x + 1), γ i ∈ (0, 1). The graph of the solution is represented on Figure 4.

Fourth case: different densities for X 0 . We now aim to see the impact of the law of the aggregate risk X 0 on the optimal individual risks.

To do so, we take d+1 = 4 and functions L i (t, x) = (1-t) α i (1+x) β i , with the particular choice α = (1.16, 1.89, 1.21, 2.92) and β = (0.52, 0.58, 0.54, 0.59).

We start with a comparison for three truncated and normalized gaussians with variance σ 2 = 0.01. More precisely, we plot on Figure 5 the graphs of the density for three different values of the average µ = 0.1, 0.5, 0.9, and the corresponding solutions. Let us note for instance that Agent 1 fully insures the others for values lower than 0.57 when the average is small (0.1). This phenomenon disappears for larger values of µ. We can also remark that Agent 5 supports a significant proportion of risk only for high values of X 0 and µ. We have drawn on Figure 6 the graphs of three other densities with the corresponding optimal risk-sharing. The densities we study are, for t ∈ [0, 1]:

• (a) -(b): f X 0 (t) = 1 (Uniform density),

• (c) -(d): f X 0 (t) = λe -λt + e -λ (Exponential-like density),

• (e) -(f ): f X 0 (t) = t 2 + 2 3 (Quadratic density, say). Fifth case: non stricly concave utilities. To close this series of numerical tests, we aim to show that our algorithm is actually able to treat concave but not strictly concave utilities. This is an important extension which is motivated by Choquet expectations and monetary risk measures Figure 8: Distance of the approximations to the numerical limit such as Expected Shortfalls. The idea is to make a small strictly concave perturbation of the utilities as explained in remark 4.2. In this example, we look at a case with 5 agents and a gaussian distribution of the aggregate risk X 0 (its density is represented on Figure 7-(a)). We then take perturbed utilities of the form L ε i (t, x) = (1 + εx β i + (1 -ε)x)(1 -t) α i , with 0 < β i < 1 and α i > 0. We then let ε → 0 to recover the (linear) Choquet expectation. We show the numerical results for ε = 1 (Figure 7-(b)) , ε = 0.05 (Figure 7-(c)) and ε = 10 -5 (Figure 7-(d)).

We have finally computed the uniform distance between the approximations and the numerical limit. If x ε denotes the numerical solution for a fixed ε, we define this quantity by ||x 10 -5 -x ε || ∞ , and we give the results on Figure 8, for some values of ε in [0.02, 0.05]. We therefore observe a numerical convergence.

  where x 2 := z 2 • y 2 and x 3 := z 3 • y 2 , is the unique solution of the sup-convolution 3 i=1

F

  x∈A J(x) with J(x) := b a (t, x(t))dt (3.1)

Figure 1 :

 1 Figure 1: Numerical convergence for a 20 agents risk-sharing problem

Figure 2 :

 2 Figure 2: Risk Sharing for 5 agents

Figure 3 :Figure 4 : 9 Figure 5 :

 3495 Figure 3: Risk Sharing for 8 agents

  Optimal Risk Sharing for the Uniform density (c) Exponential density with λ = 4 (d) Optimal Risk Sharing for the Exponential density (e) Quadratic density (f) Optimal Risk Sharing for the Quadratic density

Figure 6 :

 6 Figure 6: Plots of the solutions for three different distributions of the total risk X 0 .

( a )

 a Density of X 0 (b) ε = 1 (c) ε = 0.05 (d) ε = 10 -5

Figure 7 :

 7 Figure 7: Approximation of the solution of the problem for a linear utility

Acknowledgements

The authors wish to express their gratitude to Rose-Anne Dana for interesting discussions and helpful suggestions.

otherwise.

(3.9) and I p (t) := {j ∈ {1, ..., d} : p j (t) = min i=1,...,d p i (t)}.

For every x ∈ A, set

In the next section, we shall need the following elementary result:

Lemma 3.2. The set valued map x ∈ A → G(x) has a closed graph and x, the solution of (3.1), is characterized by the condition x ∈ G(x).

Proof. Let (x n ) n ∈ A converge to some x and z n = x n + ρ n (y n -x n ) converge to some z and let us prove that z ∈ G(x). Taking subsequences if necessary, we may assume that ρ n converges to some ρ ∈ [0, 1] and y n converges uniformly to some y ∈ A so that z = x + ρ(y -x). Using the fact that J is of class C 1 enables us to pass to the limit in

and therefore to deduce that z ∈ G(x).

It is obvious that x, the solution of (3.1), is such that x ∈ G(x). Now if x ∈ G(x) then one has either ρ = 0 or x = y for some ρ ∈ [0, 1] and some y ∈ A as in the definition of G(x). If x = y then J ′ (x) • h is maximized over A for h = x which is exactly the variational inequality characterizing the solution of (3.1). If ρ = 0 and y = x then J(x + λ(y -x)) is maximized over [0, 1] for λ = 0 so that J ′ (x) • (y -x) ≤ 0 which again means that J ′ (x) • h is maximized over A for h = x.

Algorithm and numerical results

Algorithm and convergence

The algorithm we propose to solve (3.1) is a simple optimal step gradient ascent-like method. The fact that the constraints are easy to handle relies again on the fact that linear problems over the simplex are very simple. Our algorithm is defined as follows :


Alessandro Zin

Olivier Sename
email: olivier.sename@gipsa-lab.grenoble-inp.fr

Luc Dugard
email: luc.dugard@gipsa-lab.grenoble-inp.fr

A feedback-feedforward suspension control strategy for global chassis control through anti-roll distribution

Keywords: vehicle, active suspension, yaw control, anti-roll distribution, H∞, Linear Parameter Varying systems

A new global chassis control strategy developed in two steps is proposed using the suspension actuators only. First, a robust H∞ Linear Parameter Varying (LPV) suspension controller is designed to ensure robust performance w.r.t parameter uncertainties. The interest of such a controller relies on its inherent structure which allows for performance adaptation through the variation of a scalar coefficient c0. The second step consists in adding a feedforward compensator which controls the anti-roll distribution of the suspension forces to improve the vehicle safety in the presence of longitudinal and lateral disturbances. The control scheme consists in controlling the yaw rate through a model reference control. The global chassis control then relies on the adaptation of the suspension performance (comfort/safety), according to the anti-roll distribution. The interest of the methodology is emphasized through numerical simulations on a full vehicle non linear model, including chassis and tires.

Introduction

Over the past few years, automotive engineering has been characterized by rapid growth in active systems. In particular, efforts to improve driving safety and comfort under all driving conditions have been naturally focused on the chassis behavior. During the last decade, active control of vibration isolation in automotive vehicle has known a growing interest in research [START_REF] Oustaloup | The CRONE suspension[END_REF][START_REF] Hrovat | Survey of advanced suspension developments and related optimal control applications[END_REF][START_REF] Sammier | Skyhook and H ∞ control of semi-active suspensions: Some practical aspects[END_REF]. In particular the interest of the H ∞ approach has been emphasized in the context of active suspension control of full vehicle models [START_REF] Park | An H ∞ controller for active suspensions and its robustness based on a full car model[END_REF][START_REF] Lu | Multiobjective optimal suspension control to achieve integrated ride and handling performance[END_REF][START_REF] Gaspar | Design of robust controllers for active vehicle suspension using the mixed µ synthesis[END_REF].

However active safety and comfort are known to stumble over conflicting interests for the conventional chassis. To tackle this problem, the usual approach consists of an independent optimisation of the different chassis subsystems (in particular brakes, suspension actuators, steering, etc.). However, the complex and multivariable inherent nature of the vehicle makes such an approach inoperative. The target relies today on global chassis control [START_REF] Shibahata | Progress and future direction of chassis control technology[END_REF] and on centralized control [START_REF] Valtolina | A co-ordinate approach for the control of road vehicles[END_REF] strategies, to account for these coupling phenomena between steering, suspension and braking. In [START_REF] Zin | Active comfort and handling improvement with a "3D" bicycle model[END_REF] the improvement in security is shown using steering and suspension actuators. In [START_REF] Gaspar | The design of a combined control structure to prevent the rollover of heavy vehicle[END_REF][START_REF] Gáspár | The design of an integrated control system in heavy vehicles based on an LPV method[END_REF][START_REF] Gáspár | A fault-tolerant rollover prevention system based on LPV method[END_REF]] some control strategies using braking actuators, active anti-roll bars and/or active suspensions are developed in the LPV/H ∞ context for heavy vehicles. Finally in [START_REF] Nouillant | Feedforward and CRONE feedback control strategies for automobile abs[END_REF][START_REF] D'andrea-Novel | Collaboration between braking torques and active suspension forces to control a vehicle[END_REF][START_REF] Gáspár | Toward global chassis control by integrating the brake and suspension systems[END_REF][START_REF] Poussot-Vassal | Attitude and handling improvements through gain-scheduled suspensions and brakes control[END_REF] the collaboration between the braking and suspension actuators is investigated for vehicles, and in [START_REF] Corno | On optimal motorcycle braking[END_REF] for motorcycles..

In this paper a new feedback-feedforward global chassis control strategy is proposed based on suspension actuators only. First, a robust H ∞ Linear Parameter Varying (LPV) suspension controller is designed to ensure robust performance w.r.t parameter uncertainties. The interest of such a controller also relies on its inherent structure which allows for performance adaptation through the variation of a scalar coefficient c 0 . The second step consists in an original feedforward strategy which aims at limiting the sprung mass displacements in braking or acceleration phases, and, as well, at controlling the anti-roll distribution of the suspension forces to improve the vehicle safety in the presence of longitudinal and lateral disturbances. The control scheme for vehicle safety consists in controlling the yaw rate through a model reference PI controller. The global chassis control then relies on the adaptation of the suspension performance (comfort/safety) according to the anti-roll distribution. The interest of the methodology is emphasized through numerical simulations on a full vehicle non linear model including chassis and tires.

The outline of the paper is as follows. In section 2, the full vehicle model is described. Then in section 3 the LPV/H ∞ feedback controller for active suspensions is developed. Section 4 is dedicated to the feed-forward strategy for global chassis control through anti-roll distribution. The paper ends with some concluding remarks.

Vehicle modelling

A full-car model is here described for the analysis of the dynamical car behavior, road holding and comfort. It is made of two subsystems representing the horizontal and vertical dynamics respectively, thereby extending the one developed in [START_REF] Zin | A nonlinear vehicle bicycle model for suspension and handling control studies[END_REF].

The main modelling assumptions are the use of rigid solids, neglection of longitudinal efforts, autoalignment tire moments, the kinematical effects due to the geometry of the suspensions, the gyroscopical effects on the suspended mass due to the yaw rate, the anti-roll bar effects and the aerodynamical effects. Moreover, the inertial axes of the suspended mass are assumed to cross each other at the center of gravity and are assumed to be parallel to the ground.

The horizontal model is then written as (cf figure 1) :

(1)

                             a x (t) = vx (t) -v y (t) ψ(t), m a y (t) = m vy (t) + v x (t) ψ(t) = F ytf r (.) + F ytf l (.) cos δ f (t) +F ytrr (.) + F ytrl (.) + F yd (t), I z ψ(t) = F ytf r (.) + F ytf l (.) cos δ f (t) l f -F ytrr (.) + F ytrl (.) l r + M zd (t) + F ytf r (.) -F ytf l (.) sin δ f (t) t f .
where the notations are as follows:

• m = 4 m s + 2 m usf + 2 m usr : total vehicle mass, with m s : 1/4 of the total sprung mass, m usf : 1/2 of the front unsprung mass, m usr : 1/2 of the rear unsprung mass,

• I z = I zs + 2 m usf l 2 f + 2 m usr l 2 r : moment of inertia of the vehicle around the vertical axis with I zs : moment of inertia of the sprung mass around the vertical axis (in the sprung mass frame),

• l f (resp. l r ) : distance from the front axle (resp. rear) to the center of gravity of the car, t f (resp. t r ) : distance from the wheel center to the axle center at the front (resp. rear),

• F ytij (.) : front and rear, left and right, lateral tire forces (i = {f, r}), (j = {r, l}), F yd (t) : lateral disturbance force, M zd (t) : disturbance moment of inertia around the vertical axis,

• δ v (t): steering wheel angle, δ f (t) : front steering angle,

• C : location of projection of the center of gravity of the vehicle to the ground, : vehicle attitude,

• v x (t
• ψ(t) : yaw rate,

• a x (t) (resp. a y (t)) : longitudinal acceleration (resp. lateral) of C in the vehicle frame. The tire lateral forces are calculated through the well known tire static model of Paceijka as below [START_REF] Pacejka | The magic formula tyre model[END_REF][START_REF] Pacejka | Magic formula tyre model with transient properties[END_REF][START_REF] Milliken | Race car vehicle dynamics[END_REF]:

(2)

F ytij = D. sin(C. arctan(B.α ij -E.(B.α ij -arctan(B.α ij ))))
where, B, £C, £D and E are function of the normal force F ztij and the four wheel sideslip angles are given by:

(3)

               α f r (t) = arctan vy(t)+l f ψ(t) vx(t)+t f ψ(t) -δ f (t), α f l (t) = arctan vy(t)+l f ψ(t) vx(t)-t f ψ(t) -δ f (t), α rr (t) = arctan vy(t)-lr ψ(t) vx(t)+tr ψ(t) , α rl (t) = arctan vy(t)-lr ψ(t) vx(t)-tr ψ(t) .
The vertical tire forces F ztij (.) are therefore considered as exogenous model inputs and will be provided by the vertical model.

According to figure 2, the model exogenous inputs (1) are the disturbance forces and moments F yd (t), M zd (t), and the vertical tire forces F ztij (.), the steering wheel angle δ v (t) given by the driver and the variation of the longitudinal speed vx (t) (referred to as v x st in the figure). On the other hand, the vertical model is given on figure 3, and satisfies:

(4)

                                  
4 m s zs (t) = F sf r (.) + F sf l (.) + F srr (.) +F srl (.),

I xs θ(t) = -F sf r (.) + F sf l (.) t f + -F srr (.) + F srl (.) t r +m H a y (t), I ys φ(t) = F srr (.) + F srl (.) l r -F sf r (.) + F sf l (.) l f -m H a x (t), m usij zusij (t) = -F sij (.) + F ztij ,
with:

• I xs (resp. I ys ) : moment of inertia of the sprung mass around the longitudinal (resp. lateral) axis,

• H : height of the center of gravity of the car,

• z s (t) : vertical displacement of CG,

• θ(t) : roll angle of the sprung mass,

• φ(t) : pitch angle of the sprung mass In (4) F sij (.) is the vertical force delivered by the front/rear (i = {f, r}) right/left (j = {r, l}) suspension.

Passive suspension : in this case

(5) where F ksij (.) and F csij (.) are nonlinear functions given by the static maps in figure 4 for the front and rear suspensions (spring and damper). This passive case is considered as the reference non-controlled case for comparison.

F sij (.) = F ksij (.) + F csij (.),
Active suspension: in this case (6)

F sij (.) = F ksij (.) + u csij ,
where u csij is the force provided by the active damper. It will be detailed in the next section.

The absolute position of the deflection Z def ij (t), used in the static maps, is given by:

Z def ij (t) = Z sij -Z usij + z sij (t) -z usij (t) -l ksij ,
where l ksij is the nominal spring length, Z sij and Z usij are the absolute positions of the sprung and unsprung masses, respectively. These are calculated assuming Non linear characteristics of the suspension forces.

F zd (t) = M xd (t) = M yd (t) = a x (t) = a y (t) = z 0ij (t) = z sij (t) = z usij (t) = 0 in (4)
and can be obtained, if necessary, by inversion of the characteristics F ks (.) in figure 4, i.e. :

(7) (Z s -Z us -l ks ) = F -1 ks (F 0 ms ), (Z us -l kt ) = - (F 0 ms +F 0 mus ) kt .
where F 0 ms and F 0 mus are the static loads of m s and m us . Finally the deflection positions of the four suspensions are given by : (8)

       z sf r (t) = z s (t) -l f sin φ(t) -t f sin θ(t), z sf l (t) = z s (t) -l f sin φ(t) + t f sin θ(t), z srr (t) = z s (t) + l r sin φ(t) -t r sin θ(t), z srl (t) = z s (t) + l r sin φ(t) + t r sin θ(t).
Moreover the vertical tire forces are calculated using eq. ( 9) and ( 10), i.e. ( 9)

F t (.) = F kt (.) + F ct (.),
where F kt (.) and F ct (.) satisfy the following linear equations :

(10)

           Z def t (t) = Z us -z us (t) + z 0 (t) -l kt , żdeft (t) = żus (t) -ż0 (t) , F kt (.) = -k t Z def t (t), F ct (.) = -c t żdeft (t).
Following fig. 5, the exogenous inputs of the vertical model are the ground profile z 0ij (t), the longitudinal and lateral accelerations a x (t), a y (t) of C. The "Suspension characteristics" block includes, for the case of passive suspension, the force-deflection speed static map. For controlled suspensions, it corresponds to the controlled input forces [START_REF] Gáspár | Toward global chassis control by integrating the brake and suspension systems[END_REF]. Scheme of the vertical model.

An LPV feedback controller for active suspensions

In this section, some results presented in [START_REF] Zin | An LPV/H ∞ active suspension control for global chassis technology: Design and performance analysis[END_REF] are briefly recalled. The objective is to design an LP V /H ∞ controller scheduled according to the stiffness coefficient of the suspension spring k s (.) in order to improve performance and robustness. The scheduled design controller is based on the H ∞ theory applied to polytopic systems. The theoretical results used for control design can be found in [START_REF] Apkarian | Self scheduled H ∞ control of linear parameter-varying systems: a design method[END_REF].

The control design step is carried out using a quarter-car model and then applied to each of the four suspension active dampers, as shwon in figure 6.

The suspension control-oriented model

The control scheme is given in figure 7. The signal F 0 ms is the steady stade force caused by the weight of the sprung mass, which has to be known in order to evaluate the value of the steady state deflection (Z s -Z us -l ks ) in [START_REF] Gaspar | Design of robust controllers for active vehicle suspension using the mixed µ synthesis[END_REF]. This initial measure may be delivered by a pressure sensor installed in the active shock Quarter-car suspension models (passive and active).

absorber. Then it allows to get on-line the real-time value of k s (.) used in figure 7 for each damper, i.e.

( The controller K(s, k s ) is designed using the following LPV model where

) k s (Z def ij ) = - F ks (Z def ij ) Z def ij 11
(12) ẋ(t) = A a (k s ) x(t) + B 1 d(t) + B 2 u(t), with x(t) T = żs (t), z s (t), żus (t), z us (t) , d(t) T = F zd (t), z 0 (t) ,
where

A a (k s ) =       -c0 ms -ks(.)
      B 1 =     -1/m s 0 0 0 0 k t /m us 0 0     , B 2 =     1/m s 0 1/m us 0     .
In this model, c 0 is a constant damping coefficient, while k s (.) is a varying parameter. c 0 is here the nominal damping coefficient which could be chosen as the linear approximation of the damping. Then the control-oriented model uses u cs (t) = u(t) -c 0 żdeft (t) where u corresponds to the variation of the damping force compared to the nominal passive suspension.

Therefore, assuming that k s ∈ [k sm , k sM ] it can be consider as an affine model and can be rewritten as a polytopic model.

Remark 1: Note that in the following section, c 0 will be considered as a varying parameter.

The measurement y(t) used for control design is the suspension deflection speed:

y(t) = żus (t) -żs (t).
The controlled outputs are z s (t) and u(t) and are sufficient to control the sprung masse acceleration zs (t), the suspension deflection z def (t) and the vertical displacement of the unsprung mass z us (t).

The LP V /H ∞ suspension controller

The performance specifications are as follows. For a fixed value of c 0 (t) and for all admissible values of the varying stiffness coefficient k s (.), the controller must improve the passenger comfort without deteriorating the vehicle road holding and the suspension deflection of the passive suspension for the same frozen damping value c 0 .

As in the usual H ∞ framework, the performance objectives are achieved via minimizing weighted transfer function norms. The considered LP V /H ∞ control scheme is given in Fig. 8.

The LP V /H ∞ suboptimal control problem consists in designing a stabilizing controller K(s, k s ) depending on the varying parameter k s (.) that internally stabilizes the closed-loop system and minimizes a constant γ, s.t. the closed-loop transfer H ∞ -norm between the exogenous inputs w(t) to the exogenous outputs e(t) is less or equal to γ (i.e. e(t)/w(t) ∞ ≤ γ) for any The design method LP V /H ∞ developed in [START_REF] Zin | An LPV/H ∞ active suspension control for global chassis technology: Design and performance analysis[END_REF] provides the results shown in figure 9 for c 0 = 700 N/m/s where the frequency-domain behavior of the closedloop (i.e. active) system is compared to the open-loop (i.e. passive) one for the non linear model. Non linear frequency-domain analysis of the active suspension using sinusoïdal input z0(t) and with c0(t) = 700 N/m/s.

k s ∈ [k sm , k sM ].

Remark 2:

The "frequency" response of the nonlinear dynamical model is computed as follows: first, a finite set of pure sinusoïdal signals are chosen as being the road profile. The frequency of the sinusoids varies from 0.1 to 50.0 Hz with a constant amplitude of 0.015 m. Then, at each frequency, zs (t), z s (t), z us (t) and z def (t) are recorded. Finally, the frequency response of the nonlinear system is then computed as the input and output DF T magnitudes ratio.

As seen above, the driving comfort is improved between 0.5 and 5 Hz (see zs (s)/z 0 (s) and z s (s)/z 0 (s) in figure 9) without deteriorating neither the road holding nor the suspension deflection compared to the passive case (see the transfer z us (s)/z 0 (s) and z def (s)/z 0 (s)). However, it is well known that the drawback of road vibration isolation in suspension control leads to an amplification of the induced-oscillations due to load transfers for instance during a braking phase. In the next section, a feed-forward strategy is used to reduce these oscillations and to improve vehicle safety.

About performance adaptation

In the proposed design, c 0 (t) has been chosen constant in the design step, since the aim is to design a controller which improves the comfort while keeping the road holding as good as the passive suspension. The adaptation of the suspension characteristics is provided only by changing the value of c 0 .

Indeed, for different frozen values of the parameter c 0 (t), the active suspension of the closed-loop model of Fig. 7 presents better performances than the passive model for the same frozen damping values (i.e. c s (.) = c 0 (t)).

A comparison between the closed-loop frequency responses of the nonlinear passive and active systems is shown in Fig. 10 for c(t) = {700, 2200, 3700} N/m/s and when k s (.) varies. The dash lines correspond to the passive model, the solid lines to the closed-loop active LP V model. These values represent a hard, mid and soft suspension requirements, i.e. a road holding, deflection and confort objective. For each value of c 0 (t), the frequency responses of zs (t), z s (t), z us (t), z def (t), from z 0 (t), are plotted in order to evaluate the ability of the LP V controller to improve the ride comfort, the road holding and the suspension stroke w.r.t the passive LP V model.

It can be seen in Fig. 10 that for the same frozen value of the damping coefficient c 0 (t), the active suspension presents a better ride confort without deteriorating the road holding and the suspension deflection of the passive suspension. 0n the other hand, when the magnitude c 0 (t) decreases, the passenger comfort is increased, while the road holding is deteriorated (according to the trade-off between comfort/road holding [START_REF] Hedrick | Invariant properties of automotive suspensions[END_REF]). Note also, that the peak of the transfer from the road profile z 0 (t) to the actuator forces u ca (t) does not increase when c 0 (t) increases, which of course is not the case of the LP V control signal u(t) (see Fig. 11). This emphasizes the great intrinsic interest of using c 0 (t) as a varying coefficient in the controller. While allowing performance adaptation, neither the control structure does change nor the actuator frequency gain peak is increased. Active control input for c(t) = {700, 2200, 3700} N/m/s.

Remark 3 With reference to Fig. 4, a high value of the passive damping c 0 (.) is obtained for low values (in magnitude) of the suspension velocity deflection. In fact, the damper characteristic of Fig. 4 is a typical static map of a "sporty" suspension setup: in low suspension deflection velocity (that is in general the framework of a suspension) the road holding is favoured. The trade-off between safety and comfort does not depend on the suspension deflection velocity, but on the stability of the vehicle. With the simple strategy presented in Fig. 7 the safety/comfort ability of a suspension can be changed online, depending on the stability of the vehicle.

A feed-forward strategy for global chassis control through anti-roll distribution

The need for feedforward control laws in vehicle dynamics control is today well known [START_REF] Inagaki | Development of feedforward control algorithms for active suspension[END_REF]. In this section a new strategy is proposed through an anti-roll distribution of the suspension forces in order to improve the comfort/safety trade-off.

First, the objective is to keep the sprung mass horizontal whatever the longitudinal and lateral vehicle accelerations are. It corresponds to require zs (t) = φ(t) = θ(t) = 0, in the vertical dynamical equation system [START_REF] Gáspár | A fault-tolerant rollover prevention system based on LPV method[END_REF]. If in (4) the forces F sij (.) are replaced by the feedforward signals to get u f sij (.), they should satisfy:

(13)                    0 = u f sf r (.) + u f sf l (.) + u f srr (.) + u f srl (.), 0 = -u f sf r (.) + u f sf l (.) t f + -u f srr (.) + u f srl (.) t r + m H a y (t), 0 = u f srr (.) + u f srl (.) l r - u f sf r (.) + u f sf l (.) l f -m H a x (t).
The above system includes 3 equations for 4 unknown variables u f sij (.). In order to ensure the existence of a unique solution, an equation dedicated to the anti-roll distribution is added. This distribution is represented by the scalar η(.) which corresponds to the anti-roll torque provided by the rear axle compared to the front axle, in order to compensate the induced-roll torque by m H a y (t). This is formulated below by the equation :

(14) 1 -η(.) u f sf l (.) -u f sf r (.) t f = η(.) u f srl (.) -u f srr (.) t r .
Remark 4

• Note that, in real cars, this distribution is usually imposed by the anti-roll bars and by the suspensions themselves.

• With η(.) = 0, the anti-roll behavior is obtained by the rear axle only (both forces u f sf l (.) and u f sf r (.) equals); with η(.) = 1 it is achieved by the front axle.

From equations ( 13) and ( 14) the four feedforward forces become:

               u f sf r (.) = -1 2 H m 1 l f +lr a x (t) -η(.) t f a y (t) , u f sf l (.) = -1 2 H m 1 l f +lr a x (t) + η(.) t f a y (t) , u f srr (.) = -1 2 H m -1 l f +lr a x (t) -1-η(.) tr a y (t) , u f srl (.) = -1 2 H m -1 l f +lr a x (t) + 1-η(.)
tr a y (t) .

Since we wish to control η(t) and c 0 , we need to characterize the under and oversteering, hence to choose some "reference" vehicle. Here a "neutral" car is selected as a reference car.

The reference car model

The chosen reference car model is the well known bicycle model, i.e., the horizontal vehicle model [START_REF] Apkarian | Self scheduled H ∞ control of linear parameter-varying systems: a design method[END_REF], linearized around null steering wheel angle, and null tire sideslip angles. It is assumed that the left and right parts of the vehicle are symmetric. The model then satisfies the following equations: [START_REF] Oustaloup | The CRONE suspension[END_REF] m a y (t) = m vy (t) + v x (t) ψ(t) = 2 F ytf (.) δ f (t) + 2 F ytr (.)

I z ψ(t) = 2 F ytf (.) δ f (t) l f -2 F ytr (.) l r .
where F ytf (.) and F ytr (.) are the front and rear lateral tire forces respectively. These are linear functions of the sideslip angles α i , i.e. F yti (.) = C αi α i (t) with

   α f (t) = arctan vy(t)+l f ψ(t) vx(t) -δ f (t), α r (t) = arctan vy(t)-lr ψ(t) vx (t)
.

where C αi are constant coefficients, obtained from the non linear tire characteristics (here we have chosen C α f = -8.4 × 10 4 N/rad and C αr = -6.4 × 10 4 N/rad). This model is an LPV model according to the longitudinal speed v x (t) but does not include any limitation of the attainable speed. Therefore the steering wheel angle is saturated in order to avoid too large lateral accelerations. The saturation bounds are calculated as the static gain of the transfer function from δ f (t) to a y (t), paremetrized by v x (t), from equations [START_REF] Oustaloup | The CRONE suspension[END_REF]. This concept is represented in figure 12 where the vehicle speed is 90 km/h, and where a bend is realized on a dry road considering a saturation of a y (t) at 10 m/s 2 . The simulation results are represented in figure 12 and compared with the neutral car (η = 0.5). We can observe that the neutral reference car has no overshoot in ψ(t) and gives a linear characteristic for β(t)/a y (t).

Yaw control through anti-roll distribution

In this section, the previous reference model is used to derive a yaw control strategy by controlling the distribution η(.), as shown in figure 13. The control objective is to reduce the error between the reference yaw rate and the measured one (e ψ(t) = ψ(t) -ψref (t)). The control input is the variation ∆ η (.) of the parameter η around its nominal value η = 0.5. As ∆ η (.) is between -0.5 and 0.5, a very simple PI controller with anti-windup is here chosen.

Figure 13

Yaw rate PI control strategy with anti-windup through anti-roll distribution.

To evaluate the performance of the global control scheme, the same manoeuver as in figure 12 is considered. As seen in figure 14, the overshoot is reduced thanks to the anti-roll distribution, and the response of the closed-loop system is faster than the one of the reference car (for which η = 0.5). Of course there still exists some tracking error in the transient behavior but the use of anti-roll distribution greatly improves the yaw rate behavior (even if in this case the controlled vehicle is not completely neutral).

Towards a global chassis control strategy

In section 3, a LP V /H ∞ controller has been designed for suspension control. In figure 7, the given strategy includes the variation of the suspension stiffness coefficient k s and also, the parameter c 0 (t), which is considered constant during

Figure 14

Simulation results of the yaw control strategy using anti-roll distribution the design. On the other hand, the results in the previous section have shown that ∆ η (.) represents the degree of understeering/oversteering of the vehicle.

In order to increase the vehicle stability, the suspension will be made safety or comfort oriented according to the driving situations. As shown previously, the variation of the parameter c 0 (t) modifies the vertical dynamical performances, then it will be changed on-line, from c 0 (.) = 700 N/m/s (safety oriented), to c 0 (.) = 3700 N/m/s (comfort oriented). A very simple proportional controller is here proposed for validation.

For illustration, the same manoeuver as in figure 14 is simulated, for a random road profile. As seen in figure 15 , controlling c 0 greatly reduces the oscillations of the yaw rate, compared with the case where η only is controlled. As seen in the right-lower figure, c 0 is varying according to η.

Figure 16 shows the time responses of the vertical dynamics. As expected, controlling c 0 deteriorates the vehicle comfort but in favour of the vehicle safety.

Conclusions

In this paper a new global chassis control strategy is developed using the suspension actuators only. First, a robust H ∞ controller of the suspension is designed. Its inherent structure allows for performance adaptation through the variation of a coefficient c 0 . Then, to improve the vehicle safety in the presence of longitu- dinal and lateral disturbances, a feedforward strategy is proposed, by controlling the anti-roll distribution of the suspension forces. The control scheme consists in controlling the yaw rate through a model reference control. The global chassis control then relies on the adaptation of the suspension performance (comfort/safety) according to the anti-roll distribution. This strategy needs sensors for suspension control (4 deflection sensors and the initial damping forces), acceleration sensors (longitudinal, lateral), a yaw sensor and steering wheel sensors. All these sensors are currently used now in automotive industry. Some improvements may rely on the reference car model, which needs to adapt to road conditions (dry, wet, icy ...). Some issues may concern the estimation of the road adhesion to be used in this global chassis control (see for instance [START_REF] Stephant | Avaluation of sliding mode observer for vehicle sideslip angle[END_REF]). Global chassis control: analysis of the vertical dynamics.

Figure 1

 1 Figure 1Horizontal car model

Figure 2

 2 Figure 2Scheme of the horizontal dynamical model

Figure 3

 3 Figure 3Vertical car model

Figure 4

 4 Figure 4Non linear characteristics of the suspension forces.

Figure 5

 5 Figure 5Scheme of the vertical model.

Figure 6

 6 Figure 6Quarter-car suspension models (passive and active).

Figure 7

 7 Figure 7Scheme for the LPV design

Figure 8

 8 Figure 8LP V /H∞ control design interconnection.

  Figure 9Non linear frequency-domain analysis of the active suspension using sinusoïdal input z0(t) and with c0(t) = 700 N/m/s.

  Figure 10LP V /H∞ synthesis results in frequency domain for the closed loop model (with actuator) (solid line) and the passive model (dash line) for c(t) = cs(.) = {700, 2200, 3700} N/m/s and when ks(.) varies.

  Figure 11Active control input for c(t) = {700, 2200, 3700} N/m/s.

  Figure12model analysis.

Figure 15

 15 Figure 15Global chassis control: analysis of the lateral dynamics.

  Figure 16Global chassis control: analysis of the vertical dynamics.


Luc Hillairet
email: luc.hillairet@math.univ-nantes.fr

FINITE DEFICIENCY INDICES AND UNIFORM REMAINDER IN WEYL'S LAW

teaching and research institutions in France or abroad, or from public or private research centers.

Introduction

Von Neumann theory classifies all self-adjoint extensions of a given symmetric operator A in terms of the so-called deficiency indices. Since the choice of the self-adjoint condition reflects the physics that is underlying the problem, it is natural to ask how the spectra of two different self-adjoint extensions A 0 and A 1 can differ.

In some cases, the deficiency indices are finite. This happens for instance in the following interesting settings from mathematical physics -Quantum graphs (see [START_REF] Endres | The trace formula for quantum graphs with general self adjoint boundary conditions[END_REF][START_REF] Harrison | Vacuum energy, spectral determinant and heat kernel asymptotics of graph Laplacians with general vertex matching conditions[END_REF] and section 4.1 below), -Pseudo-Laplacians, Šeba billiards (see [START_REF] De Verdière | Pseudo-laplaciens I[END_REF][START_REF] Keating | Localised eigenfunctions in Šeba billiards[END_REF], and section 4.2 below), -Manifolds with conical singularities (see [START_REF] Kokotov | Compact polyhedral surfaces of an arbitrary genus and determinants of Laplacians[END_REF][START_REF] Kirsten | Exotic expansions and pathological properties of ζ-functions on conic manifolds[END_REF]), -Hybrid manifolds (see [START_REF] Pankrashkin | Resolvent expansions on hybrid manifolds preprint[END_REF]).

In any of these settings we prove the following theorem (see section 2 and the applications for a more precise version) Theorem 1. Let A be the symmetric operator associated with one of the preceding settings. There exists a constant C such that for any self-adjoint extensions A 0 and A 1 of A we have

∀E, |N 1 (E) -N 0 (E)| ≤ C
where N i denotes the spectral counting function of A i .

This fact actually derives from [START_REF] Birman | Spectral Theory of Self-Adjoint Operators in Hilbert Space Mathematics and its Applications (Soviet Series)[END_REF] ch. 9 sec. 3 1 . Our proof is slightly different and based on the min-max principle but the underlying ideas are similar.

Motivation for this result came principally from [START_REF] Keating | Localised eigenfunctions in Šeba billiards[END_REF] and the Šeba billiard setting. In this case we can take A 0 to be the standard Dirichlet Laplace operator in R, and this theorem proves that the remainder in Weyl's law is, up to a O(1) term, uniform with respect to the location of the Delta potential. In the case of one delta potential the uniform bound can also be derived from the fact that the spectra of the pseudo-laplacian and the usual laplacian are interlaced (see [START_REF] De Verdière | Pseudo-laplaciens I[END_REF] for instance).

In constrast with [START_REF] Endres | The trace formula for quantum graphs with general self adjoint boundary conditions[END_REF][START_REF] Harrison | Vacuum energy, spectral determinant and heat kernel asymptotics of graph Laplacians with general vertex matching conditions[END_REF][START_REF] Kirsten | Exotic expansions and pathological properties of ζ-functions on conic manifolds[END_REF][START_REF] Pankrashkin | Resolvent expansions on hybrid manifolds preprint[END_REF] we consider a rather crude spectral invariant. Moreover, our result relies on the min-max principle only which is less sophisticated than the analysis performed in the former references. It should be noted, however, that our result is not a straightforward byproduct of these results and should more likely be considered as a first step. From our perspective, it 1 We are grateful to Alexander Pushnitski for indicating to us this reference. is quite interesting to have a general method allowing to get quite good hold on the spectral counting function before moving on to more complete spectral invariants such as heat kernel or resolvent estimates.

Setting and Notations

We begin by recalling some basic facts from spectral theory of self-adjoint operators as well as Von Neumann theory of self-adjoint extensions of a symmetric operator. We will use [START_REF] Reed | Methods of modern mathematical physics. I. Functional Analysis[END_REF][START_REF] Reed | Methods of modern mathematical physics. II. Fourier analysis, self-adjointness[END_REF] and [START_REF] Courant | Methods of mathematical physics[END_REF] as references.

1.1. Basic Spectral Theory. We consider a Hilbert space H with scalar product •, • and associated norm • .

On H we consider a symmetric operator A with domain dom(A) and its adjoint A * . The graph norm is defined on dom(A) by

u 2 A = u 2 + Au 2 . An operator is self-adjoint if A = A * . It has compact resolvent if the injection from dom(A) into H is compact.
The spectrum of a self-adjoint operator with compact resolvent consists in eigenvalues of finite multiplicities, that form a discrete set in R. There exists an orthonormal basis consisting of eigenvectors.

If there exists C ∈ R such that ∀u ∈ dom(A), Au ≥ C u , the operator is called semibounded.

For a semibounded self-adjoint operator with compact resolvent, the spectrum can be ordered into a non-decreasing sequence (λ n ) n∈N .

The spectral counting function is then defined by

N (E) = Card{λ n ≤ E},
and by Courant-Hilbert min-max principle (see [START_REF] Courant | Methods of mathematical physics[END_REF]) we have

(1.1) λ n = min max Au, u u 2 , u ∈ F \{0} , F ⊂ H, F vector space s.t. dim F = n .

Von Neumann Theory.

This section summarizes section X.1 of [START_REF] Reed | Methods of modern mathematical physics. II. Fourier analysis, self-adjointness[END_REF]. We define

K ± = ker(A * ∓ id) d ± = dim(K ± ).
We recall that we have the following decomposition of dom(A * ) (see Lemma in section X.1 of [START_REF] Reed | Methods of modern mathematical physics. II. Fourier analysis, self-adjointness[END_REF]).

dom(A * ) = dom( Ā) ⊕ K + ⊕ K -,
and that the decomposition is orthogonal with respect to the scalar product

•, • A * defined on dom(A * ) by ∀u, v ∈ dom(A * ), u, v A * := u, v + A * u, A * v .
We denote by π ± the orthogonal projection from dom(A * ) onto K ± and by π 0 the orthogonal projection onto dom( Ā) (so that π 0 = idπ +π -).

The following theorem (theorem X.2 and corollary of [START_REF] Reed | Methods of modern mathematical physics. II. Fourier analysis, self-adjointness[END_REF]) provides a parameterization of all self-adjoint extensions of A Theorem 2. A admits self-adjoint extensions if and only if

d + = d -.
For any self-adjoint extension A sa of A, there exists a unique isometry

U from K + onto K -such that ∀u ∈ dom(A sa ), U π + (u) = π -(u).

The theorem

Using the notations of the preceding section we have Theorem 3. Let A be a symmetric operator with equal finite deficiency indices :

d + = d -= d < ∞.
Suppose that there exists A 0 a self-adjoint extension of A such that

(i) A 0 has compact resolvent, (ii) A 0 is semibounded Then (1

) Any other self-adjoint extension also has compact resolvent and is semibounded.

(2) There exist E 0 such that, for any other self-adjoint extension A 1 the following holds

∀E ∈ R, |N 1 (E) -N 0 (E)| ≤ d
where N i (E) denotes the spectral counting function of A i .

As indicated in the introduction this result is actually already proved in [START_REF] Birman | Spectral Theory of Self-Adjoint Operators in Hilbert Space Mathematics and its Applications (Soviet Series)[END_REF] using that the difference of the resolvents (A 0z) -1 -(A 1z) -1 is finite rank for some z.

Proofs

3.1.

A 1 has compact resolvent. We prove this fact by proving the stronger lemma. Lemma 3.1. If A 0 has compact resolvent then the injection from dom(A * ) (equipped with

• A * ) into H is compact. Proof. Since A 0 has compact resolvent, the injection from (dom(A 0 ), • A 0 ) into H is compact. But dom( Ā) is closed in dom(A * ) for • A * . Since • A 0 coincide with • A * on A 0 , the injection from (dom( Ā), • A * ) is also compact. Let (u n ) n∈N be a sequence in dom(A * ) that is • A * -bounded.
We can extract a subsequence from π 0 u n that converges in H since the injection from (dom( Ā), • A * ) into H is compact. On K ± , • A * is equivalent to • and, since d ± are finite we can also extract convergent subsequences. This proves the Lemma.

The Lemma says that A * has compact resolvent and so does A 1 since A * extends A 1 and the latter is closed.

3.2.

A 1 is semibounded. (See also [START_REF] Akhiezer | Theory of linear operators in Hilbert space[END_REF] sec. 85) Let (λ k (A 0 )) k∈N denote the (ordered spectrum) of A 0 and consider n such that λ n (A 0 ) ≥ 0.

Consider F ⊂ dom(A 1 )of dimension n + d. Denote by F 1 = F ∩ dom( Ā) = F ∩ ker(idπ 0 ). Theorem 2 implies that ker(idπ 0 ) |dom(A 1 ) is of dimension d, and thus, since F is of dimension n + d, F 1 is of dimension at least n.

Moreover F 1 ⊂ dom( Ā) ⊂ dom(A 1,0 ), thus, for all u ∈ F 1 we have

A 1 u, u = A 0 u, u . Since F 1 ⊂ F it follows that max u∈F,u =0 A 1 u, u u 2 ≥ max u∈F 1 ,u =0 A 0 u, u u 2 .
Since dim F 1 ≥ n and λ n ≥ 0, it follows from the min-max principle that the right-hand side is non-negative.

We thus obtain, that for all F ⊂ dom(A 1 ) of dimension n + d we have

max u∈F,u =0 A 1 u, u u 2 ≥ 0.
This implies that A 1 has at most n+d-1 negative eigenvalues. (otherwise the subspace generated by n + d negative eigenvalues would contradict the preceding bound).

3.3.

Comparing N 0 and N 1 . According to the previous section we have that A 1 also is semibounded with compact resolvent so that we can denote by (λ k (A 1 )) k∈N its ordered spectrum. We also denote by V i n the vector space generated by the first n eigenvectors of A i . For any n, set F = V 1 n+d ∩ dom( Ā). Making the same reasoning as previously we find that F is of dimension at least n and max

u∈V n+d ,u =0 A 1 u, u u 2 ≥ max u∈F,u =0 A 0 u, u u 2 .
By definition of V n+d the left-hand-side is λ n+d (A 1 ) and, using the min-max principle, the righthand side is bounded below by λ n (A 0 ). Thus, we obtain,

∀n, λ n+d (A 1 ) ≥ λ n (A 0 ).
Observe that N 0 (E) is characterized by

λ N 0 (E) (A 0 ) ≤ E < λ N 0 (E)+1 (A 0 )
Using this and the preceding inequality we find that, for all E we have

E ≤ λ N 0 (E)+d+1 (A 1 )
An thus, for all E we have

N 1 (E) ≤ N 0 (E) + d.
Since A 0 and A 1 now play symmetric roles we also have ∀n, λ n+d (A 0 ) ≥ λ n (A 1 ).

and thus N 0 (E) ≤ N 1 (E) + d. The final claim of the theorem follows.

4. Applications 4.1. Quantum graphs. Quantum graphs are now well-studied objects from mathematical physics (see [START_REF] Kuchment | Quantum graphs: an introduction and a brief survey[END_REF] for an introduction). A very rough way of defining a (finite) quantum graph is the following.

Pick K positive real numbers (the lengths), set

H = K i = 1 L 2 (0, L i ) and D = K i = 1 C ∞ 0 (0, L i ) and define A on D by A(u 1 ⊕ u 2 ⊕ • • • ⊕ u K ) = -(u ′′ 1 ⊕ u ′′ 2 ⊕ • • • ⊕ u ′′ K )
. This operator is symmetric. Any self-adjoint extension of A is called a quantum graph.

Remark 4.1. Usually quantum graphs are constructed starting from a combinatorial graph given by its edges and vertices. This combinatoric data is actually hidden in the choice of the self-adjoint condition.

One basic question is to understand to which extent the knowledge of the spectrum determines the quantum graph (i.e. the lengths and the boundary condition).

It is known that there are isospectral quantum graphs [START_REF] Band | The isospectral fruits of representation theory: quantum graphs and drums[END_REF] and the following theorem says that, as far as counting function is concerned it is quite difficult to determine the self-adjoint condition. Theorem 4. For any quantum graph with K edges the following bound holds :

N (E) - L π E 1 2 + ≤ 3K
where L := K i=1 L i and E + := max(E, 0).

Proof. Fix K and the choice of the lengths. We choose one particular self-adjoint extension A D that consists in decoupling all the edges and putting Dirichlet boundary condition on each end of each edge. The spectrum is then easily computed and we have

spec(A D ) = K i=1 k 2 π 2 L 2 i , k ∈ N .
In particular we have (denoting by N D the counting function of the Dirichlet extension)

N D (E) = K i=1 L i π E 1 2 + ,
where [•] denotes the integer part. In particular, we have

N D (E) - L π E 1 2 + ≤ K.
We now compute the deficiency indices. A straightforward computation yields d + = d -= 2K. And thus, using the main theorem and triangular inegality we obtain that for any quantum graph

N D (E) - L π E 1 2 + ≤ 3K,
indepently of the choice of lengths.

4.2. Pseudo-Laplacian with Delta potentials. It is known that on Riemannian manifolds of dimension 2 or 3 it is possible to add so-called Delta potentials. From a spectral point of view this corresponds to choosing a finite set of points P and to consider the Riemannian Laplace operator defined on smooth functions with support in M \P. (see [START_REF] De Verdière | Pseudo-laplaciens I[END_REF] for instance) Remark 4.2. This construction is also possible starting from a bounded domain in R 2 with, say, Dirichlet boundary condition. We obtain the so-called Šeba billiards (See [START_REF] Keating | Localised eigenfunctions in Šeba billiards[END_REF] for instance)

There are several self-adjoint extensions and (a slight generalization of) Lemma of [START_REF] De Verdière | Pseudo-laplaciens I[END_REF] proves that, in this setting the deficiency indices are d := card P. Following Colin de Verdière we call any such self-adjoint extension a Pseudo-Laplacian with d Delta potentials Application of the theorem gives the following. It should be noted first that the bound depends only on the number of Delta potentials and not on their location, and second that the effect of adding Delta potentials is much smaller than the usual known remainder terms in Weyl's law for N 0 (E).

4.3.

Others. There are two other settings where finite deficiency indices occur that are worth mentioning. In both case one could apply the theorem to get a Weyl's asymptotic formula independent of the choice of the self-adjoint condition up to a O(1) term. These are (1) Manifolds with conical singularities. The common self-adjoint extension in use corresponds to Friedrichs extension and if one changes the self-adjoint extensions at the conical points (see [START_REF] Kokotov | Compact polyhedral surfaces of an arbitrary genus and determinants of Laplacians[END_REF][START_REF] Kirsten | Exotic expansions and pathological properties of ζ-functions on conic manifolds[END_REF]) then the counting function is affected only by some bounded correction.

Observe that the deficiency index associated with the Laplace operator on the cone of opening angle α is 2[ α 2π ] -1 (with [•] the integer part). (2) The so-called hybrid manifolds that are studied in [START_REF] Pankrashkin | Resolvent expansions on hybrid manifolds preprint[END_REF] which are obtained by, in some sense, grafting quantum graphs onto higher dimensional manifolds. Here again, one can compare the counting function of the chosen self-adjoint extension to the natural one which is obtained when all the parts are decoupled.

1

 1

Theorem 5 .

 5 Let M be a closed Riemannian manifold of dimension 2 or 3. Let N 0 be the counting function of the (standard) Laplace operator on M . For any pseudo-laplacian with d Delta potential the following holds |N (E) -N 0 (E)| ≤ d

Acknowledgments : We are grateful to Alexander Pushnitski, Alexey Kokotov and Jens Marklof for useful comments on the first version of this note that resulted, in particular, in great improvment of the bibliography.


Y Renier

Lal

ORBIT RECONSTRUCTION, CORRECTION AND MONITORING IN THE ATF2 EXTRACTION LINE*

The orbit in the ATF2 extraction line has to be accurately controlled to allow orbit and optics corrections to work well downstream. The Final Focus section contains points with large beta function values which amplify incoming beam jitter, and few correctors since the steering is performed using quadrupole movers, and so good orbit stability is required.

First experience monitoring the orbit and measuring transfer matrices in the extraction line during the ATF2 commissioning is described, along with relative orbit reconstruction, allowing comparisons with simulations to identify sources of variations.

INTRODUCTION

Layout of ATF-ATF2

ATF is a facility composed of a 1.3 GeV linac and a damping ring producing small emittance beams (Ex=1.2 nm.rad, Ey=12 pm.rad). This beam is then injected in ATF2 which is an Extraction (EXT) followed by a Final Focus (FF) beamline (Figure 1) [1,2].

Objectives

The first goals of ATF2 are to get reliably 37 nm vertical beam size at the interaction point (IP), with less than 10 nm position jitter. It will be the first experimental validation of the local chromaticity correction scheme, planned to be used in ILC and CLIC.

To reach these goals, a good orbit control is needed, especially in the FF sextupoles. Reconstruction of the orbit allows to monitor the parameters of the injected beam, which is the first step to implement feedback. As these studies were done at an early stage of the commissioning, it also helped to find and correct several hardware or software problems.

Principle of orbit reconstruction

Transfer matrices between the injection point and the BPMs give the influence of an input parameter's variation on the BPM readings. The reconstruction, solving equation 1, achieves the opposite : find the parameters of the beam at the injection point from the BPM readings.

(1)

TRANSFER MATRIX MEASUREMENTS

As the reconstruction relies heavily on transfer matrices, the model of the beamline was first tested, allowing at the same time a check of BPMs and corrector polarities and scales, by measuring R 12 and R 34 matrix elements between correctors and BPMs.

The measurements were made changing a corrector strength and averaging 10 readings for each setting, to obtain the coefficient between the angle introduced by the corrector and the variation of the position measured at each BPM, in order to compare with the model value. The results for the first horizontal and vertical corrector are shown in Figure 2 and3. The very good agreement (especially in the FF) for most of the BPMs gave us good confidence in the modeling of the beamline. The precision is not as good where there were too large amplitudes because of saturation in BPM readings.

RELATIVE ORBIT RECONSTRUCTION

Reconstruction of parameters at injection

The reconstruction was then tried on 500 BPM readings, changing the energy by varying the frequency of the ring cavity. The x, x', y, y', dE/E parameters at the injection point were fitted by least square minimisation solving Equation 1. Results for the horizontal plane are shown in Figure 4. One can clearly see that the steps in energy are well reconstructed. Also a correlation between x and dE/E and between x' and dE/E is visible, indicating some anomalous dispersion at that point.

In the vertical plane, natural fluctuations of the beam being 10x lower than in horizontal, reconstruction was not successful.

Scale factor determinations

The reconstructed orbit deduced from these parameters is then compared to the measurements to obtain scale factors in BPM readings by a linear fit (see Figure 5). Horizontal scale factors of all BPMs obtained in this way are shown in Figure 6. This estimation can be biased by an overall scale factor common to all BPMs, and should be calibrated, thanks to the energy fit or the transfer matrix measurements.

Most BPMs have scale factors between 0.75 and 1.5.

To evaluate the quality of the fit, the spread of residuals is computed for each BPM, with or without correcting for the scale factors, and compared to the initial spread of the BPM readings (see Figure7). For almost all BPMs, the spread of the residuals of the reconstruction is significantly smaller than the spread of the measurements.

It is especially true in the first part of the line where the dispersion function is large. That is a very good indication that the energy is well reconstructed.

One can also see a good reconstruction of x and x' at the end of the line where the the beta functions become very large, which is, this time, due to a good reconstruction of the positions at that location.

For the BPMs with non-unitary scale factors, the spread is further reduced when the reading is corrected by the found scale factor.

CONCLUSION AND PROSPECTS

It was shown that the modeling of ATF2 is reliable, allowing successful reconstruction of the horizontal positions, angles and of the energy at the injection point.

Thanks to this reconstruction, one can estimate parasitically, with just beam fluctuations, the scale factors of the BPMs.

Nevertheless, because the reconstruction suffers from the bad resolutions of the first BPMs, work is on-going on their electronics, trying to improve it. That, combined with a reconstruction optimized to exploit the sub-micron resolution of the last BPMs, should allow us to reconstruct as well the vertical plane.

Figure 1 :

 1 Figure 1: ATF and ATF2 layout.

Figure 2 :

 2 Figure 2: R 12 measured (solid line) and modeled (dashed line) for all BPMs.

____________________________________________*

  Work supported by Agence Nationale de la Recherche of the French Ministry of Research (Programme Blanc, Project ATF2-IN2P3-KEK, contract ANR-06-BLAN-0027)

Figure 3 :

 3 Figure 3: R 34 measured (solid line) and modeled (dashed line) for all BPMs.

Figure 4 :

 4 Figure 4: Evolution of x, x' and dE/E reconstructed at the injection point.

Figure 5 :

 5 Figure 5: Example of scale factor fit. Horizontal and vertical axes show the reconstructed and experimental readings, respectively.

Figure 6 :

 6 Figure 6: Horizontal scale factors of all BPMs, with estimated errors from the fit.

Figure 7 :

 7 Figure 7: Spread of initial measurements and of the residuals to the reconstruction, with and without correction for scale factor.