N° d'ordre:
lBS
THE8E
présentée li
L'UNIVERSITE PAUL SABATIER DE TOULOUSE (Sciences)
pour obtenir le
DOCTORAT de L'UNIVERSITE PAUL SABATIER
Spécialité: Calcul Scientifique
par
Sou mana
Beidi
HAMMA
ETUDE DE METHODES NUMERIQUES
D'OPTIMISATION GLOBALE
soutenue publiquement le 13 Mars 1992, devant le jury composé de:
BONNEMOY C.
Professeur. Université
d'Auvergne.
Clermont-Ferrand
HIRIART-URRUTY J.-B.
Professeur. Université Paul Sabatier de Toulouse
MALIVERT Ch.
Professeur. Université de Limoges
PHAM DINH T.
Professeur. I.N.S.A.
de Rouen
RAYMOND J.-P.
Professeur. Université Paul Sabatier de Toulouse
ROUZE M.
Chef du département Projet et Analyse C.N.E.S.•
Centre Spatial Guyanais. Kourou.
Laboratoire
d'Analyse
Numérique
Formation
Doctorale de Mathématiques Appliquées
Université
Paul
Sabatier
Travail réalisé dans le cadre d'un marché C.N.E.S.
"Etude de problèmes d'optimisation globale"
(marché N° 873 - CNES - 87 - 4937)
REMERCIEMENTS
••••••••••••••••••••••••••
C'est en sortant de l'EN.A.C. (Ecole Nationale de l'Aviation Civile) de Toulouse que
j'ai eu la chance de rencontrer le Professeur Jean-Baptiste HIRIART-URRUTY. Il m'a fait
reprendre le goût des mathématiques, m'a initié à la recherche et a dirigé mes recherches durant
tourie travail de thëse. Pour tous ses conseils, toute son aide et sa disponiblùé régulière je lui
suis infiniment reconnaissant de m'avoir guidé.
Le département Mathématiques Appliquées et Analyse Numérique du Centre National
d'Etudes Spatiales (CN.E.S.) de Toulouse est à l'origine de ce travail. Je remercie donc tous
les membres du département, en particulier Mmes C. AUBERT, P. DORIO et MM. F.
BONNEAU, G. LASSAllE-BALlER et
P. LEGENDRE.
Mr. ROUZE a été notre premier contact avec le département. Il a suivi ce travail de trés
près. Son aide et ses conseils ont été précieux. De plus il a accepté de faire le déplacement
(depuis Kourou) pour participer à mon jury. Pour tout cela je le remercie vivement.
Avec le Professeur C. BONNEMOY de l'Université d'Auvergne on a eu une
coopération trës fructueuse.ll a aussi accepté de présider mon jury. Je l'en remercie beaucoup.
Je remercie également les Professeurs C. MAL/VERT de l'Université de Limoges,
PHAM DINH Tao de l'IN.S.A. de ROUEN, J.-P. RAYMOND de l'Université Paul Sabatier
de Toulouse pour avoir accepté d'être les rapporteurs et examinateurs de ce travail et pour leurs
corrections et critiques constructives.
J'adresse ma reconnaissance à mes deux "nègres" en informatique: l'un Blanc, Mr K.
YAZDANIAN ( C.E.R.T. et Sup'Aéro de Toulouse), l'autre ... Noir, J.-C. HOCI/ON
ingénieur en Informatique de l'E.N.S.E.E.I.H.T. de Toulouse.
Je remercie aussi le Département de Mathématiques de la Faculté des Sciences de
Limoges, qui en m'accueillant depuis Octobre 1991 m'a permis de terminer la rédaction de ce
travail et de bénéficier des conseils et de /' expérience de son équipe d'optimisation dirigée par
le Professeur M. THERA.
Que tous les membres (chercheurs, techniciens et administratifs) du Laboratoire
d'Analyse Numérique
de L'université Paul Sabatier de Toulouse trouvent ici ma
reconnaissance pour m'avoir aidé a accomplir ce travail.
Je ne terminerai pas celle liste de gratitude sans remercier Thierry & Marie BELTAN et
tous mes amis Baha'is pour leur amitésincère et leur soutien moral.
Par la possesion du don idéal de la recherche scientifique,
l'homme est le produit
le plus noble de la création, le
régulateur de la nature,., Ce don est le pouvoir le plus
louable de l'homme car, par l'utilisation qui en est faite,
l'amélioration
de
la
race
humaine
s'accomplit,
et
le
développement
des
vertus
de
l'humanité
est
rendu
possible...
Abdu'l-bahà
(1844-1921)
"Les bases de J'unité du monde"
Certaines personnes dépourvues de génie personnel
sont quelquefois douées du pouvoir de le stimuler...
Anonyme
Baba,
Maman,
Permettez-moi de dédier ce travail
A celui qui, en maîtrisant la tête du serpent,
M'a permis de jouer avec la queue ...
A mon tour, à présent, de prendre
La tête du serpent pour nos enfants;
Même si cette fois-ci,
11 s'agit plutôt ... d'une pieuvre!
PRESENTATION GENERALE
Le
présent
travail
qui
porte
sur
les
méthodes
numériques
d'optimisation globale, a pour origine un sujet proposé par le C.N .E.S.
(Centre National d'Etudes Spatiales) de Toulouse "Etude de problèmes
d'optimisation globale" (marché N° 873-CNES-87-4937). Le but était de
développer des codes d'optimisation globale en vue d'applications à des
problèmes
spécifiques.
Au chapitre l , nous rappelons les problématiques de l'optimisation
locale et de l'optimisation globale d'un point de vue numérique.
Nous
nous
attachons, au chapitre 2, à faire
l'état de l'art des
méthodes et codes existants en optimisation globale, en
les
classant
suivant différents critères: leur philosophie, leur stratégie, leur garantie
de convergence, etc.
Ensuite, nous étudions de près, au chapitre 3, certaines méthodes
numériques:
le Recuit Simulé, la méthode du Tunnelier, le Polytope
Mouvant.
Dans
ce
chapitre,
nous
présentons
aussi
des
résultats
numériques de codes que
nous
avons écrits, testés et
améliorés. Ce
chapitre contient également, à la suite de
plusieurs expérimentations
numériques, des détails de choix de paramètres, conseils d'utilisation,
etc. pour certains algorithmes.
Enfin nous abordons, au chapitre 4, le problème d'optimisation de
trajectoires interplanétaires tel qu'il nous a été présenté par le C.N.E.S.
A la fin de chaque chapitre (ou
paragraphe indépendant) nous
donnons les références bibliographiques qui s'y rapportent.
En
fin
de
document,
nous
proposons
une
(proposition
de)
traduction en Français de différents termes anglais couramment utilisés
en optimisation globale numérique. C'est la traduction que nous avons
utilisée dans ce mémoire.
Une annexe propose les programmes et les détails de
plusieurs
résultats
numériques
expérimentaux.
TABLES DES MATIERES
INTRODUCTION GENERALE
4
Chapitre
1
OPTIMISATION LOCALE vs OPTIMISATION GLOBALE
1.1.
GENERALITES
7
1.1.1.
Convention de vocabulaire
8
1.1.2.
Position des Problèmes
9
1.1.3.
Solutions numériques
Il
1.1.4.
Diversités des méthodes
Il
1.1.5.
Difficultés rencontrées pour notre étude
13
1.2.
OPTIMISATION
LOCALE
1.2.1. Rappels sur la théorie de l'optimisation locale
14
1.2.2. Quelques exemples d'algorithmes efficaces
16
1.2.2.1. l'algorithme de D.F.P.
17
1.2.2.2. les algorithmes de Quasi-Newton
18
1.2.3. Transition vers l'optimisation globale
19
référen ces
20
Chapitre
2
METHODES NUMERIQUES D'OPTIMISATION GLOBALE
2.0.
Rappel
du
problème
22
2.1. Classification
des
problèmes
23
2.1.1. Problèmes avec/ou sans contraintes.
23
2.1.2.
Problèmes particuliers/généraux
23
2.2. Classification
des
méthodes
24
2.2.1. Selon l'approche
25
2.2.2. Selon la philosophie
26
2.2.3. Selon la garantie de convergence
27
2.2.4. Autres classifications
28
2.3. Choix
d'une
classirication
29
2.3.1. Méthodes avec garantie de convergence
30
2.3.2. Méthodes directes
31
2.3.3. Méthodes indirectes
38
2.4. Caractéritiques
générales
des
méthodes
39
2.4.1. Optimisation sur un compact S
39
2.4.2. Schémas généraux
40
2.5. Difficultés
Numériques
2.5.1. Critères d'arrêt - Convergence des codes
44
2.5.2. Difficultés liées aux contraintes
45
2.5.3. Propriétés minimales pour un code
d'optimisation globale
45
2.5.4. Caractérisation d'un minimum global
46
2.6. Parallélisation
48
2.6.1. Motivations
48
2.6.2. Machines parallèles
48
2.6.3. Parallélisation
49
2.6.4. Conclusion
50
2.7.
Quelques
applications
de
l'optimisation
globale
50
2.7.1. Les moindres carrés
51
2.7.2. Autres exemples
53
2.8.
Conlusion
53
références
55
Chapitre
3
EXEMPLES D'ALGORITHMES NUMERIQUES
D'OPTIMISATION GLOBALE
3.1.
LE
RECUIT
SIMULE en
optimisation
"continue" 59
3.1.1. Philosophie
60
3.1.2. Algorithme général du Recuit Simulé
61
3.1.3. Application à l'optimisation dans Rn
64
3.1.4. Améliorations du Recuit: utilisation de cycles de
Réchauffement
69
3.1.5. Résultats Numériques
81
3.1.6. Conclusion
91
références
92
3.2.
LE "POLYTOPE MOUVANT" OU SIMPLEXE DE NELDER
& MEAD
94
3.2.1
Principe
95
3.2.2 Algorithmes du Polytope Mouvant
96
3.2.2.1. Algorithme séquentiel
98
3.2.2.2. L'algorithme parallélisable
99
3.2.2.3. Description détaillée de la version séquentielle
101
3.2.3. L'algorithme à recherche adapté
107
3.2.4. Convergence de la méthode du Polytope Mouvant
110
3.2.5 Résultats Numériques
113
3.2.6. Conclusion
115
références
1 15
3.3.
LA METHODE DU TUNNELIER
1 17
3.3.1. Principe
1 17
3.3.2. Algorithme du Tunnelier
1 18
3.3.3. Choix des paramètres
1 2 1
3.3.4. Convergence de la méthode du Tunnelier
123
3.3.4. Résultats Numériques
123
références
124
Chapitre
4
OPTIMISATION DE TRAJECTOIRES INTERPLANETAIRES
4.0.
Introduction
126
4.1.
Présentation
de
la
mission
VESTA
127
4.2.
Notations
128
4.3.
Position
du
problème
129
4.4.
Résolution
du
problème
129
4.5.
Le
problème
d'optimisation
13 1
4.6. Calcul de la fonction coût
132
4.7. Code général HOPTIMA
136
4.8.
Résultats
numériques
140
4.8.1 Trajectoires initiales
140
4.8.2 Résultats
143
4.8.2.1. Exemples de résultats complets
d'optimisation locale
143
4.8.2.2. Résultats d'optimisation globale
147
4.9
Conclusion
149
références
150
Conclusion
générale
1 5 1
Glossaire
153
4
INTRODUCTION
GENERALE
Soit f: Rn --+ R une fonction objectif, continue.
La
programmation
non
linéaire ou optimisation
(locale) non
linéaire
consiste en
un ensemble de méthodes pour trouver un optimum local
(minimum ou maximum). Mais comme maximiser f revient à minimiser
-f,
nous considérerons donc
tout
problème d'optimisation comme
un
problème de minimisation c'est-à-dire
(Ploc):
trouver un point Xloc de Rn tel qu'il existe un voisinage V de
Xloc avec
f(Xloc) 1 f(X)
V X eV
(R.O.I.)
Cependant en
pratique,
il existe, en
général,
plusieurs
minimiseurs
locaux avec des valeurs de fonction différentes. Ainsi le problème de
l'optimisation globale se veut de :
(Pg):
trouver un point x* de Rn tel que f(X*) 1 f(X)
V X e Rn
(R.O.2)
Nous nous intéressons dans ce travail à la résolution numérique de ce
problème.
Pour des
raisons numériques (Le.
informatiques cf 1.1.2.) on suppose
l'existence d'un compact S
contenant x* comme point intérieur, et le
problème (Pg) devient celui que nous allons considérer pour toute la
suite, soit:
1 (P): trouver
x- de S tel que
f(X*) 1 f(x)
V X e S
(R.O.3)
L'optimisation
globale est
restée longtemps
marginale même
à
l'intérieur
du
domaine
d'activités
ou
de
recherches
dénommé
Optimisation
ou Programmation
Mathématique. II y a plusieurs
raisons à cela, l'une étant le caractère très difficile du problème. En effet,
contrairement à l'optimisation locale, on
ne dispose
pas
d'outils de
l'Analyse mathématique permettant une caractérisation "nette" (ou une
localisation) d'un
minimiseur global. Et
si une
telle
caractérisation
existe, elle doit être destinée, alors, à une implémentation numérique, et
donc
les
preuves
(théorèmes) d'optimalité
globale doivent être
non
seulement
théoriques
mais
aussi
numériques,
ce
qui
complique
davantage les choses.
Depuis quelques années cependant, on assiste à un regain d'intérêt
pour l'optimisation globale et cela est dû aux demandes de plus en plus
précises et
pressantes d'utilisateurs d'origines
très
variées. En
effet,
beaucoup de problèmes pratiques d'ingénierie se formulent comme des
problèmes d'optimisation globale (P); c'est le cas notamment en chimie,
5
biologie,
économétrie,
design
électrique,
aéronautique,
statistiques,
énergétique, ... (cf [FL090], [TOR89]). Cela suffit pour illustrer l'enjeu du
problème d'optimisation
globale car,
quand
un
utilisateur
se
rend
compte qu'il trouve des minimums locaux différents en exécutant son
code d'optimisation (locale), il est tout à fait légitime qu'Il songe à un
code qui serait capable d'en tenir compte et de donner le meilleur de
tous.
Beaucoup
d'articles,
publiés
souvent
par
des
ingénieurs,
ont
présenté des
problèmes spécifiques à leur domaine, des
propositions
d'algorithmes,
des
résultats
numériques,
inspirant
ainsi
la
théorie
mathématique de l'optimisation globale.
Jusqu'en 1978, en dehors de ces articles, peu de congrès ou livres-
revues ont fait la synthèse de l'optimisation globale, mis à part [DIXON &
SZEGO 1975] et [DIXON & SZEGO 1978]. Mais depuis, les publications,
attestant
l'intérêt
pour
cette
branche
de
l'optimisation,
se
sont
multipliées. Ceci est attesté par des "surveys" de [ARCHEITI &
SCHOEN
1984], [RINNOOY KAN & TIMMER 1986], [COLLINS & al 1988], une thèse
[TiMMER
1984]
sur
les
méthodes
stochastiques,
et
des
livres
[LAARHOVEN & al 1987] sur le Recuit, [PARDALOS & ROSEN 1987] sur
l'optimisation quadratique, [TORN & ZILINSKAS 1989] faisant l'état de
l'art sur les méthodes numériques, et [HORST & TUY
1990] sur les
méthodes déterministes. Enfin la revue nouvelle née en 1991, "Journal
of Global Optimization" qui est déjà à son numéro 3 témoigne de la
détermination
du
milieu
scientifique
à
donner
un
coup
de
pouce
supplémentaire à cette branche de l'optimisation qui, sans aucun doute,
aura des développements dans l'avenir.
Cependant, pour l'heure, force est de constater qu'on ne dispose
pas de codes d'optimisation globale "prêt-à-porter". Beaucoup de codes
sont, certes, proposés par tel ou tel auteur mais
peu ont
un champ
d'applications large. Ils
sont
souvent conçus pour résoudre de
façon
ponctuelle tel
ou
tel
problème spécifique. Et on
verra probablement
dans les années à venir, grâce aux développements plus ambitieux de la
théorie, des propositions de codes de plus en plus performants.
Dans le travail qui suit, nous verrons comment une classification
des
méthodes,
peut
servir
de
base
pour
développer
des
codes
appropriés
basés sur telle ou
telle
méthode.
Et c'est pourquoi
nous
soulignons
pour
chaque
classe
de
méthodes
sa
philosophie,
ses
hypothèses, et éventuellement ses forces et faiblesses.
Nous présentons les détails de l'implémentation de trois algorithmes:
le premier basé sur une méthode stochastique: le Recuit, les deux autres
basés sur deux méthodes déterministes, du Tunnelier et du
Polytope
Mouvant.
6
Chapitre
1
OPTIMISATION LOCALE
vs
OPTIMISATION GLOBALE
7
1.1.
Généralités
Tout le
travail qui
suit a été
essentiellement centré
sur
l'aspect "al~orithmes numériques" de l'optimisation globale. Soit donc la
résolution
numérique du problème suivant :
(P)
Minimiser f(x) où f: R n ~ R
continue.
DUjpjtjop
J J J
On dira que (Xloc.f1oc) de Rn X R. où f1oc=f(Xloc).
est un minimum local si et seulement si:
3
V un voisinage de Xloc tel que:
V X e V
f(Xloc) S f(X)
(R.L1.I.)
DUjpjljop
1.1.2
On dira que (Xg,fg) de Rn X R, où fg=f(Xg),
est un minimum global si et seulement si:
V X e Rn
f(Xg) S f(X)
(R. 1. 1.2.)
y
floc
Fg
x
Xloc
Xg
fi~ure 1.1
Xloc minimiseur local, Xg minimiseur global pour
f: R ~ R
8
1.1.1.
Convention
de
Vocabulaire
A
cause
de
la
diversité
des
vocabulaires
et
des
confusions
éventuelles que cela pourrait engendrer, nous allons convenir d'adopter
le vocabulaire qui suit:
• minimum,
maximum,
optimum
nous allons
utiliser ces trois mots et
leurs pluriels,
minimums,
maximums et optimums.
• point de
minimum
ou (état
de
minimum)
on appellera point de minimum le couple (X*, f*)
solution de (P) avec
X* E Rn
et
f* =f( X*).
• minimiseur,
maxrrmse ur,
optimiseur
on désignera par minimiseur le point X* de Rn
solution de (P).
• on appellera le min i m u m de f la valeur réelle f* . Cela nous
permettra certaines fois (cas où f est constante sur S, cas où
il y a
plusieurs minimiseurs mais avec la même valeur de fonction ... ) de dire
qu'il y a un seul minimum; sinon, les algorithmes qui prétendent passer
par tous les "minimums", s'il s'agissait des couples (X* , f(X*)), n'auraient
aucun sens car il y en aurait une infinité!
• solution
de
(P)
On désignera indifféremment solution de (P) selon les circonstances, ou
un
minimiseur, ou
un minimum ou
un point de
minimum car, bien
entendu, les trois sont équivalents.
• algorithme
général
Ce
terme
désignera
un
résumé
(en
blocs
de
programmation)
d'une
méthode,
d'une
technique ou
d'un algorithme.
On
ne
prétend
pas
y
donner
les
détails
nécessaires
pour
une
"finition"
en
programme
informatique, mais on y met en évidence la philosophie, les grandes
étapes, tout cela restant très souple de façon à résumer l'essentiel des
situations
d'implémentations
informatiques
possibles.
Chaque
fois
qu'une
théorie
mathématique
élaborée
est
liée
à
un
algorithme
on
préfèrera
le
terme
"méthode de ... " bien
que
dans
la
plupart des cas, dans notre discours, méthode numérique et algorithme
général soient confondus.
9
1.1.2.
Position
des
problèmes
Nous donnons ici dans sa généralité le problème (P) d'optimisation
globale qui nous intéresse. Pour mieux apprécier son importance, nous
formulons aussi le problème (Ploc) d'optimisation locale.
Soit
f: R n ~ R continue;
Nécessité
d'un
compact
S
contenant
un
minimiseur
comme
point
intérieur:
Minimiser f numériquement nécessite la connajssane à priori d..:.u..n
compact S de Rn sur lequel on ramènerait le problème; et ceci, que le
problème
soit
avec
ou
sans
contraintes.
En
effet,
l'optimisation
numérique (surtout) globale a besoin d'une
limitation du
champ de
recherche,
dit
aussi
de
faisabilité,
pour
plusieurs
raisons.
D'abord
numériquement on ne peut ni représenter ni parcourir convenablement
Rn de façon que l'algorithme "visite"
facilement un minimiseur. Surtout,
s'il faut en plus que l'algorithme converge en un temps convenable. Car,
comme
nous
le
verrons
au
chapitre
2,
beaucoup
de
méthodes
d'optimisation
nécessitent une partition et/ou
une discrétisation assez
convenables du domaine de recherche de façon à pouvoir le représenter
en
un certain nombre fini de points "significatifs" permettant, après
application
d'autres
procédures,
de
"visiter"
un
minimiseur
global.
D'autres méthodes ont
besoin, comme hypothèse, que
le nombre des
minimums locaux (il s'agit ici des valeurs minimales) soit fini. Et surtout,
on a besoin d'une condition suffisante d'existence d'un minimum global!
Pour toutes ces raisons et d'autres, que le problème soit avec ou sans
contraintes,
on
suppose
l'existence d'un
compact
S, contenant une
solution
comme
point
intérieur
sur
lequel
on
applique
le
code
d'optimisation choisi (cf 2.4.1). Bien sûr dans les cas où le problème est
avec contraintes, on s'assure que, sur le compact S toutes les contraintes
sont vérifiées.
En pratique ce compact S est donné, a priori par excès, par l'utilisateur
qui
connait à peu
près
un domaine du
faisable
et/ou
un domaine
contenant un minimiseur.
Problème (Ploc) d'optimisation locale
Minimiser(loc) f(X)
(Ploc)
1x e S
ou encore
10
Trouver (Xloc , floc) de Rn X R, où floc=f(Xloc),
(Ploc)
tel que:
'V X e V, f(Xloc) 5 f(X)
où V est un voisinage de Xloc.
Problème CP) d' optimisation elobale
(P)
IMinimiser(global) f(X)
x e S
ou encore
Trouver (Xopt, fopt) de Rn X R, où fopt=f(Xopt),
(P)
tel que:
1'V X e S , f(Xopt) 5 f(X)
y
y =f(X)
Fg
x
X=a
\\"'X=b
\\
Xg=b
figure1.2:
Xg minimiseur global pour f: R -t R
avec contraintes:
a i Xi b
Il
1.1.3.
Solutions
numériques
de ( P)
Théoriquement, résoudre (P) reviendrait à déterminer exactement
le
plus petit de tous les minimums locaux (ou un de ses arguments X·
au cas où
il y aurait plusieurs minimiseurs globaux). Cependant, d'un
point de vue numérique, une solution ne pourrait être obtenue qu'avec
une certaine précision donnée. Dès lors, le problème (P) sera considéré
résolu dès qu'on aura atteint un élément, pour un seuil E donné, d'un au
moins des ensembles Ax*(E). ACx*(E). ADx*(E) suivants (où X· est un
minimiseur global) :
• Ax*(E) = {x e S tel que IIx. x*1I < E}
• Arx*( E) = lx e S tel que Ir(x) - r(x*)1 < E}
• ADx*(E) = {x e S tel que D(r(x» , E}
• D()
m ( {z e S tel-.!ue r(z) , y}
ou
y -
m(S)
avec m mesure de Lebesgue.
1.1.4. Diversité
des
méthodes
Durant ces dernières années, l'optimisation globale a fait l'objet de
plusieurs études. Des articles ont été publiés sur le sujet, traduisant la
richesse des approches et des motivations.
Contrairement à ce que l'on pourrait être tenté de croire, l'optimisation
globale numérique n'a pas hérité (du moins pas
systématiquement) des
facilités des techniques des méthodes numériques d'optimisation locale.
En
effet,
ces
dernières
utilisent
pour
la
plupart
des
directions
de
descentes
(conditionnées
par
des
calculs
de
gradients
ou
de
leurs
approximations),
ce
qui
permet
de
converger
naturellement
vers
un
point de minimum local. Or, justement, l'optimisation globale, évite de
rester en de tels points. Il lui faudrait. au contraire, en échapper.
C'est pourquoi beaucoup d'approches ont été utilisées attestant (une fois
de plus!) de toute l'ingéniosité du génie humain, dans la tentative de
résoudre le problème. Car il s'agit, de façon résumée, de trouver un état
de minimum et de n'y rester que s'il est le meilleur (global).
Pour illustration de cette diversité, citons (cf ch 2):
• la méthode d'énumération
des
minimums
locaux
(local minima enumeration)
12
on
passe par tous les
rmmmurns (si cela est faisable!). On
prend le
meilleur qui est donc solution de (P).
• la méthode des initialisations
multiples
(muftistart)
on utilise plusieurs fois une technique d'optimisation locale en différents
points de S; on prend le meilleur résultat; il serait solution de (P).
• la méthode du Tunnelier
(tunneling)
on cherche (par une technique locale) un
minimum local XI; on utilise
une technique (du Tunnelier) pour éliminer tous les points qui ont une
valeur supérieure
à celle de XI (en effet ils ne pourraient être optimum
global). Puis on réoptimise pour trouver un XI+ 1 et ... ainsi de suite. La
suite des Xk ainsi construite converge vers le minimiseur global.
• les méthodes à Recherche
Aléatoire
(Random
Seareh)
on
construit une
suite (Xk)k
à
valeurs
de
fonction
décroissantes.
convergente vers
un minimiseur global en
passant de
X k à Xk+ 1; la
transition
se
faisant
grâce à des
processus
probabilistes de
façon
à
garantir (en probabilité) la "visite" d'un mimiseur global.
• la méthode du
Recuit Simulé (Simu/ated
Annea/ing)
C'est une
variante du
Random Search où la suite (Xk)k n'est pas à
valeurs de
fonction décroissantes. En effet. on accepte (sous certaines
conditions) de passer de Xk à Xk+1 même quand
f(Xk+l) ~
f(Xk). Ceci,
afin de permettre "d'al/er voir à côté si ça ne serait pas mieux". En effet.
comme on
le
verra
plus
loin, cela
permet d'échapper aux
états
de
minimums locaux et de converger vers un état de minimum global.
• la méthode des
montagnes
russes
(Steepest Ascent Descent : SAD.)
D'abord on cherche un minimum (local). Pour ne pas y rester. on cherche
cette fois-ci un maximum (local); cela permet d'échapper de
l'état de
minimum
local.
Puis
on
refait
une
minimisation
suivie
d'une
autre
maximisation et... ainsi de suite. Le plus petit minimum trouvé. à la fin,
est candidat pour être le global.
• les algorithmes
Génétiques (Genetic
A/gorithms)
Ces algorithmes sont basés tout simplement sur la théorie de l'évolution
selon DARWIN. On engendre successivement des populations NI. N2 .... et
à chaque génération. on applique la loi de l'évolution; du
moins une
adaptation
informatique:
limitation de
la population totale à un
seuil
donné par exemple. Ce qui entraîne l'élimination de certains individus
13
("inutiles"). La suite des points "sélectionnés par la nature" contient un
minimiseur global.
1.1.5. Difficultés
de
notre
étude
Une fois que nous nous sommes intéressés à l'aspect numérique de
l'optimisation
globale,
nous
nous
sommes
trouvés
confrontés
aux
problèmes
suivants:
1 -
il Y a beaucoup d'articles sur le sujet
2 -
il n'y a pas encore de code numérique "prêt-à-porter"
3 -
il faut se fixer le genre de problèmes à traiter
4 -
aucune bibliothèque scientifique ne propose de code
Cela est dû en grande partie aux difficultés rencontrées dans la
tentative
de
résoudre
numériquement
(P).
Ces
difficultés,
nous
les
discuterons, au fur et à mesure, dans les chapitres suivants.
Disons simplement qu'un bon nombre d'articles ont été publiés sur le
sujet,
plusieurs
méthodes
sont
proposées
mais
peu
de
codes
sont
disponibles. Et
souvent ces
codes
sont de portée limitée : cas
des
fonctions
L-Iipschiziennes,
cas
des
problèmes
D.C.
(différences
de
fonctions convexes), problèmes sans contraintes, avec peu de minimums
locaux, ensembles C de contraintes particuliers (polyèdres, de mesure
non nulle), etc ...
Pour notre part, après une recherche bibliographique qui nous a
fourni près de 600 références, nous nous sommes fixés comme objectif
de programmer un code qui utiliserait peu d'informations sur la fonction
f à minimiser, et qui soit basé sur une méthode stochastique. Car notre
but était de
l'appliquer à la résolution de
problèmes d'optimisation
présentés par le C.N.E.S. C'est ainsi que nous avons choisi la méthode du
Recuit Simulé.
Pour valider les codes issus de l'étude du Recuit, nous avons eu besoin
de
les
comparer
à d'autres.
Ce
qui
nous
a
valu
d'étudier et
de
programmer quelques codes basés sur d'autres méthodes (déterministes
par exemple).
14
1.2. OPTIMISAnON LOCALE
Bien que l'optimisation locale ne soit pas le thème essentiel de
notre travail. nous en présentons les grandes lignes.
En
effet,
comme
nous
le
verrons
plus
loin
(cf.
ch
2)
beaucoup
d'algorithmes d'optimisation globale utilisent des
codes d' optimisation
locale comme sous-codes. Donc l'efficacité de ces derniers affecte celle
des
premiers.
Il est utile donc de montrer qu'i1 y a dans le domaine de l'optimisation
locale
des
algorithmes
efficaces.
Comme
exemples
d 'algorithmes
d'optimisation globale utilisant des codes d'optimisation locale, citons les
algorithmes des Initialisations Multiples.
l'algorithme du Tunnelier (cf.
ch 3).
1.2.1.
Rappels
sur
la
théorie
de
l'optimisation
locale
POSITION DU PROBLEME
Soit
f: Rn -+
R continue,
On a vu au §.1.I.2
que. numériquement. pour minimiser f (que ce soit
avec
ou
sans contraintes). on
se
ramène à
la
minimisation
sur
un
compact S de Rn à cause de certaines difficultés informatiques.
Donc nous reformulons le problème d'optimisation locale sur S.
Trouver (Xloc • floc) de S X R. où f1oc=f(Xloc),
tel que:
(Ploc)
'<1 X eV. f(Xloc) s. f(X)
où V est un voisinage de Xloc.
RESOLUTION NUMERIOUE
Nous
nous intéresssons ici à l'aspect numérique de la résolution de
(Ploc).
Pour ce faire. on considère qu'on
construit itérativement une suite
XI. X2,...• Xn•... de façon à trouver un point Xloc solution de (Ploc).
15
algorithme
général
Initialisation de XQ, tk' do
Tant
que Test d'arrêt faux faire
choisir tk' dk
Xk+l= x, + tk·dk
Fin (de tant que)
~ et
dk
sont appelés respectivement pas et direction d'avancée. Pour
les calculer, l'algorithme peut utiliser des
informations sur les
valeurs
du gradient (voire du Hessien) de f au point X k courant, si toutefois ces
informations sont disponibles. C'est ainsi qu'on pourra distinguer deux
classes de méthodes:
- les méthodes du type eradjents
Elles
utilisent
les
gradients
et/ou
les
Hessiens
de
f
(ou
leurs
approximations) pour choisir des
tk et dk
appropriés, de façon que la
direction choisie soit une direction de descente.
- les méthodes
directes
Elles ne sont pas basées sur l'utilisation des gradients, mais sur les
seules évaluations de fonctions.
Les
méthodes directes sont utiles quand f n'est pas différentiable
ou quand des approximations des gradients sont difficiles à calculer.
Parmi
les
méthodes
directes,
on
trouve
beaucoup
de
méthodes
à
recherche aléatoire; on trouve aussi des méthodes dont les algorithmes
n'obéissent pas vraiment au schéma de l'algorithme général, c'est-à-dire
dont les directions d'avancée ne sont pas choisies de la façon
définie ci-
dessus; en effet, il existe même des méthodes n'utilisant pas de pas et
direction d'avancée dans le sens classique. Exemples de la recherche de
Fibonacci
où
l'algorithme
(en
dimension
1)
s'attache
à
réduire
l'intervalle
contenant
le
minimum
global
-le
seul
d' ail1eurs
par
hypothèse-, le "polytope mouvant" de Nelder et Mead où, partant de N+1
points
(pour
le
cas
N-
dimensionnel),
on
effectue
des
opérations
géométriques sur ces points et on remplace au fur et à mesure les points
inutiles
(points
de
très
grande
valeur
de
fonction),
jusqu'à
la
convergence.
16
1.2.2.
Quelques
exemples
d'algorithmes
efficaces
La programmation linéaire est une branche de l'optimisation qui a
été
suffisamment
développée.
Des
codes
efficaces
existent
pour
sa
résolution. Donc, dans ce qui suit, nous ne présenterons que des codes
pour l'optimisation non linéaire.
1.2.2.1 Algorithme général des méthodes de descente
On dira que la direction dk est une direction de descente si f( Xk+ tk.dk)
est inférieure à f(Xk) pour un pas tk donné.
Les méthodes de descente sont celles où les tk et dk sont supposés être
ceux qui donnent une direction de descente et un pas optimal. Ce qui
donne:
Initialisation de Xo• 10, do (point initial, pas. direction)
Tant
que Test d'arrêt faux faire
calcul de la direction (de descente)
d k
calcul du pas (optimal) t k •
X k + 1 = X k + tk·dk
Fin (de tant que)
(Xk , f(Xk»
est solution.
Remarque:
Ce qui différenciera les algorithmes basés sur de telles méthodes
sera
essentiellement:
- le choix de la
direction: on peut utiliser les gradients. le Hessien ....
suivant le type de problème à résoudre.
- la façon de calculer .!!Lw: on utilise souvent une autre minimisation
uni-dimensionnelle le long de la direction de descente; ce qui est alors
appelé recherche
linéaire.
- les
critères d'acceptation
des
pas
et
directions
pour éviter
non
seulement les divergences, mais aussi les lenteurs de convergence. ou
pour tenir compte d' éventuelles contraintes.
- et, bien sûr, le test d'arrêt de l'algorithme.
17
1.2,2.2 Exemple: al~orithme D.F.P (Davjdon-Flecher-PowelI)
C'est un exemple d'algorithme de Quasi-Newton (cf. 1.2.2,3.).
On note
gk = g (Xie)
le vecteur gradient de f en XIe
Hk = H (Xie)
une matrice NxN définie positive.
Njyeau
0:
Initialisation
E.
Il ,XO. HO;
Niveau
J:
poser Ilk = - Ok • gk
Trouver À.k
qui minimise ( en fonction de À.)
f(Xk+À.llk);
poser alors
Vk= U dk
Xk+1 = Xk + Vk;
calculer
gk+ 1
faire Uk = gk+1 • gk
réactualiser Hk-s1 (en utilisant Hk, Vk et Uk voir plus loin);
Niveau
2 :(Test
d' arrêt)
Si
1 V k 1 < E
ou 1 g k + 1 1< Il
STOP,
Sinon
faire k = k +1 et
aller il niveau
1.
calcul de
Hk :
Hk:+1 = Hk + Ak + Bk
avec
Ak = Vk VTk 1 < Vk, Vic>
Bk = - Hk Vie VTk Hk
1 < Vie, Hk: Uk>
<.,.> désigne le produit scalaire.
convergence
du
D,F,P.
En théorie, le D,F,P. converge en n itérations quand la fonction à
minimiser est quadratique sur Rn.
C'est un code assez satisfaisant pour des problèmes différentiables sans
contraintes, quadratiques ou non.
exemple de résultats
pour la fonction test
FI3 (cf. 3.1)
f(X) = (XI+IO X2) 2 + 5 (X3 - X4) 2 + (X2 - 2 X3) 4 + 10 (XI- X4) 4
avec Xo = (3,-1,0,1) on trouve:
f*=3.04
10. 13
X* = (1.85
10.4; -1.85 10 .5; 3,51 10.4; 3.15 10 .4)
18
1.2.2.3. Méthodes de Quasi-Newton
Nous rappelons ces méthodes car ce sont celles qui sont le plus
utilisées en optimisation locale différentiable.
Qn note
gk =g (Xk)
le vecteur gradient de f en Xk
Hk =H (Xk)
une matrice NxN définie positive.
• aleorithme
eénéral
Niyeau
0
Initialisation E, li • XO, HO;
Njyeau
1
si 1 g k 1< E STOP;
sinon poser ôk = - Hk • gk;
Niyeau
2
trouver
(en panant de Â.=I) Â.k qui minimise la fonction de Â.:
f ( Xk + Â. dk );
Niyeau
3
poser alors
Vk= Â.k dk
Xk+1 =Xk+Vk
calculer
gk+ 1;
Niyeau
4
faire Uk = gk+ 1 - gk
réactualiser Hk+ 1 (en utilisant Hk, Vk et Uk voir ci-dessous);
calcul de Hk(dans le cas général):
Hk-s I = Hk + Bk
de façon que:
(1) l'équation de Quasi-Newton, Hk+1 Uk = Vk, soit satisfaite pour tout k
(2) Bk soit symétrique définie positive pour tout k;
(3) que Bk soit minimale en un certain sens.
Ce qui laisse beaucoup de choix posibles.
RemarQue
Actuellement,
la
classe d'algorithmes de
Quasi-Newton
la
plus
populaire (et efficace) est celle des algorithmes de B.F.G.S du nom de ses
auteurs (Broyden, Fletcher, Goldfard et Shanno). Ils se caractérisent par
le calcul de Hk+ 1 suivant :
Hk+1 =Bk - [ Vk UTk Hk + Bk Uk VTk] / <V , ui> +
[ 1+ <Uk, BkUk> /<V .ue-j Vk VTk /<V .ue>
avec
Vk = Xk+1 - Xk
Uk = gk+1 - gk,
et <.,.> le produit scalaire usuel.
19
1.2.3.
Transition
vers
l'optimisation
globale
Toutes
les
bibliothèques
scientifiques
proposent
des
codes
numériques d'optimisation locale assez satisfaisants. Pour les problèmes
sans
contraintes,
les
codes
de
Quasi-Newton
donnent
de
très
bons
résultats.
D'autre part, même pour les cas où la fonction à minimiser n'est
pas
"assez"
différentiable,
il
est
possible
de
trouver
des
codes
numériques de Recherche Directe qui donnent satisfaction : le polytope
mouvant par exemple (cf. 3.3).
Il en est de même pour certains cas avec contraintes, où on trouve
satisfaction en utilisant quelques techniques de pénalités appropriées .
Enfin nous recommandons les codes de Quasi-Newton car beaucoup
de
versions ont été renforcées pour tenir compte des
problèmes
de
divergence (en adaptant des pas appropriés).
Comme
codes
efficaces de Quasi- Newton citons VA 13A de la
bibliothèque HARWELL et M2QNI du club MODULOPT.
Ces codes sont tous deux du type B.F.G.S.(cf. 1.2.2.3).
Conclusion
En optimisation locale, il existe des codes efficaces pour résoudre
pratiquement
tous
types
de
problèmes
(bien
posés)
: des
codes
de
recherche directe, des codes de recherche indirecte, avec ou sans calcul
de gradients, etc ...
Cependant, il arrive qu'un utilisateur, en exécutant plusieurs fois
de suite avec des points initiaux différents, tel ou tel code (programmé
par
lui-même
ou
pris
dans
une
bibliothèque
scientifique),
se
rende
compte qu'il trouve plusieurs minimums dont les
valeurs diffèrent de
façon significative. Il devient alors, nécessaire de faire appel à d'autres
outils pour aider à déterminer quel minimum est le plus petit. Et, de
déterminer un de ses arguments (s'il en a plusieurs) afin de trouver un
choix
de
paramètres
qui
donnent
un
critère
(valeur
de
fonction)
objectif optimal
global
sur
tous
les
choix
possibles.
C'est ce
genre de problèmes que
l'optimisation globale se propose de
résoudre.
20
1.2.4.
Références
bibliographiques
[BUN84j B. D. BUNDA Y, Basic Optimization Methods, Edward Arnold
(Publishers) Ltd. (1984).
[CEA71] J. CEA, Optimisation: théorie et algorithmes, Dunod, Paris (1971)
[DEN83]
J.E.
DENNIS and
R.B.
SCHNABEL, Numerical Methods for
Unconstrained Optimization and Nonlinear Equations, Prentice-Hall (1983).
[GIL74j P.E. GILL, & W. MURRAY, Numerical methods for Constrained
Optimization, Academie Press London & New York (1974).
[GIL81] P.E. GILL, W. MURRAY & M.H. WRIGHT, Practical Optimisation,
Academie Press London & New York(198 1).
[HOC81] W. HüCK et K. SCHlTTKüWSKl, Test Examples For Non Linear
Programming
Codes,
Lecture Notes
in Economies and
Mathematical
Systems, 187, Springer-Verlag, (1981).
[JET86] M. W. JETER, Mathematical Programming: an
introduction to
optimization, Marcel Dekker, Inc., New York and Basel (1986).
[LEM89] C. LEMARECHAL, Méthodes numériques d'optimisation, Notes de
cours D.E.A. d'automatique Université Paris IX- Dauphine (1989).
[SCH87]
K.
SCHlTTKüWSKI,
More Test Examples For Non
Linear
Programming
Codes, Lecture Notes
in Economies and
Mathernatical
Systems, 282, Springer-Verlag, (1987).
21
Chapitre 2 :
METHODES
NUMERIQUES
D'OPTIMISATION
GLOBALE
22
2.0.
Rappel
du
problème
Dans
ce
chapitre
nous
allons
présenter quelques
grands
traits
caractéristiques
de
l'optimisation
globale
numérique.
Rappelons
que
nous n'insisterons que sur l'aspect méthodes et algorithmes numériQues.
Pour ce faire nous allons d'abord nous attacher à faire une classification
des
types
de
problèmes
considérés
puis
des
types
de
méthodes
élaborées pour la résolution de ces problèmes.
Rappelons le problème général
(P): Minimiser f(x) où f: R n ~ R continue,
ce qui est équivalent à résoudre numériQuement (P) sur un
compact S
de Rn (cf. 1.1.2.), soit donc:
n
.
Trouver (Xopr , fopl) de R
X R. ou Iopt ef'(Xopt) ,
(P)
lleI que:
V X e S . f(Xopl) S f(X)
Solutions
numériques
de (P)
Théoriquement, résoudre (P)
reviendrait à déterminer exactement
le
plus petit de tous les minimums locaux (ou un de ses arguments X·
au cas où il y aurait plusieurs minimiseurs globaux). Cependant, d'un
point de vue numérique, une solution ne pourrait être obtenue qu'avec
une
certaine précision donnée. Par
conséquent, le problème (P) sera
considéré résolu dès qu'on aura atteint un élément, pour un seuil E
donné, d'un au moins des ensembles Ax.(E), Arx.(E). AOX*(E) suivants
(où X· est un minimiseur global) :
• Ax.(E) = {x E S tel que /Ix. x·1I < E}
• Arx.( E) = {x E S tel que If(x) • f(x·)1 < E}
• ADx.(E) = {x E S tel que D(f(x» , E}
. 0 ()
m ( {z E S tel que f(z) , y}
ou
y
-
m(S)
avec m mesure de Lebesgue.
23
2.1.
Classification
des
problèmes.
2.1.1 Problèmes avec contraintes, problèmes sans contraintes
Une classification naturelle des problèmes serait de
distinguer les
problèmes sans contraintes (Psc) et ceux avec contraintes (Pc).
Numériquement, il est plus aisé de résoudre les (Psc) et cela, quelle que
soit la nature de la méthode utilisée (stochastique ou déterministe). En
effet, le domaine (des variables) de recherche n'étant pas restreint, on a
moins de soucis à se faire que dans les cas (Pc) où il faut tout mettre en
œuvre pour s'assurer à chaque itération que l'on reste dans le domaine
du
"faisable" ou
sinon
s'imposer un moyen d'y
revenir
à la
limite
(méthodes de pénalités par exemple), et cela n'a d'autre effet que de
compliquer davantage les calculs! Dans tous les cas - et pour le moment
- la
majorité
des
méthodes numériques
est encore
conçue
pour les
problèmes sans contraintes. Cependant il existe des contraintes qu'on
peut qualifier de simples et qui sont très faciles à adapter dans un code
conçu initialement pour les (Psc); c'est par exemple le cas où l'ensemble
C des contraintes s'exprime comme un polyèdre.
Enfin, des études sont menées pour presque toutes les méthodes dans le
but de trouver des adaptations aux cas avec contraintes.
2.1.2 Problèmes
particuliers, problèmes
2énéraux
Par problèmes "particuliers" on entend des
problèmes ayant des
propriétés particulières qu'on utilise pour élaborer une
méthode propre
à ce
type de problèmes. Bien
entendu une telle méthode ne pourrait
traiter que de tels problèmes.
Exemples:
- les problèmes
D.C. où la fonction à rmrnrmser est sous la
forme f = g - h
avec g et h
deux fonctions convexes. Des méthodes
appropriées sont développées par [HORST et TUY 1991].
- les problèmes
à
données
fractions
(optimisation
fractionnaire) où
f = ~ avec h>O.
Une méthode est proposée par [PARDALOS et ROSENI991].
- les problèmes
Lipschitziens
où
f est
supposée
être
partout
L-Lipschitz.
Plusieurs
versions
de
la
méthode
de
Pijavskij-
Shubert sont proposées notamment dans [PI172] , [SHU72], [MLA86] etc ...
24
RemarQues
Il faut noter que la différentiabilité est très souvent exigée pour la
résolution de tels problèmes.
Il est aussi clair que ces problèmes ne représentent qu'une infime partie
des problèmes d'optimisation globale.
Tous les
autres problèmes sont classés "généraux". Parmi eux,
certains ne sont même pas différentiables (Le. la fonction f n'est pas
différentiable) voire même pas partout continus!
Evidemment, il n'y a pas de méthodes qui prétendent traiter tous
les problèmes généraux. Ainsi souvent la continuité est exigée ou bien
d'autres hypothèses (même faibles). Mais dans tous les cas la méthode
de
résolution
n'est
pas
basée
principalement
sur
des
hypothèses
particulières
de
f,
contrairement
aux
cas
des
méthodes
pour
la
résol ution de problèmes particuliers.
2.1.3. Conclusion
En résumé, on peut dresser le tableau de classification suivant:
• Problèmes sans contraintes:
• problèmes résolubles - problèmes de formes spéciales.
• problèmes généraux.
• Problèmes avec contraintes:
• problèmes de formes spéciales.
• problèmes généraux.
2.2.
Classification
des
méthodes.
En dépit de l'abondance des méthodes proposées, on peut leur
trouver des traits caractéristiques suivant l'approche, la philosophie, ou
la urantje ~ convereence.
D'où la possibilité de trouver plusieurs
classifications
suivant
le
critère
retenu.
Nous
donnons
quelques
exemples de classifications.
Pour chaque
nom
de
méthode,
nous
donnons
une
proposition
de
traduction en Français et le nom (courant) en Anglais (cf. glossaire).
25
2.2.1. Selon l'approche
On distingue deux types d'approches selon qu'elles incluent ou non
des processus probabilistes [DIXON and SZEGO 1978] et [GOMULKA
1978]:
• méthodes déterministes:
• Méthodes à Recherche par Quadrillage
(Grid Search Methods);
• Méthodes des Trajectoires
(Trajectory Methods);
• méthodes stochastiques:
• Méthodes à Recherche (linéaire) (à stratégie) Aléatoire
(Random (Iine) Search Methods);
• Méthodes de Regroupement en grappes
(Clustering Methods);
• Méthodes par Echantillonnage (Sampling Methods);
Les
méthodes stochastiques sont des
méthodes qui
utilisent des
stratégies probabilistes pour assurer la convergence vers la solution. En
général, celle convergence est assurée par une relation du type:
Soit (Xk)k
la suite des itérés, on est alors assuré que:
P ( Xk -4 X·) > 1 - E
où X· est une solution et E >0 donné.
La
convergence des
méthodes
déterministes,
contrairement aux
méthodes
stochastiques,
serait assurée de
façon
implicite
suivant
le
principe même de
la méthode. C'est-à-dire que
leur
convergence est
théoriquement
garantie
mais
seulement
grâce
à
des
hypothèses
supplémentaires sur la fonction f. C'est le cas (théoriquement du moins)
de la méthode du Tunnelier (cf. 3.2).
26
2.2.2. Selon la philosophie
Une autre façon de classer consiste
à distinguer
les
méthodes
suivant leur philosophie même [RINNOOY KAN et T1MMER 1986] :
• Partition et Recherche (Partition & Search)
• Approximation de f et Recherche (Approximation of f & Search)
• Décroissance Globale ou Amélioration permanente de f
(Global decrease or permanent improvement in f)
• Amélioration du minimum local Le. construction d'une suite de
minimiseurs à valeurs de fonctions décroissantes
tlmprovement of local minima)
• Énumération des minimums locaux
(Enumeration of local minima)
Partition & Recherche
S est
subdivisé en
sous-régions de
plus
en
plus
petites
dans
lesquelles le minimum global serail localisé;
Exemples
des
méthodes
par
Procédure
de
Séparation
et
Evaluation
Progressive (P.S.E.P) (Branch-and-Bound) [STRONGIN 1978], [PINTER
1983], [HORST 1986] ...
Aaproximation et Recherche
Dans cette philosophie, f est remplacée par une approximation de
plus en plus précise, qui est plus facile à traiter numériquement.
Exemple [SHUBERT 1972] : f est supposée L-Lipschitzienne. donc on peut
la remplacer par Gk où :
• étape k=1
GI (X) = f(X1) - LI/X - X 1 Il
• étape k :
Gk (X) = MAX [ f(X) - L 1/ X-Xk Il
Gk-I (X)]
on est alors assuré que:
Ok-I (X) :5:Gk(X):5:f(X)
pour tout X de S,
Gk (Xi) = f (Xi)
i= 1.2.....k
Xi minimiseur local trouvé à l'étape i-I
A chaque itération. l'approximation se précise davantage; et l'algorithme
s'arrête quand un élément de Af(e) est trouvé.
Méthodes de Décroissance Globale :
Ce sont des méthodes utilisant des évaluations permanentes de f et
engendrant
(par un
modèle
propre)
une
suite
de
valeurs
de
f qui
converge vers le minimum global.
27
Mlthodes d'Amélioration des minImums locaux:
Elles utilisent
des sous-méthodes "locales"(Le. des sous-méthodes
pour déterminer des minimums locaux). Leur principe est de répéter
une stratégie (propre à chaque méthode), tout en se servant des sous-
méthodes locales (qui elles, déterminent des minimums locaux),
pour
construire
une
suite
de
minimiseurs
locaux
à
valeurs
de
fonction
décroissantes.
Le
minimiseur
global
serait
théoriquement
le
dernier
local trouvé. Exemple [LEVY et MONTALVO 1985]
Méthodes d'Enumiration des minimums locaux:
Il s'agit ici, d'énumérer tous les minimums locaux (en utilisant des
sous-méthodes d'optimisation locale); le plus petit est solution. Encore
faut-il avoir les moyens (et la certitude) de les énumérer tous!
2.2.3. Selon la 2arantie de conver2ence
Torn
[TORN
1974],
influencé
par
[LEON
1966],
proposa
la
classification suivante:
• Méthodes avec garantie de convergence
(Methods with guaranteed convergence)
• Méthodes à Recherche Aveugle (Blind Search)
• Méthodes de Combinaison des Recherches Aveugle et Locale
tCombination of Blind & Local Search)
• Méthodes de Recherche non Locale (Non Local Search)
Méthodes avec garantie de
conyergence:
Elles garantissent qu'une certaine précision (prédéterminée) serait
atteinte au bout d'un nombre fini de pas. On pourrait même les appeler
aussi "Méthodes à Recherche Exhaustive".
Mlthodes à Recherche Ayeugle :
On
suppose qu'on dispose d'une
stratégie
nous
permettant de
sélectionner des
points d'intérêt : points vérifiant par
exemple
des
conditions nécessaires et/ou suffisantes d'optimalité globale, points qui
sont dans les voisinages des minimums candidats à être globaux .... ou
toute autre stratégie possible.
Dans ces cas, on peut procéder à une recherche jusqu'à "localiser" de tels
points, autrement dit, on peut faire une recherche aveugle jusqu'à les
trouver.
28
Mlthodes de Combinaison des
Recherches Aye!l~le et Locale :
Il
semble
naturel
de
faire
une
combinaison
de
la
recherche
aveugle avec
une recherche locale (voir ci-dessous pour la recherche
locale).
Méthodes de
Recherche Locale :
Les itérations se font de façon "continue" autour des derniers points.
C'est-à-dire que chaque état est pris dans
un voisinage immédiat de
l'état précédent.
Mlthodes de
Recherche non Locale :
On entend par là, toute possibilité de "sauter" un mmrmum local
et
pouvoir atteindre le global (ou du moins un meilleur point) même si
celui-ci n'est pas dans un voisinage immédiat de l'état précédent.
2.2.4. Autres classifications
Plusieurs auteurs ont proposé d'autres classifications:
[LEON 1966]
• Méthodes à Recherche Aveugle
• Méthodes de Recherche Locale
• Méthodes de Recherche non Locale.
[ARCHEITI et SCHOEN 1984]
Méthodes Déterministes:
• Méthodes par Recouvrement
(Covering Methods);
• Méthodes des Trajectoires, du Tunnelier;
Méthodes Probabilistes:
• Méthodes basées sur l'échantillonnage aléatoire
(Methods based on random sampling (Clustering methods);
• Méthodes à Recherche Aléatoire;
• Méthodes basées sur un modèle stochastique de la fonction
objectif
(Methods based on a stochastic model of the objective
function).
[ZILINSKAS 1986]
• Méthodes de transiuon d'un rmmmum local à un autre
(Sequentiel transition from one local mminimum to another);
• Méthodes à Recherche Aléatoire (Random Search methods);
• Méthodes par Recouvrement
(Covering Methods);
• Méthodes basées sur un modèle stochastique de la fonction
objectif
(Methods based on a stochastic model of the objective function).
29
2.3.
Choix
d'une
classification
On peut se demander à quoi sert une classification des types de
problèmes
et
de
méthodes.
A
notre
avis
quand
un
utilisateur
est
confronté à un
problème d'optimisation globale, la première chose à
faire est de bien cerner le problème; à savoir:
• hypothèses sur f (différentiabilité ... ),
• hypothèses sur le domaine de recherche,
• existence ou non de contraintes, quels types de contraintes,
• coût d'évaluation de la fonction (temps CPU, nombre de sous-
programmes nécessaire),
• facilité de l'évaluation (accès, formule explicite de f)
• précision dont on dispose sur les calculs,
• type de matériel utilisé,
• temps dont on dispose pour résoudre le problème,
• le problème est-il en temps réel ou pas ?
On voit bien, donc, que connaître une classification des problèmes
et des méthodes, peut faciliter la tâche de l'utilisateur et
le guider dans
son choix; ce qui lui permet de fixer ses objectifs en conséquence. CaLi!
est
préférable d'avoir
une
solution
approchée
(avec
une
précision
insuffisante) dans un délai raisonnable gu 'une solution exacte (avec la
précision souhaitée) hors délais dans certains cas contractuels.
De même,
pour un chercheur, une classification
permet de savoir à qui
va s'adresser son code et le spectre d'utilisation qu'il en aurait, et donc
cela influence l'élaboration d'un algorithme quelconque.
Pour notre part, nous avons choisi deux classifications :
•
la
première étant celle distinguant
les
méthodes
déterministes
des
stochastiques, soit:
Première classification:
Méthodes (à stratégies) déterministes;
Méthodes (à stratégies) stochastiques;
30
• la deuxième classification est celle de [TORN & ZILINSKAS 89] qui
résume à peu près l'état d'avancement actuel des méthodes.
Deuxième classification;
• Méthodes avec garantie de convergence
(Methods with
guaranteed accuracy):
• Méthodes par Recouvrement (Covering methods);
• Méthodes Directes (Direct Methods);
• Méthodes à Recherche Aléatoire
(Random Search methods);
• Méthodes de Regroupement en grappes
(Clustering methods);
• Méthodes de Descente Généralisée (Generalized Descent);
• Méthodes Indirectes ;
• Méthodes
d'approximation des surfaces de niveaux
(Methods approximating the level sets);
• Méthodes d'approximation de la
fonction objectif
(Methods approximating the objective function).
La
première classification
est
la
plus
couramment
utilisée
dans
la
littérature. Elle couvre aussi toutes les sous-classes de la deuxième.
Néanmoins, nous présentons avec un peu plus de détails la deuxième
classification, pour mieux faire ressortir la richesse et la profondeur des
philosophies de ses sous-classes.
2,3,1. Méthodes ayec
garantie de conyergence :
Il
est
naturel
de
distinguer
les
méthodes
suivant
qu'elles
garantissent ou non une
précision sur la solution numérique. En effet,
obtenir une solution avec une certaine garantie de précision implique
une recherche exhaustive sur S pour atteindre un minimum global. Or
cela n'est possible que pour les problèmes pour lesquels on a certaines
informations a priori. Une recherche exhaustive exige un effort de calcul
considérable, comparée aux méthodes ne garantissant pas qu'une telle
solution est atteinte.
2,3.1. t. Méthodes par Recouvrement
Les cas les plus simples de ces méthodes sont celles qui consistent
en la détermination de sous-régions ne pouvant contenir te minimum
global, et leur élimination pour la suite. C'est comme une procédure de
localisation qui se termine par une sous-région de taille suffisamment
31
réduite pour donner la solution avec
la précision souhaitée. On
peut
aussi les voir comme des procédures donnant au fur et à mesure des
approximations par défaut et/ou par excès des
valeurs (de
fonctions)
candidates à être la valeur minimale (globale).
Exemple:
Algorithme
de
Pijavski-Shubert (en dimension 1)
Soit
à minimiser une fonction f L-Lipschiz sur [a,b],
fi=f(Xi)
i=I, ... .k
a
~ Xi ~ b , l'échantillon de valeurs,
fOk = min ( f(Xi) i=I,... .k},
et
Ak ={ x tel que
Max (fi - L 1X - Xii) ~ fOk
alors on sait que
f(X) ~ fOk
pour tout X de Ak;
il s'en suit que l'ensemble
Ak
peut être
éliminé
pour la suite des
itérations car il ne pourrait contenir un minimiseur global.
Sans rentrer dans les détails théoriques, on imagine bien comment en
considérant les sous-régions qui restent et, en réappliquant le processus,
on éliminera au fur et à mesure toutes ceIles qui ne contiennent pas le
minimiseur global.
2.3.2. Mtthodes
Directes:
Elles n'utilisent (seulement) que des informations locales (valeurs
de fonction) pour construire un modèle convergeant vers une
solution.
2,3.2,1, Méthodes de Recherche Aléatoire
Ces méthodes ont beaucoup de succès car simples à programmer.
Leur principe consiste à utiliser des processus probabilistes pour décider
en un point comment "avancer" et un critère pour valider les résultats
au
cours des
itérations.
Notons que mathématiquement,
il
n'est pas
encore
bien
prouvé
dans
quelles
situations
exactement
on
peut
se
permettre -ou pas- d'utiliser cette notion de "hasard",
32
Exemples d'aleorithmes de Recherche Aléatojre
Dans
ce
qui
suit.
nous
allons
présenter quelques
méthodes
à
Recherche Aléatoire. La plupart de ces
méthodes consistent en deux
étapes essentielles:
étape 1:
échantillonnage (ou phase globale)
Au cours de cette étape, le domaine S est discrétisé (suivant une façon
donnée)
en
un certain
nombre
de
points.
En
ces
points
(dits
de
J'échantillon), on évalue les valeurs de la fonction à minimiser.
étape 2:
raffinement (ou phase locale)
Durant cette deuxième étape. des
processus de
raffinement
peuvent
intervenir pour sélectionner "les meilleurs" points (ou "meilleurs" sous-
ensembles) selon
certains critères propres à la méthode choisie.
A partir de
ces
points
sélectionnés
sont
appliquées
les
procédures
propres
à
la
méthode;
ce
qui
conduit
théoriquement
à cerner un
candidat au minimum global.
Méthode
de
Recherche
Aléatoire
Pure
:
C'est la plus simple des méthodes de Recherche Aléatoire.
Elle consiste seulement en une phase: la phase globale.
• algorithme de la Recherche Aléatoire Pure
Niyeau
J: Evaluer f en N points obtenus à I'aide d'une distribution
uniforme sur S;
Niyeau
2: La plus petite valeur de f trouvée est candidate pour être
le
minimum global.
• convergence de la Recherche Aléatoire Pure (ou uniforme)
Théorème: [TIMMER 1984)
Soit (Xk)k
k=I •...•N l'échantillon choisi.
si f est continu on a le résultat suivant:
f(X k) -+ f(X*)
si et seulement si l'ensemble { X k
k=I .....N } est dense dans S.
33
Commentaires
On voit bien que la convergence de la Recherche
Alëatoire
Pure
exige des hypothèses très fortes.
Ainsi, même pour assurer la convergence avec une certaine probabilité
p donnée, il faut choisir N tel que:
p := 1 - (1 - pl ) N
où pl est la probabilité, pour la discrétisation choisie, d'avoir un point
de l'échantillon dans la région d'attraction d'un minimum global.
Si on choisit p, on en tire la valeur de N nécessaire.
Méthode
des
Initialisations
Multiples
A cause de l'extrême simplicité de la Recherche
Aléatoire
Pure
mais aussi de ses faibles performances numériques, plusieurs extensions
ont été proposées, partant elles aussi d'un échantillon uniforme dans S,
mais en
rajoutant une deuxième phase de raffinement (phase locale).
Une manière simple de faire
cette deuxième phase est donnée par la
méthode des Initialisations Multiples .
• algorithme général de la méthode des Initialisations Multiples:
Niveal!
Q : Evaluer f en N points obtenus par une distribution
uniforme sur S
Niveau
1: Choisir un point de l'échantillon.
Niveau 2
Appliquer une phase locale en ce point.
Niyeau 3
Un test d'arrêt indique s'il faut s'arrêter ou
revenir au niveau
I.
A l'arrêt, le minimiseur local avec la plus petite valeur de fonction est
candidat pour être le minimiseur global...
Pour les choix possibles pour les différents niveaux voir ceux proposés
en 2.3.2.2.
34
2.3.2.2. Méthodes de Regroupement en grappes
Les
méthodes
de
Regroupement
en
grappes
sont
des
méthodes
d'optimisation
globale
utilisant pour une
grande
part
des
techniques
dites de
regroupement.
La justification est la suivante:
a) il est possible d'avoir un échantillon A (de S) de points de
concentration dans les voisinages des minimiseurs locaux de f.
b) on
peut appliquer une
certaine technique de
regroupement
aux
points de
cet échantillon de façon
à identifier les
voisinages de
minimiseurs locaux, et donc de permettre l'application des procédures
de minimisation locale.
c) l'enchaînement a) - b) peut être implémenté efficacement de
façon
à
être
complété
par
d'autres
méthodes
proposées
pour
l'optimisation globale.
Citons
comme
exemples
d'algorithmes
basés
sur
ces
méthodes
(cf.
[TOR89]):
l'algorithme de Becker et Logo,
l'algorithme de Torn,
l'algorithme de Spircu,
l'algorithme de Priee,
l'algorithme de Timmer.
exemple 1
l'algorithme de Torn
Niveau
1:
(échantillonnage)
Engendrer aléatoirement N points de S dans Rn;
Niveau
2:
(concentration)
Concentrer ces
points
autour des
minimums;
Niveau
3:
(regroupements
en
grappes)
Appliquer
une
technique
de
regroupement;
Niveau
4: (test d'arrêt)
Si une condition d'arrêt (par exemple s'il ne reste qu'un
seul
point pour chaque "grappe" ... ) aller au pas 6;
Niveau
5
(réduction
de
l'échantillon)
Retenir les meilleurs points (de chaque grappe),
Aller au Pas 2;
Niveau
6: (trouver les minimums locaux)
Déterminer les minimums locaux en partant du meilleur
point
de chaque grappe;
35
les
niveaux:
Nous allons détailler certains niveaux de cet algorithme
car
ils
sont
représentatifs
de
la
philosophie
des
méthodes
de
regroupement.
échantiJIonna~e :
Pour cette
phase
de
l'algorithme, commune
à
beaucoup de
méthodes (stochastiques surtout), il y a plusieurs possibilités dont les
deux suivantes :
1)
Engendrer
aléatoirement
(suivant
une
loi
donnée)
N
points.
2)
Utiliser une technique de subdivisions:
pour ce faire, S est divisé (suivant ses composantes par exemple) en
sous-régions. Dans
chaque sous-région, on prend des
"représentants".
Par exemple, on divise S en N sous-régions et on prend pour chacune le
"centre de gravité", ou bien on peut diviser S en N/m sous-régions et
pour chacune prendre m points.
concentration:
Nous citerons deux exemples :
1)
On retient un nombre prédéterminé de
points
ayant
les
plus petites valeurs de fonction.
Cette technique est peu utilisée car elle peut s'inclure dans la phase de
réduction.
2)
Pour
chaque
point
de
l'échantillon
initial,
appliquer
quelques pas d'un algorithme d'optimisation locale. Ceci aura pour effet,
de
faire
"descendre"
les
valeurs
de
fonction.
Les
nouveaux
points
obtenus seront, pour certains, plus
rapprochés (ou concentrés) autour
d'un
minimiseur
local
(celui
qu'on
aurait
trouvé
si
on
exécutait
l'algorithme
complètement). En
d'autres
termes,
tous
les
points
qui
étaient potentiellement des points initiaux d'un même minirniseur local,
après cette opération de quelques pas, conduiront à des points un peu
plus "rapprochés" de ce minimiseur. Ce qui, au total, revient à concentrer
des points autour des minimiseurs locaux.
Bien sûr, après cette phase on élimine un certain nombre (prédéterminé
ou à définir) de points.
technique de re2roupement
Le but est de regrouper tous les points susceptibles de donner
une même grappe. Cette grappe serait une région d'attraction d'un seul
minimiseur. Pour cela, on peut envisager diverses possibilités dont:
1) définir une relation d'appartenance à une même grappe, par
exemple
quand
deux
points
donnés
sont
séparés
d'une
distance
inférieure à un certain seuil. Il suffirait de partir de quelques points
générateurs (n'importe lesquels). On rajoute, au fur et à mesure, tous les
36
points proches d'un des points de la grappe qui grossit ainsi jusqu'à se
former
entièrement.
Evidemment,
la
technique
de
regroupement
s'achève quand tous les points ont retrouvé leur grappe!
2) ici on se fixe d'abord des centres pour chaque grappe. Ce qui
peut se faire de différentes manières, par exemple utiliser les N points
de l'échantillon pour faire l'approximation d'une fonction densité (ce qui
exige que
N soit
très grand) et utiliser les points d'intérêt de cette
dernière fonction comme centres des grappes.
Une fois que sont fixés
le nombre des grappes et
leurs centres, on
rajoute chaque point de l'échantillon à la grappe dont le centre est le
plus proche.
Cette technique est peu utilisée.
test d'arrêt:
En optimisation globale
numérique en
général, et
pour les
méthodes
stochastiques
(dont
les
méthodes
de
regroupement)
en
particulier, le test d'arrêt est une des phases les plus délicates car l'une
des plus "personnalisées" lors de l'implémentation de la version choisie
d'un algorithme (cf. 2.5.1).
Pour les méthodes de regroupement, en plus des tests classiques, on
peut retenir les suivants:
1 )
Si tous les points ne sont plus distants de plus d'un seuil
donné, alors ils appartiennent à une seule grappe qui est elle-même
réduite
au
maximum. On
peut
donc
dire
que
c'est la
grappe d'un
minimiseur
global
que
l'on
calcule
en
utilisant
(complètement)
un
algorithme d'optimisation locale.
2)
Si après d'éventuelles réductions, il ne reste qu'un point
par grappe, cela signifie que chacune de ces grappes ne contient qu'un
seul minimiseur. On peut, à ce moment-là, appliquer à chacun de ces
points (complètement) un algorithme d'optimisation locale et choisir le
plus petit minimum comme global.
réduction de l'échantillon:
A chaque boucle de l'algorithme, on réduit le dernier échantillon
réactualisé, qui, rappelons-le, est gardé en mémoire avec ses valeurs de
fonction, de façon que l'échantillon final vérifie une condition d'arrêt.
Pour cette réduction, plusieurs choix ont été proposés dont:
1)
le plus courant, qui consiste à retenir un certain nombre
de points de plus petites valeurs;
2)
garder un point sur toute paire de points qui sont séparés
par
une
distance
inférieure
au
seuil
du
regroupement;
On
pourrait
37
considérer en effet, que de tels points sont identiques numériquement
(ou du moins ont des propriétés semblables).
2.3.2.3. Méthodes de Descente Généralisée
Des d'efforts ayant été faits en optimisation locale pour élaborer des
codes satisfaisants, il est tout naturel d'envisager leur utilisation pour le
cas de l'optimisation globale.
Pour ce faire, il y a deux possibilités:
-
modifier
l'équation
différentielle
décrivant
la
trajectoire
de
descente locale;
- appliquer
plusieurs fois
un
algorithme
local
à
une
fonction
auxiliaire associée à f.
Citons comme exemples du deuxième cas,
des
méthodes où
un code
d'optimisation
locale
est
appliqué
successivement
à
des
fonctions
auxiliaires de f qui
tiennent compte de pénalités relatives à chaque
minimum local trouvé, et permettent
d'exclure ces derniers. La fonction
de
pénalités
permet d'exclure
les
minimums
locaux
déjà
rencontrés
avant d'appliquer une nouvelle minimisation locale.
Exemple de fonction de pénalités
• dans la Méthode
du
Tunnelier (Tunneling
cf. 3.2)
Soit X* un minimiseur local ,
pour trouver un X tel que f(X) :!> f(X*) on cherche plutôt un X tel que
f(X) = f(X*) ; Pour ce faire il faut trouver un zéro de la fonction T(X), dite
du Tunnelier ainsi définie :
f(X) - f(X*)
T(X)=
r
IIX - Xmll)o TI IIX - Xi*lI)J
i=1
où X*I. X*2 .....•X*r sont les derniers minimiseurs locaux trouvés.
• le terme f(X)-f(X*) assure l'élimination des points tels que
f(X) > f(X*), qui ne pourraient être minimiseurs globaux .
• le terme IIX - Xi*IIÀi permet d'éviter de choisir encore les minimiseurs
locaux (déjà rencontrés):
X*i
i= 1, ... ,r
38
• le terme IIX - Xmll~ évite à l'algorithme (de résolution de l'équation
T(X)=O) de converger vers un point stationnaire de
f(X)-f(X*). qui n'est
IIX - Xi*lIÂ.i
pas solution de T(X)=O. Xm est une constante à déterminer par une
heuristique à chaque étape de résolution de T(X)=O.
2.3,3 Méthodes
Indirectes:
Elles
utilisent,
en
plus
des
informations
locales
(valeurs
de
fonction) d'autres informations, pour construire un modèle convergeant
vers un minimiseur global.
2.3.3,1.. Mélhodes d'approximation des surfaces de niveau
Dans ces méthodes, l'optimisation est faite via des approximations
des surfaces de niveau. Nous avons peu étudié ces méthodes.
2.3.3,2" Méthodes d'approximation de la yaleur de fonction
Dans ces méthodes, l'optimisation se fait via une théorie basée sur
un modèle statistique de la fonction objectif f.
Elles sont utiles quand la
fonction objectif est très coûteuse à évaluer; par exemple, minimiser une
fonction dépendant de données sur un siècle (températures
moyennes
journalières, ... ); ainsi, si de
plus
le problème de
minimisation est à
refaire avec d'autres données (un autre siècle dans notre exemple), en
passant par un
modèle
statistique de
la
fonction
objectif basé
sur
quelques données (et non pas sur toutes), on minimise ce modèle qui, en
théorie, donne des solutions qui sont approximativement celles de f.
Les méthodes indirectes, même si des théories simples arrivent à
les justifier sont d'une implémentation numérique très difficile. En fait,
nous ne les citons ici que pour mémoire car dans notre cas, nous traitons
des fonctions définies explicitement ou dont les évaluations se font à des
coûts raisonnables.
39
2.4.
Caractéristiques
générales
des
méthodes.
Nous venons de voir que plusieurs
méthodes ont été proposées
pour trouver un minimum global, utilisant diverses philosophies.
Mais
en
dépit
de
leur
diversité,
on
peut
remarquer,
d'une
façon
générale, quelques traits caractéristiques communs.
Au
paragraphe
2.3.,
nous
avons
détaillé
délibérement
la
philosophie
profonde
de
certaines classes de
méthodes
pour
mieux
apprécier la diversité et l'ingéniosité des idées.
Cependant, on peut remarquer que toutes ces méthodes utilisent l'une
des deux grandes classes d'approche stochastique ou déterministe. C'est
pourquoi
nous
retiendrons
comme
classes
de
méthodes
ces
deux
principales classes (cf. 2.2.) :
approche
stochastique regroupant les Méthodes
Stochastiques,
approche
déterministe regroupant les Méthodes
Déterministes.
En effet, c'est en classant les méthodes numériques d'optimisation
globale de la sorte, que nous pouvons voir que pour chaque classe, on
peut
retrouver
des
schémas
généraux
auxquels
obéissent
les
algorithmes, et la nécessité numérique de travailler sur un compact S.
2.4.1. Optimisation sur un compact S
Rappelons qu'il y a deux classes de problèmes, ceux sans contraintes (on
minimise
f sur
Rn) et ceux avec
contraintes (on
minimise sur un
ensemble C dit des
contraintes).
Dans
les deux cas,
on
ramène le
problème
à
une
minimisation
sans
contraintes
sur
un
compact
S,
contenant un minimiseur global comme
point intérieur à S pour les
raisons
suivantes:
1) D'abord cela nous assure que la fonction est bornée et qu'une
solution existe. Ces deux conditions sont fondamentales pour tous les
théorèmes (ou des fois les postulats!) de convergence.
2) Pour les méthodes déterministes:
- certaines de ces méthodes partitionnent (ou quadrillent) le domaine de
recherche. Cela est facilité et simplifié sur S.
3) Pour les méthodes stochastiques:
-
pratiquement
toutes
(à
quelques
exceptions
près
dont
le
Recuit
Simulé), échantillonnent le domaine de recherche en un nombre fini de
points "représentatifs". Dans ces cas, que le problème soit avec ou sans
contraintes,
un
compact
S
comme
domaine
de
recherche
est
plus
convenable.
40
- la mesure de Lebesgue du domaine de recherche est, dans
beaucoup
de cas, utilisée pour donner la probabilité de succès des méthodes. El la
probabilité de ce succès est inversement proportionnelle à cette mesure.
Ainsi si le problème est sans contraintes, Rn ne peut pas convenir; car sa
mesure est infinie; ce qui donnerait une probabilité de succès nulle. Et
dans le cas d'un problème avec contraintes un compact S où aucune
contrainte n'est violée est préférable numériquement.
4)
Toutes
les
applications
pratiques
de
l'optimisation
globale
proviennent d'un constat, par l'utilisateur, de l'existence de
plusieurs
minimums. De ce fait, ce même utilisateur peut imposer des "limites" au
domaine
de
recherche
(ou
d'exigence).
Il
pourrait
donc
choisir
un
compact S.
2.4.2 Schémas eénéraux:
A première vue on peut dire que tout algorithme d'optimisation globale
obéit au schéma suivant:
schéma
1
phase 1: Trouver un minimum de f.
phase 2:
Décider s'il est global; dans ce cas s'arrêter;
Sinon faire la phase 3.
phase 3
Appliquer une procédure pour en échapper puis refaire
la phase 1.
Cependant,
en
distinguant
les
méthodes
suivant
les
deux
classes
d'approches
stochastiques
et
déterministes
on
peut
retrouver
des
schémas généraux plus spécifiques.
41
cas des méthodes déterministes
on peut citer quelques cas:
- avec Contions de pénalités:
schéma
2
phase 0: Minimiser f.
phase 1: Définir une fonction T auxiliaire associée à f (ayant des
minimums communs avec
1) permettant d'éliminer
tout
minimum
déjà
trouvé.
phase 2: Trouver un minimum de T.
phase 3: Décider s'il s'agit d'un minimum global dans ce cas s'arrêter;
Sinon appliquer la phase
1.
La phase
1 est la clé de tels algorithmes; voir comme illustration la
méthode du Tunnelier.
- autres cas :
schéma
3
phase 1: Décomposer le domaine en sous-régions en utilisant
éventuellement
les
propriétés de
r.
phase 2:
Eliminer
certaines
sous-régions.
phase 3:
Décider si les résultats ainsi obtenus sont satisfaisants, et
dans ce cas s'arrêter; sinon réactualiser le domaine de
recherche et
refaire l'étape
1.
42
cas des méthodes stochastiQues
Ces
méthodes, qui semblent donner plus de succès, peuvent se
décomposer en deux sous-groupes:
- le Recuit Simulé , les Recherches Aléatoires à une seule Initialisation:
schéma
4
phase 1: Trouver un itat de minimum.
phase 2:
Ne rester dans cet itat que si c'est celui d'un mimmum global;
Sinon s'en échapper (grâce aux
acceptations par
exemple
pour le Recuit cf. 3.1) et refaire la phase 1.
- les autres méthodes stochastiques (dont celles de regroupement)
schéma
S
phase 0:
(phase
globale)
Choisir un échantillon de N points.
phase 1: (réactualisatlon de
I'échantillon)
On peut éliminer des points et/ou en rajouter d'autres grâce
à des
procédures propres à chaque algorithme.
phase 2:
(phase de
raffinement)
On applique une phase comportant
éventuellement une
optimisation locale (ou quelques uns de ces pas) à certains des •
points de I'échantillon réactualisé (ou à tous).
phase 3: (critère d'arrêt)
Si le critère n'est pas vérifié retourner à la phase 2:
Sinon appliquer une
procédure complète d'optimisation
à chaque (ou certains) point(s) du dernier échantillon; le
plus
petit minimum est global.
43
Dans ce schéma général, on peut en effet regrouper pratiquement toutes
les
méthodes
stochastiques
(dites
aussi
à
stratégies
aléatoires)
en
résumant les différentes possibilités pour chaque phase:
• phase d'échantillonnnage :
• Cardinal de l'échantillon:
- un nombre fixé à l'avance,
- un nombre calculé au fur et à mesure par l'algorithme.
• Stratégie d'échantillonnage:
- points engendrés suivant une distribution quelconque,
- points engendrés suivant une distribution tenant
compte du voisinage du dernier point,
- points engendrés suivant une distribution tenant
compte de tous les points dejà engendrés.
• phase de réactualisation de l'échantillon
• nombre de points à rajouter ou à réduire de l'échantillon:
- un nombre fixé à l'avance;
- un nombre calculé au fur et à mesure par l'algorithme.
• phase de raffinement
• sans aucune recherche locale,
• recherche locale en certains points; la recherche pouvant
être une optimisation locale (ou quelques uns de ces pas);
• recherche locale en tous les points.
•
critère
d'arrêt:
cette
phase
est
délicate
pour
tous
les
algorithmes;
il est traité en 2.5.
44
2.5.
Dirficultés
Numériques.
L'implémentation
d'un
algorithme
numérique
d'optimisation
globale comporte beaucoup de difficultés. Ces difficultés, essentiellement
numériques, sont liées, de près, au type de problème qu'on veut traiter,
au
type
de
matériel, au
type
d'algorithme,
mais
aussi
au
type
de
fonction test dont on dispose pour valider le code informatique ainsi
écrit.
2.5.1 Critère d'arrêt Conyere,ence des méthodes,
C'est l'une des questions les plus délicates dans l'implémentation
d'un code. En effet, même quand des théorèmes assurent la convergence
d'un algorithme, ils ne donnent pas une caractérisation d'un minimum
global.
De ce fait,
tout
théorème ou toute propriété de convergence n'a un
intérêt que seulement pour le développement de la théorie.
En dehors de quelques cas de caractérisation d'un minimum global (voir
2.5.4), la convergence numérique d'un code est testée suivant un des
nombreux critères suivants:
• on se fixe un nombre maximal d'évaluations de f;
• on se fixe un nombre maximal d'itérations;
• pour les algorithmes stochastiques, on se fixe une certaine
probabilité à atteindre;
• si aucune amélioration (à un seuil fixé) n'est trouvée depuis un
certain temps;
• si les derniers itérés ne diffèrent pas d'un seuil donné.
Les
critères
ci-dessus
sont
adaptables
à
tout
type
d'algorithme
(stochastique
comme
déterministe).
Il
existe
d'autres
critères
de
convergence
beaucoup
plus
spécifiques
à
certaines
sous-classes
d'algorithmes.
critères propres à des ale,orithmes déterministes;
- pour les algorithmes utilisant une fonction auxiliaire T associée à f;
un critère d'arrêt plus spécifique pourrait être basé sur les résultats de
la fonction T. C'est le cas, par exemple, de l'algorithme du Tunnelier: où
l'on arrête l'algorithme dès qu'on ne trouve plus un zéro de T.
critères propres à des ale,orithmes stochastiQues:
- pour les algorithmes de regroupement, s'il ne reste plus qu'un seul
point par grappe, ou s'il ne reste plus qu'une grappe.
- pour certains algorithmes stochastiques on peut utiliser les mesures de
sous-régions d'attraction des minimiseurs locaux.
45
2.5.2. Difficultés liées aux contraintes ou à la densité des minimiseurs.
En pratique, il existe également, des difficultés spécifiques au type de
problème envisagé:
- difficultés liées aux contraintes:
algorithmes
déterministes
avec
fonction
auxiliaire:
ces
algorithmes adaptent quand
cela est
facile
(contraintes de bornes par
exemple) des
pénalités qu'on intègre dans
une
fonction
auxiliaire T
associée à f;
algorithmes stochastiques: ce sont les algorithmes qui sont les
plus sensibles au type de contraintes.
D'abord, elles ne peuvent pas traiter des contraintes de type égalités car
la mesure (de Lebesgue) de l'ensemble des contraintes est, en général.
nulle et cela incompatible avec certaines hypothèses de convergence. En
effet, comment imaginer se déplacer "au hasard" tout en restant sur des
courbes
(en
dimension
2)
des
surfaces
(en
dimension
3)....
sans
"déborder" par rapport à l'autre (ou aux autres) dimension(s).
Néanmoins
pour des
contraintes de
bornes
ou
pour des
contraintes
faciles (polyèdres, ensemble C des contraintes convexe) on peut adapter
des
sous-programmes qui
tiennent compte des
contraintes
sans
que
l'algorithme ne perde de son efficacité.
- difficultés liées à la densité des minimiseurs:
quand les minimiseurs locaux sont très proches et/ou quand l'amplitude
entre
les
différentes
valeurs
minimales
de
f,
est
faible,
beaucoup
d'algorithmes ont de sérieuses difficultés à trouver un minimum global.
2.5.3 Propriétés minimales pour un code d'optimisation elobale
Théoriquement,
un
code
d'optimisation
doit
converger
vers
un
minimum. Mais, l'expérience montre que. pour en arriver là. il Y a toute
une
heuristique.
En
effet pour les
cas
pratiques
précis,
des
codes
spécifiques sont généralement développés en connaissance de cause.
Pour le cas général, il n'y a pas encore. de code sûr. De plus. il n'y a
aucune convention sur ce qu'il faut utiliser comme test pour valider un
code qui se propose de traiter ne serait-ce qu'une classe
de fonctions.
11
n'y
a pas
longtemps encore,
les
seules
tests disponibles
étaient
académiques [DIXON & SZEGO 1975], rDIXON & SZEGO 19781 [HOC8l]. ou
[SCH87]. Récemment seulement. dans [FLOUDAS & PARDALOS 1990] on a
proposé des tests issus des applications: chimie énergétique. mécanique.
aéronautique.
De toute façon, jusqu'à présent. on ne s'est mis d'accord sur aucun test,
ni sur une performance pour valider un code d'optimisation globale.
46
Nous proposons néanmoins quelques propriétés minimales:
1) Que le code marche pour au moins une classe bien définie de
problèmes.
2) Que le code ne se fasse pas "piéger" par le premier minimum
local
rencontré.
En
d'autres
termes
le
code
doit
marcher
pour
un
problème (de la classe prévue) ne comportant que deux minimums.
3) En fin d'exécution le code doit, au pire, donner un minimum,
même si c'est un local et non un autre point (stationnaire, etc ... ).
4)
Le
code
doit
proposer
un
choix
par
défaut
de
tous
ces
paramètres de façon qu'un utilisateur puisse l'exécuter sans peine.
5)
A défaut de
trouver un
minimum
global,
on
doit pouvoir
"améliorer
une
situation
existante"
("améliorer"
un
minimum
local
trouvé par une autre méthode par exemple).
2.5.4.
Caractérisation d'un minimum global
En optimisation locale, on peut reconnaître un rmmmum local pour
les cas
où f est deux
fois différentiable. Des conditions d'optimalité
existent dans ce domaine: conditions nécessaires, conditions suffisantes,
etc.,;
De
plus,
certaines
de
ces
conditions
sont
faciles
à
adapter
numériquement dans un code: module de gradient inférieur à un seuil,
matrice hessienne définie positive, ...
Il en est tout autrement en optimisation globale. Pour le moment, et à
notre connaissance, aucun code satisfaisant basé sur une reconnaissance
analytique n'a été proposé.
Cependant, dans
le développement théorique de l'optimisation globale,
plusieurs
articles
proposent
de
résoudre
(théoriquement
pour
le
moment)
ce
problème en proposant des
conditions
nécessaires
voire
suffisantes
pour
certains
types
de
problèmes.
C'est
le
cas
dans
[HIRIART-URRUTY 1989] où les types de problèmes considérés sont:
• (Pl): Optimisation D.C.: minimiser une différence de deux fonctions
convexes sur un convexe fermé .
• (P2):
Maximisation
convexe: maximiser une fonction convexe sur
un convexe.
47
• (P3): Optimisation convexe "à rebours": rrumrmser une fonction f
convexe sous contraintes: g(x) positive et x dans C; où C est convexe et g
une fonction convexe.
Les résultats qui y sont présentés restent, rappelons-le, pour le moment
seulement
théoriques.
D'autres
caractérisations
ont
été
proposées
via
des
méthodes
numériques. C'est le cas dans [PINKUS 1968] où on utilise l'intégration
numérique pour caractériser le maximiseur global unigue sur C, de la
fonction continue f, comme égal à:
l x exp (-t f(x» dx
cl
lorsque t ~ +00
exp (-t f(x» dx
c
On a aussi la caractérisation suivante (provenant de l'étude des "grandes
déviations" en Statistiques):
avec comme hypothèses:
• f admet un maximum global unique
• C adhérence d'un ouvert borné non vide
max f(x) = Hm
log l exp ( -n f'(x) dx
lorsque n ~ +00
n
C
Mais, même pour ces deux cas, il n'est pas prouvé qu'en utilisant
l'intégration
numérique,
on
simplifie
ou
facilite
la
résolution
du
problème.
48
2.6.
Parallélisation.
2.6.1. Motiyations
Les algorithmes que nous avons cités dans ce chapitre, sont tous
prévus
pour
être
exécutés
séquentiellement,
pas
par
pas,
sur
une
machine à un seul processeur.
Actuellement, avec le développement de
l'informatique, pour certaines
applications, on peut faire exécuter par
plusieurs processeurs certains
codes, à condition qu'ils vérifient quelques règles de parallélisation.
L'optimisation
globale,
qui
est numériquement
un
problème difficile,
pourrait bénéficier des avantages de ce traitement en parallèle.
Brièvement, nous allons résumer les pour et les contre, de
l'utilisation
d'une
parallélisation .
• POUR:
• Amélioration de la vitesse de calcul
- les machines conventionnelles sont à leurs limites physiques;
- une machine avec N processeurs est moins chère que N, avec
un seul processeur chacune;
• Possibilité de résoudre des problèmes réputés complexes;
• Il y a des problèmes "parallèles" de par leur nature;
• Contraintes réelles de temps (délais ... ).
• CONTRE:
• faible utilisation (sous-utilisation) de la machine;
• Programmation plus compliquée.
2.6.2 Machines parallèles
On peut diviser les machines pour un traitement parallèle, en deux
catégories principales [SPITERI 1989]:
• les LS.D.M.:
Instruction Simple Données Multiples:
Composées
d'un
processeur
maître et de
plusieurs
processeurs dits
esclaves
faisant un type de
travail supervisé par le maître chez qui
transitent les
informations principales.
• les LM.D.M.:
Instructions Multiples Données Multiples:
49
Composées
de
processeurs
couplés
entre
eux
de
façon
à
pouvoir
communiquer entre eux (tous à la fois ou en grands sous-groupes) pour
échanger des informations de lecture, d'écriture et de messages. Dans ce
cas les processeurs ne font pas forcément le même type de travail. Ils
pourraient même
avoir besoin constamment des
résultats des
autres,
d'où la nécessité du dialogue.
Performances
Les performances d'un algorithme parallèle sont appréciées d'après les
valeurs de la vitesse-efficace, du gain de vitesse et de l'efficacité:
vitesse-efficace
V m:
Soit lm
le temps nécessaire pour une machine de m processeurs pour
traiter
parallèlement
l'algorithme;
et
Il
le temps mis par la même
machine pour traiter l'algorithme en non-parallèle;
la vitesse-efficace V m est donnée par
V m =.!.L
lm
Je gain de vitesse G m '
Gm=~
m
l'efficacité E m :
E
=
V rn
m
m . lm
2.6.3. Parallélisation
En optimisation globale, il y a des algorithmes qui sont parallélisables,
facilement, sous leur version initiale actuelle.
Ce sont essentiellement, des sous-classes d'algorithmes stochastiques:
les
Initialisations
Multiples,
les
Algorithmes
de
Regroupement,
par
exemple. En effet, selon le nombre de processeurs dont on dispose, il
suffit
de
leur
partager
l'échantillon
initial;
chacun
des
processeurs
traitera
le
domaine
qui
lui
est donné,
car
il
peut
se
passer des
informations des autres dans un premier temps.
50
Pour
plusieurs
autres
méthodes
stochastiques,
des
algorithmes
de
parallélisation sont proposés modifiant légèrement les algorithmes dans
leur version initiale.
Quant aux algorithmes déterministes, il est moins évident de leur
trouver une parallélisation. En effet, leur exécution exige une mise à jour
continuelle
des
paramètres,
des
approximations,
qui
dépendent
essentiellement des
résultats des
itérations précédentes;
donc, on
ne
peut
imaginer (dans
la
plupart des
cas)
un
processeur
pouvoir
se
charger de plus d'un pas de travail à la fois.
2.6.4. Conclusion
Les
développements
de
l'informatique
permettent
d'envisager
la
parallélisation
d'algorithmes;
cela
pourait
être
très
utile
pour
l'optimisation
globale.
Cependant,
il
faudrait
que
les
algorithmes
candidats
à
la
parallélisation
puissent
justifier
de
meilleures
performances. Par performance, on entend amélioration du nombre de
succès (cas où
un minimum est
trouvé) par
rapport à une
version
séquentielle.
Enfin
rappelons que la parallélisation en elle-même n'apporte pas une
garantie suplérnentaire de succès. Au plus, elle pourrait faciliter la tâche:
temps de calcul réduit, problème à dimension élevée, etc ...
Et de toute façon,
il faut que l'algorithme parallèle soit convergent et
vérifie les conditions minimales d'un code d'optimisation globale (cf.
2.5.3)
2.7.
Quelques
applications
de
l'optimisation
globale.
L'optimisation globale est très sollicitée dans
tous
les domaines
d'ingénierie. Cela s'explique, d'un côté, par un besoin de trouver des
paramètres qui donnent un choix optimal , et d'un autre côté, par le
constat
que
la
fonction
objectif
à
minimiser
comporte
plusieurs
minimums. On comprend donc, assez bien, pourquoi presque toutes les
branches de conception et d'ingénierie se trouvent confrontées à des
problèmes d'optimisation globale.
La plupart des articles qui proposent des codes, les testent sur des
exemples de fonctions
académiques; et la valeur du code n'était jugée
que
par
rapport à ces
tests. Quant aux
applications pratiques, elles
faisaient l'objet d'une élaboration d'un code spécifique. Les articles qui
les
présentaient.
permettaient
rarement.
à
un
autre
utilisateur,
de
pouvoir disposer de tous les paramètres de la fonction à minimiser, pour
pouvoir la tester sur un autre code. Cela fut ainsi jusqu'à la publication
51
de [FLOUDAS & PARDALOS 1990] où des
applications pratiques sont
proposées de façon explicite, par une fonction réelle et ses contraintes.
Nous
mentionnerons
dans
ce
paragraphe,
quelques
exemples
d'applications:
- les moindres carrés,
- la reconnaissance des formes
- la restauration d'images
2.7.1
Les moindres carrés
Le
problème
des
moindres
carrées
non
linéaires
est
très
souvent
rencontré en
pratique (ajustement des
paramètres d'un modèle
sur
un
ensemble de points expérimentaux). Il y a des codes disponibles pour
résoudre "localement" ce problème (dans les librairies scientifiques NAG
ou
IMSL); le choix d'une méthode dépend de données telles que: le
nombre d'inconnues, nature du problème ds moindres carrés (faibles ou
grands résidus) (cf. [ABDALLAH & BENCHEKROUN 1988]).
Que faire pour "améliorer" un minimum local? Tout d'abord il est très
difficile de llY.!Ù!. qu'on est dans un état de minimum local (lesquels
peuvent être
nombreux et rapprochés).
On
peut essayer
un
nouveau
point de départ, Le. une autre initialisation de l'algorithme ...
Cependant le problème des
moindres carrés est résolu facilement
et directement si on dispose d'un code d'optimisation globale.
En effet soit à résoudre le système d'inconnue X, suivant:
(5)
fi(X)= 0
pour i=l, ... , m
La résolution du système (S) se ramène à un problème d'optimisation
globale:
(P)
Minimiser f
m
avec f (X) =
L ( fi(X) ) 2
i=1
Si on trouve que le minimum global de f vaut zéro, alors ses arguments
sont solutions du système (S).
52
Nécessité de J'optimisation ~Iobale;
Toute solution x* de (S) est un minimiseur global de f, et f(X*)= O.
Si
l'optimisation
globale
nous
fournit
un minimiseur global
dont la
valeur de fonction est nulle, ce minimiseur est alors solution de (S).
Et, si le minimum global que l'on trouve est strictement positif, on peut
dire que le système (S) n'admet pas de solution.
Cependant,
si la
fonction f
admet des minimums locaux (dont les
valeurs
de
fonction
sont
alors
strictement
positives),
connaître
un
minimiseur local de f ne nous apprend rien sur les solutions de (S). En
effet on ne peut même pas s'en servir pour dire que le système admet
ou pas des solutions.
Conclusion:
l'optimisation
globale
permet
de
résoudre
de
façon
satisfaisante le problème des moindres carrés (Le. le système (S» : elle
permet de savoir si le problème n'apas de solution; et quand il en a, elle
trouve une solution, au moins.
11\\ ustration:
Soit le système (S) d'inconnue (X1. X2), suivant:
pour résoudre (S), on résout le problème (P) d'optimisation globale (ce
qui correspond à la fonction test Fü2 (cf. 3.1»:
f(X1,X2) - (-13+X1-2X2+5X22_X23)2+(-29+X1-14X2+X22+X23)2
sans contraintes.
cette fonction admet un rmmrnurn local
en Xlocal .. (11.4128, -0.89681)
pour
flocal ..48.98425
Xopt- (5,4)
pour fopt-O
• en utilisant le Recuit Simulé on trouve:
x* = (5 ; 4)
X* minimiseur global avec
f*= f (X*) = 0
donc X* est solution de (S);
53
en utilisant comme code d'optimisation locale le D.F.P (cf. chi) on
trouve:
XO = (lIAI; -0.89)
et
f"= f (XO) = 48.98 XO minimiseur (local)
XO n'est pas solution!
Nous
avons
utilisé
cette
technique
pour
résoudre
des
systèmes
d'équations non
linéaires, se ramenant à un problème d'optimisation
avec plusieurs minimums. Et cela avec succès pour le Recuit, car si un
minimum
global existe
il vaut
zéro;
on a donc
une
estimation du
minimum qu'on utilise dans
le test d'arrêt qui est, rappelons-le, une
étape délicate pour les
algorithmes stochastiques dont
fait
partie
le
Recuit.
2.7.2 Restauration d'imaees et Reconnaissance de formes:
Les résultats du Recuit Simulé dans Rn [BONNEMOY & HAMMA
1991a]
sont
complétés
par
des
applications
de
cette
méthode:
la
restauration
d'images
et
la
reconnaissance
de
formes
d'après
leur
luminance (cf. [BONNEMOY 199Ib)).
2.8.
Conclusion,
Les codes d'optimisation globale sont de plus en plus demandés; et
cela dans beaucoup de domaines de conception et d'ingénierie.
Plusieurs méthodes sont proposées pour répondre à celte demande: des
méthodes spécifiques à une situation donnée,
des méthodes à portée
limitée à une classe de fonctions données, des méthodes utilisant les
acquis en optimisation locale, des méthodes basées sur l'observation des
comportements naturels (sur la physique, les statistiques, l'évolution),
etc.
Malgré leur nombre important, on peut classer ces
méthodes en deux
catégories principales, selon qu'elles incluent ou pas, dans leur approche
des stratégies probabilistes: la classe des méthodes stochastiques et la
classe des
méthodes déterministes.
La convergence des algorithmes basés sur les premières méthodes est
assurée en
probabilité; celle des
algorithmes basés sur
les
méthodes
déterministes
est
implicitement
liée
aux
stratégies
déterministes
uti lisées.
Pour toutes ces classes de méthodes, il subsiste encore des problèmes
pour la validation de codes. En effet, il n'y a encore aucun test général
reconnu décisif en dehors des fonctions tests de type
académique; et
54
surtout
il
n'existe
aucune
convention
établie
sur
des
critères
de
performances à vérifier.
Les
méthodes
stochastiques
semblent
être
préférées
par
plusieurs
auteurs, car pour le moment, ce sont celles qui donnent le plus de
satisfaction.
Une de leurs faiblesses, cependant, réside dans
le
test
d'arrêt qui conditionne d'ailleurs, en grande partie le succès de tout code
stochastique (cf. 3.1 pour exemple).
Les problèmes de succès des méthodes déterministes, quant à eux, sont
beaucoup plus liés aux stratégies employées pour échapper à un état de
minimum local. En effet, ces stratégies, même si elles sont bien justifiées
théoriquement
sont
d'une
implémentation
difficile
et
surtout
leur
complexité croît, au fur et à mesure,
que se déroule l'exécution des
codes (cf. 3.2 pour exemple).
55
2.9.
Références
bibliographiques.
[ABD88] A. ABDALLAH & F. BENCHEKROUN, Méthodes et algorithmes
pour les problèmes de moindres carrés non linéaires, Mémoire de projet
D.E.S.S. "Méthodes Informatiques et Modèles Mathématiques", Université
Paul Sabatier Toulouse (1988).
[ARC84] F. ARCHETTI & F. SCHOEN, A survey on the global optimization
problem:
general
theory
and
computational
approaches, Annals of
Operations Research r, (1984).
[BOG84] P. T. BOGGS, R. H. BYRD & R. B. SCHNABEL editors, Numerical
Optimization,
section
three:
Global Optimization, Proceedings of the
SIAM Conference on Numerical Optirnization, Boulder, Colorado, June 12-
14, 1984.
[BON9Ia] C. BONNEMOY & S.B. HAMMA, La méthode du recuit simulé:
optimisation globale dans Rn, APII, Aut. Prod. Inf. Ind., vol 25, N°5,
1991, p. 477-496.
[BON91 b]
C. BONNEMOY, La méthode du recuit simulé: restauration des
images,
reconnaissance
de
surfaces, APII, Aut. Prod. Inf. Ind., vol 25,
W5, 1991, p. 497-517.
[COL87] N. E. COLLINS, R. W. EGLESE, B. L. GOLDEN, Simulated Annealing,
an Annoted Bibliography, Working Paper Revues MS/5, 87-028, College
of Business and Management, University of Maryland, (1987).
[D1X75] L.C.W. D1XON & G.P. SZEGO, (eds), Towards Global Optimization,
North-Holland, Amsterdam (1975).
[DIX78] L.C.W. D1XON & G.P. SZEGO, (eds), Towards Global Optimization 2,
North-Holland, Amsterdam (1978).
[FLOU90] C.A. FLOUDAS & P.M. PARDALOS, A Collection of Test Problems
for
constrained
Global
Optimisation
Algorithms, Lecture Notes in
Computer Science, Springer- Verlag (1990).
[FLOU91] C.A. FLOUDAS & P.M. PARDALOS editors, Recent Advances in
Global Optimization, Princeton series in computer science, (1991).
[GOM78]
J.
GOMULKA, Deterministic vs probabilistic approaches to
global optimizatlon, in: Dixon & Szegô, (eds), Toward Global Optimization
2, North- Holland, Amsterdam (1978).
,.
n
[\\V\\:fl8j
B.
1-I1r::r~\\<. 1 c..oot.:.""} Sv~J........~\\. hOt~ _
O-""...... "'-~"'-(
1
(\\... \\:\\u.\\M...... t.:u
"'lot· ~cz.. V"f A3 1 N°L n., 8~.
56
[HAM88] S.B. HAMMA, Optimisation Globale, Mémoire de D.E.A. Maths.
Appliquées, Université Paul Sabatier Toulouse (1988).
[HIR85] J.-B. HIRIART-URRUTY, Generalized Differentiability, Duality and
Optimization for problems dealing with differences of convex functions,
Lecture Notes in Economies and Mathematical Systems 256, Springer-
Verlag (1985).
[HIR89] J.-B. HIRIART-URRUTY, From Convex Optimization to Nonconvex
Optimization,
Part J:
Necessary and sufficient conditions for
global
optima lity, in: Clarke & al., Nonsmooth optimization and related topics,
Plenum(l989), p. 219-239 ..
[HOR86] R. HORST, A general class of branch-and-bound methods in
global
optimization
with
some
new
approaches
for
concave
minimization, Journ. of Optirnization Theory &
Applications 51 (1986),
p. 271-291..
[HOR90] R. HORST & H. TUY, Global
Optimisation:
Deterministic
Approaches, Springer-Verlag (1990).
[LAA87]
P.M.J. LAARHOVEN & E.H.L. AARTS, Simulated Annealing:
Theory & Applications, D. Reidel Publishing Company, Holland, 1987
[LEON66] A. LEON, A classifled bibliography on optimisation, in: A.Lavi &
T.P. Vogl (eds.), Recent advances in optimization techniques (John Wiley
1966).
[LEV85] A.V. LEVY & A. MONTAL VO, The Tunneling algorithm for the
global minimization of functions, SIAM J. Sei, Stat. Comp. 6 (1985), p. 15-
29.
[MLA86] R.H. MLADINEO, An algorithm for finding the global maximum
of a multimodal, multivariate function, Mathernatical Programming
34
(1986), p. 188-200.
[PAR87] P.M. PARDALOS & J.B. ROSEN, Constrained global optimization:
Algorithms and applications, Lecture Notes in Computer Science Vol 268,
Springer-Verlag
(1987).
[PAR91] P.M. PARDALOS & T PHILLlPS, Global Optimization of Fractional
programs, Journal of Global Optimization VoU, N°2 (1991), p. 173-182.
[PIJ72] S.A. PIJAVSKIJ, An algorithm for finding the global extremum of
a function, USSR Comput. Math. & Math. Phys. (1972), p.57-67.
57
[PINK68] M. PINKUS, A closed form solution of certain programming
problems, Operation Research 16 (1968), p. 690-694.
[PIN83]
J. PINTER, A unified approach to globally convergent one-
dimensional
optimization
algorithms,
Report
IAMI-83.5,lnst.
App\\.
Math. Inf. CNR, Milan (1983).
[RAT88] H. RATSCHEK & J. ROKNE, New computer methods for global
op t imi zat i on, Halsted Press, Ellis
Horwood Series in Math.
and
its
Applications, New York, N.Y. (1988).
[RIN86] A.H.G. RINNOY KAN and G.T. TIMMER, Global Optimization,
Report 8612/A, Erasmus University Rotterdam (1986).
[SHU72]
8.0.
SHU8ERT, A Sequential method seeking the global
maximum of a function, SIAM Journal on Numerical Analysis (1972), p.
379-388.
[SPI89]
SPITERI, Super
Calculateurs, cours de D.E.A de Maths.
Appliquées, Université Paul Sabatier, Toulouse III (1989).
[STR78]
R.G.
STRONGIN,
Numerical
methods
of
multiextremal
optimization, Nauka, Moscou (1978).
[TIM84] G.T. TIMMER, Global Optimization: a stochastic approach, Ph. D.
Thesis, Erasmus University Roterdam (1984).
[TOR74]
A.A. TORN, Global Optimization as a combination of global and
local search, Abo Akademi, HHAA A:13, Finland (1974).
[TOR89]
A.A. TORN and A. ZILINSKAS, Global Optimization, Lecture
Notes in Computer Science, Vol 350, Springer-Verlag (1989).
[ZIL86] A. ZILINSKAS, Global Optimization - Axiomatics of statiscal
models, algorithms and their applications, (Mokslas, Vilnius 1986), 166
pages en Russe.
58
chapitre
3:
EXEMPLES
D'ALGORITHMES NUMERIQUES
D'OPTIMISATION GLOBALE
59
3.1. LA METHODE DU RECUIT SIMULE 1
C. BONNEMOY·
et
S. B. HAMMA00
BESUME:
Dans ce qui suit nous présentons des résultats numériques
d'optimisation de fonctions dans
Bn,
certaines
sans
contraintes,
d'autres avec
contraintes. Nous donnons d'abord un aperçu sur les
méthodes numériques d'optimisation globale, pour situer la méthode du
Becuit Simulé. Puis nous résumons quelques résultats de l'utilisation de
cette méthode. notamment l'influence de certains paramètres: le point
initial,
la
température
initiale...
Ensuite
nous
montrons
comment
l'introduction de cycles dits de réchauffement améliore les résultats et
atténue l'influence des paramètres "de contrôle". Enfin nous comparons
le Becuit Simulé simple puis le Becuit avec les cycles de réchauffement
avec
deux
méthodes
d'optimisation
(l'une
déterministe,
l'autre
stochastique).
ABSTBACT:
One
presents
in
the
following
some
numerical
results
of
Optimization in Bn
with or without constraints. First one resumes
some results about this method practice, particulary the use of the
model's parameters such as the initial point, initial temperature.. After
that one shows that the introduction of "heating cycles" improve those
results and minimize the parameters influence. Finally one compares
the annealing methods (with or without cycles) with two others (one
deterministic, the other stochastic).
(1)
Ce travail sur le Recuit Simulé a faill'objet d'une publication dans RAIRO. APII,
vol 25, n05. 1991.
(0)
Professeur à l'Université d'Auvergne. Clermont-Ferrand 1.
(") Chercheur Labo. Analyse Numérique Univ. Paul Sabatier, Toulouse III.
Le travail de S.B.HAMMA a été financé par le CNES Toulouse ( marché N°873-CNES·87-4937).
60
3.1.0.
HISTORIQUE DES METHODES NUMERIQUES D'OPTIMISATION GLOBALE
Soit une fonction F: Rn -+ R
à minimiser;
Plusieurs approches ont été utilisées pour résoudre numériQuement
le
problème
d'Optimisation
Globale.
(voir
[TOR89]
pour
une
bibliographie assez complète). Cependant, en dépit de l'abondance des
méthodes utilisées, on peut les classer suivant deux grandes classes
d'approches: déterministes ou stochastiques.
3 1 0 1 Méthodes stochastiQues
Elles
utilisent des
processus ou propriétés
probabilistes
pour
assurer la convergence vers un minimum global. Généralement, elles
utilisent l'évaluation de la fonction F sur un échantillon aléatoirement
choisi; puis il s'ensuit des procédures sur l'échantillon. L'avantage de
ces méthodes réside dans le peu d'hypothèses qu'on exige sur la
fonction
objectiif
F. Cependant,
leur
convergence
n'est assurée
généralement que par :
lim P ( xk(n) -+ X·) - 1
lorsque n -+ +00
où X· est un point de minimum globlal et xk(n) une suite de points
convergeant vers X·.
Or il faut bien fixer n fini.
3,1.0.2. Méthodes déterministes
Elles nécessitent, souvent, des hypothèses supplémentaires
sur la
fonction objectif F: condition L-lipschitzienne, classe C1 ou C2, ce qui
limite leur champ d'applications. Pour plus de détails se référer à
[TOR89].
3.1.1.
PHILOSOPHIE
Le
nom du
Recuit
Simulé dérive
d'une analogie
avec
des
procédés
physiques
comme,
par exemple,
amener
un
système
(de
corps) à une énergie d'équilibre (E). Partant d'un état initial avec pour
énergie (E), le système est perturbé selon une
loi
aléatoire vers un
nouvel
état
(voisin)
pour
lequel
on
peut
évaluer
(facilement)
la
variation
d'énergie
(li E)
résultant
du
changement
d'états.
Si
le
changement d'états correspond à une réduction d'énergie du système,
le
nouvel
état est
alors
accepté
(validé).
Cependant,
même
si
le
61
changement a entrainé une augmentation de l'énergie du système, il
est
(éventuellement)
accepté
quand
même
avec
une
probabilité
spécifique
(qu'on
se
précise).
La
probabilité
d'acceptation
est
généralement de la forme :
P = Exp(- !! E/T)
où T est un paramètre de contrôle correspondant à une fonction de la
température dans
l'analogie citée.
Voilà en résumé la philosophie qui fait qu'un algorithme de
Recuit Simulé est apte
à converger vers un minimum global.
L'analogie,
ci-dessus
présentée,
a
été
introduite
par
[METROPOLIS et al 1953]. Une application à l'Optimisation Globale
discrète a été proposée, par la suite, par [CERNY 1982] et [KIRKP ATRIK
et al 1983].
Toutes les méthodes du type Recuit Simulé sont à classer dans
le groupe des méthodes stochastiques.
3.1.2.
ALGORITHME GENERAL DU RECUIT SIMULE
3.1.2.1
Algorithme de Métropolis
Avant de présenter l'algorithme général du Recuit Simulé, il est utile
de rappeler l'algorithme standard des méthodes
à recherche aléatoire
dont fait
partie
le
Recuit Simulé. Ce sont des
méthodes classées
stochastiques.
Alrocithme
standard
des
méthodes
de
descente
aUatoire
llJI....U
lilLl : choisir Xo initial, évaluer F(XO) ; itération k= 0
lilù:
à l'itération k+1. Engendrer un vecteur X à partir de Xk.
Si F(X) < F( Xk ) alors Xk+ 1 = X sinon Xk+ 1 = Xk
laL..J: si critère d'arrêt non vérifié faire k=k+ 1 et aller à Pas 2.
Si
critère
d'arrêt
vrai,
arrêter.
Alors
(Xk, f( X k»
est la solution
approchée.
62
Des théorèmes de convergence vers un rmnimum global ont été
établis pour ces méthodes ([BON8?]), en utilisant par exemple pour le
tirage de X une loi normale centrée sur Xk et de matrice de variance-
covariance oZ x 1n <In matrice unité de dimension (n x nj)
avec CJ2
"suffisamment
grand".
Cependant
des
difficultés
subsistent
pour
échapper aux minimums locaux, ce que fait naturellement le Recuit
Simulé.
6:
3.1 2 2
Algorithme du Recuit Simulé.
Procédure
Recuit;
Début
Initialisation de l'état x dans
Rn;
Calcul de T(O) (* initialisation de la température *)
Tant que critère d'arrêt non vérifié faire
01 :- 0 (* initialisation du compteur de répétitions *)
Répéter jusqu'à ce que 01 • N(k)
générer état y voisin de x; (* transition de x vers y *)
!1F. F(x)-F(y) ; ( * due à la transition
de x à y *)
appel
sous-programme accepter F,T(k));
Si accepte est vrai alors
x :. y (* "état y devient l'état courant *)
fin
si;
01:=01+ 1
fin (* répéter jusqu'à ce que 01 • N(k) *);
k :.. k + 1;
calculer T(k);
fin tant que (*
critère d'arrêt
"):
Fin (* procédure Recuit *);
subroutlne Accepte t\\ F,T) :
Début
Si !1F:s;; 0 alors
accepte :.. vrai; (* descente *)
sinon
A(T ,k):= Exp(-(!1 F/T( k»;
si random(O,1) <
A(T,k)
alors
accepte ;. vrai; (* remontée acceptée *)
sinon
accepte :'" faux; (* remontée refusée *)
fins i;
fin si;
Fln (* subroutine Accepte *);
N.B.
En mécanique aléatoire la fonction objectif F est appelée énergie, et le
paramètre T(k) température.
64
3,1,3.
APPLICATION AU PROBLEME D'OPTIMISATION
Soit f:
Rn
-+
R, mesurable et lnt-bomée, à minimiser sur un
ensemble de contraintes C; où C est une partie de
Rn de mesure non
nulle (ce qui interdit les contraintes de type égalité).
Soit donc:
trouver
m .. inf f(X)
Xe C ,
( P )
et pour E >0
arbitrairement petit,
trouver X de C tel que f(X) < m+E
ce qui est équivalent à :
( P' )
trouver m= inf F(X)
pour tout
X e Rn.
où F(X)=f(X)+IC(X)
avec IC(X) .. I 0 si. X e C,
+ 00 sinon.
En ce cas F n'est pas bornée supérieurement; or l'hypothèse "F bornée"
est requise pour la convergence ([ROY89]); mais on remarque que
Accepte( ~ F,T) est faux dès que Xe C mais pas Y.
Donc toute transition hors de C sera refusée par l'algorithme du
recuit; ce qui revient à dire qu'on n'utilisera que la restriction à C de
la transition.
Pour satisfaire l'hypothèse il suffit:
o
soit que f soit bornée sur C,
o
soit de restreindre la recherche
à E n C contenant au moins une
solution de
(P), ainsi f est bornée sur E n C, où E est un pavé borné
suffisamment grand de Rn.
Les transitions hors de E seront soit refusées, soit réfléchies dans E
(moins de tirages perdus cf [HAA89]).
65
3 1 3 1
IlIustratiQn
Soit à minimiser la tonction suivante:
FONCTION FOO
f(XloX2) - 1/3 (X13 + X23) - X1 2 - X22 + 8/3
SQUS contraintes :Xi ~ 0
i-1,2
r S r1 QU r2S r S r3
QU r ~ r4
avec
r.. (X1 2 + X22)112
r1 -0.25 r2-0.75
r3-1.25 et r4..1.75
L'ensemble C des contralntes n'est pas connexe.
On part du point Xo-(O,O), le minimum glQbal vaut 0 en X*..(2,2).
La IQi de tirage correspond à la IQi normale N(X,a212) avec 0-0.2.
La figure 3.1.1., ci-après. montre le tonctionnement précis de la
méthode du recuit
sur cet exemple.
La tralectclre depuis Xo jusqu'à x* (noté S) passe par tQUS les peints
de pseudo-descente (notés +) y inclus par les vrais points de descente
QU descentes absolues (notés x), et par les polnts de remontées
acceptées (notés -).
Les autres points (notés 0) sont des remontées refusées. y inclus tQUS
les points tirés hors des contraintes C.
L'écart-type 0 doit être initialement suffisamment grand pour que les
deux couronnes puissent être franchies.
On peut
améliorer la préclsfon du résultat en cholslssant un test
d'arrêt effectif (ici la solution est supposée ccnnuel) et en faisant
décroîtra 0 vers une valeur strictement positive lorsque k~ +00 •
Ainsi, avec le seul test d'arrêt k~500. soit un maximum de 500 tirages
pour obtenir une vraie descente, et avec 0 divisé par 10
à k=50 • 100 et 250. on obtient X*..(2.02,1.95)
Y*..0.002 en 19
descentes et 281 tirages.
66
Résultat de la minimisation de FOO ( voir Figure 3.1.1)
avec une précision de E-2
DONNEES INITIALES :
nombre maximal de tirages pour une vraie descente:maxtir=500
nombre maximal de vraies descentes :
maxdes=lOO
ecart-type pour le noyau de transition:
(J =0.2
températue initiale :
T(O) =1
point initial:
Xo =(0, 0)
RESULTATS
l\\. chaque vraie descente
:
notée X
(19 au total)
k compte les tirages
(281 au total)
sf compte les pseudo-descentes:
notéees +
(61 au total)
.
rm compte les remontées acceptées:
notées
(54 au total)
sq compte les remontées refusées:
notées 0
(166 au total)
ct compte
les tirages hors de C:
notées 0
(122 au total)
N" Descente
F(X)
X(l)
X(2)
k
sf
rm
ct
sq
1
2.612
0.143
0.194
11
1
0
10
10
2
2.168
0.195
0.815
78
6
8
64
64
3
2.095
0.078
0.896
1
1
0
0
0
4
2.009
0.249
0.933
11
1
0
9
10
5
1. 713
0.412
1.141
2
1
0
1
1
6
1. 627
0.975
0.723
26
6
5
7
15
7
1.143
0.529
1.792
76
20
19
27
37
8
0.746
0.922
1.955
1
1
0
0
0
9
0.644
1.030
1. 909
22
6
7
0
9
10
0.469
1.470
1.462
17
6
5
4
6
11
0.346
1. 676
1.443
1
1
0
0
0
12
0.341
1.572
1. 533
1
1
0
0
0
13
0.145
2.108
1.609
1
1
0
0
0
14
0.097
2.277
1. 880
3
1
0
0
2
15
0.064
2.216
2.117
12
2
5
0
5
16
.060
2.224
1.918
3
1
1
0
1
17
0.049
2.041
1.773
3
1
1
0
1
18
0.046
2.113
1.812
3
1
0
0
2
19
0.002
2.028
1.955
9
3
3
0
3
SOLUTION atteinte en:
x· =( 2.028,1.955)
pour F(X·)=0.002
67
..
figure 3.1.1.: détails de la minimisation de FOO.
68
3.1.3 2. Choix des paramètres
En analysant "algorithme général présenté plus
haut,
on voit
facilement ce qui différenciera les différents algorithmes du type
Recuit Simulé : choix de la fonction température T(k),
du nombre des
répétitions
N(k),
choix
des
critères
d'arrêt
(en
fonction
des
itérations, de la température, des acceptations, ft.).
Nous présentons donc les différents choix faits à travers la très
large bibliographie sur le Recuit Simulé.
a) Fonction température T(k)
Rappelons que le but de cette fonction est de permettre, grâce
à la probabilité P=Exp( - ~ F 1 T(k)), l'acceptation des remontées.
Accepter
trop
de
remontées
entraînerait
que
l'algorithme
converge très
lentement, et, au contraire, peu de
remontées
pourraient ne pas assurer la convergence vers l'optimum global.
C'est
pour cette
raison
que
plusieurs
choix
de
fonctions
température ont été proposés :
1) fonction constante
T(k) .. cte
2) fonction arithmétique T(k) .. T(k-1) - cte
3) fonction géométrique
T(k)..a(k)*T(k-1) (a souvent constant)
4) fonction inverse T(k)..
Cte/(l +d*k)
5) fonction logarithmique T(k) = CR / ( log(l +k) ), celle pour laquelle
la convergence théorique est assurée ([GEM 88]).
b) Nombre de répétitions N(k)
Ce nombre correspond au nombre maximum de tirages qu'on peut
tolérer entre deux descentes successives. Si ce nombre est petit, on
acceptera peu de remontées; donc les chances de converger vers
l'optimum global pourraient être réduites. Au contraire, s'il est grand,
la convergence pourrait être très lente. les différents choix sont:
1) constante N(k) .. cte
2) arithmétique N(k) .. N(k-1) + cte
3) géométrique N(k) = N(k-l) 1 a(k)
4) logarithmique N(k) .. CRI log(T(k))
5) puissance N(k) = (N(k-1)) (1/a(k)
69
c) Critère d'arrêt.
Un des points délicats du Recuit Simulé réside dans le choix et la
formulation d'un critère d'arrêt. La convergence de l'algorithme est
assurée en probabilité pour k -+ +00; cependant il faut bien fixer un
nombre maximal de tirages qui garantisse toutefois une convergence.
D'autre
part,
à quelles conditions sur les derniers points itérés
saurait-on qu'on a atteint l'optimum global ?
Pour répondre à toutes ces questions, plusieurs critères ont été
utilisés, séparément ou souvent ensemble. Ce sont:
1)
critère itérations : nombre maximal de tirages total fixé.
2)
critère température : T(k) s T finale donnée.
3)
critère énergie : 1Â.FI inférieure à un seuil donné .
4)
critère acceptations : si peu d'acceptations sont enregistrées
depuis un certain nombre d'itérations.
5) critère estimations : pas de meilleure estimation de "optimum pour
un certain nombre de températures successives, etc.
3.1,4,
UTILISATION DES CYCLES PE RECHAUFFEMENT
3,1 4.1, Choix des paramètres pour le code REC FOR
Les différents paramètres que nous avons choisis pour le code
REC.FOR (basé sur le recuit) sont dus à [BONal] .
a) fonction température T(k):
Nous avons choisi la fonction positive à décroissance lente
T(k) - CR 1 Log(kt + 2)
où kt est le cumul de tirages
où CR est la constante du recuit correspondant à T(O) • CR.
Pour CR, nous avons enregistré de meilleurs résultats quand CR était
voisin de 2, pour les exemples présentés .
70
b) Nombre de répétitions N{K):
Un nombre maximal constant Maxtir a été choisi entre deux
descentes consécutives.
c) Loi de tirage Gen{Xis.JQ.:
On passe de Xk à Xk+ 1 en tirant un vecteur X à partir de Xk suivant
une loi normale N(Xk, Ok) (où Ok" a2xln ' a2 variance). Bien que la loi
normale donne des résultats satisfaisants il y a bien sûr beaucoup
d'autres choix possibles (uniforme. exponentielle, ...).
d) Test
d'arrêt:
Il est basé à la fois sur:
1) les acceptations : s'il n'y a plus d'acceptations depuis un certain
nombre de tirages.
2) les estimations : si les dernières estimations de Fa pt ne
diffèrent pas de E donné.
3) les itérations : si un nombre total de tirages est atteint et/ou si
entre deux descentes consécutives, un maximum de tirages Maxtir est
atteint.
e) Problèmes avec contraintes d'optimisation
Si des contraintes existent (du type inégalités seulement), elles
sont intégrées de la façon suivante:
après chaque itération, on vérifie
grâce à un sous-programme CONTR(n,x) si Xk vérifie les contraintes.
Sinon on ne calcule même pas f(Xk), on passe directement à l'itération
suivante.
3.1.4.2.
Utilisation des cycles de réchauffement.
Le code REC.FOR. tel que décrit dans ce qui précède. qui
consiste en l'application de la procédure RECUIT en une fois à partir
d'un point initial XO, donne des résultats satisfaisants. Cependant pour
un nombre de variables du problème (~) élevé ou pour des fonctions
très "vallonnées", nous avons voulu améliorer le REC.FOR en partant du
dernier optimum trouvé comme nouveau point initial XOc après, bien
sûr. une réduction éventuelle de la température initiale TOC du cycle
précédent. Les cycles ainsi créés ont pour avantage, en repartant d'une
température relativement élevée et d'un point initial assez "affiné".
de donner l'optimum global avec plus de chances. Cependant ils ont
l'inconvénient de majorer le nombre d'évaluations de la fonction
(nombre déjà très élevé pour chaque cycle). Mais cet inconvénient est
71
déjà atténué si on fait un choix de réduction de la température initiale
(T OC pour un cycle c) de façon
que les acceptations de remontées
soient moins fréquentes quand le nombre de cycles tend vers
Ns (Ns
nombre maximal de cycles).
3 1 4 3
Algorithme du RECUIT-CYCLES
Procédure RECUIT_CYCLES;
Début
Initialiser c:=1; Ns (nombre de cycles);
Xo1(*1 er point initiaJ*);TO 1(* 1e re température initiale*);
Tant que ( c S Ns et critère d'arrêt non vérifié)
faire
RECUIT; (* appel procédure. fournit Xopt C *);
Xopt :- XoptC;
c := c-t:
TOC :- Rc * TOc-1 (* réduction de température initiale *)
XOc :- Xopt (* initialisation de Xo pour cycle suivant *)
fin (* tant que *)
(Xopt. F(Xopt))
est solution;
FIN (* procédure RecuiCcycles*)
72
Organigrammes simplifiés
a) le RECUIT_CYCLES
INITIALISATIONS
Ns (nombre de cycles)
c .. 1
X 0 1 (premier point initial)
T01
Test d' arrêt général
oui
c<Ns
r n o n _ _-l
(Xopt. F(Xopt))
est
oui
solution de (P)
73
b) Le RECUIT
INITIALISATlONS
Xo
F(XO)
T(O)=CR
. t k
etape
- - - - . - - - - - - -
ou i _~__!__--,
(X, F(X)) est
solution de (P)
appel à GEN (Xk.X)
X est engendré
à partir de xk
F(X) < F(Xk)
oui
Appel à
Accept((F(X)-F(Xk),T(k»
k.. k+1
oui __--1
mise à jour
de T(k)
74
3.1.4.5. Convergence du Recuit Simulé.
Rappelons que le problème (P) qui nous intéresse est:
Trouver x*, f* tels que
f(x*) s: f(x)
'if x E S
f(x*) .. f*
où f: S ~ R,
S compact de Rn.
La méthode du Recuit consiste à construire une suite de variables
aléatoires
XO, Xl, " , Xk, ".
de la façon qui suit:
Soit Tl, T2, ". , Tk, ".
une suite de nombres positifs tels que
Hm Tk .. 0;
k-4oo
Xo initial étant choisi, pour Xk .. x on engendre un Yk voisin avec une
probabilité
P[Yk" Y 1 Xk .. xl
..
R(x,y)
où R est une matrice de
transition telle que R(x,Y) > 0 si et seulement si y appartient au
voisinage de x.
On définit ainsi
Yk avec probabilité Pk
Xk+l"
Xk avec probabilité 1-Pk
1
avec Pk" inf{ t , exp(-(f(Yk) - f(x)) 1 Tk) }
La suite de variables aléatoires
X .. (Xk : k ~ 0) ainsi construite est
une chaîne de Markov à temps discret .
a)
Convergence dans Je cas où f S ~ R
avec
S fini
Donnons quelques définitions relatives à la suite
de variables
aléatoires (Xk : k 2: 0):
75
propriété SI : (Strong Irreducibility : forte irréductibilité).
On dira que la suite (Xk : k ~ 0) vérifie la propriété S.1. si et seulement
si:
Dans S tout état Xi est accessible depuis un état Xj c'est-A-dire
3 une suite
Xj. Xi1, Xi2•... , Xip.Xj
telle que
R(Xik,Xik+ 1) > 0 ,
VOs k < p.
Etats accessibles A une hauteur E
Pour tous Xi et Xi de S,l'état Xi est accessible depuis l'état Xi à
une hauteur E si et seulement si:
3 une suite Xi • Xi1. Xi2, ...• Xip • Xi telle que
R(XiK.Xik+1) > O.
VOs K < p
et f(Xik) s E
pour 0 S k s p.
propriété W.B : (Weak Reversibility : faible réversibilité).
On dira que la suite (Xk : k ~ 0) vérifie la propriété W.R.. si et
seulement si:
v
E e R.
v Xi, Xj e S.
les assertions suivantes
sont
équivalentes :
i) l'état Xi est accessible depuis "état Xj A une hauteur E.
ii) l'état Xj est accessible depuis l'état Xj A une hauteur E.
Minimum local, Profondeur d'un minimum.
Un état x est dit minimum local s'il n'existe pas d'état x', tel
que f(x') < f(x), accessible depuis l'etat x A une hauteur f(x).
Sa profondeur est le plus petit nombre E > 0
tel qu'un état x'. avec f(x') < f(x). soit accessible depuis l'état x à une
hauteur f(x)+E.
76
profondeur d'un Minimum global
Par définition la profondeur d'un minimum global est fixée à +00.
y
tg
---1,
x
Xlocl
Xg
Xloc2
figure 3 1 2 :
profondeur d'un minimum pour f: R -+ R
la difficulté d'échapper d'un état minimum est d'autant plus
grande que ce minimum est profond.
Soit d* la profondeur d'un minimum local (et non global) la plus
grande. Soit S* l'ensemble des points pour lesquels f prend sa valeur
minimale (globale). On a le théorème suivant :
77
THEORGME : ( cas fini, cf (HAJ88])
Si les conditions SI et WR sont satisfaites, alors les assertions
suivantes sont équivalentes
i) Iim P[Xk e SOl - 1;
k ~oo
ii) L Exp( - dO 1 Tk) _
+00.
k •
1
De plus, si la fonction température Tk est définie par Tk - C 1 Log(k+1)
alors
ii) est réalisée
si et seulement si
C ~ dO.
Remargues:
° Le choix de R symétrique se justifie par le fait que cette
condition est suffisante pour réaliser la propriété WA.
° En général, il n'est pas facile de trouver dO (ou un de ses
majorants). Aussi le choix de C peut entrainer une convergence très
lente si C est trop grand ou un risque de ne pas converger vers
l'optimum global si C est trop petit.
b) Convergence dans le cas où f: Rn ~ R
Dans ce paragraphe, nous rassemblons un certain nombre
de
résultats concernant la méthode du recuit simulé dans un espace
d'états continu (cf [ROY89), [HAA89J). Tous les résultats sont donnés
en fonction de F
fonction associée à f, qui prend en compte les
contraintes.
Probabilité de BoUzman
Soit (E, @, À) un espace mesuré polonais (espace des "états"), avec
À a-finie.
Soit F: E ~ R une fonction objectif ("énergie") à minimiser sur E, telle
que exp(-FIT) soit À-intégrable.
78
On définit la probabilité de Bollzman pour F. sur E, par:
~r - (1IZr) exp(-FIT) dÀ où Zr - JE exp (-F(x)fT) dÀ(x)
où T est une "température" >0. avec Zr fini pour au moins une valeur
T-TO.
On
suppose.
pour
l'instant,
que
F
est
essentiellement
bornée
inférieurement par m , ce qui est équivalent à:
À{xlF(x) < m} - 0
et V asü, Ea-{xlF(x) < rn-a)
est de mesure non nulle pour À.
Proposition 1
On pose E· - {xlF(x) s m},
si Fest semi-continue inférieurement. alors E· est fermé et
si ~r --? ~ faiblement quand T --? 0, alors Il est portée par E·.
Donc la probabilité de Boltzman définie par
F se concentre
sur
l'ensemble des minimums de F quand la température T décroît vers o.
La probabilité limite Il existe et peut être explicitée dans des cas
pratiques.
Ainsi:
• si E est fini, alors ~ est la loi uniforme sur E· (cf. [GEM88
)).
o
si les hypothèses du théorème de convergence du recuit H1, H2, H3 et
H4 (voir plus loin), sont satisfaites, alors (cf [HAA89J) :
a) si F a un minimum absolu x· (c'est-à-dire pour tout voisinage
W X. de x", il existe a>O tel que Â.({xlF(x) s m+a}/Wx.)-O).
alors I.l - li(x·).
b) si E est un compact de Rn muni de la topologie induite.
si À est la mesure de Lebesgue restreinte à E.
si E· - {xl' x2' ...• xr} est inclus dans P(intérieur de E),
si V xI F(x) - gj(x-xi) où les fonctions gi: Rn --? R sont homogènes de
degré di>O. avec dr-d2..... dq>dq+J>- ... dr
avec 1~:s:r,
alors
~ =(1N) t vI li(xj) , où Vj .. À{x/gj(x):!::1} et V - t vI .
j;J
j;1
o de même,
on a le résultat voisin suivant dans [ROY89]:
si À est la mesure de Lebesgue sur une partie E de Rn.
79
si Ea,,{x/F(x) < rn-a] est compact pour un a>O,
si EO 2 (x1, x2, ...• xr) inclus P(intérieur de E),
si F est deux fois différentiable en chaque point xi de EO,
alors
III converge faiblement vers Il définie par:
r
r
Il =(1/K) r kj Ô(Xj)
où kl 2 (det(F"(xj))) -1/2 et k - r kj
j:1
j:l
L'algorithme du recuit simulé
Ayant probabilisé "espace des états E par les mesures de Boltzman
III associées à F, on construit une chaine de Markov (Xk) k~ non
homogène dans E, définie par des noyaux markoviens NTk(x,dy) tels que:
V f mesurable, E [f(Xk+1) / Xk) - [NTkfJ (Xk) et que les lois successives
Pk des éléments aléatoires Xk soient définies par Pk+1 - Pk NTk.
Le résultat fondamental obtenu est que
IIPk - 11H11 ~ 0 si
Tk~ 0 "assez lentement" quand k ~ +'"'.
Ainsi on aura simulé à la limite une variable dont la loi est portée par
EO, l'ensemble des minimums de F.
Construction des novaux Nn, dynamique de Métropoljs:
° 1) On se donne un noyau Markovien N sur E. dit noyau de
"proposition". tel que >.. soit réversible par rapport à N, c'est-à-dire
V A1, A2 e@
Jd>"(X) N(X,A2) =
Jd>"(X) N(X,A1)
Al
A2
par exemple N(x,dy) - n(x,y)dy avec n symétrique.
° 2) Pour une température T on introduit le noyau de Métropolis
Nr(X,dy)
1
2
EX (y) N(x,dy) Ar(x,y) + 1EXC (y) Al(x,y) N(x.dy) +
+ Ôx(y)
J(1-Al(x,y)) Nl(x,dy)
EXC
où Ex .. ( y e E tel que F(y) S F(x)) et Al(x,y) 2 exp (-[F(y)-F(x)J+/T)
80
AT(x,y) est la probabilité "d'acceptation" de la transition de l'état x
vers l'état y, obtenue par le noyau N(x,dy), cette probabilité vaut 1 si
F(y) < F(x), sinon elle décroît exponentiellement en F(y) à T fixée, et ce
d'autant plus rapidement qua T est proche de o.
En remplaçant x par l'observation xk da Xk, et T par Tk à l'instant k,
le premier terme correspond à une descente en y toujours acceptée, le
second
à une remontée refusée avec la probabilité ATk(Xk'Y)' le
troisième à une remontée refusée avec la probabilité 1-ATk(xk'Y) (et
dans ces cas Xk+\\ - xk).
Convergence:
Proposition 2:
Si ~T est réversible pour NT alors ~T est la loi asymptotique pour le
noyau NT.
Théorème (cf [ROY89))
Si H1) (E, @, l) est un espace masuré avec la-finie,
H2) F:E ~ R est bornée,
H3) l est réversible pour le noyau de proposition N,
H4) il existe q e N", une probabilité P sur E et un réel c>O tels que
V x
Nq(x,.) 1 C PO
("contraction" de N)
ALORS
il existe une constante Cr" telle que si Tk-C,Ilog(k+2), avec Cr ~ Cr",
alors IIP k - ~Tkll --> 0 quand k ~ 00.
(Pour la norme au sens de la variation totale des mesures et, par
exemple, Cr" - q ô(F) où ô(F) - sup ( IF(y) - F(x)l' (x,y) e Ex E}).
Commentaires
Nous avons donné ces résultats de convergence du Recuit simulé
dans Rn, surtout pour justifier que pour le cas continu la convergence
da
cette
méthode
a
été
traitée,
et
que
notre
démarche
de
l'implémenter est bien justifiée. Evidemment cela ne facilite pas, pour
autant, l'implémentation,. Mais cela, on l'a vu au chapitre 2, est
courant
entre
la
théorie
d'une
méthode
et
son
implémentation
numérique.
81
3.1.5.
RESULTATS
NUMERIQUES.
Nous
présentons des résultats numériques de plusieurs essais
effectués sur une dizaine de fonctions tests choisies dans (HOC81),
[SCH8?), [BON8?), [TOR89) et d'autres sources, référencées (AUTRES);
toutes ces fonctions tests sont du type académique.
Présenter
des résultats
numériques
résultant
de l'utilisation
d'un code d'Optimisation Globale comporte plusieurs difficultés. En
effet, s'il est courant de tester un code par des fonctions connues, il
est cependant difficile de convenir d'un critère de performances :
• nombre d'évaluations de la fonction?
• vitesse de convergence?
• précision sur l'optimum?
• nombre de succès (rapport global/local)?
C'est pourquoi nous présentons à la fois des résultats sur le nombre
d'évaluations de la fonction, sur la vitesse de convergence, le nombre
de succès et une comparaison avec un autre code d'optimisation.
3,1,5,1. Fonctions tests:
Pour chaque fonction nous donnons la source bibliographique, la
dimension N, le minimum global Xopt et la valeur f(Xopt) • fopt.
FONCTION FOl :
dimension N • 2 ; Réf. [BON8?]
( différente de FOO du 3.1.3.3.)
f(X"X2) • 1/3 (X,3 + X23) - X,2 - X22 + 8/3
sous contraintes : Xi ~ 0
i-l,2
r s 0.5
ou r ~ 1
avec r.. (X,2 + X22)'/2.
Cette fonction admet un minimum local en
Xlocal .. (.5/2 112 , .5/2112)
pour
fl ocal=2.4313
Xopt = (2,2)
pour fopt -O.
FONCTION F02 :
dimension N - 2 ; Réf - [BON8?)
f(X loX2) • (-13+X,-2X2+5X22-X23)2 + (-29+X,-14X2+X22+X23)2
sans contraintes.
Cette fonction admet un minimum local
en Xlocal ..(11.4128, -0.89681)
pour flocal -48.98425
Xopt= (5,4)
pour fopt-O.
FONCTION F03 :
dimension N - 2 ; Réf - (AUTRES]
f(X,.X2) - _(X22 + X,2.2X,+X2+1)
sur [-1,1) x
[-1,1].
82
les quatre sommets de [-1.1) x
[-1,1)
sont tous des minimums
dont 3 locaux et un seul global en (-1,1):
Xopt - (-1.1)
pour fopt - -S.
FONCTION F04 :
dimension N - 2 : Réf - [SCH87)
f(Xl,X2) - X2
sous contraintes :
X2 - X12 ~ 0
Xopt - (0,0)
pour fopt-O.
FONCTION F05 :
dimension N - 3 ; Réf - [SCH87)
f(Xl,X2.X3) - - X1"X2"X3
sous contraintes: O~ Xi et 72 - Xl- 2X2• 2X3 ~ 0
Xopt - (24.12.12)
pour fopt--345S.
FONCTION FOS :
dimension N-2 ; Réf - [TOR89)
Fonction de Bran;n pour N",2
f (x.y) - a(y - bx2 + cx - d)2 + h(l-e)cosx + h
sous contraintes : ·5 s x s 10
et
0 s y ~15
avec
a-l
b-5.1I4n 2 c- 5m d&S h.. 10 e&1I8n.
Cette fonction a trois minimums globaux en :
Xopt - (-3.14159 , 12.275)
Xopt - (3.14159, 2.275)
pour fopt -0.397887.
Xopt - (9.42478 • 2.475)
FONCTION F07 :
dimension N - 5 ; Réf., [HOC81)
f(Xl. ..... X5) - 2 - 11120 (X1"X2"X3"X.."X 5)
sous contraintes:
O~ Xi s i
i-l ....5
Xopt - (1.2.3.4.5)
pour fopt.,l.
FONCTION FOS:
dimension N - 3: Réf - [AUTRES)
f(Xl.X2,X3) -(X1-1 )(X l-2)(X l-3)+X3
sous contraintes: (Xi ~ 0
i-l ....3) et ><3 s 5
et
(><32 - X1 2• X22 ~ 0) et ( X1 2 + X22 + X32 - 4
~ 0)
Xopt - ( 0,21/2,21/2)
pour
fopt - -S + 2 1/2.
FONCTION F09:
dimension N - 1 : Réf - [AUTRES)
f(X) - X4.14X3+ SOX2·70X
sans contraintes.
Cette fonction admet un minimum local en :
83
Xlocal - 5.9572
pour
flocal - 11.9576
Xopt - 0.7808
pour
fopt
-
-24.3696.
FONCTION F10 :
dimension N = 10 ; Réf - [AUTRES]
N
f(X) -
- L (X12 + CIXI)
CI-0.5
i:I
sous contraintes:
-1 S Xi S 1
i-1 •..•N.
La résolution de ce problème correspond à
la maximisation de
fonction convexe sur un polytope:
MAX (Xj2 + CIXI)
sur [-1 • 1 ] N avec
0 < CI < 1
f admet donc. 2N minimums locaux stricts tous sur les sommets
de [-1 • 1] n dont un seul global:
Xopt- ( 1, ..... , 1 )
pour
fopt - -15.
FONCTION F11 :
dimension N - 2 ; Réf - [TOR89]
Sjx-hump camel-back function
f(X,.X2) - 4 X,2 - 2.1 Xl. + 113 X,6 + X,X2 _ 4 X22 + 4 X2.
sous contraintes:
·5 S
Xi s 5
i..1,2.
Cette fonction
admet 3 paires de minimums. une pour chacune des
valeurs de fonction suivantes:
f1- -1.0316 •
f2 - -0.2154
et
f3 - 2.1042
X10pt - (0.08983 • -0.7126) ou
X20pt - (-0.08983 • 0.7126 )
pour
Icpt - -1.0316.
FONCTION F12 : dimension N - 2 ; Réf - [TOR89]
Fonction de Rastrigin
f(X,.X2) - X,2 + X22 _ COS (18X,). COS (18X2)
sous contraintes: ·1 S Xi S 1
i-1.2.
Cette fonction admet environ 50 minimums
Xopt - ( O. 0)
pour
fopt - ·2.
FONCTION F13 : dimension N-4 ; Réf - [AUTRES]
f(X) = (XJ+1O X2) 2 + 5 (X3 - X.) 2 + (X2 • 2 X3)· + 10 (XI- X.)·
sans contraintes.
Xopt - ( 0.0.0.0)
pour
fopt - O.
POINTS INITIAUX OONNES PAR LESREFERENCES:
F01
(O,O);
F02: (-15,2);
F03: (O,O);
F04: (1.1);
FOS: (10.10.10);
F06
(O,O);
F07: (2,2.2.2.2);
Foe: (0.1.2.2.1);
F09: (4):
F10
(0•..•0);
F11: (0.0);
F12: (0.75,0.75);
F13 :(3.-1.0.1).
84
3.1 52 Influence du point initial Xo
Comme (presque) la plupart des méthodes d'optimisation, le Recuit
Simulé est très sensible au choix du point initial.
Dans le tableau suivant, nous présentons quelques résultats en
fonction des points initiaux choisis. Pour chaque fonction, nous avons
pris
la valeur de CR
qui
donnait
les
meilleurs
résultats
pour
l'ensemble des points initiaux (avec 02.0.04 et Maxtir ..500):
a) influence du point initial au cours d'un
seul cycle de réchauffement
Point
Initial
(0,0)
(8,9)
(-1.2,1)
(0,2)
(15,-2)
(5,6)
fonction
fOl( CR-2)
a
a
a
L
a
f02( CR-O.S)
L
a
L
a
L
a
f06( CR-O.S)
a
L
a
a
a
fll( CR-O.S)
a
a
a
a
a signille que le minimum global est allelnl
L signifie que seul un minimum local est atteint
• signifie que le code a convergé ailleurs
CR - température Inlliale
Remarque:
Pour tous les cas où le minimum global n'a pas été atteint
(cas L et cas e) pour d'autres choix de paramètres il a été atteint.
Exemples: pour f02 avec Xo·(O,O)
si on part de CR- 2
le minimum
global est atteint; pour fOB avec Xo..(8,9) en partant de CR - 2 il est
aussi atteint.
Il faut noter aussi que tous les Cas (e) correspondent au cas où Xo est
choisi en dehors du domaine de contraintes. Dans de telles situations
un choix judicieux des paramètres est indispensable pour
y revenir
L'utilisation des cycles pourrait servir à cela. Pour ce faire, au cycle
c, on prend comme point initial XOc, le dernier point X engendré
qui
est, quelques fois en dehors des contraintes (donc pas point de
minimum). Bien sûr, on conserve durant tous les cycles, le meilleur
point de minimum.
b) pour plusieurs cycles
Au paragraphe suivant. pour chaque fonction, nous n'avons pris en
considération que la valeur de CR pour laquelle nous avons obtenu les
meilleurs
résultats.
85
Avec les cycles de réchauffement, même si la température initiale est
mal
choisie,
on
peut trouver
un
nombre de
cycles
permettant
d'atteindre le minimum global. Evidemment ce nombre est d'autant
plus petit qu'on choisit un bon CR.
Dans le tableau qui suit, nous sommes partis de CR - 2 (pour le
premier cycle) pour toutes les fonctions.
Pour chaque point initial, nous donnons le nombre minimal de cycles
nécessaires pour atteindre \\e minimum global avec
CR - 2:
dimension 2
Point Inillai
XO
(0,0)
18 , 9 )
(·1.2.1)
(15,·2)
Fonction
FOt
3
•
F02
3
3
Fil
•
8
dimension 3
Point Inillai
XO
(10.10,10)
(0.1.2.2.1)
Fonction
F05
Foe
..
dimension 10
Point Inillal
XO
XI.O 1.1,10
X5.0 Xi.·lt.
1 <> 5
Fonction
Fl0
7
3.1.5.3
Influence de la température initiale CR
Rappelons que le choix
d'une température initiale trop
grande
entraîne l'acceptation de beaucoup de remontées, donc une convergence
lente, tandis que le choix d'une trop petite valeur, risque de "geler·
trop rapidement l'algorithme vers un minimum local.
Nous présentons, pour un point initial Xo (donné par la référence)
pour un cycle, les résultats que donne le Recuit Simulé en fonction de
CR (CR - TO), pour Maxtir - 500
et a2 - 0.04.
86
a) influence de CR sur un cycle
Température Initiale CR
CR-O.5
1.5
2
8
Fonction
FOl (Kt-600
a
a
a
a
L
F02 (Kt-BOO)
L
L
a
a
F04 (Kt-600 )
a
a
a
a
a
où
a signifie que le minimum global est alleint
L seul un minimum local est atteint
• le code converge aiUeurs
Kt cumul moyen (moyenne sur les CR) des tirages de Xk ; comprend
aussi les
évaluations de Jonctions.
Remarque:
Pour CR - 8,
on ne converge pas vers le minimum global pour F01,
F02. En elfet, théoriquement on devrait y converger mais lentement.
Mais entre temps, d'autres critères d'arrêt tels que nombre maximum
de tirages, critère acceptations, ... imposent un arrêt à l'algorithme. Et
souvent cet arrêt peut intervenir à un point quelconque. C'est pourquoi,
au cas où CR est pris très grand, il faut un peu relaxer certains
critères d'arrêt.
b) influence de
CR sur plusieurs cycles
L'utilisation de cycles permet d'éviter les inconvénients cités au
paragraphe précédent. En effet même si l'algorithme s'arrête en un
point (qui
n'est pas
global), en repartant de
ce
dernier et
en
permettant un nouveau cycle de réchauffement, on finit par converger
(même si c'est lentement) vers le minimum global.
Dans le tableau suivant, sont présentés les résultats suivants: nombre
de cycles en fonction de CR pour chaque fonction.
87
Les points initiaux étant ceux donnés par les références. on choisit
encore le même maxtir, et la même variance que ci-dessus:
Tempéralure Inillale CR
CR-O.5
1.5
2
8
Fonction
FOI
(Kt-1200)
en 1 cycle
2
F02
(Kt-6000)
4
3
4
F03
(Kt-12 000)
5
5
2
2
3
F04
(Kt-600) 1
F07
(Kt-9000)
3
2
Kt cumul moyen (moyenne sur les CR des tirages de Xk : comprend aussi les évaluations
de fonctions).
Remarque:
Le monbre minimal de cycles semble croître avec les températures
initiales extrêmes CR
pour les exemples ci-dessus. Cela pourrait
s'expliquer par la remarque précédente.
3.1.5.4
Nombre d'évaluations de la fonction coût et temps moyen de
convergence
Un des points faibles du Recuit Simulé est de nécessiter un
grand nombre d'évaluations de la fonction. Cela s'explique par le fait
qu'on accepte des remontées (donc d'autres évaluations) en plus des
évaluations
nécessaires pour
n'importe
quelle
méthode du
type
Recherche Aléatoire dont fait partie le Recuit Simulé. Ce nombre
d'évaluations est,
hélas, très élevé dans le cas d'utilisation des
cycles. Néanmoins on peut le réduire en faisant une réduction Rc de
température initiale à chaque début de cycle c.
Cela entrainerait que Nt < Ns • Nc où Nt • nombre d'évaluations pour Ns
cycles, Nc • nombre total d'évaluations pour 1 cycle.
88
Nombre d'éyaluations de la fonction coût:
Sous les mêmes conditions qu'au 5.1.3
pour un cycle :
pour plusieurs cycles :
fonction évaluations
fonction
évaluations
cycles
F02
600
FOI
3 000
3
F09
900
F02
5 000
4
FIl
300
F03
10 000
5
F12
1 000
FIO
50 000
7
Kt cumul moyen (moyenne sur une dizaine d'essais; comprend aussi les évaluations de
fonctions).
Temps de convergence
La procédure RECUIT-CYCLES écrite en fortran
F77
a été testée
sur IBM PS2 XT, sur Compatible AT(avec co-processeur arithmétique),
sur SUN 3 et sur CDC (Cyber 860) .
Pour toutes ces machines, les succès ont été les mêmes; seuls le temps
de calcul et la précision (10-4 , 10.6,10.6 ) ont différé.
Pour les fonctions tests choisies, le temps d'évaluation de chaque
valeur de la fonction est peu coûteux. Néanmoins avec le nombre total
d'évaluations, le temps total mis par le code pour converger crolt.
Le temps dépend. bien sûr, de la complexité (nombre de minimums
locaux, ...) de la fonction.
Ainsi, de quelques secondes (5 à 6 pour le XT) pour un cycle pour F01.
on atteint l'ordre de l'heure pour plusieurs cycles ( avec le XT) pour F10
[qui
est
un exmple-piège en optimisation globale]. Ceci
est très
significatif surtout si on sait que dans la pratique, l'utilisateur peut
avoir des fonctions très coûteuses à évaluer.
89
3 1 5 5.
Comparaison avec la méthode du polytope mouvant" (ct. 3 2)
Comme nous l'avions dit, comparer deux méthodes d'optimisation
est très
délicat. Ainsi dans ce qui suit, nous présentons juste les
résultats moyens donnés par chacune des méthodes pour les fonctions
suivantes. Pour 10 exécutions, on a noté entre parenthèses le nombre de
fois où l'on a obtenu a, l et •
Simplexe
Recuit 1 cycle
Recuit plusieurs cycles
F01
8(2)L(4}"(4)
8(9}L(1)
8(10) en 3 cycles
F02
8(6)"(4)
8(8)L(2)
8(10) en 4 cycles
F07
8(2)"(8)
8(8)"(2)
8(9) en 3 cycles "(1)
F10
a(2)L(8)
a(3)L(7)
8(10) en 7 cycles
RemarQue1 :
Pour
FOl,
le
cas
·(1),
correspondant
au
résultat
du
Recuit
(plusieurs cycles) n'est pas, en fait un vrai échec car le code converge
en un point qui
est proche du minimum global et qui n'est ni local ni
stationnaire. Seul un problème de précision dû au fait que le point de
minimum global est sur la frontière des contraintes expliquerait cet
échec apparent. Ainsi on trouve:
X - (0.98, 1.89, 2.79, 3.89, 4.78) pour f-1.20
Alors que Xopt - (1,2,3,4,5) pour fopt-1.
Remarque 2 :
En temps de convergence, le simplexe est nettement plus
rapide.
90
3 1 5 6 Comparaison avec une méthode de descente aléatoire ([BON81))
Cette méthode PRSAN a été implantée en Pascal sur SUN et sur
IBM-PC par [BON81) . Nous le comparons avec le Recuit sur des exemples
avec minimum local. Sur 10 exécutions (les meilleures pour le Recuit) :
nombre de succès
PRSAN
Recuit 1 cycle
Recuit plusieurs cycles
FOl
a( 10)
a(9)L(1)
a(10) en 3 cycles
F02
a( 1 0)
a(8)L(2)
a(10) en 4 cycles
F07
a(8)·(2)
a(8)·(2)
a(9) en 3 cycles .( 1)
Temps de convergence ( moyennes en secondes surCompatible IBM AT)
PASAN
Recuit 1cycle
Recuit plusieurs cycles
FOI
15 S&condes
5S
6S
F02
16 S
8S
21 S
F08
13 S
10 S
28 S
RemarQue 1 :
Pour les comparaisons ci-dessus aussi bien
avec le Simplexe
Mouvant qu'avec le PRSAN, nous n'avons considéré que les 10 exécutions
sur 20 qui ont donné les meilleurs résultats pour chaque fonction. Nous
nous sommes mis à la place d'un utilisateur qui essaie d'améliorer à
chaque exécution le résultat précédent en changeant de point initial
et/ou de paramètres de façon tout à fait naturelle.
Remarque 2 :
Par ailleurs, il faut noter que le Simplexe Mouvant n'est pas conçu
pour traiter des problèmes avec contraintes. Néanmoins nous sommes
arrivés à adapter un "petit sous-programme· de contraintes dans le
Simplexe; ce qui nous a permis de traiter FOl, FOI et Fl0.
De toute façon, ni les méthodes stochastiques ni le Simplexe ne peuvent
traiter des contraintes de type égalités.
91
3.1.6. CONCLUSION
Sur
quelques exemples d'optimisation dans
Rn. avec ou sans
contraintes,
nous
avons montré
les
possibilités
offertes
par
une
méthode particulière d'optimisation à stratégie aléatoire: le
Recuit
simulé sur un ou plusieurs cycles. Nous avons délibéremment laissé de
côté les aspects théoriques (fBON91bJ) pour nous consacrer à l'étude
expérimentale du choix des paramètres:
• bien que le point initial Xo
puisse
être
théoriquement choisi
quelconque et même hors des contraintes, on observe des cas où le
mrrurnum global n'est pas atteint en un cycle (mais
il l'a été pour
plusieurs cycles).
• le choix de la température initiale est crucial pour la méthode:
Si elle est trop petite, on peut bloquer sur un minimum local; trop
grande, la convergence vers un minimum global est plus lente. On
peut tenter de "évaluer automatiquement en se basant sur
la
première descente qui devrait être atteinte avec un pourcentage de
remontées acceptées situé dans une fourchette raisonnable.
• toutes les expériences présentées l'ont été avec un paramètre de
dispersion 02 constant et petit. Le fait de le prendre constant permet de
mieux étudier l'influence des autres paramètres; le fait de le prendre
petit provoque des transitions plus locales. Mais cela a permis de
mettre en évidence, l'influence positive des cycles de réchauffement.
Ces expériences ont été aussi effectuées en faisant varier 02. ce qui
permet d'obtenir des convergences plus rapides.
Par exemple en prenant 02-4 au lieu de 0.4 pour la fonction F03, on
atteint le minimum global en 1 cycle, 6 descentes et 536 tirages au (jeu
de 21 descentes et 1116 tirages.
De même, pour F02, le minimum global est atteint en 1 cycle, 13
descentes, 469 tirages
(le minimum local étant franchi lors de la 9è me
descente), alors qu'il avait fallu 3 cycles.
• les avantages de cette méthode sont ceux de toutes les méthodes
stochastiques: simplicité de mise en oeuvre sur le plan algorithmique
et
affaiblissement conséquent des
hypothèses sur
les
problèmes
traités. Les inconvénients sont aussi les mêmes: difficulté dans le
choix des paramètres pour lesquels l'expertise humaine joue encore un
trop
grand rôle,
et surtout des temps de calcul qui peuvent être
prohibitifs pour
nombre d'applications.
• un point essentiel est que les méthodes stochastiques, et plus
précisément le Recuit qui permet d'échapper aux minimums locaux,
92
permettent de résoudre, peut-être avec une précision insuffisante, des
problèmes peu accessibles aux méthodes déterministes. A ce titre tout
atelier de génie logiciel devrait posséder cet outil.
3 1 7. REfERENCES BIBLIOGRAPHIQUES
(AZE88) R. AlENCOTT,
Simulated Annealing, Séminaire Bourbaki 1987-
1988, N° 697 p. 223-237
(BON87) C. BONNEMOY, Sur quelques aspects de l'utilisation de méthodes
déterministes en milieu stochastique et inversement, Thèse d'Etat,
Université Clermont Il, Sept 87.
(BON91b) C. BONNEMOY, La méthode du Recuit Simulé: Restauration des
images. Reconnaissance des formes, A.P.U., vol.25, N°S, p. 497-517
1991.
(CER85) V. CERNY, Thermodynamical approach to the traveling salesman
problem: an efficient solution, Journ. Opti. Theory & Appl., 45 (1985)
P.41-45.
(COL87) N. E. COLLINS, R. W. EGLESE, B. L. GOLDEN, Simulated Annealing,
an Annoted Bibliography, Working Paper Revues MS/S, 87-028, College
of Business and Management, University of Maryland, (1987).
(HAA89)
H. HAARIO & B. SAKSMAN, On the definition & convergence of
the Simulated Annealing Process in general State Space, University of
Helsinki, November 1989 ( preprint )
(HAM88)
S.B. HAMMA,
Optimisation Globale, Mémoire de DEA, 88,
Université Paul Sabatier Toulouse lit.
(HOC81] W. HOCK et K. SCHITTKOWSKI, Test Examples For Non Linear
Programming
Codes, Lecture Notes in Economics and Mathematical
Systems, 187, Springer-Verlag, (1981).
(KIR83), S. KIRKPATRICK, C.D. GELATT Jr. & P.M. VECCHI, Optimization by
simulated annealing, Science 220 (1983) p.671-680.
(LAA87) P.M.J. LAARHOVEN & E.H.L. AARTS, Simulated Annealing: Theory
& Applications, D. Reidel Publishing Company, Holland, 1987
93
[NEL65]
J. A. NELDER & MEAD R. A simplex method for funclion
minimizalion , Comput. J. 7 ( 1965) , 308-313.
[MET53] N. METROPOLlS, A. W. ROSENBLUTH, M. N. ROSENBLUTH, A. H.
TELL ER & E. TELL ER, Equation of stet« cstcutstions br fast computing
machines, The Journal of Chemical Physics 21 (1953) p. 1087-1092.
[RIN86]
A.H.G. RINNOOY KAN and G.T. TIMMER, Global Optimization,
Report 86121A, ERASMUS University, Rotterdam (1986).
[ROY89]
G. ROYER, M~thode du Recuit SimuM, Cours à l'Université d'
Orléans, Février 1989
[SCH87]
K.
SCHITTKOWSKI,
More
Test Examples For Non Unear
Programming
Codes, Lecture Notes in Economics and Mathematical
Systems, 282, Springer-Verlag, (1987).
[TOR89)
A.A. TORN and A. ZILINSKAS, Global Optimization, Lecture
Notes in Computer Science, Vol 350, Springer-Verlag.
94
3.2. LA METHODE DU POLYTOPE MOUVANT:
OUTIL D'OPTIMISATION NON LINEAIRE.
Résumé:
Totalement différente de la méthode (la plus connue) du Simplexe
en Programmation Linéaire, la méthode du Simplexe de Nelder et Mead,
appelée aussi celle du Polytope Mouvant est une méthode d'optimisation
non linéaire très utilisé dans les milieux des applications.
Ne nécessitant aucun calcul de
gradients. elle utilise uniquement
des évaluations de fonction et des opérations géométriques sur l'espace
des variables (réflexion. expansion. contraction).
On part d'un simplexe de N+ 1 sommets (pour un problème de dimension
N), on fait subir à ses sommets des opérations de réflexion. d'expansion
ou de contraction de façon à trouver des directions de descentes et de
déplacer le simplexe dans cette direction. Le simplexe, au
cours des
itérations.
varie
de
forme:
la
taille
et
les
angles.
notamment.
sont
continuellement changés d'où son nom de polytope mouvant.
En théorie elle est conçue pour trouver un point de minimum local. Dans
la
pratique,
elle
se
comporte
(parfois)
comme
une
méthode
d'optimisation globale. D'où
son intérêt.
INTRODUCTION
Soit une fonction f: Rn -+ R
à minimiser.
Nous présentons ici la méthode du Polytope (ou Simplexe) Mouvant et
les modifications que nous y avons apportées afin de la faire converger
vers un minimum global de f.
L'algorithme classique (Le.
l'original) ne prétend pas
trouver un
minimum global. D'ailleurs l'algorithme parallèle
non plus (cf.3.2.2.2.).
Cependant en
apportant les modifications dans les directions d'avancée
de
l'algorithme
classique.
non
seu lement
on
assure
une
mei lIeure
95
garantie de convergence, mais aussi, une amélioration du minimum local
qui dans certains cas, donne le global.
Et de toute façon, la méthode du Polytope Mouvant nous intéresse pour
plusieurs raisons:
• c'est
une
méthode directe.
Donc on
peut
l'utiliser pour
un
problème non différentiable comme algorithme local.
• elle pourrait être utilisée pour la phase
locale d'un algorithme
d'optimisation globale (exemple du Tunnelier).
• on peut modifier, facilement, l'algorithme classique et multiplier
les directions de recherche de façon à améliorer un minimum local.
• elle peut être utile
pour
un algorithme du
type Initialisations
Multiples (cf. 2.3.2.1.) pour un problème non différentiable.
Pour
toutes ces raisons, au
moins, elle pourrait être
utile
pour
la
résolution du
problèmes pratiques.
3,2,),
PRINCIPE
Définition et Notations
La figure géométrique composée de N+ 1 sommets dans un espace
de dimension N est appelée simplexe (de dimension N). Un simplexe sera
noté S= < XI, X2, ... ,XN+I>.
Dans ce qui suit on notera X/o Xb ,Xllies sommets respectivement de
plus petite valeur de fonction î, de plus grande valeur et de deuxième
plus grande valeur de fonction f.
X2
X3~--~Xl
fj~yre3 2.1; Exemple de Simplexe dans R2
Principe
Le
principe de
la méthode du
Polytope Mouvant (Simplexe de
Nelder et
Mead) est de considérer les
N+l
sommets d'un simplexe
(pour un problème de dimension N) puis de se servir des N sommets de
plus petites valeurs de fonction pour s'éloigner de Xb (sommet de plus
grande valeur) et d'essayer de le remplacer par un point de plus petite
96
valeur.
Eventuellement,
on
remplacerait
aussi
Xg
(deuxième
plus
grande valeur). Ce faisant on converge "naturellement" vers des états de
minimum.
En
pratique on
peut décrire l'algorithme comme composé des
5
principales étapes
suivantes:
étape 1:
On choisit un simplexe initial S avec ses N+ 1 sommets de telle
façon
que
les
vecteurs joignant
X, aux
autres
sommets
engendrent
l'espace Rn.
étape 2:
On choisit comme direction de recherche, la direction joignant
le sommet
Xh (qui a la plus grande valeur de fonction) à l'isobarycentre
X o de tous les autres sommets.
étape 3:
Si celle direction est une direction de descente, on remplace
X h par un point de (XhX o). Ce point est choisi de façon à avancer un peu
moins (contraction) ou un peu plus (expansion); et cela en fonction des
résultats des opérations géométriques sur Xh par rapport à Xo
étape 4;
On réadapte le simplexe en lui remplaçant certains sommets
en fonction des résultats de l'étape 3, de façon qu'il se "déplace" vers les
directions de descente.
étape 5:
On teste une condition d'arrêt: si elle est vérifiée on arrête
l'algorithme; sinon on repart à l'étape 2.
A l'arrêt, le dernier simplexe est un simplexe réduit autour d'un point
candidat
minimiseur.
3.2.2.
Algorithmes
du
polytope
mouvant
Nous présentons d'abord l'algorithme tel qu'il a été décrit par ses
auteurs [NELDER & MEAD 1965), puis nous présenterons un algorithme
parallèle de [TORCZON 1989) et enfin notre version.
La différence entre les trois versions réside dans le choix de la direction
de descente. L'algorithme original effectue la recherche dans une seule
direction, celle de (XhX o); l'algorithme parallèle, quant à lui, effectue les
recherches dans les N directions engendrées en faisant une réflexion de
tous les points par rapport à X"
97
Rappelons
que
ces
deux
versions
sont
toutes
deux
conçues
pour
l'optimisation locale mais, après expérimentation numérique, nous avons
constaté les
phénomènes suivants:
• la version séquentielle (l'originale) enregistre quelques échecs de
convergence (locale). Cependant il lui arrive de s'échapper de certains
états de minimum locaux et de converger ainsi vers un minimum global.
• la version parallèle converge toujours, mais vers un minimum
local; nous n'avons enregistré aucun cas où elle arrive à s'échapper d'un
minimum local pour converger vers un global.
Alors nous avons adapté un algorithme "la recherche adaptée" qui
garde les avantages de la version séquentielle (possibilité d'échapper à
certains
minimums
locaux) et
qui
poursuit
les
mêmes
ambitions
de
garantie de convergence que
la version parallèle. La recherche adaptée
est
une
version
séquentielle
mais
avec
possibilité
de
recherches
multiples à certaines phases de son exécution. En effet, le nombre de
directions est modulé de façon à permelire à l'algorithme de converger.
X2
.....
1
..... .....
1
.....
1
---
---
..........
1
---
-':>t.:-
1
.....
1
Xl
1
1
X3
1
1
1
:\\
1
direction de recherche si
1
f(X2) est la plus grande valeur de fonction.
(jure 3,2.1.
Exemples de directions de recherche.
98
3.2.2.1. L'al&orithme séquentiel
Initialiser SO = < XI, X2,,,., Xn+I>. a. P, 'Y
Calculer f .. [2, ... , fn+1
Tant
que critère d'arrêt faux faire
Trouver XI, Xg, Xh, fi, fg, fh.
Déterminer le centre Xo des opérations géométriques
faire une réflexion de Xh par rapport à X o (soit Xr son image)
si la réflexion est réussie
alors faire une
expansion
si l'expansion est réussie remplacer X h par Xe
sinon remplacer Xh par X r
sinon
faire
une
contraction partielle
si elle est réussie remplacer Xh par Xc
sinon réduire le simplexe par rapport à XI
Pour une description détaillée de cet algorithme voir § 3.2.2.3.
commentaires
Rappelons que
XI (respectivement
Xh,
Xg) est le point de plus petite
valeur
(respectivement
de
plus
grande
et
de
deuxième
plus
grande
valeur) de
fonction et que Xo est l'isobarycentre de
tous
les
points
(excepté Xh).
On dira qu'une réflexion est réussie si f r < fi; qu'une expansion est
réussie si fe < fi ; qu'une contraction est réussie si fc < fh .
Pour le test d'arrêt voir
§.3.2.4.
Enfin toutes les
opérations géométriques sont faites par rapport à un
seul centre (XO) et pour un seul point à la fois (souvent Xh).
C'est une des différences fondamentales entre la version séquentielle et
la version parallèle (cf.
TORCZON 1989». En effet, dans celle dernière les
opérations géométriques se font toutes par rapport à XI comme centre.
Et
de
plus, à chaque
itération, ce
sont tous
les
autres
sommets du
simplexe qui subissent ces opérations.
99
3.2.2,2. L'alll0rjthme parallélisable
Avec la version d'algorithme ci-dessus décrite, à chaque itération
la recherche se fail dans .LI.ll.'. seule direction qui est celle de (Xh,XO).
Une idée développée par (TORCZON 1989] est d'effectuer celle recherche
dans
~ les N directions. Ce qui peut se faire si on applique les
réflexions, expansions et contractions à tous les sommets du simplexe
par rapport au point XI.
Donc, dans celle version, on n'a besoin que de XI le point de plus petite
valeur de fonction.
Par ailleurs une opération sera dite réussie si elle améliore la valeur de
fonction par rapport f/=f(X/) (la plus petite).
soient a. 8 et ~ tels que
0 < 8 < 1 (coefficient de contraction)
a ~ 1 (coefficient de réflexion)
~ > 1 (coefficient d'expansion)
X2
X3
1
-
1
-
- -
/
-
1
--
1 ~ directions de recherche si -
1
f(X1) est la plus peri le valeur
1
de fonction.
1
1
1
fillure 3 2.3
Le simplexe initial et les images des sommets.
100
aJ&O[ilhme
parallélisable:
Initialisations: SO = < XOI, XOI , ..........• XOn+1 >, a, 9 , J.l >1.
Trouver
XOmin
=
Argmin
f ( Xi)
i=l •..,N+I
Intervertir
XOmin
et
XOI
k:=O;
Tant que critère d'arrêt faux
faire
(*REFLEXION*)
Faire pour i= 2•...•N+ 1
k+1
k
Xk'l
X
i = (1+ a)' X 1 - a
h l
Calculer
f ( X i )
h
k
Si
Min f ( X
1i ) i=l •...• N
< f ( X 1 )
(*EXPANSION*)
Faire pour i=2•...• N+ 1
k
k
Xk+ 1
Xei
= (I-J.l) X 1
+J.l
i
k
Calculer
f ( X ei)
i=2•...• N+ 1
Si Min f( X\\i) i=2. N+I < Min f( Xk+1i) i=2•...• N+I
k
X + 1. _ X k.
. 2
N 1
1 -
el
1= •••• ,
+
Sinon
(·CONTRACTION·)
Faire pour i=2•...• N+I
Xk+li
= (1- P)
Xkl
+ PXki
k+ 1
Calculer
f ( X
i )
i=2•...•N+ 1
h l
Trouver
X
min
A
.
f(Xk+li)
=
rgmm
i=l •..• N
k
Si f(Xk+lmin)
<
f ( X + Il )
k+
intervertir
X
1min
et
Xk + Il
k:=k+ 1;
fin (·de tant que")
101
commentaires:
On voit bien qu'il est possible de paralléliser l'algorithme ci-dessus.
En effet, la recherche se faisant dans les N directions, on peut imaginer
mettre un processeur sur chaque direction (ou sur un certain nombre de
directions).
3,2,2,3. Description de la version séQueOlielle
Notons que le but majeur de celle version est de se débarasser de Xh ,
point de plus grande valeur de fonction.
Nous
présentons
donc
dans
les
détails
les
différentes
phases
de
l'algorithme. Ainsi, on suivra mieux les différentes étapes de choix des
sommets à remplacer et les points qui les remplace, au fur et à mesure
du déroulement de
l'algorithme.
L'algorithme va de la
PHASE A à la PHASE G et tient compte des
résultats
partiels.
De ce
fait
nous avons
résumé
pour chaque
phase
toutes les possibilités: arrêt, sortie conditionnelle vers une
autre phase
ou passage à la phase suivante.
PHASE (A):
On pan du simplexe
S:
< XI, X2.
, Xn+l >
on évalue
ft,
f2,
, fuI.
PIIASE ( B):
Trouver
X h et fh
Tels que
fh:;; f(X h)
soit
la plus grande valeur des fi
X g et f g Tels que f g :;; frX g) soit la 2° plus grande valeur des fi
X 1 et fi
Tels que
fi:;; f(X 1) soit
la plus petite valeur des fi
soit donc S:
< XI,
, X, >.
PIIASE (C):
Calculer l'isobarycentre Xo de tous les points excepté Xh
soit:
Xi)
/
N.
102
PHASE (0):
On essaie
de s'Eloigner
de X h •
Pour cela on calcule X r l'image de X h
par une réflexion par rapport l
Xo où:
" X
X 0 Il
Xr=(I+a)XO - a
Xn
avec
r
..
Il X h
-
XOII
Evaluer
fr = f ( x. ).
PHASE (E):
On
compare
fr
à
fi',
SOUS-PHASE (El):
Si
fr <
fi
on a alors trouvé une
valeur de fonction plus petite.
La direction vers
X r semble être
bonne
d'où
l'idEe
d'une
expansion
dans
cette
direction.
On cakuleX e l'image de X r par une expansion par rapport à X 0
où:
X 0 Il
Xe = (1- "f) Xo + 1 X.
avec
1=11 Xe
-
" X r -
X 0 Il
Evaluer fe = f ( Xe ).
On
compare
fe
à
fi.
X2
Xh
Xr
-- - - -- Xe-.
XI
figure 3.2.4, Le simplexe et les images après réflexion et expansion.
103
SOUS-PHASE (KU):
Si
fe
<
fi
remplacer
X h par Xe
on a ainsi un nouveau simplexe:
S :
< Xl,
, X g, Xe >
s.atisfait
alors
ARRET
Test
d'Arrêt
1 SI non
aller à
(B).
SOUS_PHASE (E.1.2):
Si
fe
1
fi
On abandonne Xe: on s'est certainement éloigné dans une
mauvaise direction.
Néanmoins remplacer X h par X r ( rappelions qu'on a fr < fi)
Le nouveau simplexe est donc:
S =< Xl, ......., Xg , Xr >
satisfait
alors
ARRET
Test
d'Arrêt
1 si non
aller à
(B).
SOUS-PHASE (E2):
Si
fr
1
fi
mais
Si
fr
i
f g
(on a donc une amélioration par rapport aux mauvais points X h et X g)
on remplace
X h
par
X r
et le simplexe devient:
S = < Xl, ......., Xg , Xr >
s.atisrait
alors
ARRET
Test
d'Arrêt
1 SInon
aller à
(B).
SOUS-PHASE (E3):
Si
fr 1
fi
et
Si
fr >
f g
(pas d'amélioration par rapport aux mauvais points ... )
procéder à la
phase
(F).
104
PHASE ( F ):
On compare
fr
à
fh.
SOUS-PHASE (FI) :
a)
Si fr
>
fh
on procède à (F2) (i.e. une contraction).
b)
Si fr
,(
fh
on remplace
X h
par
X r
et on a le simplexe suivant
S =< Xi,
, Xg , Xr >
qu'on contracte (ensuite) en
procédant à la phase (F2).
SOUS-PHASE (F2): CONTRACTION DU SIMPLEXE
a) Si fr
>
fh
faire une contraction bis (i.e. sur le segment
[XO,Xh) pour trouver Xc où:
Xc=(I-P)XO +pXh
avec 0<1} <1
b) Si f r
> fh faire une contraction simple (i.e. sur le segment
[X O.X r] rappelons que (FI b) a été réalisé i.e. X h a été remplacé
par Xr) pour trouver Xc où:
Xc=(I-P)XO+pXr
avec 0<1} <1.
X2
Xh
Xc simple
Xr
--._- --- - - Xe-.
XI
fleure 3.2.5. Les images après réflexion, expansion et contraction.
105
PHASE ( 0 ): Evaluer fe = f ( Xc ) puis on compare fe à fh.
SOUS-PHASE (01): Si
fc
< fh
on remplace X h par X c et le simplexe devient
S :
< Xl,
, X g • Xc >.
SOUS-PHASE (02) : Si fc
l
fh
on a un échec total. On est pas
parvenu à se débarasser de
fh. On procède alors à une réduction
de la taille du simplexe:
réduction de
la taille du
simplexe
on divise par 2 la distance de chaque point à X 1; où f(X 1) est la
plus petite valeur de fonction).
On a ainsi.
Xi = 0.5· (Xi + XI)
i =1•..• N+l
i .. 1
et le simplexe réduit est
S:
< XI, 0.5· ( Xi + XI) ,....,0.5 • ( Xs + XI) >.
Evaluer fI •.....• fn+ 1
les valeurs de fonction aux sommets du nouveau simplexe;
Aller au
(0).
X2
Xh
Xl
fi&ure 3.2.6. Les sommets du simplexe après contraction.
106
Xh
X2
XI
fieure
3.2.7. Formes possibles du
simplexe et ses
éventuels sommets
après une
itération.
Le
sommet
Xh est dans tous les cas remplacé. X, n'est jamais
remplacé.
X2 n'est remplacé que quand il y a réduction du simplexe;
c'est notamment le cas où on ne parvient pas à améliorer les valeurs de
fonction par rapport aux mauvais points Xh et X g'
En partant du simplexe initial S= < X"
X2, Xh > après une itération il
devient l'un des simplexes suivants:
<X"X 2,Xe> ou <X"X2,Xr> ou <X"X2,Xc> ou <X"X2,X'c>
ou enfin < X" X'2' X'3 >.
J07
3.2 3.
Ah:orilhme à recherche adaptée
X2
- - ---
- - - --"ç
dlrection de
recherche pour l'algorithme
séquemiel.
XI
X2
X3
.... ....
1
.... .... ....
1
.... ....
/
....
1
1 .... directions de recherche si
1
f(X1) eSI la plus petlte valeur
1
de fonction pour l'algorithme parallèle,
1
1
1
fieure 3,2,8.
Nombre de directions de recherche dans R2.
108
Dans les deux cas de l'algorithme séquentiel ou parallèle le nombre
de directions est fixé à l'avance, Et ce nombre ne tient pas compte de
l'état d'avancement du processus d'optimisation, En effet il ne prend en
compte ni le nombre d'échecs sur une direction donnée ni la nature du
succès
(réflexion
ou
expansion
ou
contraction
réussie).
L'expérience
montre qu'en début d'exécution une contraction est
peu
"bienfaisante"
car elle réduit rapidement le simplexe. De même, en
fin d'exécution,
avoir
une
seule
direction
de
recherche
n'est
pas
approprié
pour
échapper à un état de minimum local.
C'est pour ces raisons que nous avons rajouté à l'algorithme séquentiel
une
procédure
Accept de détermination
du
nombre
de
directions
adaplées à chaque étape de la minimisation. Cela donne l'algorithme
suivant:
aleorithme à recherche adaptée
Initialiser SO = < X l, X2, ,Xn+1> , a, p, y
Calculer fi, f2, ... , fn+1
Tant
que critère d'arrêt faux faire
Trouver XI, Xg, Xh, fi, fg, fh.
Déterminer
Xo (isobarycentre des sommets excepté Xh)
fa i re une réflexion de Xh par rapport à X 0
si la reflexion est réussie
alors faire
une expansion
si l'expansion est réussie remplacer Xh par Xe
sinon remplacer Xh par X r
sinon
faire appel à Accept
Fin
La procédure Accept détermine, en fonction du nombre d'itérations, si
une contraction est opportune ou pas. Si la contraction est justifiée, on la
fait. Si elle est réussie tout se passe comme dans l'algorithme séquentiel;
sinon on détermine l'opportunité d'une réduction du simplexe.
Une contraction est jugée opportune si on ne trouve pas de succès sur
aucune des
autres directions supplémentaires (qu'on aura déterminées).
Il en sera de même pour la réduction du simplexe.
109
Pour le calcul du
nombre des
directions supplémentaires nous avons
uulisé deux possibilités: déterministe el aléatoire,
En plus de la direction principale qui est celle de (Xo Xh) nous retenons
un
certain
nombre,
"n_directions", de
directions
supplémentaires
sur
toutes les N directions possibles. On trouve ces N directions possibles en
déterminant les directions joignanl un sommet Xs el l'isobarycentre de
tous les autres sommets
excepté Xs. Ce sommet Xs est d'abord Xb '
ensuite Xg le sommet de deuxième plus grande valeur de fonction, puis
X g.l
le sommet de troisième plus grande valeur de fonction, puis le
sommet de quatrième plus grande valeur de fonction, et ainsi de suite...
el alors la direction est celle joignant ce point et l'isobarycenlre.
Le nombre (minimum un) de directions supplémentaires retenues
sera
"n_directions"
avec:
n fi x é
choix d é t e r m i n i s t e
n_directions =
l sup (l, n_aléatoire) choix aléatoire
n_aléaloire un nombre entier pris au hasard entre 0 el N-1.
X2
1
1
1
- -
---
1
+
direction de recherche
1
principale.
1
1
1
1
XI
1
1
1\\
:
direction de recherche supplémentaire.
fi ~y[e 3,2,9. Nombre de directions pour l'algorithme à recherche adaptée,
110
3.2.4.
Convergence
de
la
méthode
du
Polytope
Mouvant
Cette
méthode,
comme
plusieurs
méthodes
directes,
n'est
pas
basée sur
un théorème de convergence théorique établi pour tous les
types de
problèmes pour
lesquels elle est
sollicitée. L'article original
[NEL65) ne propose aucun théorème de convergence.
Le travail de [DENNIS & TORCZON
1990) (sur l'algorithme parallèle)
propose
quant
à
lui,
un
théorème
de
convergence
mais
pour
les
fonctions continument différentiables. Or, nous savons que pour de telles
fonctions
il y a des
algorithmes d'optimisation
locale
beaucoup plus
performants (car l'algorithme parallèle est un algorithme d'optimisation
locale).
Enfin
retenons
que, à ce jour, il
n'y
a pas
de
théorème de
convergence théorique établi pour le type de problèmes où l'utilisation
de cette méthode est courante et justifiée: en optimisation non linéaire
non
différentiable.
De ce fait, la convergence de la méthode est testée numériquement
par
un
test d'arrêt que
l'on
suppose
suffisant
pour
assurer que
le
simplexe final est réduit (en taille) autour d'un minimum.
Différents tests sont proposés. Nous en citons les deux suivants:
• Test de INELDER & MEAD 1965)
C'est celui présenté dans l'article original.
Le test d'arrêt est basé sur la déviation standard des
N+ 1 valeurs de
fonction aux sommets su simplexe. Si la déviation est inférieure à une
certaine valeur E 2 donnée, le simplexe donnerait alors une solution.
C'est-à-dire :
~I
si
L
S
E 2
N+I
I, f(Xi)
N+I
j-I
1 1 1
• Test de (PARKINSON & HUTCHINSON 1972]
Il est basé à la fois sur les variations d'amplitude de la fonction et sur la
proximité des sommets du simplexe.
s El
Le test est satisfait si
N+l
1
L
et
N+I
Il X k + 1 i - X k , 11 2
< E2
i= 1
où El et E2
sont deux seuils donnés, et Xki est l'image du sommet Xi à
l'itération k.
• Test de (WOODS 1985)
Il est basé sur la mesure des variations de la taille du simplexe au cours
des itérétions.
Le test est satisfait si
Max { 1- Il Xk+1i - Xkl Il,2::;; i s N+t }::;; E;
où A = Max ( l , IIXktll 1 et Xki est l'image du sommet Xi à l'itération k.
commentaires
Evidemment,
le
test
d'arrêt
choisi
pour
J'implémentation
numérique,
dans
ce
cas
aussi,
affecte
la
précision
des
résultats,
le
nombre d'évaluations de la fonction, et dans certains cas, la convergence
de la version séquentielle.
112
InirilliuDon!:
Simple.. ini,i.' < XI. Xl •• Xn>.>
PlrllrnttreS a. Il. T·
oui
non
Te"d'WI
non
oui
113
3.2.5.
Résultats
Numériques
Rappelons que nous nous sommes intéressés à celle méthode parce
qu'elle est directe et n'utilise pas de gradients de la fonction objectif. Et
c'est
pour cela
que
nous
avons envisagé
l'utiliser
pour
traiter
des
problèmes pratiques . Pour ce faire nous avions plusieurs possibilités:
• l'utiliser comme algorithme local directement;
• l'utiliser dans l'algorithme du Tunnelier (cf. 3.3), qui lui est un
algorithme global, comme sous-programme local;
• ou modifier l'algorithme de façon qu'il soit un "peu plus global".
C'est ainsi que nous avons implémenté en Pascal sur P.C. et sur
C.D.C.
les
deux
versions (séquentielle et
parallélisable) et
les
avons
expérimentés sur plusieurs fonctions tests. Les
résultats ainsi obtenus
nous ont
amené à modifier l'algorithme séquentiel pour tenir compte
des avantages et des inconvénients des deux versions.
En effet:
•
la
version
séquentielle (recherche
monodirectionnelle) échappe
parfois à certains minimums locaux, ce qui
n'est pas
le cas de l'aune
version.
• la version parallélisable (recherche multidirectionnelle) converge
toujours (échec très rare) vers un minimum (local cependanrl).
D'où
l'idée
de
l'algorithme
à
recherche
adaptée
que
nous
avons
implémenté. Cet algorithme, non
plus, n'a
pas encore de justification
théorique. Son succès dépend de beaucoup des petits coups d'ajustement
que nous y avons rajoutés. En tous les cas, dans certaines situations, il
converge
vers
un
minimum
global.
Et
ce
sont
les
résultats
de
cet
algorithme que nous avons utilisés pour une comparaison avec le Recuit
Simulé (cf. 3.2.).
1 14
Résultats de minimisation elobale de J'aleorithme (du Polytope MouvanJ)
à recherche adaptée comparés à ceux du Recuit Simulé sans cycle et
avec cycles (cL § 3.1.);
Sur une moyenne de 10 exécutions:
Simplexe
Recuit
1 cycle
Recuit
plusieurs cycles
fOI
.(9)L(I )
.(10) en 3 cycles
f02
.(8)L(2)
.(10) en 4 cycles
f07
.(9) en 3 cycles • ( 1)
flO
.(2)L(8)
.(3)L(7)
.( 10) en 7 cycles
Dans chaque cas;
a signifie que le minimum global est atteint; le nombre de fois est entre
parenthèses
L signifie qu'un minimum local est atteint
• signifie que le code n'a pas convergé.
3.2.6
Conclusion
La méthode du Polytope Mouvant est une méthode directe facile •
implémenter.
Dans
le
cas
d'un
problème
d'optimisation
locale
non
différentiable, ou dans les cas où la fonction objectif est sous la forme
d'une "boite noire", elle peut être utile car peu exigente. Cependant la
convergence n'est pas toujours garantie.
Le
succès des
algorithmes
basés sur
cette
méthode dépendent
beaucoup des choix du simplexe initial, des coefficients (de reflexion,
d'expansion et de contraction) et des tests d'arrêt.
Choix des paramètres a. Il. et Q;
Les
valeurs
des
paramètres "géométriques" avec
lesquels
nous
avons enregistré les meilleurs résultats sont les suivants;
• coefficient de réflexion
a = "
• coefficient d'expansion
Il = 2,
• coefficient de contraction
Q = 0,5.
1 15
Choix. du sÎmplexe Înitial:
Le
simplexe initial doit être impérativement lei
que
les droites
joignant le sommet X1 aux autres sommets engendrent l'espace Rn. Pour
ce faire nous avons procédé de la façon suivante:
Pour Xo un point initial donné, on construit le simplexe initial en
faisant
XI :;;Xo
Xi:;; XI + k ei-I
pour tout i=2, , N
où:
ei est le i-èrne vecteur de la base cononique (les N directions engendrant
ainsi Rn);
k un pas arbitraire.
Pour le test d'arrêt on ne peut pas désigner de test optimal car cela
dépend
essentietlernent
du
type
de
fonction
objectif:
fonction
avec
variations d'amplitude importantes, minimiseurs très
serrés, etc ...
Enfin, même si nous avons étudié de
près cette méthode, nous
n'avons pas pu l'utiliser pour les applications du C.N.E.S car le type de
problème traité au chapitre 4 n'est pas compatible avec les opérations
géométriques sur
l'espace de
variables correspondant (espace composé
de dates de survol et vitesses de survol).
rHirences
[DEN87) J. E. DENNIS, Jr & D. J. WOODS, Optimizatian an microcomputers:
the Ne/der & Mead simplex algarithm, in Arthur Wouk, editor, New
computing environments: Microcomputers in Large-scale Computing,
pages 116-122, SIAM, Philadelphia (1987).
(DEN90) J. E. DENNIS, Jr & V. TORCZON, Mu/ti-Directiona/ Search: A new
seerch
algarithm, tech. rep., Dept. of Mathematical Sciences, Rice
University, Houston, Texas 77251-1892, 1990.
(NEL65]
J. A. NELDER & R. MEAD, A slmplex method far tunction
minimizatian, Computer Journal 7, (1965), p. 308-313.
116
[PAR72] J. M. PARKINSON & D. HUTCHINSON, An investigation into the
efficiency of variants on the simplex method, in F. A. Lootsma, editor,
Numerical Methods for Non-Iinear Optimization, p. 115-135, Academic
Press, London & New York, (1972).
(SPE62) W. SPENDLEY, G. R. HEXT & F. R. HIMSWORTH, Sequential
applications of simplex designs in optimization and
evolutionary
operation, Technometrics 4, (1962), p. 441- 461.
(TORC89] V. TORCZON, MurU-Directional Search: A Direct Search
Algorithm for ParaI/el Machines, Ph. D. thesis, Rice University, Houston,
Texas, (1989).
(W0085) D. J. WOODS, An Interactive Approach for Solving Mufti-
Objective Optimization Problems, Ph. D. thesis, Rice University, (1985).
117
3.3.
La
Méthode du
Tunnelier
Résumé:
La
méthode
du
Tunnelier
(Tunneling
Method)
est
une
méthode
déterministe
qui
se
propose
de
trouver
le
minimum
global
d'une
fonction f sans contraintes en utilisant des codes d'optimisation locale.
Son principe, simple, consiste à chercher un minimiseur local et de s'en
servir pour éliminer tous les points ayant une valeur de fonction plus
grande ou égale à celle de ce minimiseur. Cela grâce à une
fonction
auxiliaire associée à f qui permet aussi de déterminer un nouveau point
initial à partir duquel on refait une minimisation locale. Chaque nouvelle
minimisation
donne
de
meilleurs
résultats
car,
auparavant
on
s'est
débarassé des points de valeur plus grande ou égale à celle du dernier
minimiseur.
L'algorithme s'arrête quand,
pour le
dernier
rninimiseur,
tout se passe comme si on avait éliminé tous
les points candidats à
donner des valeurs de fonction plus petites!
3.3.1.
Principe
Soit une fonction f: Rn ~ R
à minimiser.
La
méthode
du
Tunnelier
(Tunneling
method)
est
une
méthode
introduite
par
[MONTAL VO
19791.
généralisation
d'une
méthode
introduite par [VILKOV et ail en 1975 au cas multidimensionnel.
Le principe
consiste à engendrer une suite de minimums locaux à
valeurs de fonction décroissantes. Le dernier minimum est le minimum
global.
118
On
peut considérer le Tunneling comme une
généralisation de
la
technique dite de "technique de déflation" de [GOLDSTEIN & PRICE 1971]
pour trouver le minimum global d'un polynôme (de dimension 1)
soit f un polynôme dont on cherche un minimum global;
si X· est un minimiseur local de f
on définit
ainsi si f est de degré m alors fI sera de degré m-2;
Et si on parvient à montrer que fl(x) n'est pas négative alors X· est
un minimiseur global de f. Par contre s'il existe un XO tel que
f(XO) < f(X·),
X· n'est pas un minimiseur global, et en ce cas on pourrait repartir
de XO comme point initial pour une nouvelle minimisation locale... et
ainsi de suite. Tant que fk n'est pas positive on cherche un point
initial XO tel que
f(XO) < f(X·k)
où X·k est le dernier minimiseur trouvé, et on réoptimise,
3,3.2.
Algorithme
du
Tunnelier
La
méthode
du
Tunnelier,
qui
s'inspire
de
la
méthode
de
[GOLDSTEIN & PRICE 1971] , est constituée de deux phases principales:
phase 1 : recherche d'un minmuseur local
à partir d'un point initial Xo donné on utilise un algorithme de
minimisation locale (mais de type descente) pour trouver un minimiseur
local X·;
phase2 : phase du Tunnelier (raffinement):
le but de celle phase est de trouver un point XO • x·
tel que
f(XO) 1 f(X·) ; ce point est alors utilisé comme point de dépara pour une
nouvelle phase de
minimisation (phasel
précédente).
119
Pour trouver un X" tel que f(X") { f(X·) on cherche plutôt un X tel
que f(X) = f(X·) ; pour ce faire il faut trouver un zéro de la fonction T
dite
du Tunnelier ainsi définie :
f(X) - f(X·)
Tx·(X)=
r
IIX - Xmll~ 0 TI IIX - Xj·lI~ j
j" )
où X·). X· 2 .
. X· r sont les derniers minimiseurs locaux trouvés.
• Propriétés de la fonction
du Tunnelier T:
• le terme f(X)·f(X·) assure l'élimination des points tels que
f(X) > f(X·), qui ne pourraient être des minimiseurs globaux.
• le terme IIX - Xi·lI~ i permet d'éviter de choisir encore les points
X·j
ou tout autre minimiseur ayant pour valeur f·=f(X ·i) i= l , ....r.
• le terme IIX - XmII~ 0 évite à l'algorithme ( de résolution de T(X)=O ) de
f(X)-f(X·)
converger
vers
un
point stationnaire de
IIX _ Xi.lI~ i
qui
n'est pas
solution de
T(X)=O. En d'autres termes, les constantes Xm et ~ 0 sont
choisies de façon à éjecter tous les minimiseurs locaux inintéressants de
T
(points
où
T'(x)=O
et
T(x»O)
risquant
dattirer
une
routine
de
minimisation destinée à trouver une racine de T(x)=O.
décroissance des valeurs des minimums locaux
on a la propriété suivante:
f(X·,) ~ f(X· 2) ~ f(X· 3) ~ ... ~ f(X·g_) ~ f(X· g)
où
X·). X· 2 '
. X· g sont les minimiseurs locaux trouvés à l'étape
1.2.... et X· g un minimiseur global.
Donc à échape du Tunnelier on améliore le dernier minimum trouvé.
• condition
d'arrêt:
elle est basée sur la condition suivante:
Si on trouve X· un minimiseur tel que
V X E S Tx. ( X ) > 0,
alors
X· est minimiseur global.
120
y
x
y
y =T (X)
x
fi&urc 3.3.1
Une première optimisation locale donne X*, à partir duquel on trouve XO
(un zéro de T) ; en réoptimisant en partant de XO comme point initial on
trouvera un minimimseur de valeur plus petite, etc.
121
al&orithme &énéral du Tunnelier
Initialiser XO, k:= 2;
Trouver X·l minimum local de f
Tant
que critère d'arrêt faux fa ire
trouver les paramètres de la fonction T(X·k _, ) du Tunnelier
SI T est négative ou nulle trouver XO racine de T
trouver un minimum X·k de f en partant de XO
rin (si)
k:=k+l;
rin (tant que)
A l'arrêt X·k est un candidat minimiseur global.
3.3.3.
Choix
des
paramètres
Nous avons utilisé le choix donné dans ILEVY & MONTALVO 1980b).
Dans ce
cas,
la
recherche d'une solution de
l'équation T(X)=O est
remplacée par la recherche d'un
XO tel que T(XO) < E donné. C'est-à-dire
qu'on applique une routine de minimisation l T(X) jusqu'à obtenir un tel
XO. Remarquons que dans ces cas un minimiseur local de T qui a une
valeur de fonction supérieure à E n'est pas intéressant d'où la nécessité
du terme IIX - XmllÀ 0 au dénominateur.
Choix
des
paramètres
de
la
fonction
du
Tunnelier
f(X) - f(X·)
Tx·(X)=
r
IIX - Xm"À 0 fi IIX - X
À
j·II
j
j", t
Détermination de À~•••
Au premier passage à la phase du Tunnelier (i.e. quand l'algorithme
a trouvé X ••) on pose À 0=0 et À. '" b.
où b. est calculé de la façon
suivante:
soit
X=X.· + r
un vecteur aléatoire de norme (très) inférieur à 1; la
valeur de
b. est trouvée
itérativement en
partant de
b. =1 que
l'on
augmente jusqu'à ce qu'on ait la propriété rt.dx > 0
et la condition
d'optimalité du
premier ordre T'(x)t'dx < 0; où dx est la direction de
descente (de la roui ne qui minimise T) calculée pour b. = 1.
122
A
tous les
i-ème passages suivants à la phase du Tunnelier (on
dispose
donc
de
AI.A 2•...•Ai . I). on détermine Ajen utilisant une
heuristique:
pour
chaque
itéré
X.
on
calcule
si
les
termes
du
dénominateur
doivent
intervenir
ou
pas
en
faisant
une
sorte
de
"topographie locale" autour de X; ainsi si X est jugé proche de X* j on
pénalise par les dénominateurs pour l'éloigner de ces points indésirabes;
sinon on ne pénalise pas et on prendra A j=O. Dans le cas de pénalisation
A j est déterminé de la façon suivante:
o
si
1 X- X* d:i! 1 + C 2'
b j
si
1 X· X* j 1< 1 - C 2
(bj/2) (1 +( 1-
1X- x-, 1),C2)
sinon
avec C 2 constante (petite de l'ordre de 1. E-5).
Détermination de Ao. et de Xm
Ces termes ont pour rôle d'empêcher à la routine de minimisation de T
de
converger vers
des points stationnaires de T ayant une valeur de
fonction supérieure au seuiJ E choisi.
Pour cela si on ne trouve pas XO (que J'on cherche tel que f(XO) < E). on
regarde si la variation entre l'itéré courant Xk et Je précédent Xc (=Xk-I)
est faible ou si le gradient de T est très petit; dans ces conditions Xc
semble être un point stationnaire "indésirable" et donc on pénalise grâce
à Ao et Xm que l'on choisit comme suit.
X
= 1 Xc
si
1 Xc- Xkl<l
m
u-Xc + (I-u) Xk
sinon
où 0 < u s 1 est choisi. quand même. tel que
1x-X ml-cl;
on choisit ensuite AO pour éloigner l'algorithme de Xc. On commence par
une certaine valeur initiale que J'on incrémente par la suite jusqu'à ce
qu'une
heuristique
détermine
si
on
a
quiné
le
voisinage
du
point
indésirable.
Remarque
Ce choix de paramètres est sans doute très complexe car il nécessite
l'utilisation et la maîtrise de beaucoup d'heuristiques.
123
3.3.4.
Convergence
de
la
méthode
du
Tunnelier
C'est
une
méthode déterministe dont la propriété de
convergence
fait
partie de
sa
philosophie même. En effet sa convergence est en
théorie garantie car tant que la fonction auxiliaire T n'est pas positive il
est possible (théoriquement) de trouver un zéro de T et de réoptimiser.
Donc la seule difficulté (et
pas la moindre) est, dans la pratique, de
trouver un zéro de T ou un minimiseur X" de T tel que
f(XO) < E.
3.3.5.
Résultats
Numériques
La
Méthode du Tunnelier est une méthode -dont le principe
très simple à comprendre- lui donne assez de "charme".
En effet, il suffit de disposer d'un bon code d'optimisation locale pour
éliminer au fur et à mesure tous les minimums locaux; le dernier trouvé
est candidat minimum global.
Cela est vrai tout au moins en théorie. car en pratique pour passer d'un
minimum
local
Xi· au suivant il faut trouver XO tel que f(XO) = f(X i·)
c'est à dire il faudrait résoudre une équation du type T(X)=O.
Or,
justement. on sait que
parfois. il est
nettement plus difficile de
résoudre certaines équations que
de
les
ramener sous
la
forme
d'un
problème d'optimisation (cf. 2.7).
En somme, on remplace un problème difficile (l'optimisation
globale) par un problème parfois plus difficile (résolution d'équations)
qui. lui même. se ramène à la résolution d'un problème d'optimisation ...
Nous nous
sommes intéressés à cette méthode dans
le
but
d'en
implémenter un code pour traiter des problèmes pratiques.
Nos premières expériences ont été catastrophiques car nous n'avons eu
aucune convergence pour les
problèmes avec
plusieurs minimums. En
effet. le code C05ADF. de résolution de T = 0 que nous avons tiré de la
bibliothèque
NAG.
a toujours été
tenu
en
échec.
Alors
nous
avons
remplacé la résolution par un autre (!) problème d'optimisation de T(x)
(cf. 2.7.1) et avons tiré les conséquences suivantes:
• il est très difficile de trouver les paramètres de
la fonction du
Tunnelier
sans
qu'on
observe
des
divergences
de
l'algorithme
de
résolution d'équation dans les voisinages de certains points,
• quand les minimiseurs sont très proches, par exemple la fonction
test
FI2 (cf,3.1.5), la
phase du
Tunnelier n'est pratiquement
utilisée
qu'une seule fois. Et donc si le deuxième minimiseur n'est pas un global
il n'y a plus aucune chance de trouver un global.
124
• le choix pour la fonction du Tunnelier comme une fraction conduit
l des divergences. D'ailleurs des auteurs de travaux sur l'algorithme du
Tunnelier se penchent maintenant vers
une nouvelle fonction auxiliaire
du type T(x) = If(x) - f*] • EXP( h (X*I» où h est une fonction jouant le
même rôle que le dénominateur dans la fonction T classique.
Enfin, nous avons abandonné celle méthode pour les difficultés de son
implémentation. De plus, elle non plus, tout comme le polytope mouvant,
n'est pas adaptée au genre de problème que
nous avions en
vue (cf.
chapitre 4).
relérences
(GOL71] A.A GOLDSTEIN & J. F. PRICE, On Descent from local minima,
Mathematics of Computation 25, p. 569-574, (1971).
(LEV80al A.V. LEVY & S. GOMEZ, The Tunneling gradlent-restoration
algorithm for the global minimization of nonlinear [unctions subject to
nonlinear
inequality
constraints,
liMAS-UN AM,
Communicaciones
téchnicas, Serie Naranja: Investigaclones, N° 231 J. Sci. Star. Comp. 6
( 1980)
ILEV80b) A.V. LEVY & A. MONTALVO, A modification to the Tunneling
algorithm for the global minimization of functions, Tech. Report 240,
Univ. Nat. Auron, de Mexico (1980)
ILEV851 A.V. LEVY & A. MONTALVO, The Tunneling algorithm for the
global minimization of functions, SIAM J. Sei. Stat, Comp. 6, p. 15-29,
(1985)
(MON79] A. MONTALVO, Development of a new algorithm for the global
minimizatlon
of functions, Ph.
D. Thesis in Theorical and
Applied
Mechanics, Universidad Nacional Autonorna de Mexico (1979).
(VIK751 A. V. VIKLOV, N.P. ZHIDKOV & D.M. SHCHEDRIN, A method of
search for
the
global minimum of a function
of one
variable, J. of
Compul. Math.& Math. Phys. 15, p. 1040-1042, (1975).
Ils
Chapitre
4:
OPTIMISATION
DE
TRAJECTOIRES INTERPLANETAIRES
126
4.0
Introduction
Le travail de cette thèse ayant été proposé sur le thème "étude de
problèmes
d'optimisation
globale",
il
a
été
divisé
en
trois
phases
chronologiques:
la
première, consistait
11
faire
l'état
de
l'art
sur
les
méthodes et codes existants sur le sujet, la deuxième, 11 programmer des
codes en
vue
de
les
appliquer, en
troisième
phase,
11
des
problèmes
spéclfiques
du
C.N.E.S.
Les
résultats
des
deux
premières
phases
correspond-nt
11
ceux
présentés
aux
trois
chapitres
précédents.
La
dernière
phase,
quant
11
elle,
a
consisté
11
expérimenter
le
code
RECSIMU.FOR (cf. [BONNEMOY & HAMMA 1991a)) sur des problèmes
d'optimisation
de
trajectoires
interplanétaires.
Le
code
RECSIMU.FOR
basé sur le Recuit Simulé utilise des cycles de
réchauffement en
plus,
pour renforcer l'algorithme standard du Recuit, dans les cas difficiles.
L'optimisation
des
trajectoires
qui
nous
a été
proposée,
provient
du
constat suivant:
- dans le cadre de la mission VESTA, il est nécessaire d'optimiser les
impulsions utiles 11 la réalisation de cette mission,
- une optimisation locale a révélé l'existence de plusieurs minimums.
Donc, on est bien confronté 11 un problème d'optimisation globale.
Ainxi,
dans
un
premier
temps,
nous
avons
continué
d'optimiser
localement
les
trajectoires
pour mieux cerner le
problème (car on
ne
disposait pas de
fonction
objectif exprimée directement en
fonction de
ses variables), et choisir une
méthode d'optimisation globale adéquate.
Dans un second temps, nous avons programmé le code le RECSIMU.rOR.
Nous avons testé ce code sur des exemples tirés de
la littérature, puis
l'avons appliqué au cas du C.N.E.S.
Malheureusement,
11
la
fin
du
temps
contractuel,
nous
n'avions
pas
encore
fini
d'harmoniser
les
différents
sous-programmes
propres
au
C.N.E.S. (pour le calcul de la fonction objectif) avec le code RECSIMU.rOR.
Des propositions utiles ont été données dans (HA MMA91),
Dans ce qui suit. nous
présenterons en détail des résultats donnés par
l'optimisation
locale, puis nous parlerons des difficultés rencontées par
le
code d'optimisation globale.
127
4.1.
Présentation
de
la
mission
VESTA
Le
problème d'optimisation
auquel
nous nous sommes
intéressés
concerne I'cpumisation de trajectoires interplanétaires dans le cadre de
la mission VESTA qui poursuit deux objectifs majeurs (horizon 1994):
1) Une exploration de MARS
2)
Une exploration de ce que l'on
appelle les
"petits corps" ou
astéroïdes.
Le projet de la mission VESTA est réalisé en collaboration avec l'Union
Sov iétique.
Il est prévu de lancer une sonde composée de deux parties distinctes:
l'une de ces parties
se satellisera autour de la "planète rouge"; l'autre se
chargera de survoler le maximum d'astéroïdes possible pour une durée
du
voyage
ne dépassant pas
5 ans.
Les
collaborateurs de
l'U.R.S.S.
s'occuperont de l'étude de MARS et fourniront le lanceur. Au cours de
celle mission VESTA de
nombreuses expériences sont prévues, parmi
lesquelles on peut noter:
différentes
mesures
spectrométriques:
ultra-violet,
masse
neutre,
infra-rouge,
des pénétrateurs qui seraient largués sur les plus gros astéroïdes
visités,
une expérience de détermination de masse ponctuelle lâchée de
la sonde.
Outre la construction de la sonde et la mise au
point des différentes
expériences prévues, une des grandes phases du projet est le choix et le
calcul des
trajectoires adéquates.
Il faut, en effet, tenir compte de nombreux critères, tels que les temps
de
transfert entre deux corps à visiter (qui ne doivent pas être trop
longs) ou la taille de ces corps (plus ils sont gros, mieux c'est). N'oublions
pas que pour permettre une étude quelconque, les
vitesses relatives de
survol des corps ne doivent pas être trop élevés. La sonde ne peut, de
plus, emporter qu'une quantité limitée d'ergol.
On le voit, le choix n'est pas simple, d'autant que le nombre d'astéroïdes
recensés aujourd'hui s'élève à Il 449.
D'où la nécessité d'optimiser!
L'optimisation se situe justement au niveau de la "consommation". En
effet, afin de
recueillir le plus de
renseignements possible, on
essaie
d'emporter le maximum de matériel scientifique. Pour gagner du poids,
on minimise la quantité de carburant à embarquer, et on s'efforce d'en
consommer le moins possible pour pouvoir aller le plus loin possible.
C'est dans cette oprique qu'on est amené à utiliser l'effet de Swing-by ou
effet de "tremplin gravitationnel" des planètes [cf CNES Il; ce qui permet
de
modifier
la
direction
de
la
sonde
sans effectuer de
manœuvres,
coûteuses en
carburant.
128
4.2
Notations
: corps célestes visités
(Mi)i = I,N-I : corps fictifs
: dates de survol des corps réels
dates des
manoeuvres intermédiaires
position et vitesse du corps Pi à la date Ti
v.d
vitesse absolue de la sonde au départ de Pi
1
v.a
vitesse absolue de la sonde à l'arrivée de Pi
1
v.d
vitesse relative de la sonde au départ de Pi
1
v.a
vitesse relative de la sonde à l'arrivée de Pi
1
vitesse absolue de la sonde au départ de Mi
vitesse absolue de la sonde à l'arrivée de Mi
vitesse infinie maximale de départ de Pi
angle de déviation maximale du cône de swing-by
angle de déviation nécessaire
constante d'attraction gravitationnelle de Pi
incrément réalisé lors du survol de corps réel ou fictif
rayon de péricentre de survol de Pi
rayon de péricentre minimum de survol de Pi.
129
4.3. Position
du
problème
Le problème à résoudre est le suivant:
Elanl donnée
une
trajectoire
de
sonde
interplanétaire
visitant
les
N
corps célestes (Pi)i = l,N aux dates Ti
i=l,N,
optimiser celle trajectoire sachant que:
• les
manœuvres son! impulsionnelles,
• les manœuvres sont réalisées lors de chaque survol, avec
possibilité d'en rajouter entre deux survols,
• la fonction COÛI est égale à la somme des modules des /1v
résultant de ces manœuvres.
Les contraintes à prendre en compte sont:
• la vuesse de dé pari du premier corps ne peul pas dépasser une
valeur limite fixée (hmitarion du lanceur).
• lors du survol d'une planète, il existe une altitude minimale de
survol à ne pas dépasser: rim'
Enfin, suivant que l'on mel ou non des manœuvres entre les survols, on
obtient
deux
sortes de
problèmes,
assez
similaires,
mais
néanmoins
distincts.
4.4 Résolution
du
problème
4.4.1
Méthode
de
résulurlen
Dans les deux cas, la méthode de résolu lion choisie est celle dite des
"patched-conics".
On
considère
que
la
sonde
n'esl
soumise
qu'à
l'attraction d'un seul corps qui est:
- soit le soleil ( pratiquement tout le lemps),
- soit une planète (aux
abords de celle-ci).
La
trajectoire
est
alors
une
succession
d'arcs
d'ellipses
héliocenlriques
(I.e.
dont
le
soleil
est
un
foyer)
el
d'hyperboles
pl ané tocentriq ues.
On a alors à résoudre un nombre fini de problèmes aux deux corps.
130
4.4.2
Formulation
La formulation du problème est un peu différente suivant que l'on
considère un problème avec ou sans sans manœuvres intermédiaires.
4.4.2.1
Cas
sans
manœuvres
Intermêdlalr es
C'est le cas le plus simple. On veut trouver la trajectoire passant par
les N corps
(Pi)i = I.N aux dates Ti
i= I.N
en faisant des manœuvres
impulsionnelles au survol de ces corps.
On doit résoudre une suite de problèmes de LAMBERT de la forme:
" Trouver le transfert de Lambert allant de Pi à la date Ti jusqu'à Pi+ 1 à
la date Ti+ l ".
On a donc N-I problèmes de Lambert li résoudre.
Les variables du problème sont alors les N dates de survol ( Ti)
i= I.N .
Et le coût a pour expression :
N -1
J . IIIL1vi Il
i=1
4.4.2.2
Cas
avec
manœuvres
interm~dialres
On
rajoute
maintenant
la
possibilité d'effectuer
des
manœuvres
entre
deux
survols.
On
peut
supposer
que
ces
manœuvres
intermédiaires sont réalisées lors du survol des corps fictifs (Mi) i=I.N-1
d'attraction négligeable. On résoud alors la suite de problèmes suivants:
" Extrapoler la trajectoire du corps réel Pi à la date Ti jusqu'au corps
fictif
Mi
à la date Ti' connaissant la vitesse relative V d
i
de la sonde au
départ de
Pi;
ceci
permet de
calculer la
position de
la
manœuvre
intermédiaire,
par
simple
propagation
d'orbite.
Trouver
ensuite
le
transfert de Lambert allant du corps fictif
Mi à la date TM i au corps
réel Pi+ 1 à la date Ti+l".
Ce deuxième transfert (entre
Mi et Pi+! ) est de la même forme que
ceux du.4.4.2.1.;
On obtient ainsi ( N-I) problèmes de Lambert à
résoudre et (N-I) extrapolations.
13 1
variables du problème
Dans ces conditions les variables du problèmes sont les suivantes:
• les
N dates de survol
• les
N-I
dates de manœuvres
• les
3(N-1) composantes des vitesses relatives de départ des
corps survolés.
On a ainsi 5N-4 variables.
Le coût s'exprime sous la forme
N -1
N·I
J.
Lllâvi Il + L"âVMi Il
i=1
i=1
4.5.
Optimisation
Dans les deux cas de figure developpés ci-dessus, l'optimisation est
de la même forme et consiste à minimiser la fonction coût J sous les
contraintes
suivantes:
• vitesse de départ du premier corps inférieure ou égale à la vitesse
maximale fixée:
d
vl
SVIM
• altitude de survol d'une planète supérieure ou égale à l'altitude
minimale fixée:
rI
SrI m
Pour l'optimisation nous avons utilisé les codes suivants
• VA I3A, de la bibliothèque HARWELL, basé sur une méthode de
gradients.
Elle nécessite le calcul de la fonction coût et de ses dérivées
partielles par rapport aux variables du problème .
• M2QN 1 du même genre que le précédent, fourni par le club
MODULOPT (C. LEMARECHAL, I.N.R.I.A Rocquencourt).
132
• RECSIMU basé sur le Recuit Simulé que nous avons écrit
spécialement pour résoudre ce problème. Ce code n'utilise
aucun calcul de dérivées.
4.6.
Calcul
de
la
roncHon coût
Pour utiliser les codes d'optimisation locale VAI3A et
M2QNI
on
a
besoin:
- du calcul de la fonction coût,
- du calcul des dérivées partielles.
Les calculs de dérivées partielles sont intégrés numériquement dans les
sous-programmes
qui
calculent
également,
à
différentes
étapes,
la
fonction
coût.
En
effet,
en
même
temps
que
les
résolutions
des
problèmes
de
Lambert,
les
sous-programmes
appropriés
calculent
les
dérivées partielles correspondantes voir
[CNES 2).
calcul de la fonction coût:
C'est la somme totale de tous les modules des incréments 11 v .
Les
incréments
réalisés
lors
des
survols
sont
de
trois
natures
différentes:
• incrément tenant compte de la contrainte sur la vitesse de départ
du premier corps,
• incrément de swing-by effectué lors du survol des planètes.
• incrément réalisé lors du survol de corps d'attraction
négligeable: astéroïdes, comètes, corps fictifs.
4.6.1.
Incrément
relatlf
il
la
contrainte
On veut traiter la contrainte:
133
Le principe est le suivant : si la vitesse de la sonde au dé part du premier
corps
est
trop
grande, on
la corrige
au
moyen
d'une
manœuvre
supplémenraire.
L'Incrément de vitesse a pour expression:
--+
."...
Av 1 • 0
si
--+
--+d
Av 1 •
v 1
si
4.6.1,
Cas de
survol d'une
planète
Il ,'agil de
prendre en compte de
manière analytique
la
contrainte
d'ahilude minimale de survol :
ri 1 ri m
L'angle de déviation maximale du cône de swing-by,
8i M• est donné par:
8 M • 2 Arc sin ( --_.::....---:---
i
1 +
L'angle de déviation nécessaire. 8i• s'exprime par
--+ a
--+ d
v i · v i
8 i = Arc cos (
a
d
v i
vi
134
La manœuvre
optimale de swing-by peut être effectuée dans deux cas:
• si via> vid , la manœuvre a lieu, avant le survol de la planète: il
faut diminuer la vitesse d'arrivée pour que
son
module soit égal au
module de la vitesse de départ pendant le survol, et fournir la rotation
qu'il fallait.
La déviation maximale 9iM
est calculée à partir de vid .
• si v,a 1 v,d
la manœuvre a lieu après le survol de la planète: il
1
l '
faut
augmenter
la
vitesse d'arrivée
et
fournir
l'éventuel
surplus' de
rotation.
La déviation maximale ai M
est atteinte pour via.
135
Dans le calcul de l'incrément de swing-by, deux cas sont à distinguer:
• si
8i 1 8im alors
• si
8i > 8im alors
4.&.3.
Cas
de
survol
d'un
corps
d'attraction
négligeable
(réel
ou
fictif)
l'incrément est donné par la différence vectorielle des vitesses absolues
d'arrivée et de départ de la sonde :
soit encore
136
4.7. Code Général I10PTIMA
4.7.1
le
programme
général
Rappelons
que
le
but
du
programme
VESTA
est
l'étude
d'une
mission
spatiale (horizon
1994)
pour la visite de
gros
astéroïdes
se
trouvant entre les
planètes
MARS
et JUPITER.
Après
une
étude de
faisabilité, qui a montré plusieurs possibilités pour une telle mission, il
s'agit d'étudier de près les trajectoires ainsi trouvées, pour choisir en fin
de
compte les
meilleures. Ce qui
conduit à les
optimiser en
tenant
compte de facteurs coûts et de contraintes supplémentaires.
Il est prouvé que pour Je survol de N corps, en notant:
Tj
i= 1,...,N
les N dates de manoeuvres lors du survol des corps
T"] i= 1,..,N-I les N-I
dates de manoeuvres intermédiaires
vrjd i=I, ..,3(N-l) les 3(N-l) vitesses relatives de départ des corps i,
ces
5N-4
paramètres
suffisent
à
déterminer
entièrement
les
trajectoires de survol de tous les corps i.
Pour l'optimisation des
trajectoires,
les
variables
à considérer seront
donc ces paramètres qui serviront ainsi à déterminer la fonction coOt et
à traiter mathématiquement le problème.
En définitive le problème posé revient à résoudre
M I N J
(P)
{
sous contraintes C
où
M
1)
J=L Il li v i Il
j",1
M = N . I
dans le cas s an s manœuvres intermédiaires
avec
{M=2N-2 dans le cas avec manœuvres intermédiaires
2)
les contraintes C pouvant être :
- d'altitude minimale de survol
de vitesse maximale (de départ par exemple)
. de vitesse minimale (de survol).
137
Le programme général HOPTIMA répond à cet objectif.
Descriplion des sous-pro~rammes de HOPTIMA
Le programme général comporte près de 50 sous-programmes écrits
en fortran 77, sans compter le code d'cptimisarion. Comme cela est dit
en
4.4., les
calculs des
trajectoires, de
leurs COÛIS et
des
dérivées
partielles,
sont
faits
après
résolution
de
différentes
équations
de
la
mécanique spatiale: problèmes aux deux corps, problèmes de
Lambert,
El c'est donc toul naturel que ce soit les sous-progammes résolvant ces
problèmes qui donnent aussi les paramètres pour calculer les valeurs de
la fonction COÛI el de ses dérivées partielles,
Enfin, notons que pour noIre part, nous
n'avons d'Intervention qu'au
niveau de l'initialisation des trajectoires et de leur optimisation, tout
le
resle est Iraité par les sous-programmes fournis par le C.N.E.S. .
On distingue les différems types de sous-programmes suivants:
• ceux nécessaires à la résolu lion du problème de LAMBERT el qui
fournissent donc les paramètres de
la trajectoire courante pour chaque
étape de I'optimisation;
• ceux
calculant la fonction
• ceux déterminant les gradients à chaque étape;
• ceux appelant les codes d'optimisarion;
• les codes d'cptimisation proprement dits VA13A, M2QNl, RECSIMU ...
Nous avons exécuté le code général en fortran 77 sur C.D.C. (Cyber
860) avec le système NOSNE, au C.t.C.T. (Centre Interuniversitaire de
Calcul de Toulouse).
138
Schématiquement il est ainsi constitué:
Sous-programmes
Sous-programmes
appelés
principaux
CONsr
f-+
CONSTANlES
1
•
lBIXJN
-+
DONNEES
Données
nécessaires
•
la
résotution
du
probtëme
d 'oplimisalion
~
INIT
--..
INITIALISATION
,
lnirialisation
des
données
COlITUT
Calcul de la fonction coQt en
~
CAl.CUL DES 6V
donc de la fonction COÙI
unlisant dcs
méthodes de
calcul de mécanique
sp ati a le
!
VALMIN ~
OPTIMISATION
On
minimise la fonClion coût
de la fonction COÙI
en faisant appel • un code
doptimisation
(locale
ou
globale
selon
disponibilitl!)
~
EDITION
des
RESULTATS
139
4.7.2.
Entrées/Sorties
de
HOPTIMA
ENTREES
En entrée, on donne une
trajectoire
Initiale de
la manière suivante:
- on lit plusieurs fichiers:
un fichier DONNEES donne :
• les N corps à survoler (dont la terre)
• les N dates de leur survol
• les altitudes minimales de survol
un fichier PAROPTIMMA_DON donne:
• la masse initiale de la sonde
• les coefficients de pondération sur les dates et vitesses
un fichier UNF8599_DlR donne:
• pour les corps à survoler les informations:
nature du corps (planète ou astéroïde)
différents paramètres (rayon au soleil, taille ... )
orbites
- en interactif on peut choisir:
•
de figer certaines dates et/ou certaines vitesses
•
le code d'optimisation
• de tester le calcul des dérivées par différences finies
• de réoptimiser (avec le même code ou avec un autre)
SORTIE
En sortie, on a une trajectoire optimisée complète avec la possibilité de
pouvoir sortir les détails de:
• dates précises (à l'heure près) de survol des N corps
• dates de
manœuvres intermédiaires
• détails sur les  V : Vitesses arrivée/départ
• somme totale des  V
(i.e. J optimisée)
• masse totale optimisée
On peut aussi suivre les itérations et ainsi voir:
• les valeurs de J
• éventuellement les
gradients
• et bien d'autres paramètres suivant le code d'optimisation choisi.
140
4.8.
Résultats
Numériques
Dans ce qui suit nous
présentons quelques résultats numériques de
['optimisation de trajectoires avec les codes VAI3A, M2QNI et RECSIMU.
Les données initiales qui fournissent
le choix des corps et les dates
de survol déterminent une trajectoire initiale qu'il faut optimiser.
4.8.1
Trajectoires
initiales
Tous les numéros négatifs de corps signifient qu'il s'agit d'astëroîdes.
DONNEES J.
Numéro
Nom du corps
date de survol
commentaires
corps
3
Terre
29 Il 1996
Altitude mini.
départ:
110 km
4
Mars
28 06 1997
Altitude mini.
survol: 10 km
-46
Hestia
26 05 1998
-2942
1932BG
04 04 2000
-2599
Veseli
04 Il 2000
-80
Sapho
10 04 2001
141
DONNEES 2
Numéro
Nom du corps
date de survol
commentaires
corps
3
Terre
26 Il 1996
Altitude mini.
départ: 110 km
4
Mars
07 07 1997
Altitude mini.
survol: 10 km
-46
Hestia
25 05 1998
-2942
1932BG
08042000
-2599
Veseli
04 Il 2000
-80
Sapho
07 04 2001
DONNEES 3
Numéro
Nom du corps
date de survol
commentaires
corps
3
Terre
03 12 1996
Altitude mini.
départ: 110 km
4
Mars
04 07 1997
Altitude mini.
survol: 10 km
-46
Hestia
19 05 1998
-2942
1932BG
05042000
-2599
Veseli
06 Il 2000
-3216
190RB
0508 2002
142
DONNEES 4
Numéro
Nom du corps
date de survol
commentaires
corps
3
Terre
03 12 1996
Altitude mini.
dénart: 110 km
4
Mars
04 07 1997
Altitude mini.
survol: 10 km
-46
Hestia
1905 1998
-2942
1932BO
05042000
-2599
Veseli
06 Il 2000
-80
Sapho
05 08 2002
DONNEES S
Numéro
Nom du corps
date de survol
commentaires
corps
3
Terre
01 Il 1994
Altitude mini.
départ: 110 km
4
Mars
1008 1995
Altitude mini.
survol: 10 km
-4152
10 12 1996
4
Mars
15 08 1997
Altitude
'
,
mini,
survol: 10 km
-2087
05 05 1998
-1
03 03 1999
143
4.8.2.
Résultats
Tous les résultats qui
suivent sont obtenus après optimisation prenant
en compte à la fois
:
• l'introduction de
manœuvres intermédiaires
• l'optimisation en masse (voir (CNES 3]
).
Par ailleurs on a considéré les valeurs suivantes:
• masse initiale
=
2000. 0000
Kg
• masse larguée
=
500.0000 Kg
• gisp
=
2.7860 Km!s
• Dv de freinage
=
.1000 Km/s
4.8.2.1.
Exemples
de
résultats
complets
d'optimisation
locale
Pour chacun des résultats nous présentons:
• un tableau: donnant les dates de survols et les dates de maœuvres
intermédiaires,
• la masse optimisée: correspondant à la masse sonde plus la masse
de
la charge
utile; celle masse est
donnée
pour
les
deux
types de
résolutions, à savoir optimisation avec ou sans contraintes,
• l'incrément total de vitesse: c'est la somme de tous les modules
d'incréments; on donne les détails de celle valeur en la décomposant en
deux parties l'une avant le largage du pénétrateur (dans nos exemples
c'est toujours le troisième corps), et l'autre après le largage.
144
Trajectoire optimisée la partir de DONNEES 1
(par
M2QNI)
Numero Nom du
date de survol
date de manœuvre inter-
Corps
Corps
médiaire
3
lERRE
02 12 1996
à
2111
06 03 1997
à
0911
4
MARS
04 07 1997
à
09 "
06 12 1997
à
23"
-46
HESTIA
21 05 1998
à
22U
26 05 1999
à
0211
-2942
1932RB
05 04 2000
à
03"
13 07 2000
à
0711
-2599
VESELI
05 II 2000
à
om
12 01 2001
à
1711
-80
SAPHO
10 04 2001
à
02H
M I ·
. é
7 16
sans contraintes
1
. 1 Kg
asse tota e optJmls e :
716.5 Kg
avec contraintes
Incrément de vitesse total : 1.743 Km/s
.778
avant le survol du 3 e m e c o r p s
dont: 1 .965
après le largage du pénétrateur
Trajectoire optimisée la partir de DONNEES 2
(par
M2QNI)
Numero Nom du
date de survol
date de manœuvre inter-
Corps
Corps
médiaire
3
lERRE
01 12 1996
à
23"
17 03 1997
à
0911
4
MARS
04 07 1997
à
14 12 1997
à
1311
07 "
-46
HESTIA
22 05 1998
à
0011
26 05 1999
à
1511
-2942
1932RB
05 04 2000
à
03"
20 07 2000
à
2211
-2599
VESELI
05 II 2000
à
O6H
18 01 2001
à
2211
-80
SAPHO
10 04 2001
à
0211
M
1
irnisé
1715.9Kg
sans contraintes
asse tota e oplImls e:
7 16.3 Kg
avec contraintes
Incrément de vilesse total : 1.743 Km/s
.778
avant le survol du 3 e m e corps
dont: 1 .965
après le largage du pénétrateur
145
Trajectoire optimisée à partir de DONNEES 3
(par M2QNI)
Numero Nom du
date de survol
date de manœuvre inter-
Corps
Corps
médiaire
3
TERRE
03 12 1996 à
22H
16 03 1997
à
02H
4
MARS
29 06 1997 à O9H
29 10 1997
à
OOH
-46
HESTIA
23 05 1998 à
IOH
03 06 1999
à
OOH
-2942
1932RB
05042000
à
12H
04 07 2000
à
I3H
-2599
VESELI
06 Il 2000
à
02H
22 10 2001
à
07H
- 3216
I90RB
05 08 2002
à
19H
M I ' . é
1 850.7 Kg
sans contraintes
asse tota e oplJmls e:
850.9 Kg
avec c o nrr a in te s
Incrément de vitesse lotal : 1.258 Km/s
.788
avant le survol du 3 e m e corps
dont: 1 .470 après le largage du pé n étr at eur
Trajectoire optimisée à partir de
DONNEES 4 (par VA I3A
après 7
réoptimisations)
Numero Nom du
date de survol
date de manœuvre inter-
Corps
Corps
médiaire
3
TERRE
04 12 1996 à
22H
18 12 1996
à
08H
4
MARS
03 07 1997 à O9H
05 07 1997
à
19H
-46
HESTIA
23 05 1998 à
18H
04 07 1999
à
O4H
-2942
1932RB
05 04 2000
à
08H
07 04 2000
à
08H
-2599
VESELI
05 Il 2000
à
22H
07 Il 2000
à
22H
- 80
SAPHO
09 04 2001
à
21H
M I '
. é
1 706.9 Kg
sans c on rr ai nte s
asse Iota e oplJmls e: 707.2 Kg
avec c o nrr ain re s
locn!meot de vitesse total : 1.780 Km/s
.776
avant le survol du 3 em e corps
dont: 1 1.004 après le largage du pénétrateur
146
Trajectoire optimisée è partir de DONNEES 5
(par VAI3A)
Numero Nom du
date de survol
date de manœuvre inter-
Corps
Corps
médiaire
3
lERRE
30 10 1994 à 01"
Il 01 1995
à
03"
4
MARS
24 08 1995 à
19 "
16 05 1996
à
0711
-4152
13 12 1996 à
19"
09 03 1997
à
1811
4
19 08 1997 à
18"
02 12 1997 à 06"
-2087
10 05 1998 à
22"
22 08 1998
à
1411
-1
03 03 1999 à
15"
10 3 9 .0 Kg
sans contraintes
Masse totale opti misée: 1 1039.4 Kg
avec contraintes
Incrément de vitesse total: .947 Km/s
.216 avant le survol du 3 e m e corps
dont: 1 .731 après le largage du pénétrateur
147
4.8.1.1.
Résultats
d'optimisation
globale
Nous
exposons
ici,
brièvement,
nos
remarques
à
la
suite
de
l'utilisation du code RECSIMU.fOR pour la résolution du problème (P).
Pour plus de détails voir [UAMMA91).
OjWculth
repcoptrées
Dar
le
code
A
la
fin
du
temps
contractuel,
nous
n'avions
pas
fini
encore
d'harmoniser
les
différents
paramètres
entre
le
code
d'optimisation
globale et
les autres sous-programmes du
code général HOPTIMA. En
effet,
pendant
longtemps
nous
avons été
confrontés
à
des
problèmes
d'overflow. Cela est dû, en grande partie, aux difficultés de la méthode
du Recuit Simulé à s'adapter au type de résolution adoptée dans le code
général. On pourrait résumer ces difficultés comme suit:
• plateaU!
Si on résume le calcul de la fonction coût 1 on pourrait dire qu'elle
est de la forme:
JI
Si
Xe
CI
12
Si
Xe
C2
1 = {
Jk
Si
Xe
C k
où CI, C2, ... , Ck sont des sous-ensembles de R26.
Autrement dit,
J est une ronction à "plateaux". Dans certains
cas en
modifiant
légèrement les variables (ce qui est inévitable pour la bonne
marche
du
Recuit),
la
fonction
coût
reste
constante
dans
un
"bon"
voisinage (sur CI
par exemple). Le
Recuit, dans ces cas, a du mal à
trouver de vraies descentes; au plus arrive-t-il
à trouver des descentes
relatives! et donc au bout d'un certain nombre d'itérations (qu'on a pris
de plus en plus grand) on aueint le maximum de tirages autorisés entre
deux
vraies descentes.
Ces plateaux proviennent de troncatures faites dans le calcul des
11 V. Elles ont été nécessaires poui
intégrer les contraintes et permettre
de ce fait l'utilisation (dans un premier temps par le C.N.E.S.) de VA I3A
qui est un code d'optimisation sans
contraintes. En effet, pour VA BA
le problème est résolu comme s'il était sans contraintes. Car, en pratique,
grâce aux troncatures, pour tout jeu de variables (soumis par VA 13A) il
était
possible
(dans
le
code
général) de
lui
calculer
une
valeur
de
fonction.
Malheuresernent pour le Recuit cela est fatal.
148
• tirue au hasard dans Rz..ti
Une autre difficulté pour le Recuit, est causée par
le tirage au
hasard autour de Xk (pour passer de Xk à Xk+ 1) dans R 26 . En effet cela
devient "pénible" de sortir des états de minimums locaux, surtout très
três difficile d'échapper aux "plateaux" dans
R 26.
Au cours de ce tirage, on est souvent amené à générer des variables Xi
qui correspondent à des dates incohérentes : dates antérieures à la date
de départ choisie (antérieure de beaucoup même, parfois!), dates avant
l'origine du calendrier ..Julien-CNES.....
et cela expliquerait en
partie
certains
overflows
et/ou
underflows
dans
les
sous-programmes
de
résolution du problème aux deux corps.
• beaucoup de descentes relatives
En
observant
en
détail
quelques
résultats
du
Recuit,
on
peut
remarquer que la méthode
marche même si elle est piégée par tout ce
qu'on
vient
de
dire.
En
effet,
on
remarque
pour
l'optimisation
de
DONNEES 2
- beaucoup de descentes relatives qui confirment qu'on a eu
beaucoup de remontées acceptées. Or, bien entendu, la force du Recuit
réside dans cette acceptation de remontées.
- une
convergence
autour d'un
point
tel
que
Y=O.
Pour
atteindre ce
point il a fallu
beaucoup de descentes (malheureusement
relatives) et
une précision de E-9. Mais en fin d'exécution, YO étant
meilleur (-1.47), c'est comme si
le Recuit s'était fait
piéger par
un
minimum local qui, de surcroît, est pire que Je point initial!
• peu de vraies
descentes
On a observé quand même quelques vraies descentes en jouant sur
la matrice de covariance intiale VOAMO. En effet en prenant
.
1 0.4
si i = indice de date
VGAMO(I) = 0.04
si i = indice de vitesse
on est arrivé à avoir quelques vraies descentes avant que le code
se
fasse piéger à nouveau par un plateau.
149
4.9.
Conclusion
le
problème
pratique
étudié
dans
ce
c hapitre
est
bien
un
problème
d'opu mlsatton
globale.
En
effel,
les
résultats
d'expérimentation
numérique
ont
montré
l'existence
de
plusieurs
minimums.
Pour lui trouver des solutions globales. nous avons choisi la méthode du
Recuit
Simulé. En
effet, après avoir étudié
plusieurs
possibilités de
méthodes, seule la méthode du Recuit nous avait semblé être adaptée à
la situation dont les impératifs étaient:
• on ne disposait pas de formule explicite de la fonction objectif en
fonction de ses variables;
• l'espace des variables était mixte: réels (vitesses) el entiers (dates).
• pour l'optimiseur (le code) c'est comme si on ne disposait que d'une
"boîte noire" pour le calcul de fonction: pour un jeu d'entrée, on récupère
une valeur de fonction à la sortie;
• l'optimisation consistait, en gros. à considérer une trajectoire initiale,
la modifier légèrement pour trouver une trajectoire optimale; en effet, il
y
a
quelques
jours
seulement
de
différence
entre
les
dates
d'une
trajectoire initiale et celle de la trajectoire optimisée;
• la différentiabilité de la fonction J (par rapport aux dates) n'est pas
partout
garantie!
•
la
dimension
du
problème
(5N·4)
eSI
élevé:
jusqu'à
36
pour
seulement 3 à 4 asteoïdes à visiter en repassant par Mars;
Or le Recuit, même s'il n'est pas plus efficace que d'autres méthodes
pour les problèmes de
grande dimension, est,
néanmoins.
très
adapté
aux
autres
impératifs ci-dessus cités.
En
effet.
le
Recuit
Simulé
ne
nécessite, pour son exécution, que des
valeurs de fonction pour un jeu
de variables X=(X l , .... XN). Au cours de celte exécution, l'état (Le (X ,
f(X»
est
perturbé
vers
un élat
voisin dont
on
calcule la
valeur de
fonction, etc ... Et cela permet à l'algorithme de construire un modèle qui
converge vers un minimum global.
C'est pourquoi nous avons programmé le code RECSIMU.FOR et l'avons
renforcé par des cycles de réchauffement (cf. [BONNEMOY & HAMMA
199Ia)).
Ce code, malheureusement. n'a pas donné les résultats escomptés dans
les délais (cf.
[HAMMA
1991)) pour des raisons
que
nous
pensons
d'harmonie avec les sous-programmes propres au calcul des trajectoires.
de leurs paramètres et de leurs coûts.
1.50
RHérences
du
Chapitre
4
[CNES
IJ
HOPTlMA: optimisation de trajectoires interplanétaires avec
manœuvres
intermédiaires,
CT!DTI/ MS! MN! 283 , Toulouse 06 - Il -1986
[CNES 2] Optimisation de trajectoires interplanétaires avec contraintes
supplémentaires. Application au projet VESTA,
CT/DTI! MS! MN! 321 , Toulouse 16 - 12 - 1986
[CNES 3J VESTA: Modifications du programme HOPTIMA d'optimisation
de
trajectoires
interplanétaires,
CT/DTI/ MS! AS! 194 , Toulouse 07
- 07 - 1988
[CNES 4] Test d'optimalité et animation de trajectoires interplanétaires ,
CT!DTI! MS! MN! 127 .
[CNES .5J VESTA détermination de trajectoires créneau 1994,
CT!DTI! MS! MN! 106 , Toulouse 06 - 05 - 1986
[DONNEMOY & HAMMA 1991a) C. DONNEMOY & S.D. HAMMA, La
méthode du recuit simulé: optimisation globale dans Rn, APII, Aut, Prod.
Inr. Ind.,vol 25, W5, 1991, p. 477-496.
[HAMMA 1991) S.D. HAMMA, Etude de Problèmes d'optimisation de
trajectoires interplanétaires, Rapport final de la phase D, Marché C.N .E.S
- 873 - CNES - 87 - 4937, Juillet 1991.
15 1
CONCLUSION GENERALE
Le
travail présenté dans ce
mémoire a été réalisé dans le cadre
d'un
marché
contractuel
établi
entre
la
division
Mathématiques
Spatiales
du
Centre
National
d'Etudes
Spatiales
de
Toulouse
et
le
Laboratoire
d'Analyse
Numérique
de
l'Université
Paul
Sabatier.
Ce
marché portait sur "l'étude de problèmes d'optimisation globale".
Dans l'état actuel de nos recherches nous pensons avoir fait l'état
de l'an sur les méthodes numériques existentes. En effet, suite à notre
étude, nous pouvons souligner dans ce qui suit, en conclusion, quelques
grands
traits.
• Existence de deux grandes classes d'approches pour la résolution
numérique du
problème d'optimisation globale:
- l'une stochastique utilisant des stratégies probabilistes pour construire
un modèle convergeant vers une solution,
- l'autre déterministe.
•
La
diversité
des
philosophies
au
sein
d'une
même
classe
d'approche.
•
La
popularité
des
méthodes
stochastiques
par
rapport
aux
déterministes;
en
effet
elles
sont d'implémentation
relativement
plus
facile et ne
sont pas très exigentes sur la
fonction objectif (sur des
hypothèses
supplémentaires)
comme
c'est
le
cas,
malheureusement,
pour les méthodes déterministes.
• Des solutions ont été souvent apportées, de façon locale, à tel ou
tel problème spécifique par un algorithme conçu en conséquence et qui
ne marche que pour ce genre de cas isolé.
•
Aucune
bibliothèque
scientifique
ne
propose
de
codes
d'optimisation globale. Les quelques rares embryons de
codes à portée
un peu plus élargie sont d'accès très difficile (confidentialité?).
• Beaucoup de
difficultés persistent encore:
prise en
compte de
contraintes, problème à dimension élevée, etc ...
Numériquement,
notre
contribution
a
été
de
montrer
que
pour
certaines méthodes que nous avons étudiées de près (cf. chapitre 3), il
est toujours possible d'en améliorer les
algorithmes pour trouver des
implémentations dont les succès sont meilleurs. C'est ainsi que pour le
Recuit nous avons introduit des cycles de réchauffement pour améliorer
le pourcentage de succès, el pour le Polytope Mouvant (qui n'est pas un
152
code
d'optimisation
globale)
en
introduisant
la
recherche
adaptée
on
arrive
à
converger
"un
peu
plus"
vers
une
solution
globale.
Pour
beaucoup d'autres
méthodes nous avons des versions de
codes qui
ne
nous
ont
pas
encore
donné
satisfaction
mais
pour
lesquelles
nous
envisageons de
continuer à travailler pour "forcer la
convergence":
le
Tunnelier, les Initialisations Multiples, etc ...
rn raison du marché que nous avions avec le C.N.E.S., nous avons
centré
tout
notre
travail
sur
la
réalisation
d'un
code
pour
traiter
le
problème présenté au chapitre 4. El c'est pour cela que nous avons un
peu plus insisté sur l'implémenlation des méthodes et
leur convergence
numérique
et
avons
mis
un
peu
de
côté
les
aspects
théoriques:
théorèmes
de
convergence
et
leur
preuve
théorique.
Ces
derniers
aspects font
partie de
nos
projets futurs. Nous avons donc spécialement
amélioré
le
Recuit
Simulé
pour traiter
le
problème
d'optimisation
de
trajectoires interplanétaires présenté par le C.N.E.S.
Pour le
valider il
nous a été indispensable de
le comparer à d'autres codes. En
pins des
codes (PRSAN et PRSMD2) que nous a fournis le Professeur BONNEMOY
de
Clermont-Ferrand l, nous en avons programmé quelques autres. Et à
ce stade, nous avons pu apprécier les difficultés que posent les con nits
entre la convergence théorique d'une méthode et ses succès numériques.
C'est le cas de la méthode du Tunnelier, des Initialisations Multiples, de
la méthode de PilA VSKIJ et bien d'autres dont nous avons des versions
de
programmes non
satisfaisants pour le momemt.
Enfin. concluons que "avancement de
nos
recheches sont, à l'état
actuel, tels que nous puissions dire que nous disposons de suffisamment
d'éléments
pour
envisager
l'élaboration,
dans
l'avenir,
d'un
petit
"package" informatique d'optimisation globale proposant quelques codes
stochastiques
et
déterministes
basés
sur
tous
les
embryons
de
codes
dont
nous
disposons.
Ce
package
dans
un
premier
temps
inclura
le
Recuit
Simulé, le
Tunnelier,
le
Polytope
Mouvant et
les
Initialistions
Multiples.
Ce
genre
de
logiciel
aura
pour
ambition
principale
d'améliorer
un
minimum
local
dans
un
voisinage
de
plus en
plus
large.
Puis
nous
essayerons
de
voir dans
quelles
mesures
on
peut
réaliser
des
codes
numériques
basés
sur
des
études
déjà
faites:
par
exemple
dans
[JIIRIART-URRUTY
1989J
des
conditions
nécessaires
et
suffisantes
d'optimalité
globale
sont
données
pour
certains
types
de
problèmes
(d.c.).
Tout cela, en
pins de
la justification théorique de
la convergence
du Polytope Mouvant vers un minimum global, fera partie de nos futurs
projets de
recherche.
153
GLOSSAIRE
Traduction en Français de termes anglais
Local
minima
enumeratlon : énumération des minimums locaux.
Tunneling Method
: Méthodes du Tunnelier.
Simulated
Annealing: Recuit Simulé.
Multlstart : Initialisations multiples.
S.A.D.: Seteepest Ascent Descent : Montagnes
Russes.
Genetlc
A Igorithms : Algorithmes Génétiques.
Grid
Search Methods : Méthodes de Recherche par Quadrillage.
Trajectory
Methods : Méthodes des Trajectoires.
Determinlstlc
Methods : Méthodes Déterministes.
Stochastic
Methods : Méthodes Stochastiques.
Random
(line)
Search
Methods : Méthodes à Recherche(linéaire)
Aléatoire ou Méthodes à Stratégie Aléatoire.
Sampling
Methods
Méthodes par Echantillonnage.
Partition
&
Search
Partition & Recherche.
Approximation of f & Search : Approximation de r & Recherche.
Global
Decrease
or
Permanent
Inprovement
ln
f
:
Décroissance
Globale ou Amélioration permanente de r.
Methods
with
guaranteeed
convergence
Méthodes avec garantie
de
convergence.
Blind Search : Recherche Aveugle.
Combination of Blind &
Search : Combinaison de Recherches
Aveugle et Locale.
Non local Search : Recherche Non Locale.
154
Coverlng
Methods
Méthodes par Recouvrement.
Methods
based
on
Random
Sampling: Méthodes basées
sur
l'échantillonnage Aléatoire.
Methods
based
on
a
stechasrle
model
or
the
objective
runction:
Méthodes basées sur un modèle stochastique de la fonction objectif.
Sequential
Transition
rrom
one
local
minimum
te
another:
Méthodes de transition d'un minimum local à un autre.
Coverlng
Methods : Méthodes par Recouvrement.
Methods
with
guaranteeed
aecuracj : Méthodes avec garantie de
précision.
Direct
Methods : Méthodes Directes.
Clusterlng
Method. : Méthodes de Regroupement en Grappes.
Cluster: Grappe.
Generalized
Descent : Descente Généralisée.
Metbod. apprexlmatlng th. L•••I Set.
:
Méthod. . d'approximat lon
\\.
des surfaces de niveau.
f)
Methods
approxlmatlng
the
Objective
Functlon
Méthodes
d'approximation de la fonction objectif.
Pure Random
Search : Recherche Aléatoire Pure.
Hranch & Round : P.S.E.P.:
Procédure de Séparation et Evaluation
Progressi ve.
RESUME
Le
présent
travail
qui
porte
sur
les
méthodes
numériques
d'optimisation globale, a pour origine lin sujet proposé par le C.N.E.S.
(Centre National d'Etudes Spatiales) de Toulouse "Etude de problèmes
d'optimisation globale" (marché N° 873-CNES-87-4937). Le but était de
développer des codes d'optimisation globale eil vile d'applications à des
problèmes
spécifiques.
Au chapitre I , nous rappelons les
'"
problématiques de J'optimisation
locale ct de l'optimisation globale d'un point de vue numérique.
..
Nous nous
attachons, au chapitre 2, à faire
l'état de
l'urt : des
méthodes et codes existants en optimisation globale, en les
classant
suivant différents critères: leur philosophie, leur stratégie, leur garantie
de convergence, etc.
'
Ensuite, nous étudions de près, au chapitre 3, certaines méthodes
numériques:
le Recuit Simulé, la méthode du Tunuelier, le Polytope
Mouvant.
Dans
ce
chapitre,
nous
présentons
aussi
des
résultats
numériques de codes que
nous
avons
écrits, testés et améliorés. Ce
chapitre contient également, à
la suite de plusieurs expérimentations
numériques, des détai ls de choix de paramètres, conseils d'utilisation,
etc. pour certains algorithmes.
Enfin nous abordons, au chapitre 4, le problème d'optimisation de
trajectoires interplanétaires tel qu'il nous a été présenté par le C.N.E.S.
Il s'agit de la trajectoire d'une sonde visitant de gros astéroïdes entre les
planètes Mars ct Jupiter. Afin de recueilJ ir le plus de renseignements
possible, on cherche à emporter le maximum de matériels scientifiques.
Pour
gagner
cil!
poids,
on
cherche
à
embarquer
le
moins
d'ergol
(carburant) possible, puis on essaie d'en consommer le moins possible
au cours de
la mission pour pouvoir aller le plus loin
possible. La
trajectoire
d'une
telle
sonde,
pour
remplir
sa
mlsslon ,
avec
des
hypothèses
de
survol
(vitesses
minimales
ou
maximales
de
survol,
altitudes minimales ou maximales de survol, etc.) sur
les
différents,
corps,
peut
être
considérée
comme
une
suite
de
manœuvres
impulsionnellcs: pour corriger la trajectoire, pour freiner, pour accélérer,
etc. Ainsi, outre l'utilisation de l'effet Swing-by (011 effet dc "tremplin
gravitationnel") qui permet de modifier la direction de la sondé 'sans
effectuer
de
manœuvres
coûteuses,
il
s'agira
pour
minimiser
la
consommation d'ergol, de minimiser la somme
totale des modules des
incréments de vitesses réalisés lors des différentes manœuvres.
MOTS CLES:
Oplimisalion
(;tohale,
Recuit
Simulé,
Cycles
de
Réchaufr~'~
tope
Mouvanl,
Méthode
du
Tunnelier,
Trajectoires
int
.