Pierre Bauguion

Amal Benhamiche

Matthieu Chardy

Yannick Carlinet

Bruno Kauffmann

Ruby Krishnaswamy

Luca Muscariello

Philippe Olivier

Nancy Perrot

Alain Simonian

Christian Tanguy

Vincent Angilella

Paul Beaujean

Mathieu Besson

Antoine Glo- Rieux

Claudio Imbrenda

Thibaut Lefebvre

Wuyang Li

Jakub Mareček

Léonce Mekinda

Yassine Naghmouchi

Felipe Olmos

Bobby Ung

Yuhui Wang

etc., thank you Xuan Zeng

for the wonderful discussions about new things, coffee breaks, dinners and jokes. Your optimism, kindness and perseverance made my PhD experience unique and unforgettable

Résumé

Contexte technologique

Le Cloud Computing, ou informatique dans le nuage, est une nouvelle tendance porteuse dans les technologies de l'information et de la communication. La virtualisation et l'utilisation des réseaux sont les technologies clés pour mettre en place une infrastructure de cloud computing. Le marché du cloud computing, qui permet de payer à l'usage de ressources de calcul et/ou de stockage connait une croissance forte. Dans ce paradigme, les serveurs de calcul ou de stockage sont partagés par plusieurs utilisateurs/clients de manière sécurisée au travers de l'utilisation de machines virtuelles (VM). L'intérêt pour les clients est la réduction de leurs investissements dans des ressources informatiques physiques en payant les ressources suivant leur utilisation. Pour l'Opérateur, le partage d'une infrastructure commune permet de mutualiser les coûts, ce qui se reflète dans des prix compétitifs. L'économie d'échelle permet aussi à l'Opérateur de réduire ses coûts et sa capacité à utiliser un nombre minimum de ressources physiques pour assurer un service à un grand nombre de clients, est essentiel pour lui assurer du profit.

Un problème central pour mettre en place une infrastructure de cloud computing est la gestion des data centers qui abritent les machines physiques. On considère ici que les data centers sont constitués d'un certain nombre de machines physiques constituées en groupes. L'organisation des data centers suit une architecture en trois tiers suivant la hiérarchie des switchs du réseau (accès, agrégation et core), architecture la plus couramment utilisée de nos jours, voir [START_REF] Headquarters | Cisco data center infrastructure 2.5 design guide[END_REF]. Typiquement, un utilisateur demande à exécuter une application qui consiste en une demande d'un nombre donné de machines virtuelles avec des demandes de ressources. La première phase du déploiement d'une application est de décider de l'affectation des machines virtuelles constituant l'application à un certain nombre de machines physiques. Généralement, ce placement se fait en optimisant un objectif (maximiser le nombre de machines physiques libres) sous des contraintes exprimant les besoins des utilisateurs et les contraintes des administrateurs du data center. Cependant, la définition de ce problème évolue avec les évolutions technologiques (possibilité de migrer les machines virtuelles dans le data center sans interrompre le service, méthodes avancées de contrôle de la consommation énergétique, etc), et l'amélioration des pratiques. De nouvelles préoccupations continuent d'apparaître, notamment sur l'évolution des fonctionnalités comme les demandes de VLANs (Virtual Local Area Networks) pour supporter un groupe de machines virtuelles, Le cloud computing a quatre composantes principales: utilisateurs, allocations des ressources, machines virtuelles et machine physiques.

Généralement, ce placement se fait en optimisant un objectif (maximiser le nombre de machines physiques libres) sous des contraintes exprimant les besoins des utilisateurs et les contraintes des administrateurs du data center. Cependant, la définition de ce problème évolue avec les évolutions technologiques (possibilité de migrer les machines virtuelles dans le data center sans interrompre le service, méthodes avancées de contrôle de la consommation énergétique, etc), et l'amélioration des pratiques.De nouvelles préoccupations continuent d'apparaître, notamment sur l'évolution des fonctionnalités comme les demandes de VLANs (Virtual Local Area Networks) pour supporter un groupe de machines virtuelles, Le cloud computing a quatre composantes principales: utilisateurs, allocations des ressources, machines virtuelles et machine physiques.

L'utilisateur soumet une demande de service qui contient les demandes détaillées de ressources (puissance de calul, mémoire, etc.). Ces ressources sont ensuite provisionnées en forme de machines virtuelles en temps réel suivant les disponibilités et un accord sur les niveaux de service selon le SLA (Service Level Agreement). Enfin, les machines virtuelles sont allouées à un ensemble de machines physiques qui exécutent la demande. Ce processus prend aussi en compte les facteurs économiques. Ce dernier concerne surtout l'opérateur dont l'objectif est de maximiser son revenu ou bien minimiser son coût opérationnel. Par exemple, réduire le nombre de serveurs physiques utilisés permet d'éteindre ceux qui ne sont pas utilisés, ou bien rentabiliser au maximum l'usage de ses serveurs physiques. Ici nous nous intéressons en particulier à l'allocation des machines virtuelles aux machines physiques. En partant d'un placement statique, un modèle devra être construit pour allouer un ensemble donné de machines virtuelles à des machines physiques. Ce placement devra prendre en compte les ressources de calcul et la capacité du réseau (mais pas les temps de calcul). L'objectif du problème est de minimiser un coût d'affectation et de routage en prenant en compte les contraintes matérielles.

Modèle mathématique

Le problème précédemment décrit se modélise sous la forme du programme mixte en nombre entier avec des contraintes quadratiques. Nous introduisons la notation suivante pour le développement de formulations mathématiques.

• R ensemble des requêtes virtuelles.

• H = (S, E) H graphe physique.

• G r = (V r , L r ) le graphe de la requête virtuelle r.

• S ensemble des serveurs.

• E ensemble des arcs du graphe H.

• V r ensemble des machines virtuelles de la requête r.

• L r ensemble des arcs de G r .

• C ri demande en ressources de calcul CPU de la VM i de V r .

• M ri mémoire requise pour VM i de V r .

• C k capacité de calcul du serveur k.

• M k mémoire disponible au serveur k.

• F k coût fixe du serveur k.

• A k coût additionnel linéaire induit par la charge CPU du serveur k.

• F rij demande de débit du lien virtuel (i, j).

• B e capacité de transmission (bande passante) du lien e.

• W e coût fixe du lien e.

• P kp le plus court chemin pour k -p, (k, p) ∈ S × S : k = p.

• x ri k 1 si le VM i de la requête r est affectée au serveur k.

• θ k 1 si le serveur k est utilisé.

• φ e 1 si le lien e est utilisé.

Le modèle résultant est donné ci-dessous:

min k∈S F k θ k + k∈S A k r∈R i∈V r c ri x ri k + e∈E W e φ e . (P)
s.t.

k∈S

x ri k = 1, ∀r ∈ R, i ∈ V r , (AC)

i∈V r

x ri k ≤ θ k , ∀r ∈ R, ∀k ∈ S, (LC)

r∈R i∈V r c ri x ri k ≤ C k θ k , ∀k ∈ S, (
θ k , φ e , x ri k ∈ {0, 1}, ∀r ∈ R, i ∈ V r , k ∈ S, e ∈ E.
L'objectif est de minimiser le coût total additif composé de trois termes: les coûts fixes encourus par commutation sur les serveurs, le coût supplémentaire provenant de la charge CPU, et les coûts fixes de l'utilisation des liens. Nous modélisons le coût supplémentaire induit par la charge CPU comme une fonction linéaire de la charge tandis que la mémoire est indépendante de la charge. Les Contraintes (AC) signifient que chaque machine virtuelle doit être affectée à un seul serveur.

(LC) modélise le fait que les machines virtuelles d'une même requête sont généralement affectées à des serveurs différents pour des raisons pratiques (sécurité, fiabilité). Les contraintes (KP, KP') expriment le fait que la mémoire et le CPU disponibles au niveau de chaque serveur sont limités. Les contraintes (QC) sont liées à la limitation de la bande passante. Ce problème est donc clairement un programme en 0,1 à contraintes quadratiques que nous nous proposons de résoudre.

Notons également que les paramètres du problème peuvent être incertains. Il est en effet difficile de déterminer à l'avance les débits requis pour chaque demande. Ceci nous incite donc à envisager les techniques d'optimisation robuste.

Outils mathématiques

Nous nous intéressons donc aux programmes mathématiques du type (1) min c T 0 x s.t.

q i (x) ≤ 0, i = 1, . . . , m,

x ∈ X ⊂ R n .

(1) Où l'ensemble q i (x), X sont donné par q i (x) = x T Q i x+c T i x+d i (i = 0, . . . , m), X = x ∈ R n : Ax ≤ b, x j ∈ Z, j ∈ J et J ⊆ {0, . . . , n}. Les c i sont des vecteurs de dimension n, et les d i sont des scalaires.

On suppose que K := X ∩{x : q i (x) ≤ 0 (i = 1, . . . , m)} est compact et non-vide. Ces programmes sont notés MIQCP. Il n'est pas trivial de relaxer ces problèmes MIQCP sous la forme de programmes convexes dans l'espace d'origine. En passant par l'idée du lifting (élevant le problème dans un espace de plus grande dimension), on peut reformuler [START_REF] Acemoglu | Opinion dynamics and learning in social networks[END_REF] en ajoutant des variables supplémentaires donnant ainsi [START_REF] Adams | A hierarchy of relaxations leading to the convex hull representation for general discrete optimization problems[END_REF]. min c T 0 x s.t.

Q i , X ≤ 0, i = 1, . . . , m,

x ∈ X ,

X = xx T , (2)
Où X est une matrice réelle de taille n × n. Notons que la contrainte X = xx T n'est pas convexe. Une approche classique consiste à relâcher cette contrainte et rajouter des inégalités valides pour renforcer la relaxation convexe. Ceci inclut l'approche de McCormick [START_REF] Mccormick | Computability of global solutions to factorable nonconvex programs: Part i-convex underestimating problems[END_REF] où la contrainte X ij = x i x j sur le domaine rectangulaire x i ∈ [l i , u i ], x j ∈ [l j , u j ] est remplacée par les inégalités max{u i x j + u j x i -u i u j , l i x j + l j x i -l i l j } ≤ X ij ≤ min{l i x j + u j x i -l i u j , l j x i + u i x j -u i l j } Une autre linéarisation a été proposée par Glover [START_REF] Glover | Improved linear integer programming formulations of nonlinear integer problems[END_REF] utilisant de l'ordre de O(n) variables. Citons également la technique de reformulation-linéarisation de Sherali et Adams [START_REF] Sherali | A hierarchy of relaxations between the continuous and convex hull representations for zero-one programming problems[END_REF]. D'autres inégalités peuvent être obtenues dynamiquement par les techniques de programmation disjonctive [START_REF] Balas | A lift-and-project cutting plane algorithm for mixed 0-1 programs[END_REF][START_REF] Balas | Disjunctive programming: Properties of the convex hull of feasible points[END_REF].

Des relaxations non-linéaires peuvent être également considérées. Il est en effet possible de relaxer la contrainte X = xx T en la remplaçant par: Xxx T 0. Cette relaxation peut être également obtenue par relaxation lagrangienne. Une autre technique consiste à convexifier les fonctions quadratiques. En effet, en considérant une matrice C i semi-définie-positive, la contrainte

q i (x) = q i (x) = x T Q i x + c T i x + d i ≤ 0 est relâchée grâce à xx T , C i + c T i x + Q i -C i , X + d i ≤ 0, X xx T .
Il est facile de prouver que déterminer la meilleure convexification se ramène à la résolution d'un programme semi-défini-positif (X -xx T 0) [START_REF] Billionnet | Exact quadratic convex reformulations of mixed-integer quadratically constrained problems[END_REF]. Une autre méthode pour obtenir des bornes pour le problème d'origine consiste à relâcher lagrangiènement certaines contraintes en les intégrant dans la fonction objectif et leur associant un ensemble de multiplicateurs exprimant le poids donné à ces contraintes de relâchées. Il s'agit ensuite de résoudre le problème dual qui consiste à trouver les meilleurs multiplicateurs maximisant la fonction objectif pénalisée. Un compromis est à trouver entre la difficulté de résoudre le problème d'évaluation (pour des multiplicateurs données) et la qualité de la borne.

D'autres paradigmes ont été récemment proposés, tel que la programmation copositive. Il s'agit de remplacer la contrainte Xxx T = 0 par une contrainte du type: S(x, X) ∈ CP, Où S(x, X) = 1 x T x X et CP est le cône des matrices complètement positives qui sont les matrices qui s'écrivent sous la forme BB T où B est une matrice réelle non négative. Le cône dual du cône des matrices complètement positives est le cône des matrices copositives qui sont les matrices X vérifiant v T Xv ≥ 0 pour tout vecteur v non négatif. Un autre paradigme utilisé dans le contexte de programmation polynomiale (incluant donc la programmation quadratique) est l'approche SoS (somme de carrés) de Lasserre [START_REF] Lasserre | Global optimization with polynomials and the problem of moments[END_REF][START_REF] Lasserre | Semidefinite programming vs. lp relaxations for polynomial programming[END_REF] et Parrilo [START_REF] Parrilo | Structured semidefinite programs and semialgebraic geometry methods in robustness and optimization[END_REF][START_REF] Parrilo | Semidefinite programming relaxations for semialgebraic problems[END_REF]. Il s'agit de remplacer la contrainte de positivité d'un polynôme par une contrainte imposant que ce polynôme est la somme de carrés de polynômes. Or si on se restreint à des sommes de carrés de polynômes de degré inférieur à une constante, on peut déterminer si un polynôme s'écrit sous la forme d'une telle somme en résolvant un programme semi-défini-positif.

L'autre outil mathématique que nous utilisons dans cette thèse pour traiter l'incertitude est l'optimisation robuste. Rappelons que les paramètres de la demande peuvent être incertains. Il est en effet difficile de déterminer à l'avance les débits requis pour chaque demande. Ceci nous incite donc à envisager les techniques d'optimisation robuste. Les variables non ajustables sont parfois interprétées comme variables du type"ici et maintenant", tandis que les variables ajustables peuvent être considérés comme des variables "attend et regarde". Cette contrepartie robuste ci-dessus est généralement appelé contrepartie robuste entièrement ajustable (FARC). FARC est parfois appelé la contrepartie robuste dynamique puisque v dépend de ξ. FARC peut être considérée comme un problème d'optimisation à deux étages où u sont les variables de premier étage et v sont les variables de deuxième étape.

Pour obtenir un problème d'optimisation traitable, plusieurs approches ont été proposées. Si nous limitons v à une fonction constante par rapport aux paramètres incertains, FARC devient simplement la contrepartie robuste standard statique notée SRC. Si les variables ajustables v dépendent affinement de ξ:

v(ξ) = w + Wξ, ξ ≡ [U, b] ∈ Ξ,
où w et W sont à optimiser, on obtient l'approche (AARC). Si v est exprimé en un polynôme en ξ ayant une degré ne dépassant pas une constante fixe, le problème de la contrepartie solide est alors lié au test de positivité d'un polynôme. Ensuite, on peut utiliser des techniques connexes, par exemple, la hiérarchie sos mentionné précédemment pour résoudre le problème. Il y a aussi d'autres approches, telles que, les politiques affines par morceaux, les approches multistatiques [START_REF] Ben-Ameur | Robust routing in communication networks[END_REF], l'adaptabilité finie [START_REF] Bertsimas | Finite adaptability in multistage linear optimization[END_REF].

Même si un grand nombre de propositions en optimisation robuste sont apparues, il y a encore des défis à relever. Tout d'abord, à notre connaissances, aucune approche n'est suffisamment générale pour englober la robustesse statique, la robustesse affine et robustesse entièrement ajustable ou dynamique décrites précédemment. Deuxièmement, comme observé dans [START_REF] Bertsimas | A hierarchy of near-optimal policies for multistage adaptive optimization[END_REF], il n'y a pas de façon systématique pour jouer sur le compromis entre la performance des politiques qui en découlent et de la complexité de la politique. Troisièmement, les paramètres incertains d'un problème d'optimisation peuvent être parfois difficiles à observer. Dans plusieurs applications, seul un sous-ensemble de ces paramètres ou certains agrégats d'entre eux peut être observé.

(programme linéaire à variables mixtes). Nous avons implémenté le relaxation de Mc-Cormick ainsi que la relaxation de Glover. Les résultats numériques montrent que, même si le modèle basé sur la relaxation de Glover a moins variables que celui basé sur la relaxation de McCormick, ses performances de calcul ne sont pas beaucoup mieux que ce dernier. La qualité de la relaxation semble être l'élément le plus déterminant. On remarque que nous pouvons appliquer les inégalités RLT pour renforcer le modèle basé sur la relaxation de McCormick et de supprimer un certain nombre de contraintes redondantes. D'un autre côté, il peut être compliqué d'appliquer RLTs au modèle basé sur la linéarisation de Alors que les inégalités valides sont souvent très efficaces pour le renforcement de la relaxation linéaire, il y a des cas où l'écart de relaxation peut être de 20%(ce qui peut rendre le branch-and-bound coûteux en temps de calcul). Ainsi, nous avons proposer des relaxations lagrangiennes pour décomposer le problème résolvant ainsi le sous-problème relâché de manière efficace avec des bornes inférieures de bonne qualité. En particulier, nous proposons une décomposition qui conduit à une séquence de sous problèmes associés à chaque demande, chaque serveur, et chaque lien. Nous donnons en outre un résultat géométrique sur la force de la borne lagrangienne montrant qu'elle est généralement beaucoup plus forte que la borne de la relaxation continue. Une hiérarchie de relaxation est également proposée en considérant une séquence de couverture de l'ensemble de la demande. Nous pouvons, par exemple, diviser l'ensemble des demandes (requêtes) en certains sous-ensembles disjoints et décomposer le problème en fonction de ce partitionnement. Les bornes de la relaxation lagrangienne résultante est meilleure que celle basée sur une simple décomposition des demandes. Pour résoudre le problème dual, nous utilisons l'algorithme in-out proposé par Ben-Ameur et Neto [START_REF] Ben-Ameur | Acceleration of cutting-plane and column generation algorithms: Applications to network design[END_REF] avec un bon choix de des paramètres, notamment le premier point intérieur.

Nos expériences numériques montrent que: (1) Pour les petites instances de problèmes que CPLEX peut résoudre à l'optimum en un temps raisonnable, nous montrons que la décomposition lagrangienne avec in-out peut fournir rapidement des bornes quasioptimales. [START_REF] Adams | A hierarchy of relaxations leading to the convex hull representation for general discrete optimization problems[END_REF] Pour les instances où CPLEX nécessite des heures de temps CPU pour les traiter, le schéma de décomposition lagrangienne avec in-out peut produire des bornes de bonne qualité en un temps moindre. De plus, on remarque que, en moyenne, la décomposition lagrangienne est environ dix fois plus rapide que CPLEX (avec ses paramètres standards) pour atteindre la même borne. La hiérarchie proposée est très attrayante pour les instances de problèmes avec un grand nombre de demandes, que CPLEX 12.6.3 ne peut résoudre dans la pratique.

Ensuite, nous introduisons une nouvelle formulation induite par des questions de symétries du problème. Cette formulation permet de réduire considérablement le nombre de termes bilinéaires dans le modèle, et comme prévu, semble plus efficace que les modèles précédents. On peut traiter le modèle résultant comme une reformulation de l'agrégation du modèle original. Pour les instances de problèmes qui ne peuvent pas bénéficier directement de cette nouvelle formulation, nous proposons une heuristique pour générer des solutions possibles et les bornes supérieures correspondantes.

Notez que, dans la pratique, l'arrivée des machines virtuelles est dynamique et imprévisible. Le modèle (P) ne peut guère fonctionner dans ce cas. A la fin du Chapitre 3, nous proposons un modèle qui optimise certaines politiques d'affectation.

Dans le chapitre 4, nous discutons de quelques approches pour l'approximation des enveloppes convexes et enveloppes concaves des fonctions bilinéaires

f (x) = i<j Q ij x i x j , (1 ≤ i < j ≤ n),
en particulier sur les hypercubes en raison de leur généralité. Une propriété de la fonction bilinéaire f dûe à Sherali [START_REF] Sherali | Convex envelopes of multilinear functions over a unit hypercube and over special discrete sets[END_REF] et Rikun [START_REF] Rikun | A convex envelope formula for multilinear functions[END_REF] est que son enveloppe convexe et concave sur un hypercube sont sommet-polyédriques, à savoir les enveloppes de f sur un hypercube coïncident avec les enveloppes de sa restriction aux sommets de l'hypercube. Cela nous permet de définir l'enveloppe convexe de f comme suit. min

λ i 2 n i=1 λ i f (v i ) : 2 n i=1 λ i v i = x, 2 n i=1 λ i = 1, λ i ≥ 0, i = 1, . . . , 2 n , où v i (i = 1, . . . , 2
n ) sont les sommets d'un hypercube à n dimensions. Cependant, comme le nombre de points extrêmes d'un hypercube est exponentiel, la caractérisation complète de ces enveloppes polyédriques est en général difficile.

Nous examinons tout d'abord un estimateur basé sur un programme SDP, où on considère l'intersection d'une contrainte de SDP avec un ensemble d'inégalités quadratiques valides pour l'hypercube. Notre contribution est liée aux travaux d' Anstreicher [START_REF] Anstreicher | On convex relaxations for quadratically constrained quadratic programming[END_REF]. Nous montrons que les deux variantes suivantes de formulations pour approcher l'enveloppe convexe des fonctions bilinéaires sont équivalentes.

min 1 2 Q, X s.t. X ii ≤ x i , (i = 1, . . . , n), max{x i + x j -1, 0} ≤ X ij ≤ min{x i , x j }, i < j, S(x, X) 0, min 1 2 Q, X s.t. X ii = x i , (i = 1, . . . , n), max{x i + x j -1, 0} ≤ X ij ≤ min{x i , x j }, i < j, S(x, X) 0
Ensuite, nous proposons une nouvelle approche en considérant un outil prédéfini pour approcher les enveloppes convexe et concave de f sur un hypercube. Plus précisément, pour un hypercube à n dimensions, nous définissons une couverture comme étant une ensemble C tel que l'ensemble des indices {1, . . . , n} est couvert par set C. Ensuite, le polytope associé à la couverture C est explicitement caractérisé par

C C = ⎧ ⎨ ⎩ (x, (X ij ) i<j ) ∈ [0, 1] n × [0, 1] n(n-1) 2 : X I = 2 |I| l=1 λ l I V l I , x I = 2 |I| l=1 λ l I v l I , 2 |I| k=1 λ k I = 1, λ I ≥ 0, I ∈ C. ⎫ ⎬ ⎭ , où V l I = v l I (v l I )
T , et v l I est un sommet de l'hypercube de dimension |I|. Un estimateur convexe de fonction bilinéaire f pourrait être donné par :

1 2 min (x,(X ij ) i<j )∈C C Q, X .
Cette approche conduit à une reformulation linéaire et nous sommes en mesure d'établir des liens avec diverses inégalités et la méthode reformulation-linéarisation-technique (RLT) proposée dans [START_REF] Adams | A hierarchy of relaxations leading to the convex hull representation for general discrete optimization problems[END_REF]. Enfin, on montre que l'enveloppe convexe (resp. Concave) d'une fonction bilinéaire f est affinement équivalente à l'enveloppe concave (resp. Convexe) d'une fonction enveloppe polyédrique. Plus précisément, nous considérons la fonction polyédrique

φ : [0, 1] n x → φ(x) = i<j Q ij |x i -x j |.

Nous avons montré que

f (x) = 1 2 i<j Q ij (x i + x j ) - 1 2 q φ(x), ∀x ∈ [0, 1] n , q f (x) = 1 2 i<j Q ij (x i + x j ) - 1 2 φ(x), ∀x ∈ [0, 1] n ,
où q f ( f ) est l'enveloppe convexe (concave) d'une fonction f . Avec ces résultats, nous montrons que si Q ij ≥ 0, l'enveloppe concave d'une fonction bilinéaire f est entièrement caractérisé par la relaxation de McCormick.

Nos résultats numériques montrent que dans l'ensemble, les estimateurs basés sur la couverture sont généralement performants.

Après les résultats numériques, nous fournissons quelques idées pour une nouvelle approche pour approximer l'enveloppe convexe à un point x prédéfini lorsque tous les coefficients de f ne sont pas négatifs. L'approche est motivée par la propriété que l'enveloppe concave de la fonction φ est sommet-polyédrique sur toutes les polytopes arbitraires si les coefficients Q ij (i < j) sont non négatifs. Nous considérons un polytope qui couvre un hypercube et a moins de 2 n sommets. Ensuite, nous montrons que l'enveloppe convexe f sur ce polytope (qui peut être obtenu facilement) est une borne inférieure de l'enveloppe convexe de f . D'autres études sont nécessaires pour trouver un moyen raisonnable de définir des polytopes appropriés.

Dans le chapitre 5, nous considérons les programmes linéaires impliquant des paramètres incertains et proposons une nouvelle contrepartie robuste traitable pour approcher la stratégie dynamique (FARC). Pour introduire le nouveau paradigme, nous introduisons deux outils, à savoir, la matrice de l'ombre et l'ensemble de pôles. Une matrice de l'ombre pourrait être utilisée pour représenter le fait que les actions ou décisions sont généralement prises sur des informations incomplètes de l'incertitude. Si la matrice d'ombre est une matrice d'identité dont la taille est le nombre de paramètres incertains, nous disons que nous avons une mesure complète de l'incertitude. Désignons par P ∈ R n 0 ×dim(Ξ) une matrice d'ombre, où n 0 est la dimension de l'ombre (à savoir, l'information partielle résultant) et dim(Ξ) est la dimension de l'ensemble d'incertitude Ξ. L'information résultante partielle est définie par

Ξ P := PΞ ≡ {Pξ, ξ ∈ Ξ}.
Un élément clé de l'approche multipolaire est un ensemble fini de pôles, qui sont des vecteurs dans l'espace image de la matrice d'ombre. On note Ω un tel pôle-set. Nous disons que Ω est un ensemble de pôles Ξ P si pour tout ξ ∈ Ξ, Pξ appartient à l'enveloppe convexe de Ω (une combinaison convexe de pôles) notée conv Ω. Étant donné un ensemble Ξ P , une collection d'ensembles de pôles Ξ P est définie comme

F Ξ P := {Ω : Ξ P ⊆ conv Ω} .
Évidemment, les points extrêmes de Ξ P forment un possible ensemble de pôles, à savoir, ext(Ξ P ) ∈ F Ξ P .

Pour un Ω ∈ F Ξ P fixé, pour chaque ξ ∈ Ξ un coefficient positif λ ξ ω est associé à chaque pôle ω. Par définition d'un ensemble de pôles, le système suivant a une solution

ω∈Ω λ ξ ω ω = Pξ, ω∈Ω λ ξ ω = 1, λ ξ ω ≥ 0, ω ∈ Ω.
Soit Λ ξ l'ensemble des vecteurs de poids Λ ξ satisfaisant le système ci-dessus pour une ξ ∈ Ξ donné. Dans le paradigme considéré, chaque pôle est associé à une action de recours, et l'action de recours lorsque ξ ∈ Ξ est considéré, est approchée par une combinaison convexe des opérations de recours associées aux pôles. Plus précisément, soit v ω l'action de recours associée au pôle ω dans le système ci-dessus. Nous avons besoin des variables ajustables v(ξ) afin de se limiter à

v(ξ) = ω∈Ω λ ξ ω v ω ,
Où λ ξ ∈ Λ ξ . Nous pouvons facilement présenter la contrepartie robuste multipolaire définie par

Π Ξ (P, Ω) = min u,v c T u s.t. Uu + V ω∈Ω λ ξ ω v ω ≤ b, ξ ∈ Ξ, λ ξ ∈ Λ ξ .
La contrepartie robuste multipolaire cherche une paire de solution non-ajustable u et un ensemble d'actions de recours liées aux pôles v ω , ω ∈ Ω tel que la fonction objectif est minimisée lorsque les paramètres incertains varient dans l'ensemble d'incertitude Ξ. En bref, étant donné Ξ, l'approche robuste multipolaire peut être vue comme une fonction d'ensemble d'un Ω pôle-set et une matrice d'ombre P . On note la fonction par Π Ξ : R n 0 ×dim(Ξ) × F Ξ P (P, Ω) → Π Ξ (P, Ω) ∈ R et on désigne par Π Ξ (P, Ω) la valeur robuste par rapport à (P, Ω). En outre, nous appelons (u, v) solution multipolaire. Nous avons ensuite prouvé que la contrepartie robuste multipolaire généralise la contrepartie robuste statique, contrepartie robuste affine ajustable, ainsi que contrepartie robuste entièrement ajustable. Par ailleurs, dans le cadre de l'optimisation robuste multipolaire, on montre que si l'ensemble d'incertitude Ξ est un simplex, la contrepartie robuste affine ajustable est équivalente à la contrepartie robuste entièrement ajustable dans le sens où les valeurs des fonction objectif sont égales à l'optimum.

Nous montrons également que la contrepartie robuste multipolaire est traitable soit par une procédure de génération de coupe ou d'une formulation compacte. Un résultat important est que, étant donné une matrice d'ombre P, la fonction Π Ξ est monotone croissante par rapport à l'inclusion de l'enveloppe convexe des ensembles de pôles. En outre, nous montrons que l'approche multipolaire peut générer une séquence de bornes supérieure et une séquence de bornes inférieures en même temps et les deux séquences convergent vers la valeur robuste des FARC sous certaines hypothèses modérées.

Comme indiqué précédemment, l'approche multipolaire est basée sur la définition d'un ensemble de pôles. Pour le construire, nous commençons par générer un simplex puis nous calculons la meilleure transformation homothétique de ce simplex lui permettant de contenir un ensemble convexe donné. Un algorithme efficace est proposé pour calculer cet ensemble homothétique. En tant que conséquence directe de notre algorithme, nous offrons une preuve simple des résultats géométriques de [START_REF] Nevskii | Properties of axial diameters of a simplex[END_REF] liés aux hypercubes. Les ensembles de pôles obtenus après cette transformation homothétique sont alors améliorés en utilisant une procédure de raffinement basée sur la génération de coupes qui est utilisée pour mettre à jour l'ensembles de pôles.

Pour illustrer numériquement l'approche multipolaire, un problème de lobbying est considéré où un lobby vise à minimiser le budget nécessaire pour convaincre un ensemble d'électeurs en tenant compte d'un modèle de dynamique d'opinion raisonnable sous une certaine incertitude. Dans nos expériences numériques, nous utilisons notre contrepartie robuste multipolaire pour modéliser le problème. Nous considérons des hypercubes ou des boules comme ensemble d'incertitude. Généralement, la contrepartie robuste multipolaire ferme de manière significative l'écart entre contrepartie robuste affine et la contrepartie robuste entièrement ajustable, en particulier pour les incertitudes du type hypercube. En outre, comme on pouvait s'y attendre, la valeur robuste diminue à mesure que plus d'informations est utilisée dans la contrepartie robuste multipolaire. Enfin, nous avons illustré le bénéfice de l'adaptabilité.

Dans l'ensemble, Résumons les principaux résultats de cette thèse comme suit. Many real life problems are characterized by making decisions with current information to achieve certain objectives. Mathematical programming has been developed as a successful tool to model and solve a wide range of such problems. However, many seemingly easy problems remain challenging. And some easy problems such as linear programs can be difficult in the face of uncertain inputs. Motivated by a telecommunication problem where assignment decisions have to be made such that the cloud resources are utilized in a minimum-cost way, we employ several mathematical programming tools to solve the problem efficiently and develop new tools for general theoretical problems. In brief, our work can be summarized as follows.

• An exact formulation and several reformulations on the cloud virtual machine mapping problem. We use several valid inequalities to strengthen the exact formulation, thereby accelerating the solution procedure significantly. In addition, an effective Lagrangian decomposition and a symmetry-induced formulation are proposed.

• A couple of new perspectives on the construction of convex and concave envelopes for bilinear optimization over a hypercube. We establish several theoretical connections between different techniques and provide a novel approach.

• A novel paradigm on general linear problems with uncertain parameters. We propose a hierarchical and convergent framework of adjustable robust optimizationmultipolar robust approach, which generalizes notions of static robustness, affinely adjustable robustness, fully adjustable robustness and fills the gaps in-between.

Our goal is to provide a mapping solution of cloud virtual machines to cloud operators such that the resource of data center networks are utilized economically and efficiently in a cloud computing environment. This easy-to-state problem however is hard to solve even for small size problem instances as it is a nonconvex quadratically constrained problem. Then, several questions arise. How to improve the scalability of computation? How to model the problem in the face of uncertainty? Can we generalize some techniques or provide new perspectives to some theoretical problems? This thesis is planned to address these questions.

Thesis motivation

This thesis is motivated by an optimization problem in the context of cloud computing where a number of different sets of virtual machines need to be assigned to servers in a date-center network. In this section, we present the problem and some background information. Then, we review some existing work and point out several challenges.

Mapping of cloud virtual machines

Cloud computing paradigm emerges as a combination of many advanced technologies, e.g., hardware virtualization, Web services, distributed computing and date center automation. It significantly improves the scalability and flexibility of the utilization of network resources while cutting the IT cost as a pay-as-you-go service.

In a cloud computing environment, clients of telecommunication operators demand cloud resources for different applications. Cloud resources primarily include memory and CPU, which are usually packed as virtual machines (VMs). To meet practical requirements (e.g., data transmission, live migration), communications between virtual machines are required. Therefore, a request from a single client can be described as a set of virtual machines with associated network throughput.

For each demand from a client, cloud operators organize their activities in a top-down process. It involves pricing, accounting, resource allocation and mapping procedures. In the mapping procedure, cloud operators map the virtual requests to physical network machines by the technology of virtualization, which vertically scales the usage of hardware resources (CPU and memory).

Today, with an increasing number of demands on cloud services, the utilization of server virtualization in Infrastructure as a Service (IaaS) and Platform as a Service (PaaS) grows rapidly. Therefore, improving the scalability and reducing the power consumption of data center networks becomes a practical issue to cloud operators. Several solutions are available to alleviate this situation. For example, one can change the existing network architecture and routing protocol. Alternatively, one may consider optimizing the placement of VMs on the hosts of virtualised servers under multiple constraints while minimizing the power consumption. This thesis focuses on the assignment of virtual resource to physical servers. In general, a Virtual Network Function (VNF) is responsible for the implementation of functionalities of multiple VMs (a request) on top of the hardware networking infrastructures. Possible mapping solutions are limited to certain hard constraints, for example, resource capacity constraints, traffic routing constraints. The solution to this problem is how to map the cloud resources and which servers and links should be used. Our objective is to minimize the power consumption cost incurred by the mapping decisions, which will also be called mapping cost. And henceforth, we will call this problem mapping problem.

To conclude this section, we illustrate the mapping problem by an example in Figure 1.1. In the example, we have two virtual requests which are also two VNFs. The left-hand-side one has three VMs and the right-hand-side one has two VMs. Also, VMs within the same request communicate with each other. One feasible solution is that for

Problem background

In this section, we provide some basic information on virtual requests.

As Cloud Computing is built on the virtualization technologies, VMs are important components. After a customer's request is received, a set of VMs will be allocated on servers to execute a specific program, even on different operating systems. Without special restrictions, a server can usually run multiple VMs simultaneously. During operating, one VM is isolated from the other, but some communication throughput will be generated between VMs within the same request. Note also that in practice, each request involves a small number of VMs.

Apart from the basic requirements, some additional features are required for advanced computing services. In order to improve the utility of data center resource, VMs should be dynamically started or stopped, and sometimes live migration should be conducted, i.e. move the VM from one physical to another. To ensure the service availability, VMs of the same request should be separately allocated.

In general, VMs are usually categorized by the type, e.g., small, medium, large, xlarge [START_REF] Wang | The impact of virtualization on network performance of amazon ec2 data center[END_REF] (see Table 1.1). In addition, one may estimate the distribution of the throughput based on real statistics. It has been observed that in practice, while 80% around 800 Kbytes/min, 4% of them have a rate 10 times higher (see [START_REF] Meng | Improving the scalability of data center networks with traffic-aware virtual machine placement[END_REF] for example). Furthermore, one may need to be aware that all the VMs communications between two corresponding servers (which host VMs) should route on a single path, as multipath routing may cause discrepancies among arrivals of data at destination.

Related work and challenges

Virtualization technology enables the emergence of cloud computing as a flexible and ondemand service. It has been observed that virtualization can also impact the computation and communication performance of cloud services. These impacts are greatly due to the placement of virtual machines in a virtualization based data-center network, which has attracted a considerable amount of studies. Those studies differ from each other in terms of the objective function, resource associated constraints, and solution methods. We review in this section several investigations related to the mapping problem.

Google in 2012 proposed a challenge organized by the French Operational Research and Decision Aid Society (ROADEF) and the European Operational Research society (EURO), where a set of VMs assigned to a set of machines needs to be reassigned to minimize the assignment cost while balancing the usage of physical machines under several resource constraints. As reported in [START_REF] Murat Afsar | Machine reassignment problem: the roadef/euro challenge 2012[END_REF], the proposal takes into account capacity constraints of servers regarding CPU, memory, storage as well as location constraints. However, it does not include bandwidth constraints respecting the throughput requirements among VMs. Exact formulations of this problem are mixed integer linear programs. To deal with large scale problems, different heuristics were proposed. Similar work has been conducted in [START_REF] Espling | Modeling and placement of structured cloud services[END_REF], where authors constructed a structural representation of virtual services. The root of the structure holds all the information associated with the structure and each child is associated with VM type information and constraints. However, the problem formulation does not concern the constraint of bandwidth incurred by communication flows. One of the specific challenges on the collaborations among data centers is that a local cloud site has no information on the specifics of remote sites (i.e. the amount of available capacity). To address the scalability of data centers, authors in [START_REF] Meng | Improving the scalability of data center networks with traffic-aware virtual machine placement[END_REF] introduced the traffic-aware virtual machine placement model where only bandwidth constraints are considered and it results in a Quadratic Assignment Problem (QAP). Rather than solving the problem exactly, a two-tier approximate algorithm was proposed. Furthermore, traffic patterns and network architectures are taken into account as some impacts of the scalability.

In the context of virtual network embedding (VNE), similar mapping problems are also considered. Houidi et al. [START_REF] Houidi | Virtual network provisioning across multiple substrate networks[END_REF] proposed a general approach to virtual network provision among multiple substrate networks. The virtual provision process is characterized by four stages: resource matching, splitting, embedding and binding. In the phase of VNE, a mixed integer linear programming model was developed and an exact algorithm was proposed for the solution. Their global objective is to maximize the acceptance ratio. They also considered the dynamic case, where the arrival of VMs is dynamic and unpredictable over time. Heuristics on sequential request processing and VM embedding with parallel request processing where discussed. Note that due to the scalability of cloud computing, a graph splitting procedure is not necessarily considered in the cloud environment. Authors in [START_REF] Papagianni | On the optimal allocation of virtual resources in cloud computing networks[END_REF] presented a unified resource allocation framework for IaaS clouds based on VNE techniques. They applied the methodology proposed in [START_REF] Chowdhury | Virtual network embedding with coordinated node and link mapping[END_REF] to correlate the node and link mapping phase, leading to a mixed integer program. The formulation allows splittable traffic flows (i.e., multiple routing paths) for any O-D pair. Authors considered two phases solution procedure: node mapping and link mapping. In node mapping phase, they applied random rounding techniques [START_REF] Raghavan | Randomized rounding: A technique for provably good algorithms and algorithmic proofs[END_REF] to correlate flow variables and binary variables. A potential substrate node which maximizes the products of flow variables and binary variables is then selected. Decision on the mapping of virtual links are made by solving a Multi-Commodity Network Flow (MCNF) problem.

A dynamic placement model for clustered web applications was proposed in [START_REF] Karve | Dynamic placement for clustered web applications[END_REF].

In their proposal, a placement controller decides to start or stop servers as needed to control the number and placement of application instances, which reduces overheads (e.g., configuration file loading), but it does not allow combinations of applications to be deployed and executed on application servers. They considered multiple objectives, i.e., maximize the acceptance ratio, minimize the placement changes, and balance the loads of servers, while respecting some linear knapsack constraints. To solve the problem, a heuristic method was proposed.

In the face of uncertain demands of VMs, authors in [START_REF] Chaisiri | Optimal virtual machine placement across multiple cloud providers[END_REF] presented a model to place VMs across multiple cloud providers, where they assumed four components in the cloud environment: user, VM repository, cloud provider, and cloud broker. The algorithm was designed for cloud brokers to allocate VMs onto cloud providers and it assumes that there are three phases provisioning resources: reservation, utilization and on-demand plans. Whenever the demand exceeds the amount of reservation resources, additional prices (usually greater than the sum of the corresponding reservation and utilization cost) are incurred. They used stochastic integer programming techniques (a two-stage recursive formulation) to minimize the total cost. The first stage variable is defined as the number of VM reservation, and the second stage variables are defined as the numbers of VMs in utilization phrase and on-demand phase. The probability distributions of uncertain parameters of demands and prices are given. Instead of considering bandwidth constraints explicitly, authors associated with each VM a bandwidth cost. mapping problem. Second, parameters of demands are usually uncertain, for example, the communication throughput between VMs fluctuates largely.

Content of the PhD 1.3.1 The scope

In this thesis, we address challenges presented at the end of last section. Point by point, let us identify the focus of our thesis.

• As mentioned earlier, the mapping problem is a 0-1 quadratically constrained problem. Its combinatorial structure is rather complex as it contains the quadratic assignment and the quadratic knapsack problems. The tractability of this problem relies heavily on the development of techniques on Mixed-Integer Quadratically Constrained Programming (MIQCP). However, as mentioned in [START_REF] Burer | The milp road to miqcp[END_REF][START_REF] Burer | On nonconvex quadratic programming with box constraints[END_REF][START_REF] Burer | Non-convex mixed-integer nonlinear programming: A survey[END_REF], even though significant progress has been made, the "breakthrough" results are yet to come and many fundamental problems have not been handled. Indeed, standard linear relaxations bring large gaps, while techniques such as semi-definite relaxation is computationally costly in the branch-and-bound procedure. Thus developing effective relaxation techniques and schemes for the deterministic model is the central task of the mapping problem.

• As MIQCP covers not only the mapping problem, but also a wide range of other applications and research problems, for example, quadratic assignment problems [START_REF] Hahn | Lower bounds for the quadratic assignment problem based upon a dual formulation[END_REF][START_REF] Wang | Optimal mapping of cloud virtual machines[END_REF], pooling and blend [START_REF] Foulds | A bilinear approach to the pooling problem †[END_REF][START_REF] Audet | Pooling problem: Alternate formulations and solution methods[END_REF][START_REF] Misener | Apogee: Global optimization of standard, generalized, and extended pooling problems via linear and logarithmic partitioning schemes[END_REF], max-cut problems [START_REF] Karp | Reducibility among combinatorial problems[END_REF], we are also interested in revisiting general relaxation techniques and establish some theoretical results for MIQCP problems.

• Motivated by the mapping problem with uncertain demands, we devote efforts to developing a general novel paradigm for linear programs with uncertain parameter in the framework of robust optimization. Linear programs with uncertain parameters can be very difficult when we consider adaptability of recourse actions. Thus, relaxation techniques need to be proposed. Moreover, in practice, the uncertain parameters of an optimization problem can be sometimes difficult to observe. So how to incorporate sufficient uncertainty to a model while allowing decisions to be made upon incomplete information? And is there a systematical way to make a trade-off between the computational complexity and the quality of the solution?

Outline and contributions

The remainder of this thesis is organized in five chapters. Chapter 2 provides an overview on relaxation techniques of MIQCPs. Chapter 3 and Chapter 4 focus on relaxations techniques on MIQCPs from an industrial perspective and a research perspective. Chapter 5 is devoted for developing novel approaches dealing with uncertainty in the framework of robust optimization.

In Chapter 2, we review some relaxation techniques on MIQCP. These techniques include general algorithmic frameworks, polyhedral relaxations, nonlinear relaxations, Lagrangian relaxation, as well as copositive programming and moment-SoS hierarchies. Derivations on some important results are given.

In Chapter 3, we firstly present several ideas to reformulate the mapping problem to get the problem solved efficiently. Further, a Lagrangian-based decomposition is proposed, reducing the problem into a number of subproblems of small sizes, which are then solved by a cutting plane algorithm. A set of experiments is carried out by implementing the branch-and-cut procedure and the accelerated cutting plane algorithm (a.k.a. in-out algorithm) with CPLEX 12.6.3 routines. The numerical results show that several global optima and lower bounds of good quality can be found in a reasonable time by the proposed approaches. Finally, we proposed a reformulation exploiting symmetries of virtual requests, which can handle some larger instances.

Chapter 4 discusses a couple of approaches to approximating the convex and concave envelopes of bilinear functions, particularly over hypercubes due to its generality. The first approach is based on a semi-definite program. The second approach considers some predefined set covers of a hypercube and leads to a linear program. Then we establish a connection between the convex envelope of a bilinear function and the concave envelope of a polyhedral function. Numerical experiments are conducted to compare the two approaches. As a perspective, a novel approach to approximate the envelopes is discussed and illustrated.

In Chapter 5, we study linear programs involving uncertain parameters and propose a new tractable robust counterpart which contains and generalizes several other models including the existing Affinely Adjustable Robust Counterpart and the Fully Adjustable Robust Counterpart. It consists in selecting a set of poles whose convex hull contains some projection of the uncertainty set, and computing a recourse strategy for each data scenario as a convex combination of some optimized recourses (one for each pole). We show that the proposed multipolar robust counterpart is tractable and its complexity is controllable. Further, we show that under some mild assumptions, two sequences of upper and lower bounds converge to the optimal value of the fully adjustable robust counterpart. To illustrate the approach, a robust problem related to lobbying under some uncertain opinions of authorities is studied. Several numerical experiments are carried out showing the advantages of the proposed robustness framework and evaluating the benefit of adaptability.

Finally, in Chapter 6, we conclude the thesis by evaluating the results and giving some directions of future research.

Introduction

A Mixed-Integer Quadratically Constrained Program (MIQCP) is an optimization problem of the following form

min c T 0 x s.t. q i (x) ≤ 0, i = 1, . . . , m, x ∈ X ⊂ R n , (2.1)
where andc i (i = 1, . . . , m) are n-dimensional vectors, d i (i = 1, . . . , m) are scalars. We do not assume any convexity property of functions q i but explicitly assume that X is a compact set in an n-dimensional real space. We further assume that the constraint set denoted by K := X ∩ {x : q i (x) ≤ 0 (i = 1, . . . , m)} is nonempty and compact, so that c T 0 x attains its minimum at some x ∈ K.

q i (x) = x T Q i x + c T i x + d i (i = 1, . . . , m), X = {x ∈ R n : Ax ≤ b, x j ∈ Z, j ∈ J},
Our assumption on set X can be justified by the following reasons. First, it is known that the decision version of a general MIQCP without restriction on the topological property of X is undecidable. Jeroslow [START_REF] Jeroslow | There cannot be any algorithm for integer programming with quadratic constraints[END_REF] showed that minimizing a linear form over quadratic constraints in integer variables is not computable by a recursive function (i.e., no computing device can be programmed to compute the optimums of all problems in the problem class.). One way to address this issue is to ensure the compactness of set X . Second, for most of practical optimization applications, the compactness of X usually holds 1 .

MIQCP constitutes a very general family of problems and there are several ways to categorize them. Given a problem instance ((q) m i=1 , X ), we can categorize it by the convexity, linearity, and integrality of the problem. • By integrality. If X does not contain integrality constraints, (2.1) becomes a continuous quadratically constrained problem.

Note that formulation (2.1) also includes optimization problems with quadratic objective function as one can move the objective function to the constraint set by introducing an additional variable. In addition, despite the fact that polynomial optimization contains the family of MIQCP problems, we notice that any polynomial optimization problem might be reduced in the form of (2.1) by introducing auxiliary variables.

Since the objective function of (2.1) is linear, the global optima of (2.1) lie among extreme points of the constraint set. Henceforth, we will focus on the constraint set of (2.1). It has been commonly acknowledged that convexity plays a fundamental role in mathematical programming in the tractability of an optimization problem. Rockafellar in [START_REF] Rockafellar | Lagrange multipliers and optimality[END_REF] pointed out that "In fact the watershed in optimization isn't between linearity and nonlinearity, but convexity and non-convexity.". However, the definition of convex optimization varies in the literature. By saying a convex problem in the form of (2.1), Recockafallar [START_REF] Rockafellar | Lagrange multipliers and optimality[END_REF] refers to cases where X is convex and q i (i = 1, . . . , m) are convex relative to set X (or equivalently, K is convex) as a convex optimization problem, while Boyd and Vanderghe [START_REF] Boyd | Convex optimization[END_REF] refers to the situation that q i (i = 1, . . . , m) are convex functions and the algebraic representation of X is convex. Following Lasserre [START_REF] Lasserre | On convex optimization without convex representation[END_REF], we say that (2.1) is a convex optimization problem if K is convex and we say that (2.1) is a convex program (or has a convex representation) if q i (i = 1, . . . , m) and the algebraic representation of X are convex.

We remark that the development of convex programming, including paradigms (e.g., linear programming, conic programming) and algorithms, has brought success in solving convex MIQCPs. But it appears that breakthroughs for nonconvex MIQCPs are yet to come.

Outline. In this chapter, we will review convex relaxation techniques handling a nonconvex MIQCP. In Section 2.3, we present two general algorithmic frameworks dealing with nonconvex MIQCP along with several definitions and basic concepts. Then we discuss both classical and recent techniques on polyhedral relaxation in Section 2.4.

In Section 2.5, we discuss several semidefinite relaxations and their equivalent reformulations, followed by dynamic strategies of adding valid inequalities. We also mention several quadratic convexification techniques and their connection with the standard semidefinite relaxation. Section 2.6 revisits the classic Lagrangian relaxation from a geometric view, which sheds light on the strength of the associated Lagrangian bound. Finally, we presents two novel paradigms briefly.

Notation. Throughout this chapter, for any matrix C, C T denotes its transpose and C ij represents its i th row j th column component. For a set S, conv S stands for its convex hull; Int S refers to its interior and cl denotes its closure. If S is finite, |S| stands for its cardinality. For a general function f : S → R, we denote by q f (resp. f ) the convex (resp. concave) envelope of f over a n dimensional hypercube. The inner product between two matrices A, B ∈ R m×n is denoted by A, B . For a symmetric matrix X ∈ S n and vector x ∈ R n , S(x, X) refers to 1 x T x X in S n+1 . S n + represents the cone of n × n real positive semidefinite matrices.

General algorithmic frameworks

It is known that a nonconvex MIQCP is much harder to solve than a convex one. For a convex MIQCP, dropping the integrality constraints usually leads to a tractable convex program. Methods that address convex MIQCPs include, for example, generalized Benders' decomposition [START_REF] Geoffrion | Generalized benders decomposition[END_REF], outer approximation [START_REF] Duran | An outer-approximation algorithm for a class of mixed-integer nonlinear programs[END_REF], Non-Linear Programming (NLP) or Linear Programming (LP) based Branch-and-Bound (BB) algorithms [START_REF] Quesada | An lp/nlp based branch and bound algorithm for convex minlp optimization problems[END_REF]. We refer readers to excellent surveys [START_REF] Bonami | Algorithms and software for convex mixed integer nonlinear programs[END_REF][START_REF] Grossmann | Review of nonlinear mixed-integer and disjunctive programming techniques[END_REF][START_REF] Belotti | Mixed-integer nonlinear optimization[END_REF] for details. In contrast, for a nonconvex MIQCP, dropping the integrality constraints leads to a nonconvex program, which is essentially a global optimization problem. To address this issue, one usually considers replacing the nonconvex functions q i (i = 1, . . . , m) with their convex underestimators.

Definition 2.1 For a nonconvex function f : X → R, its convex underestimator denoted by f over X is a convex function such that f (x) ≤ f (x), ∀x ∈ X . Its convex envelope is defined as the supremum functional among the set of all convex underestimators of f . Likewise, a concave overestimator of f denoted by f over X is a concave function such that f (x) ≤ f (x), ∀x ∈ X and its concave envelope is the infimum functional among the set of all concave overestimators.

Another methodology is to find a convex relaxation of set K. Relations between these two methodologies can usually be established by Lagrangian duality (see, for example, [START_REF] Lemaréchal | Semidefinite relaxations and Lagrangian duality with application to combinatorial optimization[END_REF][START_REF] Malick | The spherical constraint in boolean quadratic programs[END_REF]). Details about these techniques will be discussed in subsequent sections and we now briefly review two general algorithmic frameworks for exact solutions.

In general, the replacements above will produce a mixed-integer convex program. Following [START_REF] Liberti | Reformulations in mathematical programming: Definitions and systematics[END_REF], we call it an auxiliary problem of (2.1). For the following discussion, it is useful to present definitions below. Definition 2.2 [START_REF] Liberti | Reformulations in mathematical programming: Definitions and systematics[END_REF] Let P be auxiliary problem of (2.1) and F(P ) refers to the feasible region of P . If K F(P ), we say that P is a relaxation of (2.1). If problem P preserves local (resp. global) optimality information of (2.1),i.e., there is a mapping from set F(P ) to set K, such that its restriction the set of all local (resp. global) optima of P is surjective with respect to local (resp.global) optima of (2.1), then we say that P is a local (resp. global) reformulation of (2.1). If P is a local and global reformulation of (2.1), we say that P is an exact reformulation of (2.1).

Note that a relaxation of P might also be an exact reformulation of the original problem (2.1) even though the feasible region F(P ) is strictly larger than that of the original problem. A family of such relaxations has been explicitly studied and characterized by Ben-Ameur and Neto [START_REF] Ameur | A geometric characterization of "optimalityequivalent" relaxations[END_REF] and they call them optimality-equivalent relaxations.

Two general algorithmic frameworks are often employed depending on the type of the auxiliary problem P . If P is a global reformulation of (2.1), then we can employ any method dealing with convex Mixed-Integer Non-linear Programs (MINLPs) to solve P and then recover the optima of the original problem (2.1) by certain mapping function. Methods reformulating a nonconvex MIQCP to a convex MIQCP include, but not limited to, conversions to MILP [START_REF] Beale | Special facilities in a general mathematical programming system for non-convex problems using ordered sets of variables[END_REF], quadratic convex reformulations [START_REF] Billionnet | Improving the performance of standard solvers for quadratic 0-1 programs by a tight convex reformulation: The qcr method[END_REF][START_REF] Billionnet | Extending the QCR method to the case of general mixed integer programs[END_REF], copositive reformulation [START_REF] Burer | On the copositive representation of binary and continuous nonconvex quadratic programs[END_REF]. The efficiency of the solution procedure is heavily influenced by the strength of reformulations.

If P is a convex relaxation of (2.1), one can rely on global optimization methods, e.g., spatial branch-and-bound, branch and reduce. Spatial Branch-and-Bound (sBB) was first proposed by McCormick [START_REF] Mccormick | Computability of global solutions to factorable nonconvex programs: Part i-convex underestimating problems[END_REF]. Like standard branch-and-bound methods, sBB method recursively partitions the searching space in a convergent way, relying on evaluating a sequence of subproblems. Specifically, sBB first solves the convex relaxation P , which provides a lower bound. Then, it seeks to solve the (nonconvex) nonlinear problem with a local solver. If the resulting solution is feasible to (2.1), then it provides an upper bound. If the gap between upper bound and lower bound is within some predefined tolerance, the algorithm terminates and outputs the optimal solution. If not, it partitions the domain of a variable, leading to two new subproblems. Each subproblem is solved to get upper and lower bound repeatedly. In general, techniques tightening the bound of variables helps to improve the strength of convex relaxations.

Apart from the strength of convex estimators, the convergence rate also depends on the partition strategy and bounding techniques of a sBB algorithm. For example, Linderoth [START_REF] Linderoth | A simplicial branch-and-bound algorithm for solving quadratically constrained quadratic programs[END_REF] proposed a simplicial sBB algorithm for solving QCQPs over a box, where searching space is recursively partitioned into triangles and rectangles. The author showed theoretically and computationally, that the convergence of the proposed sBB algorithm is generally faster than algorithms that partition a searching space to rectangles.

Another well-known sBB algorithm in global optimization community is the so-called "αBB" method, which was proposed by Floudas and co-workers [START_REF] Androulakis | αbb: A global optimization method for general constrained nonconvex problems[END_REF]. For nonconvex MIQCPs, the principle idea of "αBB" method is to construct convex underestimator of nonconvex functions (q i ) by adding certain quadratic terms with multipliers α, such that the resulting Hessian matrix Q(α) is positive semidefinite. Investigations on the convergence of αBB algorithms and McCormick's sBB algorithm with respect to the diameter of searching region have been discussed in [START_REF] Androulakis | αbb: A global optimization method for general constrained nonconvex problems[END_REF][START_REF] Bompadre | Convergence rate of mccormick relaxations[END_REF].

We remark that the both standard BB and sBB frameworks heavily rely on the strength of relaxations. In what follows, we will discuss different convex relaxations to problem (2.1).

It is usually not trivial to relax a nonconvex MIQCP to a convex program in the original space. By the idea of lifting, one can reformulate problem (2.1) in a higher dimensional space at the expense of additional variables as follows:

min c T 0 x s.t. Q i , X ≤ 0, i = 1, . . . , m, x ∈ X , X = xx T , (2.2)
where X is an n × n real symmetric matrix. Note that X = xx T is nonconvex in variables (x, X). Obtaining a relaxation of (2.2) can usually be done by relaxing the constraint X = xx T with various valid inequalities. Of-course, one can also project valid inequalities to original space and solve a relaxed problem in the original space (see, e.g., [START_REF] Saxena | Convex relaxations of non-convex mixed integer quadratically constrained programs: extended formulations[END_REF][START_REF] Saxena | Convex relaxations of non-convex mixed integer quadratically constrained programs: projected formulations[END_REF]).

In what follows, we review several well-known methods to generate valid inequalities of different forms.

Polyhedral relaxations

The most common methods of converting a nonconvex MIQCP to a convex one might be polyhedral relaxation. Among those, we mention the following.

• McCormick in [START_REF] Mccormick | Computability of global solutions to factorable nonconvex programs: Part i-convex underestimating problems[END_REF] proposed linear relaxations of quadratic term X ij = x i x j over the rectangular domain

x i ∈ [l i , u i ],
x j ∈ [l j , u j ] by the following inequalities

max{u i x j + u j x i -u i u j , l i x j + l j x i -l i l j } ≤ X ij ≤ min{l i x j + u j x i -l i u j , l j x i + u i x j -u i l j }
These inequalities are referred to as McCormick inequalities in the literature. Note that it takes O(n 2 ) variables and constraints. McCormick also pointed out that the constraints above characterize the exact convex and concave envelope of X ij = x i x j over rectangles. One can also find detailed proof in [START_REF] Al-Khayyal | Jointly Constrained Biconvex Programming[END_REF].

• Glover in [START_REF] Glover | Improved linear integer programming formulations of nonlinear integer problems[END_REF] proposed to linearize the product of a binary variable and a bounded variable with O(n) variables. For example, let us linearize

n i=1 n j=1
x i D ij x j , where (x i , x j ) ∈ {0, 1} 2 and D is a square matrix. We introduce n new variables, w i = x i j D ij x j and add 4n constraints, i.e., for each i, we have:

D - i x i ≤ w i ≤ D + i x i , j D ij x j -D + i (1 -x i ) ≤ w i ≤ j D ij x j -D - i (1 -x i ),
where

D - i = max{ j D - ij x j }, D + i = max{ j D ij x j }.
• Sherali and Adams in [START_REF] Sherali | A hierarchy of relaxations between the continuous and convex hull representations for zero-one programming problems[END_REF] proposed the so-called reformulation-linearization techniques (RLTs) to reformulate general non-convex quadratic problems. One may treat RLTs as an unified way to generate valid inequalities in variables (x, X) by some basic algebraic manipulations (i.e., addition and multiplication) on inequalities in original space. For example, let us generate RLT inequalities for problem (2.2) by exploiting set X .

First, for interval constraints x ∈ [l, u], we multiple pairs of nonnegative terms (u i -x i ), (x i -l i ), (u j -x j ) and (x j -l j ) to get inequalities:

(u i -x i )(u j -x j ) ≥ 0, (u i -x i )(x j -l j ) ≥ 0, (x i -l i )(u j -x j ) ≥ 0, (x i -l i )(x j -l j ) ≥ 0.
Replacing x i x j with variables X ij leads to McCormick constraints mentioned above. Second, for linear equality a i x = b i , we multiply any x j at both sides to obtain quadratic equalities involving x i x j terms; For linear inequality of the form a i x ≤ b i , we can multiply both sides (x j -l j ) and (u j -x j ) respectively. In addition, one can multiple pairs of constrains within set X to get constraints involving quadratic terms.

Later Adams and Sherali [START_REF] Adams | A hierarchy of relaxations leading to the convex hull representation for general discrete optimization problems[END_REF] showed that a hierarchy of increasingly stronger relaxations can be obtained by generating higher order polynomials constraints. In the light of representations of polynomials in real algebraic geometry, Lasserre [START_REF] Lasserre | Semidefinite programming vs. lp relaxations for polynomial programming[END_REF] showed that RLT hierarchy falls exactly within the type of linear relaxation hierarchies based on Krivine-Stengle's positivity certificate. Liberti and Pantelides [START_REF] Liberti | An exact reformulation algorithm for large nonconvex nlps involving bilinear terms[END_REF] showed that one can sometimes eliminate unnecessary bilinear terms in MIQCP by augmenting a subset of RLT constraints.

• For unconstrained 0-1 quadratic programs, there has been a considerable amount of investigations devoted to generating valid inequalities. Most of them are attempts to approximate 0-1 polytopes for classical combinatorial problems. Let us mention cut polytope [START_REF] Deza | Geometry of cuts and metric embeddings[END_REF], boolean quadratic polytope [START_REF] Padberg | The boolean quadric polytope: some characteristics, facets and relatives[END_REF] (a.k.a. correlation polytope [START_REF] Deza | Geometry of cuts and metric embeddings[END_REF]) and quadratic assignment polytope [START_REF] Kaibel | Polyhedral combinatorics of the quadratic assignment problem[END_REF][START_REF] Saito | A study of the quadratic semi-assignment polytope[END_REF]. We refer readers to [START_REF] Deza | Geometry of cuts and metric embeddings[END_REF] for details. Recently, some interesting results have been established between Boolean Quadratic Polytope (BQP) and the polytope of Quadratic Problem over a Box (QPB), where

QPB = conv (x, y) ∈ R n × R n(n+1) 2 : y ij = x i x j , 1 ≤ i ≤ j ≤ n, x ∈ [0, 1] n , BQP = conv (x, y) ∈ R n × R n(n-1) 2 : y ij = x i x j , 1 ≤ i < j ≤ n, x ∈ {0, 1} n .
Burer and Letchford [START_REF] Burer | On nonconvex quadratic programming with box constraints[END_REF] proved that the projection of QPB (by ignoring elements of y ii , (i = 1, . . . , n)) is exactly BQP, which implies that all valid inequalities in variables

y ij (1 ≤ i < j ≤ n) of BQP is valid for QPB.
There are also many other techniques to generate strong linear valid inequalities in a specific context. For example, authors in [START_REF] Meier | A compact linearisation of euclidean single allocation hub location problems[END_REF] proposed a novel linearization technique for hub location problems, where they assumed that the unit transport cost between hubs is proportional to the distance between them. Specifically, they consider linearizing y ij = k,l d kl x ik x jl , where d kl represents the distance between locations k and l, x ik , x jl are 0-1 decision variables assigning hub i, j to location k, l respectively. They decomposed the distance d kl as the Euclidean norm of two vectors associated with location k and l, which can then be bounded by linear valid inequalities in a static or dynamic way.

Disjunctive cuts

A family of linear cuts can also be generated dynamically in the framework of disjunctive programming. Classical disjunctive programming [START_REF] Balas | Disjunctive programming[END_REF] has been studied intensively in the context of MILPs [START_REF] Balas | Disjunctive programming: Properties of the convex hull of feasible points[END_REF][START_REF] Balas | A lift-and-project cutting plane algorithm for mixed 0-1 programs[END_REF]. It was then generalized by Grossmann and co-workers for general MINLPs [START_REF] Grossmann | Generalized disjunctive programming: A framework for formulation and alternative algorithms for minlp optimization[END_REF][START_REF] Grossmann | Generalized convex disjunctive programming: Nonlinear convex hull relaxation[END_REF]. Following investigations in [START_REF] Balas | Disjunctive programming: Properties of the convex hull of feasible points[END_REF][START_REF] Belotti | Disjunctive cuts for nonconvex minlp[END_REF], we review ideas from disjunctive programming and present a procedure to generate (linear) disjunctive cuts.

Disjunctions are logical operations (composed by basic logical operators "or" ) that return true if at least one of their operands is true. Classical disjunctive programming [START_REF] Balas | Disjunctive programming: Properties of the convex hull of feasible points[END_REF] requires a linear relaxation and a disjunction of (2.1). Let us denote by LP a polyhedral relaxation of K. Regarding problem (2.1), a disjunction can be written as

F = q k=1 F k , and F k = x ∈ LP : H k x ≤ h k (k = 1, . . . , q)
where each

H k is a m k × n real matrix, h k ∈ R m k , m k is the number of linear inequalities.
We say that this disjunction is satisfied if there exists k such that F k = ∅. For any x ∈ K, disjunction F is satisfied. A single problem may has many disjunctions. For example, a standard BB procedure partitions the searching space by imposing integer disjunctions on integer variable x, x i ≤ α ∨ x i ≥ α + 1, α ∈ Z in a finite number of steps [START_REF] Belotti | Disjunctive cuts for nonconvex minlp[END_REF]. For a nonconvex MIQCP, disjunctions are usually taken among non-convex constraints. We postpone related techniques in Section 2.5.1.

Given a current point x ∈ LP, the central task in disjunctive programming is to generate a cut π T x ≤ π 0 for cl conv F that is violated by x. This cut is referred to as disjunctive cut. Assume that LP := x ∈ R n + : Bx ≤ β with B ∈ R m×n and β ∈ R n . Balas [START_REF] Balas | Disjunctive programming: Properties of the convex hull of feasible points[END_REF] showed that the cut identified by (π, π 0 ) is valid for cl conv F if and only if π belongs to the following set

Π := y ∈ R n : ∃(u k , v k ) ≥ 0, s.t. y ≤ (B) T u k + (H k ) T v k , (u k ) T β + (v k ) T h k ≤ π 0 , (k = 1, . . . , q)
Observe that Π is related to the reverse polar of F . Several properties of reverse polar and its application in deriving deep disjunctive cuts were given in [START_REF] Cadoux | Computing deep facet-defining disjunctive cuts for mixed-integer programming[END_REF]. Given Π, an automatic procedure for generating a cut that is maximally violated relies on solving a so-called Cut Generating Linear Programming (CGLP) problem [START_REF] Balas | A lift-and-project cutting plane algorithm for mixed 0-1 programs[END_REF]:

max (π,π 0 ):π∈Π π T x -π 0 (2.3)
If the optimum of (2.3) is nonpositive, x ∈ F , otherwise, we generate a disjunctive cut.

As remarked by Balas [START_REF] Balas | Disjunctive programming: Properties of the convex hull of feasible points[END_REF] and Saxena et al. [START_REF] Saxena | Convex relaxations of non-convex mixed integer quadratically constrained programs: extended formulations[END_REF], a normalization constraint of the form q k=1 ((u k ) T ξ + (v k ) T ξ k ) = 1 is usually added to the constraint set Π to ensure the strength and numerical stability of the resulting cuts. For details of implementations, we refer readers to [START_REF] Balas | Disjunctive programming: Properties of the convex hull of feasible points[END_REF][START_REF] Belotti | Disjunctive cuts for nonconvex minlp[END_REF][START_REF] Saxena | Convex relaxations of non-convex mixed integer quadratically constrained programs: extended formulations[END_REF][START_REF] Saxena | Convex relaxations of non-convex mixed integer quadratically constrained programs: projected formulations[END_REF].

Nonlinear relaxations

In addition to polyhedral relaxations, two methodologies are frequently used to relax (2.1), namely, semidefinite relaxations and quadratic convexifications.

Semidefinite relaxations

It is well-known that semidefinite programming (SDP) can be used to derive strong relaxations for a wide range of optimization problems. Most SDP relaxations consider relaxing X = xx T to Xxx T 0, which is convex with respect to (x, X). Equivalently, one can replace the constraint X = xx T with S(x, X) 0 and rank(S(x, X)) = 1 (see, e.g., [START_REF] Vandenberghe | Semidefinite programming[END_REF]), then drop the rank-one constraint. This standard relaxation is often referred to as Shor's relaxation. It is easy to verify that this standard SDP relaxation can be derived by Lagrangian duality on problem (2.1)(see, e.g., [START_REF] Lemaréchal | Semidefinite relaxations and Lagrangian duality with application to combinatorial optimization[END_REF][START_REF] Lemaréchal | Lagrangian relaxation[END_REF]).

One can strengthen the SDP relaxations with valid inequalities, e.g., triangle inequalities, McCormick inequalities, and RLT inequalities. Anstreicher [START_REF] Anstreicher | Semidefinite programming versus the reformulationlinearization technique for nonconvex quadratically constrained quadratic programming[END_REF] showed the benefits of such combinations in terms of bound strength. Later, Anstreicher and Burer [START_REF] Anstreicher | Computable representations for convex hulls of low-dimensional quadratic forms[END_REF] showed that the combination of SDP relaxation and RLT inequalities describes the convex hull of (x, X) :

X = xx T , x ∈ [l, u] ⊂ R n if n ≤ 2.
Buchheim and Wiegele [START_REF] Buchheim | Semidefinite relaxations for non-convex quadratic mixed-integer programming[END_REF] investigated SDP relaxations for unconstrained mixedinteger quadratic programs whose constraint set is of the form X 1 × • • • × X n . In addition to the standard SDP relaxation, they propose to separate points (x i , X ii ) from each conv{(x i , x ii ) : x i ∈ X i } in a branch-and-bound scheme.

Solving large scale SDP problem is much more expensive than solving LPs or SOCPs with the same number of variables as SDP. On the one hand, special purpose solvers for conic programming cannot handle arbitrary convex constraints; on the other hand, general purpose nonlinear solvers (such as Ipopt [START_REF] Wächter | On the implementation of an interiorpoint filter line-search algorithm for large-scale nonlinear programming[END_REF]) are not designed to handle conic constraints [START_REF] Saxena | Convex relaxations of non-convex mixed integer quadratically constrained programs: extended formulations[END_REF]. Alternative to adding SDP constraint explicitly, one can also add convex quadratic constraints or linear constraints in an iterative way. Kim and Kojima [START_REF] Kim | Second order cone programming relaxation of nonconvex quadratic optimization problems[END_REF] proposed to rewrite constraint Xxx T 0 as convex quadratic constraints:

x T Cx ≤ C, X , ∀C ∈ S n + , (2.4)
or as linear constraints S(x, X), Z ≥ 0, ∀Z ∈ S n+1 + .

(2.5)

In addition, they showed that for a given C ∈ S n + , convex quadratic inequality (2.4) is stronger than linear inequality (2.5) with matrix

Z = β b/2 b/2 C , where β ∈ R, b ∈ R n . Indeed, we have 0 ≤ (1, x T )Z 1 x = xCx T + b T x + β ≤ C, X + b T x + β = S(x, X), Z
where the first inequality follows from the positive semi-definiteness of Z ∈ S n+1 + and the second inequality comes from (2.4). Authors in [START_REF] Kim | Second order cone programming relaxation of nonconvex quadratic optimization problems[END_REF] then decompose the matrix [START_REF] Horst | Dc programming: overview[END_REF], where

Q i into Q + i , Q - i by the concept of difference-of-convex programming
Q + i = n i=1:λ i >0 λ i v i v T i , Q - i = n i=1:λ i <0 λ i v i v T i
, and v i is the eigenvector corresponding to eigen-

value λ i . Therefore Q + i , -Q - i ∈ S n
+ and can be used to replace C in (2.4). Of-course, there are other ways to choose a proper psd matrix C. For example, Saxena et al. [START_REF] Saxena | Convex relaxations of non-convex mixed integer quadratically constrained programs: extended formulations[END_REF][START_REF] Saxena | Convex relaxations of non-convex mixed integer quadratically constrained programs: projected formulations[END_REF] add valid convex quadratic inequalities of the form (2.4) with C = vv T in an iterative way, where v is the eigenvector corresponding to negative eigenvalue of matrix ( Xxx T ) at a given point (x, X). Moreover, they consider adding constraints enforcing Xxx T 0 (which is nonconvex) in the spirit of disjunctive programming. Specifically, with eigenvector v, it holds that (x T v) 2 ≥ v T Xv, which reduce to w 2 ≥ z by a change of variables. Given point ( ŵ, ẑ) and assume that w ∈ [l, u], a disjunction is taken as

F 1 ∨ F 2 , where F 1 = {l ≤ w ≤ ŵ, z ≤ w(l + ŵ) -l ŵ} , and F 2 = { ŵ ≤ w ≤ u, z ≤ w(u + ŵ) -u ŵ}.
Then, one can generate a disjunctive cut by solving problem (2.3). In addition, authors [START_REF] Saxena | Convex relaxations of non-convex mixed integer quadratically constrained programs: extended formulations[END_REF] showed that adding disjunctive cuts sequentially reveals efforts of convexifications of K.

Quadratic convexifications

A number of quadratic relaxation techniques, e.g., α-BB method [START_REF] Adjiman | A global optimization method, αbb, for general twice-differentiable constrained nlps -i. theoretical advances[END_REF], the smallest eigenvalue method [START_REF] Hammer | Some remarks on quadratic programming with 0-1 variables[END_REF], difference-of-convex programming [START_REF] Horst | Dc programming: overview[END_REF], and Quadratic Convex Reformulation (QCR) [START_REF] Billionnet | Improving the performance of standard solvers for quadratic 0-1 programs by a tight convex reformulation: The qcr method[END_REF][START_REF] Billionnet | Extending the QCR method to the case of general mixed integer programs[END_REF][START_REF] Billionnet | Exact quadratic convex reformulations of mixed-integer quadratically constrained problems[END_REF], have been proposed in the literature. Among those methods, let us present QCR method [START_REF] Billionnet | Exact quadratic convex reformulations of mixed-integer quadratically constrained problems[END_REF] in detail. Assuming that x ∈ X ⊂ Z n , the so-called QCR method has two phases, the convexification phase by solving an SDP and the standard BB solution phase relying on the convex quadratic program as the output of convexification phase. The key distinction between QCR and other approaches lies in the unconventional way (to reformulate the original problem) of solving an SDP. Specifically, they rewrite each q i with a psd matrix

C i xx T , C i + c T i x + Q i -C i , X + d i ≤ 0 (2.6)
and X = xx T . Therefore, the resulting formulation is an exact reformulation of (2.2). Authors in [START_REF] Billionnet | Exact quadratic convex reformulations of mixed-integer quadratically constrained problems[END_REF] then relax X = xx T by some linear inequalities leading to a convex program.

To get the best possible choice of (C i ) m i=1 , they solve an SDP problem. Interestingly, the SDP problem is simply Shor's relaxation strengthened by linear inequalities.

To conclude this section, let us mention that SDP relaxations and quadratic convexifications are in fact two methodologies. On the one hand, SDP relaxations can be looked at as attempts to the convexification of X = xx T over set X ; on the other hand, quadratic convexifications can be regarded as efforts in constructing the convex underestimator for each q i (x) over X . Anstreicher [START_REF] Anstreicher | On convex relaxations for quadratically constrained quadratic programming[END_REF] showed that given the same information (valid inequalities), the latter is usually dominated by the former. We remark that the continuous relaxation of the quadratic convex program obtained by the QCR method gets exactly the same optimum of an continuous SDP program strengthened by linear inequalities.

Lagrangian relaxation

Lagrangian relaxation is a powerful tool in optimization, even nonconvex. As advocated by Lemaréchal, it is an essential method to generate bounds and bring new perspective on standard relaxations, e.g., SDP relaxation, to combinatorial problems (see, e.g., [START_REF] Lemaréchal | Lagrangian relaxation[END_REF]108,[START_REF] Lemaréchal | Semidefinite relaxations and Lagrangian duality with application to combinatorial optimization[END_REF]). Designing a proper Lagrangian relaxation (or decomposition) scheme is not a trivial task as one would like to have an easier problem (compared with the original (2.1)) while ensuring strong bounds. We present a couple of geometric results on the strength of Lagrangian relaxation in the literature. Let us assume that q i (i = 1, . . . , m) are complicated constraints in problem (2.1). Relaxing these constraints with a dual vector λ ∈ R n leads to the following lagrange function

L (x, λ) = c T 0 x + m i=1 λ i q i (x)
and the corresponding dual function is

Ψ(λ) :=    inf x∈X L(x, λ) λ ≥ 0, -∞, otherwise.
Recall the assumption that X is compact, therefore the minimum is attainable.

Let us introduce variables γ ∈ R m and move the nonlinear components q i (i = 1, . . . , m) to the constraint set, i.e.,

Ψ(λ) := min c T 0 x + m i=1 λ i γ i : (x, γ) ∈ S (2.7)
where

S = (x, γ) ∈ R n+m+1 : x ∈ X , q i (x) ≤ γ i (i = 1, . . . , m) .
As the objective function is linear, we can safely replace S with its convex hull conv S. Therefore, we have

Ψ(λ) = min c T 0 x + m i=1 λ i γ i : (x, γ) ∈ conv S . (2.8)
Now observe that (2.8) is exactly the dual objective function of

min c T 0 x s.t. (x, γ) ∈ conv S, γ ≤ 0.
(2.9) Clearly, Ψ is bounded. If Slater's condition holds that ∃(x, γ) ∈ conv S such that γ < 0, we can safely define max λ≥0 Ψ(λ) (2.10)

By strong duality, the optimal value of (2.10) is equal to the optimal value of (2.9). Lemaréchal [108] has shown that the constraint set of (2.9) is equivalent to S ∩ conv S, where S = {(x, γ) ∈ R n+m : γ = 0}. This reads exactly ∃(x, 0) ∈ conv S. Now notice that the original problem (2.1) is equivalent to

min c T 0 x : (x, 0) ∈ conv {(x, 0) ∈ S} (2.11)
Clearly, we can see that the duality gap between (2.10) and (2.1) is introduced due to the general relation that

conv {(x, 0) ∈ S} ⊆ {(x, 0) ∈ conv S} .
In the light of this geometric result, one can have an intuition on the quality of the bounds obtained by a Lagrangian relaxation scheme. Moreover, we can show that the dual bound provided by any Lagrangian relaxation scheme is no worse than the optimum of the continuous relaxation of the primal problem, provided that q i (i = 1, . . . , m) are convex. To formalize this result, let us denoted by X the polyhedral set obtained by dropping the integrality constraints of X .

Theorem 1 If q i (i = 1, . . . , m) are convex functions, the optimal value of (2.10) is greater than or equal to the optimum of the continuous relaxation of problem (2.1) (by replacing X with X ).

Proof: Let C be the constraint set of continuous relaxation of (2.1), i.e.,

C = (x, 0) : x ∈ X , q i (x) ≤ 0 (i = 1, . . . , m)
It is sufficient to show that {(x, 0) ∈ conv S} ⊆ C. Let (x, 0) ∈ conv S, then there exists

(n + 2) points (x k , γ k ) ∈ S, such that (x, 0) = n+2 k=1 λ k (x k , γ k ), where λ ∈ R n+2 + and n+2 k=1 λ k = 1.
Obviously, it holds that x ∈ X . For each i = 1, . . . , n, we have

q i (x) ≤ n+2 k=1 λ k q(x k ) ≤ n+2 k=1 λ k γ k = 0,
where the first inequality is due to the convexity of function q i , and the second inequality comes from the fact that (x k , γ k ) ∈ S. Therefore (x, 0) ∈ C, which completes the proof.

We now consider a special case where q i (x) (i = 1, . . . , m) are linear functions. With little trouble, one can derive the following result.

Theorem 2 [108] If all functions q i , 1, . . . , m are linear functions, the dual optimal value of Ψ(λ) is equal to the optimal value of the following primal problem

min c T 0 x : x ∈ conv X , q i (x) ≤ 0, (i = 1, . . . , m }.
This geometric interpretation of Lagrangian relaxation will be used to design a Lagrangian decomposition scheme in Chapter 3.

Recent paradigms

Recall that relaxing (2.1) amounts to approximating conv (x, X) : X = xx T , x ∈ X . In convex optimization, there exists an important subclass called conic optimization. It optimizes a linear form over the intersection of an affine subspace and a convex cone.

Looking at the aforementioned techniques from the perspective of conic optimization, we may realize that we were actually using polyhedral cone, second-order cone, and psd cone to strengthen our relaxations. A natural question arises: can we find some other convex cones that are "smaller"?

To answer this question, we briefly review two recently developed frameworks, namely, copositive programming and moment-SoS hierarchy that have been used to solve nonconvex MIQCPs.

Copositive programming

Copositive programming deals with two convex cones: the cone of copositive matrices and the cone of completely positive matrices, which are dual to each other. An n × n real symmetric matrix X is copositive if its quadratic form v T Xv is nonnegative for all v ∈ R n + . X is a completely positive matrix if it can be factorized as BB T , where

B ∈ R n×k + .
It is easy to verify that both cones are closed and convex. The cone of completely positive matrices is smaller than the cone of psd matrices, and thus one may expect stronger relaxations introduced by enforcing X being completely positive. Unfortunately, it has been recently proved that checking if a matrix belongs to the completely positive cone is NP-hard [START_REF] Dickinson | On the computational complexity of membership problems for the completely positive cone and its dual[END_REF]. And it is well-known that the membership problem of copositive matrix cone is co-NP-complete [START_REF] Murty | Some np-complete problems in quadratic and nonlinear programming[END_REF]. We refer readers to excellent surveys [START_REF] Hiriart-Urruty | A variational approach to copositive matrices[END_REF][START_REF] Dür | Copositive programming-a survey[END_REF] for details about copositive matrices and copositive programming.

Burer [START_REF] Burer | On the copositive representation of binary and continuous nonconvex quadratic programs[END_REF] showed a quite general result that every quadratic problem with linear and binary constraints has a copositive representation. More precisely, he showed that if a quadratic binary problem

min x T Qx + 2c T x s.t. c T i x = d i (i = 1, . . . , m), x ≥ 0, x j ∈ {0, 1} (j ∈ J), (2.12
) satisfies some key condition, i.e., c T i x i = d i (i = 1, . . . , m) and x ≥ 0 implies x j ≤ 1 (j ∈ J), then it can be rewritten as the following completely positive problem:

min Q, X + 2c T x s.t. c T i x = d i (i = 1, . . . , m), c i c T i , X = d 2 i (i = 1, . . . , m), x j = X jj (j ∈ J) S(x, X) ∈ CP .
where CP represents the cone of completely positive matrices and J ⊆ {1, . . . , n}. As pointed in [START_REF] Burer | On the copositive representation of binary and continuous nonconvex quadratic programs[END_REF], this key condition can be easily satisfied by augmenting the constraint x j + s j = 1 with slack variable s j ≥ 0 for all j ∈ J. Note that (2.12) encompasses a number of NP-hard problems, e.g., all 0-1 linear integer programs, all unconstrained 0-1 quadratic programs. It is still unknown if a general MIQCP problem of the form (2.1) has a copositive representation or not.

Based on the definition of completely positive matrix, one can approximate the completely positive matrix cone from outside to any precision using a sequence of polyhedralsemidefinite relaxations [START_REF] Parrilo | Structured semidefinite programs and semialgebraic geometry methods in robustness and optimization[END_REF]. For example, when n ≤ 4, it is known that the so-called doubly nonnegative matrix cone (the intersection of SDP cone and componentwise nonnegative symmetric matrix cone) characterizes the cone of completely positive matrices exactly [START_REF] Maxfield | On the matrix equation x'x= a[END_REF]. For more general cases, several approaches have been proposed to approximate the copositive cone and its dual by testing positivity of polynomials. Some of those are based on moment-SoS hierarchies, which will be reviewed subsequently.

Moment-SoS hierarchies

In this section, we briefly review the well-known relaxation techniques, moment-SoS hierarchies due to Lasserre [START_REF] Lasserre | Global optimization with polynomials and the problem of moments[END_REF][START_REF] Lasserre | Semidefinite programming vs. lp relaxations for polynomial programming[END_REF], Parrilo [START_REF] Parrilo | Structured semidefinite programs and semialgebraic geometry methods in robustness and optimization[END_REF][START_REF] Parrilo | Semidefinite programming relaxations for semialgebraic problems[END_REF]. For more details, we refer readers to the survey of Laurent [START_REF] Laurent | Sums of squares, moment matrices and optimization over polynomials[END_REF].

Let us focus on problems of the form (2.1) and introduce c(x (2.13)

) = c T x. c ∈ R[x]
where the second minimum is taken over all the probability measures µ supported by K (i.e., K dµ = 1). Indeed, it holds that K c dµ ≥ k p * dµ = p * , ∀µ ∈ M(K). On the other hand, for any x ∈ K, let µ be the Dirac measure at any point x ∈ K, i.e., µ := δ x , so we will have

K c dµ = c(x) showing that p * ≥ min µ∈M(K) c dµ. Further,
observe that (2.13) is an infinite-dimensional LP problem. Nevertheless, we have that

K c dµ = α c α K x α dµ = α c α y α , where y α = K x α dµ(i = 1, . . . , m) is called moment of order α for measure µ supported by K. Let y 0 = K dµ = 1. Problem (2.13
) then becomes a finite dimensional LP problem:

p * = min α c α y α : y 0 = 1, y ∈ M d (2.14)
where

M d = y α : ∃µ ∈ M(K) s.t. y α = K x α dµ, ∀α ∈ N n d . It is easy to verify that M d is a convex cone. So problem (2.14) is a convex problem.
Next, consider the dual of (2.13)

d * := sup γ∈R {γ : c(x) -γ ≥ 0, ∀x ∈ K} . (2.15)
Now, let us introduce the cone of nonegative polynomial with degree d (in our case, d = 1) over

K : d = {f ∈ R[x] d : f (x) ≥ 0 ∀x ∈ K}. Therefore (2.15) is equivalent to d * = sup γ∈R {γ : c -γ ∈ d } . (2.16) So (2.16
) is also a convex problem. In fact, M d and d form a dual pair of convex cones. Assume that the compact set X has an algebraic representation {x : g j (x) ≤ 0 (j = 1, . . . , k)}, then we can use some powerful results from real algebraic geometry to characterize cone d and M d . Let us introduce the quadratic module Q(q 1 , . . . , q m , g 1 , . . . , g k ) generated by q 1 , . . . , q m , g 1 , . . . , g k :

Q(q 1 , . . . , q m , g 1 , . . . , g k ) :=    f ∈ R[x] : f = - m i=1 σ i q i - k j=1 σ j+m g j + σ 0   
where (σ) m+k i=0 are sums-of-squares (SoS) polynomials in variables x. A polynomial f is a sums of squares of polynomials if it can be written as f = l j=1 u 2 j for some u 1 , . . . , u l ∈ R[x]. A nice result is that any nonnegative quadratic polynomial is a sum of squares, which follows from celebrated results of Hilbert (see, e.g., [START_REF] Laurent | Sums of squares, moment matrices and optimization over polynomials[END_REF]).

Putinar's Positivstellensatz (See [START_REF] Laurent | Sums of squares, moment matrices and optimization over polynomials[END_REF][START_REF] Lasserre | Global optimization with polynomials and the problem of moments[END_REF]) provides a positivity certificate for a polynomial over set K.

Theorem 3 [START_REF] Putinar | Positive polynomials on compact semi-algebraic sets[END_REF] Under the assumption that ∃M > 0,

s.t. M - n i=1 x 2 i ∈ Q(q 1 , . . . , q m , g 1 , . . . , g k ), for f ∈ R[x], if f > 0 on K, then f ∈ Q(q 1 , . . . , q m , g 1 , . . . , g k ).
If we replace the constraint (c -γ) ∈ d with Putinar's positivity certificate (c -γ) ∈ Q(q 1 , . . . , q m , g 1 , . . . , g k ), we get a lower bound of the dual problem (2.16). In addition, the membership problem of the quadratic module reduces to solving an SDP, provided that the degree of SoS is bounded.

u 2 j = z 2 d ( j u T j u j )z d . Therefore σ is an SoS polynomial if and only if ∃L ∈ S n + , σ = z d Lz d .
The size of the matrix L is n+d d . The strength of the relaxation is increased as the bounded degree number increases. Similar arguments can also be derived for the primal side (see [START_REF] Lasserre | Global optimization with polynomials and the problem of moments[END_REF][START_REF] Parrilo | Structured semidefinite programs and semialgebraic geometry methods in robustness and optimization[END_REF]). Consequently, the sequence of the primal side relaxations is called moment hierarchy and the sequence of the dual side relaxations is called SoS hierarchy. Notice that one can use different positivity certificates to replace Putinar's Positivstellensatz, which may reduce to different forms of convex problems, e.g., LP, SOCPs [START_REF] Lasserre | Moments, positive polynomials and their applications[END_REF].

Chapter 3

Optimal Mapping of Cloud Virtual Machines

Summary

One of the challenges of cloud computing is to assign virtual machines to physical machines optimally and efficiently. The aim of cloud operators is to minimize the mapping cost while respecting constraints regarding location, assignment and capacity. After establishing the exact MIQCP model, we propose several ideas to reformulate the problem in order to accelerate solution procedure. Further, a Lagrangian-based decomposition is proposed, reducing the problem to a number of subproblems that are significantly easier to handle and maintaining a theoretical guarantee on the quality of bounds. Also, a reformulation exploiting symmetries of virtual requests is presented to reduce the number of bilinear terms involved in the formulation. Numerical experiments are conducted showing the effectiveness of the proposed approaches.

Introduction

Since the background and a brief history of the mapping problem have been introduced in Chapter 1, we focus on the solution procedure in this chapter. The organization of this chapter and its contributions are summarized below.

In Section 3.3, we provide an exact formulation which falls within the realm of MIQCPs. To solve the problem, we propose to reformulate the problem via classical linearization techniques, thereby handling the problem in the framework of MIP. Then we employ the Reformulation-Linearization-Technique (RLT) and various valid inequalities to strengthen the model. Some numerical experiments are conducted to show the effectiveness and limitations of these methods. Section 3.4 considers a Lagrangian decomposition scheme based on each virtual request to achieve lower bounds of good quality for relatively larger problem instances. We show, both theoretically and numerically, these bounds are generally much stronger than those obtained via continuous relaxations. In Section 3.5, we introduce a novel formulation that exploits symmetries of each virtual request and allows a reduction of bilinear terms. Numerical experiments suggest that the resultant model is more scalable than the aforementioned ones. To leverage the usage of the new formulation, a heuristic is proposed.

Followed by concluding remarks, we discuss a possible formulation that optimizes certain mapping polices as a future research direction to capture the dynamic arrival of virtual requests.

A general model

Remember that a virtual request consists of a set of virtual machines and virtual communications between them. Therefore, we may represent each virtual request as a directed graph. Our goal is to map such graphs to a physical network. Henceforth, we will use the following notation to construct the mathematical model.

• Sets -R set of virtual request.
-H = (S, E) graph of a physical network.

-G r = (V r , L r ) a graph of virtual network for request r ∈ R.

-S set of servers in the physical network.

-E set of directed edges in the physical network.

-V r set of VMs of request r.

-L r set of directed virtual links of request r.

• Parameters

-c ri required CPU of VM i ∈ V r .
m ri required memory of VM i ∈ V r .

-C k CPU capacity of server k.

-M k memory capacity of server k.

-F k fixed cost of server k ∈ S.

-A k additional cost of server k imposed from CPU loads.

f rij required throughput associated with logical link (i, j) ∈ L r .

-B e bandwidth of edge e ∈ E.

-W e fixed cost of edge e ∈ E.

-

P kp shortest k -p path, (k, p) ∈ S × S : k = p. • Variables -x ri k ∈ {0, 1} 1 if VM i of request r is mapped to server k. -θ k ∈ {0, 1} 1 if server k is used (switched on).
φ e ∈ {0, 1} 1 if edge e is used (switched on).

To help readers bear in mind the notation, we remark that for a parameter or a variable, its subscripts (if it has) are associated with physical resources and its superscripts (if it has) are associated with virtual resources.

In computer networks, the traffic between an O-D pair is usually routed on a single path. As will be shown soon, the resulting formulation contains rather complicated (nonconvex quadratic) constraints. Hence, we make the following assumption in order to reduce the complexity of the problem.

Assumption 1

The traffic between any O-D pair is routed on the shortest path in terms of hops.

With the above notations, we construct the exact mathematical model of the mapping problem as follows.

min k∈S F k θ k + k∈S A k r∈R i∈V r c ri x ri k + e∈E W e φ e . (P) s.t. k∈S x ri k = 1, r ∈ R, i ∈ V r , (AC) i∈V r x ri k ≤ θ k , r ∈ R, k ∈ S, (LC) r∈R i∈V r c ri x ri k ≤ C k θ k , k ∈ S, (KP) r∈R i∈V r m ri x ri k ≤ M k θ k , k ∈ S, (KP') r∈R k,p∈S: k =p,e∈P kp (i,j)∈L r f rij x ri k x rj p ≤ B e φ e , e ∈ E, (QC) θ k , φ e , x ri k ∈ {0, 1}, r ∈ R, i ∈ V r , k ∈ S, e ∈ E. (BC)
In what follows, we refer to (P) as the above model with the interpretation below.

• The objective is to minimize the total cost, which is additively composed of three terms: the fixed cost incurred by switching on servers, the additional cost coming from the CPU load, and the fixed cost from the usage of links. We model the additional cost induced by CPU load as a linear function to represent the fact that CPU is usually categorized as load dependent resource, while memory is load independent [START_REF] Karve | Dynamic placement for clustered web applications[END_REF].

• Constraints (AC) mean that each virtual machine must be mapped to a single server. Constraints (LC) model the fact that virtual machines are usually mapped separately in a cloud environment due to some practical issues, e.g., security, reliability.

• Constraints (KP, KP') ensure that for each server, the aggregated required CPU, memory resource cannot exceed its limits. Constraints (QC) emphasis the fact that for each edge, the aggregated throughputs on the edge cannot exceed the bandwidth.

Before the solution procedure, we analyze the difficulties and structure of the problem briefly. From the perspective of mathematical programming, the combination of integrality constraints and bilinear constraints (nonconvex) makes the problem rather difficult.

In terms of the problem structure, it shares features of several well-known combinatorial problems, e.g., the multiple knapsack problem and quadratic assignment problem. Authors in [START_REF] Amaldi | On the computational complexity of the virtual network embedding problem[END_REF] have proved that the mapping problem of the form (P) is strongly NP-hard even if |R| = 1 by showing a polynomial time reduction from the maximum stable set problem.

McCormick inequalities and Glover's linearization

Model (P) is a 0-1 non-convex quadratically constrained problem. A fundamental idea to deal with such problems is lifting it to a higher dimensional space [START_REF] Burer | Non-convex mixed-integer nonlinear programming: A survey[END_REF][START_REF] Burer | The milp road to miqcp[END_REF]. By introducing new variables y rij kp and enforcing y rij kp = x ri k x rj p for each (r, i, j, k, p) : i = j, k = p, we lift the problem to a higher dimensional space and thus results in a MIP problem. This comes at a price of introducing non-convex equations. Simple convex relaxations can be achieved by linearization techniques. In our thesis, we use two linearization techniques:

McCormick inequalities [START_REF] Mccormick | Computability of global solutions to factorable nonconvex programs: Part i-convex underestimating problems[END_REF] and Glover's linearization [START_REF] Glover | Improved linear integer programming formulations of nonlinear integer problems[END_REF].

A direct use of McCormick's inequalities produces an approximation for the nonconvexity, i.e., for each (r, i, j, k, p) : i = j, k = p, we add four inequalities,

x ri k + x rj p -1 ≤ y rij kp , (3.1a)
y rij kp ≤ x ri k , (3.1b)
y rij kp ≤ x rj p , (3.1c)
y rij kp ≥ 0. (3.1d)
With this relaxation, we can solve problem (P) by a MIP model:

P MC : {P} ∩ {(
= x ri k j:j =i x rj p f rij by constraints 0 ≤ w ri kp ≤ max j:j =i {f rij }x ri k , (3.2a) j:j =i x rj p f rij + (x ri k -1) max j:j =i {f rij } ≤ w ri kp ≤ j:j =i f rij x rj p . ( 3

Polyhedral relaxations

In this section, we propose several valid linear inequalities to strengthen (P MC ). One of the most straightforward and powerful methods is to add valid linear inequalities. An important subset of those can be derived by the RLTs which is known to strengthening non-convex discrete and continuous formulations [START_REF] Sherali | A hierarchy of relaxations between the continuous and convex hull representations for zero-one programming problems[END_REF]. We apply the RLT to constraints (AC), which produces

r∈R |V r |(|V r | -1)|S| additional constraints, k∈S:k =p y rji pk = x rj p , r ∈ R, i, j ∈ V r : i = j, p ∈ S. (AC RLT )
Let us denote the formulation in the sequel as

P RLT1 : {(P)} ∩ {(AC RLT ), (3.1d)} . (P RLT1 )
Note that constraints (3.1a-3.1c) are implied by (AC RLT ) and model (P RLT1 ) is indeed stronger than formulation (P MC ). Let F (P RLT1 ) be the feasible region of (P RLT1 ) and F (P MC ) be the feasible region of (P MC ).

(y rji p k -x rj p ) ≤ y rij kp -x rj k ≤ y rij kp .
The first equality and the first inequality follows from (AC RLT ). The second equality comes from (AC).

Optimal Mapping of Cloud Virtual Machines

Of-course we can also apply the RLTs to (LC), which adds 2 r∈R |V r ||S|(|S| -1) constraints.

i∈V r :i =j

y rij kp -x rj p θ k ≤ 0, r ∈ R, j ∈ V r , (k, p) ∈ S 2 : k = p, (LC RLT ) i∈V r x ri k - i∈V r :i =j y rij kp -θ k + x rj p θ k ≤ 0, r ∈ R, j ∈ V r , (k, p) ∈ S 2 : k = p. (LC RLT )
Notice that we can simply replace the left-hand-side constraints with their convex envelopes which are readily identified by McCormick inequalities (see Section 2.4 for details). Specifically, we replace x rj p θ k in (LC RLT ) with min θ k , x rj p . Likewise, we replace

x rij p θ k in (LC RLT ) with max 0, θ k + x rj p -1 . P RLT2 : {(P)} ∩ {(AC RLT ), (LC RLT ), (LC RLT ), (LC), (3.1d)} . (P RLT2 )
Clearly, the continuous relaxation of (P RLT1 ) is stronger than that of (P MC ), but is weaker than that of (P RLT2 ). In practice, it entails a trade-off between the strength of relaxations and their sizes to determine a modest model that is robust to the inputs. Apart from linear valid inequalities coming from the RLT, some strong valid inequalities can be derived by exploiting the problem structure. Specifically, for each r ∈ R, (i, j) ∈ L r , if there is a required throughput between i and j, i.e., f rij > 0, then for each link e, we have

(k,p)∈S 2 :k =p,e∈P kp y rij kp ≤ φ e . (Cut)
The number of such additional constraints is just O(|R||E|), but surprisingly, these simple constraints can often accelerate the solution process of a standard MIP solver. Similarly, for each r ∈ R, (k, p) ∈ S 2 : e ∈ P kp , if there exists some throughput mapped to P kp , then for each link e, we have

(i,j)∈L r y rij kp ≤ φ e . (Cut')
Notice that the number of such constraints is O(|R||S| 2 |E|), a quartic polynomial in the number of physical servers, so it might be computationally beneficial to add them dynamically.

To conclude this section, we present two MIP models for the solution procedure of problem (P).

Numerical experiments

In this section, we evaluate the effectiveness and limitations of the formulations by conducting some numerical experiments on formulations (P MC ), (P RLT1 ), (P RLT2 ), (P V1 ) and (P V2 ). Regarding the problem instances, we assume that all the VMs in a request communicate with each other and fix the number of VMs for each request to be 5. The small physical network instance with 8 nodes is generated randomly. Other instances are taken from SND library [START_REF] Orlowski | SNDlib 1.0-Survivable Network Design Library[END_REF]. All the models are coded in C++ and solved via CPLEX 12.6.3 on a Mac with Inter Core i5 clocked at 2.7 GHz and with 8 GB of RAM.

As CPLEX 12.6.3 can handle MIQCPs where quadratic terms are multiplications of binary variables, we also use it to solve (P) directly as a benchmark to (P MC ). Results are summarized in Table 3.1. For each problem instance, the size of the physical network (number of servers, number of edges), the number of requests, the number of binary variables, as well as CPU time and the number of branching nodes are reported. As might be expected, the computational time increases exponentially for both formulations as the size of the problem instances increase. Clearly, the simple McCormick relaxation based model outperforms the routines of CPLEX solver for this mapping problem. Let us now focus on the performance of models (P RLT1 ), (P RLT2 ), (P V1 ) and (P V2 ). To show the effectiveness of RLT-based and additional valid inequalities, we first solve the continuous relaxations of different formulations. Then we solve each problem using standard CPLEX 12.6.3 Branch-and-Cut with default settings. Results are summarized in Table 3.2.

Continuous relaxation figures in the table show that for each problem instance, the strength of formulations increases as more valid inequalities are added, which is within our expectation. However notice that the differences w.r.t. the optimality gap between (P RLT1 ) and that of (P RLT2 ) is negligible while the computational cost of solving (P RLT2 ) can be ten times higher than that of solving (P RLT1 ). This is probably because the former has 2 r∈R |V r ||S|(|S| -1) more constraints than the latter. Another observation is that adding (Cut),(Cut') to models (P RLT1 ) and (P RLT2 ) reduces around 40% gap, which confirms the necessity of exploring additional valid inequalities in addition to the RLT based ones.

To balance the strength of (P V1 ) and (P V2 ) and their problem sizes, we put constraints (Cut'), (LC RLT ), (LC RLT ) to the user-cut pool of CPLEX. The B&B results imply that compromises have to be made between the strength of some valid inequalities and their computational costs. For example, both (P RLT2 ) and (P V2 ) are stronger (cf. Columns '#Nodes') than (P MC ), but they are also computationally more costly. On the other hand, model (P RLT1 ) and (P V1 ), that are slightly weaker than model (P RLT2 ) and (P V2 ), are computationally five times as efficient as model (P MC ). For the largest instance, (P RLT1 ) reduces the computation time from 13hs to 3.74hs, and (P V1 ) performs best if we add constraints (Cut') in an iterative way. In addition, the results also suggest two factors that influence the computational performance: the number of virtual requests and the size of the physical network. For a fixed number of requests, the computation time increases gradually as the size of the physical network increases; on the other hand, given a fixed physical network, the computation time rises significantly as the number of requests grows. Generally, the size of virtual requests might affect computational performance heavier than that of a physical network. Overall, the computational experiments indicate that model (P V1 ) performs best among all models.

A Lagrangian decomposition

The previous numerical results show that the gap incurred by some convex relaxation of the mapping programs is typically large. While several proposed linear valid inequalities are quite effective in strengthening the linear relaxation, there are cases where the gap can be 20%. The goal of this section is to use Lagrangian relaxation to decompose the problem thereby solving the relaxed problem efficiently along with lower bounds of good quality. Of-course, there are plenty of ways to decompose the problem and we have tried many. In this section, we present a request based decomposition scheme, which has a nice interpretation in the context of cloud computing. It also guarantees that the resultant lower bound dominates the optimum of the continuous relaxation of formulation (P V1 ) which is numerically the strongest one among the aforementioned models. Numerical experiments suggest that the proposed scheme can computationally be attractive for large scale problems. Let us recall (P V1 ) below.

min k∈S F k θ k + k∈S A k r∈R i∈V r c ri x ri k + e∈E W e φ e (P V1 ) s.t. k∈S x ri k = 1, r ∈ R, i ∈ V r , (AC) p∈S:p =k y rij kp = x ri k , r ∈ R, i, j ∈ L r , k ∈ S, (AC RLT ) i∈V r x ri k ≤ θ k , r ∈ R, k ∈ S, (LC)
r∈R i∈V r c ri x ri k ≤ M k θ k , k ∈ S, (KP) r∈R i∈V r m ri x ri k ≤ C k θ k , k ∈ S, (KP') r∈R i,j∈V r : i =j k,p∈S: k =p,e∈P kp f rij y rij kp ≤ B e φ e , e ∈ E, (QC)
(k,p)∈S 2 :k =p,e∈P kp y rij kp ≤ φ e , e ∈ E, r ∈ R, i, j ∈ V r : f rij > 0, ( Cut
) (i,j)∈L r y rij kp ≤ φ e , e ∈ P kp , k = p ∈ S, r ∈ R, (Cut') y rij kp ≥ 0, r ∈ R, i, j ∈ V r : i = j, k = p ∈ S, (3.1d)
θ k , φ e , x ri k ∈ {0, 1}, r ∈ R, i ∈ V r , k ∈ S, e ∈ E.

A request based decomposition

We first present a decomposition leading to a sequence of subproblems associated with each request, each sever, and each link. To this end, we disaggregate constraints (KP), (KP'), (QC) by reformulating them with a handful of auxiliary variables with corresponding interpretations.

• w r k reserved CPU for request r on server k.

• z r k reserved Memory for request r on server k.

• κ r e reserved bandwidth for request r on link e.

The equivalent counterpart of (KP), (KP'), (QC) is then

i∈Vr c ri x ri k ≤ w r k , r ∈ R, k ∈ S, (3.5
)

i∈Vr m ri x ri k ≤ z r k , r ∈ R, k ∈ S, (3.6)
i,j∈Vr:

i =j k,p∈S: k =p,e∈P kp f rij y rij kp ≤ κ r e , r ∈ R, e ∈ E, (3.7
)

r∈R w r k ≤ C k θ k , k ∈ S, λ ∈ R |S| + (3.8) r∈R z r k ≤ M k θ k , k ∈ S, µ ∈ R |S| + (3.9) r∈R κ r e ≤ B e φ e , e ∈ E. σ ∈ R |E| + (3.10)
To make the problem separable by request while ensuring the strong lower bounds, we proceed the following. First, we copy variables θ and φ by introducing following constraints:

θ r k = θ k , r ∈ R, k ∈ S, η ∈ R |R|×|S| (3.11) φ r e = φ e , r ∈ R, e ∈ E. ζ ∈ R |R|×|E| (3.12)
Then, we replace θ k with θ r k in constraints (LC). Likewise, we replace φ e with φ r e in constraints (QC). Finally, we denote the resulting reformulation of (P V1 ) as Proof: It holds that the Proj (x,y,θ,φ) F(P v1 ) ⊆ F(P V1 ). The reverse also holds by letting

w r k = i∈Vr c ri x ri k ≤ w r k , z r k = i∈Vr m ri x ri k , κ r e = i,j∈Vr: i =j k,p∈S: k =p,e∈P kp f rij y rij kp for each r ∈ R, k ∈ S, e ∈ E.
Relaxing the five constraints with associated Lagrangian multipliers λ, µ, σ, η, ζ leads to the Lagrangian over variables (λ, µ, σ, η, ζ; w, z, κ, θ, x, φ). For ease of notation, let

v * = (λ, µ, σ, η, ζ), v = (w, z, κ, θ, x, φ). The Lagrangian is then L (v * , v) = r∈R k∈S A k i∈Vr c ri x ri k + w r k λ k + µ k z r k + θ r k η r k + e∈E (σ e κ r e + ζ r e φ r e ) + k F k -C k λ k -M k µ k - r∈R η r k θ k + e W e -B e σ e - r∈R ζ r e φ e .
The Lagrangian problem reduces to |R| + |S| + |E| subproblems, i.e., for each r ∈ R, we obtain L r (v * ) by solving

min k∈S A k i∈Vr x ri k a ri + λ k w r k + µ k z r k + θ r k η r k + e∈E (σ e κ r e + ζ r e φ r e ) (Sub r ) s.t. k∈S x ri k = 1, i ∈ V r , p∈S:p =k y rij kp = x ri k , i, j ∈ L r , k ∈ S, i∈V r x ri k ≤ θ r k , k ∈ S, i∈Vr c ri x ri k ≤ w r k , k ∈ S, i∈Vr m ri x ri k ≤ z r k , k ∈ S, i,j∈Vr: i =j k,p∈S: k =p,e∈P kp f rij y rij kp ≤ κ r e , e ∈ E, (k,p)∈S 2 :k =p,e∈P kp y rij kp ≤ φ r e , e ∈ E, i, j ∈ V r : f rij > 0, (i,j)∈L r y rij kp ≤ φ r e , e ∈ P kp , k = p ∈ S, y rij kp ≥ 0, i, j ∈ V r : i = j, k = p ∈ S, θ r k , φ r e , x ri k ∈ {0, 1}, i ∈ V r , k ∈ S, e ∈ E.
Note that one can even relax the integrality of θ r and φ r to reduce the complexity of solving a subproblem associated with each request. For each k ∈ S, and e ∈ E, we solve the following problems analytically by checking the coefficients.

L k (λ, µ, η) = min θ k ∈{0,1} F k -C k λ k -M k µ k - r∈R η r k θ k , L e (σ, ζ) = min φe∈{0,1} W e -B e σ e - r∈R ζ r e φ e .
Consequently, the dual objective denoted by Ψ(v * ) is

Ψ(v * ) = r∈R L r (v * ) + k∈S L k (λ, µ, η) + e∈E L e (σ, ζ).
And the dual problem is then

max v * Ψ(v * ). (3.13)
Before the solution procedure of (3.13), we verify that evaluating Ψ(v * ) is relatively easier than solving (P V1 ) itself. In particular, when the number of requests is large such that CPLEX cannot handle (due to the increasing number of binary variables), this decomposition scheme might be attractive. In addition, the lower bound achieved by this decomposition is stronger than the bound obtained by solving the continuous relaxation of (P V1 ). We formalize this result below.

Theorem 4 The optimum of (3.13) dominates that of continuous relaxation of (P V1 ).

Proof: According to Theorem 2, the optimum of (3.13) is equivalent to the optimum of minimization of the objective function of (P V1 ) over S, where × r∈R S r represents the Cartesian product of (S r ) r∈R , where S r represent the constraint set of each subproblem (Sub r ). It also holds that conv{ × r∈R S r } = × conv{S r }. More- over, it holds that S is "smaller" than the continuous relaxation of F(P v1 ), which completes the proof in conjunction with Lemma 3.2.

S = conv{ × r∈R S r } {(v,
Theorem 4 also indicates a hierarchy of request based Lagrangian relaxations. We may consider a sequence of covers of set {S r } r∈R . For example, we can partition set {S r } r∈R into disjoint subsets and each subset is associated with two requests. The resulting Lagrangian relaxation bound will be larger than the one based on single request decomposition. When the size of each subset is |R|, we get problem (P v1 ).

The solution procedure

In this section, we develop a solution procedure of (3.13). Ben-Ameur and Neto [START_REF] Ben-Ameur | Acceleration of cutting-plane and column generation algorithms: Applications to network design[END_REF] proposed the so-called in-out algorithm for convex programs, which is treated as an acceleration of cutting-plane and column generation algorithms. We employ this algorithm to solve a Lagrangian dual problem. Let us recall several features of problem (3.13).

• The objective function Ψ(v * ) is additively composed of functions L r , L k , L e and all functions are concave functions in variables v * .

• For each input v * , we have evaluations of functions L r , L k , L e and access to their subgradients at v * .

Therefore, we are able to construct the following master problem by approximating each concave function with a sequence of subgradients and evaluations.

max r∈R γ r + k∈S γ k + e∈E γ e (M) s.t. γ r ≤ L r (v * ) + k∈S w rn k (λ k -λ n k ) + k∈S z rn k (µ k -µ n k ) + k∈S θ rn k (η r k -η rn k ) + e∈E κ rn e (σ e -σ n e ) + e∈E φ rn e (ζ r e -ζ rn e ), r ∈ R, n = 1, . . . , N, γ k ≤ L k (λ, µ, η) -α k θ n k (λ k -λ n k ) -β k θ n k (µ k -µ n k ) -θ n k r∈R (η r k -η rn k ), k ∈ S,
γ k + e∈E γ e ≤ U λ k ,µ k , σ e ≥ 0, r ∈ R, k ∈ S, e ∈ E.
where U is an upper bound of Ψ(v * ). A trivial choice of U can be an upper bound of the primal cost function, e.g., k∈S

F k + e∈E W e φ e + k∈S A k r∈R max i∈V r c ri .
To use the in-out algorithm, we need to choose a feasible point of the dual problem (3.13). A key observation is that any v * defined before is feasible to (3.13). For each v * , let us denote by

L(v * ) ∈ R |R|+|S|+|E| the vector ((L r (v * )) r∈R , (L k (λ, µ, η)) k∈S , (L e (σ, ζ)) e∈E ).
Correspondingly, any (v * , L(v * )) is a feasible point to the master problem (M). The in-out algorithm usually outperforms the standard cutting-plane (or column generation) algorithms since it allows relatively flexible choose of the subsequent evaluation points. Following Ben-Ameur and Neto [START_REF] Ben-Ameur | Acceleration of cutting-plane and column generation algorithms: Applications to network design[END_REF], let us call the evaluation point separation point, the optimal solution of (M) out-point, and call the feasible point in-point. The separation point is chosen as a convex combination of the in-point and out-point. Next, we state the in-out Algorithm when applied to problem (3.13).

The choice of the initial in-point: Even though each subproblem (Sub r ) (a MILP) is easier than problem (P V1 ), it is still computationally costly. This motivates us to reduce the number of iterations of Algorithm 1. To this end, we may need to select a good inpoint with some cost. Specifically, we first solve the continuous relaxation of (P v1 ) to optimality and take the dual multipliers corresponding to constraints (3.8), (3.9), (3.10), (3.11) and (3.12) respectively as the initial values of v * in . Then we solve the subproblems with v * in ∈ R |R|+|S|+|E| , which outputs γ in . In the sequel, we shall denote by v * 0 in and γ 0 in the respective initial values of v * in and γ in . We show that at the initialization procedure, we obtain a lower bound that dominates the continuous relaxation of (P V1 ).

Lemma 3.3 The value of r∈R γ r0 in + k∈S γ k0 in + e∈E γ e0
in dominates the optimal objective value of the continuous relaxation of (P V1 ) Proof: First, by Lemma 3.2 and strong duality, the optimum of the continuous relaxation of (P V1 ) is equal to the minimization of L(v r0 in , •) over continuous relaxation of (P v1 ). Second, observe that all the integrality constraints are not relaxed in subproblems associated with L r , L k , L e . Combining the above observations completes the proof.

(v * sep , γ sep ) = α(v * out , γ out ) + (1 -α)(v * in , γ in ) evaluate function L r , L k , L e with input v * sep if r∈R L r (v * sep ) + k∈S L k (v * sep ) + e∈E L e (v * sep ) < r∈R γ r sep + k∈S γ k sep + e∈E γ e sep then add cuts at point v * sep by concavity if r∈R L r (v * sep ) + k∈S L k (v * sep ) + e∈E L e (v * sep ) > r∈R γ r in + k∈S γ k in + e∈E γ e in then (v * in , γ in ) ←-(v * sep , L(v * sep )) end solve (M) else end (v * in , γ in ) ←-(v * sep , L(v * sep )
) end Remark 3.1 Our numerical results show that the lower bound at the initialization procedure, i.e., r∈R

γ r0 in + k∈S γ k0 in + e∈E γ e0
in , is generally close to the dual optimum, i.e, the optimal objective value of (3.13). Therefore, for computational benefits, one can terminate Algorithm 1 at the early stage of the algorithms or even after the initialization procedure.

Numerical experiments

In this section, we carry out some numerical results to evaluate the effectiveness of the proposed decomposition scheme. All test instances have a moderate number of small request (i.e., each request has 5 VMs). Model (P V1 ) is used as the benchmark to the request based decomposition. Algorithm 1 is implemented in C++. Numerical results are summarized in Table 3.3.

For each problem instance, we solve model (P V1 ) exactly and solve Lagrangian dual problem (3.13) by Algorithm 1. MILPs are solved by CPLEX12.63. To keep the size of model (P V1 ) moderate, we put valid inequalities of (Cut') to the user-cut pool of CPLEX. Optimality gap is quantified by gap = Exact-v Exact , where v is replaced by the optimum of the continuous relaxation of (P V1 ) and the lower bound obtained from Algorithm (1) respectively. Another observation is that as the number of requests increases, the strength of formulation (P V1 ) and that of Lagrangian decomposition scheme appear weaker. On the one hand, it implies that the convex relaxation of (P V1 ) is weaker for instances with relatively larger number of requests. On the other hand, it reveals the fact that the relaxed coupling constraints (3.8) and (3.9), (3.10) become more important in (P V1 ) as the number of requests goes up. For problem instances having more than 6 requests, CPLEX cannot solve the to optimality within the 10 hours time limit. We use the proposed hierarchy in Section 3.4.1 to generate a sequence of strong bounds in a reasonable time. More precisely, we consider two levels of Lagrangian decomposition. Level 1 refers to the one where the problem (P V1 ) is decomposed by each single request, each server, and each link. Level 2 refers to the one where we decompose the problem by each pair of requests, each server, and each link.

Results in

Due to the lack of exact optimal values, we measure the strength of the above two decomposition schemes by the closed gap defined below closed gap = Lagrangian bound -continuous relaxation value upper bound -continuous relaxation value where the upper bound is the one provided by CPLEX after 10 hours branch-and-cut. The numerical results are summarized in Table 3.3. We highlight the following observations. First, as might be expected, the bound provided by Level 2 decomposition is stronger than the one obtained by Level 1. Second, it appears that as the number of requests increases, Level 2 decomposition does not perform significantly better than the Level 1 decomposition. Finally, it might be interesting to integrate the request based Lagrangian decomposition scheme into the branch and cut procedure so that the exact solution procedure could be further accelerated.

A symmetry-induced model

In this section, we reformulate the initial model (P) by exploiting the symmetric structure of the virtual requests. To this end, we introduce the definition below. Definition 3.1 For each request r, a virtual machine i is equivalent to a virtual machine j if and only if the following conditions hold.

• The frequencies of CPU resource of VM i and j are equal.

• The capacity of memory resource of VM i and j are equal.

• The throughput requirements of logical links between VM i and VM j are equal, i.e., f rij = f rji .

• For each VM l where l = i, j, the throughput requirements of logical links (i, l) and (j, l) are same, i.e., f ril = f rjl and f rli = f rlj .

It is clear that the equivalence is well defined as it is a reflective, symmetric, transitive relation. This definition allows a partition of a request (a set of virtual machines) V r . Each member of a partition consists of equivalent VMs. We may interpret a member of a partition as a VM type. A type t is identified by the frequency of CPU c t , the memory capacity m t , and the throughput f t . We emphasize here that the type we defined here is also associated with the throughput, whereas in practice, the term "type" is generally characterized by CPU, memory (see section 2.4). For the sake of further development, we introduce additional parameters and variables below:

• T r : set of VM types required in virtual request r ∈ R.

• V rt : set of VMs of type t ∈ T r in request r, i.e., V r = t∈T r V rt .

• x rt k : 1 if we mapping a type t VM of request r to server k, 0 otherwise.

It is clear that for each request, feasible solutions are equivalent among decision variables associated with same VM type under permutation in the sense that their corresponding objective values are equal. As will be shown below, this simple feature allows us to deduce an equivalent model with fewer binary variables and constraints. We may also exploit symmetries among edges to reduce the number of (QC). It is obvious that in the context of formulation (P), two edges e and e are equivalent if and only if W e = W e , B e = B e , and they are related to the same set of servers in terms of shortest path routing. As a result, redundant edges do not need to be considered.

Consequently, the objective function becomes

k∈S F k θ k + k∈S A k r∈R t∈T r c t x rt k + e∈E W e φ e . (3.14)
The assignment constraint (AC) is rewritten as

k∈S x rt k = V rt , r ∈ R, t ∈ T r . (3.15)
The location constraint (LC) is replaced with

t∈T r x rt k ≤ θ k , r ∈ R, k ∈ S. (3.16)
The knapsack constraint associated with CPU and the one associated with memory constraints can be rewritten as

r∈R t∈T r c t x rt k ≤ C k θ k , k ∈ S, (3.17
)

r∈R t∈T r m t x rt k ≤ M k θ k , k ∈ S. (3.18)
The bilinear type knapsack constraint can then be replaced by Proposition 1 Formulation (SP) is equivalent to (P) in a sense that for any feasible solution of (P), there exits a feasible solution of (SP) such that their objective values are equal.

As stated before, bilinear terms and binary variables make (P) and (SP) hard to solve. We compare the sizes of these two problems in terms of the numbers of binary variables and constraints. Notice that the number of variables associated with linear constraints and the number of bilinear constraints in both formulations are same. So we just report in Table 3.5 the number of bilinear variables and bilinear terms and linear constrains.

It is obvious that if |T r | < |V r |, formulation (SP) is probably attractive. In general, Table 3.5: Problem sizes of (P) and (SP)

Model Bilinear variables Linear constraints Bilinear terms

(P) r∈R |V r ||S| r∈R |V r | + (2 + |R|)|S| |S||V r |(|S| -1)(|V r | -1)/2 (SP) r∈R |T r ||S| r∈R |T r | + (2 + |R|)|S| |S||T r | 2 (|S| -1)/2
this condition holds in the mapping problem. First, in a virtualized environment, the combination of CPU and memory for virtual machines is designed in an optimal way such that the computing resource of servers can be allocated efficiently. Second, it is really tough for clients to quantify exact the requirements of CPU, memory, and throughputs. Thus cloud operators usually design several VM types. For example, Amazon Elastic Cloud Computing (EC2) mainly provides "small" VMs configured with 1.7 GB memory and 1 compute unit [START_REF] Wang | The impact of virtualization on network performance of amazon ec2 data center[END_REF]. Third, even if throughput between VMs are not strictly categorized, there exists strong correlations between VM types and throughputs among them, i.e., the "larger" (regarding CPU and memory) two VMs are, the more throughputs will be needed in-between.

The partition of a request can be done before the solution procedure. First, we partition the request graph by characteristics of CPU and memory. Then for each subset of the partition, we further partition it based on the throughput values.

To conclude this section, we remark that one may employ all the relaxation techniques proposed previously for (P) apply to model (SP) directly for resolution. Yet, implications among constraints need to be examined carefully. For example, if one introduce RLT constraints based on assignment constraints (3.15), the McCormick constraints cannot be eliminated as (P RLT1 ). Strong inequalities of the form (Cut) cannot be employed either. In brief, we will use the following formulation to solve (SP).

SP V1 : {(SP)} ∩ {(AC RLT ), (3.1a) -(3.1d), (Cut')} . (SP V1 )
Note also that the strength of formulation (SP V1 ) is usually weaker than the general formulation (P V1 ).

A heuristic

For model (SP), it is clear that when |T r | = 1, the benefits can be significant and when |T r | = |V r |, the two models are equivalent. However, for problem instances not satisfying the equivalence conditions defined in 3.1, the advantages of model (SP) can hardly manifest. To address this issue, we propose a heuristic to obtain a feasible solution with an upper bound.

Assume that |T r | is large and we want to find a feasible solution of (P) in a reasonable time. Clearly the VMs can be easily categorized regarding the requirements of CPU and memory, thus we elaborate on partitioning VMs by throughput values. Intuitively, VMs with similar f rij will be grouped in the same set. Meanwhile, we want to minimize the amount of throughput among two subsets. To this aim, we consider minimizing the normalized cut (that minimizes the multi-cuts while balancing sizes of each subgraph) for a virtual request graph G r . The problem itself is NP-hard, but we might as well to solve it approximately usng the well-known shi-malik algorithm [START_REF] Shi | Normalized cuts and image segmentation[END_REF]. To be specific, we proceed the following. For VMs in the same subset, every throughput value is augmented to the maximum; for VMs among two different subsets, each mutual throughput value is shifted to the largest one. This setting allows us to make use of model SP, which outputs a feasible solution and an upper bound. Intuitively, the upper becomes stronger as the number of subsets decreases. When the number of subsets is equal to the number of VMs in G r , we get the global optimum of the mapping problem.

Numerical experiments

In this section, we test the effectiveness of model (SP) numerically. We choose to test on problem instances of a single request with 10, 20, 30, 40 VMs respectively. Model (SP) and the aforementioned heuristic coded in C++.

We carry out experiments on two sets of problem instances. For all problem instances, VMs are categorized by CPU and Memory. In addition, for the first set of problem instances, throughput among a same category of VMs are assumed to be the same, whereas for the second set of problem instances, throughput values are usually different.

Numerical results of the first problem instance set are summarized in Table 3.6. Note that the optimum of (SP V1 ) is exactly the optimum of (P V1 ). As expected, the symmetryinduced model (SP V1 ) performs significantly better than (P V1 ) for all the tested problem instances. This is mainly due to the reduction of binary variables and bilinear terms. However, as the size of the physical graph increases, the solution time of (SP V1 ) also increases exponentially. For example, for 10 VMs, it takes around 400 seconds for instance [START_REF] Ben-Ameur | Acceleration of cutting-plane and column generation algorithms: Applications to network design[END_REF][START_REF] Kall | Stochastic programming[END_REF] while it just takes 1.22 seconds for instance [START_REF] Anstreicher | Computable representations for convex hulls of low-dimensional quadratic forms[END_REF][START_REF] Ben-Tal | Adjustable robust solutions of uncertain linear programs[END_REF]. Now let us focus on the second problem instance set. Following the observation from [START_REF] Meng | Improving the scalability of data center networks with traffic-aware virtual machine placement[END_REF] that 80% of throughput values are quite low and 4% of them have a 10 times higher rate, we generate throughput values with a probability 0.8 being lower than 2 Mbps/min and with a probability 0.2 lying between 10 Mbps/min and 30 Mbps/min. We employ the heuristic proposed in Section 3.5.1 to partition virtual request G r to N subsets. Results are summarized in Table 3 . Exact values are computed using (P V1 ). We stress two observations:

• For each problem instance, as might be expected, the quality of upper bounds become better as the size of the partition increases. • As the size of the problem instance increases, the gap becomes larger and the influence of the partition size N becomes more significant. To conclude this section, we remark that there is still a large number of symmetries in model (SP V1 ). For example, in a physical network, two servers k, p might be equivalent in the sense that (θ k , θ p ) = (0, 1) or (θ k , θ p ) = (1, 0) leads to the same mapping cost. One may introduce symmetry breaking constraints statically to the model or dynamically during a branching-and-cut procedure.

A perspective

In practice, the arrival of VMs is dynamic and unpredictable. The static mapping model (P) can hardly work in that case. In this section, we propose a model that optimizes certain mapping policies as a perspective.

In a cloud computing environment, requests might be categorized. For example, let us consider the cloud service of Platform as a Service (PaaS), where platforms might be Linux, Windows, Unix, etc.. Let set R be set of request types, and let x r denote the mapping polices for requests of type r. Specifically, we may interpret x ri k as the probability that VM i will be assigned to sever k. The model can be developed simply by relaxing the integrality constraints of x in model (P). The objective function turns out to be the expectation of future cost. When requests of type r arrive, cloud provider will map x ri k percent of VM i to server k. Note that unlike a MILP or a convex MIQCP, relaxing integrality constraints of x does not lead to an easier problem. In fact, it is computationally more challenging than problem (P). For problem (P), one can get the optimal solution and objective value using standard finite branch and bound algorithms. For this model where continuous variables are involved in bilinear constraints (QC), one may need to consider spatial branch and bound framework for two reasons. First, the optimal solution to a given relaxation is not necessarily feasible. Second, the convex underestimators of bilinear constraints may need to be updated in each subdivision of the search region. Many related fundamental problems, such as, convex and concave envelopes of a general bilinear function over a hypercube need to be investigated.

Conclusion

We proposed a model for the optimal mapping of virtual machines in a cloud environment, by taking into account link capacity constraints, location constraints, and knapsack constraints regarding CPU and memory. We then strengthened the model by employing different levels of RLT and adding strong valid inequalities. Computational results confirm the power of valid inequalities, especially those that exploit the problem structure. To handle relatively larger problem instances, we proposed a request based Lagrangian decomposition scheme. We show, both theoretically and numerically, that this decomposition can produce bounds greatly stronger than continuous relaxation bounds. In addition, a novel model exploiting symmetries of a request was discussed and analyzed. Finally, we pointed out that it would be interesting to associate each request with certain mapping policy and investigate related techniques to solve the resulting model.

Appendix: numerical comparison of model (P MC ) and model (P GL )

We provide numerical results on some instances to compare the differences of two different models (P MC ) and (P GL ). Convex and concave envelopes of bilinear functions

Summary

Convex and concave envelopes of nonconvex functions play an important role in global optimization. We discuss in this chapter a couple of approaches to approximating the convex and concave envelopes of bilinear functions, especially over hypercubes due to its generality. The first approach is based on a semidefinite reformulation. The second approach considers some predefined set covers of a hypercube and leads to a linear program.

Then we establish a connection between the convex envelope of a bilinear function and the concave envelope of a polyhedral function. Numerical experiments are conducted to compare the two approaches. As a possible extension, a novel approach to approximate the envelopes is discussed and illustrated.

Introduction

There has been a great amount of studies devoted to developing convex underestimators and concave overestimators of nonlinear functions f (x) over various polyhedral sets. One of the most important motivations in optimization community is that, one can reformulate a complicated usually non-convex problem as an easier problem with convex representations of the objective function and constraints. Such convex reformulations can then be solved repeatedly in branch-and-bound algorithms, where searching spaces are refined in a convergent way. In general, the computational efficiency of branch-and-bound is greatly influenced by the strength of convex reformulations.

Among such investigations, the construction and approximations of convex and concave envelopes of bilinear functions over boxes draws much attention. This is mainly due to the following reasons. First, many important problems, for example, quadratic assignment problems [START_REF] Hahn | Lower bounds for the quadratic assignment problem based upon a dual formulation[END_REF][START_REF] Wang | Optimal mapping of cloud virtual machines[END_REF], pooling and blend [START_REF] Foulds | A bilinear approach to the pooling problem †[END_REF][START_REF] Audet | Pooling problem: Alternate formulations and solution methods[END_REF][START_REF] Misener | Apogee: Global optimization of standard, generalized, and extended pooling problems via linear and logarithmic partitioning schemes[END_REF] involve bilinear terms. Second, boxes are naturally outer approximations of polytopes and therefore the estimators of convex and concave envelopes of f over a box are also valid estimators for function f over polytopes belonging to the box. In addition, most branch-and-bound algorithms partition the searching space by divisions of the feasible region into boxes, which also motivates investigations in deriving strong convex and concave estimators over boxes. It is well known that a n-dimensional box is simply linear transformations of a hypercube. Therefore, we consider in this chapter the convex (resp. concave) envelope of a bilinear function (4.1) over a hypercube H n = [0, 1] n :

f (x) = i<j Q ij x i x j (4.1)
where

Q ij ∈ R, 1 ≤ i < j ≤ n.
There are two main research directions on the construction of envelopes of f over H n . The first direction is motivated by the definition that the convex envelope of any bounded f over any non-empty compact set P is a function

q f : P x → q f (x) = inf {u : (x, u) ∈ conv epi f } ,
where epi f = {(x, u) : f (x) ≤ u} and conv of a set represents its convex hull. Geometrically, one can interpret the convex envelope of f as the bottom of the convex hull of the epigraph of f over P . Since P is compact, the definition is equivalent to

min λ i ⎧ ⎨ ⎩ v i ∈P λ i f (v i ) : λ i v i = x, v i ∈P λ i = 1, λ i ≥ 0 ⎫ ⎬ ⎭ . ( 4.2)
The concave envelope f (x) of f over P is defined symmetrically by replacing minimization to maximization. When P = H n , a nice property of bilinear function f due to Sherali [START_REF] Sherali | Convex envelopes of multilinear functions over a unit hypercube and over special discrete sets[END_REF] and Rikun [START_REF] Rikun | A convex envelope formula for multilinear functions[END_REF] is that its convex and concave envelopes over a hypercube H n are vertex polyhedral, i.e., the envelopes of f on H n coincide with the envelopes of its restriction to the vertices of H n . This property allows us to simplify the definition (4.2) to the following,

q f (x) = min λ i 2 n i=1 λ i f (v i ) : 2 n i=1 λ i v i = x, 2 n i=1 λ i = 1, λ i ≥ 0, i = 1, . . . , 2 n , ( 4.3)
where v i , i = 1, . . . , 2 n are vertices of H n . However, as the number of extreme points of a hypercube is exponential in the dimension, it is known that a full characterization of such polyhedral envelopes is in general difficult. The characterizations of convex and concave envelopes in lower dimensional hypercubes have been studied in [START_REF] Meyer | Convex envelopes for edge-concave functions[END_REF]. The second research line can be treated as the dual of first research direction since the convex (resp. concave) envelope is interpreted as the pointwise supreme (resp. infimum) of an affine under-(resp. over-) estimator of f over P . For example, the convex envelope is

q f (x) = max α,γ γ : α T (x -v) + γ ≤ f (v), ∀ v ∈ P , (4.4)
where (α, γ) ∈ R n+1 defines the supporting hyperplane of the q f at point x. When P = H n , one can replace set P with its extreme points due to the vertex polyhedral property of the convex envelope of f . However, the number of extreme points is exponential. A cutting-plane algorithm was proposed in [START_REF] Bao | Multiterm polyhedral relaxations for nonconvex, quadratically constrained quadratic programs[END_REF] to find a facet of the convex envelope by separating the supporting function (α, γ) ∈ R n+1 at a pre-defined point. One of the most frequently used relaxations to approximate the convex and concave envelopes of f (x) over a hypercube due to McCormick [START_REF] Mccormick | Computability of global solutions to factorable nonconvex programs: Part i-convex underestimating problems[END_REF] is the following inequalities max{x i + x j -1, 0} ≤ x i x j ≤ min{x i , x j }, i < j.

(4.5)

We will refer to these inequalities as McCormick inequalities or McCormick relaxation.

McCormick also pointed out that (4.5) characterize the exact convex and concave envelope over rectangles. One can also find detailed proof in [START_REF] Al-Khayyal | Jointly Constrained Biconvex Programming[END_REF]. However, for higher dimensional (n ≥ 3) cases, term-wise McCormick relaxation, in general, do not provide a characterization of the convex and concave envelope over a hypercube. The triangle inequalities proposed by Padberg [START_REF] Padberg | The boolean quadric polytope: some characteristics, facets and relatives[END_REF] in concert with McCormick inequalities provide an explicit characterization of convex and concave envelopes of f (x) over a cube (n = 3). As noticed by [START_REF] Meyer | Convex envelopes for edge-concave functions[END_REF], the number of facets defining the convex and concave envelopes is indicated by the number of simplices induced by the triangulations of a hypercube and unfortunately, the number grows super-exponentially in the number of dimension. Locatelli and Schoen [START_REF] Locatelli | On convex envelopes for bivariate functions over polytopes[END_REF] utilized the definition (4.4) and derived explicit convex envelopes for some general bivariate functions over various polytopes in R 2 .

A number of papers are also dedicated to giving explicit characterizations of convex and concave envelopes over special polytopes. For example, Locatelli and Schoen [START_REF] Locatelli | On convex envelopes for bivariate functions over polytopes[END_REF] used the definition (4.4) and derived explicit convex envelopes for some general bivariate functions over various polytopes in R 2 . Sherali and Alameddine derived the convex envelope of a bivariate function over the so-called D-polytope. Linderoth [START_REF] Linderoth | A simplicial branch-and-bound algorithm for solving quadratically constrained quadratic programs[END_REF] analytically derived the explicit functions for the disjoint triangular regions of a rectangle via disjunction x ≤ y or x ≥ y. Tawarmalani et al. [START_REF] Tawarmalani | Explicit convex and concave envelopes through polyhedral subdivisions[END_REF] derived convex and concave envelopes for several nonlinear functions via detecting certain polyhedral divisions of a hyper-rectangle.

To deal with the difficulty in (4.3) and (4.4), one natural approach is to consider extended formulations with moderate sizes thereby obtaining strong convex and concave estimators. However, as observed by Luedtke et al. [START_REF] Luedtke | Some results on the strength of relaxations of multilinear functions[END_REF], extended formulations do not always bring improvements compared with McCormick relaxations (4.5). A comprehensive study on the strength of McCormick relaxations of multilinear functions over a hypercube is provided in [START_REF] Luedtke | Some results on the strength of relaxations of multilinear functions[END_REF]. In particular, it is shown that if the coefficients of a bilinear function f (x) defined as (4.1) are nonnegative, i.e., Q i,j ≥ 0, ∀ i < j, the ratio of the gap between McCormick's overestimator and underestimator constructed via (4.5) over the true gap of the convex and concave envelopes is always less than 2. As proposed by Lasserre and Thanh [START_REF] Lasserre | Convex underestimators of polynomials[END_REF], one can also construct a convex polynomial p d ∈ R[x] with degree no greater than d that underestimates f . Based on Putinar's Positivstellensatz [START_REF] Putinar | Positive polynomials on compact semi-algebraic sets[END_REF], p d can be approximated as close as desired by solving a hierarchy of semidefinite programs.

We investigate in this chapter the strength of different approaches on approximating the convex and concave envelopes of the bilinear f over H n . The outline and contributions are summarized below.

Outline and contributions. In Section 4.3, we firstly review a semidefinite program (SDP) based estimator, which is based on the intersection of a SDP constraint and Mc-Cormick inequalities. This work has been proposed by Anstreicher [START_REF] Anstreicher | On convex relaxations for quadratically constrained quadratic programming[END_REF]. Our contribution lies in some theoretical results on two variants of formulations on bilinear functions as a special case of general quadratic functions over a box. As a consequence, an alternative proof to a result from Burer and Letchford [START_REF] Burer | On nonconvex quadratic programming with box constraints[END_REF] is provided. Then we propose a novel approach by considering a predefined tool to approximate the convex and concave envelope of f over H n . This approach leads to a linear reformulation and we are able to establish connections with various inequalities and the reformulation-linearization-technique (RLT) method proposed in [START_REF] Adams | A hierarchy of relaxations leading to the convex hull representation for general discrete optimization problems[END_REF]. Finally, we show that the convex (resp. concave) envelope of a bilinear function f is affinely equivalent to the concave (resp. convex) envelope of a polyhedral function. In Section 4.4, we report a set of numerical results on a number of problem instances to compare these two approaches. Following the numerical results, we propose an interesting research direction on approximating the convex envelope valid over a hypercube when all coefficients of f are nonnegative. Concluding remarks follow in Section 4.6.

Notation. Throughout this chapter, for any matrix C, C T denotes its transpose and C ij represents its i th row j th column component. C I,J represents the submatrix composed by set I of columns of C and set J of rows of C. When I = J, we simply use C I in place of C I,I . For a set S, we use conv S to represent its convex hull. If S is finite, we use |S| to represent its cardinality. We use H n to represent an n-dimensional hypercube, i.e., H n = [0, 1] n . For a general function f : H n → R, we denote by q f (resp. f ) the convex (resp. concave) envelope of f over an n-dimensional hypercube. The inner product between two matrices A, B ∈ R m×n are denoted by A, B . For a square matrix A ∈ R n×n , diag(A) ∈ R n represents the vector of diagonal elements. Conversely, for a vector a ∈ R n , Diag(a) ∈ S n stands for the n × n diagonal matrix filled by vector a. For a symmetric matrix X ∈ S n and vector x ∈ R n , we use S(x, X) to represent the matrix 1 x T x X in S n+1 . We use vector τ (X) = (X ij ) 1≤i≤j≤n to represent a vector that is the upper triangular part of the matrix X and τ + (X) = (X ij ) 1≤i<j≤n to represent the strictly upper triangular part.

Main results

In this section, we first review a SDP based convex (resp. concave) underestimator (resp. overestimator) and establish some connections between two formulations. Then, we propose a novel way to construct estimators based on a pre-defined cover set. Finally, we conclude this section by making a connection with the convex and concave envelopes of a polyhedral function, which implies some results in the literature.

SDP based estimators

Given a vector x, we are interested in constructing a convex under-estimator and a concave over-estimator that are valid over a hypercube. Towards this aim, we consider the following assumption.

Assumption 2 A set of valid inequalities involving quadratic terms is available. For each k ∈ K,

x T A k x + b T k x + d k ≤ 0, ∀x ∈ H n , (4.6)
where K is an index set.

Note that Assumption 2 in general holds for bilinear optimization problems over any compact set. For example, such valid inequalities can be McCormick's inequalities (4.5), Padberg's triangle inequalities [START_REF] Padberg | The boolean quadric polytope: some characteristics, facets and relatives[END_REF], RLT based inequalities [START_REF] Adams | A hierarchy of relaxations leading to the convex hull representation for general discrete optimization problems[END_REF]. The principal idea is to use these valid inequalities to construct strong estimators. Let us start by constructing the convex underestimator. We associate nonnegative multipliers α ∈ R |K| + with the set of inequalities. Then we add the negative terms

k∈K α k (x T A k x + b T k x + d k ) to f (x).
For ease of presentation, we define

f α (x) = x T Q(α)x + c(α) T x + p(α), where Q(α) = 1 2 Q + k∈K α k A k , c(α) = k∈K α k b k and p(α) = k∈K α k d k .
Recall that f α is a convex underestimator of f over H n if and only if f α ≤ f and f α is convex over the hypercube H n . Obviously, it holds that f α (x) ≤ f (x), ∀x ∈ H n . The convexity of f α can be ensured by restricting Q(α) to be positive semidefinite. Let us denote by A the set of feasible choices of α:

A := α ∈ R n + : Q(α) ∈ S n + .
We ensure the non-emptiness of set A by explicitly adding valid inequalities x 2 i ≤ x i (i = 1, . . . , n) to inequality system (4.6). As a consequence, it also follows that A has a nonempty interior. The strongest convex underestimator of the form f α can be obtained by solving the following problem sup α∈A f α (x).

(4.7)

Problem (4.7) is then strictly feasible. The dual problem of (4.7) reads inf

X 1 2 Q, X + xx T s.t. A k , X + xx T + b T k x + d k ≤ 0, ∀k ∈ K, X ∈ S n + . (4.8)
The feasible region of (4.8) is nonempty (e.g., take X as the null matrix of size n). In addition, as (4.7) is strictly feasible, the dual optimum is attained and is equal to the primal optimum of (4.7).

Remark 4.1 Note that in model (4.7), we just exploit the information of a set of valid inequalities in K and we restrict Q(α) ∈ S n + . Other than using more valid inequalities, we may replace the SDP matrix cone with "larger" convex cones in geometric sense e.g, copositive cone, in the hope that the underestimator values will be increased. From its dual model (4.8), we see that if we replace the SDP matrix cone with a "smaller" convex cone, e.g., completely positive matrix cone, the underestimator values will be higher. However, checking the membership of copositive cone is co-NP-complete [START_REF] Murty | Some np-complete problems in quadratic and nonlinear programming[END_REF] and checking if a matrix is completely positive cone is NP-hard [START_REF] Dickinson | On the computational complexity of membership problems for the completely positive cone and its dual[END_REF]. One can approximate the copositive matrix cone by doubly nonnegative matrix cone, i.e., the intersection of SDP cone and componentwise nonnegative symmetric matrix cone. For more details, we refer readers to [START_REF] Anstreicher | On convex relaxations for quadratically constrained quadratic programming[END_REF][START_REF] Hiriart-Urruty | A variational approach to copositive matrices[END_REF].

By change of variables

X = X + xx T , one can rewrite (4.8) as min X 1 2 Q, X : (x, X ) ∈ K , ( 4.9)
where

K := (x, X ) ∈ H n × S n : A k , X + b T k x + d k ≤ 0, ∀k ∈ K, S(x, X ) 0 .
Notice that when x is not fixed, problem (4.9) is exactly the Shor's relaxation in conjunction with valid inequalities. Formally, we summarize the foregoing as follows.

Theorem 5 The strength of the SDP relaxation min (x,X )∈K 1 2 Q, X is equivalent to pointwise supremum of quadratic convexifications of quadratic functions defined by (4.7).

Theorem 5 clearly implies that any quadratic convexification based underestimator is dominated by SDP-based estimator. This result in fact coincides with Theorem 2 in [START_REF] Anstreicher | On convex relaxations for quadratically constrained quadratic programming[END_REF], where they claim that the well-known αBB method (c.f., [START_REF] Adjiman | A global optimization method, αbb, for general twice-differentiable constrained nlps -i. theoretical advances[END_REF][START_REF] Skjäl | A generalization of the classical αbb convex underestimation via diagonal and nondiagonal quadratic terms[END_REF]) are dominated by SDP based relaxation in conjunction with valid inequalities. Recall that convex underestimators of a bilinear function f are usually used to compute lower bounds of good quality over a hypercube for a nonconvex problem involving f . Therefore, for practical applications, it might be worth computing a sequence of (α) k i=1 by solving (4.7) approximately. Then a sequence of underestimators (f α i ) k i=1 can be obtained. A formula of the convex underestimator f of a bilinear function f over a hypercube H n is then

f (x) = max α 1 ,...,α k f α i (x).
By definition, K is closed, convex, and bounded. Let us consider a SDP based convex underestimator defined by q f sdp (x) := min

X∈F (x) 1 2 Q, X (4.10)
where

F (x) := {X ∈ S n : (x, X) ∈ K, X ii ≤ x i (i = 1, . . . , n)}.
For the sake of clarity, X ii ≤ x i (i = 1, . . . , n) are explicitly written out despite the redundancy w.r.t. set K. The SDP based concave overestimator f sdp can be defined symmetrically. An important observation is that values of f sdp (x), q f sdp (x) do not change if we replace

X ii ≤ x i with X ii = x i . To see this, let us introduce set F (x) := {X ∈ S n : (x, X) ∈ K, X ii = x i (i = 1, . . . , n)} and define F(x) := τ (X) ∈ R n(n+1) 2 : X ∈ F (x) and F (x) := τ + (X) ∈ R n(n-1) 2
: X ∈ F (x) . We show below that by considering the n(n -1)/2 components of τ (X) corresponding with the strictly upper triangular part of X we get a member of F (x).

Theorem 6 For any

x ∈ H n , Proj F(x) = F (x).
Proof: It is clear that F (x) ⊆ Proj F(x). We now show that the reverse also holds. Let (X ij ) i≤j be any point of F(x). It holds that Xxx T 0. Since X ii ≤ x i , it holds that Diag(xdiag(X)) + Xxx T 0. By Schur's lemma, it means that X + Diag(xdiag(X)) ∈ F (x) and therefore (X ij ) i<j ∈ F (x). Hence Proj F(x) ⊆ F (x). Theorem 6 immediately yields the desired result.

Corollary 4.1 If the objective function is bilinear, i.e., diag(Q) = 0, then for each x ∈ H n , the following equalities hold

min X∈F (x) Q, X = min X∈F (x) Q, X and max X∈F (x) Q, X = max X∈F (x) Q, X .
Proof: Since the diagonal associated terms vanish, the support of the objective function is essentially Proj F(x), which is exactly F (x) by Theorem 6. Now let us consider set QPB n and the polytope BQP n that are explicitly studied in [START_REF] Burer | On nonconvex quadratic programming with box constraints[END_REF], where

QPB n = conv (x, y) ∈ R n × R n(n+1) 2 : y ij = x i x j , 1 ≤ i ≤ j ≤ n, x ∈ H n , BQP n = conv (x, y) ∈ R n × R n(n-1) 2 : y ij = x i x j , 1 ≤ i < j ≤ n, x ∈ {0, 1} n .
The following result due to Burer and Letchford [START_REF] Burer | On nonconvex quadratic programming with box constraints[END_REF] is a consequence of Theorem 6.

Corollary 4.2 Proj x,(y ij ) i<j QPB n = BQP n Proof: It is trivial that BQP n ⊆ Proj x,(y ij ) i<j QPB n .
We show that the reverse also holds. Let all valid inequalities of QPB n be included in set K, then we have

QPB n = {(x, τ (X)) : X ∈ F (x), x ∈ H n }. By Theorem 6, we have Proj x,(y ij ) i<j QPB n = (x, τ + (X)) : X ∈ F (x), x ∈ H n . Let (x, τ + (X)) ∈ Proj x,(y ij ) i<j QPB n , then there exist some points (x k ) with associated convex combination coefficients λ k ≥ 0, k λ k = 1 such that x = k λ k x k and X = k λ k x k (x k ) T .
In addition, X ii = x i (i = 1, . . . , n) holds. These facts lead to x k i ∈ {0, 1}, ∀λ k > 0 (i = 1, . . . , n). Therefore the extreme points of Proj x,(y ij ) i<j QPB n are binary, showing that (x, τ + (X)) ∈ BQP n .

Cover based estimators

In this section, we present a novel way to construct convex and concave estimators via linear programming. To formalize the idea, we introduce definitions below. It is clear that {S} is a cover of itself. Any partition of S is a cover. To highlight the fact that BQP n defined above is associated with a set J ⊆ {1, . . . n}, we will use BQP J . We can lift BQP J to a higher dimensional space as

C J = (x, y) ∈ H n(n+1)/2 : (x J , y J ) ∈ BQP J .
q f C (x) = 1 2 min (x,τ + (X))∈C C Q, X . (4.11)
Note that C C can be explicitly described by

C C = (x, τ + (X)) ∈ H n(n+1)/2 : X c = 2 |c| l=1 λ l c V l c , x c = 2 |c| l=1 λ l c v l c , 2 |I| k=1 λ k c = 1, λ c ≥ 0, c ∈ C.
where

V l c = v l c (v l I ) T and v l c is a vertex of H |c| .
To improve the bound obtained from (4.11), one might consider a partition P as a cover of {1, . . . , n} and add RLT inequalities to bound X I,J where I, J ∈ P, i.e.,

X I,J = 2 |I| k=1 2 |J| p=1 γ kp V kp I,J , 2 |I| k=1 γ kp = λ p J , p = 1, . . . , 2 |J| , 2 |J| p=1 γ kp = λ k I , k = 1, . . . , 2 |I| , γ kp ≥ 0, k = 1, . . . , 2 |I| , p = 1, . . . , 2 |J| , ( 4.12)
where V k,p I,J = v k I (v p J ) T , ∀k, p.

Proposition 2 For any I, J ∈ P, BQP I∪J is characterized by The concave overestimator f C (x) can also be obtained in the same way by maximizing the objective function. Also, note that one can make a trade-off between the quality of estimations and computational performance via the choice of the cover C.

S = (x, τ + (X)) ∈ R |I∪J| × R ( |I∪J|
In practice, one can surely combine this cover-based estimators with valid inequalities to strengthen the formulation. To avoid adding redundant cuts, we now study the strength of cover-based estimators by establishing connections between various valid inequalities and the choice of cover sets.

Proposition 3 For BQP n , various sets of inequalities in the same dimensional space of BQP n can be identified by proper choices of cover set C of the index set {1, . . . , n} in the sense that the cover-induced polytope C C defined before can be characterized by some inequalities. For example, McCormick inequalities [START_REF] Mccormick | Computability of global solutions to factorable nonconvex programs: Part i-convex underestimating problems[END_REF][START_REF] Al-Khayyal | Jointly Constrained Biconvex Programming[END_REF] are identified by a collection of all two subsets. Triangle inequalities [START_REF] Padberg | The boolean quadric polytope: some characteristics, facets and relatives[END_REF] joint with McCormick inequalities are identified by three subsets. When the graph induced by Q ij , i < j has no K 4 minors, triangle inequalities and McCormick inequalities are equivalent to the characterization of BQP n [START_REF] Padberg | The boolean quadric polytope: some characteristics, facets and relatives[END_REF][START_REF] Deza | Geometry of cuts and metric embeddings[END_REF].

As mentioned before, we may choose a partition of the index set as a cover. Intuitively, a partition should be chosen by considering the coefficients of the bilinear function and the size of each subset. To balance these two criteria, we may consider minimizing normalized cut defined in [START_REF] Shi | Normalized cuts and image segmentation[END_REF]. Since it is NP-hard to solve the problem, one may consider the algorithm of [START_REF] Shi | Normalized cuts and image segmentation[END_REF] to solve it approximately. Alternatively, one can consider multiple random partitions of the index set at the same time. However, according to our numerical experiments, the benefits are not significant.

Extensions

In this section, we show that the convex (rep. concave) envelope of a bilinear function f (x) = i<j Q ij x i x j is affinely equivalent to the concave (rep. convex) envelope of a polyhedral function. This connection implies some well-known results in the literature. Lemma 4.2 For each pair of (x i , x j ) ∈ H 2 , it holds that

x i x j ≤ 1 2 (x i + x j -|x i -x j |)
with equality at binary points.

Let us define the function below:

φ : H n x → φ(x) = i<j Q ij |x i -x j |. (4.13)
It is clear that function φ(x) is polyhedral. In addition, it enjoys the following property.

Lemma 4.3 The convex hull of {(x, φ(x)) : x ∈ H n } is solely determined by the set of points {(x, φ(x)) : x ∈ {0, 1} n }, i.e., conv{(x, φ(x)) : x ∈ H n } = conv{(x, φ(x)) : x ∈ {0, 1} n } Proof: By definition, conv{(x, φ(x)) : x ∈ {0, 1} n } ⊆ conv{(x, φ(x)) : x ∈ H n }. Let (x, φ(x)
) be an extreme point of conv{(x, φ(x)) : x ∈ H n }. Suppose that there exists xk ∈ (0, 1) for some 1 ≤ k ≤ n, and define S k = {j : xj = xk , 1 ≤ j ≤ n}. Let be a sufficiently small quantity. Then, define x + by

x + j = ⎧ ⎨ ⎩ xj j / ∈ S k xj + j ∈ S k , j = 1, . . . , n.
Similarly, define x -by

x - j = ⎧ ⎨ ⎩ xj j / ∈ S k xj - j ∈ S k , j = 1, . . . , n,
and φ(x + ) , φ(x -) are

φ(x + ) = φ(x) - i∈S k ,j / ∈S k : x i <x j Q ij + i∈S k ,j / ∈S k : x i >x j Q ij φ(x -) = φ(x) + i∈S k ,j / ∈S k : x i <x j Q ij - i∈S k ,j / ∈S k : x i >x j Q ij . This implies that (x, φ(x)) = 1 2 (x + , φ(x + )) + 1 2 (x -, φ(x -)), which contradicts the as- sumption that (x, φ(x)) is an extreme point. Thus ext{conv{(x, φ(x)) : x ∈ H n }} = {(x, φ(x)) : x ∈ {0, 1} n }.
Note that Lemma 4.3 holds regardless of the signs of Q ij , ∀i < j. With above lemmas, we arrive at the following result.

Theorem 7 Given a bilinear function f

(x) = i<j Q ij x i x j , its convex (resp. concave) en- velope is affinely equivalent to the concave (resp.convex) envelope of φ(x) = i<j Q ij |x i -x j |, i.e., f (x) = 1 2 i<j Q ij (x i + x j ) - 1 2 q φ(x), ∀x ∈ H n , q f (x) = 1 2 i<j Q ij (x i + x j ) - 1 2 φ(x), ∀x ∈ H n . Proof: It follows from Lemma 4.2 that f (x) = 1 2 i<j Q ij (x i + x j ) - 1 2 φ(x), ∀x ∈ {0, 1} n .
As mentioned in Section 4.2, convex and concave envelopes of bilinear functions are vertex polyhedral. Lemma 4.3 show that envelopes φ are also vertex polyhedral. By the definition of (4.3), one can get the desired result. Note that Theorem 7 indicates that we can approximate concave and convex envelops of f by approximating the convex and concave envelope of φ over H n . Conversely, all techniques on approximating convex and concave envelopes of f also apply to function φ over a hypercube. In addition, this relation implies the following. Proof:

If Q ij ≥ 0, φ(x) is convex. Then it holds that for any x ∈ H n , f (x) = 1 2 i<j Q ij (x i + x j ) - 1 2 i<j Q ij |x i -x j | = i<j Q ij min {x i , x j } ,
which completes the proof. Note that this result has been proved independently in [START_REF] Luedtke | Some results on the strength of relaxations of multilinear functions[END_REF] and [START_REF] Meyer | Convex envelopes for edge-concave functions[END_REF]. We may regard the above result as an alternative proof.

Corollary 4.4 [START_REF] Mccormick | Computability of global solutions to factorable nonconvex programs: Part i-convex underestimating problems[END_REF] When n = 2, McCormick inequalities characterize the convex and concave envelopes of a bilinear function f (x) exactly.

Proof: Let n = 2, for each x ∈ H 2 , it is easy to see that

φ(x) = ⎧ ⎨ ⎩ Q 12 min{x i + x j , 2 -x i -x j } Q 12 ≥ 0, Q 12 max{x i -x j , x j -x i }, Q 12 ≤ 0. and q φ(x) = ⎧ ⎨ ⎩ Q 12 max{x i -x j , x j -x i }, Q 12 ≥ 0 Q 12 min{x i + x j , 2 -x i -x j }, Q 12 ≤ 0.
By Theorem 7, we have

q f (x) = ⎧ ⎨ ⎩ q 12 max {0, x 1 + x 2 -1} , q 12 ≥ 0, q 12 min {x 1 , x 2 } q 12 ≤ 0, and
f (x) = ⎧ ⎨ ⎩ q 12 min {x 1 , x 2 } , q 12 ≥ 0, q 12 max {0, x 1 + x 2 -1} q 12 ≤ 0.
which is exactly the McCormick's under (over-) estimator.

Numerical experiments

In this section, we compare the strength of the proposed SDP based estimators q f sdp , f sdp and cover based estimators q f C , f C by evaluating their values, as well as their gaps on a set of problem instances. q f sdp is defined in consistent with (4.10), where constraints in F (x) are

X ii ≤ x i , i = 1, . . . , n, X ij ≥ max {x i + x j -1, 0} , 1 ≤ i < j ≤ n X ij ≤ min {x i , x j } , 1 ≤ i < j ≤ n, S(x, X) 0.
Likewise, f sdp is defined symmetrically. The convex hull gap at a pre-defined point x is quantified as

gap(x) = f (x) -q f (x),
where we respectively use f sdp (x), q f sdp (x) and f C (x), q f C (x) in place of f (x), q f (x) for SDP based estimators and cover-based estimators. The instances are randomly generated following the rules below.

1. The dimension of the cube n belongs to {10, 20, 30} and the density of the matrix Q denoted by ρ belong to {0.5, 0.75, 1}.

2. For nonnegative Q, any nonzero component equals 1. For mixed signed entries of Q, 1 3 of the nonnegative components are -1 and the others are 1.

3. We randomly generate 500 test points within a hypercube with respect to different dimensions.

4. For the cover based estimators, we use a partition P as a cover. The partition of each considered index set {1, . . . , n} has been done beforehand using Shi-Malik algorithm [START_REF] Shi | Normalized cuts and image segmentation[END_REF]. The size of partition is 4 for dimension 10, 7 for dimension 20, and 10 for dimension 30.

All formulations are formulated by YALMIP [START_REF] Lofberg | Yalmip: A toolbox for modeling and optimization in matlab[END_REF] and problem instances are solved by CPLEX 12.62 with default settings on a Mac with Inter Core i5 clocked at 2.7 GHz and with 8 GB of RAM.

Table 4.1 exhibits the associated approximated values for mixed coefficients cases. For each type of estimators, the evaluations of approximated convex and concave estimators are displayed as the averages of values of 500 test points. We report several observations below.

• Both SDP-based estimators and cover-based estimators perform significantly better than McCormick's inequalities.

• Values of both SDP-based estimators and cover-based estimators are close to values of the convex and concave envelopes. Even when the density is 1, the cover-based estimators can be close to values returned by the SDP-based estimators.

• Even though the strength of SDP-based estimators is generally stronger than that of cover-based estimators. There are cases that cover-based estimators dominate the SDP-based estimators, e.g., row (10, 0.75). Of-course, the strength of cover-based estimators generally depends on the choice of the cover set.

Recall from Corollary 4.3 that the concave envelopes of bilinear functions over hypercubes are characterized by McCormick overestimators if all bilinear coefficients are nonnegative. Hence we omit values associated with concave overestimators in Table 4.2 and summarize observations below.

• The proposed cover-based convex estimators are evidently stronger than McCormick's estimators.

• If the matrix Q is not dense, e.g., 0.75, 0.5, the cover-based convex estimators can be competitive with the SDP-based convex estimators.

• If the matrix Q is very dense, the SDP-based convex estimators exhibit distinct advantages over other estimators.

q f MC f MC gap q f SDP f SDP gap q f C f C gap q f f

A perspective

As a future research direction, we briefly discuss a novel approach to construct convex underestimator when all the coefficients of a bilinear function f are nonnegative. Recall the concept of pole-set introduced in [START_REF] Ben-Ameur | Multipolar robust optimization[END_REF].

Definition 4.3 [START_REF] Ben-Ameur | Multipolar robust optimization[END_REF] Given a convex set S, a finite set Ω is called a pole-set of S if and only if any element v ∈ S belongs to the convex hull of Ω (a convex combination of poles) denoted by conv Ω.

Obviously, extreme points of S form a pole-set of itself. The extreme points of any simplex enclosing S also form a pole-set of S. The polyhedral function φ defined in (4.13) has a desired property below.

Lemma 4. [START_REF] Adjiman | A global optimization method, αbb, for general twice-differentiable constrained nlps -i. theoretical advances[END_REF] The concave envelope of φ(x) is vertex polyhedral over arbitrary polytopes if

Q ij ≥ 0, ∀i < j.
Now assume that a pole-set Ω of a hypercube H n is available, we are then able to construct a convex underestimator of f via constructing a concave overestimator of φ over a hypercube.

Theorem 8

Let Ω be a pole-set of a hypercube H n , then if

Q ij ≥ 0, ∀i < j, for each x ∈ H n , q f (x) ≥ 1 2 i<j Q ij (x i + x j ) - 1 2 φ Ω (x), (4.14)
where

φ Ω (x) = max ω∈Ω λ ω φ(ω) : λ i ≥ 0, ω∈Ω λ ω = 1, ω∈Ω λ ω ω = x .
Proof: By definition, it follows that

φ Hn (x) ≤ φ Ω (x). Since Q ij ≥ 0, ∀i < j, Lemma 4.4 implies φ Ω (x) = max ω∈Ω λ ω φ(ω) : λ i ≥ 0, ω∈Ω λ ω = 1, ω∈Ω λ ω ω = x .
It completes the proof in conjunction with Theorem 7.

One can see that when the pole-set Ω is the set of extreme points of H n , the concave envelope of φ(x) is obtained and therefore we get the exact convex envelope of f (x). Ofcourse, the quality of the approximation relies on the choice of pole-set Ω. An investigation on the construction of pole-sets has been carried out in [START_REF] Ben-Ameur | Multipolar robust optimization[END_REF]. We illustrate here the method with an example.

Example 1 Let us approximate the convex envelope values of a bilinear function f (x)

with

Q ij = 1, 1 ≤ i < j ≤ 5 over a 5-dimensional cube.
In general, we need to use 32 extreme points to get the optimality certificate. Let us take a pole-set

Ω = (1, x) : x ∈ {0, 1} 4 (0, x) : x ∈ ext{x ∈ R 4 : x -0.5 1 ≤ 2} ,
where ext represent the set of extreme points of a set. It is clear that the pole-set has 24 points, where 16 of them come from extreme points of a 4-dimensional cube and 8 of them are members of a norm-1 ball centered at 0.5 with radius 2.

For point x = (0.5, 0.5, 0.5, 0.5, 0.5), the exact value of convex envelope is q f (x) = 2, the value of McCormick relaxation gives q f Mc (x) = 0 and our pole-set based convex estimator returns q f Ω (x) = 1.0625.

Conclusion

In this chapter, we have discussed two different ways to approximate the convex and concave envelopes of a general bilinear function over hypercubes. The first approach aims at finding the best estimators of the bilinear function in a convex quadratic form at a pre-specified point, which turns out to be a semidefinite program parameterized at the point. We then established a connection between quadratic convexification based techniques and the standard semidefinite relaxation for bilinear optimization problems. We also demonstrated that two seemingly different formulations are in fact equivalent by establishing a geometric connection between their feasible regions. The second approach focused on approximating the boolean quadratic polytope by adding fully-described BQPs of lower dimensions.

In addition, some pole-set based estimators have been proposed as a perspective. Further studies are needed to find a reasonable way to define the pole-sets.

Chapter 5

Multipolar Robust Optimization

Summary

We consider linear programs involving uncertain parameters and propose a new tractable robust counterpart which contains and generalizes several other models including the existing Affinely Adjustable Robust Counterpart and the Fully Adjustable Robust Counterpart. It consists in selecting a set of poles whose convex hull contains some projection of the uncertainty set, and computing a recourse strategy for each data scenario as a convex combination of some optimized recourses (one for each pole). We show that the proposed multipolar robust counterpart is tractable and its complexity is controllable. Further, we show that under some mild assumptions, two sequences of upper and lower bounds converge to the optimal value of the fully adjustable robust counterpart. To illustrate the approach, a robust problem related to lobbying under some uncertain opinions of authorities is studied. Several numerical experiments are carried out showing the advantages of the proposed robustness framework and evaluating the benefit of adaptability.

Introduction

Uncertainty in optimization parameters arises in many applications due to the difficulty to measure data or because of their variability. To deal with uncertainty, there are mainly two approaches: stochastic optimization and robust optimization. In the first case, some probabilistic assumptions are made about the uncertain data [START_REF] Birge | Introduction to stochastic programming[END_REF][START_REF] Dantzig | Linear programming under uncertainty[END_REF][START_REF] Kall | Stochastic programming[END_REF]. One is then interested in computing a solution optimizing some moments of random variables depending on the data. Another variant, known as chance constrained programming [START_REF] Charnes | Chance-constrained programming[END_REF], consists in imposing that some constraints are satisfied only with some probability.

Robust optimization is a more recent approach dealing with uncertainty. It does not require specifications of the exact distribution of the problem's parameters. Roughly speaking, uncertain data are assumed to belong to a known compact set, called uncertainty set, and we aim at finding a solution that is immunized against all possible realizations in the uncertainty set. An early contribution related to robust optimization is the work of Soyster [START_REF] Soyster | Technical note-convex programming with set-inclusive constraints and applications to inexact linear programming[END_REF] followed by intensive investigations in the last 20 years starting with [START_REF] Ben-Tal | Robust convex optimization[END_REF][START_REF] Ei-Ghaoui | Robust solutions to least-square problems to uncertain data matrices[END_REF] in the context of convex optimization and the book [START_REF] Kouvelis | Robust discrete optimization and its applications[END_REF] dealing with discrete optimization. Almost at the same time, and in an independent way, a lot of work was initiated in [START_REF] Fingerhut | Designing least-cost nonblocking broadband networks[END_REF] and [START_REF] Duffield | A flexible model for resource management in virtual private networks[END_REF] on robust optimization in communication networks dealing with uncertain traffic matrix, see [START_REF] Ben-Ameur | Robust routing in communication networks[END_REF] for a survey.

Robust optimization and stochastic programming are related in numerous ways. For example, using some knowledge about the distribution of uncertain data, it is sometimes possible to define an uncertainty set in such a way that the robust solution is an approximated solution of a chance constrained problem (see, e.g., [START_REF] Ben-Tal | Robust optimization[END_REF][START_REF] Bertsimas | Theory and applications of robust optimization[END_REF][START_REF] Klopfenstein | A robust approach to the chance-constrained knapsack problem[END_REF] for details and references). An approach combining robust optimization and stochastic programming consists in computing solutions that are distributionally robust where the distribution of parameters is assumed to vary within some set (for example, when the mean and the covariance matrix are known) (see, e.g., [START_REF] Ghaoui | Worst-case value-at-risk and robust portfolio optimization: A conic programming approach[END_REF][START_REF] Goh | Distributionally robust optimization and its tractable approximations[END_REF]).

The definition of the uncertainty set is a critical issue since a bad choice might lead to very expensive solutions. One way to alleviate overconservatism of the robust approach is to assume that a subset of the decision variables are adjustable on the realization of the uncertain data. Let us for example consider the following linear problem

min c T x s.t. Ax ≤ b, x ∈ R n ,
involving uncertain parameters. We assume that x is partitioned as x = (u, v), where u represents the non-adjustable and v the adjustable variables. The robust counterpart of this uncertain problem under consideration reads min

u,v c T u s.t. Uu + Vv(ξ) ≤ b, [U, b] ∈ Ξ, (FARC)
where the uncertain parameters are U ∈ R m×n and b ∈ R m , while V and c are assumed to be known. We denote by ξ ≡ [U, b] ∈ Ξ the uncertain parameters belonging to the uncertainty set Ξ assumed to be compact, convex and with a non-empty interior . ξ will be considered as a vector in the rest of this chapter.

The non-adjustable variables are sometimes interpreted "here and now" variables, while the adjustable ones can be seen as "wait and see" variables. This robust counterpart above is generally called fully-adjustable robust counterpart (FARC). FARC is sometimes called the dynamic robust counterpart since v depends on ξ. FARC can be seen as a two-stage optimization problem where u are the first-stage variables and v are the second-stage variables.

If variables v are also static, then FARC simply becomes the standard static robust counterpart denoted by SRC. Cases where FARC and SRC are equivalent have been pointed out in [START_REF] Ben-Tal | Adjustable robust solutions of uncertain linear programs[END_REF] where it is shown that adaptability does not lead to any improvement in the constraint-wise uncertainty case. Still, FARC is generally much less conservative than SRC. In other words, there is generally some benefit of adaptability. Solving FARC is, unfortunately difficult in general cases as shown by many authors [START_REF] Ben-Tal | Adjustable robust solutions of uncertain linear programs[END_REF][START_REF] Chekuri | Hardness of robust network design[END_REF][START_REF] Minoux | Robust network optimization under polyhedral demand uncertainty is np-hard[END_REF]. Another concern related to FARC is the inherent difficulty of implementing the solution v(ξ) in a practical way.

To get a tractable optimization problem and also to alleviate some overconservatism of SRC, an affinely adjustable approach was proposed in [START_REF] Ben-Tal | Adjustable robust solutions of uncertain linear programs[END_REF], where the adjustable variables v are not fully adaptable (dynamic), but are assumed to depend on the uncertain data ξ in an affine way:

v(ξ) = w + Wξ, ξ ≡ [U, b] ∈ Ξ,
where w and the elements of matrix W are new decision variables (a.k.a. affine decision rules). The induced formulation is called affine-adjustable robust counterpart (AARC). An affine approach was also independently proposed for network optimization problems where the traffic matrix is supposed to be uncertain and the way how traffic is splitted through network's paths is optimized [START_REF] Kerivin | New economical virtual private networks[END_REF][START_REF] Ben-Ameur | Routing of uncertain traffic demands[END_REF][START_REF] Ouorou | A model for robust capacity planning for telecommunications networks under demand uncertainty[END_REF]. Further developments appeared in [START_REF] Babonneau | Robust capacity assignment solutions for telecommunications networks with uncertain demands[END_REF][START_REF] Ouorou | Tractable approximations to a robust capacity assignment model in telecommunications under demand uncertainty[END_REF].

Applying affine decision rules naturally leads to less expensive solutions than those obtained by the static approach. The performance gap quantified by the difference between optimum of AARC and the optimum of FARC was discussed for robust linear problems with right-hand-side uncertainty in [START_REF] Bertsimas | On the power and limitations of affine policies in two-stage adaptive optimization[END_REF][START_REF] Bertsimas | Optimality of affine policies in multistage robust optimization[END_REF]. One of the results of [START_REF] Bertsimas | On the power and limitations of affine policies in two-stage adaptive optimization[END_REF] states that AARC is equivalent to FARC when the uncertainty set is a simplex. Some tight approximation bounds relating the optimum of AARC to that of FARC in the right-hand-side uncertainty case are also given there.

Related investigation on problems with some special uncertainty sets (integer sublattices of the unit hypercube) are discussed in [START_REF] Iancu | Supermodularity and affine policies in dynamic robust optimization[END_REF], where they provide sufficient conditions such that the associated affinely adjustable decision rules lead to exact optimum of FARC.

The affine approach is related to the well-known linear or first-order decision rules used in the context of multi-stage stochastic optimization [START_REF] Garstka | On decision rules in stochastic programming[END_REF]. Linear decisions rules were also used in [START_REF] Kuhn | Primal and dual linear decision rules in stochastic and robust optimization[END_REF] in the context of stochastic programming not only to get upper bounds (as done above) but also to get lower bounds by properly approximating the dual problem using linear decision rule.

As observed in [START_REF] Chen | Uncertain linear programs: Extended affinely adjustable robust counterparts[END_REF], even though AARC has been successfully applied to several problems, its performance might be unsatisfactory under situations where the adjustable variables exhibit high nonlinearity in terms of the uncertain parameters. This led to some extensions of the affine approach in [START_REF] Chen | A linear decision-based approximation approach to stochastic programming[END_REF][START_REF] Chen | Uncertain linear programs: Extended affinely adjustable robust counterparts[END_REF] (see also references therein) by reparametrizing the uncertainties and then applying the affinely adjustable approach. Roughly speaking, a new set of variables is introduced (for example the positive and the negative parts of the original uncertainties), and the adjustable variables are assumed to affinely depend on the new set of parameters. A similar idea is also proposed in [START_REF] Bertsimas | Optimality of affine policies in multistage robust optimization[END_REF] in the context of one-dimensional constrained multistage robust optimization.

Other extensions of affine decision rules have been proposed in literature. In [START_REF] Bertsimas | A hierarchy of near-optimal policies for multistage adaptive optimization[END_REF] polynomial recourse actions are considered where v is expressed as a polynomial in uncertainty parameters with degree no larger than a fixed constant. The complexity of the robust counterpart problem is then related to testing the positivity of a polynomial. Using some recent results in algebraic geometry stating that under mild conditions, a positive polynomial can be expressed as a sum of squares (not a priori bounded), the robust counterpart is approximated by considering sums of squares of degree no larger than a fixed constant. As a sum of squares can be represented by a semidefinite programming [START_REF] Lasserre | Global optimization with polynomials and the problem of moments[END_REF], the proposed robust counterpart can be efficiently handled [START_REF] Bertsimas | A hierarchy of near-optimal policies for multistage adaptive optimization[END_REF].

Another robust approach dealing with uncertainty, termed as multi-static approach in [START_REF] Ben-Ameur | Robust routing in communication networks[END_REF], was proposed and studied in [START_REF] Ben-Ameur | Between fully dynamic routing and robust stable routing[END_REF][START_REF] Ben-Ameur | Robust routing and optimal partitioning of a traffic demand polytope[END_REF][START_REF] Zotkiewicz | More adaptive robust stable routing[END_REF]. It consists in partitioning the uncertainty set Ξ into a finite number of subsets Ξ 1 , . . . , Ξ p and using a recourse action v i for each subset Ξ i . In other words, if ξ ∈ Ξ i , then we take v = v i . The recourse actions v i are of course subject to optimization. A quite close idea is proposed in [START_REF] Bertsimas | Finite adaptability in multistage linear optimization[END_REF], where it was called finite adaptability. The performance of finite adaptability in a fairly general class of multi-stage stochastic and adaptive optimization problems was investigated in [START_REF] Bertsimas | A geometric characterization of the power of finite adaptability in multistage stochastic and adaptive optimization[END_REF].

One can also combine finite adaptability and the affinely-adjustable approach by partitioning the uncertainty set into some subsets and considering some optimized specific affine decision rules for each subset. This was also considered in [START_REF] Ben-Ameur | Between fully dynamic routing and robust stable routing[END_REF][START_REF] Ben-Ameur | Robust routing in communication networks[END_REF] in the context of network design problems. This type of adaptability might also be called piecewise-affine adaptability. Piecewise-affine rules were also considered in several other papers such as [START_REF] Bemporad | Min-max control of constrained uncertain discrete-time linear systems[END_REF][START_REF] Goh | Distributionally robust optimization and its tractable approximations[END_REF].

We mention that dynamic programming approaches have been successfully designed to solve certain robust optimization problems. For example, authors in [START_REF] Klopfenstein | A robust approach to the chance-constrained knapsack problem[END_REF] proposed a dynamic programming approach to solve the chance-constrained robust kapsack problem iteratively. Very recently, Agra et al. [START_REF] Agra | A dynamic programming approach for a class of robust optimization problems[END_REF] decomposed a general robust optimization model whose dynamic programming value function is convex and separable, to a master problem and a separation problem. To solve the separation problem, they proposed a couple of dynamic programming approaches, which allows a FPTAS.

While a great number of proposals in robust optimization have appeared, there are still challenges. First, to the best of our knowledges, none are general enough to encompass static robustness, affinely adjustable robustness and fully adjustable robustness. Second, as observed in [START_REF] Bertsimas | A hierarchy of near-optimal policies for multistage adaptive optimization[END_REF], there is no systematic way to influence the trade-off between the performance of the resulting policies and the computational complexity required to obtain them. Third, the uncertain parameters of an optimization problem can be sometimes difficult to observe. In several applications, only a subset of such parameters or some aggregates of them can be observed.

The objective of our work is to provide a framework addressing those challenges at the same time. Our contributions are four-fold:

1. A novel approach. We propose a hierarchical and convergent framework of adjustable robust optimization -multipolar robust approach, which generalizes notions of static robustness, affinely adjustable robustness, fully adjustable robustness and fill the gaps in-between. As a byproduct, a new way to look at the affine adaptability is proposed. The result of [START_REF] Bertsimas | On the power and limitations of affine policies in two-stage adaptive optimization[END_REF] stating that affine rules are optimal when the uncertainty set is a simplex is also obtained as a consequence of the multipolar approach.

2. A comprehensive analysis. We show that the multipolar robust counterpart is tractable by either a cut generation procedure or a compact formulation. Further, we prove that the multipolar approach can generate a sequence of upper bounds and a sequence of lower bounds at the same time and both sequences converge to the robust value of FARC under some mild assumptions.

3. A general constructive algorithm of pole-sets. The multipolar approach is based on some tools related to the uncertainty set, that we term as pole-sets. For their construction, we start with a simplex and then compute the best homothetic transformation of this simplex to allow it to enclose a given convex set. An efficient algorithm is proposed to compute such homothetic set. As a byproduct, we provide a very simple proof of the geometric results of [START_REF] Nevskii | Properties of axial diameters of a simplex[END_REF] related to hypercubes. The pole-sets obtained after this homothetic transformation are then improved using a tightening procedure.

4. An application. To numerically illustrate the multipolar approach, a lobbying problem is considered where a lobby aims to minimize the budget needed to convince a set of voters taking into account a reasonable opinion dynamics model under some uncertainty. The benefit of adaptability is clearly shown for this problem.

Outline. In Section 5.3, we present the concept and ingredients of multipolar robust optimization and show that static robustness, affinely adjustable robustness, fully adjustable robustness are special cases of multipolar robust framework. In Section 5.4, we discuss the tractability, the monotonicity and the convergence of the proposed approach. A simple illustrative example is described in Section 5.5. In Section 5.6, we propose algorithms for pole-set generation. Section 5.7 is dedicated to a numerical example on a lobbying problem under several uncertainty scenarios. Finally, concluding remarks follow in Section 5.8.

Notation. Throughout this chapter, we use Ξ to represent a compact convex uncertainty set and ξ to denote a member of Ξ. We use I to denote the identity matrix. Vectors and matrices are marked in bold, and their scalar components are presented in italic. Given any matrix C, C T denotes its transpose. We also use [C, D] to denote the matrix where C and D are concatenated by columns assuming they have the same number of rows. Similarly, (C, D) denotes the matrix obtained by row concatenation of two matrices C and D when they have the same number of columns. Observe that v = (v 1 , . . . , v n ) is then a vector and [v 1 , . . . , v n ] = v T . We use δ ij to represent the Kronecker's delta function, where δ ij = 1 if i = j, 0 otherwise. For a set S ∈ R n , we use ext(S) to represent the set of its extreme points, conv S to represent its convex hull and dim(S) to denote its dimension. If S is finite, we use |S| to represent its cardinality. We also use the standard notation for usual norms: • ∞ for the infinity norm, • 1 for the Manhattan norm and • 2 for the Euclidean norm.

The multipolar robust optimization concept

In this section, we introduce the main ingredients of multipolar robustness and then setup the multipolar robust counterpart as a novel approximation of FARC.

Shadow matrix. Like other robust approaches, multipolar approach is also based on an uncertainty set Ξ. In addition, we consider a matrix associated with certain operations on the uncertain information, which can be data aggregation, filtering, and selection. Note that these operations can either be natural or artificial. Natural operations are induced by the difficulty of measurements or shortage of data. For example, in communication networks, traffic flows are usually observed in an aggregated manner (the consequence of aggregating uncertain demands from multiple origin-destination pairs). Nevertheless, adjustable recourse actions should be implemented based the observed partial information. On the other hand, artificial operations can be certain techniques to control the complexity of the multipolar robust counterpart, as explained in the concluding remarks of this section. We call the associated matrix of an operation shadow matrix since the operation either reduces the size of the multipolar robust counterpart or is a direct consequence of observations. We use P ∈ R n 0 ×dim(Ξ) to denote a shadow matrix, where n 0 is the dimension of the shadow (i.e., the resulting partial information) and dim(Ξ) is the dimension of the uncertainty set Ξ. The resulting partial information is defined by Ξ P := PΞ ≡ {Pξ, ξ ∈ Ξ}.

(5.1)

When P is identity matrix, we have a complete measure of uncertainty.

We will assume that P is full row rank matrix. Consequently, Ξ P is also compact, convex and has a non-empty interior.

Pole-set. A key component of the multipolar approach is a finite set of poles, which are given vectors in the range space of the shadow matrix. We denote by Ω such a poleset. We say that Ω is a pole-set of Ξ P iff for any ξ ∈ Ξ, Pξ belongs to the convex hull of Ω (a convex combination of poles) denoted by conv Ω. Given a set Ξ P , a collection of pole-sets of Ξ P is defined as

F Ξ P := {Ω : Ξ P ⊆ conv Ω} .
(5.2)

Obviously, extreme points of Ξ P form a pole-set, i.e., ext(Ξ P ) ∈ F Ξ P . Multipolar robust counterpart. We now setup the multipolar robust counterpart w.r.t. an uncertainty set Ξ, a shadow matrix P, a pole-set Ω ∈ F Ξ P . For each ξ ∈ Ξ, we consider a weight λ ξ ω for each pole ω in Ω. Then, for each scenario ξ ∈ Ξ, the following system has a solution

ω∈Ω λ ξ ω ω = Pξ, ω∈Ω λ ξ ω = 1, λ ξ ω ≥ 0, ω ∈ Ω.
(5.3) Let Λ ξ be the set of weight vectors λ ξ satisfying the above system for a given ξ ∈ Ξ.

In the considered paradigm, each pole is associated with a recourse action, and the recourse action in the presence of ξ ∈ Ξ is approximated by a convex combination of the recourse actions associated with the poles. Specifically, let vector v ω be the recourse action associated with pole ω in the above system. We require the adjustable variables v(ξ) to be restricted to

v(ξ) = ω∈Ω λ ξ ω v ω , (5.4)
where λ ξ ∈ Λ ξ . We can readily present the multipolar robust counterpart defined by

Π Ξ (P, Ω) = min u,v c T u (MRC) s.t. Uu + V ω∈Ω λ ξ ω v ω ≤ b, ξ ∈ Ξ, λ ξ ∈ Λ ξ . (5.5)
Following the spirits of robust optimization, the multipolar robust counterpart (MRC) seeks a pair of non-adjustable solution u and a set of recourse actions related to poles v ω , ω ∈ Ω such that the objective function is minimized while hedging against the uncertainty set Ξ. In brief, given Ξ, the multipolar robust approach can bee seen as a set function of a pole-set Ω and a shadow matrix P. We denote the function by Π Ξ : R n 0 ×dim(Ξ) × F Ξ P (P, Ω) → Π Ξ (P, Ω) ∈ R and call Π Ξ (P, Ω) multipolar robust value w.r.t. (P, Ω). Also, we call (u, v) multipolar solution.

To conclude this section, we add few remarks on the concept of shadow matrix and pole-set to clarify the motivation behind these ingredients.

• Note that by (5.5), the solution is protected against the considered uncertainty Ξ.

Neither a shadow matrix P nor a pole-set Ω changes the uncertainty set, so P and Ω are not used to approximate the uncertainty set.

• Observe that 1 + dim(Ξ P ) ≤ |Ω|, so we can use the shadow matrix P to reduce the number of recourse actions and therefore the number of variables of (5.5). Reducing the number of poles leads to an MRC which is easier to solve as will be shown in Section 5.4.

• In several applications, after data is revealed, the adjustable variables should be quickly chosen and used. This is fortunately easy to do in the multipolar robust framework since the only thing to do is to find the coefficients λ ξ ω and use them to combine the already computed recourse vectors v ω , ω ∈ Ω.

Special cases

We show in this section that MRC generalizes SRC, AARC, and FARC by different settings of pole-sets and recourse actions associated with poles.

First, we show that SRC is a special case of MRC. Imposing v ω = v ω for any pair of ω and ω belonging to Ω leads to v(ξ) = v ω , ∀ξ ∈ Ξ, which means that the recourse action is static. Another way to get SRC is to impose that P is a null matrix having one row (relaxing in this case the full row rank constraint related to P ) and Ω contains just the zero vector.

Second, we show that FARC is a special case of MRC. Let Ω be the set of extreme points of Ξ and P = I. Then conv Ω = Ξ, that is for ξ ∈ Ξ, there exists λ ξ ≥ 0 such that ω∈Ω λ ξ ω = 1 and ω∈Ω λ ξ ω ω = ξ. By linearity of inequalities (5.5), imposing that Uu + Vv ω ≤ b for each extreme point ω ∈ Ω is necessary and sufficient to ensure the satisfaction of all inequalities (5.5) for each ξ ∈ Ξ. We get here the fully adjustable case representing the best that we can obtain for this problem since it is equivalent to assuming that v can vary with no restrictions. Note that if the number of extreme points of Ξ is limited, then the robust optimization counterpart can be efficiently solved. However, if the number of extreme points of Ξ is non-polynomial, the problem is generally difficult (as already mentioned in Section 5.2, see for example [START_REF] Ben-Tal | Adjustable robust solutions of uncertain linear programs[END_REF][START_REF] Chekuri | Hardness of robust network design[END_REF][START_REF] Minoux | Robust network optimization under polyhedral demand uncertainty is np-hard[END_REF]). Third, we show that AARC can also be generalized by MRC by proving the following theorem.

Theorem 9 Let Ω ∈ F Ξ P such that |Ω| = 1 + dim(PΞ). Then the optimal solution of the corresponding MRC problem is exactly the best solution that is affine in Pξ.

Proof: Since PΞ has non-empty interior, PΞ ⊂ conv Ω and |Ω| = 1 + dim(PΞ), the elements of Ω are affinely independent. Let d = dim(PΞ) and assume Ω = {ω (1) , . . . , ω (d+1) }. The shadow matrix P is here the identity matrix. Consider matrix D obtained by taking vectors ω (i) as columns and adding a final line containing only coefficients equal to 1.

D =       ω (1) 1 . . . ω (d+1) 1 . . . . . . . . . ω (1) d
. . . ω

(n+1) d 1 . . . 1       .
Observe that D is a non-singular square matrix of size (d + 1). Given any ξ, there are unique coefficients λ ξ ω such that Pξ = ω∈Ω λ ξ ω ω and

ω∈Ω λ ξ ω = 1.
This can be written as Pξ, 1 = Dλ ξ where λ ξ is the vector whose components are the λ ξ ω for ω ∈ Ω. This immediately implies that λ ξ = D -1 Pξ, 1 . Using E to denote the matrix whose columns are the recourse vectors v ω , equation (5.4) becomes v = ED -1 Pξ, 1 . This clearly implies that v affinely depends on Pξ.

Let us now consider any affine policy w + WPξ. As shown above, the recourse vector v provided by the multipolar approach is given by ED -1 Pξ, 1 . By taking E = [W, w]D, we get v = w + WPξ. In other words, any recourse policy that is affine in Pξ can be obtained through the multipolar approach. When P = I, we get the desired result below.

Corollary 5.1 The affinely adjustable approach is a special case of the multipolar approach. It corresponds to any set of (dim Ξ + 1) affinely independent poles, in multipolar robust optimization when P = I.

The following corollary is also immediate. Corollary 5.2 If the uncertainty set Ξ is a simplex, then the affinely adjustable robust counterpart is equivalent to the fully adjustable robust counterpart in the sense that their objective values are equal.

Proof: Taking all the vertices of the simplex uncertainty set as the set of poles in multipolar robust approach leads to the optimum of FARC. By Corollary 5.1, this poleset corresponds to affine adjustable approach, which completes the proof. Corollary 5.2 has been presented in [START_REF] Bertsimas | On the power and limitations of affine policies in two-stage adaptive optimization[END_REF] in the special case of right-hand-side uncertainty, so we may treat the result here as an alternative proof using the framework of multipolar approach.

Analysis

In this section, we first analyze the tractability of the multipolar robust counterpart MRC. Then, we show that the proposed framework can generate a monotonic sequence converging to the fully adjustable robust value of FARC. In fact, we will simultaneously generate a lower and an upper bound both converging to the optimal value of FARC under some mild assumptions.

Tractability

In this section, we show that MRC is computationally tractable. It can be solved either by cut generation or using a compact reformulation.

First, a cutting plane algorithm for solving MRC may be devised as follows. Assume that |Ω| is finite and has a reasonable size. Given a solution (u, v), we have to check if there exists a pair of ξ ∈ Ξ and λ ξ ∈ Λ ξ violating the constraints of MRC. This can be done by checking the sign of the optimum of each i th problem max λ,ξ

U i u + V i ω∈Ω λ ξ ω v ω -b i (5.6a) s.t.
λ ξ ω ≥ 0, ω ∈ Ω, (5.6d) ξ ∈ Ξ, (5.6e)
where U i and V i are the i th rows of U and V. If it is positive, then constraint

U i u + V i ω∈Ω λξ ω v ω ≤ b i , (5.7)
needs to be added to the restricted problem, where ( λξ , ξ) solves (5.6). Problem (5.6) can generally be solved easily when Ξ is polyhedral or ellipsoidal. In these cases, by equivalence of separation and optimization [START_REF] Grötschel | Geometric algorithms and combinatorial optimization[END_REF], the multipolar robust optimization counterpart problem can also be solved in polynomial time if the number of poles |Ω| is polynomially bounded. Second, we may solve MRC by duality. It is sometimes possible for several kinds of convex uncertainty sets to write a strong dual of (5.6) leading to an extended reformulation of MRC. This holds for example if Ξ is a polytope defined by a limited number of constraints, i.e., Ξ :

= {ξ ≡ [U, b] : Cξ ≤ d}, where C = [C 1 , . . . , C m ], C i ∈ R n d ×(n+1) , d ∈ R n d
and ξ is expressed as a column vector of size (n + 1) × m. ξ contains m blocks of size n + 1 vectors: the i th block contains U T i followed by b i .

By strong duality, the constraints of the multipolar robust counterpart MRC w.r.t. Ξ can be replaced with a polynomial number of inequalities. For each i, the inequalities

U i u + V i ω∈Ω λ ξ ω v ω ≤ b i , ξ ∈ Ξ, λ ξ ∈ Λ ξ are replaced with d T η i + V i v ω -ω T σ i ≤ 0, ω ∈ Ω, C T j η i -P T j σ i = δ ij u, -1 , j = 1, . . . , m, η i ∈ R n d + , σ i ∈ R n 0 , (5.8)
where the shadow matrix P = [P 1 , . . . , P j , . . . , P m ], P j ∈ R n 0 ×(n+1) , j = 1, . . . , m, δ ij is Kronecker's delta function.

When Ξ is ellipsoidal, i.e., Ξ := {ξ : Fξ 2 ≤ 1}, the multipolar robust counterpart can be represented by a second order cone program. Then for each i, the i th constraint of MRC is replaced with

η i 2 + V i v ω -ω T σ i ≤ 0, F T η i -L i = 0, η i ∈ R nq , σ i ∈ R n 0 .
(5.9)

where L i = L i1 , . . . , L im , L ij = δ ij u, -1 + P T j σ i , j = 1, . . . , m, n q is number of rows of matrix F. For sake of completeness, a proof of (5.9) is provided in Appendix.

Monotonicity

We show in this section that the function Π Ξ (P, •) is monotonic w.r.t. a partial order defined on F Ξ P when the shadow matrix P is fixed.

Given an uncertainty set Ξ, we now define a partial order over the collection of its pole-sets F Ξ P denoted by F Ξ P . We set members of F Ξ P ordered by the inclusion of their convex hulls, i.e., for any Ω , Ω ∈ F Ξ P , Ω F Ξ P Ω ⇐⇒ conv Ω ⊆ conv Ω.

(5.10)

The next theorem emphasizes the fact that the function Π Ξ (P, •) is monotonic regarding the partial order F Ξ P for each fixed P ∈ R n 0 ×dim(Ξ) . In other words, the multipolar value gets smaller when Ω is smaller w.r.t. F Ξ P .

Theorem 10 Given P ∈ R n 0 ×dim(Ξ) , for any Ω , Ω ∈ F Ξ P , if Ω F Ξ P Ω, then we have Π Ξ (P, Ω ) ≤ Π Ξ (P, Ω).

Proof: If (u, (v ω ) ω∈Ω ) is an optimal solution of MRC, then a feasible solution, when the set of poles is defined by Ω , is given as follows. Each ω ∈ Ω writes as a convex combination of the poles of Ω:

ω = ω∈Ω λ ω ω ω. Let v ω = ω∈Ω λ ω ω v ω .
The solution given by (u, (v ω ) ω ∈Ω ) is clearly feasible for MRC w.r.t. the set of poles defined by Ω , which completes the proof. Theorem 10 not only implies that the smaller the Ω w.r.t. F Ξ P , the lower the multipolar robust value, but also implies that for a given P, Π Ξ (P, Ω) is minimum if Ω = ext(Ξ P ).

Given two pole-sets Ω, Ω ∈ F Ξ P , Theorem 10 also indicates that: first, for a fixed shadow matrix P, if |Ω| > |Ω |, then Π Ξ (P, Ω) is not necessarily less than Π Ξ (P, Ω ); second, the function Π Ξ (P, •) is not strictly monotonically increasing. For example, let S, S ∈ F Ξ P , S F Ξ P S and their convex hulls are simplices. By Theorem 9, Π Ξ (P, S ) = Π Ξ (P, S) while by Theorem 10, Π Ξ (P, S ) ≤ Π Ξ (P, S), which illustrates the second point. Now take any pole-set Ω whose cardinality is strictly greater than 1 + dim(Ξ P ), such that S F Ξ P Ω F Ξ P S, then Π Ξ (P, S ) = Π Ξ (P, Ω) = Π Ξ (P, S), which illustrates the first point.

Observe also that when P = I, any pole-set whose convex hull contains Ξ is contained in a simplex. This immediately implies that the optimal value of AARC represents the worst that can be obtained by the multipolar approach.

Convergence

The aim of this section is to show that under some mild assumptions, using the multipolar framework, one can simultaneously compute a sequence of upper bounds and a sequence of lower bounds converging to Π Ξ (I, ext (Ξ)), the optimal robust value of FARC. Throughout this section the shadow matrix is the identity matrix. (5.13)

We obviously have Ξ z Ω ⊆ Ξ. Let (u * , (v * zω ) ω∈Ω ) be the optimal solution of the MRC problem related to Ξ z Ω . Due to the definition of Ξ z Ω , MRC and FARC are equivalent. Moreover, from Ξ z Ω ⊆ Ξ, we get that

c T u * = Π Ξ z Ω (I, ext(Ξ z Ω )) ≤ Π Ξ (I, ext(Ξ)).

An illustrative example

To illustrate the multipolar concept, we present a simple example, which had been previously studied in [START_REF] Chen | Uncertain linear programs: Extended affinely adjustable robust counterparts[END_REF] and is as follows:

min u s.t. ∀ξ ∈ R n , ξ 1 ≤ 1, ∃v, v i ≥ ξ i , v i ≥ -ξ i , i = 1, . . . , n u ≥ n i=1 v i . (5.19)
Observe that u is here the unique first-stage (non adjustable) variable. On the other hand v i , for each i = 1, . . . , n, are second-stage (adjustable) variables. The uncertainty set is given by Ξ ≡ {ξ ∈ R n , ξ 1 ≤ 1}. As noticed in [START_REF] Chen | Uncertain linear programs: Extended affinely adjustable robust counterparts[END_REF], an optimal fully adjustable solution is given by u = 1 and v i = ξ i 1 , whereas the optimal affinely adjustable solution requires that u = n. In other words, the affine approach does not lead to any improvement compared to the static approach.

Following the paradigm of multipolar approach in Section 5.3, let us take Pξ = (ξ 1 , . . . , ξ n 0 ), where n 0 ∈ N, n 0 ≤ n. In other words, the shadow matrix P limits the dimension of Ξ to n 0 by leaving ξ i as they are for i ≤ n 0 and disregarding the other components for i > n 0 . Let Ω ⊆ R n 0 be the set of poles containing for i = 1, . . . , n 0 , vectors φ i = (0, . . . , 0, 1, 0, . . . , 0) and φ i = -φ i , whose components are 0 except the i th component. Hence Ω contains 2n 0 poles and Ξ P = conv Ω. Given any ξ ⊆ Ξ, let λ φ i and λ φ i be the convex combination coefficients such that

Pξ = n 0 i=1 λ φ i φ i + λ φ i φ i . The equation can be transformed to Pξ = n 0 i=1 φ i λ φ i -λ φ i ;
thus these coefficients should satisfy the equations λ φ i -λ φ i = ξ i for 1 ≤ i ≤ n 0 . Let v φ i (resp. v φ i ) be the recourse vector associated with pole φ i (resp. φ i ). These vectors belong to R n . In the considered example, inequalities (5.7) are equivalent to the following set of inequalities:

n 0 i=1 λ φ i v φ i + λ φ i v φ i ≥ (|ξ 1 |, |ξ 2 |, . . . , |ξ n |) , (5.20a) u ≥ n 0 i=1 λ φ i v φ i + λ φ i v φ i 1 .
(5.20b) Let us take v φ i = v φ i = (0, . . . , 0, 1, 0, . . . , 0, 1, . . . , 1), where the first n 0 components are 0 except the i th component, which is equal to 1, while the last (n -n 0 ) components are equal to 1.

Observe that the last (n -n 0 ) components of the vector

n 0 i=1 λ φ i v φ i + λ φ i v φ i are equal to 1. Moreover, for 1 ≤ i ≤ n 0 , we have λ φ i + λ φ i ≥ |λ φ i -λ φ i | = |ξ i |.
This clearly implies that inequalities (5.20a) are satisfied. In addition, inequality (5.20b) leads to u ≥ n 0 i=1 λ φ i + λ φ i (1 + n -n 0 ), where (1 + n -n 0 ) is the L 1 norm of each recourse vector v φ i . Consequently, u ≥ 1+n-n 0 . Since we are minimizing u, we get u = 1+n-n 0 . The cost decreases when n 0 increases. When n 0 is equal to 1, we get a static solution,

The monotonicity of multipolar robust values w.r.t. pole-sets might suggest using minimum volume circumscribed polytopes. Considering the Nikodym distance (related to volumes) instead of the Hausdorff distance, the same kind of results can be obtained [START_REF] Bronstein | Approximation of convex sets by polytopes[END_REF]. One might be interested in a minimum volume simplex containing a convex set Ξ P . We know for example that if Ξ P is the hypercube H n 0 , then a minimum volume circumscribed simplex has a volume equal to n 0 n 0 n 0 ! [START_REF] Lassak | Parallelotopes of maximum volume in a simplex[END_REF]. If Ξ P is the unit ball, then a minimum volume simplex containing the ball is a regular simplex whose volume is

n n 0 /2 0 (n 0 +1) (n 0 +1)/2
n 0 ! [START_REF] Bronstein | Approximation of convex sets by polytopes[END_REF] and whose dihedral angle is arccos( 1 n 0 ) [START_REF] Parks | An elementary calculation of the dihedral angle of the regular n-simplex[END_REF]. It is also known that a minimum volume simplex enclosing Ξ P satisfies the centroid property: the centroid of each facet of this simplex should be in Ξ P [START_REF] Klee | Facet-centroids and volume minimization[END_REF]. A polynomial-time algorithm to find such a minimum volume simplex enclosing a set of points in R 3 is given in [START_REF] Zhou | Algorithms for a minimum volume enclosing simplex in three dimensions[END_REF]. However, it is generally unknown how to solve the problem in higher dimensions [START_REF] Gritzmann | On the complexity of some basic problems in computational convexity[END_REF].

As observed by [START_REF] Bronstein | Approximation of convex sets by polytopes[END_REF], most constructive algorithms were generally proposed for low dimensional cases (2 or 3). For more general cases, constructive algorithms of circumscribed polytopes such as the algorithm of [START_REF] Kamenov | Algorithm of approaching polyhedrons[END_REF] are generally based on the addition of inequalities without controlling the number of vertices of the circumscribed polytope. This can hardly accommodate the need of multipolar framework since we want to control the complexity of MRC by limiting the number of poles.

Note also that we are required to construct the pole-set of Ξ P in a reasonable time. Algorithms checking whether each extreme point of Ξ P belongs to the convex hull of Ω fail to work, since the number of extreme points of a polytope can be exponential or even infinite.

The rest of this section is organized as follows. First, we describe a general algorithm to construct a simplex enclosing Ξ P . The resulting simplex is guaranteed to be smallest in the sense that it cannot be shrinked. Then, a project-and-cut based tightening procedure is proposed to construct pole-sets that are closer to Ξ P .

Generation of a circumscribed simplex

In this section, we describe a general algorithm for the construction of a circumscribed simplex of Ξ P . Specifically, we first randomly generate a set of (n 0 + 1) affinely independent points, whose convex hull forms a simplex S. Then, we compute the best homothetic transformation of S such that the resulting simplex contains Ξ P .

We denote by the ω (i) , i = 1, . . . , (n 0 + 1) the (n 0 + 1) affinely independent points. Then the n 0 -simplex set can be expressed as {x : Dλ = (x, 1), λ ≥ 0}, where

D =       ω (1) 1 . . . ω (n 0 +1) 1 . . . . . . . . . ω (1) n 0 . . . ω (n 0 +1) n 0 1 . . . 1       .
Since the (n 0 + 1) points are affinely independent, matrix D is invertible; therefore, λ j , j = 1, . . . , (n 0 + 1), can be expressed as an affine function of x; the coefficients of the Observe that having l i(n 0 +1) (and σ * as a consequence), we can get the translate t through the linear system n 0 k=1 t k l ik = σl i(n 0 +1) -l i(n 0 +1) , i = 1, . . . , (n 0 + 1). Multiplying by D, we get that (t 1 , ..., t n 0 , 0) = σ(0, ..., 0, 1) -D(l 1(n 0 +1) , ..., l (n 0 +1)(n 0 +1) ) which leads

to t = n 0 +1 i=1 z i ω (i) .
According to (5.24), the smallest σ * is -

n 0 +1 i=1 z i . We should however check if σ * ≥ 0.
This holds since by considering any x ∈ Ξ P , one can write that

n 0 +1 i=1 z i ≤ n 0 +1 i=1 n 0 j=1 l ij x j = n 0 j=1 x j n 0 +1 i=1 l ij =0,
where the last equality is based on (5.23).

Since the matrix D -1 σ,t does not exist when σ * = 0, we have to study this special case. It is clear that σ * = 0 if and only if Ξ P is a single point. Observe that in this case, we necessarily have

n 0 +1 i=1 z i = 0 since z i = n 0 j=1 l ij x j , where x is the single point of Ξ P . Then formula σ * = - n 0 +1 i=1 z i is still valid and t * = x = n 0 +1 i=1 z i ω (i) also occurs.
Note that values of z i , i = 1, . . . , (n 0 +1) defined in Theorem 12 can easily be computed for any Ξ P since we only have to minimize a linear function over a convex set.

As a special case, pole-sets of a hypercube are of great use in multipolar robust approach. First, hypercubes are one of the most common uncertainty sets in many applications. Second, general box sets of the form {x : x ∈ [l, u] n 0 ⊆ R n 0 } are simply affine transformations of a hypercube, so the pole-sets of a hypercube also apply to boxes with some simple transformations.

Corollary 5.4 If Ξ P is a hypercube, then σ * and t * are given by:

t * = n+1 i=1 n 0 j=1 min{0, l ij }ω (i) and σ * = 1 2 n 0 +1 i=1 n 0 j=1 |l ij |.
Proof: If Ξ P is a hypercube, by Theorem 12, we have

z i = n 0 j=1 min{0, l ij }, i = 1, . . . , n 0 .
By (5.23), we have -

n 0 +1 i=1 n 0 j=1 min{0, l ij } = 1 2 n 0 +1 i=1 n 0 j=1
|l ij |, which completes the proof.

According to Corollary 5.4, we have a closed formula for the homothetic translation for a n 0 -simplex S containing the n 0 -hypercube, i.e,

x σ,t = 1 2 n 0 +1 i=1 n 0 j=1 |l ij |x + n 0 +1 i=1 ω (i) n 0 j=1 min {0, l ij } , ∀x ∈ S. (5.25)
Note that the value σ * presented in Corollary 5.4 has been given in [START_REF] Nevskii | Properties of axial diameters of a simplex[END_REF] but the proof here is much simpler. To sum up the foregoing, we present a general algorithm for the generation of a circumscribed simplex as follows.

ω 0 ω 1 ω 2 ω 3 ω 1 ω 2 ω 3
1. Generate (n 0 + 1) affinely independent points (ω (i) ) i=n 0 +1 i=1 .

2. Compute σ * and t * by Theorem 12 and output the σ * ω (i) + t * , i = 1, . . . , (n 0 + 1).

A tightening procedure

In this section, we propose a general procedure to construct pole-sets of good quality by tightening a given pole-set.

The procedure is the following: among the vertices of Ω select the farthest one in L 2 sense from Ξ P and compute the projection of this vertex on Ξ P . Then we consider the hyperplane separating this vertex from Ξ P (containing the projection) and compute the extreme points of the intersection of this hyperplane with conv Ω. These extreme points are added to Ω while the vertex that has been projected is removed from Ω. Figure 5.1 illustrates a tightening procedure of a 3-D simplex covering H 3 . The procedure is repeated until the cardinality of Ω reaches some fixed upper bound. Details are given below:

1. Assume Ω = ω (k) , k ∈ I . For each k ∈ I, compute the distance between ω (k) and Ξ P . Let z k be the projection of ω (k) on Ξ P . z k can be usually expressed in a closed form. For example, in the ball case, we have

z k = ω (k) ω (k) 2
, while in the hypercube case we get

z k i = ω (k) i if ω (k) i ∈ [0, 1], z k i = 1 if ω (k) i ≥ 1 and z k i = 0 if not.
The distance between ω (k) and Ξ P is then given by ω (k)z k 2 . Let ω (k 0 ) be the vertex of Ω maximizing the distance from Ξ P :

ω (k 0 ) = arg max ω (k) ∈Ω ω (k) -z k 2 .
2. Let α = ω (k 0 )z k 0 and let B ω (k 0 ) , α 2 be the ball of radius α 2 centered at ω (k 0 ) . Since B ω (k 0 ) , α 2 ∩ Ξ P = {z k 0 } and B ω (k 0 ) , α 2 and Ξ P are convex, there is a hyperplane separating them. This hyperplane, denoted by h(ω (k 0 ) ), is here uniquely defined since it contains z k 0 and is orthogonal to α. It is then given by h ω (k 0 ) = x : (x -z k 0 ) T α = 0 . We use h -ω (k 0 ) = x : (x -z k 0 ) T α < 0 and h + ω (k 0 ) = x : (x -z k 0 ) T α ≥ 0 to respectively denote the inner and outer half spaces.

3. Now partition the vertices ω (k) k∈I into two disjoint sets: Ω -and Ω + , where Ω -= ω (k) : ω (k) ∈ h -ω (k 0 ) and Ω + = ω (k) : ω (k) ∈ h + ω (k 0 ) . Then consider the set of vertices Ω obtained as intersections between the hyperplane h(ω (k 0 ) ) and the set of lines (ω (i) , ω (j) ) where ω (i) ∈ Ω + and ω (j) ∈ Ω -: Ω = ω (i) ∈Ω + ,ω (j) ∈Ω -h(ω (k 0 ) )∩ (ω (i) , ω (j) ). The number of such intersections is of course less than |Ω -|×|Ω + |. Also note that we need to remove redundant points from Ω if they are convex combinations of other points of Ω . Finally update Ω by deleting Ω + and adding Ω : Ω = Ω -∪ Ω . If cardinality of Ω is still under a prescribed upper bound, the procedure is repeated.

To conclude this section, we might add that it is sometimes more efficient to start with a pole-set having more than (n 0 + 1) poles. Assume, for example, that Ξ P is the unit ball {x ∈ R n 0 : x 2 ≤ 1}. Then one can consider a 2n-pole-set where poles are the 2n 0 extreme points of x ∈ R n 0 : x 1 ≤ √ n 0 . Of course, 2n-pole-sets can also be easily generated for many other convex sets.

A numerical example: the lobbying problem

Let us consider a lobbying problem where a set of voters (for example, legislators or members of regulatory agencies) have to take some decisions. The opinion of each voter depends on the opinion of some authorities. Authority's opinions are generally uncertain.

A lobby would like to ensure that an important decision will be unanimously approved by all voters. The lobby will spend some effort (energy, money, etc.) to convince each voter, while the total lobbying budget is minimized. Assume that there are m voters and n authorities. The opinion of voter i is given by n j=1 Q ij ξ j where Q ij is an estimated number belonging to [-1, 1] and ξ j represents the uncertain opinion of authority j. If Q ij is close to 1, then j has a big impact on i, while Q ij = 0 means that i does not care about j, while a negative value of Q ij can be interpreted as a negative effect (i.e., when j recommends something, i is inclined to have an opposite opinion). We assume here that the lobby would be satisfied if n j=1 Q ij ξ j ≤ 0 for each voter. Since this might not occur for some voters, some effort modeled here by v i (ξ) can be made by the lobby to convince them. The total effort is quantified by m i=1 r i v i (ξ) where r i is a unit effort price corresponding with voter i. The problem can be formulated as follows, min u (5.26) where Q ∈ [-1, 1] m×n , Ξ is the convex uncertainty set and u is the budget that has to be secured by the lobby. The lobby problem is related to opinion dynamics in social networks (see [START_REF] Acemoglu | Opinion dynamics and learning in social networks[END_REF] and the references therein). Notice that interactions between voters are also possible since the set of authorities might include the set of voters as a subset.

s.t. m i=1 r i v i (ξ) ≤ u, ξ ∈ Ξ, Qξ ≤ v(ξ), ξ ∈ Ξ, v(ξ) ≥ 0, ξ ∈ Ξ,
To illustrate the multipolar robust approach, we consider here two different uncertainty sets: the hypercube H n and a unit volume ball B n . The numbers r i are assumed to be equal to 1. Specializing (5.8) to H n , we get the following formulation for MRC associated with a shadow matrix P and a feasible pole-set Ω ∈ F P Hn , Π Hn (P, Ω) = min

u,vω u s.t. m i=1 v i ω + n j=1 σ j ω j + n j=1 β j ≤ u, ω ∈ Ω, n j=1 η T ij ω j + n j=1 α ij -v i ω ≤ 0, ω ∈ Ω, i = 1, . . . , m, n j=1 µ ij -V i ω - n j=1 δ ij ω j ≤ 0, ω ∈ Ω, i = 1, . . . , m, β + P T σ ≥ 0, α ij + P T η i ≥ Q ij , i = 1, . . . , m, j = 1, . . . , n, µ i -P T τ i ≥ 0, i = 1, . . . , m, α ij , β j , µ ij ≥ 0, i = 1, . . . , m, j = 1, . . . , n, σ, η i , τ i ∈ R n 0 , i = 1, . . . , m.
(5.27)

As said above, we also consider the case where the uncertainty set is a unit volume ball B n , whose center is ξ = ( 1 2 , . . . , 1 2 , . . . , 1 2 ) and radius ρ = ( Γ(n/2+1) π n/2 ) 1 n . According to (5.9), the dual of the multipolar robust counterpart w.r.t. B n , writes

Π Bn (P, Ω) = min u,vω u s.t. ρ P T σ 2 + 1 T v ω -ω T σ ≤ u, ω ∈ Ω, ρ Q T i + P T η i 2 -v i ω -ω T η i + Q i ξ ≤ 0, ω ∈ Ω, i = 1, . . . , m, ρ P T µ i 2 -v i ω -ω T µ i ≤ 0, ω ∈ Ω, i = 1, . . . , m, σ, η i , µ i ∈ R n 0 , i = 1, . . . , m, (5.28)
where Ω ∈ F P Bn .

As stated in Section 5.3.1, the fully adjustable robust value w.r.t. hypercube H n can be achieved by simply taking P = I and Ω = ext(H n ). The problem looks more complex in the ball case since the number of extreme points is infinite. Assume again that P is identity. Given ξ ∈ B n , the optimal solution of v i is max{0, Q i ξ} for each i = 1, . . . , m. Denote by P the power set of the index set {1,. . . ,m }. We partition the ball B n into a family of disjoint subsets by a set valued mapping S : P → 2 Bn , i.e., for each J ∈ P,

S(J) := ξ ∈ B n : Q i ξ ≥ 0, i ∈ J, Q j ξ ≤ 0, j ∈ J ,
where J = {1, . . . , m} \ J. Therefore, the fully adjustable robust program writes

Π * Bn = max J∈P:S(J) =∅ max ξ∈S(J) i∈J Q i ξ.
(5.29)

Notice that (5.29) takes an exponential number (in the number of constraints m) of seconder order cone programs to obtain the fully adjustable robust value Π * Bn . We show that (5.29) is equivalent to a much simpler problem. We show that it is also a lower bound of Π * Bn . Let J max ∈ P be a subset for which the maximum is achieved:

J max = arg max J∈P ρ i∈J Q i 2 + i∈J Q i ξ . Take ξ ∈ B n such that ξ = ξ + ρ i∈Jmax Q i i∈Jmax Q i 2 . Let K ∈ P such that Q i ξ ≥ 0, i ∈ K, Q i ξ ≤ 0, i ∈ {1, . . . , m} \ K.
We have that

ρ i∈Jmax Q i 2 + i∈Jmax Q i ξ = i∈Jmax Q i ξ ≤ i∈K Q i ξ ≤ max ξ∈S(K) i∈K Q i ξ ≤ Π *

Bn

where the first equality follows from the choice of ξ. The second inequality is due to J max ∩ K ⊆ K, while the third inequality is from the fact that ξ belongs to S(K) by the definition of K. The last inequality is a direct consequence of (5.29).

Although problem (5.30) is easier than problem (5.29), it is still computationally costly when the number of constraints m is large. Notice that (5.30) can also be seen as a integer quadratic program that can be approximated by semidefinite programming and solved using standard quadratic programming tools. We will not elaborate more on this since this falls out of the scope of this thesis.

Numerical experiments

The problem instances are randomly generated following the rules below.

1. Let m ∈ {10, 20, 30, 40, 50}, n ∈ {5, 9, 10, 12, 15, 20, 30}.

2. Generate the components of Q uniformly over [-1, 1].

3. We build four different sizes of pole-sets for each considered hypercube by the circumscribed simplex generation algorithm and the tightening procedure described in Section 5.6. As a result, for a hypercube H n , (Ω i ) is a monotonic sequence w.r.t. the set inclusion of their convex hulls, i.e., Ω i Hn Ω j for all i > j. Table 5.1 displays the cardinality of different pole-sets of H n . The number of vertices of H n is also provided in the last column. 4. As an illustration of multipolar robust approach for smooth convex uncertainty sets, we generate pole-sets of ball B n as well. In Table 5.2, Ω 0 denotes a first pole-set whose convex hull is a simplex, while Ω 1 is the 2n-pole-set defined at the end of Section 5.6. Starting from Ω 1 and applying the tightening procedure, we get the pole-sets Ω i , i = 2, 3, 4 as outputs. The cardinality of pole-sets are shown in Table 5.2. The pole-sets had been readily generated before the solution procedure. Compact formulations (5.27) and (5.28) are modeled by YALMIP [START_REF] Lofberg | Yalmip: A toolbox for modeling and optimization in matlab[END_REF] and all the problem instances are solved by the Linux version of CPLEX 12.5 with default settings on a Dell E6400 laptop with Intel Core(TM)2 Duo CPU clocked at 2.53 GHz and with 4 GB of RAM. We evaluate our multipolar approach from different measuring: the impact of pole-sets, the impact of the shadow matrix P, and the benefit of adaptability.

The influence of pole-sets

Recall that multipolar robust approach closes to some extend the gap between affine robust value and the fully adjustable value. To test the impact of pole-sets, we fix P = I. For relatively lower dimensional cases, we report the multipolar robust values w.r.t. different pole-sets, and compute the percentage of the closed gap induced by the multipolar robust approach Π Ξ (Ω 0 )-Π Ξ (Ω i ) Π Ξ (Ω 0 )-Π * Ξ × 100. For higher dimensional cases, where FARC can hardly be solved in a reasonable time, we report the multipolar robust values w.r.t. different pole-sets.

Results related to hypercubes are presented in Table 5.3. The closed gap percentages are given within parentheses. Overall, these results appear encouraging as indicated by the closed gaps. Observe that when the uncertainty set is fixed, the multipolar robust values in general get lower as the pole-set Ω gets smaller. Also, we report the computing time for higher dimensional instances associated with hypercube uncertainty set in Table 5.4. While the computational time compared with the affine robust approach scales in magnitude, the complexity of multipolar robust approach is controlled by the choice of pole-sets. In particular, in higher dimensional cases, where fully adjustable robust values are difficult to obtain, the robust multipolar approach brings lower cost (compared with affine approach) in a reasonable time.

A sequence of lower bounds can also be generated as stated in Corollary 5.3 of Sec- i=0 onto the surface of hypercubes. The obtained lower bounds are denoted by Π Hn (Γ i ). Note that conv Ω ⊆ Ω does not necessarily lead to conv Γ ⊆ Γ or conv Γ ⊇ Γ. Thus it may happen that Π Hn (Γ i ) ≥ Π Hn (Γ i+1 ). The results are summarized in Table 5.5, where the best lower bound for each problem instance is marked in bold. balls and hypercubes are the same. Second, the closed gaps by multipolar robust approach on robust problems with ball uncertainty sets might be less significant than that with hypercube uncertainty sets. As might be expected, larger poles-sets are required for balls compared to hypercubes. Third, despite the limitations, multipolar approach closes around 30% of the optimality gap. In particular, it appears compelling when the number of constraints are large, while the dimension of the uncertainty set is small. The lower bounds obtained in the ball case are reported in Table 5.7. Interestingly, the observed sequences of lower bounds associated with ball B n are monotonically increasing and their best bounds in general are close to the fully adjustable robust value.

The impact of the shadow matrix

To investigate the impact of the shadow matrix on the robust value of the robust problem, we conduct some experiments on problem instances w.r.t. hypercube uncertainty sets.

The shadow matrices considered here are simply projection matrices on lower subspaces. The results are displayed in Table 5.8, where the uncertainty set is a hypercube and several projections are considered (on H 5 , H 7 , H 10 , and H 12 ). The pole-set considered is the set of extreme points of the projected set. As might be expected, the robust value decreases as more information is employed in MRC. The benefit of adaptability

To illustrate the concept of benefit of adaptability in the framework of multipolar robust approach, we compute the multipolar robust values of problem (5.26) with different proportions of adjustable variables v. We allow the first θm components of v to be adaptable to the realization of ξ, while keeping the remaining m -θm variables independent of the realization of ξ, where θ ∈ [0, 1]. Note that when θ = 0, we get the static case SRC. The results are summarized in Table 5.9 and we emphasize here two observations:

1. For each problem instance, as the adaptability ratio θ increases, the robust value decreases significantly, which is reasonable both in theory and practice.

2. As the adaptability ratio θ increases, the influence of pole-sets on the robust value increases. For example, when the adaptability ratio θ = 0.25, the multipolar robust values of all problem instances remain the same with different pole-sets except instance [START_REF] Belotti | Disjunctive cuts for nonconvex minlp[END_REF][START_REF] Belotti | Disjunctive cuts for nonconvex minlp[END_REF]. When the adaptability increases, the robust values of more instances improve as the better pole-sets are used, which can be clearly seen when θ = 0.75 and θ = 1.

Conclusion

In this chapter, we have presented a novel approach to handle uncertainty in optimization problems called the multipolar robust approach, which is based on a set of poles Table 5.9: The benefit of adaptability (hypercube) (m,n), Ω 1 θ = 0.25 θ = 0.5 θ = 0.75 θ = 1 of the proposed approach in terms of the robust values, the complexity, and the benefit of adaptability. In general, the results appear encouraging.

It would be interesting to investigate further the performance of the multipolar robust approach on other problems. A systematic study of good approximations of convex bodies by enclosing polytopes with a limited number of extreme points should help to alleviate overconservatism and get closer to the optimal fully adaptable robust value. One can also put more focus on the approximation of convex bodies from inside using, for example, maximum volume inscribed polytopes to get better lower bounds of the fully adjustable robust value.

While the approach was proposed in the context of a two-stage optimization problem, it can be adapted to multistage optimization. Multipolar decision rules can also be considered in stochastic programming. The multipolar approach might also be combined with finite adaptability or multi-static robustness by partitioning the uncertainty set into several subsets and considering some multipolar decision rules for each subset.

Appendix: the derivation of (5.9)

We derive the compact formulation (5.9) of Section 5.4.1 w.r.t. an ellipsoidal uncertainty set defined by Ξ := {ξ : Fξ 2 ≤ 1}.

For each i th constraint, MRC requires the optimum of the following problem nonpositive. max ξ,s,λ≥0

U i u -b i + ω∈Ω λ ξ ω V i v ω s.t. s 2 ≤ 1, Fξ = s, η i ∈ R nq ω∈Ω λ ξ ω ω = Pξ, σ i ∈ R n 0 ω∈Ω λ ξ ω = 1, τ i ∈ R
where ξ ≡ [U, b] and η i , σ i , τ i are dual multipliers corresponding to each group of constraints. Consider the corresponding Lagrangian

L (λ, ξ, s, η i , τ i , σ i ) = U i u -b i + ω∈Ω λ ξ ω V i v ω + η T i (s -Fξ) + σ T i ω∈Ω λ ξ ω ω -Pξ + τ i ω∈Ω λ ξ ω -1 .
The dual function is then max λ,ξ, s 2 ≤1

L (λ, ξ, η i , τ i , σ i ). Setting the derivative w.r.t. λ, ξ leads to the dual constrains

V i v ω + τ i -ω T σ i ≤ 0, F T η i -L i = 0, (5.31)
where P = [P 1 , . . . , P k , . . . , P m ] , L i = L i1 , . . . , L im , L ij = δ ij u, -1 + P T j σ i , j = 1, . . . , m. The dual objective is min η i ,τ i ,σ i η i 2 -τ i .

(5.32)

By duality, the optimum of the above dual problem is equal to the optimum of the problem. Thus restricting the non-positivity of the primal optimum can be equivalently represented as

η i 2 + V i v ω -ω T σ i ≤ 0, F T η i -L i = 0, η i ∈ R nq , σ i ∈ R n 0 .
(5.33)

Chapter 6

Conclusions

To conclude this thesis, we briefly summarize our main contributions, and outline some directions of future research.

Evaluation

This thesis is motivated by a mapping problem where a set of virtual machines need to be assigned to servers with multiple 0-1 quadratic constraints. Two main challenges exist. First, to the best of our knowledge, efficient exact solution procedures or approximation schemes with performance guarantees of bounds do not appear in the literature. Second, the demands of computing resources and networking resources are highly uncertain. Few investigations have been conducted to incorporate uncertainty properly in the context of robust optimization. It is well-known that both problems are difficult. Thus this thesis aims at proposing and applying novel relaxation techniques to handle MIQCPs and robust programs.

In Chapter 3, an exact model dealing with static mapping procedure is provided. Relaxations involving certain standard techniques such as RLT inequalities are proposed. To further strengthen the formulation, two types of strong cuts are proposed by exploiting the problem structure, which numerically prove quite effective. Then, we have tried a number of decompositions based on Lagrangian relaxations with a hope that the Lagrangian dual bound is much stronger (compared with the bound provided by the continuous relaxation of a strong formulation) while the solution procedure is efficient. A request based decomposition scheme is then discussed in Section 3.4. We also provided a geometric result on the strength of the decomposition scheme, which implies a novel hierarchy of relaxations by grouping the set of requests. Our numerical results show that the proposed hierarchy can handle relatively larger problem instances much more efficiently than the standard branch-and-bound routines of CPLEX12.6.3. Finally, we proposed a reformulation exploiting symmetries of virtual requests, which scales the proposed algorithms greatly.

Chapter 4 discusses a couple of relaxation approaches for bilinear optimization problems over a hypercube from the perspective of convex and concave envelopes. A connection between quadratic convexification techniques and Shor's relaxation strengthened by RLT inequalities was established. We showed that two seemingly different formulations are in fact same, which generalizes a related result from Burer and Letchford [START_REF] Burer | On nonconvex quadratic programming with box constraints[END_REF]. Then we proposed a novel relaxation approach by considering a predefined tool to approximate the convex and concave envelope of f over a hypercube. This approach leads to a linear reformulation and we also established its connections with various inequalities and the RLT method proposed in [START_REF] Adams | A hierarchy of relaxations leading to the convex hull representation for general discrete optimization problems[END_REF].

To address the issues of uncertainty, we proposed a new tractable robust paradigm called multipolar robust optimization in Chapter 5. It provides a new way to model the recourse action for a fairly general linear program with uncertain parameters. It generalizes notions of static robustness, affinely adjustable robustness, fully adjustable robustness and fills the gaps in-between. The result of [START_REF] Bertsimas | On the power and limitations of affine policies in two-stage adaptive optimization[END_REF] stating that affine rules are optimal when the uncertainty set is simplex in the special case of right-hand-side uncertainty is generalized as a consequence of the multipolar approach. Moreover, we show that the multipolar robust counterpart is tractable by either a cut generation procedure or a compact formulation. Further, we prove that the multipolar approach can generate a sequence of upper bounds and a sequence of lower bounds at the same time and both sequences converge to the robust value of fully adjustable robust counterpart under some mild assumptions. The multipolar approach is based on some tools related to the uncertainty set, that we call as pole-sets. For their construction, we start with a simplex and then compute the best homothetic transformation of this simplex to allow it to enclose a given convex set. An efficient algorithm is proposed to compute such homothetic set. Several numerical experiments are carried out showing the advantages of the proposed robustness framework.

Perspectives

As future research directions, it would be interesting to investigate the following topics, which have been left due to the lack of time.

Regarding the mapping problem discussed in Chapter 3, we have the following perspectives. First, we may integrate the proposed decomposition hierarchy into the branchand-bound procedure of CPLEX with a hope of accelerating the solution procedure further. We might also consider combining the symmetry-induced model and decomposition hierarchy together to deal with much larger problem instances. Second, it might be worthwhile developing a spatial branch-and-bound scheme to solve the extended mapping model discussed in Section 3.6. We can also consider certain convex and concave envelops over some polytopes associated with assignment constraints and location constraints. Third, a future project could be applying the multipolar robust optimization paradigm to the mapping problem with uncertain parameters associated with virtual requests.

In Chapter 4, we have proposed a perspective on approximating the convex and concave envelopes using the idea of pole-set, which covers an n-dimensional hypercube with much less extreme points. A valuable direction might be to combine the idea of poleset and certain valid inequalities for QPB polytope to derive strong estimators for general quadratic forms rather than bilinear forms. And, it might be interesting to integrate ideas of cover-based estimators and pole-set based estimators in order to find the most expressive linear formulation to approximate the estimators of bilinear functions.

To extend the work of multipolar robust optimization, we may investigate the approximation of convex bodies from inside using, for example, maximum volume inscribed polytopes, to get better lower bounds of the fully adjustable robust value. We can also extend our work to model recourse actions that are binary or general integer variables.

All these topics will be pursued in future research.

  Nous nous intéressons à des problèmes linéaires incertains min c T x s.t. Ax ≤ b, x ∈ R n . Nous supposons que x est divisée en x = (u, v), où u représente les variables non ajustables et v les variables ajustables. La contrepartie robuste de ce problème incertain à l'étude s'écrit min u,v c T u s.t. Uu + Vv(ξ) ≤ b, [U, b] ∈ Ξ, (FARC) où les paramètres incertains sont U ∈ R m×n et b ∈ R m alors que V et c sont supposées être connues. On note par ξ ≡ [U, b] ∈ Ξ les paramètres incertains appartenant à l'ensemble d'incertitude fixé Ξ supposé être compact, convexe et avec un intérieur non vide. ξ sera considéré comme un vecteur.

4. 1

 1 An evaluation of convex and concave estimators (mixed entries) . . . . . . 4.2 An evaluation of convex estimators (nonnegative coefficients) . . . . . . . 5.1 The pole-sets of hypercubes . . . . . . . . . . . . . . . . . . . . . . . . . . 5.2 The pole-sets of the unit volume ball . . . . . . . . . . . . . . . . . . . . . 5.3 The multipolar robust values with different pole-sets (hypercube uncertainty sets) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.4 Computing time (in seconds) . . . . . . . . . . . . . . . . . . . . . . . . . 5.5 Lower bounds related to hypercubes . . . . . . . . . . . . . . . . . . . . . 5.6 The multipolar robust values with different pole-sets (ball) . . . . . . . . . 5.7 The lower bounds in the ball case . . . . . . . . . . . . . . . . . . . . . . . 5.8 Impact of the shadow matrix . . . . . . . . . . . . . . . . . . . . . . . . . 5.9 The benefit of adaptability (hypercube) . . . . . . . . . . . . . . . . . . .

Figure 1 . 1 :

 11 Figure 1.1: An illustration of the mapping procedure

  is a polynomial (with degree 1) in variables x over the real field. We represent a polynomial c by c(x) = α c α x α , where α ∈ N n d := α ∈ N n : n i=1 α i ≤ d and x α stands for the monomial x α i 1 . . . x αn n . We identify c by its coefficients (c α ) α∈N n d . As observed by Lasserre [101], it holds that p * := min x∈K c(x) = min µ∈M(K) K c dµ.

  Indeed, let us denote by z d := (x α = x α 1 1 . . . x αn n : n i=1 α i ≤ d) the vector that contains all monomials of degree at most d. Then for any u j ∈ R[x] d (a polynomial with degree at most d), we have that u j (x) = u T j z d and thus l j=1

P

  V1 : {(P RLT1 )} ∩ {(Cut), (Cut')} , (P V1 ) P V2 : {(P RLT2 )} ∩ {(Cut), (Cut')} .(P V2 )

PLemma 3 . 2

 32 v1 : {(P V1 )} ∩ {(3.8), (3.9), (3.10), (3.11), (3.12)} . (P v1 ) The projection of the continuous relaxation of F(P v1 ) in variables (x, y, θ, φ) is exactly the continuous relaxation of set F(P V1 ).

  y) : (3.8), (3.9), (3.10), (3.11), (3.12)} .

n = 1 ,

 1 . . . , N, γ e ≤ L e (σ, ζ) -D e φ n e (σ e -σ n e ) -φ n e r∈R (ζ r e -ζ rn e ), e ∈ E, n = 1, . . . , N, r∈R γ r + k∈S

Algorithm 1 :

 1 In-out AlgorithmLet (v * in , γ in ), (v * out , γ out ) be the initial in-point and out-point choose an α ∈ (0

Definition 4 . 1 A

 41 cover of a set S is a collection C of subsets of S whose union is the whole set S, i.e., S = c∈C c.

Definition 4 . 2

 42 The polytope induced by a cover C of {1, . . . , n} is defined asC C = c∈C C c .Lemma 4.1 Any cover C-induced polytope is an outer approximation of the BQP n , i.e., BQP n ⊆ C C . Moreover, the equality holds if {1, . . . , n} ∈ C. Proof: Let (x, y) be an arbitrary point in BQP n . It is easy to verify that (x, y) satisfies constraints of C c , for all c ∈ C. Thus it holds that (x, y) ∈ C C and therefore BQP n ⊆ C C . If {1, . . . , n} ∈ C, it holds that C C ⊆ BQP n , which leads to BQP n = C C . By Lemma 4.1, the convex underestimator of a bilinear function f can be constructed as

Corollary 4 . 3

 43 If Q ij ≥ 0, the concave envelope of a bilinear function f is characterized by McCormick over-estimator.

Definition 5 . 1 ω -ξ 2 .

 512 Let Ω ∈ F Ξ be a pole-set of a non-empty set Ξ, the distance function between them is defined asd (Ω, Ξ) = max ω∈Ω min ξ∈ΞThe distance function is well-defined since Ω and Ξ are closed and bounded. It characterizes the furtherest distance between pole-set Ω and the uncertainty set Ξ. This distance is nothing other than the well-known Hausdorff distance.Let Ω ∈ F Ξ such that d (Ω, Ξ) = . For each ω ∈ Ω, we have d (ω, Ξ) ≤ . Let z ω be the projection of ω on Ξ, i.e.,z ω = argmin ξ∈Ξ d (ω, ξ) , e ω = ω -z ω ,(5.11)where e ω 2 ≤ . For each ξ ∈ Ξ, consider convex combination coefficients β subscripts to avoid confusion:ξ ≡ [U ξ , b ξ ], E ≡ [U E , b E ] and z ω ≡ [U zω , b zω ]. We define the convex set Ξ z Ω = conv {z ω : ω ∈ Ω} .

Figure 5 . 1 :

 51 Figure 5.1: An illustration: (Left) Pole ω 0 is replaced with poles (ω i ) 3 i=1 . (Right) The updated convex hull of the new set of poles.

Lemma 5 . 2 Q

 52 Program (5.29) is equivalent to i ξ is an upper bound of Π * Bn since it is obtained by relaxing the constraints ξ ∈ S(J).

  KP)

	m ri x ri k ≤ M k θ k ,	∀k ∈ S,	(KP')
	r∈R i∈V r
		f rij x ri k x rj
	r∈R k,p∈S:	(i,j)∈L r
	k =p,e∈P kp

p ≤ B e φ e , ∀e ∈ E, (QC)

  Glover. Par conséquent, nous allons nous concentrer sur la linéarisation McCormick. Ensuite, nous utilisons la RLTs et nous ajoutons des inégalités valides pour renforcer le modèle McCormick. Notre premier résultat montre que les inégalités RLT peuvent éliminer de nombreuses inégalités McCormick redondantes. En outre, nous avons proposé deux types d'inégalités valides en exploitant la structure du problème. Dans l'ensemble, plusieurs formulations sont proposées et ils sont codés en C ++ et résolus par CPLEX 12.6.3. Leurs performances de calcul sont évaluées sur certains cas de problème. Les résultats numériques montrent que la combinaison des types d'inégalités valides peut accélérer sensiblement la résolution du problème. Nous renvoyons le lecteur au Tableau 3.2 pour plus de détails.

•

  Une formulation exacte et plusieurs reformulations pour le problème d'affectation de machines virtuelles dans le cloud. Nous utilisons plusieurs inégalités valides pour renforcer la formulation exacte, accélérant ainsi l'algorithme de résolution de manière significative. En outre, une décomposition lagrangiennes efficace et une formulation induite par symétrie sont proposées. Un nouveau paradigme sur les problèmes linéaires généraux avec des paramètres incertains. Nous proposons une hiérarchie convergente de problèmes d'optimisation robuste -approche robuste multipolaire, qui généralise les notions de robustesse statique, robustesse affinement ajustable, et la robustesse entièrement ajustable. Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2 Thesis motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2.1 Mapping of cloud virtual machines . . . . . . . . . . . . . . . . . . Virtual machine types . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.1 Numerical results for P and P MC . . . . . . . . . . . . . . . . . . . . . . . 3.2 Numerical results for different formulations . . . . . . . . . . . . . . . . . 3.3 Evaluation of Lagrangian decomposition . . . . . . . . . . . . . . . . . . .

	• Deux approches sont proposées pour la construction d'enveloppes convexes et con-caves pour l'optimisation bilinéaire sur un hypercube. Nous établissons plusieurs connexions théoriques entre différentes techniques et nous discutons d'autres exten-sions possibles. 1 Introduction • Contents 1.1

1.2.2 Problem background . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2.3 Related work and challenges . . . . . . . . . . . . . . . . . . . . . 1.3 Content of the PhD . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.3.1 The scope . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.3.2 Outline and contributions . . . . . . . . . . . . . . . . . . . . . . . 2 Literature Review 2.1 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.3 General algorithmic frameworks . . . . . . . . . . . . . . . . . . . . . . . . 2.4 Polyhedral relaxations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.4.1 Disjunctive cuts . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.5 Nonlinear relaxations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.5.1 Semidefinite relaxations . . . . . . . . . . . . . . . . . . . . . . . . 2.5.2 Quadratic convexifications . . . . . . . . . . . . . . . . . . . . . . . 2.6 Lagrangian relaxation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.7 Recent paradigms . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.7.1 Copositive programming . . . . . . . . . . . . . . . . . . . . . . . . 2.7.2 Moment-SoS hierarchies . . . . . . . . . . . . . . . . . . . . . . . . 3 Optimal Mapping of Cloud Virtual Machines 3.1 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3 A general model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3.1 McCormick inequalities and Glover's linearization . . . . . . . . . 3.3.2 Polyhedral relaxations . . . . . . . . . . . . . . . . . . . . . . . . . 3.3.3 Numerical experiments . . . . . . . . . . . . . . . . . . . . . . . . . 3.4 A Lagrangian decomposition . . . . . . . . . . . . . . . . . . . . . . . . . 3.4.1 A request based decomposition . . . . . . . . . . . . . . . . . . . . 3.4.2 The solution procedure . . . . . . . . . . . . . . . . . . . . . . . . i List of Tables 1.1 3.4 Evaluation of the Lagrangian decomposition hierarchy . . . . . . . . . . . 3.5 Problem sizes of (P) and (SP) . . . . . . . . . . . . . . . . . . . . . . . . . 3.6 Performances of (SP V1 ) and (P V1 ) . . . . . . . . . . . . . . . . . . . . . . 3.7 Upper bounds of (P) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.8 Continuous relaxation statistics of (P MC ) and (P GL ) . . . . . . . . . . . . 3.9 Branch and Cut statistics of (P MC ) and (P GL ) . . . . . . . . . . . . . . .

Table 1 .

 1

		1: Virtual machine types
		CPU (core) Memory (GB)
	Small	1	1
	Medium	2	2
	Large	4	4
	XLarge	8	16

  GL ) outperforms McCormick's relaxation for some small instances. However for most instances, even though model (P GL ) introduces much less continuous variables than McCormick's relaxation, its overall computational performance are generally no better than those of model (P MC ). Thus we believe that the main difficulty is not the number of variables and constraints. Besides, we can apply RLT inequalities to strengthen (P MC ) and remove a number of redundant constraints (see Lemma (3.1)), whereas it appears complicated to apply RLTs to model (P GL ). Hence, we will focus on techniques strengthening the McCormick linearization. One may refer to Appendix of this chapter for some numerical results on the comparison between two models.

	.2b)
	Consequently, constraints (QC) become,

(k,p):k =p,e∈P kp r∈R,i∈Vr w ri kp ≤ D e φ e , e ∈ E. (3.3) The resultant model is min {f (θ, φ, x) : (AC), (LC), (KP), (KP'), (BC), (3.2a), (3.2b), (3.3)} . (P GL ) Our numerical experiments indicate that model (P

Table 3 .

 3 1: Numerical results for P and P MC

				Branch-and-Cut statistics
	(|S|, |E|) #Req. #Bin.	P		P MC
				#CPU(s) #Nodes	#CPU(s) #Nodes
	(8, 20)	1	68	3.33	88	0.03	0
	(8, 20)	2	108	10.00	1704	0.45	229
	(12, 30)	1	102	371	53305	1.18	16
	(12, 30)	2	162	-	-	134	8897
	(12, 30)	4	282	-	-	4180	10766

-Time limit: 36000 Seconds (10hs)

Table 3 . 2 :

 32 Numerical results for different formulations

	relaxation statistics
	Continuous

Table 3 .

 3 [START_REF] Adams | A hierarchy of relaxations leading to the convex hull representation for general discrete optimization problems[END_REF] show that the request based decomposition provides strong lower bounds efficiently in the following senses.1. For small problem instances that CPLEX can solve to optimality in a reasonable time, Algorithm 1 can provide near optimality bound quickly.2. For instances that CPLEX takes hours of CPU time (e.g. instances with 5 and 6 requests), Algorithm 1 can produce strong bounds in hundreds of seconds. Additionally, Algorithm 1 is generally about ten times faster than the standard CPLEX branch-and-cut algorithm in reaching the same optimality gap (in the last column of Table3.3.

			Table 3.3: Evaluation of Lagrangian decomposition
	(|S|, |E|) #Req. #Bin.	P V1		Decomposition 3.13
				#CPU(s)	Root gap(%)	#CPU(s) Gap(%)
	(12, 30)	2	162	15.40	9.70	21.38	0
	(12, 30)	3	222	65.08	9.26	60.92	0.55
	(12, 30)	4	282	125.78	8.52	93.23	1.85
	(12, 30)	5	342	3050.21	20.15	139.78	12.72
	(12, 30)	6	402	32492.00 (9hs)	24.09	278.56	16.21

Table 3 .

 3

			4: Evaluation of the Lagrangian decomposition hierarchy
	(|S|, |E|) #Req. #Bin.		Level 1	Level 2
				#CPU(s) Closed gap(%)	#CPU(s) Closed gap(%)
	(12, 30)	8	522	535.12	38.91	907.42	45.01
	(12, 30)	10	642	612.71	32.99	1534.88	35.63
	(12, 30)	20	1284	1924.31	18.29	3247.95	20.12
	-Time limit: 36000 Seconds (10hs)

  .7. The optimality gap induced by an upper bound is quantified by Gap = Upper bound -Exact value

	Exact value

Table 3 .

 3

			6: Performances of (SP V1 ) and (P V1 )
	(|S|, |E|) #VMs.	(SP V1 )		(P V1 )
			Cpu(Sec.) #Nodes #Cpu(sec.) #Nodes
	(12, 30)	10	1.22	0	13.68	0
	(25, 90)	10	424	121	1216	30068
	(25, 90)	20	1080	150	3680	73808
	(40, 178)	10	522	0	1024	25197
	(40, 178)	20	824	121	13428(3.7hs) 310124
	(40, 178)	30	4780	875	-	-
	(40, 178)	40	4410	1950	-	-
	Time limit: 36000 Seconds (10hs).

Table 3 .

 3

	(12, 30)	10	2	2.72	0	9.22
	(12, 30)	10	4	3.67	0	4.21
	(25, 90)	10	2	739	121	11.69
	(25, 90)	10	4	1280	187	8.90
	(25, 90)	20	4	2182	478	23.70
	(25, 90)	20	7	4821	157	15.84
	(25, 90)	25	4	1328	952	24.33
	(25, 90)	25	10	3120	1253	13.89

7: Upper bounds of (P) (|S|, |E|) #VMs. N #Cpu(sec.) #Nodes Gap(%)

Table 3 .

 3 8: Continuous relaxation statistics of (P MC ) and (P GL )

	(|S|, |E|) #Req.	P MC		P GL
			Cpu(Sec.) Gap(100%) #Cpu(sec.) Gap(100%)
	(8,20)	1	0.00	63.9	0.00	63.9
	(8,20)	2	0.02	60.8	0.01	60.8
	(12,30)	1	0.02	79.6	0.01	79.6
	(12,30)	2	0.07	74.7	0.06	74.7
	(12,30)	4	0.38	65.88	0.64	65.88

Table 3 .

 3

			P MC		P GL
			Cpu(Sec.) #Nodes #Cpu(sec.) #Nodes
	(8,20)	1	0.29	0	0.53	0
	(8,20)	2	0.45	229	0.38	251
	(12,30)	1	1.18	16	0.52	284
	(12,30)	2	134	8897	129	9016
	(12,30)	4	4180	10766	3966	14750

9: Branch and Cut statistics of (P MC ) and (P GL ) (|S|, |E|) #Req.

  Then (x I , x J ) can be written as a finite convex combination of binary points of H |I|+|J|

	by imposing γ kp = λ k I λ p J . Also notice that	k=1 2 |I|	2 |J|
	Remark 4.2 Alternative to adding constraint set (4.12) to formulation (4.11) associated
	with partition P, one can also generate a new cover set whose members are all possible
	unions of two subsets of P, thereby obtaining an equivalent formulation. One can check
	that both formulations have the same number of constraints and variables. For example

2 ) : (4.12), C I ∩ C J Proof: It holds that BQP I∪J ⊆ S since all the inequalities of S are valid inequalities to BQP I∪J . So it suffices to show that the reverse also holds. Suppose (x, τ(X)) ∈ S. j=1 γ kp = 1. Therefore (x I∪J , τ + (X)) can be written as a finite convex combination of vertices of BQP I∪J ending the proof. if P = {p 1 , p 2 , p 3 }, the newly generated cover set can be {p 1 ∪ p 2 , p 1 ∪ p 3 , p 2 ∪ p 3 }.

Table 4 .

 4 1: An evaluation of convex and concave estimators (mixed entries)

	MC-estimators	SDP-estimators	cover-estimators	Envelopes
	(n, Density)

Table 5 .

 5 1: The pole-sets of hypercubes

	Hypercube |Ω 0 |	|Ω 1 |	|Ω 2 |	|Ω 3 |	#ext.
	H 9	10	32	162	387	512
	H 10	11	36	112	322	1,024
	H 12	13	44	144	449	4,046
	H 15	16	56	192	353	32,768
	H 20	21	76	144	514	1,048,576
	H 30	31	60	116	432	1,073,741,824

Table 5 .

 5 2: The pole-sets of the unit volume ball

	Unit ball	|Ω 0 |	|Ω 1 |	|Ω 2 |	|Ω 3 |	|Ω 4 |
	B 5	6	10	118	218	308
	B 9	10	18	62	152	352
	B 10	11	20	72	132	374
	B 12	13	24	88	164	478

Table 5 .

 5 3: The multipolar robust values with different pole-sets (hypercube uncertainty sets)

	(m,n)	static	affine/ΠH n (Ω0)	ΠH n (Ω1)	ΠH n (Ω2)	ΠH n (Ω3)	Π * Hn
	(10,9)	24.84	12.42	12.18(9.45 )	10.55(73.62)	10.16(88.98)	9.88
	(10,10) 25.50	12.75	11.53(58.65)	10.96(86.06)	10.70(98.56)	10.67
	(10,12) 30.66	15.33	14.63(35.18)	13.71(81.41)	13.43(95.48)	13.34
	(20,9)	50.75	25.37	23.82(46.69)	22.08(99.10)	22.06(99.70)	22.05
	(20,10) 50.88	25.44	23.56(16.95)	20.74(42.38)	18.58(61.86)	14.35
	(20,12) 59.79	29.89	27.54(23.64)	25.40(45.17)	23.58(63.48)	19.95
	(10,15) 35.81	17.90	16.91	15.98	15.29	-
	(10,20) 50.88	25.44	24.82	24.35	23.33	-
	(10,30) 64.32	32.16	31.70	31.14	30.22	-
	(20,15) 82.49	41.25	39.36	36.20	34.87	-
	(20,20) 99.28	49.64	47.17	45.66	40.80	-
	(20,30) 157.20	78.60	77.82	76.83	73.77	-

Table 5 .

 5 4: Computing time (in seconds) (m,n) static affine/Π Hn (Ω 0 ) Π Hn (Ω 1 ) Π Hn (Ω 2 ) Π Hn (Ω 3 ) All we need to do is to generate a sequence of (Γ i ) i=3 i=0 by projecting the pole-sets (Ω i ) i=3

	(10,15)	0.00	0.01	0.20	1.49	6.31
	(10,20)	0.00	0.03	0.87	2.66	27.54
	(10,30)	0.00	0.04	1.15	4.53	39.12
	(20,15)	0.00	0.03	0.68	10.43	34.98
	(20,20)	0.00	0.07	4.65	16.41	152.79
	(20,30)	0.00	0.14	3.91	15.40	220.48
	tion 5.4.3.

Table 5 .

 5 5: Lower bounds related to hypercubes (m,n) Π Hn (Γ 0 ) Π Hn (Γ 1 ) Π Hn (Γ 2 ) Π Hn (Γ 3 ) Π *

	Hn

Table 5 .

 5 6: The multipolar robust values with different pole-sets (ball) .54(20.92) 19.51(22.88) 19.45(26.80) 19.36(32.68) 18.33 (20,10) 33.12 17.19 16.42(18.08) 15.62(36.85) 15.57(38.03) 15.36(42.96) 12.93 (20,12) 39.85 20.93 20.05(17.25) 19.96(19.02) 19.90(20.20) 19.81(21.96)

	(m,n)	static	ΠB n (Ω0)	ΠB n (Ω1)	ΠB n (Ω2)	ΠB n (Ω3)	ΠB n (Ω4)	Π * Bn
	(10,9)	17.27	9.43	9.28(11.11)	9.21(16.30)	9.11(23.70)	9.01(31.11)	8.08
	(10,10) 16.09	9.21	8.98(23.00)	8.95(26.00)	8.94(27.00)	8.89(32.00)	8.21
	(10,12) 19.64	10.93 10.76(20.99) 10.76(20.99) 10.74(23.46) 10.70(28.40)	10.12
	(20,9)	35.87	19.86 1915.83
	(30,5)	19.51	11.00	10.39	9.22	9.14	9.03	-
	(40,5)	37.57	21.63	21.11	20.71	20.63	20.55	-
	(50,5)	38.14	20.94	20.06	19.33	19.14	19.02	-

Table 5 .

 5 7: The lower bounds in the ball case

	(m,n)	Π Bn (Γ 0 )	Π Bn (Γ 1 )	Π Bn (Γ 2 )	Π Bn (Γ 3 )	Π Bn (Γ 4 )	Π * Bn
	(10,9)	6.60	6.70	6.70	6.71	7.67	8.08
	(10,10)	6.07	6.28	6.93	7.44	7.60	8.21
	(10,12)	7.26	7.51	8.01	8.08	8.44	10.12
	(20,9)	16.04	16.17	16.89	16.89	16.89	18.33
	(20,10)	12.03	12.03	12.03	12.03	12.03	12.93
	(20,12)	12.35	13.41	13.67	13.67	13.67	15.83
	(30,5)	7.48	7.48	8.11	8.11	8.34	-
	(40,5)	16.54	17.31	19.27	19.27	19.87	-
	(50,5)	15.53	15.53	17.22	17.22	17.80	-

Table 5 .

 5

			8: Impact of the shadow matrix
	(m,n)	H 5	H 7	H 10	H 12
	(10,30)	56.95	51.50	45.59	43.29
	(10,50)	109.84	105.68	100.69	96.11
	(10,70)	162.14	159.52	156.12	154.53
	(10,100)	246.98	244.13	239.34	237.43
	(20,30)	133.35	126.59	117.63	111.53
	(20,50)	233.73	224.82	213.04	203.89
	(20,70)	345.56	340.41	328.49	319.79
	(20,100)	462.51	453.26	439.71	431.39

Special cases where unbounded variables are necessary, e.g.,[START_REF] Liberti | Mathematical programming: Turing completeness and applications to software analysis[END_REF], will not be considered in model (2.1).

-Time limit: 54000 Seconds (15hs).

We will also assume that there is a positive number µ such that (u * , 1) 2 ≤ µ. This assumption generally holds. For example, if the cost vector c is positive and variables u are non-negative, then c T u * ≤ Π Ξ (I, ext(Ξ)) implies that (u * , 1) 2 is upper-bounded. The number µ does not depend on .

Assumption 3 There exists a constant number µ such that (u * , 1) 2 ≤ µ for any Ξ z Ω ⊆ Ξ and any optimal solution u * of the FARC problem related to Ξ . Lemma 5.1 Under Assumption 3, for each ξ ∈ Ξ, (U E u * -b E ) is bounded from above by µ1, where 1 is an all-ones vector.

Proof: The result follows from Cauchy-Schwartz inequality applied to each row of E ≡

Let δ be a small positive number and let

We will assume that for some small number δ, the static robust counterpart problem SRC is still solvable.

Assumption 4 There exists a static robust solution (u δ , v δ ) w.r.t. uncertainty set Ξ δ .

Theorem 11 Under Assumptions 3 and 4, for each pole-set

where c T u * and c T u δ are respectively fully adjustable robust cost w.r.t. Ξ z Ω and static cost w.r.t. Ξ δ .

Proof: Assume that the optimal solution of FARC w.r.t. uncertainty set

Let us show that (û, vω ) is a feasible solution of the MRC problem related to Ξ and Ω.

For any ξ ≡ [U, b] ∈ Ξ, by (5.12), one can write:

where (5.16) follows from the fact that u * , v * zω satisfies constraint Uu

The robust cost incurred by (û, (v

Moreover, the corresponding sequence of sets Ξ z Ω i defined in (5.13) satisfies

Proof:

and we know that

Consequently,

hold in the limit at the same time.

while the optimal fully adjustable solution is obtained when n 0 = n. Finally, taking 1 < n 0 < n, we obtain a compromise between the simplicity of the static approach and the efficiency of the fully adjustable solution. As mentioned earlier, such a compromise cannot be obtained for this example with the affinely adjustable approach.

Consider now a slightly changed example with the uncertainty set being the nonpolyhedral set defined by Ξ := {ξ ∈ R n : ξ 2 ≤ 1}. The rest of the problem remains as in (5.19); thus the new problem can be formulated as follows:

(5.21)

Observe that since the new uncertainty set Ξ contains the previous one based on L 1 norm, the optimal value of (5.21) is greater than or equal to that of (5.19). Optimal solutions based on either the static approach or the affine approach still incur a cost of n while the optimal fully adjustable solution has a cost of √ n. Let us now consider the multipolar approach, where P is still defined by P ξ = (ξ 1 , . . . , ξ n 0 ). Let us choose the following set of poles:

. One can easily show that Ξ P ⊆ conv(Ω). Moreover, by taking v √ n 0 φ i = v√ n 0 φ i = (0, . . . , 0, √ n 0 , 0, . . . , 0, 1, . . . , 1), where the first n 0 components are 0 except the i th component, which is equal to √ n 0 , while the last n -n 0 components are equal to 1, we get a solution of the multipolar robust counterpart with u = √ n 0 + n -n 0 . Similarly to the previous case, when n 0 is equal to 1, we get a static solution, while the optimal fully adjustable solution is obtained when n 0 = n.

The construction of pole-sets

We know from Section 5.4 that the multipolar robust value converges to a fully adjustable robust value when the distance between Ω and Ξ P gets close to 0, and P = I. We also proved the monotonicity of multipolar robust value w.r.t. the inclusion of conv Ω. Therefore, the objective of this section is to find a pole-set Ω ∈ F P as close to Ξ P as possible, while minimizing the number of poles. This is clearly related to the theory of approximation of convex sets by polytopes.

A considerable amount of work has been done in this area. A recent survey of relevant results is given in [START_REF] Bronstein | Approximation of convex sets by polytopes[END_REF]. It is proved in [START_REF] Bronshteyn | The approximation of convex sets by polyhedra[END_REF][START_REF] Dudley | Metric entropy of some classes of sets with differentiable boundaries[END_REF] that given a convex body Ξ P ∈ R n 0 , there exists a polytope

where d H denotes the Hausdorff distance and k(Ξ P ) is a constant only depending on Ξ P . More precise approximations are obtained in dimension 2, where we can ensure the existence of

2n sin π n where l is the length of the boundary of Ξ P . Moreover, if the boundary of Ξ P is two-times smooth, then an explicit asymptotic result is known about the distance between Ξ P and the set of circumscribed polytopes having n vertices: the closest polytope

where k(Ξ P ) is a constant depending on n 0 and the Gaussian curvature of the boundary of Ξ P [START_REF] Böröczky | Approximation of general smooth convex bodies[END_REF]. affine function λ i are the components of the ith row of D -1 , i.e., λ i (x) = n 0 j=1 l ij x j + l i(n 0 +1) , i = 1, . . . , (n 0 + 1).

(5.22)

Note that λ i (x) ≥ 0, i = 1, . . . , (n 0 + 1) iff x belongs to the n-simplex. Let σ σ , T t be the associate matrices for the operations of scaling with factor σ > 0 and translation t ∈ R n 0 . Thus the associated matrix D σ,t of the simplex with homothetic transformation on the simplex S is D σ,t = T t σ σ D, where

Its corresponding inverse is then

.

Let σ * be the smallest scaling factor σ such that a translate of σS contains Ξ P . The translate used when σ = σ * is denoted by t * .

Theorem 12 σ * and t * are given by: t * = n 0 +1 i=1 z i ω (i) and σ * = - Proof: Assume that the homothetic copy of S given by σS + t contains Ξ P . Then the coefficients λ i (x) defined in (5.22) should be nonnegative for any point x ∈ Ξ P . Considering the matrix D -1 σ,t defined above and computing the minimum values of λ i (x), i = 1, . . . , (n 0 + 1), we get

For ease of notation, we express this as l i(n 0 +1) + z i ≥ 0, i = 1, . . . , (n 0 + 1), where l i(n 0 +1) = σl i(n 0 +1) -n 0 k=1 t k l ik , i = 1, . . . , (n 0 + 1).

Since the matrix (l) i,j=1,...,(n 0 +1) is the inverse of matrix D, we have