Memoire Online - Contribution à l'analyse et la synthèse des systèmes d'ordre fractionnaire par la représentation d'état

La modélisation du comportement dynamiquedes systèmesdans esquels se déroulent des phénomènes physiques â caractère diffusif, par des modèlesutilisant adérivation entière classique donne des modèles entier de dimensiontrès élevée.Lutilisation dea dérivation d'ordre non entier par contre, permet d'obtenirdesmodèlesn'utilisant que très peu de paramètres grâce â lordre nonentier Néanmoins, a richessentroduite par cet ordre de dérivation rendle calcul fractionnairetrès complexe, c'este cas notamment de l'identiification des systèmesEn effet, si pour les systèmes entiers, e problème deanon linéarité due aux pôles du modèlea été résolu par les algorithmesd'identiification bien établis tels que l'algorithme de Marquardt-Levenberg [0] et'algorithme de Levy33],63] l'ordre de dérivation rend le problème plus complexe encore.Les méthodes d'optimisation Heuristiques, telles que les algorithmes génétiques 23], 29]ou 'optimisationpar essaims de particules (OEP) [16], [33], sont d'un intérêt particulier pour résoudre cegenrede problème puisqu'elles ne nécessitent aucune information autreque a fonction optimiser elle-même. C'est dans ce cadre que sinscrit le travailprésenté dans ce présent chapitre. L'OEP est associée â l'algorithme didentiification Vector Fitting" VF)22],22] pour

développer un nouvel algorithme didentificationdes systèmes non entiersdansedomaine fréquentiel. Cet algorithme fonctionnede manièrehiérarchisée dansun niveau supérieur, en supposant connus les paramètres, lOEP permet d'optimiser l'ordre non entier et dans un niveau inférieur, l'ordre non entier étant connu, l'algorithme VectorFitting permet d'optimiser les paramètres du modèle en utilisant la forme pôlessrésidus.

Le premier paragraphe est consacré à la présentationde 'algorithme d'identiication des systèmes entiers dans le domaine fréquentiel.Dans le deuxième paragraphe on présente le principe de l'optimisation par essaim particulaire ainsi que sesdééfinitions debase.La combinaison des deux algorithmes permettant lidentiification des systèmesnon entiers dans le domaine fréquentiel est développée dans le troisième paragraphe.Les exemples numériques d'application sont également présentés dans es paragraphes3, 4 et 5.

Etant donné un système G dont le comportement fréquentiel est connu dans un intervalle de fréquences donnéLe problème est didentiifier ce système parun modèle pôle résidus de la forme :

ri et pi sont respectivement les résidus et les pôlesde G(s), d et h sont deux nombres réels h = 0 lorsque le modèle doit être propre et d = 0 lorsqu'il doit être strictement propre. n étant la dimension de G(s). L'objectif est d'identifier ces paramètres en utilisant a méthode des moindres carré même si le problème est non linéaireà causedes pôles de G(s). L'algorithme "Vector Fitting" permet alorsde résoudre ce problèmendirectement et itérativement enintroduisant deux autres fonctions de transfert

? ?

H(s) = (ó(s)G(s)) Pn i=1 s-ãj + õ + s æ

(4.2)

Àj

ó(s) Pn i=1 s-ãj + 1

lever le problème de non linéarité, on suppose également quees pôlesãi, appelés les pôles d'initialisation, sont connus au début de chaque itération.

4.1.1 Identification des paramètres de H(s) et ó(s) En substituant l'expression de ó(s) dans celle de H(s) on obtient :

_Xn
i=1

Qi
s _ ãi

( _Xn )

Ài

+ y + s æ _ G(s) = G(s) (4.3)

s _ ãi

i=1

Au début de chaque itération, Les pôles ãi étant connus (égaux à ceux calculés dans l'itération précédente)cette équation est linéairevis à visdes coefficients%i, y, æ et Ài. Elle peut donc être mise sous la formelinéaire

Fk et gk sont les k`eme lignes de F et g respectivement. Le vecteur ayant (2n+2) éléments, la résolution de l'équation (4.5) nécessite ( 2n + 2) mesures. Il faut noter également que la matrice F et le vecteur g sont complexes alors que le vecteursolution doit être réel. L'équation (4.5) doit alors être écrite sous a forme

FR et FC sont les parties réelles etimaginaires de la matrice F. gR et gC sont les parties réelles et imaginaires du vecteur g.

L'équation (4.6) devient ainsi une équation linéaire dont es coefficients sont tous réels. Néanmoins, sa résolution nécessite lutilisationde la pseudo nversion de matrice puisque F n'est plus une matrice carrée. La solution est dans ce cas donnée par :

Les deux fonctions de transfert H(s) et ó(s) étant déterminées, pour calculer les pôles pi de G(s), H(s) et ó(s) doivent être écrites sous la forme

_fln+1

%i i=1 (s--àzi)

H(s) = Lan fln

i=1 s--yi + V + s æ = i=1(s--yi)

_fln (4.8)

Ài i₌1(s--zi)

ó( s) = Lan i=1 s--yi + 1 = fln i=1(s--yi)

G(s) est alors donnée par

H(s)

G(s) = ó(s)

_fIn+1

(ó(s)G(s))

ó(s) = fIn i=1 (s - àzi)

i=1(s - zi) (4.9)

Cette équation montre que les pôles de G(s) ne sont autre que les zéros de ó(s) et que les pôles d'initialisation ãi sont simplifiés lors de la division de H(s) par ó(s). De plus, les zéros de ó(s) étant les pôles de son transfert inverse 1/ó(s), au lieu de calculer les zéros de ó(s), il est plus indiqué de calculer les pôles de 1/ó(s) qui ne sont autre que les valeurs propres de la matrice système du modèle détatcorrespondant.

xÿ =Ax+Bu
y=Cx+Du

(4.10)

A est une matrice diagonale dont les éléments sont les pôles ãi de ó(s), B un vecteur
colonne unité, C est un vecteur ligne dont les élémentssont les résidus de ó(s) et D = 1.

Pour déterminer le modèle d'état de 1/ó(s) à partir de celui de ó(s), il suffit que l'entrée u(t) de ó(s) devienne la sortie de 1/ó(s) et que la sortie y(t) de ó(s) devienne l'entrée de 1/ó(s). Ainsi, à partir de l'équation de sortie du modèle détat de ó(s) (y = Cx + Du), on a:

En remplaçant cette expression dans léquation dynamique et entenant comptedu fait que D = 1, on obtient :

A la fin de chaque itération, silors du calcul des nouveaux pôles de G(s), on trouve des pôles instables, ils doivent être transformés en pôles stablesorsqu'on souhaitemposer un modèle G(s) stable. Dans le cas des systèmes entiers, ceci peut êtreobtenu en changeant le signe de la partie réelle des pôles sans modifier eur dynamique.Danse cas des systèmes non entiers commensurables, le changement designe de apartie réelle d'un pôle dans le plan complexe s^á (a étant l'ordre commensurable compris entre 0 et 1) le rend en effet stable mais modifie également sa dynamique.Eneffet, un pôle situé dans la partie stable du demi plan droit du plan complexe donne une dynamique oscillatoire alors que s'il est situé dansle demi plan gaucheil donneraitunedynamique apériodique. Il est donc nécessaire de maintenir le pôle dans la partie stabledudemiplan droitdu plan complexe. La figure (4.1)illustre la manièrede procéder orsque es pôles sont complexes. Dans le cas où les pôles sont réelsil suffit de changerleurs signes comme danse cas entier.

Si on note les pôles instables par s1,2 instable et les pôles stables par s1,2 stable respectivement, la relation qui exprimeles pôles stables en fonctiondes pôlesnstables est donnée par:

En effet, le déplacement des pôles doit se faire de sorte à mainteniremême module des pôles et l'argument des pôles stables doit être symétrique par rapport à a droite de pente a^ð₂, de l'argument des pôles instablesUn simple calcul géométrique permet alors de déduire la relation (4.14)

Après avoir calculé les pôles de G(s), les autres paramètres (les n résidus ri et les deux paramètres d et h) sont calculés par une autre équationlinéairede la forme F è = g définie par :

Les étapes principales de l'algorithme "Vector Fitting" sont récapitulées danslorgaa nigramme de la figure (4.2)

Pour arrêter l'itération de lalgorithme, deuxcritères ddarrêt sont utilisés. Le premier limite le nombre maximal d'itérations, empéchant ainsi llalgorithme d'itérerndéfiniment. Le deuxième critère d'arrêt est le critère à minimiser constitué parlerreur entrees valeurs des réponses fréquentielles produites par e modèledentifié et cellesdonnéespar le système. Il est donné par :

G^*(sk) étant les donnés du systèmeG(sk) les données générées parle modèle identifié et Nf est le nombre de de fréquences utilisées.

L'optimisation par essaim particulaire (OEP), comme es algorithmes génétiques, est une méthode d'optimisation heuristique baséesur lasimulationducomportement collec tif des êtres vivants tels que des oiseaux ou des poissons.Cette méthode d'optimisation, inventé par l'électricien Ebenhart R.et le socio-psychologue Kennedy .en 1995, [33] s'appuie nottament sur un modèle développé par le biologiste ReynoldC.W72] permettant de simuler le déplacement d'un groupe doiseaux. Cette méthode se base sur a collaboraa tion des individus d'un même essaim en essayant de maintenirconstante adistante qui les sépare afin d'éviter de se chevaucher lorsquils changent de direction.

Cette technique est souvent décrite comme une sorted'algorithme évolutionnaire avec une population d'individus (les particules) dans laquelle, àchaque pasde temps, es "meilleurs" (selon un critère prédéfini) sont plus au moins mitésparesautres.Un autre aspect essentielle, propre à cette technique, est l'existence d'une mémoire quedoit posséé der chaque élément de l'essaim lui permettant de se souvenir de sa meilleure performance et celle transmise par ses congénères. De plus, Les individus de 'essaimtravaillent en collaboration (en s'échangeant des informations) et non pas en compétition comme dans les algorithmes génétiques par exemple.

Pour expliquer le principe de cet algorithme appliqué pour résoudreun problème de minimisation ou de maximisationconsidérons le problèmed'optimisation

L'espace de recherche représente lespacede variationdes paramètres(xi) à optimiser, il est délimité par les valeurs minimales et maximales de ces paramètres.Le nombre de paramètres à optimiser d constitue la dimension de l'espace de recherche.

Une particule, également appelée "élément de lassaim" représente une solution poo tentielle au problème à optimiser (4.15) Elle est constituée parune combinaison donnée des paramètres à optimiser (xi). Comme la particule est ammenée à évoluerelle est représentée, dans l'espace de recherche, par une position X(i,j). (i = 1, 2, · · · , M et j = 1, 2, · · · , d. i est le rang de la particule danslessaim et j le rang du paramètre x(j) qui compose le i`eme individu. M est la taille de l'essaim et d la dimension de l'espace de recherche.

A chaque paramètre x(j) est associée la vitesse dévolution v(j). La vitesse d'évolution de la i`eme particule est alors définie par V(i, j). La meilleure position déjà occupée par la i`eme particule est représentée par pbest(i) :

pbest(i, j) étant la valeur du paramètre x(j) correspondant à la meilleure position occupée par la i`eme particule. On lui associe également la valeur de safiteness F itpbest(i).

La meilleure position déjà occupée par la meilleure particule de 'essaim est représenn tée par gbest :

gbest(j) est la valeur du paramètre x(j) correspondant à la meilleure position occupée par la meilleure particule de l'essaimOn lui assicie aussi la valeur de safitnessFit gbest.

Les trois éléments fondamentaux pour calculer le déplacement d'une particule, d'une position à l'autre, sont décrites par la figure (4..3).

~ La particule se déplace selon sa vitesse propre(elle se déplace selon sonntuition) (flèche 1).

~ Elle se déplace vers la meilleure position quelle a déà occupée.On ditqu'elle a tendance à retourner vers la position de sa meilleureperformance elle se déplace selon sa propre expérience) (flèche 2).

~ Elle se déplace également versla position de la meilleure performancedéjà trouvée par une autre particule delessaim. (elle atendance a faireconfiance à 'information transmise par les autres particules) (flèche 3).

Un coefficient de confiance est alors associé à chacune de ces troisvitesses.Ainsi, a particule ne rejoint aucune des trois positions précédentes mais sedéplace vers une nouu velle position qui est la combinaisonlinéaire de ces trois positions.

La vitesse V(i,j) et la position X(i,j) de chaque paramètre sont alors mises à our à chaque itération par :

V (i, j) = w V (i, j) + c1 rand1 _pbest(i) - X(i, j) + c2 rand2 _gbest - X(i, j) X(i,j) = X(i,j) + V(i,j)

(4.21) Pour être plus précis, on devrait représenter a vitesse eta position de chaqueparamètre à l'itération k par respectivement V'

on a préféré garderla notation usuellement utilisée pour ne pas surchargeres variables.

c1 et c2 sont deux constantes d'accélération, elles caractérisent a capacité deapartii cule à chercher dans un autre endroit de lespacederecherche oubienà affiner sa recherche à l'endroit où elle se trouveEn général on choisitc1 et c2 telles que c1 +c2 <4 [16]. rand1 et rand2 sont deux nombres aléatoires compris entre 0 et 1. Les coefficients de confiance de la particule en sa propre expérience et la confiance qu'elledonne à'information transmise par les autres particules sont ainsi générés aléatoirement à chaquetération. La pondéé ration w change à chaque itérationAu début de la recherche, on lui donne une valeur assez grande pour accélérer la recherche avec des variationsde a position assez grandes (recherche approximative) Ensuite, au fur et à mesureque aparticules'approche dea meilleure solution de l'essaim, cette pondération devient plus petite afin de permettre d'affiner la recherche de la position optimale. On peut utiliser 'expression suivante pour déterminer les valeurs de cette pondération16].

iter : est le rang de l'itération actuelle. itermax : est le nombre maximum d'itération. wmax : est la valeur initiale de la pondération, on la prend généralement égale à 0.9. wmin : est la valeur finale de la pondération elle est comprise entre 0.3 et 0.4 [16]. L'organigramme de la figure (4.4), montre les étapes de lalgorithme d'optimisation par essaim particulaire.

gorithmes d'optimisationitératifsstochastiques. Linitialisation dea position et dea vitesse de chaque paramètre de chaque particuleest obtenue par

( )

X(i, j) = X(j)min + X(j)max - X(j)_mi_nrand

( ) (4.23)

V (i, j) = V _(j)min + V _(j)max - V (j)_mi_nrand

où : X(j)min ^et X(j)max sont les valeurs limites du paramètre x(j). V(j)min = 0 et V(j)max = 1. rand est un nombre aléatoire compris entre 0 et 1.

Initialement, chaque particule a sa propre vitesse et sa propre positionimitée dans l'espace de recherche. A chaqueitération, toutes les particules changent eur vitesse ete déplacent selon les équations (4.19) Certaines peuvent alors se déplacerhors de'espace de recherche. Pour éviter ce problème, on assigne à la particule sortantea valeur du point de frontière le plus proche La vitesse de la particuleconcernée est alors annulée pour l'empécher de se déplacer à la prochaine itération.En général, Lorsquea solution trouvée par l'algorithme se trouve sur la limitede 'espace de recherche, cela signiie que quelques limites des paramètres ne sont pas correctes. Le mécanisme de con~nnement des particules est donné par :

????

if X(i,j) > X(j)max = X(i,j) = X(j)max if X(i,j) <X(j)min = X(i,j) = X(j)min V(i,j) = 0

(4.24)

Dans ce paragraphe, l'optimisation par essaim particulaire est associée à 'algorithme d'identification "Vector Fitting" pour développer un nouvel algorithme d'identification des systèmes non entiers dansle domaine fréquentiel.

pour lequel il faut déterminer lordre non entier a ainsi que les paramètres ri, pi, d et h. La dimension n du modèle étant préalablement fixé

Le problème étant fortement nonlinéaire àcause des pôesqui sont au dénominateur de la fonction de transfert et de lordre nonentier qui esta puissance de'opérateur de Laplace, pour le résoudre on propose de procéder comme suit.

~ Dans une première étape, on suppose connus les paramètres du modèleri, pi, d et h), l'ordre non entier a est déterminé en utilisant l'optimisation paressaim particulaire.

~ Dans une seconde étape, a étant déterminé, à l'aide du changement de variable

dont les paramètres sont déterminés en utilisant l'algorithme "Vector Fitting"

L'algorithme d'identification global du modèle non entier basculede a première à a seconde étape et vise et versa itérativement usqu'àce que 'un desdeux critères d'arrrt est vérifié. Le premier critère étant lerreur quadratique entrees données générées pare modèle et les mesures effectuées sur lesystème.Le deuxième critère étant enombremaxii mal d'itérations. L'organigramme de la figure (4.5) lluste e principe de cet algorithme.

FIGURE 4.5: Algorithme hybride d'identification utilisant simultanément VectorFitting" et l'Optimisation par Essaim de Particules

4.4.1 Utilisation de l'algorithme "Vector Fitting" pour l'approxii mation d'un dérivateur non entier

L'approximation du dérivateur non entier par un modèle entier a fait'obbet de pluu sieurs travaux puisqu'il consitue lélément principal dans lasimulation des systèmesnon entiers. On utilise dans ce qui suit lalgorithme "Vector Fitting" pour proposer une autre méthode d'approximation du dérivateur non entier basée sur une technique d'identiicaa tion. Pour montrer l'intérêt de cette méthode, on compare ses résultats à ceux obtenus avec la méthode d'approximation CRONE. On considèrepour ce faire, edérivateurs^0.6 qui est d'abord approximé àl'aide dun transfert rationnelde dimension n = 10 dans la bande de fréquences [10^-5 10⁺⁵]. Les paramètres d'un transfert entierdu même dimension sont ensuite identifiés à partir des données générées par edérivateurs^0.6 dans la même bande de fréquences.

FIGURE 4.6: Position des pôles et zéros des modèles entiers qui approximent s^0.6.

La figure (4.6) représente les pôles et zéros des deuxtransferts rationnels.les axes0 dB sont décalés pour mieux représenter la position des pôleset zéros de chaque fonction de transfert entière). la figure (4.7) illustrees diagrammes de Bode des deuxmodèles entiers qui approximent le dérivateur s^0.6.

Ces résultats montrent que les pôles et zéros du modèle entier obtenus à 'aide de l'algorithme "Vector Fitting" sont distribués de a même manière récursive que ceux du modèle obtenu par la méthode CRONE (figure 4.6) Cette courbe montreégalement que la bande de fréquences de validité de lapproximation obtenue enutlisant 'algorithme "Vector Fitting" est pluslarge que celle utilisée par a méthode CRONE.Commee montre la figure (4.7). Lapproximationde s^0.6 par l'algorithme "Vector Fitting" estdonc meilleure.

FIGURE 4.7: Diagrammes de Bode des modèles entiers qui approximent s^0.6 dans la bande de fréquences [10_⁵ 10⁺⁵] pour n = 10.

Un vecteur contenant 500 valeurs des pulsations logarithmiquementréparties dans a bande de fréquences [10_⁶ 10⁺⁶] a été généré, pour chacune delle on calcule le module en dB et la phase en radian de Gf rac(s).

Trois valeurs différentes de n ont été choisies (n = 3, n = 5 et n = 8). La taille de l'essaim M est choisie égale à 20. La dimension de l'espace de recherche étant égal à 1, l'espace de recherche est délimité par les valeurs 0 et 1. La valeur de l'itération maximale est égale à 40.

Pour chacune des valeurs de n, l'algorithme trouve la bonne valeur de a et les bonnes valeurs des paramètres. Lerreur quadratique de a différence desmodules estégale à 3.36 10_⁶ dB.

FIGURE 4.8: Diagrammes de Bode du système (4.28) et de son modèledentifié 4.29

Pour montrer la robustesse de lalgorithme d'identification, es valeursdu module et de la phase générés par le modèle non entier (4.28) ont été perturbées eneur ajoutant des valeurs aléatoires comprises entre -30% et +30% de leurs valeurs. Les résultats obtenus sont illustrés dans la figure (4.8) lerreur quadratique est1, 12 10_² dB, l'ordre non entier est a = 0, 495 et le modèle identifié est donné par

Il faut noter que les pôles complexes, qui caractérise e comportement dynamique du système sont correctement identifiés, ce qui nest pasle casdespôles et éros entiers. ela justifie la différence entre les diagrammes de Bode en hautes fréquences.

Contribution à l'analyse et la synthèse des systèmes d'ordre fractionnaire par la représentation d'état

Chapitre 4