Roger Prud'homme

Kwassi Anani

EVAPORATION D'UNE PASTILLE EQUIVALENTE VAPORIZATION OF AN EQUIVALENT PASTILLE

Keywords: pastille, évaporation, oscillations harmoniques, facteur de réponse, perturbations de température pastille, evaporation, harmonic oscillations, response factor, temperature perturbation Liste des symboles

On examine la réponse dynamique d'une goutte en évaporation en forme de pastille sous l'effet de petites oscillations acoustiques. La pastille est alimentée de façon continue avec un coefficient d'échange thermique h et il est tenu compte des effets de la convection de la chaleur et de la conduction sur le facteur de réponse en fréquence du système et sur le champ thermique dans la gouttelette. L'effet du coefficient d'échange thermique h est étudié.

' , f f  f ˆ    G h H , k l L,  n m, M M  M  N p Pe q p, Q L Q 

Introduction

On examine la réponse dynamique d'une goutte en évaporation en forme de pastille (un petit cylindre de liquide appelé "pastille", Figure 1) sous l'effet de petites oscillations acoustiques. Dans ce type d'études inspirées par le modèle de Heidmann i pour les gouttes sphériques ii , la pastille est alimentée de façon continueadiabatique ou isotherme -et il est tenu compte des effets de la convection de la chaleur et de la conduction sur le facteur de réponse en fréquence du système et sur le champ thermique dans la gouttelette. Dans le présent travail l'alimentation de la pastille se fait avec un coefficient d'échange thermique et n'est plus limité aux cas limites adiabatique et isotherme.

La pastille équivalente a été introduite dans une modélisation simplifiée iii des instabilités HF de la combustion de gouttes. Il s'agit en effet de remplacer la résolution d'un problème instationnaire à symétrie sphérique de variables r et t par un problème 1D instationnaire à variable x et t.

L'avantage est que les équations aux petites perturbations de la pastille soumise à un champ acoustique sont linéaires à coefficients constants et présentent donc des solutions analytiques.

Comparée à la goutte sphérique de diamètre moyen s r 2 , une pastille équivalente a le même volume et s'évapore uniquement sur ses surfaces de base comme le montre la figure 1. ).

Pour la phase gazeuse, nous ferons les mêmes hypothèses que dans les travaux précédents : mélange constitué de la vapeur de la goutte F et du gaz brûlé A, sachant que la combustion complète a lieu loin de la goutte dans la chambre de combustion (indice c). L'état stabilisé de référence est celui où les grandeurs f ont la valeur locale f , alors que dans l'état perturbé de manière harmonique par une onde acoustique stationnaire,

f f f    , 1 , '     f f f f f , et     t i e f t f  x x , ' 
. On suppose la goutte positionnée en un ventre de pression (qui est aussi un noeud de vitesse) d'amplitude maximale c p ˆ. On montre que l'évolution est gouvernée, au niveau de la surface de la goutte par les deux équations :

    S c L c S T p a M Q p T b iu iu M , 1 ˆ           [1] Dans le système [1],   , , , , c S L p T Q M
sont respectivement le débit masse d'évaporation, le flux de chaleur pénétrant dans la goutte (n'ayant pas servi à l'évaporation), la température de surface, la pression de chambre et la chaleur latente d'évaporation. Les coefficients intervenant dans ces équations dépendent des valeurs des quantités thermodynamiques et chimiques de référence du mélange non perturbé dans la chambre et au niveau de la surface d'évaporation. On trouve leurs expressions dans les articles cités et dans l'annexe A4.3 de l'ouvrage de Prud'homme vi . Les équations du système [1] seront couplées à celle de la phase liquide lors de la résolution finale.

Equations du problème

La phase liquide est supposée incompressible et non dilatable. Ses coefficients thermodynamiques sont supposés constants.

Les notations sont les suivantes : -Masse moyenne de la goutte cylindrique -Temps de séjour moyen (ou durée de vie) du liquide

L S v v x M M   
 .

-Temps caractéristique de conduction thermique

L S T x   2 ~ (avec L L L L c k    ).
-Nombre de Péclet dans le liquide

v T L   Pe  (son inverse  est utilisé dans certains calculs)
Les équations de base de la phase liquide sont :

-Bilan de masse de la pastille :

M M dt dM    .
-

Bilan local d'énergie interne : 0 2 2          x T x T v t T L L  Les conditions limites s'écrivent : -En S x x  (surface libre)     L x l L S S l Q x T k S t T t x T S     , , . -En 0  x (au fond du cylindre où a lieu l'alimentation)   S l l T T H x T     0 . Pour une alimentation adiabatique 0  H et si l'alimentation est isotherme   H .

Analyse linéaire de la phase liquide

La fonction  

L Pe u G ,
On pose, comme il a été indiqué dans l'introduction :

f f f    et 1 , '     f f f f f . L'équation de l'énergie devient : 0 ' ' ' 2 2          x T x T v t T L L  et la condition au niveau de la surface libre s'écrit : L x l L c Q x T k r S     ' 2  avec S c r r 2  . En posant     t i e x f t x f  , ' 
, on obtient :

  0 ~2    T S v S x s x s i    , soit encore :     0 3 2    u Pe i x s Pe x s L S L S  avec T u   3  , 1   en présence de convection et O  
en négligeant la convection. Cette équation caractéristique a deux racines   s s et . On exprime alors la perturbation de température avec les quantités :

H x h x x x s n x s m S S S S       , , ,  :           h n e h m e h n e h m e T T n m n m S l           ˆ [2]
La condition à la surface libre s'écrit

    L S S L x l S L L Pe u G x T T k S x T T k S Q S , ˆ     et l'on a :   n m n m L e h n e h m e h n n e h m m Pe u G        1 1 , [3]

Résolution

Les coefficients m et n sont reliés aux quantités   , , u par l'enchainement suivant :

t it R n t it R m u R u u t L L L L L L                                          [4] On déduit : c S p Z iu iu b T 1 1 1 ˆ          
, avec la fonction de transfert complexe

Z :     L L u G B u G A iu iu Z Pe , 3 Pe , 3 1       [5]
Rappelons que le facteur de réponse N est défini à partir de :

           t V t V dV dt t V p dV dt t V p t V q N , 2 , , ' , ' , ' [6]
Où p' est la perturbation incidente (la perturbation de pression de chambre c p' ) et q' la réponse du système (ici la perturbation de débit de la goutte ' M  ) à cette perturbation (Critère de Rayleigh iv ).

Le facteur de réponse réduit est alors égal à la partie réelle de

Z :   Z N    [7]
En définissant la perturbation de température du liquide dépendant du temps par :

c réd l p b T 
, on obtient :

        h n e h m e h n e h m e G B B A T n m n m réd l             3 ˆ [8]
Et la perturbation thermique réelle à l'intérieur de la pastille s'écrit :

                                   h n e h m
         [11]
et donne une évaluation de la profondeur de pénétration de l'onde acoustique dans la goutte lorsque cette quantité est faible.

3 Ce développement limité n'a de valeur que pour les faibles profondeurs de pénétration   .

Quelques résultats

Sans convection thermique

Avec convection thermique Fond adiabatique Pour des valeurs plus petites de PeL, la température de la surface de la pastille pénètre mieux l'intérieur de la pastille, sans atteindre la situation de la figure 3a, vu que le coefficient d'échange h n'est pas nul cette fois.

0 , 0   h  0 , 1   h  Coefficient d'échange 0 , 0   h  0 , 0   h  Fond isotherme   h , 0  ∞   h , 1  ∞ Tableau

Facteur de réponse

Conclusion

Cette étude permet d'affiner des résultats relatifs aux effets de l'alimentation en liquide d'une goutte moyenne dans un moteur en introduisant un coefficient d'échange. En effet les études précédentes ne prenaient en considération que les cas limites de l'injection adiabatique et isotherme. On voit donc que l'extension de la configuration de Heidmann en tenant compte du champ thermique intérieur à la goutte moyenne supposée représenter la suspension en combustion dans un moteur à ergol liquide nécessite quelques précautions. La vérification expérimentale s'impose bien entendu. Mais l'introduction du coefficient d'échange introduit une nouvelle souplesse au modèle rendant possible une meilleure correspondance à la réalité. D'autre part, des aspects physiques portant sur cette goutte moyenne peuvent aider à mieux comprendre le comportement de l'écoulement diphasique lui-même. Il en est ainsi pour le temps caractéristique v  , qui est à la fois le temps de séjour du liquide entre le point d'alimentation et la surface d'évaporation, et le temps de vie moyen des gouttelettes de la suspension en combustion. L'intervention du nombre de Péclet de la phase liquide a certainement un sens pour l'écoulement de la suspension de même que les paramètres d'alimentation de la goutte moyenne.

Il reste à adapter ce modèle au cas de l'interface d'évaporation sphérique.

Les calculs effectués dans cet article et les tracés de courbes résultantes ont été réalisés avec les programmes Maple de tableaux 2 et 3.

Le modèle de la goutte moyenne a été utilisé par ailleurs pour une application de la thermodynamique irréversible étendue par Prud'homme vi (pp.282-288).

Il a servi de référence pour le perfectionnement du code CEDRE de l'ONERA ii . Des travaux récents ont été publiés sur la goutte sphérique moyenne par Anani et al. vii . D'autres travaux plus généraux complèteront utilement la documentation sur le sujet des instabilités en présence de gouttes en évaporation et en combustion viii,ix,x,xi .

2 Figure 1 . 2 .Figure 2 .

 2122 Figure 1. Configuration géométrique d'une pastille équivalente à une sphère de diamètre

  La courbe des variations de température coupe l'axe des  pour :

1 .

 1 Les six types de configuration en présence d'un coefficient d'échange au fond d'alimentation. Le tableau1 donne les différentes configurations rendues possibles avec notre modélisation. par Anani v . Nous allons examiner les influences de coefficients d'échange de valeurs finies non nulles sur les courbes du facteur de réponse réduit  N et sur les profils de perturbations thermiques réduits réd lT ' à l'intérieur de la pastille.

  4.1. Perturbations thermiquesLes deux paramètres ayant une influence sur la conduction de la chaleur sont le nombre de Péclet PeL entre le haut et le bas de la pastille, et le coefficient d'échange h au bas de la pastille au niveau où pénètre le liquide d'alimentation qui est à la température S T . On vérifie bien que PeL contrôle la pénétration de la perturbation thermique à partir de la surface, alors que h a une action décisive au fond.

Figure 3 .Figure 4 .

 34 Figure 3. Influence du coefficient d'échange (h=0 ; 10 ; 1000) pour des valeurs identiques des autres paramètres

  Si l'on prend en compte la convection interne au liquide (ce que nous avons fait en section 4.1 où nous avons pris (fournie par la valeur uC de l'intersection avec l'axe N=0) décroit d'abord en fonction de  puis croit jusqu'à une valeur limite pour  infini. Cette valeur limite correspond à la conduction thermique infiniment rapide dans le liquide. Pour h=0,1, la courbe 1000   disparaît du demi-plan supérieur pour passer dans les N<0 (stabilité). Puis seules les courbes 1 positives (donc des zones d'instabilité) . Enfin, pour les valeurs de h supérieures (h =10 et h =1000, ce qui signifie une alimentation isotherme), on constate la stabilisation complète ; en particulier la stabilité s'est installée pour les valeurs de 5   .

Figure 5 .

 5 Figure 5. Courbes donnant le facteur de réponse réduit

Figure 6 .

 6 Figure 6. Différences entre les champs de température du liquide pour les deux conditions en 0 x  :

  Nous porterons notre attention sur les conditions d'alimentation de la goutte en forme de pastille. Le débit masse injecté est toujours le débit moyen d'évaporation M  , mais les conditions thermiques peuvent être différentes. Nous les caractérisons par un simple coefficient d'échange

								H . Jusqu'à présent seuls ont été étudiées
	l'injection adiabatique (	H		0	) et l'injection isotherme (	H		

Sorbonne Université , CNRS, Institut Jean Le Rond d'Alembert, UMR 7190, F-75005 Paris, France -e-mail : roger.prud_homme@upmc.fr

Université de Lomé, Faculty of Sciences, Department of Mathematics, 02 BP 1515, Lomé, Togo -e-mail : kanani@univ-lome.tg


Paolo Viappiani
email: paolo.viappiani@lip6.fr

Positional Scoring Rules with Uncertain Weights

Keywords: Scoring rules, rank aggregation, possible winners, minimax regret, stochastic dominance, convex sequences

Positional scoring rules are frequently used for rank aggregation (for example in social choice and in sports). These rules are highly sensitive to the weights associated to positions: depending on the weights, a different winner may be selected. In this paper we explicitly consider the role of weight uncertainty in both the case of monotone decreasing weights and of convex decreasing weights. First we discuss the problem of finding possible winners (candidates that may win for a feasible instantiation of the weights) based on previous works that established a connection with the notion of stochastic dominance. Second, we adopt decision-theoretic methods (minimax regret, maximum advantage, expected value) to pick a winner based on the weight uncertainty and we provide a characterization of these methods. Finally, we show some applications of our methodology in real datasets.

Introduction

In many contexts it is necessary to aggregate several rankings and either pick a winner or determine an output ranking. For example, rank aggregation emerges in recommender systems and social choice (preference aggregation), in information retrieval (aggregation of the output of several search engines), in sports (aggregation of the performance in several races into a single score). Positional scoring rules are frequently used due to their simplicity; in addition they satisfy a number of interesting properties (one of the most prominent results is that of Young [START_REF] Young | Social choice scoring functions[END_REF]: a symmetric social choice function is continuous and consistent if and only if it is a scoring rule). Among positional scoring rules, Borda is a well-known method whose properties have been studied in depth [START_REF] Young | An axiomatization of borda's rule[END_REF][START_REF] Fishburn | Borda's rule, positional voting, and Condorcet's simple majority principle[END_REF].

A scoring rule assigns a scores to each candidate based on the rank obtained in each ranking. The output of a scoring rule crucially depends on the weights (attached to ranks); in general several different winners are possible with different weights. In this article we consider methods to generate a winner using positional scoring rules under weight uncertainty. We assume non increasing weights; we also discuss the case of scoring rules where weights constitute a convex sequence.

Work supported by the ANR project Cocorico-CoDec (ANR-14-CE24-0007-01).

We assume n agents express preferences in the form of rankings involving a set of m candidates A = {a, b, c, . . .}; rankings are assumed to be linear orders (complete, transitive, asymmetric and irreflexive binary relation). Positional scoring rules discriminate between candidates by assigning a fixed score to each rank. Scoring rules assign a score to each alternative based on its rank distribution; let v x j be the number of times alternative x was ranked in the j-th position. Note that m j=1 v x j = n for each x ∈ A and x∈A v x j = n for each j = 1, . . . , m. A scoring rule specifies the vector of weights w 1 , . . . , w m (also called scoring vector) to be assigned to each position. The score obtained by a candidate according to weight vector w = (w 1 , . . . , w m ) is s w (x) = m j=1 w j v x j . These total scores can be used to pick a winner or to rank the alternatives. For example plurality is obtained by setting w 1 = 1 and w j = 0 for all j ∈ {2, . . . , m}.

By choosing a particular w, it is possible to specify some preferences on which kind of aggregation is desired, by giving more or less weight to the first positions compared to the positions that came afterwards in the ranking. First of all, we assume that not all weights are null, otherwise the alternatives are not discriminated (degenerated scoring rule). A natural hypothesis1 , that we adopt here, is to require that the sequence of weights is non-increasing: w i ≥ w i+1 for all i ∈ {1, . . . , m -1}; intuitively, in a ranking an alternative is (weakly) preferred to the alternatives that comes afterwards in the ranking.

A scoring rule is invariant to affine positive transformation of the scoring vector, that means that ranking the w = αw + β with α > 0 and arbitrary β gives the same output. Therefore, with no loss of generality, we let w 1 = 1 and w m = 0 (therefore we have m -2 degrees of freedom). Given this assumption the Borda rule is given by setting w j = m-j m-1 . Moreover, it is often (but not always) assumed in practice that the positional weights constitute a convex sequence, meaning that the difference between the first and the second weight is not less than the difference between the second and of the third, and so on. In such a case the weights need to satisfy the following constraint, for each i between 1 and m -2:

w i -w i+1 ≥ w i+1 -w i+2 ⇐⇒ w i -2w i+1 + w i+2 ≥ 0. (1)
Note that Borda and plurality are convex; furthermore convexity is often satisfied by the weights used when combining ranks in sports, races and other situations (e.g. formula one world championship, alpine skiing world cup).

We argue that setting a precise vector of weights for a scoring rule can be seen arbitrary; indeed the decision of which weights to use is critical since different weights often lead to different winners. Therefore there is interest to reason not just with a fixed vector w but about possible sets of parameters.

Notation We introduce some of the notation that we will use in the paper. We use [[m]] to denote the set {1, . . . , m}. Given two vectors v 1 and v 2 we write

v 1 v 2 iff v 1 j ≥ v 2 j for all components j. If the inequalities are strict, v 1 j > v 2 j
for all components j, then we write v 1 v2 . We will use W D to denote the set of scoring vectors with non-increasing weights and W C to denote the set of non-increasing (i.e. weakly decreasing) scoring vectors whose weights constitute a convex sequence 2 (with our boundary assumptions, w 1 = 1 and w m = 0):

W D = (w 1 , . . . , w m ) 1 = w 1 ≥ w 2 ≥ . . . ≥ w m-1 ≥ w m = 0 , (2)
W C = (w 1 , . . . , w m ) w ∈ W D ∧ w i -2w i+1 + w i+2 ≥ 0 ∀i ∈ [[m-2]] . (3)

Dominance and Possible Winners

In this section we discuss dominance relations and possible winners. The first step is to reformulate scoring rules in terms of cumulative ranks. This is useful in order to establish dominance relations between alternatives in the context of scoring rules. Dominance [START_REF] Stein | A stochastic dominance analysis of ranked voting systems with scoring[END_REF] holds between alternative x and alternative y iff the former has higher score than the latter alternative for any possible scoring rule; y is then said to be dominated. A rational decision maker will then never choose a dominated alternative; since, no matter how the weights are defined, there is another alternative that is at least as good (in case of weak dominance) or strictly better (strict dominance). We will then refine undominated alternatives in order to identify which candidates are possible winners.

Reformulation using Cumulative Standings

We consider cumulative standings that represent the fraction of times that a candidate was ranked over a certain point. Cumulative standings are defined as the cumulative sum of the rank vectors, starting from the first position. The vector

V x = (V x 1 , . . . , V x m-1 ) is such that V x j = j l=1 v x l
is the number of times that alternative x has been ranked in position j or better. Note that we consider in V x only m -1 components: indeed V x m , the number of cumulative standings in the last place, would always equal to the number of voters n.

We now define the vector δ = (δ 1 , . . . , δ m-1 ), dubbed differential weights, as the vector of the differences between two successive positional weights of the original scoring vector:

δ j = w j -w j+1 ∀j ∈ [[m-1]].
Remember that we assumed that w 1 = 1 and w m = 0; this assumption means that we have δ 1 = 1-w 2 and δ m-1 = w m-1 . The score obtained by an alternative

x can now be expressed in function of δ and the cumulative standings V :

s δ (x) = m-1 j=1 δ j V x j
where we use the subscript δ to underline the dependency on the differential weight. Decreasing weights in the original formulation correspond to positive differential weights. Note that the original rank vectors can be expressed in terms of cumulative standings: v j = V j -V j-1 ; similarly the original weights can be recovered from the differential weights: w j = m-1 l=j δ l . The requirement that the first weight should be equal to one can therefore be written as: m j=1 δ j = 1 (and this indirectly bounds all other weights to be lower than 1 due to monotonicity). The score of an alternative can therefore be seen as a convex combination of the cumulative ranks.

We now derive a reformulation that is useful for scoring vectors that constitute convex sequences. To do so we need to introduce the vector of cumulative of the cumulative standings:

V x j = j l=1 V x l = j l=1 l o=1 v x o = j l=1 (j -l + 1)v x l j ∈ [[m]]
We now define a new vector of parameters φ that represents the convexity of the weights w:

φ j = δ j -δ j+1 = w j -2w j+1 + w j+2 j ∈ [[m -2]] φ m-1 = δ m-1 = w m-1 . (4)
Given that we assumed w 1 = 1, we have φ 1 = 1 -2w 2 + w 3 and, since w m = 0, we have

φ m-2 = δ m-2 -δ m-1 = w m-2 -2w m-1 .
Stating that φ j ≥ 0 for all components j is equivalent to require the weight vector w to be convex. The score s w (x) of an alternative x under a scoring rule with weight vector w can now be expressed as m-1 j=1 φ j V x j with φ obtained from w using Equation 4. The constraint that all weights should be bounded and the highest weight, w 1 has value equal to one, becomes, in terms of convex weights, in assuming that m-1 l=1 l φ l = 1. We now show, in the following table, how some common scoring rules are expressed in terms of differential weights and in terms of convex weights.

) δ = (1, 0, . . . , 0 m-2 ) φ = (1, 0, . . . , 0 m-2 ) k-approval w = (1, . . . , 1 k , 0, . . . , 0 m-k ) δ = (0, . . . , 0 k-1 , 1, 0, . . . , 0 m-k-1 ) φ = (0, . . . , 0 k-2 , -1, 1, 0, . . . , 0) m-k-2 Borda w = (1, m-2 m-1 , m-3 m-1 , . . . , 0) δ = ( 1 m-1 , . . . , 1 m-1 ) φ = (0, . . . , 0 m-2 , 1 m-1 ) top-k Borda w = (1, k-2 k-1 , . . . , 1 k-1 , 0, . . . , 0 m-k ) δ = ( 1 k , . . . , 1 k k , 0, . . . , 0 m-k-1 ) φ = (0, . . . , 0 k-1 , 1 k , 0, . . . , 0 n-k-1 )
We highlight the following observations:

k-approval is not convex when k ≥ 2. The k-approval score of alternative x is exactly V x k .

we call top-K Borda the scoring rule based on Borda restricted to the top k positions. The score of an alternative x with respect to top-k Borda is

V x k k .

Dominance relations

The usefulness of the reformulations presented in Section 2.1 is that they can be used to discriminate candidates according to dominance relations, that allow to identify candidates that are less preferred than another one for any feasible scoring vector. Note that dominance only gives us with a partial order, so it is usually not enough to unambiguously define a winner.

Non-increasing weights When dealing with scoring vectors in W D , i.e. with non increasing sequences of weights the set of possible scores obtained by a candidate x, with cumulative ranks V x , are given by

m-1 j=1 δ j V x j δ 1 ≥ 0, . . . , δ m ≥ 0 ∧ m j=1 δ j = 1 .
Basically, all convex combinations of the components of V x are possible. Since all elements of δ are non-negative, we can compare the cumulative ranks of two candidates componentwise to check if a dominance relation exists.

Proposition 1 [START_REF] Stein | A stochastic dominance analysis of ranked voting systems with scoring[END_REF] If V x V y then the score x is necessarily at least as good than y for any scoring rule with non-increasing weights (and x is necessarily strictly better than x when V x V y ):

-V x V y =⇒ s w (x) ≥ s w (y) ∀w ∈ W D , -V x V y =⇒ s w (x) > s w (y) ∀w ∈ W D .
We say that x weakly dominates y in the first case, and that x strongly dominates y in the second case.

The previous statement can be seen as a form of first-order stochastic dominance.

A candidate is said dominated if there exists another candidate that dominates the former.

Convex weights Assuming a non-increasing convex sequence (the scoring vector belongs to W C ), the space of possible scores associated to an alternative is given by

m-1 j=1 φ j V j φ 1 ≥ 0, . . . , φ m ≥ 0 ∧ m-1 j=1 j φ j = 1 .
All φ j are non-negative since the sequence is convex. If each element of the vector V x is at least as big as the corresponding element of V y , than x has at least the same score of y for any scoring rule with convex weight (and the analogous relation holds with strict inequalities).

Proposition 2 [START_REF] Stein | A stochastic dominance analysis of ranked voting systems with scoring[END_REF] If V x V y , then x is at least as good than y, for any scoring rule with a convex sequence of decreasing weights (and x is necessarily strictly better than x when V x V y ):

-V x V y =⇒ ∀w ∈ W C s w (x) ≥ s w (y), -V x V y =⇒ ∀w ∈ W C s w (x) > s w (y).
Again, the first case is referred as weak dominance and the second as strong dominance. This is akin to second order stochastic dominance, but considering convex and not concave utility.

Example 1. Consider the following numeric example.The first table reports the distribution of the ranks, the second the cumulative ranks and the third the double cumulative ranks. -When considering monotone weights, one can establish dominance by pairwise comparisons of rows in the table of cumulative ranks; for instance a weakly-dominates c since

Candidate v • 1 v • 2 v • 3 v • 4 a 2 2 2 2 b 0 6 2 0 c 2 0 4 2 d 4 0 0 4 Candidate V • 1 V • 2 V • 3 a 2 4 6 b 0 6 8 c 2 2 6 d 4 4 4 Candidate V • 1 V • 2 V •
V a 1 = V c 1 , V a 2 > V c 2 , and V a 3 = V c 3 .
The set of weakly-undominated candidates is then {a, b, d}. No strong domination holds.

-When considering convex weights, now d weakly dominates a; moreover d strongly dominates c since

V d j > V c j for all j ∈ [[m-1]
]. The set of weakly undominated candidates is {b, d}, while the set of strictly undominated is {a, b, d}.

Possible Winners

In the following we present the notions of possible and necessary winners under different assumptions about the scoring vector. The possible winners are all those candidates that may be winners under a realization of the weights; a necessary winner (if it exists) is a winner under any possible weight. Possible winners (and necessary winners) need to be undominated according to the relations described in the previous section. However it is important to note that there might be undominated alternatives that may not be a possible winner.

In what follows, let let W be the set of feasible weights; we focus on W being either W C or W D . Definition 1. A candidate x is a possible co-winner iff there is a weight vector w such that the score of x under w is higher or equal than the score of all other candidates.

∃w ∈ W : s w (x) ≥ s w (y) ∀y ∈ A

If the above formula holds with a strict inequality sign, x is a possible winner.

Definition 2. A candidate x is a necessary co-winner iff for all weight vectors w ∈ W the score of x under w is higher or equal than the score of all other candidates. ∀w ∈ W : s w (x) ≥ s w (y) ∀y ∈ A

If the above holds with a strict inequality sign, x is a necessary winner.

Possible winners are a subset of maximal elements of the dominance relations seen in the previous sections; indeed a dominated candidate (in either sense, weak or strong) cannot be a possible winner for any scoring rule. A candidate that is weakly dominated by another candidate but it is not strongly dominated may be a possible co-winner. Moreover, note that if there is only a single candidate that is undominated, then it is a necessary winner. We show that it may occur that an undominated alternative is not a possible winner with an example.

Example 2. Assume the rank distributions (n = 12, m = 3) associated to three candidates a, b, c presented in the following table (to be read as follows: a is ranked first 4 time, second 3 times and 5 times last). Can a be a winner for some values of w 2 ? The answer is no: for a to be better than b, we need w 2 ≥ 2 3 but for a to be better than c, w 2 should be less than 1 3 . The cumulative ranks are V a = (4, 7), V b = (6, 6), V c = (2, 11) and no pairwise domination holds between a, b, or c. Therefore a is neither dominated by b or c but is not a possible winner (and not even a co-winner). This kind of reasoning has been discussed in [START_REF] Benabbou | Incremental weight elicitation for multiobjective state space search[END_REF], dealing with combinatorial problems; we now present similar techniques for computing possible winners.

v • 1 v • 2 v •

Computation of Possible Winners

We test whether candidate x is a possible winner by finding the maximum difference between its score and that of the best ranked alternative other than x. If this is positive, then candidate x is a possible winner, if it is zero is only a possible co-winner; otherwise it is not a possible winner. Formally, define the maximum advantage or margin MA(x) of x by: MA(x) = max

where W ∈ {W D , W C } is the set of possible weight vectors (either the class of non-increasing or the class of convex weights). This can be achieved by the following optimization: max Note that when considering convex weights, we reduce the feasibility region and therefore the maximum margin will be less. Candidates b, d are possible winners in both cases, since they are associated with a positive maximum advantage. Note that even if a was undominated in W D , it is not a possible winner; in W D a and c are possible co-winners for scoring weight w = (1, 0.5, 0.5, 0).

Z,w Z (7) s.t. Z ≤ m j=1 w j v x j - m j=1 w j v y j ∀y ∈ A-{x} (8)
w j ≥ w j+1 (9)
w j -2w j+1 + w j+2 ≥ 0 ∀j ∈ [[m -2]] (10)
w 1 = 1; w m = 0 (11)
. Candidate v • 1 v • 2 v • 3 v •

Aggregation of scoring rules with uncertain weights

In general there may be many possible winners, it is therefore often necessary to have a method to pick a single winner. We now discuss how to adapt classic criteria for decision-making under uncertainty for the case of scoring rules with uncertain weights.

Maximin and maximax Selecting the alternative according to maximin (criterion that picks the alternative whose worst score is highest) corresponds to setting each weight to zero, except w 1 that is equal to 1 by assumption; therefore this case corresponds to plurality. Instead maximizing the maximum possible score3 (that means using maximax decision rule) corresponds to setting each weight to one, except w m that is set to 0 by assumption; this case corresponds to (m-1)approval for the case of non-increasing scoring rules. For convex rules, maximax is attained by Borda. Of course, these methods are trivial and do not actually support the idea that we should possibly consider a variety of scoring vectors (in W m or in W c depending on the case) for evaluating different candidates.

Minimax regret Minimax regret is a robust decision criterion classically used for optimization under uncertainty [START_REF] Savage | The Foundations of Statistics[END_REF][START_REF] Kouvelis | Robust Discrete Optimization and Its Applications[END_REF]; it has been more recently advocated to be used in decision-making where the uncertainty is over utility values [START_REF] Boutilier | Constraint-based Optimization and Utility Elicitation using the Minimax Decision Criterion[END_REF][START_REF] Salo | Preference ratios in multiattribute evaluation (PRIME)-elicitation and decision procedures under incomplete information[END_REF]. In the context of voting, it has been proposed by Lu and Boutilier [START_REF] Lu | Robust approximation and incremental elicitation in voting protocols[END_REF] as a way to deal with partially elicited rankings (but they assume that the social choice function is known precisely).

We consider now using minimax regret to identify a candidate (to be declared winner) in face of uncertainty over the values of the scoring vector w 1 , . . . , w m . The idea is to associate each alternative with the maximum loss incurred in terms of score points (with respect to the "true" winner) assuming that an adversary can freely set the scoring vector; the minimax-regret alternative is the one that minimizes such loss. The max regret of alternative x is:

MR(x) = max
w j v x j ( 13
)
where W is either W D (decreasing weights) or W C (convex weights), depending on the context. We then pick the alternative with minimum max regret: arg min x∈A Z MR (x); whatever the weight in W , the picked alternatives is at most MR(x) score points from optimality.

We now show that we can characterize the minimax regret alternative in a way that it is not necessary to solve an optimization problem. As before, we slice out our analysis dealing with 1) non-decreasing weights and 2) convex sequences. Before presenting our results note that the winner according to kapproval is then just the alternative with highest value V k (number of times x was ranked in position k or better); let

V * k = max x∈A V x k be such value. Similarly let V * k = max x∈A V
x k be the maximum value of the double cumulative ranks for a given position j.

Proposition 3 In the case of non-increasing weights:

MR(x) = m-1 max j=1 V * j -V x j ; (14)
instead, assuming convex weights:

MR(x) = m-1 max j=1 V * j -V x j j . (15)
Proposition 3 allows 4 us to compute max regret in an efficient way without solving an optimization problem. While in the original formulation of maximum regret (Equation 13) the maximization is over a continuum of values (all feasible scoring vectors w), Equations 14 and 15 compute the maximum among a fixed number of alternatives (in the two cases, monotone and convex).

For non-increasing scoring vectors, note that that, since V * j is equivalent to the best score according to j-approval, the term V * j -V x j is the loss occurred to x when considering j-approval; max regret can be seen as the maximum loss occurred by x with respect to the family of k-approval voting rules. The minimax regret optimal alternative is then the candidate that is the least far away from the optimal score attained with any k-approval voting rule. Example 4. We provide an example of max regret computation using the first case of Proposition 3. One needs to consider the table of the cumulative standings; the max-regret of a given candidate is the the maximum shortfall between the values in the candidate's row compared to the best value (in bold below) in each column. Since we wish to minimize the maximum regret MR, the optimal alternative w.r.t. minimax regret is therefore a when considering non-increasing weights. Now, note that, since V x j j represent the score associated to the top-j Borda rule, Equation 15 basically states that the computation of max regret when w lies in W C is equivalent to consider the loss occurred by x with respect to the family of top-k Borda aggregators (that includes plurality and Borda as a special case).

Candidate v • 1 v • 2 v • 3 v • 4 a 2 2 2 2 b 0 6 2 0 c 2 0 4 2 d 4 0 0 4 Candidate V • 1 V • 2 V •
Example 5. We consider again the running example and compute the max regret values when considering that the scoring vector is a convex sequence, using Equation 15. First, we compute the double cumulative distribution of the ranks; we then divide the second column by two, the third by three, etc. Then, for computing the maximum regret of a, we consider the maximum, among columns, between the value of a and the best (bold) values, i.e. max{4-2, 4-3, 4.66-4} = 2. The minimax regret optimal candidate in this case is therefore d. Expected score Assuming a distribution over the weights, it is possible to rank alternatives by the expected score. Let P (w) be the such distribution. By the linearity of expectation, we can sort alternatives by using a scoring rules using the expected values of the uncertain weights:

Candidate

V • 1 1 V • 2 2 V •
E w∼P (w) [s w (x)] = m j=1 E w∼P (w) [w j ]v x j .
Therefore it means that sorting the candidates by expected score -under distribution P (w) -is equivalent to using a scoring rule whose weights are given by E[w j ], the expectations of the weights. The choice of the probability distribution over the weights is critical; we observe that if we assume an uniform distribution over non-increasing weights, sorting by expectation is equivalent to Borda.

Proposition 4 Assume w uniformly distributed in W D (non-increasing weights)

E w∼U [s w (x)] = s Borda (x)
For convex sequences, we don't have a closed formula for the expected weights in W C . We rely on Monte Carlo methods based on Gibbs sampling to derive numerical values for the expected weights. In Figure 1 we show the weights obtained by using expectation assuming an uniform distribution over all possible convex sequences. In the plot we compare the expected weights with that (normalized) used in official F1 car races (that adopt a fixed convex sequence). The F1 weights (dashed lines) can be considered rather "steep" since they only award points to a small number of top drivers while often there are about 15-20 drivers (cfr with expected weights with m = 22, in light blue); but if one wishes to award points only to the best 8 or 9 drivers, then the F1 weights are actually less "steep" than the ones obtained by expectation over convex sequence (in the plot consider the line representing m = 9, in violet).

Numerical Tests

We consider the rankings of the F1 race championship from 1961 to 2008 5 . In this context, a candidate is a driver and each ranking is the result of a race. In Table 1 below, we show the number of undominated alternatives and the number of possible winners in the case of monotone weights and that of convex weights. Several alternatives are Pareto dominated (e.g. a driver is Pareto dominated if there is another one that is better ranked in all races); the number of Pareto candidates is shown on the fourth column (for instance 16 m w1 w2 w3 w4 w5 w6 w7 w8 w9 w10 3 1.00 0.25 0.00 -------4 1.00 0.39 0.12 0.00 ------5 1.00 0.47 0.21 0.07 0.00 -----6 1.00 0.55 0.29 0.13 0.04 0.00 ----7 1.00 0.57 0.32 0.17 0.08 0.02 0.00 ---8 1.00 0.61 0.39 0.23 0.13 0.06 0.02 0.00 --9 1.00 0.64 0.42 0.27 0.16 0.09 0.04 0.01 0.00 -10 1.00 0.69 0.48 0.33 0.21 0.13 0.07 0.03 0.01 0.00 We then analyze (using the propositions of Section 2) which drivers are undominated, and which are possible winners (Section 2.3). For monotone scoring vectors, there are often several possible winners but never more than 8. In just 5 years we have a necessary winner (1962,1963,1991,1993,2002), meaning that setting the scoring vector correctly is critical. When assuming convex scoring vectors, the number of possible winners is considerably reduced (and there is a necessary winner in several circumstances). Note that, although not frequently, there are sometimes undominated candidates that are not possible winners.

In Table 2 (top) we analyze how often the different rank aggregation methods (including ranking with the F1 point system and k-approval) disagree6 on picking the winner; while most of the time they pick the same winners, there is a considerable difference between. There is a considerable disagreement between Borda and the expected score assuming convexity (more than 1 out of 4 times the winner is different), and as well between Borda and convex MR. Since scoring rules can also be used to produce rankings, we also compare (Table 2 bottom) the rankings obtained by the different methods (e.g. sorting by ascending MR scores) using the Spearman distance.

Conclusions

In this article we have considered preference aggregation in the context of social choice when the scoring vectors are not defined apriori. The fact that using scoring rules with different weights may produce very different outputs have been noticed before; some of these methods are reviewed in [START_REF] Llamazares | Preference aggregation and DEA: an analysis of the methods proposed to discriminate efficient candidates[END_REF]. In this paper we first considered the problem of discriminating which alternatives are possible winners, and the considered aggregation methods motivated by decision-theory (in particular, minimax regret and expected score).

In this paper we take the assumption that a committee specifies a generic settings for the scoring vector, that is using monotone or convex sequences; we believe that this is quite reasonable if one wish to use our proposed methods as social choice function. For minimax regret, the computation is quite easy using the cumulative (and double cumulative) ranks, therefore one could adopt such technique as an aggregation method. One could of course add additional desiderata about the weights or even consider elicitation protocols -in that case one would need to compute regret using optimization tools [START_REF] Viappiani | Regret-based optimal recommendation sets in conversational recommender systems[END_REF].

Note that the provided methods could be also used with an informative goal: for instance, for inspecting which candidates are possible winners (letting the decision maker discriminate between such alternatives in a second moment). The fact that the expectation of a uniformly distributed decreasing sequence gives the weights that are linearly spaced (as Borda), can be seen as an additional theoretical justification for using Borda.

  As usual we let w 1 = 1 and w 3 = 0, so the only free parameter is w 2 . The score of a is s(a) = 4 + 3w 2 ; the score of b is s(b) = 6 and the score of c is s(c) = 2 + 9w 2 .

Fig. 1 .

 1 Fig. 1. Table (on the left): weights obtained by uniformly sampling a convex sequence of m positions; plot (right): sampled weights with the weights used by F1 races.

Table 1 .

 1 drivers out of 54 are Pareto in 1961, while 17 out 22 in 2008). Number

				Monotone		Convex
	dataset m n Pareto #Und. #P.W. #Und. #P.W.
	1961	54 8	16	3	3	2	2
	1962	51 9	12	1	1	1	1
	1963	54 10	8	1	1	1	1
	1964	42 10	17	6	5	4	4
	1965	52 10	13	2	2	2	2
	1966	43 9	20	4	4	2	2
	1967	47 11	14	4	4	3	3
	1968	44 12	22	4	4	2	2
	1969	42 11	15	4	4	1	1
	1970	43 13	21	8	8	5	4
	1971	50 11	21	4	4	1	1
	1972	44 12	25	7	7	1	1
	1973	44 15	19	4	3	1	1
	1974	62 15	23	3	3	2	2
	1975	52 14	22	6	6	1	1
	1976	57 16	23	6	6	3	3
	1977	61 17	26	5	5	5	4
	1978	49 16	26	6	6	1	1
	1979	36 15	20	3	3	3	3
	1980	41 14	23	2	2	1	1
	1981	40 15	22	5	5	3	3
	1982	40 16	26	6	6	4	4
	1983	35 15	26	3	3	2	2
	1984	35 16	24	3	3	2	2
	1985	36 16	24	3	3	1	1
	1986	32 16	20	4	4	2	2
	1987	32 16	25	8	8	2	2
	1988	36 16	16	2	2	2	2
	1989	47 16	30	5	5	2	2
	1990	40 16	20	5	5	3	3
	1991	41 16	19	1	1	1	1
	1992	39 16	23	5	5	1	1
	1993	35 16	24	1	1	1	1
	1994	46 16	24	4	4	2	2
	1995	35 17	26	6	6	1	1
	1996	24 16	20	2	2	1	1
	1997	28 17	22	4	4	2	2
	1998	23 16	17	3	3	2	2
	1999	24 16	18	2	2	2	2
	2000	23 17	21	4	4	2	2
	2001	26 17	19	2	2	1	1
	2002	23 17	5	1	1	1	1
	2003	24 16	16	3	3	2	2
	2004	25 18	16	2	2	1	1
	2005	27 19	19	2	2	1	1
	2006	27 18	19	2	2	1	1
	2007	26 17	17	4	4	3	3
	2008	22 18	17	4	4	2	2
				of undominated items
	and possible winners in actual rank data.
	(Legend: #Und.=number of undominated,
	#P.W.:=number of possible winners,
	m=number of drivers, n=number of races).

This hypothesis is removed in Goldsmith et al.[START_REF] Goldsmith | Voting with rank dependent scoring rules[END_REF] where the authors allow preferences for intermediate positions.

There is some redundancy in the constraints: it is enough to assume convexity and wm-1 ≥ 0 to ensure that the sequence is not increasing.

Several authors, including[START_REF] Cook | A data envelopment model for aggregating preference rankings[END_REF], have proposed to take a similar optimistic approach, although the way the feasible set is defined makes the resulting rules quite different.

Obtained from the PREFLIB data repository (http://www.preflib.org/).

To handle the case of ties (when a method returns multiple winners) we compute disagreement as the cardinality of the symmetric set difference normalized by the cardinality of the union.

Disagreement

Monotone Convex F1 [START_REF] Goldsmith | Voting with rank dependent scoring rules[END_REF]


Raphael Cohen
email: raphael.cohen@gatech.edu

Eric Feron
email: feron@gatech.edu

Pierre-Loïc Garoche

Eric Feron Dutton-Ducoffe

Credible Autocoding of The Ellipsoid Algorithm Solving Second-Order Cone Programs

The efficiency of modern optimization methods, coupled with increasing computational resources, has led to the possibility of real-time optimization algorithms acting in guidance of systems. Unfortunately, those algorithms are still seen as new and obscure and are not considered as a viable option for safety critical roles. This paper deals with the formal verification of convex optimization algorithms. Additionally, we demonstrate how theoretical proofs of real-time convex optimization algorithms can be used to describe functional properties at the code level, thereby making it accessible for the formal methods community. In seeking zero-bug software, we use the Credible Autocoding framework. We focused our attention on the Ellipsoid Algorithm solving second-order cone programs (SOCP). The paper also considers floatingpoint errors and gives a framework to numerically validate the method.

I. INTRODUCTION

Formal verification of optimization algorithms used online within control systems is the sole focus of this research. Recently, such algorithms have been used online with great success for the guidance of dynamical systems, including, autonomous cars [START_REF] Falcone | Predictive active steering control for autonomous vehicle systems[END_REF] and reusable rockets [START_REF] Blackmore | Lossless convexification of control constraints for a class of nonlinear optimal control problems[END_REF]. Thus, powerful algorithms solving optimization problems are already used online, have been embedded on board, and yet still lack the level of qualification required by civil aircraft or manned rocket flight. Automatic code generation for solving convex optimization problems has already been done [START_REF] Mattingley | Automatic code generation for real-time convex optimization[END_REF], but does not include the use of formal methods. Likewise, work within the field of model predictive control already exists where numerical properties of algorithms are being evaluated [START_REF] Patrinos | A dual gradient-projection algorithm for model predictive control in fixed-point arithmetic[END_REF]. Nevertheless, this work is only valid for Quadratic Programming and using fixed-point numbers. As well, no formal verification is performed. On the other hand, some contributions already have been made concerning formal verification of control systems [START_REF] Feron | From control systems to control software[END_REF], [START_REF] Wang | From Design to Implementation: An Automated, Credible Autocoding Chain for Control Systems[END_REF], but focuses on formal verification and code generation for linear control system and typical pierre-loic.garoche@onera.fr feedback control techniques. Research has also been made toward the verification of numerical optimization algorithms [START_REF] Wang | Credible autocoding of convex optimization algorithms[END_REF], yet it remains theoretical and no proof was performed. Work already exists about formal verification for convex optimization algorithms [START_REF] Cohen | Formal verification for embedded implementation of convex optimization algorithms[END_REF]. This work remains unfortunately very basic. Only a formalization of annotations at code level are presented and no proof was actually performed with the use of formal methods. As well, only an initial formulation of a numerical analysis is given, without concrete result. The paper is structured as follows: at first, Section II presents backgrounds for convex optimization and axiomatic semantics. Section III focuses on the axiomatization of an optimization problem and the formal verification of the ellipsoid method. Section IV presents a modification of the original ellipsoid method in order to avoid ill-conditioned ellipsoids. Following this, a floating-point analysis is presented in Section V. Finally, Section VI concludes.

II. PRELIMINARIES

A. Second-Order Cone Programming

Optimization algorithms solve a constrained optimization problem, defined by an objective function, the cost function, and a set of constraints to be satisfied:

min f o (x) s.t. f i (x) ≤ b i for i ∈ [1, m]
This problem searches for x ∈ R n , the optimization variable, minimizing f o ∈ R n → R, the objective function, while satisfying constraints f i ∈ R n → R, with associated bound b i . A subclass of these problems can be efficiently solved: convex problems. In these cases, the functions f o and f i are required to be convex [START_REF] Boyd | Convex optimization[END_REF]. Here, we only present a specific subset of convex optimization problems: Second-Order Cone Programs. For x ∈ R n , a SOCP in standard form can be written as:

min f T x s.t. A i x + b i 2 ≤ c T i x + d i for i ∈ [1 , m] With: f ∈ R n , A i ∈ R ni×n , b i ∈ R ni , c i ∈ R n , d i ∈ R.
Because we are focusing on SOCP, this work does not include semidefinite programs (SDP). In control systems, SDP's are mostly used off-line, checking beforehand system's stability [START_REF] Boyd | Linear matrix inequalities in system and control theory[END_REF]. Thus, SOCP represents a trade off between being general and useful for the control community.

B. Axiomatic Semantic and Hoare Logic

Semantics of programs express their behavior. Here, we specify a program using axiomatic semantics. In this case, the semantics can be defined in an incomplete way, as a set of projective statements, ie. observations. This idea was formalized by [START_REF] Floyd | Assigning meanings to programs[END_REF] and then [START_REF] Hoare | An axiomatic basis for computer programming[END_REF] as a way to specify the expected behavior of a program through preand post-condition, or assume-guarantee contracts.

Hoare Logic: A piece of code C is axiomatically described by a pair of formulas (P, Q) such that if P holds before executing C, then Q should be valid after its execution. This pair acts as a contract for the function and (P, C, Q) is called a Hoare triple. In our case we are interested in specifying, at code level, algorithm specific properties such as the convergence or feasibility. Software frameworks, such as the Frama-C platform [START_REF] Cuoq | Frama-c: a software analysis perspective[END_REF], provide means to annotate a source code with these contracts, and tools to reason about these formal specifications. For the C language, ACSL [START_REF] Baudin | ACSL: ANSI/ISO C Specification Language[END_REF] (ANSI C Specification Language), can be used as source comments to specify function contracts, or local annotations such as loop invariants. In this work, we used the WP Frama-C plugin, and the SMT (Satisfiability modulo theories) solver Alt-Ergo to prove properties at code level.

Linear Algebra-based Specification: ACSL also provides means to enrich underlying logic by defining types, functions, and predicates. In its basic version, ACSL contains no predicate related to linear algebra. It is however possible to introduce such notions, supporting the expression of more advanced properties. Figure 1 presents the definition of new ACSL types and formalization of matrix addition. We also wrote a library for operators used in convex optimization, defining norm, grad, scalar product, det, etc. Figure 1 presents as well the definition of the vector two norm.

III. AUTOMATIC ANNOTATED CODE GENERATION

A. The Ellipsoid Method Solving SOCP

Let us now recall the main steps of the algorithm detailed in [START_REF] Nemirovski | Introduction to Linear Optimization[END_REF]. In the following, we denote E k = Ell(B k , c k ), the ellipsoid computed by the algorithm at ACSL 1 / * @ 2 type matrix; 3 type vector; 4 logic vector vec_of_3_scalar(double * x) 5 reads x[0..2]; 6 logic real vector_select(vector x, 7 integer i); 8 logic integer vector_length(vector x); 9 logic vector vec_add(vector x, vector y); 10 axiom vec_add_length:

Ell(B, c) = {Bu + c : u T u ≤ 1} (1)
Ellipsoid cut: We start the algorithm with an ellipsoid containing the feasible set X, and therefore the optimal point x * . We iterate by transforming the current ellipsoid E k into a smaller volume ellipsoid E k+1 that also contains x * . Given an ellipsoid E k of center c k , we find a hyperplane containing c k that cuts E k in half, such that one half is known not to contain x * . Finding such a hyperplane is called the oracle separation step, cf. [START_REF] Nemirovski | Introduction to Linear Optimization[END_REF]. In our SOCP setting, this cutting hyperplane is obtained by taking the gradient of either a violated constraint either the cost function. Then, we define the ellipsoid E k+1 by the minimal volume ellipsoid containing the half ellipsoid Êk that is known to contain x * that is computed thanks to the Equations ( 2), ( 3) and ( 4). In addition to that, we know an upper bound, γ, of the ratio of Vol(E k+1 ) to Vol(E k ) (see Property (1)). Figure 2 illustrates such ellipsoids cuts.

c k+1 = c k -1/(n + 1)B k p (2) and B k+1 = n √ n 2 -1 B k + n n + 1 - n √ n 2 -1 (B k p)p T (3) with: Fig. 2. Ellipsoid Cut p = B T k e/ e T B k B T k e. ( 4
)
Property 1: Let k ≥ 0, by construction:

Vol(E k+1 ) ≤ exp{-1/(2 • (n + 1))} • Vol(E k )
Please find the proof of this property in [16].

Hypotheses: In order to know the number of steps required for the algorithm to return an -optimal solution, three scalars and a point x c ∈ R n are needed:

• a radius R such that X ⊂ B R (x c ) • a scalar r such that B r (x c ) ⊂ X • and V such that max x∈X f o (x) -min x∈X f o (x) ≤ V .
The main result can be stated as:

Theorem 1: Let us assume that X is bounded, not empty and such that R, r and V are known. Then, for all > 0, the algorithm, using N iterations, will return x, satisfying:

f o (x) ≤ f o (x * ) + and x ∈ X ( -solution)
With N = 2n(n + 1) log R r V , n being the dimension of the optimization problem. This result, when applied to LP, is historically at the origin of the proof of the polynomial solvability of linear programs. Its proof can be found at [START_REF] Nemirovski | Introduction to Linear Optimization[END_REF], [START_REF] Leonid G Khachiyan | Polynomial algorithms in linear programming[END_REF].

B. Semantics of SOCP

To axiomatize an optimization problem we intend to see it, independently of the method used to solve it, as a pure mathematical object. Our goal is to axiomatize it with enough properties so that when coupling it with the ACSL annotated C code implementation of a solving algorithm, the resulting code would be formally verifiable. Using ACSL, we define a new type and a high level function, providing the possibility to create objects of the type "optim" (Figure 3). When applying a method to solve an actual optimization problem, many concepts are crucial. The work here is to highlight those concepts and write a complete enough library so that when adding the specifications and the semantics of the implementation of a method to this library, the code would be formally verifiable. The concepts of feasibility and optimality are being axiomatized. For instance, please find in Figure 4 the axiomatization of a constraint computation and the feasibility predicate definition.

C. Annotating the Ellipsoid Algorithm

In order to annotate the algorithm and produce a globally formally verifiable code, we adopt a specific technique. Each function will be written in separated files. For each function, two files are generated: a header and a body file. On the header file (.h), the header of the function and the ACSL contract will be written and on the body file (.c) appears the ACSL annotated C code implementation of the function (updateEllipsoid.c, getGrad.c, ini-tializationEllipsoid.c, ellipsoidMethod.c, etc). Figure 5 displays ACSL contract for the ellipsoid method. To C Code + ACSL 1 / * @ 2 requires \forall vector x,y; keep things simpler, we only presented annotations in the case where we were working with linear programs (only defined by matrix A, vectors b, c) and not with second-order cone programs. As well, we simplified the annotations in order to make it more readable.

Progress: So far, all the matrix-based properties have been proven (addition, multiplication, scalar product, . . . ). As well, the volume related triples are also proved by the SMT solver Alt-Ergo. The constraint and cutting hyperplane vector computations are proved as well. The loop invariant stating the optimality of c best is also demonstrated. On the other hand, covering ellipsoids property and the -optimality of returned solutions is not quite proved by itself for now. These properties being mathematically complex, further work is required in order to have them proved by the SMT solver.

IV. BOUNDING THE CONDITION NUMBER

Bounding the condition number of B is fundamental and represents the main argument of the algorithm numerical stability. Unfortunately, for the original algorithm, no reasonable bound on k(B) can be found. Therefore, we slightly modified the ellipsoid algorithm to make it able to correct the current ellipsoid E k in the case where its condition number had become too high (ellipsoid too flat). That way we can control the condition number of B. Also, this correcting step, when it occurs, does not break the convergence of the algorithm and its semantics described in Section III-A. First let us define the condition number of a matrix:

Definition 1: For a non-singular matrix A of R n×n , we define the condition number of A as the scalar:

k(A) = A • A -1 = σ max (A)/σ min (A) (5)
Where σ i are the singular values of the matrix.

A. Bounding the Singular Values

When updating B i by the usual formulas of the ellipsoid algorithm, B i evolves according to a) Minimum Half Axis: First, we claim that if σ min (B) is less than than r /V then the algorithm has already found an -solution (we can thus stop the algorithm and return the current best point found). The scalar being the wanted precision and the scalars r and V being defined in Section III-A (proof in [16]).

B i+1 = B i D i , where n -1 singular values of D i are n/ √ n 2 -1,
b) Maximum Half Axis: When the largest singular value of B i is less than 2R √ n + 1, we carry out a step as in the basic ellipsoid method. When it is greater, we take some time to "correct" B i in such a way that E + i is a localizer along with E i , specifically E i ∩X ⊂ E + i ∩X, and on the top of it:

• The volume of E + i is at most γ times the volume of E i ; • The largest singular value of B + i is at most 2R √ n + 1.
For this, let us define σ = σ max (B i ) > 2R √ n + 1 and let e o being the corresponding direction. We then consider the matrix G such that:

G = diag n/(n + 1), √ n + 1/σ, . . . , √ n + 1/σ
We conclude then this case by performing the below update on B i and c i :

B i+1 = B i • G and c i+1 = c i -(e T o c i ) • e o (6)
Please find in figure 6 an illustration of such a correction. The unit ball being the feasible set. Hence, we conclude from this that, throughout the execution of the code we have:

σ min (B i ) ≥ σ min = 1/2 • r /V σ max (B i ) ≤ σ max = 4R √ n + 1

B. Corresponding Bounds

By having a lower bound on σ min (B) and a upper bound on σ max (B), we conclude that throughout the execution of the program:

k(B i ) ≤ 2 1/2 • 2R √ n + 1 r /V = 8RV √ n + 1 r (7)
and,

B i 2 = σ max (B) ≤ 4R √ n + 1. (8)
At each iteration we know that x * ∈ Ell(B i , c i ) Which implies that:

c i ≤ R + x c + B i . (9)
V. FLOATING-POINT CONSIDERATIONS Let F denotes the set of all floating-point numbers and R the set of reals. We use standard notation for rounding error analysis [START_REF] Siegfried | Error estimation of floating-point summation and dot product[END_REF]. We write the relative rounding error unit u and the underflow unit eta. For IEEE 754 double precision (binary64) we have u= 2 -53 and eta= 2 -1074 . We present in this section an analysis targeting the numerical properties of the Ellipsoid Algorithm. Contributions already have been made concerning finiteprecision calculations within the ellipsoid method [START_REF] Leonid G Khachiyan | Polynomial algorithms in linear programming[END_REF]. However, this work only shows that it is possible to compute approximate solutions without giving exact bounds, remains very theoretical and only applied to linear programming. Also, the analysis performed considers abstract finite-precision numbers and floating-points are not mentioned. Thanks to the analysis performed in this section, using the IEEE standard for floating-point arithmetic and knowing exactly how the errors are being propagated, we would be able to check a posteriori the correctness of the analysis using static analyzers [START_REF] Putot | Static analysis-based validation of floating-point computations[END_REF].

Throughout this section, we work on the modified version of the algorithm presented in Section IV and use the presented bounds.

A. Problem Formulation

To take into account the uncertainties on the variables due to floating-point rounding, we modify the algorithm to make it more robust. For this, we choose to evaluate those uncertainties and conclude on a coefficient λ that represents by how much we are going to widen the ellipsoid E k at each iteration (see Figure 7). Within this algorithm, we focus our attention on the update formulas (2), ( 3) and ( 4). Let us assume we have B ∈ F n×n , p ∈ F n , c ∈ F n . We want to find λ ≥ 1 ∈ R such that:

Ell B + , c + ⊂ Ell λ • fl(B + ), fl(c + ) . (10)
Beforehand, we state an equivalent condition. Its proof can be found in [16].

Lemma 1: [Equivalent Condition] Ell B + , c + ⊂ Ell λ • fl(B + ), fl(c + ) ⇐⇒ fl(B + ) -1 B + u + c + -fl(c + ) ≤ λ , ∀ u ∈ B 1 (0).
∆ B = fl(B + ) -B + ; ∆ c = fl(c + ) -c + ∆ B -1 = (fl(B + )) -1 -(B + ) -1
and assume that after performing the floating-point analysis we found E B and E c such that:

|(∆ B ) i,j | ≤ E B and |(∆ c ) i | ≤ E c ∀i, j ∈ [1, n],
We dedicated Section V-B to the computation of E c and E B . Additionally, we need to compute a number E B -1 such that:

|(∆ B -1 ) i,j | ≤ E B -1 ∀i, j ∈ [1, n].
The quantity (B + ) -1 is not used explicitly in the algorithm and its floating-point error could not be evaluated by numerically analyzing the computer instructions. Instead, we will use perturbation matrix theory [START_REF] El | Inversion error, condition number, and approximate inverses of uncertain matrices[END_REF] and the theorem 2 below: Theorem 2: Let A be a non-singular matrix of R n×n and ∆A a small perturbation of A. Then, from [START_REF] El | Inversion error, condition number, and approximate inverses of uncertain matrices[END_REF], we know that,

(A + ∆A) -1 -A -1 A -1 ≤ k(A) ∆A A (11)
This will give us a lower bound on E B -1 given E B , the norm of B and its condition number. The result is stated in the following lemma (see proof in

[16]). Lemma 2: [Widening - Analytical Sufficient Condition] 1 + nE B σ max + (σ min + nE B ) √ nE c σ 2 min ≤ λ =⇒ Ell B + , c + ⊂ Ell λ • fl(B + ), fl(c + ) .
After founding such a λ, we would like to know whether the algorithm is still converging. As well, because the method's proof lies in the fact that the final ellipsoid has a small enough volume, this correction will have an impact of the number of iterations. Lemma 3 addresses those issues. Again, its proof can be found in [16].

Lemma 3:

[Convergent Widening Coefficient] Let n ∈ N, n ≥ 2.
The algorithm implementing the widened ellipsoids, with coefficient λ will converge if:

λ < exp{1/(n(n + 1))} (12)
In that case, if N denotes the original number of iteration needed, the algorithm implementing the widened ellipsoids will require:

N λ = N/ 1 -n(n + 1) log(λ) iterations (13)

B. Floating-Point Rounding of Elementary Transformations

In this section, we express the floating-point errors taking place when performing the update formulas ( 2) and [START_REF] Mattingley | Automatic code generation for real-time convex optimization[END_REF]. For this, we present first the error analysis for basic operations appearing in the algorithm. We present an axiomatization of SOCP using the specification language ACSL. In addition to that, annotations for numerical algorithms solving those problems were proposed. We focused our attention on the ellipsoid method. We presented how the process of code generation and code verification for a given optimization problem can be automated. We show how to propagate the errors due to floating-point calculations through the operations performed by the program. As it was said earlier, the proof is not finalized yet and further work is therefore required to finalize the ACSL proof of the algorithm. As well, we only presented the correctness of the optimization process for a given and fixed problem. However, when using this technique for control and hence online, we repeat this process for different initializations. Therefore, one future work includes extracting proof of convergence concerning the online utilization of those algorithms, and finally annotating the code accordingly. This would guarantee a completely sound and bug-free implementation.

18 0Fig. 1 .

 181 Fig. 1. Linear Algebra-based ACSL Specification the k -th iteration. Throughout the paper, we denote the Ellipsoid Ell(B, c) by the set:

ACSL 1 /Fig. 3 .Fig. 4 .

 134 Fig. 3. ACSL Optim Type Definition

3 isFeasible(Fig. 5 .

 35 Fig. 5. ellipsoidMethod ACSL Function Contract

  and one singular value is n/(n + 1). It follows that at a single step the largest and the smallest singular values of B i can change by a factor from [1/2, 2]. Let us argue now that one can bound the singular values of the matrix B i throughout the execution of the program.

Fig. 6 .

 6 Fig. 6. Correcting the Ellipsoid

Fig. 7 .

 7 Fig. 7. Ellipsoid Widening

  Rounding of a Real. Let z ∈ R z = fl(z) = z + δ + η with |δ| < u and |η| < eta/2 Product and Addition of Floating-Points. Let a, b ∈ F. fl(a × b) = (a × b)(1 + 2 ) + η 2 fl(a + b) = (a + b)(1 + 1 )with:| 1 | < u, | 2 | < u, |η 2 | < eta and 2 η 2 = 0 Reals-Floats Product. Let z ∈ R and a ∈ F, |fl fl(z) • a -z • a| ≤ |z||a| • u + |a| • 2u(1 + u)Following the same technique, we know evaluate the uncertainties taking place on the computations of c + and B + . ., . denotes the scalar product.Errors on c + and B + : From Equations (2) and (3), we can see that for each component of c and B, the program performs floating-point operations of the form: fl c + fl fl(z) • fl a, b and fl fl fl(z 1 ) • d + fl fl(z 2 ) • fl fl a, b • c . With: a, b ∈ F n , c, d ∈ F and z 1 , z 2 ∈ R. Hence, by propagating the errors through the elementary operations, we found:E c ≤ u • 16n 2 + 16n + 3 • B + c(14)E B ≤ u B • n 2 1 -nu + 2 |β| + n + 2|α| + 1(15)VI. CONCLUSION

ACKNOWLEDGMENTS

This work was partially supported by project NSF CPS SORTIES under grant 1446758. The authors would also like to deeply thank Pierre Roux and Arkadi Nemirovski for their help and participation in this work.