Algorithmes d'optimisation de la population : Essaim de Particules (OEP ou PSO en anglais)

      Elles forment des essaims distincts pour flotter dans les rayons du soleil

                                        ou pour suivre les nuages d'orage. Il est possible qu'ils tirent leur énergie des décharges atmosphériques.

Mais au moment du danger ou plus largement, lors d'un changement soudain qui menace leur existence, ils s'unissent...

Stanisław Lem "L'Invincible"

Sommaire

1. Introduction
2. Principes de l'algorithme
3. Implémentation classique
4. Version modifiée
5. Résultats des tests
   

1. Introduction

Nous sommes nombreux à avoir lu le merveilleux best-seller de science-fiction de Stanisław Lem, "L'Invincible". L'une des premières descriptions de l'intelligence "en essaim" est née précisément avec la publication de ce roman de science-fiction. L'histoire est celle de robots ayant survécu sans contrôle centralisé. Les spécimens les plus simples et les plus nombreux ont ainsi survécu, plutôt que les plus complexes, les plus intelligents et les plus puissants.

Au cours de milliers d'années de microévolution, ces machines ont appris à faire face à des concurrents qui les surpassaient en intelligence et en disponibilité d'énergie. Ils ont dû se battre non seulement contre d'autres robots, mais aussi contre le monde vivant de la planète. Les éléments fantastiques de cette œuvre peuvent être comparés de manière fiable à l'évolution et à la nature elle-même.

Depuis les temps les plus anciens, les hommes s'intéressent au comportement des animaux au sein d'un groupe (appelé comportement en essaim) : comment fonctionnent les oiseaux lorsqu'un troupeau se dirige vers des pays chauds, comment un essaim d'abeilles produit de la nourriture, comment une colonie de fourmis survit tout en créant des structures complexes, comment les poissons se comportent en essaim et pourquoi leur comportement est autant synchronisé. L'organisation des individus dans la société, qui présente certains des aspects d'un organisme intégral bien coordonné, encourage de nouvelles idées dans le domaine de l'optimisation algorithmique.

L'intelligence en essaim décrit la simulation du comportement collectif d'un système auto-organisé. Il existe dans le monde un nombre relativement important d'algorithmes de ce type. Dans la version canonique, écrite en 1995 par J.Kennedy et R.Eberhart, un modèle sous-jacent à cette méthode a été obtenu en simplifiant le modèle de Reynolds. À la suite de cette simplification, les individus distincts de la population ont commencé à apparaître comme des objets séparés n'ayant pas de taille, mais une certaine vitesse.

En raison de l'extrême similitude avec les particules matérielles, les objets simples qui en résultent sont appelés particules. Leur population est appelée essaim. À chaque instant (à chaque itération), les particules ont une position et un vecteur vitesse dans l'espace. Pour chaque position de la particule, la valeur correspondante de la fonction objective est calculée. Sur cette base, selon certaines règles, la particule modifie sa position et sa vitesse dans l'espace de recherche. Lors de la détermination de la position suivante d’une particule, les informations relatives à la meilleure position parmi toutes les autres particules voisines, correspondant aux tâches de la fonction d'aptitude, sont prises en compte.

Exemples d'algorithmes en essaim :

• Méthode de l'essaim de particules (PSO)
• Algorithme des fourmis
  
• Algorithme des abeilles
  
• Système immunitaire artificiel
  
• Algorithme du loup gris
  
• Algorithme des chauves-souris
  
• Algorithme de recherche par gravité
  
• Algorithme d'altruisme
  
• et bien d'autres...

Le passage de la modélisation du comportement collectif à l'optimisation collective repose sur l'idée biologique suivante : les organismes s'unissent en colonies pour survivre et améliorer leurs conditions de vie. Chaque organisme de la colonie a, en moyenne, de meilleures chances de survie dans la lutte contre les prédateurs, la colonie peut rechercher, traiter et stocker la nourriture plus efficacement que les organismes individuels, etc. En d'autres termes, toute colonie d'organismes, tout au long de son existence, résout divers problèmes d'optimisation avec plus ou moins d'efficacité, par exemple, pour maximiser la quantité de nourriture tout en minimisant les pertes dues aux prédateurs. Ces considérations ont servi de base à l'élaboration de diverses méthodes d'optimisation mathématique.

L'essaim de particules est l'un des algorithmes d'optimisation les plus célèbres et les plus populaires depuis sa création. De nombreux auteurs de ses différentes implémentations affirment que l'algorithme est très efficace pour optimiser des fonctions complexes avec de nombreux arguments, et qu'il est même adapté à l'apprentissage des réseaux neuronaux.

Dans cet article, j'essaierai de déterminer si l'algorithme est réellement efficace pour résoudre des problèmes complexes. Dans la version classique de l'algorithme, ainsi que dans beaucoup de ses modifications, il existe des limitations importantes liées au fait que la fonction optimisée doit être lisse et continue, ce qui signifie qu'elle est totalement inadaptée aux fonctions discrètes. Dans la lignée de cette série d'articles, tous les algorithmes examinés seront modifiés afin d'éliminer les lacunes, ou du moins de faire en sorte que les algorithmes fonctionnent au moins sur le plan purement technique. En d'autres termes, tous les algorithmes doivent être indifférents au caractère lisse des fonctions (comme dans les problèmes des traders) et n'avoir aucune restriction sur ses paramètres.

2. Principes de l'algorithme

Si l'article précédent était une introduction au monde de l'optimisation, il ne couvrait pas le principe d'interaction du programme principal (que ce soit un EA, un script ou un indicateur) avec le cœur de l'algorithme d'optimisation. Il est important d'y prêter attention, car dans tous les cas, un lecteur attentif se posera la question de savoir pourquoi les algorithmes et les exemples sont écrits de cette manière. Les versions actuelles des algorithmes d'optimisation utilisent la fonction d'aptitude comme un objet externe, alors que l'algorithme est le principal programme exécutable.

La figure 1 ci-dessous montre un diagramme dans lequel l'algorithme envoie les paramètres optimisés à la fonction d'aptitude, et obtient en retour la valeur d'aptitude (critère d'évaluation). Ce système n'est pas adapté à la construction d'un programme destiné à résoudre les problèmes des utilisateurs, y compris des traders. Mais pourquoi est-ce gênant ? Par exemple, on ne peut pas effectuer un test sur tout l’historique.

Fig. 1 : Interaction entre l'algorithme d’optimisation (PSO) et la fonction d'aptitude

La structure de la figure 2 est beaucoup plus efficace. L'algorithme d'optimisation n'est pas un programme indépendant, mais un module séparé ou une sorte de "boîte noire". Ce module comporte les paramètres "min", "max" et "step" pour chaque argument optimisé. Le programme MQL reçoit des arguments optimisés sur demande et renvoie les valeurs d'aptitude (i.e. les valeurs de la fonction d'aptitude). Cette structure permet de construire des solutions très flexibles, allant de l'utilisation de l'optimisation automatique dans les Expert Advisors à l'écriture d'un gestionnaire d'optimisation personnalisé.

Fig. 2 : Interaction entre le programme MQL et l’algorithme d’optimisation (PSO)

Il important également de mentionner que l'organisation des appels aux méthodes des algorithmes d'optimisation (bloc MQL dans la figure 2) peut être représentée par un schéma général qui est le même pour tous les algorithmes d'optimisation (AO) :

Initialisation_АО_0

Cycle d'itérations (époques)
{
1) Méthode_АО_1
2) Obtention des valeurs d’aptitude (fitness) pour chaque variante des paramètres optimisés
3) Méthode_АО_2
}

Nous pouvons ainsi constater que seules 3 méthodes publiques sont utilisées : Initialisation_АО_0, Méthode_АО_1 et Méthode_АО_2. Tout cela suffit pour organiser le processus d'optimisation dans les projets des utilisateurs, quelle que soit leur complexité.

Le flux de travail de l'Algorithme de l’Essaim de Particules (PSO) est illustré à la figure 3 ci-dessous et comprend les étapes logiques suivantes :

1. Génération de particules aléatoires (1ère itération)
2. Calcul de la valeur d'aptitude pour chaque particule
3. Calcul de la valeur d'aptitude pour toutes les particules en général
4. Réglage de la vitesse des particules
5. Stop, ou passage à l'étape 2
   
6. Fin du programme.

Fig. 3 : Flux de travail de l'Algorithme de l’Essai de Particules (PSO)

Examinons maintenant l'Algorithme de l'Essaim de Particules plus en détail.

Le système d'intelligence en essaim se compose de nombreuses particules qui interagissent entre elles et avec leur environnement. Chaque particule suit des règles simples, bien qu'il n'y ait pas de système centralisé de contrôle du comportement, pour dire à chacune d'entre elles ce qu'elle doit faire. Les interactions locales et aléatoires entre elles conduisent à l'émergence d'un comportement de groupe intelligent qui n'est pas contrôlé par les individus.
Si l'on fait l'analogie avec un troupeau, on peut dire que toutes les particules doivent accomplir des tâches simples :

• éviter de toucher d'autres particules
• ajuster leur vitesse en fonction de celle des particules environnantes
• essayer de maintenir une distance assez faible entre elles et l'environnement

L'algorithme PSO (Essaim de Particules) commence par l’initialisation de la population. La 2ème étape consiste à calculer les valeurs d'aptitude de chaque particule. Les meilleurs scores individuels et globaux sont ensuite mis à jour, puis la vitesse et la position des particules sont actualisées. Lors de l'utilisation du PSO, une solution possible à un problème d'optimisation numérique est représentée par la position de la particule. Chaque particule a également une vitesse actuelle qui reflète sa magnitude absolue et sa direction vers une nouvelle solution/position supposée meilleure.

La particule stocke aussi sa valeur d'aptitude, sa position actuelle, sa meilleure position connue (une position précédente avec la meilleure aptitude connue) et l'aptitude de la meilleure position connue. Les étapes 2 à 4 sont répétées jusqu'à ce que la condition d'achèvement soit remplie. Lors de la 1ère itération, toutes les particules sont dispersées pour trouver la meilleure solution (exploration). Chaque particule est évaluée. Les meilleures solutions pour la topologie du voisinage sont trouvées. Les meilleures particules personnelles et globales pour chaque membre de l'essaim sont mises à jour. La convergence est ensuite obtenue en attirant toutes les particules vers la particule ayant la meilleure solution. 

Bien que l'algorithme PSO soit assez simple à la base, il est important de bien le comprendre pour pouvoir modifier le code de cet article suivant nos besoins. PSO est un algorithme à processus itératif. À chaque itération de la boucle principale, la vitesse actuelle de chaque particule est d'abord mise à jour. La vitesse actuelle de la particule, ses informations locales et les informations globales de l'essaim sont prises en compte. La position de chaque particule est ensuite mise à jour en utilisant la valeur de la nouvelle vitesse de cette particule.

Mathématiquement, ces 2 équations de mise à jour des coordonnées des particules se représentent comme suit :

v(t+1) = w * v(t) + c1 * rp * (p(t) - x(t)) + (c2 * rg * (g(t) - x(t))

x(t+1) = x(t) + v(t+1)

La mise à jour de la position est en fait beaucoup plus simple que les équations proposées. La 1ère équation permet de mettre à jour la vitesse de la particule.

Le terme v(t+1) désigne la vitesse à l'instant t+1. La nouvelle vitesse dépend de 3 termes :

• w * v(t) — Le facteur w est la fraction de poids de l'inertie. Il s’agit simplement d’une constante. v(t) est la vitesse actuelle à l'instant t.
  
  
• c1 * rp * (p(t) - x(t)) — Le facteur c1 est une constante appelée fraction de poids cognitive (autrement appelée personnelle ou locale). Le multiplicateur rp est une variable aléatoire comprise dans l'intervalle [0, 1]. Le vecteur p(t) représente la meilleure position de la particule trouvée jusqu'à présent, et le vecteur x(t) représente la position actuelle de la particule.
  
  
• la mise à jour de la vitesse : c2 * rg * (g(t) - x(t) — Le facteur c2 est une constante appelée fraction de poids sociale (ou globale). Le multiplicateur rg est une variable aléatoire comprise dans l'intervalle [0, 1]. La valeur du vecteur g(t) est la meilleure position connue trouvée jusqu'à présent pour l'une des particules de l'essaim. Une fois la nouvelle vitesse v(t+1) déterminée, elle est utilisée pour calculer la nouvelle position de la particule x(t+1).

3. Implémentation classique

Une unité logique décrivant un ensemble de coordonnées dans l'espace (paramètres optimisés) est une particule, qui peut être représentée comme une structure, où c[] sont les coordonnées de la particule, cB[] sont les meilleures coordonnées de la particule pour toutes les itérations, v[] est la vitesse pour chacune des coordonnées de la particule, ff est la valeur actuelle de l'aptitude de la particule et ffB est la meilleure valeur de l'aptitude de la particule pour toutes les itérations. Dans le constructeur de la structure de particules, nous initialisons les valeurs de ff et ffB en utilisant la valeur minimale possible représentée par le type "double", puisque l'algorithme est conçu pour trouver le maximum de la fonction (pour trouver le minimum, il suffit d'ajouter un signe "-" devant la valeur d’aptitude résultante).

//——————————————————————————————————————————————————————————————————————————————
struct S_Particles
{
  public:
    double c  []; //coordinates
    double cB []; //best coordinates
    double v  []; //velocity

    double ff;    //the value of the fitness function
    double ffB;   //best value fitness function

    S_Particles ()
    {
      ff  = -DBL_MAX;
      ffB = -DBL_MAX;
    }
};
//——————————————————————————————————————————————————————————————————————————————

Mettons en œuvre l'algorithme PSO sous la forme d'une classe ne comportant que 3 méthodes publiques : InitPS(), Preparation() et Dwelling()(Initialization_АО_0, Method_АО_1 et Method_АО_2). Dans les méthodes privées, GenerateRNDparticles () et ParticleMovement () sont uniques pour le PSO, les autres ont déjà été examinées dans l'article précédent. Le tableau de structures p [] est un essaim de particules. En plus du fait que chaque particule possède des valeurs d'aptitude, ses propres coordonnées et les meilleures coordonnées, l'essaim dans son ensemble possède les meilleures coordonnées cB et la meilleure valeur d'aptitude ffB.

//——————————————————————————————————————————————————————————————————————————————
class C_AO_PSO
{
  public:
  //----------------------------------------------------------------------------
  S_Particles p    []; //particles
  double rangeMax  []; //maximum search range
  double rangeMin  []; //manimum search range
  double rangeStep []; //step search
  double cB        []; //best coordinates
  double ffB;          //FF of the best coordinates

  void InitPS (const int    params,       //number of opt. parameters
               const int    size,         //swarm size
               const double inertiaP,     //inertia
               const double selfBoostP,   //boost
               const double groupBoostP); //group boost

  void Preparation ();
  void Dwelling ();

  private:
  //----------------------------------------------------------------------------
  int swarmSize; //swarm size
  int parameters;//number of optimized parameters

  double inertia;
  double selfBoost;
  double groupBoost;
  bool   dwelling;

  void   GenerateRNDparticles ();
  void   ParticleMovement     ();
  double SeInDiSp             (double in, double inMin, double inMax, double step);
  double RNDfromCI            (double min, double max);
};
//——————————————————————————————————————————————————————————————————————————————

La méthode InitPS() est destinée à initialiser l'algorithme avant le début de l'optimisation(Initialization_АО_0). Les valeurs des arguments de la méthode sont affectées aux membres privés de la méthode, de même que la taille est affectée à l'essaim et aux paramètres internes de chaque particule de l'essaim.

//——————————————————————————————————————————————————————————————————————————————
void C_AO_PSO::InitPS (const int    paramsP,
                       const int    sizeP,
                       const double inertiaP,
                       const double selfBoostP,
                       const double groupBoostP)
{
  ffB = -DBL_MAX;

  parameters = paramsP;
  swarmSize  = sizeP;

  ArrayResize (rangeMax,  parameters);
  ArrayResize (rangeMin,  parameters);
  ArrayResize (rangeStep, parameters);

  dwelling = false;

  inertia    = inertiaP;
  selfBoost  = selfBoostP;
  groupBoost = groupBoostP;

  ArrayResize (p, swarmSize);

  for (int i = 0; i < swarmSize; i++)
  {
    ArrayResize (p [i].c,  parameters);
    ArrayResize (p [i].cB, parameters);
    ArrayResize (p [i].v,  parameters);
  }

  ArrayResize (cB, parameters);
}
//——————————————————————————————————————————————————————————————————————————————

La méthode Préparation () est d'abord appelée à chaque itération (époque)(Method_АО_1). La méthode est simple mais très importante. Selon que la méthode sera appelée à la 1ère époque ou aux époques suivantes (déterminées par le flag dwelling), la valeur de l'aptitude de l'essaim sera réinitialisée et une population aléatoire sera créée, ou les particules se déplaceront vers de nouvelles coordonnées.

//——————————————————————————————————————————————————————————————————————————————
void C_AO_PSO::Preparation ()
{
  if (!dwelling)
  {
    ffB = -DBL_MAX;
    GenerateRNDparticles ();
    dwelling = true;
  }
  else ParticleMovement ();
}
//——————————————————————————————————————————————————————————————————————————————

Une population d'essaims aléatoires est générée dans la méthode GenerateRNDparticles (). Les particules ont des coordonnées aléatoires dans l'intervalle spécifié pour chacune d'entre elles, et une vitesse aléatoire pour chacune des coordonnées.

//——————————————————————————————————————————————————————————————————————————————
void C_AO_PSO::GenerateRNDparticles ()
{
  for (int s = 0; s < swarmSize; s++)
  {
    for (int k = 0; k < parameters; k++)
    {
      p [s].c  [k] = RNDfromCI (rangeMin [k], rangeMax [k]);
      p [s].c  [k] = SeInDiSp (p [s].c [k], rangeMin [k], rangeMax [k], rangeStep [k]);
      p [s].cB [k] = p [s].c [k];
      p [s].v  [k] = RNDfromCI (0.0, (rangeMax [k] - rangeMin [k]) * 0.5);
    }
  }
}
//——————————————————————————————————————————————————————————————————————————————

La méthode ParticleMovement() déclenche l'algorithme de déplacement des particules vers de nouvelles positions. Pour cela, il est nécessaire de calculer la vitesse pour chaque coordonnée avec l'équation ci-dessus. Je ne sais pas pourquoi le terme "vitesse" est utilisé, car il s'agit essentiellement d'une valeur de déplacement ou, en d'autres termes, de la différence entre l'endroit où se trouve la particule à ce moment-là et l'endroit où elle devrait se déplacer. Après avoir calculé cette différence pour chacune des coordonnées, il suffit de l'ajouter aux valeurs actuelles. Vérifiez ensuite s'il est inacceptable de dépasser les limites min/max des paramètres optimisés (pour une particule, il s'agit de coordonnées) avec un pas donné.

//——————————————————————————————————————————————————————————————————————————————
void C_AO_PSO::ParticleMovement ()
{
  double rp;       //random component of particle movement
  double rg;
  double velocity;
  double posit;
  double positBest;
  double groupBest;

  for (int i = 0; i < swarmSize; i++)
  {
    for (int k = 0; k < parameters; k++)
    {
      rp = RNDfromCI (0.0, 1.0);
      rg = RNDfromCI (0.0, 1.0);
      
      velocity  = p [i].v  [k];
      posit     = p [i].c  [k];
      positBest = p [i].cB [k];
      groupBest = cB [k];

      p [i].v [k] = inertia * velocity + selfBoost * rp * (positBest - posit) + groupBoost * rg * (groupBest - posit);
      p [i].c [k] = posit + p [i].v [k];

      p [i].c [k] = SeInDiSp (p [i].c [k], rangeMin [k], rangeMax [k], rangeStep [k]);
    }
  }
}
//——————————————————————————————————————————————————————————————————————————————

La méthode Dwelling() est la 3ème méthode publique de l'algorithme utilisé pour l'optimisation (Method_АО_2). L'objectif de cette méthode est de mettre à jour les meilleures coordonnées et les valeurs d'aptitude de chaque particule par rapport à ses performances précédentes, mais aussi de mettre à jour l'aptitude de l'essaim et les meilleures coordonnées de l'essaim si nécessaire. La méthode est appelée après avoir obtenu les valeurs d'aptitude dans la boucle d'itération.

//——————————————————————————————————————————————————————————————————————————————
void C_AO_PSO::Dwelling ()
{
  for (int i = 0; i < swarmSize; i++)
  {
    //remember the best position for the particle
    if (p [i].ff > p [i].ffB)
    {
      p [i].ffB = p [i].ff;
      for (int k = 0; k < parameters; k++) p [i].cB [k] = p [i].c [k];
    }

    if (p [i].ff > ffB)
    {
      ffB = p [i].ff;
      for (int k = 0; k < parameters; k++) cB [k] = p [i].c [k];
    }
  }
}
//——————————————————————————————————————————————————————————————————————————————

La fonction de discrétisation du nombre de type "double" avec un pas donné dans l'intervalle spécifié.

//——————————————————————————————————————————————————————————————————————————————
// Choice in discrete space
double C_AO_PSO::SeInDiSp (double in, double inMin, double inMax, double step)
{
  if (in <= inMin) return (inMin);
  if (in >= inMax) return (inMax);
  if (step == 0.0) return (in);
  else return (inMin + step * (double)MathRound ((in - inMin) / step));
}
//——————————————————————————————————————————————————————————————————————————————

Fonction permettant d'obtenir un nombre de type "double" aléatoire dans un intervalle spécifié.

//——————————————————————————————————————————————————————————————————————————————
// Random number generator in the custom interval
double C_AO_PSO::RNDfromCI (double min, double max)
{
  if (min == max) return (min);
  double Min, Max;
  if (min > max)
  {
    Min = max;
    Max = min;
  }
  else
  {
    Min = min;
    Max = max;
  }
  return (double(Min + ((Max - Min) * (double)MathRand () / 32767.0)));
}
//——————————————————————————————————————————————————————————————————————————————

La théorie est terminée. Passons maintenant à la pratique.

Étant donné que j'utilise la même structure de construction des algorithmes que dans le 1er article de la série (et que je continuerai à le faire à l'avenir), décrite dans la figure 2, il ne sera pas difficile pour nous de relier l'algorithme au banc d'essai.

Lors de l'exécution des tests, nous verrons des animations similaires à celles présentées ci-dessous. Grâce à cela, on voit clairement comment se comporte un essaim de particules. L'essaim se comporte vraiment comme un essaim dans la nature. Nous pouvons voir sur la carte thermique de la fonction qu’ils se déplace sous la forme d'un nuage dense.

Comme précédemment, le cercle noir représente l'optimum global (max) de la fonction, et le point noir représente les meilleures coordonnées moyennes de l'algorithme de recherche obtenues au moment de l'itération en cours. Laissez-moi vous expliquer d'où viennent les valeurs moyennes. La carte thermique est bidimensionnelle en coordonnées. La fonction optimisée peut comprendre des centaines de variables (mesures). La moyenne des résultats est donc calculée en fonction des coordonnées.

  PSO sur la fonction de test Skin.

  PSO sur la fonction de test Forest.

  PSO sur la fonction de test Megacity.

Comme vous pouvez le voir dans l'animation, les tests ont montré que le PSO s'adapte assez bien à la première fonction lisse, mais seulement lorsqu'il s'agit d'optimiser 2 variables. Avec l'augmentation de la dimension de l'espace de recherche, l'efficacité de l'algorithme diminue fortement. Ceci est particulièrement visible pour la 2ème et la 3ème fonction discrète. Les résultats sont nettement moins bons que ceux de l'algorithme aléatoire décrit dans l'article précédent. Nous reviendrons sur les résultats et nous en discuterons en détail lors de la constitution d'un tableau comparatif des résultats.

En voyant comment le célèbre algorithme PSO perd honteusement au profit d'un algorithme aléatoire, on peut avoir envie de lui donner une seconde chance. Dans la section suivante, nous essayerons de modifier l'algorithme PSO.

4. Version modifiée

Pour moi, les faiblesses de PSO sont les suivantes :

1) Chacune des coordonnées est nécessairement modifiée avec une probabilité presque égale à 1. Cela signifie que chaque particule de l'essaim, à chaque itération, oscille, au mieux, quelque part dans l'extremum local de la région globale, et au pire, elle n'atteindra jamais exactement le point de l'extremum global. Cela implique une caractéristique de l'algorithme : le blocage dans des extrema locaux, c'est-à-dire une mauvaise convergence.

2) L'algorithme ne fonctionne pas bien avec les fonctions discrètes, en partie à cause du premier défaut. Les coordonnées des particules sautent aux "zones" les plus proches de la surface de la fonction à optimiser, sans pouvoir étudier en détail le voisinage d'un extremum local.

3) Faible capacité de l'algorithme à explorer de nouveaux domaines. L'essaim se concentre quelque part au même endroit sans essayer de s'échapper du "trou" local. Certains auteurs mentionnent des tentatives de mise en œuvre d'une modification à plusieurs essaims. Mais la question de l'optimisation de fonctions complexes à plusieurs variables reste ouverte, car le principe de la distance mutuelle n'est pas clair. Il faut non seulement respecter le principe du mouvement conjoint, mais aussi la possibilité d'une répulsion mutuelle. Dans le cas contraire, une telle mise en œuvre n'a aucun intérêt car les essaims individuels convergeront simplement vers une zone ou un point donné. En même temps, l'optimisation de fonctions simples à une ou deux variables est réalisée par les méthodes les plus simples et avec une excellente convergence.

Que pouvons-nous faire pour améliorer l'algorithme ?

Il est évident, mais pas nécessairement vrai, que nous devons transmettre les meilleures coordonnées individuelles aux particules à partir d'autres particules. Plus les coordonnées globales de la particule "donneuse" sont bonnes, plus la probabilité de passer les coordonnées est élevée. L'évolution de la probabilité de choisir une particule est schématisée à la figure 4. Nous générons un nombre aléatoire compris entre 0 et 1, puis transformons le nombre obtenu à l'aide d'une fonction parabolique. Nous l'adaptons ensuite à la plage des numéros de série des particules de l'essaim, comprise entre 0 et <Taille de l’Essaim>-1. Nous devons pour cela introduire un paramètre supplémentaire pour PSOm (l'algorithme modifié) : la probabilité de copier la coordonnée. Nous devons également trier l'essaim de sorte que plus la particule est bonne, plus elle est proche de l'indice 0.

Fig. 4 : Probabilité de sélection des particules décalées

Modifions légèrement la méthode ParticleMovement(). Nous générons un nombre aléatoire entre 0 et 1. Si le nombre s'avère supérieur au paramètre 'copy', nous effectuerons les opérations habituelles avec la particule décrites en détail ci-dessus. Sinon nous copierons la coordonnée d'une autre particule avec l'index choisi selon la règle illustrée graphiquement à la figure 4.

//——————————————————————————————————————————————————————————————————————————————
void C_AO_PSOm::ParticleMovement ()
{
  double rp;       //random component of particle movement
  double rg;
  double velocity;
  double posit;
  double positBest;
  double groupBest;

  for (int i = 0; i < swarmSize; i++)
  {
    for (int k = 0; k < parameters; k++)
    {
      rp = RNDfromCI (0.0, 1.0);
      rg = RNDfromCI (0.0, 1.0);

      double rC = RNDfromCI (0.0, 1.0);

      if (rC > copy)
      {
        velocity  = p [i].v  [k];
        posit     = p [i].c  [k];
        positBest = p [i].cB [k];
        groupBest = cB [k];

        p [i].v [k] = inertia * velocity + selfBoost * rp * (positBest - posit) + groupBoost * rg * (groupBest - posit);
        p [i].c [k] = posit + p [i].v [k];

        p [i].c [k] = SeInDiSp (p [i].c [k], rangeMin [k], rangeMax [k], rangeStep [k]);
      }
      else p [i].c [k] = p [GetPartcileAdress ()].cB [k];
    }
  }
}
//——————————————————————————————————————————————————————————————————————————————

La méthode Dwelling() doit également être modifiée. Ajoutons l'appel à la fonction de tri SortParticles().

//——————————————————————————————————————————————————————————————————————————————
void C_AO_PSOm::Dwelling ()
{
  for (int i = 0; i < swarmSize; i++)
  {
    //remember the best position for the particle
    if (p [i].ff > p [i].ffB)
    {
      p [i].ffB = p [i].ff;
      for (int k = 0; k < parameters; k++) p [i].cB [k] = p [i].c [k];
    }

    if (p [i].ff > ffB)
    {
      ffB = p [i].ff;
      for (int k = 0; k < parameters; k++) cB [k] = p [i].c [k];
    }
  }

  SortParticles ();
}
//——————————————————————————————————————————————————————————————————————————————

La fonction GetParticleAdress() permet de choisir l'adresse d'une particule avec une probabilité décalée vers la meilleure particule.

//——————————————————————————————————————————————————————————————————————————————
//shift of probability in the smaller party (to an index 0)
int C_AO_PSOm::GetParticleAdress ()
{
  double x = RNDfromCI (-1.0, 0.0);
  x = x * x;
  x = Scale (x, 0.0, 1.0, 0, swarmSize - 1);
  x = SeInDiSp (x, 0, swarmSize - 1, 1);
  return ((int)x);
}
//——————————————————————————————————————————————————————————————————————————————

La fonction SortParticles() est un tri à bulles conventionnel.

//——————————————————————————————————————————————————————————————————————————————
//Sorting of particles
void C_AO_PSOm::SortParticles ()
{
  //----------------------------------------------------------------------------
  int   cnt = 1;
  int   t0 = 0;
  double t1 = 0.0;
  //----------------------------------------------------------------------------

  // We will put indexes in the temporary array
  for (int i = 0; i < swarmSize; i++)
  {
    ind [i] = i;
    val [i] = p [i].ffB; //ffPop [i];
  }

  while (cnt > 0)
  {

    cnt = 0;
    for (int i = 0; i < swarmSize - 1; i++)
    {
      if (val [i] < val [i + 1])
      {
        t0 = ind [i + 1];
        t1 = val [i + 1];
        ind [i + 1] = ind [i];
        val [i + 1] = val [i];
        ind [i] = t0;
        val [i] = t1;

        cnt++;
      }
    }
  }

  // On the received indexes create the sorted temporary population
  for (int u = 0; u < swarmSize; u++) pT [u] = p [ind [u]];

  // Copy the sorted array back
  for (int u = 0; u < swarmSize; u++) p [u] = pT [u];
}
//——————————————————————————————————————————————————————————————————————————————

Fonction permettant de mettre à l'échelle un nombre d'une plage numérique à une autre.

//——————————————————————————————————————————————————————————————————————————————
double C_AO_PSOm::Scale (double In, double InMIN, double InMAX, double OutMIN, double OutMAX)
{
  if (OutMIN == OutMAX) return (OutMIN);
  if (InMIN == InMAX) return (double((OutMIN + OutMAX) / 2.0));
  else
  {
    if (In < InMIN) return (OutMIN);
    if (In > InMAX) return (OutMAX);
    return (((In - InMIN) * (OutMAX - OutMIN) / (InMAX - InMIN)) + OutMIN);
  }
}
//——————————————————————————————————————————————————————————————————————————————

5. Résultats des tests

Résumons enfin les résultats de la recherche actuelle. Malheureusement, l'algorithme PSO ne s'est pas justifié, même si j'espérerais de bons résultats. Les études ont montré sa faible convergence pour les fonctions discrètes avec ruptures et avec un grand nombre d'arguments. Une tentative de modification de l'algorithme s'est avérée infructueuse, les résultats s'avérant encore plus mauvais que ceux de la mise en œuvre classique.

L'augmentation du paramètre de copie à des valeurs proches de 0,8 permet toujours d'obtenir une convergence instantanée, mais uniquement pour la première fonction lisse des tests et avec seulement 2 arguments. Pour d'autres tests, ce paramètre ne fait qu'aggraver encore plus les résultats. L'implémentation classique de PSO n'a réussi à montrer quelque chose d'intéressant que sur la fonction Skin avec 1.000 arguments. D'autres tests se sont révélés mauvais.

Le résultat final du test est de 0,47695 pour l'algorithme classique et de 0,45144 pour l'algorithme modifié. Les résultats sont présentés dans le tableau ci-dessous. Le banc d'essai nous permet de sélectionner le nombre d'essais dans les paramètres (la valeur par défaut est de 5). Les lecteurs peuvent ainsi obtenir des résultats statistiquement plus précis en augmentant ce paramètre si la puissance de calcul le permet.

AO	Exécutions	Skin			Forest			Megacity (discrète)			Résultat final
		2 paramètres (1 F)	40 paramètres (20 F)	1.000 paramètres (500 F)	2 paramètres (1 F)	40 paramètres (20 F)	1.000 paramètres (500 F)	2 paramètres (1 F)	40 paramètres (20 F)	1.000 paramètres (500 F)
RND	1.000	0,98744	0,61852	0,49408	0,89582	0,19645	0,14042	0,77333	0,19000	0,14283	0,51254
	10.000	0,99977	0,69448	0,50188	0,98181	0,24433	0,14042	0,88000	0,20133	0,14283
PSO	1.000	0,98436	0,72243	0,65483	0,71122	0,15603	0,08727	0,53333	0,08000	0,04085	0,47695
	10.000	0,99836	0,72329	0,65483	0,97615	0,19640	0,09219	0,82667	0,10600	0,04085
PSOm	1.000	0,96678	0,64727	0,57654	0,80616	0,13388	0,06800	0,53333	0,08067	0,04211	0,45144
	10.000	0,99505	0,64986	0,57654	0,90401	0,18194	0,07104	0,74667	0,10400	0,04211

En résumé, PSO tend à conduire à une convergence rapide et prématurée aux points optimaux moyens, en plus d'une convergence lente dans la zone de recherche affinée (avec une faible capacité de recherche locale). En d'autres termes, l'algorithme est très faible et n'est pas adapté à l'optimisation des fonctions complexes, et plus encore de fonctions discrètes, en particulier de fonctions à arguments multiples.

Plusieurs approches peuvent être utilisées pour améliorer le PSO en général. La taille de l'essaim est l'un des facteurs importants. Une taille d'essaim plus importante peut augmenter la probabilité d'une convergence plus rapide et plus précise. En pratique, il existe cependant souvent un seuil du nombre d'exécutions de la fonction d'aptitude, et l'augmentation de la taille de l'essaim ne fera que réduire proportionnellement le nombre d'époques, et donc la possibilité d'évolution de l'essaim.

La deuxième approche consiste à trouver un équilibre entre l'exploration et l'exploitation. Au début des itérations, le degré élevé d'exploration donne une grande chance de trouver une solution proche de l'optimum global. Et à la fin de l'itération, un degré élevé d'exploitation donne à la particule une chance de trouver la solution la plus précise dans la zone prometteuse. La pré-optimisation avant l'étape de l'essaim par d'autres méthodes est une autre technique pouvant être utilisée pour améliorer les performances sous-jacentes du PSO, assez couramment utilisée de nos jours. L'attribution de tâches ou d'objectifs différents à chaque sous-groupe peut également accroître l'efficacité du PSO dans les problèmes multitâches. 

Une autre approche encore pour améliorer les performances du PSO consiste à établir les éléments constitutifs de l'équation de vitesse (contrôle dynamique de la vitesse). Une telle approche devrait envoyer les particules dans des directions différentes et conduire à une convergence plus rapide vers l'optimum global, mais ce n'est qu'une hypothèse.

Avantages :

1. L'algorithme est très simple et peu exigeant en termes de ressources informatiques. Le code fonctionne très rapidement, puisque l'implémentation classique ne trie même pas les tableaux.
2. L'algorithme donne de bons résultats avec les fonctions lisses. Jusqu'à présent, le PSO est le leader incontesté dans le domaine de l'optimisation des fonctions lisses à arguments multiples.
   

Inconvénients :

1. La faible qualité de l'étude de la fonction optimisée. L'algorithme ne peut pas être utilisé pour résoudre des problèmes nécessitant un ensemble de plusieurs solutions uniques.
2. Vitesse et précision de convergence faibles
3. Inadaptation à l'optimisation des fonctions discrètes
   
4. Absence quasi totale d'évolutivité