Hybridation d’algorithmes évolutionnaires et de
méthodes d’intervalles pour l’optimisation de problèmes difficiles

Optimisation par essaims particulaires 

L’OEP [Kennedy 95] s’inspire du déplacement des groupes d’oiseaux ou des bancs de poissons. Les solutions candidates (appelées particules) se déplacent dans l’espace de recherche en modifiant leur position et leur vitesse selon des règles simples : le mouvement de chaque particule est influencé par sa meilleure position connue et par la meilleure solution connue de l’essaim dans l’espace de recherche (algorithme 3). Les positions des particules améliorées localement permettent à l’essaim de se déplacer globalement vers de bonnes solutions.

Formules de mise à jour 

La vitesse d’une particule V (k+1) à la génération k+1 est fonction de sa position courante X(k) , sa vitesse courante V (k) , sa meilleure position connue Xl et la meilleure position connue Xg de l’essaim : V (k+1) = ωV (k) + α >(Xl − X (k) ) + β >(Xg − X (k) ) (2.6) où ω est un hyperparamètre (facteur d’inertie) décidé par l’utilisateur, et α et β sont des vecteurs tirés dans [0, 1]n avec une probabilité uniforme. La position de la particule est ensuite simplement mise à jour : X (k+1) = X (k) + V (k+1) (2.7) 

Voisinage

 L’OEP standard est dite « global best », car chaque particule est attirée par la meilleure position connue Xg de l’essaim. Ce phénomène de biais est souvent la source d’une convergence prématurée vers des minima locaux. Une alternative, dite « local best », consiste à guider une particule vers la meilleure position connue d’un sous-ensemble de particules voisines (par exemple, les k plus proches).

 Algorithme à évolution différentielle 

L’algorithme à ED est parmi les AE les plus simples et les plus performants [Storn 97]. Initialement conçue pour les problèmes sans contraintes à variables continues, l’ED a été étendue aux problèmes mixtes et aux problèmes contraints. Sa robustesse et son nombre faible d’hyperparamètres en font un outil de choix pour la résolution de problèmes difficiles ; l’algorithme s’est notamment illustré dans des problèmes d’entraînement de réseaux de neurones [Slowik 08], de conception aérodynamique [Rogalsky 00], de décision multicritère, d’approximation polynomiale et d’ordonnancement de tâches. Contrairement aux AG dont les nouveaux individus proviennent des opérateurs de croisement et de mutation, l’ED n’est pas inspirée de l’évolution naturelle mais génère de nouveaux individus via des opérations géométriques. Basée sur le concept de différence de vecteurs, l’ED combine avec une certaine probabilité les composantes d’individus existants pour former de nouveaux individus (algorithme 4). Algorithme 4 Evolution différentielle fonction EvolutionDifférentielle(f : fonction objectif, NP : taille de la population, W : facteur d’amplitude, CR : taux de croisement) P ← population initiale aléatoire répéter P 0 ← ∅ . population temporaire pour x ∈ P faire (u, v, w) ← ChoixParents(x, P) y ← Croisement(x, u, v, w, W, CR) . génère un nouvel individu si f(y) < f(x) alors P 0 ← P 0 ∪ {y} . y remplace x sinon P 0 ← P 0 ∪ {x} . x est conservé fin si fin pour P ← P 0 . la population temporaire remplace la population jusqu’à critère d’arrêt vérifié renvoyer meilleur individu de P fin fonction 

 Opérateur de croisement quaternaire 

Notons NP la taille de la population, W > 0 le facteur d’amplitude et CR ∈ [0, 1] le taux de croisement. A chaque génération, NP nouveaux individus sont générés : pour chaque individu x = (x1, . . . , xn), trois autres individus u = (u1, . . . , un) (appelé individu de base), v = (v1, . . . , vn) et w = (w1, . . . , wn), tous différents et différents de x, sont aléatoirement choisis dans la population. Les composantes yi (i ∈ {1, . . . , n}) du nouvel individu y = (y1, . . . , yn) sont alors calculées de la manière suivante : yi =    ui + W × (vi − wi) si i = R ou ri < CR xi sinon (2.8) où ri est tiré uniformément dans [0, 1] et R est un indice aléatoire tiré uniformément dans {1, . . . , n} pour tout x, garantissant qu’au moins une composante de y diffère de celle de x. y remplace alors x dans la population s’il améliore la fonction objectif. L’opérateur de sélection est donc réduit à une sélection élitiste. La figure 2.5 illustre le croisement entre les individus x, u (individu de base), v et w en deux dimensions. La fonction objectif est représentée en courbes de niveaux. La différence v − w donne la direction de déplacement (une approximation de la direction opposée au gradient) dans laquelle u est translaté. minimum w v u y x1 x2 W × (v − w) v − w f x Figure 2.5 – Croisement quaternaire de l’évolution différentielle 

 Choix de l’individu de base 

Dans un AG, la probabilité de sélection d’un individu est généralement proportionnelle à son évaluation. Dans l’ED, la sélection de l’individu de base u est équiprobable parmi tous les individus de la population. Deux variantes ont été proposées par [Price 06] pour garantir que tous les individus de la population courante jouent le rôle d’individu de base 24 une et une seule fois par génération. Les indices des individus de base dans la population peuvent être obtenus : 1. par permutation aléatoire dans {1, . . . , NP} ; 2. en tirant un offset dans {1, . . . , NP −1} avec une probabilité uniforme et en l’ajoutant (modulo NP) à l’indice de chaque individu x. 

Contraintes de bornes 

Il arrive régulièrement que certaines composantes du nouvel individu y se trouvent en dehors de l’espace de recherche D. Deux alternatives sont alors possibles : – pénaliser la fonction objectif : un terme constant ou dépendant du nombre et de la magnitude des violations est ajouté à la fonction objectif. Lorsque les individus générés ont une forte tendance à violer les contraintes de bornes, cette approche présente un risque de convergence lente ; – générer une nouvelle composante à l’intérieur du domaine : saturer la borne du domaine, réinitialiser aléatoirement la composante dans le domaine entier ou sur le segment liant la composante de base ui à la borne du domaine Di (méthode du rebond [Price 06]) : yi =    ui + ω(Di − ui) si yi > Di ui + ω(Di − ui) si yi < Di (2.9) où ω est tiré avec une probabilité uniforme dans [0, 1]

Gestion directe des contraintes 

La gestion directe des contraintes consiste à séparer les évaluations de la fonction objectif et des contraintes du problème : à chaque individu est associé un tableau contenant la valeur de la fonction objectif et l’évaluation de chacune des contraintes. Le choix de l’individu de base u satisfait alors une des règles suivantes : – u appartient au domaine réalisable, contrairement à x; – u et x appartiennent au domaine réalisable, et f(u) < f(x); – u et x n’appartiennent pas au domaine réalisable, mais u ne viole aucune contrainte davantage que x. 2.5 Critères d’arrêt On distingue deux grandes catégories de critères d’arrêt : – un critère statique est généralement basé sur les ressources matérielles disponibles (temps CPU, nombre d’itérations ou d’évaluations de la fonction objectif) et connues a priori ; 25 – un critère dynamique fait référence à la qualité de la solution (suffisamment proche d’un optimum connu a priori) ou à la fin de convergence (nombre d’itérations consécutives sans améliorer la meilleure solution connue).

Analyse par intervalles 

L’analyse par intervalles est la branche de l’analyse numérique dédiée à l’encadrement des erreurs d’arrondis. Les méthodes d’intervalles permettent de calculer un minorant et un majorant rigoureux d’une fonction sur un intervalle, même en présence d’arrondis. Elles constituent donc une approche priviliégiée pour l’optimisation globale fiable. Le calcul par intervalles est décrit dans la section 3.1. Le concept de fonction d’inclusion est introduit dans la section 3.2. Les algorithmes de branch and bound par intervalles, dédiés à l’optimisation de problèmes continus, sont présentés dans la section 3.3. La section 3.4 mentionne les différents modes de différentiation automatique. 27 3.1 Calcul par intervalles L’arithmétique flottante (AF) est une méthode de représentation approchée des nombres réels sur des machines discrètes [Goldberg 91]. Un nombre flottant x est représenté par son signe, sa mantisse (la partie fractionnaire de x) et son exposant. Les floating point units ou unités de calcul en virgule flottante (FPU) présentes dans les processeurs manipulent des nombres dont la mantisse est de taille fixe, ce qui conduit à des erreurs d’approximation lorsque les réels ne sont pas représentables exactement en machine. Par exemple, la constante π arrondie à 3 décimales est soit 3.141, soit 3.142, et la valeur exacte se situe quelque part dans l’intervalle [3.141, 3.142]. La norme IEEE-754 fixant la représentation des réels et le comportement des opérations de base en virgule flottante a depuis 1985 été adoptée par la plupart des FPU. Elle spécifie quatre modes d’arrondi : au plus proche, vers zéro, vers +∞ et vers −∞. Le format IEEE double précision (64 bits) fournit une précision avancée (15 à 17 décimales). Néanmoins, l’accumulation de nombreux arrondis sur des problèmes numériquement instables est la source de résultats erronés (voir exemple 3). Exemple 3 (Accumulation des erreurs d’arrondis) Une accumulation catastrophique d’erreurs d’arrondis a été illustrée par [Rump 88]. Considérons la fonction : f(x, y) = 333.75y 6 + x 2 (11x 2 y 2 − y 6 − 121y 4 − 2) + 5.5y 8 + x 2y (3.1) En évaluant f(77617, 33096), on obtient 1.172603 en simple précision et 1.1726039400531 en double précision. Pourtant, la valeur exacte est − 54767 66192 = −0.827396. 3.1.1 Contrôle des arrondis La thèse de doctorat de [Moore 66] a posé les fondements du calcul par intervalles : l’idée de Moore est de représenter chaque résultat intermédiaire d’un calcul par un intervalle contenant le résultat exact. Un réel x représentable en machine est alors remplacé par l’intervalle dégénéré [x, x], tandis qu’un réel non représentable exactement en virgule flottante (par exemple 0.1) est encadré rigoureusement par un intervalle. L’arithmétique d’intervalles (AI) étend l’arithmétique réelle aux intervalles de manière rigoureuse. Les fonctions élémentaires (+, −, ×, /, log, exp, etc.) sont implémentées en arrondissant les calculs vers l’extérieur : la borne gauche du résultat est calculée avec un mode d’arrondi vers −∞ et la borne droite avec un mode d’arrondi vers +∞. La valeur exacte est alors numériquement garantie d’appartenir à l’intervalle obtenu par AI. De nombreuses bibliothèques logicielles implémentent l’AI : Profil/BIAS [Knüppel 94] (développée en C++ à l’Université de Technologie de Hamburg), Gaol [Goualard 03] (implémentation C++ d’opérateurs de programmation par contraintes sur intervalles), Boost [Brönnimann 06] (template C++), MPFI [Revol 02] (bibliothèque multi-précision en C et C++), Sun [Microsystems 01] (Fortran 95 et C++) et Filib [Lerch 01]. Nous avons récemment implémenté une bibliothèque de calcul par intervalles dans le langage fonctionnel OCaml [Alliot 12b]. Par souci de performance, des routines de bas niveau (C et assembleur) permettent un contrôle fin des arrondis. 

Arithmétique d’intervalles 

Nous adoptons les notations suivantes : – R dénote l’ensemble des réels ; – F dénote l’ensemble des flottants ; – un intervalle X = [X, X] à bornes flottantes définit l’ensemble : X := {x ∈ R | X ≤ x ≤ X} (3.2) – un intervalle est dit dégénéré lorsque X = X ; – I représente l’ensemble des intervalles à bornes flottantes : I := {[X, X] | (X, X) ∈ F 2 ∧ X ≤ X} (3.3) – pour un ensemble D ⊂ R, I(D) dénote l’ensemble des intervalles inclus dans D. La définition s’étend au cas multivarié ; – int(X) désigne l’intérieur d’un intervalle X non dégénéré, c’est-à-dire l’ensemble : int(X) := {x ∈ R | X < x < X} (3.4) – m(X) := 1 2 (X + X) désigne le milieu de l’intervalle X ; – w(X) := X − X désigne la largeur de l’intervalle X ; – une boîte X = (X1, . . . , Xn) est un produit cartésien d’intervalles ; – m(X) := (m(X1), . . . , m(Xn)) désigne le milieu de la boîte X ; – w(X) = maxi=1…n w(Xi) désigne la largeur de la boîte X ; – (X, Y ) désigne l’enveloppe convexe de X et Y , soit le plus petit intervalle de I contenant X et Y , et calculable par la machine. Par la suite, nous notons en majuscules les quantités intervalles et en gras les vecteurs. Un intervalle est donc noté X, une boîte X et un vecteur de réels x. L’extension à l’AI des opérateurs binaires  ∈ {+, −, ×, /} fournit le plus petit intervalle contenant l’image directe : X  Y = {x  y | x ∈ X ∧ y ∈ Y } (3.5) On peut calculer explicitement les bornes gauche et droite du résultat en fonction des bornes des opérandes : [a, b] + [c, d] =[a + c, b + d] [a, b] − [c, d] =[a − d, b + c] [a, b] × [c, d] =[min(ac, ad, bc, bd), max(ac, ad, bc, bd)] 1 [a, b] =[1 b , 1 a ] si 0 ∈/ [a, b] [a, b] [c, d] =[a, b] × 1 [c, d] si 0 ∈/ [c, d] (3.6) 29 Rappelons que les calculs sur intervalles (par exemple l’addition notée de manière abusive [a, b] + [c, d] = [a + c, b + d]) doivent être effectués en arrondissant la borne gauche vers −∞ et la borne droite vers +∞. La plupart des fonctions élémentaires unaires sont étendues aux intervalles en utilisant leurs propriétés de monotonie : exp([a, b]) = [exp(a), exp(b)] log([a, b]) = [log(a), log(b)] (a > 0) » [a, b] = [√ a, √ b] (a ≥ 0) (3.7) Les fonctions non monotones sur l’intervalle considéré (puissances paires, fonctions trigonométriques) sont néanmoins monotones par morceaux et font l’objet d’une étude au cas par cas. L’AI possède des propriétés plus faibles que l’arithmétique réelle : la soustraction n’est pas l’opération inverse de l’addition sur I, de même que la division n’est pas l’opération inverse de la multiplication (exemple 4). En outre, la propriété de distributivité de la multiplication par rapport à l’addition n’est pas vérifiée sur I; seule la propriété plus faible de sous-distributivité est valable : ∀(X, Y, Z) ∈ I 3 , X(Y + Z) ⊂ XY + XZ (3.8) Exemple 4 (Non-inversibilité de l’addition et de la multiplication) Soient X = [2, 4] et Y = [3, 5]. Alors : W = X + Y = [2, 4] + [3, 5] = [5, 9] W − X = [5, 9] − [2, 4] = [1, 7] ) Y Z = X × Y = [2, 4] × [3, 5] = [6, 20] Z X = [6, 20] [2, 4] = [1.5, 10] ) Y (3.9) L’AI étendue [Hanson 68, Kahan 68, Hansen 92] généralise l’AI aux intervalles à bornes infinies, ce qui permet notamment de décrire la division par un intervalle contenant zéro : 1 [a, b] =    ∅ si a = b = 0 [ 1 b , 1 a ] si 0 < a ou b < 0 [ 1 b , +∞] si a = 0 [−∞, 1 a ] si b = 0 [−∞, 1 a ] ∪ [ 1 b , +∞] si a < 0 < b (3.10) ou les calculs à bornes infinies : » [4, +∞] = [2, +∞] (3.11) 30 3.2 Fonctions d’inclusion Les propriétés de conservativité du calcul par intervalles permettent de construire des extensions aux intervalles (définition 12) des fonctions réelles décomposables (définition 10). Définition 10 (Fonction décomposable) Une fonction est décomposable (factorable function) lorsqu’elle s’exprime récursivement comme une composition de termes élémentaires (opérateurs, fonctions élémentaires ou variables). Définition 11 (Image directe) Soient f : D ⊂ R n → R et X ∈ I(D). On note f(X) l’image directe de f sur X : f(X) := {f(x) | x ∈ X} (3.12) Définition 12 (Fonction d’inclusion) Soient f : D ⊂ R n → R et F : I(D) → I. On dit que F est une fonction d’inclusion (ou extension aux intervalles) de f si : ∀X ⊂ I(D), f(X) ⊂ F(X) (Conservativité) ∀(X,Y ) ∈ I(D) 2 , X ⊂ Y =⇒ F(X) ⊂ F(Y ) (Isotonicité d’inclusion) (3.13) Les fonctions d’inclusion d’une fonction réelle ne sont pas uniques : on peut généralement construire des fonctions d’inclusion dont les ordres de convergence (définition 13) sont différents, c’est-à-dire dont la surestimation décroît à des vitesses différentes avec la taille de l’intervalle. L’évaluation d’expressions syntaxiquement équivalentes en arithmétique réelle peut fournir des encadrements plus ou moins précis en AI. Ce point est discuté dans la section 3.2.1. Définition 13 (Ordre de convergence) Soient f : D ⊂ R n → R et F une fonction d’inclusion de f. On dit que F est d’ordre de convergence α (α > 0) lorsque, pour tout X ∈ I(D) : w(F(X)) − w(f(X)) = O(w(X) α ) (3.14) w(F(X)) − w(f(X)) représente l’erreur de surestimation de l’image de f sur X. L’extension aux intervalles la plus simple à construire pour une fonction décomposable est la fonction d’inclusion naturelle (définition 14). Elle a un ordre de convergence linéaire, c’est-à-dire que l’erreur de surestimation varie linéairement avec la taille des domaines des variables. Définition 14 (Fonction d’inclusion naturelle) Soit f : D ⊂ R n → R. Sa fonction d’inclusion naturelle FN : I(D) → I est obtenue en remplaçant dans l’expression de f chacune des variables par son domaine et chacune des opérations élémentaires par son extension aux intervalles.