13.13. Algorithme de calcul « ParticleSwarmOptimization »

13.13.1. Description

Cet algorithme réalise une estimation de l’état d’un système par minimisation d’une fonctionnelle d’écart J en utilisant une méthode évolutionnaire d’essaim particulaire. C’est une méthode qui n’utilise pas les dérivées de la fonctionnelle d’écart. Elle est basée sur l’évolution d’une population (appelée « essaim ») d’états (chaque individu étant appelé une « particule » ou un insecte). Elle entre dans la même catégorie que les Algorithme de calcul « DerivativeFreeOptimization », Algorithme de calcul « DifferentialEvolution » ou Algorithme de calcul « TabuSearch ».

C’est une méthode d’optimisation mono-objectif, permettant la recherche du minimum global d’une fonctionnelle d’erreur J quelconque de type L^1, L^2 ou L^{\infty}, avec ou sans pondérations, comme décrit dans la section pour Approfondir l’estimation d’état par des méthodes d’optimisation. La fonctionnelle d’erreur par défaut est celle de moindres carrés pondérés augmentés, classiquement utilisée en assimilation de données.

Il existe diverses variantes de cet algorithme. On propose ici les formulations stables et robustes suivantes :

  • « CanonicalPSO » (Canonical Particule Swarm Optimisation, voir [ZambranoBigiarini13]), algorithme classique dit « canonique » d’essaim particulaire, robuste et définissant une référence des algorithmes d’essaims particulaires,

  • « OGCR » (Simple Particule Swarm Optimisation), algorithme simplifié d’essaim particulaire sans bornes sur les insectes ou les vitesses, déconseillé car peu robuste, mais parfois beaucoup plus rapide,

  • « SPSO-2011 » ou « SPSO-2011-AIS » (Standard Particle Swarm Optimisation 2011, voir [ZambranoBigiarini13]), algorithme de référence 2011 d’essaim particulaire, robuste, performant et défini comme une référence des algorithmes d’essaims particulaires. Cet algorithme est parfois appelé « \omega-PSO » ou « Inertia PSO » car il intègre une contribution dite d’inertie, ou encore appelé « AIS » (pour « Asynchronous Iteration Strategy ») ou « APSO » (pour « Advanced Particle Swarm Optimisation ») car il intègre la mise à jour évolutive des meilleurs éléments, conduisant à une convergence intrinsèquement améliorée de l’algorithme.

  • « SPSO-2011-SIS » (Standard Particle Swarm Optimisation 2011 with Synchronous Iteration Strategy), très similaire à l’algorithme de référence 2011 et avec une mise à jour synchrone, appelée « SIS », des particules.

  • « SPSO-2011-PSIS » (Standard Particle Swarm Optimisation 2011 with Parallel Synchronous Iteration Strategy), similaire à l’algorithme « SPSO-2011-SIS » avec mise à jour synchrone et parallélisation, appelée « PSIS », des particules.

Voici quelques suggestions pratiques pour une utilisation efficace de ces algorithmes :

  • La variante recommandée de cet algorithme est le « SPSO-2011 » même si l’algorithme « CanonicalPSO » reste par défaut le plus robuste. Dans le cas où l’évaluation de l’état peut être réalisé en parallèle, on peut utiliser l’algorithme « SPSO-2011-PSIS » même si sa convergence est parfois un peu moins performante.

  • Le nombre de particules ou d’insectes usuellement recommandé varie entre 40 et 100 selon l’algorithme, à peu près indépendamment de la dimension de l’espace des états. En général, les meilleurs performances sont obtenues pour des populations de 70 à 500 particules. Même si la valeur par défaut de ce paramètre de base provient d’une expérience étendue sur ces algorithmes, il est recommandé de l’adapter à la difficulté des problèmes traités.

  • Le nombre recommandé de générations, lors de l’évolution de la population, est souvent de l’ordre de 50, mais il peut facilement varier entre 25 et 500.

  • Le nombre maximal d’évaluation de la fonction de simulation doit usuellement être limité entre quelques milliers et quelques dizaines de milliers de fois la dimension de l’espace des états.

  • La fonctionnelle d’erreur décroît usuellement par pallier (donc avec une progression nulle de la valeur de fonctionnelle à chaque génération lorsque l’on reste dans le palier), rendant non recommandé un arrêt sur critère de décroissance de la fonction-coût. Il est normalement plus judicieux d’adapter le nombre d’itérations ou de générations pour accélérer la convergence des algorithmes.

  • Si le problème est contraint, il faut définir les bornes des variables (par la variable « Bounds »). Si le problème est totalement non contraint, il est indispensable de définir des bornes d’incrément (par la variable « BoxBounds ») pour circonscrire la recherche optimale de manière utile. De manière similaire, si le problème est partiellement contraint, il est recommandé (mais pas indispensable) de définir des bornes d’incrément. Dans le cas où ces bornes d’incréments ne sont pas définies, ce sont les bornes des variables qui seront utilisées comme bornes d’incréments.

Ces conseils sont à utiliser comme des indications expérimentales, et pas comme des prescriptions, car ils sont à apprécier ou à adapter selon la physique de chaque problème que l’on traite.

Le décompte du nombre d’évaluations de la fonction à simuler lors de cet algorithme est déterministe, à savoir le « nombre d’itérations ou de générations » multiplié par le « nombre d’individus de la population ». Avec les valeurs par défaut, il faut entre 40x50=2000 et 100*50=5000 évaluations par défaut. C’est pour cette raison que cet algorithme est usuellement intéressant lorsque la dimension de l’espace des états est grande, ou que les non-linéarités de la simulation rendent compliqué, ou invalide, l’évaluation du gradient de la fonctionnelle par approximation numérique. Mais il est aussi nécessaire que le calcul de la fonction à simuler ne soit pas trop coûteuse pour éviter une temps d’optimisation rédhibitoire.

13.13.2. Commandes requises et optionnelles

Les commandes générales requises, disponibles en édition dans l’interface graphique ou textuelle, sont les suivantes :

Background

Vecteur. La variable désigne le vecteur d’ébauche ou d’initialisation, usuellement noté \mathbf{x}^b. Sa valeur est définie comme un objet de type « Vector » ou « VectorSerie ». Sa disponibilité en sortie est conditionnée par le booléen « Stored » associé en entrée.

BackgroundError

Matrice. La variable désigne la matrice de covariance des erreurs d’ébauche, usuellement notée \mathbf{B}. Sa valeur est définie comme un objet de type « Matrix », de type « ScalarSparseMatrix », ou de type « DiagonalSparseMatrix », comme décrit en détail dans la section Conditions requises pour décrire des matrices de covariance. Sa disponibilité en sortie est conditionnée par le booléen « Stored » associé en entrée.

Observation

Liste de vecteurs. La variable désigne le vecteur d’observation utilisé en assimilation de données ou en optimisation, et usuellement noté \mathbf{y}^o. Sa valeur est définie comme un objet de type « Vector » si c’est une unique observation (temporelle ou pas) ou « VectorSerie » si c’est une succession d’observations. Sa disponibilité en sortie est conditionnée par le booléen « Stored » associé en entrée.

ObservationError

Matrice. La variable désigne la matrice de covariance des erreurs a priori d’ébauche, usuellement notée \mathbf{R}. Cette matrice est définie comme un objet de type « Matrix », de type « ScalarSparseMatrix », ou de type « DiagonalSparseMatrix », comme décrit en détail dans la section Conditions requises pour décrire des matrices de covariance. Sa disponibilité en sortie est conditionnée par le booléen « Stored » associé en entrée.

ObservationOperator

Opérateur. La variable désigne l’opérateur d’observation, usuellement noté H, qui transforme les paramètres d’entrée \mathbf{x} en résultats \mathbf{y} qui sont à comparer aux observations \mathbf{y}^o. Sa valeur est définie comme un objet de type « Function » ou de type « Matrix ». Dans le cas du type « Function », différentes formes fonctionnelles peuvent être utilisées, comme décrit dans la section Conditions requises pour les fonctions décrivant un opérateur. Si un contrôle U est inclus dans le modèle d’observation, l’opérateur doit être appliqué à une paire (X,U).

Les commandes optionnelles générales, disponibles en édition dans l’interface graphique ou textuelle, sont indiquées dans la Liste des commandes et mots-clés pour un cas d’assimilation de données ou d’optimisation. De plus, les paramètres de la commande « AlgorithmParameters » permettent d’indiquer les options particulières, décrites ci-après, de l’algorithme. On se reportera à la Description des options d’un algorithme par « AlgorithmParameters » pour le bon usage de cette commande.

Les options sont les suivantes :

Bounds

Liste de paires de valeurs réelles. Cette clé permet de définir des paires de bornes supérieure et inférieure pour chaque variable d’état optimisée. Les bornes doivent être données par une liste de liste de paires de bornes inférieure/supérieure pour chaque variable, avec une valeur None chaque fois qu’il n’y a pas de borne. Les bornes peuvent toujours être spécifiées, mais seuls les optimiseurs sous contraintes les prennent en compte.

Exemple : {"Bounds":[[2.,5.],[1.e-2,10.],[-30.,None],[None,None]]}

BoxBounds

Liste de paires de valeurs réelles. Cette clé permet de définir des paires de bornes supérieure et inférieure pour chaque incrément de variable d’état optimisée (et non pas chaque variable d’état elle-même, dont les bornes peuvent être indiquées par la variable « Bounds »). Les bornes d’incréments doivent être données par une liste de liste de paires de bornes inférieure/supérieure pour chaque incrément de variable, avec une valeur None chaque fois qu’il n’y a pas de borne. Cette clé est requise uniquement s’il n’y a pas de bornes de paramètres, et il n’y a pas de valeurs par défaut.

Exemple : {"BoxBounds":[[-0.5,0.5], [0.01,2.], [0.,None], [None,None]]}

CognitiveAcceleration

Valeur réelle. Cette clé indique le taux de rappel vers la meilleure valeur connue précédemment de l’insecte courant. C’est une valeur réelle positive. Le défaut est à peu près de 1/2+ln(2) et il est recommandé de l’adapter, plutôt en le réduisant, au cas physique qui est en traitement.

Exemple : {"CognitiveAcceleration":1.19315}

InertiaWeight

Valeur réelle. Cette clé indique la part de la vitesse de l’essaim qui est imposée à l’insecte, dite « poids de l’inertie ». C’est une valeur réelle comprise entre 0 et 1. Le défaut est de à peu près 1/(2*ln(2)) et il est recommandé de l’adapter au cas physique qui est en traitement.

Exemple : {"InertiaWeight":0.72135}

InitializationPoint

Vecteur. La variable désigne un vecteur à utiliser comme l’état initial autour duquel démarre un algorithme itératif. Par défaut, cet état initial n’a pas besoin d’être fourni et il est égal à l’ébauche \mathbf{x}^b. Sa valeur doit permettre de construire un vecteur de taille identique à l’ébauche. Dans le cas où il est fourni, il ne remplace l’ébauche que pour l’initialisation.

Exemple : {"InitializationPoint":[1, 2, 3, 4, 5]}

MaximumNumberOfFunctionEvaluations

Valeur entière. Cette clé indique le nombre maximum d’évaluations possibles de la fonctionnelle à optimiser. Le défaut est de 15000, qui est une limite arbitraire. Il est ainsi recommandé d’adapter ce paramètre aux besoins pour des problèmes réels. Pour certains optimiseurs, le nombre effectif d’évaluations à l’arrêt peut être légèrement différent de la limite à cause d’exigences de déroulement interne de l’algorithme.

Exemple : {"MaximumNumberOfFunctionEvaluations":50}

MaximumNumberOfIterations

Valeur entière. Cette clé indique le nombre maximum d’itérations interne possibles en optimisation itérative. Le défaut est 50, qui est une limite arbitraire. Il est ainsi recommandé d’adapter ce paramètre aux besoins pour des problèmes réels.

Exemple : {"MaximumNumberOfIterations":50}

NumberOfInsects

Valeur entière. Cette clé indique le nombre d’insectes ou de particules dans l’essaim. La valeur par défaut est 100, qui est une valeur par défaut usuelle pour cet algorithme.

Exemple : {"NumberOfInsects":100}

QualityCriterion

Nom prédéfini. Cette clé indique le critère de qualité, qui est minimisé pour trouver l’estimation optimale de l’état. Le défaut est le critère usuel de l’assimilation de données nommé « DA », qui est le critère de moindres carrés pondérés augmentés. Le critère possible est dans la liste suivante, dans laquelle les noms équivalents sont indiqués par un signe « <=> » : [« AugmentedWeightedLeastSquares » <=> « AWLS » <=> « DA », « WeightedLeastSquares » <=> « WLS », « LeastSquares » <=> « LS » <=> « L2 », « AbsoluteValue » <=> « L1 », « MaximumError » <=> « ME » <=> « Linf »]. On pourra se reporter à la section pour Approfondir l’estimation d’état par des méthodes d’optimisation afin de disposer de la définition détaillée de ces critères de qualité.

Exemple : {"QualityCriterion":"DA"}

SetSeed

Valeur entière. Cette clé permet de donner un nombre entier pour fixer la graine du générateur aléatoire utilisé dans l’algorithme. Par défaut, la graine est laissée non initialisée, et elle utilise ainsi l’initialisation par défaut de l’ordinateur, qui varie donc à chaque étude. Pour assurer la reproductibilité de résultats impliquant des tirages aléatoires, il est fortement conseillé d’initialiser la graine. Une valeur simple est par exemple 123456789. Il est conseillé de mettre un entier à plus de 6 ou 7 chiffres pour bien initialiser le générateur aléatoire.

Exemple : {"SetSeed":123456789}

SocialAcceleration

Valeur réelle. Cette clé indique le taux de rappel vers le meilleur insecte du voisinage de l’insecte courant, qui est par défaut l’essaim complet. C’est une valeur réelle positive. Le défaut est à peu près de 1/2+ln(2)=1.19315 et il est recommandé de l’adapter, plutôt en le réduisant, au cas physique qui est en traitement.

Exemple : {"SocialAcceleration":1.19315}

StoreSupplementaryCalculations

Liste de noms. Cette liste indique les noms des variables supplémentaires, qui peuvent être disponibles au cours du déroulement ou à la fin de l’algorithme, si elles sont initialement demandées par l’utilisateur. Leur disponibilité implique, potentiellement, des calculs ou du stockage coûteux. La valeur par défaut est donc une liste vide, aucune de ces variables n’étant calculée et stockée par défaut (sauf les variables inconditionnelles). Les noms possibles pour les variables supplémentaires sont dans la liste suivante (la description détaillée de chaque variable nommée est donnée dans la suite de cette documentation par algorithme spécifique, dans la sous-partie « Informations et variables disponibles à la fin de l’algorithme ») : [ « Analysis », « BMA », « CostFunctionJ », « CostFunctionJb », « CostFunctionJo », « CurrentIterationNumber », « CurrentState », « Innovation », « InternalCostFunctionJ », « InternalCostFunctionJb », « InternalCostFunctionJo », « InternalStates », « OMA », « OMB », « SimulatedObservationAtBackground », « SimulatedObservationAtCurrentState », « SimulatedObservationAtOptimum », ].

Exemple : {"StoreSupplementaryCalculations":["CurrentState", "Residu"]}

SwarmTopology

Nom prédéfini. Cette clé indique la manière dont les particules (ou insectes) se communiquent des informations lors de l’évolution de l’essaim particulaire. La méthode la plus classique consiste à échanger des informations entre toutes les particules (nommée « gbest » ou « FullyConnectedNeighborhood »). Mais il est souvent plus efficace d’échanger des informations sur un voisinage réduit, comme dans la méthode classique « lbest » (ou « RingNeighborhoodWithRadius1 ») échangeant des informations avec les deux particules voisines dans l’ordre de numérotation (la précédente et la suivante), ou la méthode « RingNeighborhoodWithRadius2 » échangeant avec les 4 voisins (les deux précédents et les deux suivants). Une variante de voisinage réduit consiste à échanger avec 3 voisins (méthode « AdaptativeRandomWith3Neighbors ») ou 5 voisins (méthode « AdaptativeRandomWith5Neighbors ») choisis aléatoirement (la particule pouvant être tirée plusieurs fois). La valeur par défaut est « FullyConnectedNeighborhood », et il est conseillé de la changer avec prudence en fonction des propriétés du système physique simulé. La topologie de communication possible est à choisir dans la liste suivante, dans laquelle les noms équivalents sont indiqués par un signe « <=> » : [« FullyConnectedNeighborhood » <=> « FullyConnectedNeighbourhood » <=> « gbest », « RingNeighborhoodWithRadius1 » <=> « RingNeighbourhoodWithRadius1 » <=> « lbest », « RingNeighborhoodWithRadius2 » <=> « RingNeighbourhoodWithRadius2 », « AdaptativeRandomWith3Neighbors » <=> « AdaptativeRandomWith3Neighbours » <=> « abest », « AdaptativeRandomWith5Neighbors » <=> « AdaptativeRandomWith5Neighbours »].

Exemple : {"SwarmTopology":"FullyConnectedNeighborhood"}

Variant

Nom prédéfini. Cette clé permet de choisir l’une des variantes possibles pour l’algorithme principal. La variante par défaut est la formulation « CanonicalPSO » d’origine, et les choix possibles sont « CanonicalPSO » (Canonical Particule Swarm Optimisation), « OGCR » (Simple Particule Swarm Optimisation), « SPSO-2011 » (Standard Standard Particle Swarm Optimisation 2011).

Il est conseillé d’essayer la variante « CanonicalPSO » avec une centaine de particules pour une performance robuste, et de réduire le nombre de particules à une quarantaine pour toutes les variantes autres que la formulation « CanonicalPSO » originale.

Exemple : {"Variant":"CanonicalPSO"}

VelocityClampingFactor

Valeur réelle. Cette clé indique le taux d’atténuation de la vitesse de groupe dans la mise à jour pour chaque insecte, utile pour éviter l’explosion de l’essaim, c’est-à-dire une croissance incontrôlée de la vitesse des insectes. C’est une valeur réelle comprise entre 0+ et 1. Le défaut est de 0.3.

Exemple : {"VelocityClampingFactor":0.3}

13.13.3. Informations et variables disponibles à la fin de l’algorithme

En sortie, après exécution de l’algorithme, on dispose d’informations et de variables issues du calcul. La description des Variables et informations disponibles en sortie indique la manière de les obtenir par la méthode nommée get, de la variable « ADD » du post-processing en interface graphique, ou du cas en interface textuelle. Les variables d’entrée, mises à disposition de l’utilisateur en sortie pour faciliter l’écriture des procédures de post-processing, sont décrites dans l”Inventaire des informations potentiellement disponibles en sortie.

Sorties permanentes (non conditionnelles)

Les sorties non conditionnelles de l’algorithme sont les suivantes :

Analysis

Liste de vecteurs. Chaque élément de cette variable est un état optimal \mathbf{x}^* en optimisation, une interpolation ou une analyse \mathbf{x}^a en assimilation de données.

Exemple : xa = ADD.get("Analysis")[-1]

CostFunctionJ

Liste de valeurs. Chaque élément est une valeur de fonctionnelle d’écart J choisie.

Exemple : J = ADD.get("CostFunctionJ")[:]

CostFunctionJb

Liste de valeurs. Chaque élément est une valeur de fonctionnelle d’écart J^b, c’est-à-dire de la partie écart à l’ébauche. Si cette partie n’existe pas dans la fonctionnelle, sa valeur est nulle.

Exemple : Jb = ADD.get("CostFunctionJb")[:]

CostFunctionJo

Liste de valeurs. Chaque élément est une valeur de fonctionnelle d’écart J^o, c’est-à-dire de la partie écart à l’observation.

Exemple : Jo = ADD.get("CostFunctionJo")[:]

Ensemble des sorties à la demande (conditionnelles ou non)

L’ensemble des sorties (conditionnelles ou non) de l’algorithme, classées par ordre alphabétique, est le suivant :

Analysis

Liste de vecteurs. Chaque élément de cette variable est un état optimal \mathbf{x}^* en optimisation, une interpolation ou une analyse \mathbf{x}^a en assimilation de données.

Exemple : xa = ADD.get("Analysis")[-1]

BMA

Liste de vecteurs. Chaque élément est un vecteur d’écart entre l’ébauche et l’état optimal.

Exemple : bma = ADD.get("BMA")[-1]

CostFunctionJ

Liste de valeurs. Chaque élément est une valeur de fonctionnelle d’écart J choisie.

Exemple : J = ADD.get("CostFunctionJ")[:]

CostFunctionJb

Liste de valeurs. Chaque élément est une valeur de fonctionnelle d’écart J^b, c’est-à-dire de la partie écart à l’ébauche. Si cette partie n’existe pas dans la fonctionnelle, sa valeur est nulle.

Exemple : Jb = ADD.get("CostFunctionJb")[:]

CostFunctionJo

Liste de valeurs. Chaque élément est une valeur de fonctionnelle d’écart J^o, c’est-à-dire de la partie écart à l’observation.

Exemple : Jo = ADD.get("CostFunctionJo")[:]

CurrentIterationNumber

Liste d’entiers. Chaque élément est l’index d’itération courant au cours du déroulement itératif de l’algorithme utilisé. Il y a une valeur d’index d’itération par pas d’assimilation correspondant à un état observé.

Exemple : cin = ADD.get("CurrentIterationNumber")[-1]

CurrentState

Liste de vecteurs. Chaque élément est un vecteur d’état courant utilisé au cours du déroulement itératif de l’algorithme utilisé.

Exemple : xs = ADD.get("CurrentState")[:]

Innovation

Liste de vecteurs. Chaque élément est un vecteur d’innovation, qui est en statique l’écart de l’optimum à l’ébauche, et en dynamique l’incrément d’évolution.

Exemple : d = ADD.get("Innovation")[-1]

OMA

Liste de vecteurs. Chaque élément est un vecteur d’écart entre l’observation et l’état optimal dans l’espace des observations.

Exemple : oma = ADD.get("OMA")[-1]

OMB

Liste de vecteurs. Chaque élément est un vecteur d’écart entre l’observation et l’état d’ébauche dans l’espace des observations.

Exemple : omb = ADD.get("OMB")[-1]

SimulatedObservationAtBackground

Liste de vecteurs. Chaque élément est un vecteur d’observation simulé par l’opérateur d’observation à partir de l’ébauche \mathbf{x}^b. C’est la prévision à partir de l’ébauche, elle est parfois appelée « Dry ».

Exemple : hxb = ADD.get("SimulatedObservationAtBackground")[-1]

SimulatedObservationAtCurrentState

Liste de vecteurs. Chaque élément est un vecteur d’observation simulé par l’opérateur d’observation à partir de l’état courant, c’est-à-dire dans l’espace des observations.

Exemple : hxs = ADD.get("SimulatedObservationAtCurrentState")[-1]

SimulatedObservationAtOptimum

Liste de vecteurs. Chaque élément est un vecteur d’observation obtenu par l’opérateur d’observation à partir de la simulation d’analyse ou d’état optimal \mathbf{x}^a. C’est l’observation de la prévision à partir de l’analyse ou de l’état optimal, et elle est parfois appelée « Forecast ».

Exemple : hxa = ADD.get("SimulatedObservationAtOptimum")[-1]

13.13.4. Exemples d’utilisation en Python (TUI)

Voici un ou des exemples très simple d’usage de l’algorithme proposé et de ses paramètres, écrit en [DocR] Interface textuelle pour l’utilisateur (TUI/API). De plus, lorsque c’est possible, les informations indiquées en entrée permettent aussi de définir un cas équivalent en interface graphique [DocR] Interface graphique pour l’utilisateur (GUI/EFICAS).

Cet exemple décrit le recalage des paramètres \mathbf{x} d’un modèle d’observation H quadratique. Ce modèle est représenté ici comme une fonction nommée QuadFunction. Cette fonction accepte en entrée le vecteur de coefficients \mathbf{x}, et fournit en sortie le vecteur \mathbf{y} d’évaluation du modèle quadratique aux points de contrôle internes prédéfinis dans le modèle. Le calage s’effectue sur la base d’un jeu initial de coefficients (état d’ébauche désigné par Xb dans l’exemple), et avec l’information \mathbf{y}^o (désignée par Yobs dans l’exemple) de 5 mesures obtenues à ces mêmes points de contrôle internes. On se place en expériences jumelles (voir Pour tester une chaîne d’assimilation de données : les expériences jumelles) et les mesures sont parfaites. On privilégie les observations au détriment de l’ébauche par l’indication d’une très importante variance d’erreur d’ébauche, ici de 10^{6}.

L’ajustement s’effectue en affichant des résultats intermédiaires lors de l’optimisation itérative.

# -*- coding: utf-8 -*-
#
from numpy import array, ravel
def QuadFunction( coefficients ):
    """
    Simulation quadratique aux points x : y = a x^2 + b x + c
    """
    a, b, c = list(ravel(coefficients))
    x_points = (-5, 0, 1, 3, 10)
    y_points = []
    for x in x_points:
        y_points.append( a*x*x + b*x + c )
    return array(y_points)
#
Xb   = array([1., 1., 1.])
Yobs = array([57, 2, 3, 17, 192])
#
NumberOfInsects = 40
#
print("Résolution du problème de calage")
print("--------------------------------")
print("")
from adao import adaoBuilder
case = adaoBuilder.New()
case.setBackground( Vector = Xb, Stored=True )
case.setBackgroundError( ScalarSparseMatrix = 1.e6 )
case.setObservation( Vector = Yobs, Stored=True )
case.setObservationError( ScalarSparseMatrix = 1. )
case.setObservationOperator( OneFunction = QuadFunction )
case.setAlgorithmParameters(
    Algorithm='ParticleSwarmOptimization',
    Parameters={
        'NumberOfInsects':NumberOfInsects,
        'MaximumNumberOfIterations': 20,
        'StoreSupplementaryCalculations': [
            'CurrentState',
            ],
        'Bounds':[[0,5],[-2,2],[0,5]],
        'SetSeed':123456789,
        },
    )
case.setObserver(
    Info="  État intermédiaire en itération courante :",
    Template='ValuePrinter',
    Variable='CurrentState',
    )
case.execute()
print("")
#
#-------------------------------------------------------------------------------
#
print("Calage de %i coefficients pour une forme quadratique 1D sur %i mesures"%(
    len(case.get('Background')),
    len(case.get('Observation')),
    ))
print("--------------------------------------------------------------------")
print("")
print("Vecteur d'observation.............:", ravel(case.get('Observation')))
print("État d'ébauche a priori...........:", ravel(case.get('Background')))
print("")
print("Coefficients théoriques attendus..:", ravel((2,-1,2)))
print("")
print("Nombre d'itérations...............:", len(case.get('CurrentState')))
print("Nombre de simulations.............:", NumberOfInsects*len(case.get('CurrentState')))
print("Coefficients résultants du calage.:", ravel(case.get('Analysis')[-1]))
#
Xa = case.get('Analysis')[-1]
import matplotlib.pyplot as plt
plt.rcParams['figure.figsize'] = (10, 4)
#
plt.figure()
plt.plot((-5,0,1,3,10),QuadFunction(Xb),'b-',label="Simulation à l'ébauche")
plt.plot((-5,0,1,3,10),Yobs,            'kX',label='Observation',markersize=10)
plt.plot((-5,0,1,3,10),QuadFunction(Xa),'r-',label="Simulation à l'optimum")
plt.legend()
plt.title('Calage de coefficients', fontweight='bold')
plt.xlabel('Coordonnée arbitraire')
plt.ylabel('Observations')
plt.savefig("simple_ParticleSwarmOptimization1.png")

Le résultat de son exécution est le suivant :

Résolution du problème de calage
--------------------------------

  État intermédiaire en itération courante : [ 1.76770856 -1.2054263   1.22625259]
  État intermédiaire en itération courante : [ 2.04776518 -1.17449716  3.14493347]
  État intermédiaire en itération courante : [ 2.04776518 -1.17449716  3.14493347]
  État intermédiaire en itération courante : [ 2.01933291 -1.          3.15162067]
  État intermédiaire en itération courante : [ 1.96384202 -0.74855119  3.40642058]
  État intermédiaire en itération courante : [ 1.96384202 -0.74855119  3.40642058]
  État intermédiaire en itération courante : [ 1.96384202 -0.74855119  3.40642058]
  État intermédiaire en itération courante : [ 1.96384202 -0.74855119  3.40642058]
  État intermédiaire en itération courante : [ 1.95417745 -0.73191939  3.14451887]
  État intermédiaire en itération courante : [ 1.96217646 -0.7883895   2.91127919]
  État intermédiaire en itération courante : [ 1.97610485 -1.00254825  2.89582746]
  État intermédiaire en itération courante : [ 2.0007262  -1.06443275  2.69825603]
  État intermédiaire en itération courante : [ 1.9934285  -1.02432071  2.39823319]
  État intermédiaire en itération courante : [ 1.9934285  -1.02432071  2.39823319]
  État intermédiaire en itération courante : [ 1.9942533  -0.99256953  2.30174702]
  État intermédiaire en itération courante : [ 1.9942533  -0.99256953  2.30174702]
  État intermédiaire en itération courante : [ 1.99742923 -0.99796085  2.1278678 ]
  État intermédiaire en itération courante : [ 1.99742923 -0.99796085  2.1278678 ]
  État intermédiaire en itération courante : [ 1.99742923 -0.99796085  2.1278678 ]
  État intermédiaire en itération courante : [ 1.99742923 -0.99796085  2.1278678 ]
  État intermédiaire en itération courante : [ 2.00166149 -1.0012696   2.02137857]

Calage de 3 coefficients pour une forme quadratique 1D sur 5 mesures
--------------------------------------------------------------------

Vecteur d'observation.............: [ 57.   2.   3.  17. 192.]
État d'ébauche a priori...........: [1. 1. 1.]

Coefficients théoriques attendus..: [ 2 -1  2]

Nombre d'itérations...............: 21
Nombre de simulations.............: 840
Coefficients résultants du calage.: [ 2.00166149 -1.0012696   2.02137857]

Les graphiques illustrant le résultat de son exécution sont les suivants :

_images/simple_ParticleSwarmOptimization1.png