Identification et commande des systèmes non linéaires : Utilisation des modèles de type NARMA Brahim Tlili To cite this version: Brahim Tlili. Identification et commande des systèmes non linéaires : Utilisation des modèles de type NARMA. Automatic. Ecole Nationale d’Ingénieurs de Tunis, 2008. French. <tel00589735v1> HAL Id: tel-00589735 https://tel.archives-ouvertes.fr/tel-00589735v1 Submitted on 1 May 2011 (v1), last revised 9 Jan 2011 (v2) HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers. L’archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d’enseignement et de recherche français ou étrangers, des laboratoires publics ou privés. U NIVERSITE DE T UNIS EL M ANAR ECOLE NATIONALE D’INGENIEURS DE TUNIS THESE Pr é se n t é e pou r obt e n ir le Diplôme de Doctorat S PECIALITE G ENIE E LECTRIQUE Pa r Tlili Brahim ___________________________________________________________________________ IDENTIFICATION ET COMMANDE DES SYSTEMES NON LINEAIRES : UTILISATION DES MODELES DE TYPE NARMA ___________________________________________________________________________ Soutenue le 29 juillet 2008 devant le jury composé de : Mme. MOUFIDA KSOURI Présidente MM. RIDHA BEN ABDENNOUR Rapporteur FARHAT FNAIECH Rapporteur MEKKI KSOURI Examinateur, directeur de thèse FAOUZI BOUANI Examinateur TABLE DES MATIERES Introduction générale Chapitre 1 Contexte et état de l’art 8 11 1.1. Introduction ............................................................................................................ 12 1.2. Sur l’identification et la modélisation .................................................................... 12 1.3. Sur la commande prédictive à base de modèle ...................................................... 18 1.3.1. Un aperçu historique de la commande prédictive ....................................... 18 1.3.2. Principe de la commande prédictive .......................................................... 19 1.4. Conclusion.............................................................................................................. 22 Chapitre 2 ldentification des systèmes non linéaires monovariables 23 2.1. Introduction ............................................................................................................ 24 2.2. Le modèle NARMA ............................................................................................... 25 2.3. Identification du modèle NARMA......................................................................... 26 2.3.1. Méthode basée sur l’algorithme de Kortmann ............................................. 27 2.3.2. Méthode basée sur les Algorithmes Génétiques Binaires ............................ 38 2.3.3. Méthode basée sur les Réseaux de Neurones Artificiels............................. 45 2.4. Résultats de simulation........................................................................................... 50 2.5. Conclusion.............................................................................................................. 61 Chapitre 3 La commande prédictive non linéaire 63 3.1. Introduction ............................................................................................................ 64 3.2. Calcul du prédicteur à j - pas ................................................................................ 65 3.3. Le critère de performance ..................................................................................... 66 3.4. Optimisation du critère........................................................................................... 67 1 3.4.1. Généralités sur les méthodes d’optimisation................................................ 67 3.4.2. La méthode de Nelder-Mead ....................................................................... 71 3.4.3. La méthode de Rosenbrock .......................................................................... 76 3.4.4. Améliorations des performances des méthodes d’optimisation................... 79 3.5. La loi de commande .............................................................................................. 86 3.6. Résultats de simulation .......................................................................................... 87 3.7. Conclusion............................................................................................................ 100 Chapitre 4 Identification et commande prédictive des systèmes non linéaires multivariables 101 4.1. Introduction .......................................................................................................... 102 4.2. Le modèle NARMA (cas multivariable).............................................................. 103 4.3. Identification du modèle NARMA....................................................................... 105 4.4. Commande prédictive multivariable .................................................................... 109 4.4.1. Calcul du prédicteur .................................................................................. 109 4.4.2. Le critère de performance .......................................................................... 109 4.4.3. Optimisation du critère............................................................................... 110 4.5. Résultats de simulation......................................................................................... 113 4.6. Conclusion............................................................................................................ 124 Conclusion générale 125 Bibliographie 127 Annexes 136 2 ABREVIATIONS MCR Méthode des Moindres Carrées Récursifs Méthode AGB Méthode d’identification basée sur l’Algorithme Génétique Binaire Méthode K Méthode d’identification basée sur l’algorithme de Kortmann Méthode RN Méthode d’identification basée sur les Réseaux de Neurones Artificiels Méthode NM Méthode d’optimisation de Nelder-Mead Méthode ROS Méthode d’optimisation de Rosenbrock AIC Akaike Information Criterion BIC Bayes’s Information Criterion CFON Constrained First Objective Next FPE Final Prediction Error GPC Generalized Predictive Control MBC Model Based Control MBPC Model Based Predictive Control NARMA Nonlinear AutoRegressive Moving Average NMPC Nonlinear Model Predictive Control OVF Overal F_test 3 NOTATIONS u Commande appliquée au système y Sortie du système yc Consigne ŷ Sortie prédite à l’aide du modèle q Degré de non linéarité du modèle Hc Horizon de commande Hp Horizon de prédiction sur la sortie λ Coefficient de pondération ΔU Vecteur des incréments de commande J, J1, J2, J3, J4 Critères de performances n1 Nombre de sorties du système multivariable m1 Nombre d’entrées du système multivariable U Vecteur des entrées Y Vecteur des sorties Φ Θ̂ T Vecteur d’observation Vecteur des paramètres transposé yi Valeur moyenne de la sortie i λj Coefficient de pondération sur la commande (pour l’entrée j) Hcj Horizon de commande sur l’entrée j Hpj Horizon de prédiction sur la sortie j R²ajust Coefficient de corrélation ajusté 2 Rmult Coefficient de corrélation multiple 2 Rtot Coefficient de corrélation total PFi Critères d’information partiels NM-V0 et ROS-V0 Versions originales de Nelder-Mead et de Rosenbrock NM-V1, NM-V2 et NM-V3 Versions de Nelder-Mead avec améliorations ROS-V1, ROS-V2 et ROS-V3 Versions de Rosenbrock avec améliorations 4 LISTE DES FIGURES Figure 1-1 : Modèle de Wiener Figure 1-2 : Modèle de Hammerstein Figure 1-3 : Modèle de Wiener-Hammerstein (LNL) Figure 1-4 : Modèle de Hammerstein-Wiener (NLN) Figure 1-5 : Représentation temporelle du principe de la commande prédictive Figure 1-6 : Schéma de principe d’une commande prédictive à base de modèle Figure 2-1 : Organigramme de l’algorithme de régression par pas échelonné modifié Figure 2-2 : Organigramme de l’identification du modèle NARMA avec l’algorithme génétique binaire Figure 2-3 : Schéma descriptif d’un réseau de neurone non bouclé à une seule couche cachée Figure 2-4 : Réseau de neurones à une seule couche cachée Figure 2-5 : Structure du superviseur basé sur l’algorithme génétique Figure 2-6 : Les évolutions de la séquence d’entrée, de la sortie du système et des sorties des modèles considérés (système 1) Figure 2-7 : Les évolutions de la séquence d’entrée, de la sortie du système et des sorties des modèles considérés (système 2) Figure 2-8 : Les évolutions des sorties (système et modèles) Figure 2-9(a) : Evolutions des sorties (système et modèles) et des erreurs de validation (système 2, Rv=2%) Figure 2-9(b) : Evolutions des sorties (système et modèles) et des erreurs de validation (système 2, Rv=10%) Figure 3-1 : Schéma de principe du prédicteur Figure 3-2 : Les méthodes d’optimisation les plus importantes Figure 3-3 : Les mouvements possibles dans la méthode du polytope de Nelder-Mead dans ℜ 2 Figure 3-4 : Initialisation des algorithmes d’optimisation versions (NM-V1 et NM-V2 ou ROS-V1 et ROS-V2) Figure 3-5 : Initialisation des algorithmes d’optimisation versions (NM-V3 ou ROS-V3) 5 Figure 3-6 : Initialisation des algorithmes d’optimisation (approche CFON) dans ℜ Figure 3-7 : Initialisation des algorithmes d’optimisation versions (NM-V1 et NM-V2 ou ROS-V1 et ROS-V2) Figure 3-8 : Initialisation des algorithmes d’optimisation versions (NM-V3 ou ROS-V3) Figure 3-9 : Sortie du système, commande et incréments de la commande (NM-V0) Figure 3-10 : Sortie du système, commande et incréments de la commande (ROS-V0) Figure 3-11 : Sortie du système, commande et incréments de la commande (NM-V1) Figure 3-12 : Sortie du système, commande et incréments de la commande (ROS-V1) Figure 3-13 : Sortie du système, commande et incréments de la commande (NM-V2) Figure 3-14 : Sortie du système, commande et incréments de la commande (ROS-V2) Figure 3-15 : Sortie du système, commande et incréments de la commande (NM-V3) Figure 3-16 : Sortie du système, commande et incréments de la commande (ROS-V3) Figure 3-17 : Sortie système, commande et ses incréments avec | ∆u | ≤ 0.1 (NM-V3) Figure 3-18 : Sortie système, commande et ses incréments avec | ∆u | ≤ 0.1 (ROS-V3) Figure 3-19 : Sortie système (S2), commande et ses incréments (méthode NM-V3) Figure 3-20 : Sortie système (S2), commande et ses incréments (méthode ROS-V3) Figure 3-21 : Sortie système (S2), commande et ses incréments (méthode NM-V3) avec bruit Figure 3-22 : Sortie système (S2), commande et ses incréments (méthode ROS-V3) avec bruit Figure 4-1 : Réseau de neurones à une seule couche cachée (cas d’un système multivariable) Figure 4-2 : Structure du superviseur basé sur l'algorithme génétique Figure 4-3(a) : Procédé de régulation du niveau d’eau à trois réservoirs interconnectés Figure 4-3(b) : Schéma du système à trois réservoirs (système multivariable 3) Figure 4-4 : Evolutions des séquences d’entrées (u1, u2), des sorties (système et modèle) et des erreurs de modélisation Figure 4-5 : Sortie y1, commande et incréments de commande (NM-V3) Figure 4-6 : Sortie y2, commande et incréments de commande (NM-V3) Figure 4-7 : Sorties y1, commande et incréments de commande (ROS-V3) Figure 4-8 : Sortie y2, commande et incréments de commande (ROS-V3) 6 LISTE DES TABLEAUX Tableau 2-1 : Paramètres des algorithmes Tableau 2-2 : Structures de l’équation simulée et des modèles identifiés (système 1) Tableau 2-3 : Paramètres du modèle NARMA pour une variance du bruit égale à 7.1 10-3 (système 1) Tableau 2-4 : Les valeurs de R²mult pour différentes valeurs de la variance du bruit (système 1) Tableau 2-5 : Structures des modèles identifiés (système 2) Tableau 2-6 : Variance des erreurs de validation Tableau 2-7 : Evolution de R²mult en fonction des paramètres m, n et q (système 3) Tableau 2-8 : Les valeurs des paramètres estimés, de R²mult et de temps de calcul (système 3) Tableau 2-9 : Evolution de R²mult en fonction des paramètres m, n et q (système 4) Tableau 2-10 : Les valeurs de R²mult et du temps de calcul (système 4) Tableau 3-1 : Résultats de l’optimisation avec versions originales (méthode NM et ROS) Tableau 3-2 : Résultats de l’optimisation avec réinitialisation (méthode NM) Tableau 3-3 : Résultats de l’optimisation avec réinitialisation (méthode ROS) Tableau 3-4 : Critères de performances Tableau 3-5 : Temps de simulation selon la méthode d'initialisation Tableau 3-6 : Critères de performances (système S2) Tableau 4-1 : Paramètres de réglage des algorithmes Tableau 4-2 : Coefficients de la fonction d’activation polynomiale pour les trois systèmes identifiés Tableau 4-3 : Termes et paramètres estimés pour les trois systèmes multivariables Tableau 4-4 : Valeurs de la somme des erreurs quadratiques normalisées entre le modèle trouvé et le modèle rectifié (système 3) Tableau 4-5 : Les valeurs de R²tot et du temps de calcul (système multivariable) Tableau 4-6 : Paramètres du modèle nominal et du modèle incertain Tableau 4-7 : Valeurs des critères de performances et du temps de calcul moyen (système multivariable 1) 7 Chapitre 1 Introduction générale Dans ce travail, on s’intéresse à la commande des systèmes non linéaires. L’approche utilisée est basée sur la commande prédictive à base de modèle MBPC (Model Based Predictive Control). En effet, ce type de commande présente plusieurs avantages. Comme illustré dans les travaux de Clarke [Clarke, 1988], la commande prédictive s’avère être une structure suffisamment complète proposée pour résoudre un problème de commande très général. Elle peut également conduire à un système asservi stable pour des paramètres de réglage donnés. De plus, cette stratégie a montré son efficacité, sa flexibilité et son succès dans des applications industrielles, même pour des systèmes à faible période d’échantillonnage. Cependant, la mise en œuvre de la commande prédictive à base de modèle nécessite : - La définition d’un modèle numérique qui peut être défini comme un ensemble d’équations mathématiques, permettant de représenter au mieux le comportement du système réel afin de prédire le comportement futur de celui-ci. Cette particularité permet de classer la commande prédictive dans la grande famille des commandes à base de modèles MBC (Model Based Control). - Le calcul d’une séquence de commandes futures qui est obtenue par la minimisation d’un critère de performance quadratique, sur un horizon fini, comportant un terme lié aux erreurs de prédiction futures (écarts entre la sortie prédite du système et la consigne) et un terme lié à la commande. Dans la pratique plusieurs types de modèles sont utilisés pour la synthèse de la loi de commande, parmi lesquels les modèles de type entrée-sortie. Dans la famille des modèles de type entrée-sortie, on distingue plusieurs types tels que le modèle de Volterra [Boyd, 1985], [Hazem, 2005], [Hernando, 1988], le modèle de Wiener [Westwick, 1996] [Wigren, 1994], [Zhu, 1999], le modèle de Hammerstein [Bai, 2008], [Bai, 2003], [Bai, 1990], [Greblicki, 2000], [Narendre, 1966], les modèles cascades Wiener-Hammerstein [Bai, 2002], [Bauer, 2002], [Hunter, 1986], [Zhu, 2002] et le modèle NARMA (Nonlinear AutoRegressive Moving Average) [Billings, 1981], [Leontaritis, 1985a]. L’utilisation d’un modèle de type 8 Chapitre 1 NARMA permet de réduire d’une manière considérable le nombre de termes, et par conséquent le nombre de paramètres du modèle. Ceci, conduit à la réduction du temps de calcul de la loi de commande comparé aux régulateurs prédictifs basés sur les autres types de modèles. Ces avantages nous ont poussés à l’exploitation de ce modèle dans le calcul de la loi de commande. Dans ce travail, nous présentons d’abord la méthode de Kortmann [Kortmann, 1988], puis nous proposons deux autres méthodes basées sur les approches non conventionnelles pour l’identification du modèle NARMA. La loi de commande est obtenue par minimisation d’un critère de performance non convexe. La minimisation du critère de performance est généralement effectuée par une approche classique tel que la méthode du gradient, qui nécessite le calcul de la dérivée de ce dernier [Filali, 2002], [Kenneth, 1999], [Yonghong, 1996]. L’emploi de cette méthode a souvent été couronné de succès, puisqu’elle dispose de nombreuses qualités (preuve de convergence, faible nombre d'évaluations nécessaires, faible dépendance vis-à-vis de la dimension du problème, prise en compte des contraintes). Néanmoins, plusieurs difficultés limitent son usage. Pour certains problèmes, la fonction de coût est non différentiable par nature ou présente des discontinuités (par exemple les problèmes de type Min Max), il s’en suit que le calcul du gradient est délicat (différences finies, coûteux, entaché d’erreur, mise en oeuvre délicate pour des problèmes complexes ou présence de minima locaux). Dès lors, cette approche ne peut pas être employée sans aménagement spécifique, l’utilisation des méthodes sans gradient se révèle alors plus adéquate. Ce mémoire comporte quatre chapitres, il est organisé comme suit : Le premier chapitre est un chapitre d’introduction qui présente la problématique générale de la commande prédictive à base de modèle. Nous présentons dans un premier temps les modèles de type entrée-sortie les plus connus, en particulier le modèle NARMA. Ensuite, nous nous intéressons à la stratégie et au principe de la commande prédictive non linéaire. Le deuxième chapitre sera consacré à l’identification paramétrique et structurelle du modèle NARMA. Trois approches sont utilisées pour l’identification du modèle de type NARMA, la première exploite l’algorithme de régression par pas échelonné, proposé par Kortmann et Unbehauen [Kortmann, 1988]. La deuxième méthode est basée sur l’algorithme génétique avec sa représentation binaire, et la troisième approche constitue une combinaison 9 Chapitre 1 des réseaux de neurones artificiels à fonction d’activation polynomiale avec l’algorithme génétique sous sa représentation réelle. Dans le troisième chapitre, nous nous intéresserons à la synthèse d’un régulateur prédictif sous contraintes qui est basé sur le modèle NARMA. La minimisation du critère de performance est effectuée avec deux méthodes qui n’utilisent pas la dérivée. La première méthode est basée sur l’algorithme de Nelder-Mead. La deuxième approche utilise la méthode de Rosenbrock qui est également une méthode d’optimisation numérique. Une fonction de pénalité a été employée pour le traitement des contraintes. Le quatrième chapitre traite les systèmes multivariables non linéaires. Dans ce cadre, nous proposons l’extension de la troisième approche utilisée dans le chapitre 2, pour l’identification du modèle NARMA dans le cas multivariable. Cette approche a été validée par une application pratique sur un procédé réel, qui présente un système à trois réservoirs d’eau interconnectés. Ensuite, nous proposons l’extension des méthodes d’optimisation présentées dans le chapitre 3, pour traiter la commande prédictive non linéaire des systèmes multivariables. Une conclusion générale ainsi que les perspectives de ces travaux clôturent ce rapport de synthèse. 10 Chapitre 1 Chapitre 1 Contexte et état de l’art ___________________________________________________________________ Chapitre 1 Contexte et état de l’art 1.1. Introduction ............................................................................................................ 12 1.2. Sur l’identification et la modélisation .................................................................... 12 1.3. Sur la commande prédictive à base de modèle ...................................................... 18 1.3.1. Un aperçu historique de la commande prédictive ....................................... 18 1.3.2. Principe de la commande prédictive .......................................................... 19 1.4. Conclusion.............................................................................................................. 22 ______________________________________________________________________ 11 Chapitre 1 Chapitre 1 Contexte et état de l’art 1.1. Introduction L’objectif de ce chapitre est de situer ce travail dans son contexte scientifique qu’est la commande des systèmes non linéaires. Dans ce cadre, la commande prédictive non linéaire, en utilisant les modèles de type NARMA, est développée. La stratégie de commande à base de modèle (MBC : Model Based Control) nécessite un modèle qui décrit le comportement dynamique du système. L’utilisation d’un modèle linéaire ou linéarisé autour d’un point de fonctionnement est souhaitable lorsqu’il s’agit de l’identification des systèmes simples, puisque ce type de modèle possède peu de paramètres tout en restant facile à déterminer. Malheureusement la complexité des systèmes industriels limite l’utilisation des modèles linéaires et plusieurs approches de modélisation des systèmes non linéaires ont été développées. Dans ce travail, on s’intéressera à la commande prédictive non linéaires sous contraintes en utilisant les modèles de type NARMA. 1.2. Sur l’identification et la modélisation Dans la pratique, deux démarches sont possibles pour l’obtention du modèle d’un système. La première approche est analytique et la deuxième est expérimentale [Ben Abdennour, 2001]. La modélisation analytique utilise les équations physico-chimiques qui régissent le comportement dynamique du système. La modélisation expérimentale traite le système comme une boîte noire et utilise les mesures expérimentales pour extraire le modèle. Bien que la modélisation analytique offre une description précise facilitant la compréhension du comportement des procédés, elle est souvent complexe et peu utilisée pour la commande. Les automaticiens utilisent surtout la modélisation expérimentale car elle présente une grande souplesse lors de la synthèse des régulateurs. Cette dernière modélisation comporte les modèles de type entrée-sortie, appelés aussi modèles de fichiers [Ben Abdennour, 2001], [Landau, 1988]. Dans [Giannakis, 2001] les auteurs ont donné une liste bibliographique, regroupant plus de 1400 références, qui porte sur l’identification des systèmes non linéaires. Cette liste présente la majorité des travaux publiés dans la période 1979-2000. 12 Chapitre 1 Les modèles non linéaires de types entrée/sortie existent sous plusieurs formes, parmi lesquels on distingue le modèle de Volterra [Hazem, 2004], [Hazem, 2005], [Hernando, 1988], le modèle de Hammerstein, le modèle de Wiener [Haber, 1990] et le modèle NARMA [Billings, 1981], [Leontaritis, 1985a]. a- Modèle de Volterra : Le modèle de Volterra appartient à la classe des modèles polynomiaux non récursifs. La relation d’entrée-sortie qui lui est associée est : y (k ) = ∑ ∞ i =1 ∑ ∑ ∞ j1 = 0 ∞ j2 = 0 ... ∑ hi ( j1 , ..., j N ) ∏ u (k − jl ) ∞ N j N =0 l =1 (1-1) où u, y et hi sont respectivement l’entrée, la sortie et le noyau de Volterra d’ordre N. Les jl, l = 1,…,N sont des entiers. Un noyau d’ordre N peut être vu comme une généralisation à l’ordre N de la notion de réponse impulsionnelle bien connue dans le cas des systèmes linéaires. En appelant M la mémoire du système, on peut considérer le modèle tronqué, à temps discret, à l’ordre N suivant : y (k ) = ∑ ∑ ∑ N M −1 M −1 i =1 j1 = 0 j2 = 0 ∑ h ( j , ... , j ) ∏ u (k − j ) M −1 ... jN =0 N i 1 N l =1 (1-2) l b- Modèles de Wiener et de Hammerstein Ces modèles sont obtenus en mettant en cascade un modèle dynamique linéaire et un bloc non linéaire statique. Lorsque la non-linéarité statique suit le bloc dynamique linéaire, le modèle est dit modèle de Wiener (figure 1-1); dans le cas contraire, on parle de modèle de Hammerstein (figure 1-2). Modèle de Wiener u Bloc linéaire dynamique w Non linéarité statique h y g(.) Figure 1-1 : Modèle de Wiener 13 Chapitre 1 Si la non-linéarité statique est de forme polynomiale d’ordre N et le bloc linéaire dynamique représente une réponse impulsionelle finie FIR (Finite Impulse Response), la relation d’entrée-sortie de ce modèle sera comme suit : w( k ) = ∑ h( j1 ) u (k − j1 ) M −1 (1-3) j1 = 0 et y (k ) = g [w(k )] = ∑ ci w i (k ) N (1-4) i =0 Modèle de Hammerstein Le modèle de Hammerstein peut être vu comme étant le dual du modèle de Wiener composé des mêmes blocs mais placés dans un ordre inversé. On a : u w Non linéarité statique Bloc linéaire dynamique y h g(.) Figure 1-2 : Modèle de Hammerstein Comme dans le cas du modèle de Wiener, si g(.) est un polynôme et h est une FIR, alors la relation de la sortie du modèle sera comme suit : w(k ) = g [u (k )] = ∑ ci u i (k ) N (1-5) i =0 et y (k ) = ∑ h( j ) w(k − j ) M −1 j1 = 0 1 1 (1-6) 14 Chapitre 1 Modèles cascades A partir des combinaisons des modèles de Wiener et de Hammerstein, on distingue d’autres types de modèles. Le modèle Wiener-Hammerstein (LNL, figure 1-3) est défini par un bloc dynamique linéaire en cascade avec un bloc non linéaire statique suivi par un autre bloc dynamique linéaire et le modèle Hammerstein-Wiener (NLN, figure 1-4) est défini par un bloc non linéaire statique en cascade avec un bloc dynamique linéaire suivi par un autre bloc non linéaire statique. u Bloc linéaire dynamique w1 h1 Non linéarité statique w2 Bloc linéaire dynamique y h2 g(.) Figure 1-3 : Modèle de Wiener-Hammerstein (LNL) u Non linéarité statique g1(.) w1 Bloc linéaire dynamique w2 y Non linéarité statique h g2(.) Figure 1-4 : Modèle de Hammerstein-Wiener (NLN) 15 Chapitre 1 c- Modèle NARMA La structure dynamique d’un modèle NARMA opérant dans un environnement déterministe, est régie par l’équation mathématique discrète suivante [Billings, 1981], [Kortmann, 1988] : yˆ (k ) = y + ∑ bi u (k − i ) + ∑∑ bij u ( k − i )u (k − j ) m m i =1 m i =1 j =i + ... + ∑ ...∑∑ bi...l u (k − i)...u (k − l ) m m i =1 m v = p l =v + ∑ ai y (k − i ) + ∑∑ aij y (k − i ) y (k − j ) n n i =1 n i =1 j =i + ... + ∑ ...∑∑ ai...l y (k − i )... y (k − l ) + ∑∑ cij u (k − i) y (k − j ) n n i =1 n v = p l =v m n (1-7) i =1 j =1 + ... + ∑ ...∑∑ ci...l u (k − i )...u (k − v) y (k − l ) m m i =1 n v =1 l =1 + ... + ∑∑ ...∑ cij ...l u (k − i ) y (k − j )... y (k − l ). m n i =1 j =1 y n l =1 représente la valeur moyenne. Les paramètres m et n représentent, respectivement, l’ordre de la régression sur l’entrée u(k) et l’ordre de la régression sur la sortie y(k). Comme on peut le constater, avec la relation ci-dessus, le modèle NARMA est un modèle non linéaire qui présente l’avantage d’être linéaire par rapport à ses paramètres, la non linéarité s’exprimant par rapport à l’entrée et à la sortie du système. L’implication pratique de cette propriété est que plusieurs des résultats obtenus pour des problèmes d’estimation paramétrique linéaire peuvent être repris pour ce type de modèles, moyennant quelques aménagements appropriés. Dans [Hernando, 1988], [Hazem, 2004] et [Hazem, 2005], les auteurs ont proposé des méthodes pour l’identification des modèles de Volterra. Cependant, ce type de modèle nécessite un nombre élevé des paramètres et par conséquent le temps de calcul de la loi de commande sera élevé. Le modèle NARMA fournit une représentation unifiée pour une large classe des systèmes non linéaires [Leontaritis, 1985a]. Billings et Leontaritis [Billings, 1981] 16 Chapitre 1 ont prouvé que plusieurs modèles bien connus tels que le modèle Hammerstein, de Wiener et les modèles bilinéaires sont des cas spéciaux du modèle NARMA. Malgré le nombre important de termes existants dans la forme générale du modèle NARMA, ce dernier nécessite un nombre limité parmi les anciennes mesures et les couples croisés (entrée-entrée, entrée-sortie et sortie-sortie) dans l’expression finale de la sortie [Liu, 2003]. Plusieurs approches ont été développées pour la détermination des modèles de type NARMA. L’algorithme de régression par pas échelonné, qui est basé sur des tests statistiques, assure l’estimation des paramètres ainsi que la sélection des termes qui forment le modèle NARMA [Kortmann, 1988]. L’exploitation des réseaux de neurones artificiels à fonction d’activation polynomiale a été également développée dans les travaux de [Ki, 1997]. Cependant, dans ce papier aucune approche n’a été proposée pour la détermination des coefficients de la fonction d’activation des neurones. Dans [Sheng, 2001] et [Sheng, 2003], les auteurs ont utilisé la géométrie affine et la minimisation de la distance «hyper surface» pour l’estimation des paramètres du modèle NARMA. Cette approche qui est basée sur la combinaison de l’algorithme OPS (Optimal Parameter Search) avec la méthode TLS (Total Least Squares), permet d’estimer les paramètres du modèle NARMA même en présence d’un bruit significatif. Cependant, cette approche nécessite un temps de calcul élevé et la solution analytique ne peut pas être calculée pour les systèmes de dimension supérieure à trois [Sheng, 2003]. Dans [Ruano, 2003], l’optimisation multi objectif basée sur les algorithmes génétiques MOGA (Multi Objective Genetic Algorithm) a été exploitée pour la détermination du modèle NARMA. Dans ce cas, l’algorithme MOGA a été utilisé pour réaliser un compromis entre plusieurs critères de performances. Ces critères sont exprimés en fonction de l’ordre du modèle, de l’erreur de modélisation et du nombre de termes du modèle. Dans le cas linéaire, l’ordre des régressions sur l’entrée et sur la sortie du modèle ARIMA a été également déterminé à l’aide des algorithmes génétiques à codage binaire [Ong, 2005]. Les modèles NARMA qui sont non linéaires par rapport aux anciennes valeurs de la sortie mais linéaires par rapport à la valeur actuelle de l’entrée ont été également développés dans [Narendra, 1997] et [Adetona, 2004]. Ces modèles, obtenus en utilisant un développement en série de Taylor des modèles neuronaux, sont utiles pour la synthèse des régulateurs. En effet, la loi de commande peut être obtenue par minimisation d’un critère de performance convexe. 17 Chapitre 1 1.3. Sur la commande prédictive à base de modèle 1.3.1. Un aperçu historique de la commande prédictive La commande prédictive (MPC) est née d'un besoin réel dans le monde industriel. Un besoin de système de régulation capable de performances plus élevées que les contrôleurs classiques, tel que le régulateur de type PID, tout en respectant des contraintes de fonctionnement et de production toujours plus sévères. Historiquement, cette approche a commencé à donner ses premiers résultats théoriques et pratiques à la fin des années 1970. Sa mise en œuvre pratique dans l’industrie a vu le jour dans le domaine pétrolier et pétrochimique, notamment avec les travaux pionnier de Richalet en 1976. Les publications de Richalet et al. [Richalet, 1976] [Richalet, 1978] présentaient la commande prédictive dite heuristique MPHC (Model Predictive Heuristic Control), qui fut connue plus tard sous le nom de commande algorithmique MAC (Model Algorithmic Control) et qui utilise la réponse impulsionnelle du système, tronquée aux N premiers échantillons [Camacho, 1999]. Cette approche c'est vite étendu à d'autres industries grâce à ses succès incontestables dans l'industrie pétrolière sérieusement éprouvée par des raisons économiques. En effet, le but économique recherché par des commandes évoluées, et par MPC en particulier, et de réduire les coûts de production en ramenant le point de régulation (setpoint), à des niveaux proches des contraintes de fonctionnement et/ou de production. Ceci est obtenue en réduisant la variance donnée par un contrôleur classique, tel que PID, en utilisant une commande évoluée, e.g. MPC, le point de régulation est ramené vers un niveau plus bas engendrant un coup de production moindre. En 1980 apparaît la commande matricielle dynamique DMC (Dynamic Matrix Control) présentée par Cutler et Ramaker [Cutler, 1980], qui utilise explicitement le modèle de la réponse indicielle du système pour prédire l’effet sur la sortie des commandes futures. Cellesci étaient calculées par la minimisation de l’erreur prédite qui était répétée à chaque période d’échantillonnage avec les dernières mesures fournies par le processus. Ces formulations étaient heuristiques et algorithmiques et tiraient parti du potentiel croissant des ordinateurs de l’époque. La commande GPC (Generalized Predictive Control) développée par Clarke et al. en 1987 [Clarke, 1987], qui applique des idées de la commande GMV (Generalized Minimum Variance ) [Clarke, 1979] est sans doute la plus populaire actuellement. Dans les travaux de Morari [Morari, 1994], la commande MPC a également été formulée dans l’espace d’état. 18 Chapitre 1 Ceci permet non seulement d’utiliser des théorèmes bien connus de la théorie de la représentation d’état, mais facilite aussi la généralisation de la commande à des cas plus complexes tels que les systèmes avec des perturbations stochastiques et du bruit dans les variables mesurées. Le principe de l’horizon fuyant, l’une des idées centrales de la commande MPC, fut quant à lui proposé en 1963 par Propoi [Propoi, 1963] dans le cadre du « retour optimal en boucle ouverte », et a été largement repris ensuite dans les années soixante-dix. 1.3.2. Principe de la commande prédictive La commande prédictive MPC (Model Predictive Control) est une méthode relativement récente, elle n’a connu un réel essor dans l’industrie que depuis le milieu des années 80. Actuellement, cette approche représente l’une des techniques de commandes les plus avancées en automatique. En effet, sa formulation intègre des concepts tirés de la commande optimale, de la commande stochastique, de la commande de systèmes à temps morts, de la commande avec modèle interne, de la commande multivariable et prend en compte les références futures lorsqu’elles sont connues d’avance. En outre, des contraintes sur la commande et sur la sortie peuvent être considérées en pratique avec ce type de commande. Cette méthode se repose sur les idées suivantes : - utilisation d’un modèle pour prédire les sorties du procédé à des instants futurs [Borne, 1997] [Clarke, 1989], - calcul de la séquence des commandes à appliquer au système de façon à minimiser un critère à horizon fini portant sur l’écart entre la sortie prédite et la sortie future désirée, - à chaque instant d’échantillonnage, l’horizon de prédiction est déplacé vers le futur, et seule la première valeur des commandes calculées est effectivement appliquée au système (notion d’horizon fuyant). La stratégie de commande prédictive MPC, appelée aussi commande à horizon glissant ou fuyant, est similaire à la stratégie utilisée pour le pilotage d’un avion. Le pilote connaît la trajectoire de référence désirée sur un horizon de commande fini (son champ visuel), et en prenant en compte les caractéristiques de son avion, il décide quelles actions (accélérer, décélérer, monter ou abaisser en altitude) doit prendre afin de suivre la trajectoire désirée. 19 Chapitre 1 Seule la première action est exécutée à chaque instant, et la procédure est répétée à nouveau pour les prochaines actions. Le principe général de MPC, illustré par la figure 1-5 dans le domaine temporel, se caractérise par les éléments suivants : - à l’aide du modèle du processus, on prédit à chaque instant k les sorties futures, sur un horizon déterminé de taille Hp (horizon de prédiction), - calcul de la suite de commandes à appliquer au système, de façon à minimiser un critère quadratique sur un horizon fini. Le critère quadratique dépend des erreurs de prédiction futures qui représentent l’écart entre la sortie prédite du système et la trajectoire de référence future et de l’énergie de la commande. Cette optimisation a pour objectif de garder le processus aussi proche que possible de la trajectoire de référence. Ainsi, une séquence de commandes futures sera générée par le calculateur. Des hypothèses peuvent être faites sur la structure de la loi de commande future, cette dernière sera constante à partir d’un instant donné Hc (horizon de commande), - parmi la séquence de commande calculée juste le premier élément u(k) sera appliqué à l’entrée du système, tous les autres éléments de la séquence sont omis, - répétition de la procédure complète de calcul à chaque période d’échantillonnage, selon le principe de l’horizon fuyant. Consigne yc ŷ Sortie prédite sortie y k Passé k+1 commandes futures u k+Hc-1 k+Hp temps Futur Figure 1-5 : Représentation temporelle du principe de la commande prédictive 20 Chapitre 1 La figure 1-6 représente un schéma bloc simplifié qui caractérise la commande prédictive à base de modèle. Le pointillé correspond aux calculs effectués par le calculateur. Avec yc(k) consigne. y(k) sortie du processus. yˆ ( k ) sortie prédite à l’aide du modèle. u(k) commande appliquée au système. CNA : Convertisseur Numérique Analogique. CAN : Convertisseur Analogique Numérique. yc(k), …,yc(k+Hp) Régulateur u (k) u(t) CNA Procédé Modèle du processus y(t) y(k) CAN yˆ ( k ), L , yˆ ( k + Hp ) Figure 1-6 : Schéma de principe d’une commande prédictive à base de modèle Aujourd’hui, la commande prédictive représente l’une des techniques modernes de commande parmi les plus utilisées dans la pratique industrielle [Boucher, 1996]. Cette technique, présente un certain nombre d’avantages par rapport aux autres méthodes de commande. Parmi lesquels, on distingue : - son principe très intuitif et le réglage relativement facile de ses paramètres qui la rendent accessible même aux personnes ayant des connaissances limitées en automatique, - elle peut être utilisée pour commander une grande variété de processus, ceux avec des dynamiques simples ou complexes, par exemple les systèmes à retard, à phases non minimales ou instables, - elle traite, également, les systèmes multivariables, 21 Chapitre 1 - elle est capable intrinsèquement de compenser les retards et les temps morts, - le traitement de contraintes sur le système à commander peut être inclus systématiquement dans la définition du correcteur, - elle est très utile lorsque les consignes ou trajectoires à suivre sont connues à l’avance (ce qui est le cas dans certains processus industriels). La plupart des algorithmes de commande prédictive reposent sur des modèles internes linéaires. La raison de cette approche linéaire se fonde dans la maîtrise complète de la théorie des systèmes linéaires. En effet, un modèle linéaire utilisé avec une approche prédictive donne un résultat analytique, c'est-à-dire, une loi de commande linéaire sous forme d'une équation mathématique. Un tel résultat permet le calcul de la commande future en des temps très courts, ceci avait une importance non négligeable quand les calculateurs n'offraient pas la rapidité de calcul désirée. Toutefois, les procédés à contrôler se révéleront de plus en plus non linéaires. Les méthodes linéaires arrivent donc à leurs limites et l’utilisation des modèles non linéaires, et par conséquent la commande prédictive non linéaire NMPC (Nonlinear Model Predictive Control) commence à s’imposer. Grâce à ses concepts intuitifs et aux bons résultats obtenus, la commande prédictive a été implantée dans un grand nombre d’applications industrielles. Ces applications industrielles ont toutes un dénominateur commun : la connaissance de la trajectoire à suivre par le système dans le futur, au moins sur un certain horizon. En 1997, Qin et Badgwell [Qin, 1997] recenseront plus de 2000 applications industrielles commandées avec MPC, dont 600 procédés commandés avec NMPC. 1.4. Conclusion L’objet de ce chapitre a été de situer la problématique de la mise en place d’un correcteur prédictif à base de modèle décrivant le type de modèle utilisé, à savoir le modèle NARMA. Nous avons ainsi présenté les différents modèles (Volterra, Wiener, Hammerstein) et nous avons donné le principe de la commande prédictive. L’objet du chapitre suivant est la présentation de la structure du modèle NARMA ainsi que la présentation des méthodes heuristiques proposées pour l’identification de ce type de modèle, sur lequel se base notre commande prédictive. 22 Chapitre 2 Chapitre 2 Identification des systèmes non linéaires monovariables __________________________________________________________________________ Chapitre 2 ldentification des systèmes non linéaires monovariables 2.1. Introduction ............................................................................................................ 24 2.2. Le modèle NARMA .............................................................................................. 25 2.3. Identification du modèle NARMA......................................................................... 26 2.3.1. Méthode basée sur l’algorithme de Kortmann ............................................. 27 2.3.2. Méthode basée sur les Algorithmes Génétiques Binaires ............................ 38 2.3.3. Méthode basée sur les Réseaux de Neurones Artificiels............................. 44 2.4. Résultats de simulation........................................................................................... 49 2.5. Conclusion.............................................................................................................. 60 ___________________________________________________________________________ 23 Chapitre 2 Chapitre 2 Identification des systèmes non linéaires monovariables 2.1. Introduction L’utilisation d’un modèle de type NARMA permet non seulement de réduire le nombre de termes, et par conséquent le nombre de paramètres du modèle mais il permet aussi la modélisation des systèmes avec un ordre relativement faible [Liu, 2003]. Dans ce chapitre, on s’intéressera à l’identification paramétrique et structurelle des modèles de type NARMA des systèmes monovariables. D’abord, on présentera l’algorithme de régression par pas échelonné [Kortmann, 1988]. Ensuite, on proposera deux nouvelles méthodes pour la détermination du modèle NARMA. La première méthode utilise l’algorithme génétique avec sa représentation binaire. Dans ce cas, les individus de la population forment les termes du modèle alors que les paramètres sont estimés avec la méthode des moindres carrées récursifs. Des critères sont utilisés afin d’assurer la convergence de la population d’une génération à une autre vers le vecteur d’observation souhaitable. La seconde approche est basée sur la combinaison du réseau de neurones artificiels à fonction d’activation polynomiale avec l’algorithme génétique sous sa représentation réelle. Dans ce cas, les individus de la population de l’algorithme génétique représentent les coefficients du polynôme de la fonction d’activation. Les paramètres du modèle NARMA sont alors estimés à partir des pondérations des connexions du réseau neuronal après la fin de la phase d’apprentissage. Le chapitre est organisé comme suit : Le deuxième paragraphe sera consacré pour présenter la structure du modèle NARMA. Le troisième paragraphe est réservé au développement des trois approches utilisées pour l’estimation structurelle et paramétrique du modèle NARMA. Dans le quatrième paragraphe, on présente les résultats de simulation. La conclusion est donnée dans le dernier paragraphe. 24 Chapitre 2 2.2. Le modèle NARMA On considère la classe des systèmes non linéaires monovariables qui peuvent être modélisés par la relation suivante : y (k ) = g ( X (k )) + e(k ) avec : (2-1) X (k ) = [ y (k − 1),..., y (k − n), u (k − 1),...u (k − m)] et g une fonction non linéaire supposée inconnue. Les paramètres m et n représentent, respectivement, l’ordre de la régression sur l’entrée u(k) et l’ordre de la régression sur la sortie y(k), e(k) est un bruit blanc de moyenne nulle et de variance finie. Dans ce travail, l’objectif de la modélisation consiste à caractériser le comportement du système par un modèle général de type NARMA. En effet, la sortie du modèle général yˆ (k ) (équation 1-7) peut s’écrire sous la forme vectorielle suivante [Kortmann, 1988] : ˆ T Φ (k ) yˆ ( k ) = Θ (2-2) avec Θ̂ T et Φ (k ) désignent, respectivement, le vecteur des paramètres et le vecteur d’observation. = [y θ1 θ2 ... θ r ]T (2-3) Φ ( k ) = [1 v1 v2 ... v r ]T (2-4) ˆ Θ où v1 (k ) = u(k-1 ), ..., vm (k ) = u (k − m) , vm+1 (k ) = y(k-1 ) ,…, vn + m(k) = y(k- n) , (2-5) vn+ m +1(k) = u 2(k-1 ) vn+ m+ 2(k) = u(k-1 ) ⋅ u(k- 2 ) ,…, , vr(k) = y q(k-n) Le paramètre q représente le degré de non linéarité du modèle. Les paramètres m, n et q caractérisent le modèle NARMA. Puisque la sortie du modèle est linéaire par rapport aux paramètres, on peut alors utiliser la méthode des Moindres Carrées Récursifs (MCR) pour 25 Chapitre 2 estimer ces paramètres [Landau, 1988]. Le nombre de termes intervenant dans l’expression de la sortie du modèle yˆ (k ) est donné par la relation suivante [Kortmann, 1988] : r= (n + m + q)! −1 q! (n + m)! (2-6) Par conséquent, le nombre de combinaisons possibles qu’on peut avoir est : w = 2r −1 (2-7) Pour m, n et q donnés, il existe un nombre élevé de combinaisons possibles. Par exemple pour un cas simple (m=n=q=2), on aura 14 termes et par conséquent 16383 combinaisons. Il est donc nécessaire d’utiliser une procédure permettant la sélection d’un modèle parmi l’ensemble des modèles possibles qui représente le comportement du système réel. Si on désigne par vi (k ) le ième terme du modèle considéré pour (n = m = q = 2) , alors les 14 termes possibles seront comme suit : v 1 (k ) = u (k − 1 ) v 2 (k ) = u(k - 2) v3 (k ) = y(k - 1) v 8 (k ) = u (k-1) y (k-1) v 9 (k ) = u (k-1) y (k − 2) v10 (k ) = u (k-2) y (k-1) v 4 (k ) = y(k - 2) v11 (k ) = u(k - 2) ⋅ y(k - 2) v6 (k ) = u(k - 1) ⋅ u(k - 2) v13 (k ) = y(k - 1) ⋅ y(k - 2) v5 (k ) = u 2 (k - 1) v7 (k ) = u 2 (k - 2) (2-8) v12 (k ) = y 2 (k - 1) v14 (k ) = y 2 (k - 2) 2.3. Identification du modèle NARMA La détermination du modèle d’un système revient à déterminer les termes qui forment ce modèle ainsi que les paramètres correspondants. On présente par la suite trois méthodes différentes pour l’identification des modèles de type NARMA. Ces méthodes exploitent hors ligne un fichier formé par les mesures de l’entrée et la sortie du système : {y(k), u(k), k=1, …, N}, N est le nombre total des mesures. 26 Chapitre 2 2.3.1. Méthode basée sur l’algorithme de Kortmann Cette méthode exploite l’algorithme de régression par pas échelonné modifié, proposé initialement par Kortmann et Unbehauen [Kortmann, 1988]. Cet algorithme, basé sur les corrélations des mesures entrée-sortie, permet à la fois de déterminer les termes et d’estimer les paramètres correspondants d’un modèle. Il comporte plusieurs critères d’informations ainsi que des tests statistiques, afin de donner des informations sur la qualité du modèle identifié [Akaike, 1974], [Labrousse, 1977] et [Pukkila, 1988]. Cet algorithme comporte huit étapes, ces étapes sont résumées comme suit : Etape 1 : Initialisation de l’algorithme Etape 2 : Calcul des coefficients de corrélation normalisés de tous les termes du vecteur d’observation avec la sortie du système. Etape 3 : Sélection d’un terme à ajouter à l’expression de la sortie du modèle. Etape 4 : Estimation des paramètres à l’aide de la méthode de (MCR) ; et calcul des critères de performances (OVF, R², FPE, AIC, LILC, BIC). Etape 5 : Réalisation d’un test qui permet d’accepter provisoirement le dernier terme sélectionné, de le rejeter ou de quitter la procédure. Etape 6 : Calcul des critères d’informations partiels (PFi, PEi, AICi, LILCi, BICi). Etape 7 : Réalisation d’un test qui permet la sélection des termes à garder dans l’expression du modèle actuel et des termes qui seront rejetés. Etape 8 : Calcul des coefficients de corrélation partiels normalisés et retour à l’étape 3. L’algorithme de régression par pas échelonné modifié suppose que les paramètres m, n et q du modèle NARMA sont connus. Ceci permet de calculer, dans l’étape 2, les coefficients de corrélation normalisés de tous les termes du vecteur d’observation. Ensuite, l’analyse d’un certain nombre de critères de performances permet de juger la qualité du modèle trouvé et par conséquent, l’arrêt de la procédure de recherche du modèle. 27 Chapitre 2 Pour expliquer le fonctionnement de cet algorithme, on présente par la suite les tests statistiques et les critères d’informations utilisés dans les étapes 4 et 6. 2.3.1.1. Tests statistiques et critères d’informations [Akaike, 1974], [Haber, 1990] A. Tests statistiques a. Critère de l’erreur de prédiction finale Ce critère noté FPE (Final Prediction Error), est défini par : FPE N + n = V ( Θˆ ) N − n (2-9) avec N est le nombre d’observations, Θ̂ est le vecteur de paramètres de dimension n et V ( Θˆ ) est l’estimation de la variance du résidu. ˆ)= 1 V (Θ N ∑ε N k =1 2 (k ) (2-10) avec ε (k) est l’erreur de prédiction ( ε ( k ) = y ( k ) − yˆ ( k ) ). Dans notre cas, n = r+1. b. Le coefficient de corrélation total. ∑ε N Il est défini par : R 2 tot =1− 2 k =1 N ∑ k =1 (k ) (2-11) 2 y (k ) avec N est le nombre de mesures. 2 La valeur de Rtot peut être utilisée pour juger la qualité du modèle obtenu. En effet, 2 lorsque la valeur de Rtot est très proche de l’unité, alors le modèle estimé caractérise convenablement le comportement du système considéré. 28 Chapitre 2 c. Le coefficient de corrélation multiple. ∑ [ yˆ (k) N = 2 R mult Ce critère est défini par : k =1 N ∑ [y(k) k =1 - y] 2 100 % - y] (2-12) 2 2 Rmult mesure la proportion de la variation totale de la valeur moyenne y par rapport à la sortie du modèle. ( d. Le coefficient de corrélation ajusté. 2 = 1 − 1 − R mult 2 R ajust Ce critère est défini par : ) NN −− 1v (2-13) 2 2 , prend en considération le critère Rmult v : désigne le nombre de termes du modèle. Rajust tout en pénalisant les modèles qui possèdent un nombre élevé de termes. B. Critères d’information [Landau, 1988], [Sheng, 2003]. La plupart des critères d’information utilisés pour déterminer le nombre de paramètres d’un modèle sont exprimés par : = N Ln ( Γ (v) σ 2 ) + v g(N) (2-14) avec : - σ 2 : l’estimé du maximum de vraisemblance ou son approximation par la variance du résidu qui est obtenu en utilisant l’estimateur des moindres carrés récursifs. σ 2 est défini par : ∑ N σ²= k =1 ε ∑ ( y (k ) − N 2 N (k ) = k =1 N yˆ ( k )) 2 (2-15) - N est le nombre d’observations - v est le nombre de paramètres estimés. - g (N) : est un terme de pénalité. 29 Chapitre 2 g (N) augmente quand le nombre de paramètres dans le modèle augmente, par contre N Ln ( σ 2 ) a tendance à décroître quand le nombre de paramètres croît. Le modèle optimal est choisi pour la valeur de v minimisant le critère Γ (v) . a. Le critère d’information d’Akaike (AIC) [Akaike, 1974] Le critère d’information d’Akaike (Akaike Information Criterion) AIC est défini en prenant g (N) = 2 dans la relation (2-14) : AIC (v) = N Ln (σ 2 ) + 2 v (2-16) b. Le critère d’information de Bayes (BIC) : Le critère d’information de Bayes (Bayes’s Information Criterion) BIC est défini lorsque g(N) = Ln (N ) dans la relation (2-14) : BIC (v) = N Ln (σ 2 ) + v Ln(N) (2-17) c. La Loi du Critère Logarithmique Itératif ( test LILC ) La Loi du Critère Logarithmique Itératif (Khinchin’s Law of Iterated Logarithm Criterion) LILC est définie en fixant g(N) = cLn[Ln(N)] dans la relation (2-14) : LILC (v) = N Ln (σ 2 ) + v c Ln[ Ln(N)] (2-18) où c est une constante positive (c ≥ 2 ) . 2.3.1.2. L’algorithme de régression par pas échelonné On présente par la suite les différentes étapes de l’algorithme de régression par pas échelonné modifié. L’algorithme, décrit par l’organigramme de la figure 2-1, comporte huit étapes [Kortmann, 1988]. 30 Chapitre 2 Initialisation des paramètres de l’algorithme : m, n, q, Xα, R2, v=1, j=0 Calcul des coefficients de corrélation - Sélection du terme le plus corrélé avec la sortie y(k) - Incrémentation de v : v=v+1 - Estimation paramétrique - Calcul de ε(k), F_test, R2, AIC, BIC et LILC Test 1 Test 2 v=v-1 Arrêt j=2 Sélection du 2eme terme le mieux corrélé avec la sortie y(k) Calcul de F_test partiel et des critères d’informations partiels j=j+1 Décision Tous les termes retenus Dernier terme rejeté Calcul des coefficients de corrélation de tous les termes Figure 2-1 : Organigramme de l’algorithme de régression par pas échelonné modifié Etape 1 (Initialisation): Au cours de cette étapes, on spécifie les valeurs maximales de l’ordre de la régression sur l’entrée (m) et l’ordre de la régression sur la sortie (n) ainsi que le degré de non linéarité (q). Ensuite, on calcule le nombre de termes possibles r donné par l’équation (2-6) et on forme le 31 Chapitre 2 vecteur d’observation Φ (k ) . On initialise également le nombre de paramètres : v = 1 (valeur moyenne y ) et on spécifie le seuil de signification α et on donne la valeur X α correspondante (en utilisant les abaques de Fisher donnés dans l’annexe1). Etape 2 On calcule les coefficients de corrélation normalisés de tous les termes du vecteur d’observation avec la sortie y (k ) ∑ v (k) . y(k) N ρi = N k =1 i ∑ v (k) ∑ y(k) N - 2 ⎡⎡ N 2 ⎤ ⎤ ⎡N ⎢⎢ N ∑ vi (k) - ⎢∑ vi (k)⎥ ⎥ ⎢⎣⎢⎣ k =1 ⎦ ⎥⎦ ⎣ k =1 k =1 N i k =1 ⎡ ⎢ N ∑ y 2(k) ⎢⎣ k =1 N ⎤ ⎡ - ⎢∑ y(k)⎥ ⎦ ⎣ k =1 N 2 ⎤⎤ ⎥⎥ ⎥⎦ ⎥⎦ i = [1,...,r ] (2-19) avec N le nombre d’observation et vi (k) les termes du modèle général donnés par l’équation (2-5). Etape 3 Au niveau de cette étape, on détermine le plus grand coefficient de corrélation ( ρ opt ) permettant de sélectionner le terme optimal à ajouter au modèle. ρ opt = max ( ρ i ), vopt = vi Ensuite, on incrémente le nombre des paramètres de 1, ( v = v + 1 ). Exemple : On considère un modèle NARMA avec n = m = q = 2 , on a 14 termes possibles {v1 (k) ,..., v14 (k) } donnés par la relation (2-8). Si ρ12 est le coefficient de corrélation maximal. Le terme sélectionné par l’algorithme est v opt (k) = v12 (k) = y 2 (k-1 ) (2-20) Le modèle actuel est alors : yˆ (k) = Θˆ T Φ(k) (2-21) 32 Chapitre 2 avec : Φ(k) = [ 1 y²(k-1 ) ]T (2-22) et Θ T = [ y θ12 ] (2-23) Etape 4 Estimer les paramètres (valeur moyenne incluse) en utilisant l’algorithme des moindres carrés récursifs et calculer les résidus : ε(k) = y(k) - yˆ (k) ; k = 1, ... , N (2-24) Déterminer le F_Test global ( Overal F_test ) OVF [Bariani, 1988] : OVF = 2 ⎤ ⎡ N ⎢ ∑ [ yˆ ( k ) - y ] ⎥ ⎢⎣ k = 1 ⎥⎦ 2 ⎡ N ⎤ ⎢ ∑ [ yˆ ( k ) - y ( k ) ] ⎥ ⎢⎣ k = 1 ⎥⎦ (v − 1 ) (2-25) (N - v) Calculer le coefficient de corrélation multiple. Calculer les critères d’informations définis comme suit : a) L’erreur de prédiction finale (Final Prediction Error ) FPE ⎡1 FPE (v) = N Ln ⎢ ⎣N ∑ε N 2 k =1 ⎤ ⎡N + v⎤ (k) ⎥ + N Ln ⎢ ⎣ N - v ⎥⎦ ⎦ (2-26) b) Le critère d’information d’Akaike AIC. c) Le critère d’information de Bayes BIC. d) La loi du critère logarithmique itératif LILC (en prenant c=2 dans l’équation (2-18)) ⎡1 LILC(v) = N Ln ⎢ ⎣N ∑ ε (k) ⎥⎦ N 2 k =1 ⎤ + 2 v Ln [Ln(N) ] (2-27) 33 Chapitre 2 Etape 5 Cette étape contient deux tests [Boujmil, 1991] : Test 1 : 2 est supérieure ou égale à une valeur fixée a Si la valeur du coefficient de corrélation Rmult priori, on arrête la procédure. Test 2 : Tester si le F_Test est significatif en comparant sa valeur à une valeur fixée a priori (Xα). Comparer 2 Rmult , FPE(v), AIC(v), BIC(v), LILC(v) à leurs valeurs précédentes (valeurs à v–1) : 2 2 (v ) < Rmult (v − 1) la dernière variable sélectionnée est rejetée. Si Rmult ou si FPE(v) (ou respectivement AIC(v), BIC(v), LILC(v)) est supérieur à FPE(v-1) (respectivement AIC(v-1), BIC(v-1), LILC(v-1)), la dernière variable sélectionnée est rejetée. La procédure passe à l’étape 8 si la valeur du F_Test est non significative ou si la dernière variable sélectionnée est rejetée. Remarque : Dans la version originale de Kortmann [Kortmann, 1988], la procédure s’arrête si le F_test est non significatif ou si les critères d’informations sont supérieurs à ceux du modèle précédent. Dans [Boujmil, 1991], l’auteur a proposé d’effectuer deux tests assurant plus de chance à l’algorithme pour trouver un meilleur modèle. Etape 6 : Calculer pour chaque terme du modèle actuel (sauf la valeur moyenne) le F_Test partiel : ⎤ ⎡N ⎤ ⎡N ˆ [ y (k) y ]² − − ⎢∑[yˆ(k) − y]²⎥ ⎥ ⎢∑ ⎦v−1 variables , ( i = 1, ... , v-1) ⎣k=1 ⎦v variables ⎣k=1 PFi = N ⎛⎡ ⎞ ⎤ ⎜⎢∑[yˆ(k) − y(k)]2 ⎥ ⎟ (N− v) ⎜ k=1 ⎟ ⎦v variables⎠ ⎝⎣ v> 2 (2-28) et les critères d’informations partielles : AICi, BICi, FPEi et LILCi , i = 1, …,v-1, v > 2, en utilisant les équations (2-16), (2-17), (2-26) et (2-27). 34 Chapitre 2 Dans l’équation (2-28), la première somme du numérateur ainsi que la somme des résidus du dénominateur sont calculée à l’étape 4. La deuxième somme du numérateur de l’équation (2-28) et les critères d’informations partielles se calculent comme suit : Rejeter la première variable du modèle, estimer les paramètres restant du modèle et calculer : ∑ N k =1 [ yˆ (k) - y ] 2 (2-29) v- 1 variables et la somme des résidus : ∑ ε(k) N k =1 2 (2-30) v-1 variables avec les valeurs estimées, calculer le PFi et les critères d’informations partielles. Introduire le premier paramètre et rejeter le second. Pour ce nouveau modèle estimer les paramètres et calculer le PFi et les critères d’information partiels. Refaire cette opération jusqu’au dernier paramètre du modèle. Ces valeurs caractérisent la mesure de la contribution du terme considéré dans le modèle encore inconnu. Exemple On considère l’exemple de l’étape 3. Soit y ²(k-1 ) est le premier terme sélectionné. Si u(k-2 ) est le dernier terme sélectionné, le modèle actuel sera alors : yˆ (k) = θ12 y ²(k-1 ) + θ 2 u(k- 2 ) + ε(k). (2-31) Pour calculer le PF1 et les critères d’informations partiels, on procède comme suit : on rejette le terme y ²(k-1 ) d’où le nouveau modèle : yˆ (k) = θ 2' u(k- 2 ) + ε 1 (k) - (2-32) on estime θ 2' et on calcule yˆ 1 ( k ) = θ 2' u(k - 2) et ε 1 (k). puis on calcule le F_Test partiel PF1 donné par l’équation (2-28) et les différents critères d’informations partiels AIC1, BIC1, FPE1 et LILC1, donnés par les équations (2-16), (2-17), (2-26) et (2-27), 35 Chapitre 2 - on introduit le terme y²(k-1) et on rejette u(k-2 ) , ainsi le nouveau modèle sera : ' yˆ (k) = θ12 y ²(k-1 ) + ε 2(k). - (2-33) ' y ²(k-1 ) et ε 2 (k) puis on calcule On estime le paramètre θ12' et on calcule yˆ 2 (k) = θ12 le F_Test partiel correspondant PF2 et les différents critères d’informations partiels correspondants AIC2, BIC2, FPE2 et LILC2. Etape 7 : Comparer la plus petite valeur de AICi (respectivement BICi, FPEi et LILCi) à la valeur de AICv (respectivement, BICv, FPEv et LILCv) calculée à l’étape 4. si ou min ( AICi ) < AICv min (FPEi ) < FPEv i i ou min (BICi ) < BICv ou min (LILCi ) < LILCv i i alors la variable correspondante est rejetée. - De plus, tester si la valeur du F_Test partiel est significative ou non en comparant sa valeur à la valeur Xα fixée a priori (valeur de la distribution de Fisher donnée à l’étape 1). Si PFi < Xα, la variable correspondante est rejetée. La variable est retenue dans le cas contraire. Si toutes les variables sont retenues, la procédure passe à l’étape 8. Dans le cas où une variable est rejetée, on décrémente le nombre des paramètres v (v=v-1) et on passe à l’étape 4 pour refaire l’estimation paramétrique et l’évaluation du modèle. Si la dernière variable est rejetée, on choisit alors le deuxième terme de plus grand coefficient de corrélation partiel et on passe à l’étape 4. Par ailleurs, si la dernière variable sélectionnée est rejetée deux fois, la procédure s’arrête. Etape 8 : Calculer les coefficients de corrélation partiels normalisés de tous les termes restants : 36 Chapitre 2 N ∑ εi′(k) ε(k) N ρi = k =1 ∑εi′(k) N k =1 ∑ ε(k) N k =1 N N N ⎤ ⎡ N 2 2 ⎤ ⎡ 2 ′ ′ N ( ε (k)) [ ε (k) ] N ε (k) [ ε(k) ]2 ⎥ ∑ ∑ i ⎥⎢ ∑ ⎢ ∑ i k =1 k =1 ⎦ ⎦ ⎣ k =1 ⎣ k =1 (2-34) où ε(k), k = 1, ... , N est calculé déjà à l’étape 4 et ε ′(k ) est obtenu comme le montre l’exemple suivant : Exemple On reprend le modèle présenté dans l’exemple de l’étape 6. Le modèle actuel est donné par : yˆ (k) = θ12 y²(k-1 ) + θ 2 u(k- 2 ) + ε(k) . (2-35) Les termes restant dont on veut calculer leurs coefficients de corrélation partiels sont : v1 (k ) = u(k - 1) v 3 (k ) = y(k - 1) v 8 (k ) = u(k - 1) ⋅ y(k - 1) v 9 (k ) = u(k - 1) ⋅ y(k - 2) v 4 (k ) = y(k - 2) v10 (k ) = u(k - 2) y(k - 1) v 6 (k ) = u(k - 1) ⋅ u(k - 2) v13 (k ) = y(k - 1) y(k - 2) v5 (k ) = u 2 (k - 1) v 7 (k ) = u 2 (k - 2) v11 (k ) = u(k - 2) y(k - 2) (2-36) v14 (k ) = y 2 (k - 2) Soit v1 (k ) le premier terme que l’on veut calculer son coefficient de corrélation. On considère alors le nouveau modèle ayant pour sortie v1 (k ) et pour entrée le modèle actuel : v1 (k ) = u (k − 1) = θ '12 y²(k - 1) + θ ' 2 u(k - 2) + ε '1 (k) (2-37) On estime les paramètres et on calcule ε '1 (k) par simulation d’où : ε 1′ (k) = u(k - 1) - [θ12′ ⋅ y²(k - 1) + θ 2′ ⋅ u(k - 2)] On répète le même travail pour les autres termes : v i (k), (2-38) i = 3, L ,14 et i ≠ 12. 37 Chapitre 2 Commentaires : La stratégie de cet algorithme consiste à construire le modèle en commençant par le terme le mieux corrélé avec la sortie. Puis, on ajoute par la suite des nouveaux termes jusqu'à l’obtention du modèle final. 2 , Cet algorithme utilise plusieurs tests statistiques, des critères de performances (OVF, Rmult FPE, AIC, BIC, LILC, …), des critères d’information partiels (FPEi, AICi, BICi, LILCi) ainsi que des coefficients de corrélation et des coefficients de corrélation partiels normalisés. Bien que l’adoption de tous ces critères et de ces tests statistiques offre une robustesse à cet algorithme, la difficulté lors de leurs programmations reste contraignante. L’implémentation de l’algorithme nécessite beaucoup de précaution pour ne pas faire des confusions ou des erreurs, surtout dans les étapes 6 et 8 lors du calcul des critères d’information partiels et des coefficients de corrélation partiels. De même cette difficulté se présente dans l’étape 7 pour décider quels sont les termes à rejeter de la structure actuelle. C’est pour ces raisons qu’on a choisi des exemples simples, expliquant les différentes étapes de cet algorithme. La complexité de cette méthode et son limite aux cas des systèmes monovariables nous ont incités à concevoir d’autres méthodes qui utilisent d’autres approches, moins complexes, et qui peuvent êtres extensibles pour traiter les systèmes non linéaires multivariables. Nous proposons dans ce qui suit deux nouvelles approches pour l’identification des modèles de type NARMA. 2.3.2. Méthode basée sur les Algorithmes Génétiques Binaires 2.3.2.1. Présentation de l’algorithme génétique L’algorithme génétique est appliqué à de nombreux problèmes d’optimisation [Goldberg, 1989]. Il a été mis au point dans les années 70 par l’équipe de John Holland à l’Université du Michigan. Après être demeuré mal connu pendant des années, l’intérêt pour cet algorithme s’est accru récemment. L’avantage de cette méthode d’optimisation est qu’elle trouve une solution qui s’approche de la solution optimale, et ce même pour des problèmes fortement non linéaires [Goldberg, 1989]. Cet algorithme donne de bonnes performances même lorsque le nombre de variables est très élevé. De nombreuses variantes quant aux opérateurs et à l’encodage utilisés existent dans la littérature [Goldberg, 1989], [Man, 1996]. L’algorithme 38 Chapitre 2 génétique peut fonctionner également, soit avec une représentation réelle soit avec une représentation binaire des variables [Goldberg, 1989]. Les algorithmes génétiques sont inspirés du phénomène de la sélection naturelle qui permet à la nature de créer avec succès des organismes adaptés à leur environnement [Goldberg, 1989]. Les organismes qui sont le moins adaptés disparaissent peu à peu tandis que ceux qui survivent et qui sont mieux adaptés peuvent se reproduire. Les descendants sont semblables à leurs parents, mais pas identiques. Ils retiennent une partie de leurs parents. Lorsque l’environnement change, les espèces s’adaptent graduellement. Les mécanismes qui gèrent la recherche d’un optimum dans un algorithme génétique sont déduits de ce phénomène. Des opérateurs génétiques servent à faire évoluer artificiellement un ensemble de vecteurs. Par analogie, l’ensemble des vecteurs forme une population et chaque vecteur constitue un individu. Dans les algorithmes génétiques simples, trois opérateurs principaux sont utilisés : la reproduction, la mutation et le croisement. Les opérateurs modifient de façon successive la population, ce qui permet d’explorer l’espace de recherche afin de trouver un optimum. Une fonction d’évaluation sert à évaluer l’adaptation des individus. Cette fonction évalue l’adéquation de l’individu à la solution recherchée. Plus la valeur de la fonction évaluée est élevée, plus cet individu est adapté. L’application des trois opérateurs sur une population d’individus représente une itération ou une génération. On présente ci-dessous les principaux opérateurs de traitement de l’algorithme génétique. • Reproduction : La reproduction est une opération selon laquelle les individus d’une population sont copiés dans la population de la génération suivante selon la fonction d’évaluation. Plus la valeur obtenue par la fonction d’évaluation est grande, plus la probabilité d’être copiée dans la prochaine population est élevée. La sélection étant aléatoire, les meilleurs ont tendance à être copiés plusieurs fois dans la prochaine population proportionnellement à l’importance de la valeur de la fonction d’évaluation des parents dans la population. Ceux qui sont moins adaptés sont probablement abandonnés et non présents dans la population suivante. Cette phase a pour but d’accroître l’importance des bonnes solutions. Elle tend à augmenter l’uniformité de la population. • Croisement : C’est un opérateur qui permet de mélanger les caractéristiques de deux individus dans la recherche de nouvelles solutions. Le principe consiste à utiliser l’information contenue dans deux solutions afin d’en produire deux nouvelles. Dans le cas du 39 Chapitre 2 croisement uniforme les couples sont créés en sélectionnant deux individus de façon aléatoire. Ensuite, les caractéristiques des individus sont mélangées de façon aléatoire. • Mutation : La mutation consiste à modifier aléatoirement la valeur d’un individu. La probabilité que cette modification ait lieu est en réalité très faible. Cet opérateur permet de changer la direction de l’exploration et ainsi de ne pas garder un optimum local. • Fonction d’évaluation : Cette fonction évalue l’adéquation de l’individu à la solution recherchée. Dans un problème de maximisation, plus la valeur de la fonction est élevée, plus cet individu est adapté. L’initialisation des individus de la population initiale se fait de façon aléatoire. L’algorithme s’arrête lorsqu’il y a convergence ou après un nombre d’itérations (ou générations) prédéfini. Les paramètres à ajuster pour l’algorithme génétique sont : - Le nombre maximal de générations : qui spécifie le nombre d’itérations. - Le nombre d’individus formant une population : il représente en quelque sorte la mémoire de l’algorithme, plus ce nombre est grand, plus l’algorithme a tendance à garder de bonnes parties de solutions. - La probabilité de croisement : elle spécifie le taux avec lequel les deux parents (individus) sont mélangés lors du croisement. - La probabilité de mutation : cette probabilité contrôle la variation aléatoire des individus et dicte alors l’exploration de l’espace de recherche. - Le nombre d’évaluations requises (Ne) pour une recherche par algorithme génétique correspond au produit du nombre de générations effectué (Ng) avec le nombre d’individus (Ni). Ne = Ng Ni (2-39) Dans le paragraphe suivant, on présente l’identification des paramètres du modèle NARMA en utilisant un algorithme génétique binaire. 2.3.2.2. Méthode basée sur les algorithmes génétiques binaires Dans un algorithme génétique binaire simple, chaque individu est représenté par un vecteur binaire qui désigne la variable de la fonction à optimiser. Ce vecteur binaire sera codé 40 Chapitre 2 pour donner une valeur réelle dans l’intervalle espace de solutions [Goldberg, 1989]. La souplesse que présente ce type de codage permet de résoudre plusieurs types de problème [Fnaiech, 2000]. Cet algorithme a été utilisé depuis le début de la dernière décennie, pour l’identification des systèmes de type entrée-sortie [Li, 1993], [Li, 1994], [Luh, 1998], [yang, 2000] et [Madár, 2005]. Les méthodes présentées exploitent différents types de codages. Dans [Li, 1993], les auteurs ont proposé un codage binaire où chaque bit représente un terme du modèle à déterminer. Avec cette méthode, le nombre de termes dans un individu peut dépasser le nombre maximal de termes possibles, ce qui engendre l’annulation de la solution trouvée et la création d’une nouvelle génération. On note aussi que la condition d’arrêt dépend de l’écart entre la valeur du meilleur fitness actuel et celle de l’ancien meilleur fitness. Si cet écart est inférieur à une certaine précision, la procédure sera arrêtée et le modèle sera donné après avoir éliminé les termes dont les coefficients sont proches de zéro. L’application de cette condition d’arrêt et la possibilité d’avoir des solutions non adaptées au cours du fonctionnement de l’algorithme, engendrent un temps très important lors de l’identification. Le codage proposé dans [Luh, 1998], consiste à diminuer le nombre de termes possibles dans le modèle en divisant chaque individu en plusieurs sous vecteurs binaires et en conservant les termes significatifs inchangés. Cette approche est basée sur les résultats de [Leontaritis, 1985b] qui montrent que 10 termes sont généralement suffisants pour la modélisation d’un système, avec une légère détérioration du modèle. Puis, dans [yang, 2000] les auteurs ont utilisé une nouvelle technique de codage dite hiérarchique. Ce codage permet de classer les chromosomes sous formes de couches. Dans chaque couche on trouve les termes qui représentent le même degré de non linéarité, ainsi que leur nombre. La couche d’ordre 0 contient la valeur constante (valeur moyenne), la couche d’ordre 1 représente les termes simples (les régressions linéaires sur l’entrée et sur la sortie). Ensuite, on trouve les différentes combinaisons (les doublés puis les triplés, etc…) jusqu’à l’atteinte de l’ordre du polynôme. La méthode présentée par [Madár, 2005], consiste à mettre les individus de chaque population de l’algorithme génétique sous forme arborescente. Les branches de chaque individu présentent les termes possibles du modèle à définir. Ainsi, l’algorithme OLS (Orthogonal Least Squares) consiste à donner la contribution de chaque branche et donc permet d’éliminer les termes qui ne sont pas significatifs dans le modèle. On note que dans les méthodes [Li, 1993] et [yang, 2000], la valeur moyenne est considérée dans le codage binaire, d’où la possibilité d’absence de ce terme de la structure du 41 Chapitre 2 modèle. Dès lors, si on veut garder ce terme constant dans l’équation du modèle, le fonctionnement de l’algorithme génétique sera forcé à chaque itération. Cependant, dans le travail de [Madár, 2005], ce terme constant ne figure pas dans la structure arborescente. Notre méthode, contrairement aux autres méthodes, permet de considérer le terme constant dans l’évaluation du fitness, au cours de l’estimation des paramètres, sans l’avoir intégré dans le codage de l’algorithme génétique. En plus, cette approche permet d’avoir plusieurs structures différentes présentant des performances proches. Cet avantage permet le choix d’une structure optimale en faisant un compromis entre la complexité du modèle et sa précision. Dans ce travail, le vecteur binaire est l’image du vecteur d’observation Φ (k ) . Une ligne de la population de l’algorithme génétique, notée Ipop, ainsi que le vecteur d’observation sont donnés par : Ipop = [1 0 Φ ( k ) = [1 v1 1 v2 ... 1] ... v r ] (2-40) (2-41) L’analyse de Ipop permet de sélectionner les termes intervenants dans l’expression de la sortie du modèle. Seuls les termes dont la valeur du bit correspondant est à ‘’1‘’ logique sont sélectionnés, les autres ne seront pas considérés dans l’expression de la sortie du modèle. Chaque individu comporte autant de valeurs booléennes que de termes possibles du modèle NARMA. L’exemple suivant explique le codage du vecteur binaire. Dans le cas où (m=n=q=2), on a 14 termes possibles. On considère alors une représentation binaire de l’algorithme génétique où chaque individu comporte 14 bits, chaque bit représente un terme du vecteur d’observation. Un individu peut être présenté par le vecteur suivant : [1 1 0 1 1 0 0 0 0 0 0 1 0 0] (2-42) Dans ce vecteur, les bits 1, 2, 4, 5 et 12 sont à ‘’1’’, les autres bits sont à ‘’0’’. Ceci correspond au modèle défini par la structure suivante : yˆ (k ) = y + b1 u(k −1) + b2 u(k − 2) + a2 y(k − 2) + b11 u 2 (k −1) + a11 y 2 (k −1) (2-43) Ensuite, la procédure de MCR est utilisée pour estimer les paramètres du modèle considéré (y compris la valeur moyenne). Pour chaque individu, on calcule la fonction fitness afin de l’utiliser dans la phase de sélection. La fonction d’évaluation (fitness) dépend de 42 Chapitre 2 l’erreur entre la sortie du modèle obtenu et la sortie du système. Cette fonction est définie par la relation suivante : fitness = 2 avec J 1 = ∑ ( y ( k ) − yˆ (k ) ) N 1 1+ J 1 k =1 (2-44) Avec cet algorithme, on utilise les opérateurs génétiques traditionnels (sélection, croisement et mutation) comme le montre la figure 2-2. Les caractéristiques de l’algorithme génétique considéré sont les suivantes : - La population initiale : la population initiale est choisie d’une façon aléatoire. La taille de chaque individu est égale à la taille du vecteur d’observation (sans considérer la valeur moyenne). - La sélection : la sélection consiste à choisir les individus qui sont caractérisés par les meilleurs fitness. La probabilité de sélection des individus est proportionnelle à leurs fitness et elle est donnée par : Pi = ∑ fitness fitness i (2-45) max _pop j =1 j avec fitnessi est le fitness de l’individu i et max_pop est le nombre total d’individus. Dans ce travail, on a utilisé la sélection basée sur la roue de loterie [Goldberg, 1989]. - Le croisement : l’objectif du croisement est l’échange d’information entre deux parents pour former deux nouveaux individus. Dans le présent travail, nous avons opté pour le croisement à un cite avec une probabilité Pc. - La mutation : la mutation est un processus aléatoire qui consiste à modifier l’état d’un bit de 0 à 1 ou vice versa. La probabilité de mutation Pm est généralement choisie de valeur très faible. - Condition d’arrêt : l’algorithme génétique évolue d’une génération à une autre jusqu’à une condition d’arrêt est atteinte. Dans la littérature, il existe plusieurs critères qui peuvent être utilisés comme condition d’arrêt de l’algorithme génétique. On cite par exemple l’utilisation d’un nombre maximal de génération, le test d’uniformité des individus de la population après un certain nombre de génération 43 Chapitre 2 ou bien la non amélioration du fitness. Dans notre travail, on a considéré le nombre maximal de génération de l’algorithme génétique comme condition d’arrêt. La solution obtenue, avec cette approche, détermine le modèle identifié du système. Le meilleur modèle est celui qui donne l’erreur de modélisation la plus faible. Cependant, avec cette approche, on peut avoir plusieurs solutions donc plusieurs modèles puisque l’algorithme génétique est une méthode non déterministe. Initialisation de l'algorithme Création de la population initiale Evaluation des chromosomes de la population - Sélection des termes : Φ (k ) - Estimation des paramètres par MCR : Θ̂ ˆ T Φ (k ) yˆ (k ) = Θ J1 = ∑ ( y(k ) − yˆ (k ))2 - Calcul de la sortie du modèle : N - Calcul du critère: k =1 - Calcul du fitness : 1 1 + J1 Sélection Croisement et mutation Insertion des nouveaux chromosomes dans la population suivante. Condition d'arrêt Arrêt et renvoi du meilleur chromosome produit Figure 2-2 : Organigramme de l’identification du modèle NARMA avec l’algorithme génétique binaire 44 Chapitre 2 2.3.3. Méthode basée sur les Réseaux de Neurones Artificiels 2.3.3.1. Réseau de neurones artificiels à rétro-propagation de l’erreur Les réseaux de neurones sont des fonctions mathématiques complexes et non linéaires. Ces réseaux sont composés de couches, qui sont elles mêmes composées d’éléments appelés neurones. Chaque neurone a pour entrée les sorties de chacun des neurones de la couche précédente. Ces entrées sont pondérées par les poids, pour ensuite être combinées et traitées par une fonction d’activation. L’objectif de l’apprentissage est d’obtenir ces poids. L’ajustement des poids se fait de façon à réduire l’erreur quadratique entre les valeurs obtenues en sortie et les valeurs désirées. Le gradient de l’erreur est utilisé à cette fin. Le processus se poursuit jusqu’à la convergence, c’est-à-dire tant que l’erreur totale n’atteint pas une valeur inférieure à un seuil. x1 ∑ f x2 x3 xn Couche d’entrée ŷ ∑ f ∑ f Couche cachée ∑ Couche de sortie Figure 2-3 : Schéma descriptif d’un Réseau de neurone non bouclé à une seule couche cachée f étant la fonction d’activation des neurones. Il existe plusieurs fonctions d’activation utilisées en pratique tel que la fonction tangente hyperbolique (tanh) et les fonctions (signe, sigmoïde, polynomiale, …). La fonction polynomiale est utilisée dans notre cas. 45 Chapitre 2 2.3.3.2. Méthode d’identification utilisant les réseaux de neurones artificiels et l’algorithme génétique Les réseaux de neurones à trois couches illustrés à la figure 2-3, sont utilisés dans ce travail. Soit le vecteur d’entrée : X ( k ) = [ y(k-1 ) ⋅ ⋅ ⋅ y(k-n) u ( k − 1) ⋅ ⋅ ⋅ u ( k − m ) ] (2-46) m et n désignes respectivement l’ordre de la régression sur l’entrée et l’ordre de la régression sur la sortie. Soient W1 et W2 respectivement la matrice des pondérations entre la couche d’entrée et la couche cachée et le vecteur des pondérations entre la couche cachée et la couche de sortie. Les paramètres du modèle NARMA peuvent être estimés en utilisant un réseau de neurones à fonction d’activation polynomiale [Ki, 1997]. La structure du réseau de neurones est de type «feed-forward» avec une seule couche cachée. La topologie de ce réseau est représentée dans la figure 2-4. Le réseau neuronal reçoit en entrée les éléments du vecteur X (k ) et délivre à sa sortie une estimation de la sortie actuelle. ∑/f y(k-1) ŷ(k) y(k-n) ∑ u(k-1) ∑/f u(k-m) W1 Couche d'entrée W2 Couche cachée Couche de sortie Figure 2-4 : Réseau de neurones à une seule couche cachée 46 Chapitre 2 Les sorties de la première couche sont données par : S (k ) = f (W 1 X T (k )) = [s1 ( k ) ... s M ( k )]T , (2-47) f étant une fonction d’activation polynomiale. La relation suivante présente une fonction d’activation polynomiale d’ordre q : f ( xi ) = f 0 + f1 xi + f 2 xi2 + ... + f q xi q (2-48) En tenant compte des poids des connexions entre la couche d’entrée et la couche cachée, la quantité xi sera exprimée par : xi = ∑ m j =1 w1ji u (k − j ) + ∑ w1ji y(k − j + m), i ∈ [1, M ] m+n j = m +1 (2-49) avec M le nombre de neurones de la couche cachée. Le choix de M dépend de la complexité du système. En effet, un nombre faible de neurones cachés peut conduire à un modèle non précis alors que le choix d’un nombre élevé de neurones cachés peut rendre le modèle très compliqué et le temps de calcul sera plus important. En général, le nombre de neurones cachés est choisi de façon à réaliser un compromis entre la précision et la complexité du modèle estimé. L’équation de la sortie du réseau sera alors: yˆ (k ) = ∑ wi2 f ( xi ) = W 2 S (k ) M i =1 (2-50) Les dimensions de W1 et W2 sont fonction de m, n et du nombre de neurones dans la couche cachée. Pour un nombre de neurones dans la couche cachée égal à (M) les dimensions de W1 et W2 seront comme suit : dim (W1)= (M, m+n) et dim (W2) = (M, 1). (2-51) En remplaçant la quantité f ( xi ) de l’équation (2-48) par sa valeur dans l’équation (2-50), on obtient : 47 Chapitre 2 ( + w (f ) x ) yˆ ( k ) = w12 f 0 + f1 x1 + ... + f q x1q 2 2 + f1 x2 + ... + f q ( ) q 2 + ... + wM2 f 0 + f1 x M + ... + f q x Mq . 0 (2-52) En remplaçant (xi) par son expression donnée par l’équation (2-49), l’expression de la sortie du réseau pour un degré de non linéarité q = 2 sera comme suit : yˆ ( k ) = f 0 + f1 + f2 + f2 + f2 ∑ M i =1 wi2 + f1 ∑ (w w m+n j = m +1 m j =1 l =1 2 1 1 1 j1wl1 ∑ ∑ (w w m+n m+n j = m +1 l = m +1 m+n j =1 l = m +1 )y(k − j + m) ) )u (k − j ) 2 1 + ... + wM w jM w1lM u ( k − j )u ( k − l ) 2 1 1 2 1 1 j1 wl1 + ... + wM w jM ∑ ∑ (w w m j =1 2 1 2 1 1 j1 + ... + wM w jM 2 1 2 1 1 j1 + ... + wM w jM ∑∑ (w w m ∑ (w w m 2 1 1 1 j1wl1 ) (2-53) w1lM y ( k − j + m ) y ( k − l + m ) ) 2 1 + ... + wM w jM w1lM u ( k − j ) y ( k − l + m ) Les pondérations (W1 et W2) des connexions du réseau de neurones sont déterminées avec l’algorithme de rétro propagation du gradient. La déduction des paramètres du modèle NARMA se fait à partir des pondérations (W1 et W2) et des paramètres du polynôme considéré (f0, f1,…, fq). En effet, par identification de l’équation (2-53) avec l’équation (1-7), on trouve les coefficients des termes du modèle NARMA [Ki, 1997]. Exemple : pour q = 2 et (m, n) quelconques, on aura le polynôme d’ordre 2 suivant : f ( x i ) = f 0 + f 1 x i + f 2 x i2 (2-54) Les paramètres du modèle NARMA seront alors obtenues comme suit : 48 Chapitre 2 y = ∑ ai = bi = a ij = M f 0 w 2j j =1 ∑w 2 j f 1 w ij1 i = 1 ,...,m 2 j f 1 w ij1 i = m + 1 ,...,m + n M j =1 ∑w M j =1 ∑w i = j, i = 1 ,...,m i ≠ j, i = 1 ,...,m ; j i = j, i = m + 1 ,...,m + n M j =1 2 j f 2 w ij1 w 1jl a ij = 2 ∑ w 2j f 2 w ij1 w 1jl M b ij = j =1 ∑w M j =1 2 j f 2 w ij1 w 1jl b ij = 2 ∑ w 2j f 2 w ij1 w 1jl M j =1 c ij = 2 ∑ w 2j f 2 w ij1 w 1jl i ≠ j, M j =1 = i + 1 ,...,m i = m + 1 ,...,m + n ; j = i + 1 ,...,n i = 1 ,...,m ; j = m + 1 ,...,m + n (2-55) La qualité du modèle obtenu, avec cette méthode, dépend du choix adéquat des coefficients du polynôme (f0, f1,…, fq). Dans [Ki, 1997], on ne retrouve aucune indication concernant le choix de ces coefficients. Dans le but de réduire le nombre d’essais, on propose l’incorporation d’une procédure basée sur l’algorithme génétique réel avec l’algorithme de rétro propagation du gradient (figure 2-5). u(k) y(k) Procédé - Modèle RNA + yˆ ( k ) W1, W2 f0, f1,…,fq Algorithme d’adaptation Algorithme Génétique Figure 2-5 : Structure du superviseur basé sur l’algorithme génétique Chaque individu de la population de l’algorithme génétique caractérise les coefficients du polynôme de la fonction d’activation des neurones. La population initiale de l’algorithme 49 Chapitre 2 génétique est choisie d’une manière aléatoire à partir d’un intervalle fixé. Chaque génération de l’algorithme génétique consiste à prendre les coefficients donnés par l’individu de la population et à réaliser l’apprentissage du réseau neuronal avec l’algorithme de rétro propagation. Pour chaque individu, on calcule une fonction d’évaluation (i.e. : la somme des erreurs quadratiques entre la sortie du système et celle du réseau de neurones). Ensuite, les opérateurs génétiques (sélection, croisement et mutation) sont utilisés pour former les individus de la nouvelle population. Dans ce travail, on a considéré la sélection basée sur la roue de loterie, le croisement arithmétique et la mutation uniforme. Les nouveaux individus (e1 et e2) sont obtenus à l’aide des parents (p1 et p2) en utilisant les deux relations suivantes : e1 = α c p1 + (1 − α c ) p 2 e2 = (1 − α c ) p1 + α c p 2 , α c ∈ [0, 1] (2-56) (2-57) La mutation consiste à modifier aléatoirement la valeur de l’individu en respectant son intervalle de variation. Les étapes de l’algorithme génétique seront activées jusqu’à une condition d’arrêt est atteinte. Dans notre travail, on a considéré le nombre maximal de génération de l’algorithme génétique comme condition d’arrêt. La population initiale de l’algorithme génétique est choisie d’une manière aléatoire. De même, les valeurs des matrices de pondération initiales du réseau de neurones sont choisies aléatoirement dans l’intervalle [-1,1]. Le paragraphe suivant comporte les résultats de simulation pour les trois méthodes présentées précédemment. 2.4. Résultats de simulation Le choix des paramètres (m, n et q) a été fixé après le calcul du coefficient de corrélation total R²tot, qui est exprimé par l’équation (2-11). Le coefficient de corrélation total est calculé en utilisant la méthode de Kortmann. Les paramètres (m, n et q) qui donnent la plus grande valeur de R²tot sont alors utilisés pour caractériser le modèle NARMA. Trois critères d’informations sont utilisés (les coefficients de corrélations : total R²tot, multiple R²mult et ajustable R²ajust) pour décider si le modèle obtenu est acceptable ou non [Haber, 1990]. En 50 Chapitre 2 effet, si l’un de ces critères est proche de ‘’1’’ alors le modèle est acceptable [Haber, 1990]. Les paramètres de réglage utilisés dans les simulations sont résumés dans le tableau 2-1. Algorithme génétique binaire Algorithme génétique réel Réseau de neurones artificiels - nombre de génération 40 40 - - nombre d’individus 50 100 - - probabilité de croisement 0.4 0.8 - - probabilité de mutation 0.05 0.1 - - nombre de passages - - 3000 - coef. d’apprentissage - - 5.10-3 Tableau 2-1 : Paramètres des algorithmes Au cours de la phase d’identification, le signal d’entrée utilisé est un bruit blanc de moyenne nulle et de variance égale à un. On a considéré la notation suivante : Méthode K : méthode basée sur l’algorithme de Kortmann. Méthode AGB : méthode basée sur l’Algorithme Génétique Binaire. Méthode RN : méthode basée sur l’utilisation des Réseaux de Neurones Artificielles. On a considéré quatre systèmes non linéaires et monovariables. Système 1 : le premier système simulé est caractérisé par l’équation aux différences suivante : y ( k ) = 0 . 3528 u ( k − 1) + 0 . 9844 y ( k − 1) + 0 . 1892 u 2 ( k − 1) + 0 . 6235 u ( k − 1) y ( k − 1) − 0 . 283 y 2 ( k − 1). (2-58) Dans ce cas, le modèle NARMA estimé est caractérisé par 14 termes (n=m=q=2). Le tableau 2-2 comporte les termes sélectionnés ainsi que les paramètres estimés avec les trois méthodes. Les coefficients du polynôme de la fonction d’activation identifiés par l’algorithme génétique sont les suivants : f0 = -4.6758; f1=14.5596; f2=15. D’après les résultats du tableau 2-2, on constate que les modèles identifiés sont pratiquement identiques. En effet, les termes des modèles ainsi que les paramètres correspondants sont identifiés avec succès par les trois méthodes. Les valeurs des coefficients de corrélation (R²tot, R²mult et R²ajust) sont toutes proches de un. Ceci signifie que les modèles 51 Chapitre 2 identifiés avec les différentes méthodes caractérisent convenablement le comportement du système, par conséquent, ils sont tous acceptables. Pour valider les modèles trouvés, on a comparé les réponses des modèles avec la réponse du système considéré pour une séquence du signal d’entrée différente de la séquence utilisée en identification. La figure 2-6, présente la séquence d’entrée et les sorties des modèles identifiés. Cette figure montre que la sortie du système et les sorties des modèles trouvés sont presque confondues. Les modèles trouvés avec les différentes méthodes N° Termes possibles Modèle Simulé (termes et paramètres) 1 u(k-1) x 2 u(k-2) 3 y(k-1) 4 y(k-2) 5 (u(k-1))² 6 u(k-1)u(k-2) 0 7 (u(k-2))² 0 8 u(k-1)y(k-1) 9 u(k-1)y(k-2) 0.0001 10 u(k-2)y(k-1) -0.0001 11 u(k-2)y(k-2) 0 12 (y(k-1))² 13 y(k-1)y(k-2) 0 14 (y(k-2))² 0 15 y 0.3528 K AGB RN 0.3528 0.3527 0.3528 0 x 0.9844 0.9843 0.9842 0.9845 0 x x x 0.1892 0.6235 -0.283 0 0.1889 0.6235 -0.2831 0.1879 0.6236 -0.2830 0.1892 0.6234 -0.283 0.0005 0.0021 0 R²tot 0.9997 0.9918 0.9996 R²mult 1.0075 1.0418 0.9973 R²ajust 1.0076 1.0422 0.9972 Tableau 2-2 : Structures de l’équation simulée et des modèles identifiés (système 1) 52 Chapitre 2 ent r ée u( k) 1 0.5 0 -0.5 -1 0 50 100 150 200 250 300 350 400 Evolut ions des sor t ies ( syst èm e et m odèles ident ifiés) 2 y yK yAGB 1 yRN 0 -1 0 50 100 150 200 250 300 350 400 It ér at ions Figure 2-6 : Les évolutions de la séquence d’entrée, de la sortie du système et des sorties des modèles considérés (système 1) Afin d’étudier l’influence de bruit sur la qualité de l’estimation, on a ajouté à la sortie du système un bruit blanc de moyenne nulle et de variance finie. Le tableau 2-3 comporte les paramètres estimés par les trois méthodes dans le cas où la variance de bruit est égale à 7.1 10-3. Le tableau 2-4 présente les valeurs du coefficient de corrélation multiple R²mult pour différentes valeurs de la variance de bruit. L’analyse des résultats donnés sur les tableaux 2-3 et 2-4 permet de conclure que la présence de bruit de mesure entraîne une légère dégradation de la qualité des modèles trouvés. En effet, les paramètres estimés sont biaisés et les modèles sont moins précis puisque les valeurs de R²mult diminuent au fur et à mesure que la variance de bruit augmente. 53 Chapitre 2 Les modèles trouvés avec les différentes méthodes Numéro du paramètre K AGB RN 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 0.3579 0.3552 0.9286 0.9961 0.1806 0.1876 0.5951 0.6374 -0.0208 -0.0082 -0.2928 -0.2807 -0.0176 0.012 0.0021 0.3597 0.0300 0.9357 0.0165 0.1932 0.0108 -0.0002 0.6130 0.0125 0.0414 0.0012 -0.3074 0.0449 -0.0011 -0.0366 Tableau 2-3 : Paramètres du modèle NARMA pour une variance du bruit égale à 7.1 10-3 (système 1) R²mult Variance du bruit K AGB RN 3.1 10-3 0.9837 0.9895 0.9854 7.1 10-3 0.9709 0.9814 0.9737 13 10-3 0.9516 0.9626 0.9545 Tableau 2-4 : Les valeurs de R²mult pour différentes valeurs de la variance du bruit (système 1) Système 2 : le fonctionnement du deuxième système est décrit par la relation suivante [Narendra, 1997] : y (k ) = y ( k − 1) 1 + y ( k − 1) 2 + u 2 ( k − 1) (2-59) Le choix (n=m=q=2) permet d’avoir la plus grande valeur du coefficient de corrélation total R²tot. Les différents termes possibles du modèle NARMA sont présentés dans le tableau 2-5. Les coefficients du polynôme de la fonction d’activation identifiés par l’algorithme génétique réel sont les suivants : f0=0.4388; f1=5.8 et f2=7.41. D’après les résultats obtenus on constate que : - le modèle identifié par la méthode de Kortmann, comporte seulement cinq termes, - la méthode utilisant l’algorithme génétique binaire permet d’avoir trois représentations possibles, donc trois modèles possibles malgré leurs ressemblances. 54 Chapitre 2 La première et la deuxième représentation comportent 9 termes et la troisième comporte 10 termes. Si on écarte les termes dont la valeur absolue du paramètre correspondant est inférieure à (0.01) on aura alors un modèle qui comporte seulement sept termes, - la troisième méthode suppose l’existence de tous les termes possibles du modèle NARMA. La sélection des termes du modèle choisi est basée sur les valeurs des paramètres correspondants. Ces paramètres, contrairement aux autres méthodes ne sont pas estimés par la méthode de MCR. Dans cet exemple, si on écarte les termes dont les paramètres sont inférieurs à (0.001) en valeur absolue, on aura alors un modèle qui comporte seulement sept termes, - les valeurs des coefficients de corrélation (R²tot, R²mult et R²ajust) sont toutes proches de un. Ceci signifie que les modèles identifiés avec les différentes méthodes caractérisent convenablement le comportement du système, par conséquent, ils sont tous acceptables. Les modèles identifiés avec les différentes méthodes Termes K AGB RN Termes et paramètres estimés u(k-1) -0.0121 -0.0118 -0.0119 0.0034 u(k-2) -0.0165 -0.0148 -0.0172 -0.0004 0.3559 0.3490 0.3495 -0.0146 0 -0.0032 -0.0034 -0.0045 0.9988 0.9991 0.9994 0.9993 1.0034 -0.3923 -0.2316 -0.2246 -0.2254 0 0.0136 0.0133 0.0133 -0.003 y(k-1) 0.5223 y(k-2) (u(k-1))² u(k-1)u(k-2) (u(k-2))² -0.0027 u(k-1)y(k-1) u(k-1)y(k-2) -0.0023 u(k-2)y(k-1) 0.0082 u(k-2)y(k-2) (y(k-1))² -0.0392 0.0080 0.0081 0 0.0016 0 0.0024 0 -0.0388 -0.0389 -0.0388 0 y(k-1)y(k-2) 0 (y(k-2))² 0 y 0.1897 0.2632 0.2690 0.2690 0.4418 R²tot 0.9999 0.9965 0.9965 0.9965 0.9972 R²mult 0.9864 0.9844 0.9844 0.9843 0.9922 R²ajust 0.9863 0.9841 0.9840 0.9840 0.9917 Tableau 2-5 : Structures des modèles identifiés (système 2) 55 Chapitre 2 Pour valider les modèles trouvés, on a comparé la réponse du système avec les réponses des modèles. Dans ce cas, la séquence d’entrée est un bruit blanc de moyenne nulle et de variance finie. La figure 2-7 présente la séquence d’entrée, la sortie du système et les sorties des modèles trouvés. La figure 2-8 présente les différentes sorties entre les itérations 19 et 41. On remarque que les sorties des modèles identifiés sont pratiquement confondues avec la sortie du système. Figure 2-7 : Les évolutions de la séquence d’entrée, de la sortie du système et des sorties des modèles considérés (système 2) Figure 2-8 : Les évolutions des sorties (système et modèles). Pour vérifier la robustesse des différents modèles identifiés, on a ajouté à la sortie du système un bruit de variance finie. Le rapport entre la variance du bruit et l’amplitude du 56 Chapitre 2 signal d’entrée est noté Rv. Le signal d’entrée utilisé dans ce cas est décrit par la relation suivante : ⎧1.6 ⎪1.6 exp((k −100) / 5) ⎪⎪ u(k ) = ⎨0.013 k − 2.5373 ⎪ ⎪1.6 sin(π + k ) ⎪⎩ 2 5π pour k = 1:100 pour k = 101: 200 pour k = 201:300 (2-60) pour k = 301: 400 La figure 2-9(a) représente l’évolution des sorties du système et des modèles considérés ainsi que les erreurs de modélisation pour un rapport Rv=2%. La figure 2-9(b) représente les mêmes résultats pour Rv=10%. Le tableau 2-6 présente les variances de l’erreur de validation pour chaque méthode présentée. Var. Err. (K) Var. Err. (AGB) Var. Err. (RN) Rv=2% 0.02 0.0205 0.0193 Rv=10% 0.1628 0.1650 0.1621 Tableau 2-6 : Variance des erreurs de validation On remarque, d’après les figures 2-9(a) et 2-9(b), que même avec un rapport Rv qui atteint 10%, la sortie du modèle ne diverge pas par rapport à la sortie du système. Le tableau 2-6 montre que les variances des erreurs, pour les trois méthodes, sont faibles. Ce qui valide la robustesse des modèles identifiés. Néanmoins, la variance de l’erreur obtenue par la troisième méthode (RN) est moins faible que les variances obtenues par les deux autres méthodes. Evolution des sorties (système et modèles) 4 système K AGB RNR 3 2 1 0 0 50 100 150 200 250 300 350 400 Evolution des erreurs de validation 0.4 K AGB RNR 0.2 0 -0.2 -0.4 0 50 100 150 200 itérations 250 300 350 400 Figure 2-9(a) : Evolutions des sorties (système et modèles) et des erreurs de validation (système 2, Rv =2%) 57 Chapitre 2 Evolutions des sorties (système et modèles) 4 Système K AGB RNR 3 2 1 0 -1 0 50 100 150 200 250 300 350 400 Evolutions des erreurs de validation 1 AGB AGB RNR 0.5 0 -0.5 -1 0 50 100 150 200 itérations 250 300 350 400 Figure 2-9(b) : Evolutions des sorties (système et modèles) et des erreurs de validation (système 2, Rv =10%) Système 3 : le troisième système est décrit par [Chen, 1999] : x1 (k + 1) = 0.1x1 (k ) + 2 u (k ) + x 2 (k ) 1 + (u (k ) + x 2 (k ) )2 ⎛ ⎞ u ²(k ) ⎟⎟ x 2 (k + 1) = 0.1x 2 (k ) + u ( k )⎜⎜1 + ⎝ 1 + (x1 ²(k ) + x 2 ²(k ) ) ⎠ y ( k ) = x1 (k ) + x 2 (k ). (2-61) Le tableau 2-7 comporte les valeurs de R²mult en fonction des paramètres m, n et q. On constate que la plus grande valeur de R²mult est obtenue lorsque n=2 et m=q=3. Par conséquent, le modèle NARMA considéré est caractérisé par 55 termes. m n q 2 2 3 4 1 2 2 3 3 3 2 1 2 2 2 2 3 3 R²mult 87.5018 87.5018 87.5018 87.5018 91.1783 96.1299 m n q 2 3 3 4 5 5 2 2 4 5 4 5 3 3 3 3 3 3 R²mult 96.0043 96.5732 95.4492 96.4712 94.9868 95.4942 Tableau 2-7 : Evolution de R²mult en fonction des paramètres m, n et q (système 3) 58 Chapitre 2 Le tableau 2-8 comporte les termes sélectionnés, les paramètres estimés, les valeurs de R²tot ainsi que le temps demandé par chaque méthode. Les coefficients du polynôme de la fonction d’activation identifiés par l’algorithme génétique réel sont les suivants : f0=3.35; f1=3.43; f2=-4.83 et f3=-2.05. D’après le tableau 2-8, on constate que la méthode de Kortmann a conduit à un modèle NARMA caractérisé par un nombre réduit des paramètres. De plus, cette méthode nécessite un temps très faible par rapport aux deux autres approches. Le nombre de termes du modèle obtenu par la méthode (RN) peut être réduit en supprimant les termes ayant des paramètres de valeurs très faibles. Termes possibles N° 0 y K AGB RN -0.0444 -0.0600 0.1295 N° Termes possibles K 28 u²(k-1) y(k-1) 1 u(k-1) 2.1089 2.2439 1.0432 29 u²(k-1) y(k-2) -0.071 2 u(k-2) 1.2887 1.0852 -0.0046 30 u(k-1) u(k-2) y(k-1) -0.715 3 u(k-3) 0.2361 0.1517 0.0166 31 u(k-1) u(k-2) y(k-2) 4 y(k-1) 5 y(k-2) 6 u²(k-1) 7 8 RN 0.0222 -0.0103 -0.0685 -0.0341 0.1493 -0.0413 -0.0289 0.2333 32 u(k-1) u(k-3) y(k-1) -0.0488 0.2299 33 u(k-1) u(k-3) y(k-2) -0.0441 0.1078 0.0368 34 u²(k-2) y(k-1) 0.4282 u(k-1) u(k-2) -0.0938 -0.0172 35 u²(k-2) y(k-2) -0.0041 -0.0261 u(k-1) u(k-3) 0.0407 -0.0012 36 u(k-2) u(k-3) y(k-1) 0.0032 -0.2515 9 u²(k-2) 0.3288 -0.3361 37 u(k-2) u(k-3) y(k-2) 0.0134 -0.0144 10 u(k-2) u(k-3) -0.3308 0.1543 38 u²(k-3) y(k-1) -0.0067 0.0169 11 u²(k-3) -0.6170 0.0568 39 u²(k-3) y(k-2) -0.0089 -0.0415 12 u3(k-1) -0.0344 40 u(k-1) y²(k-1) 13 u²(k-1) u(k-2) -0.0056 41 u(k-2) y²(k-1) -0.2774 14 u²(k-1) u(k-3) -0.0364 0.0176 42 u(k-3) y²(k-1) 0.1519 15 u(k-1) u²(k-2) 0.0267 -0.0704 43 u(k-1) y(k-1) y(k-2) -0.0001 16 u(k-1) u(k-2) u(k-3) 0.0488 44 u(k-2) y(k-1) y(k-2) 0.1989 17 u(k-1) u²(k-3) 18 u3(k-1) 19 u²(k-2) u(k-3) 20 u(k-2) u²(k-3) 21 u3(k-3) 22 u(k-1) y(k-1) 23 24 25 u(k-2) y(k-2) 0.0357 AGB 0.1337 0.3172 -0.269 -0.0865 0.0005 45 u(k-3) y(k-1) y(k-2) -0.1228 46 u(k-1) y²(k-2) 0.0092 -0.0371 0.0524 0.1222 47 u(k-2) y²(k-2) 0.0230 0.0738 -0.0119 48 u(k-3) y²(k-2) -0.0043 0.0049 -0.0055 49 y²(k-1) -0.0011 -0.1593 -0.2168 0.1097 50 y(k-1) y(k-2) u(k-1) y(k-2) -0.0346 0.0032 51 y²(k-2) u(k-2) y(k-1) 0.0018 0.3310 52 y3(k-1) -0.0536 -0.0852 53 y²(k-1) y(k-2) -0.0094 26 u(k-3) y(k-1) 27 u(k-3) y(k-2) -0.0904 0.1465 0.0192 0.7010 0.0032 -0.0340 -0.1431 54 y(k-1) y²(k-2) -0.0769 -0.0565 55 y3(k-2) -0.0148 R²mult 0.9657 0.9751 0.9740 Temps (sec) 9.5505 61.501 2665.09 Tableau 2-8 : Les valeurs des paramètres estimés, de R²mult et du temps de calcul (système 3) 59 Chapitre 2 Système 4 : le quatrième système est décrit par [Wen, 2004]: y (k ) = y ( k − 1) y ( k − 2 ) ( y ( k − 1) + 2 . 5 ) + u ( k − 1) 1 + y ²( k − 1 ) + y ²( k − 2 ) (2-62) Dans ce cas, le modèle NARMA estimé est caractérisé par 219 termes (n=5, m=4 et q=3). En effet, ce choix des paramètres m, n et q a conduit à la plus grande valeur de R²mult (tableau 2-9). Les modèles trouvés par les méthodes (K) et (AGB) sont donnés par les équations (2-63) et (2-64). Les paramètres du modèle NARMA estimés par la méthode (RN) sont donnés dans l’annexe 2. Les coefficients du polynôme de la fonction d’activation identifiés par l’algorithme génétique réel sont les suivants : f0=0.41; f1=-3.00; f2=-1.59 et f3=-0.47. m n q R²mult m n q R²mult 222 322 232 332 432 442 113 123 213 223 88.4192 88.6387 90.0294 90.3396 90.3396 88.8969 69.7119 89.8282 79.6418 92.9585 323 233 333 343 433 443 93.9618 94.5875 95.9938 95.0294 96.0186 96.0007 453 96.5307 543 553 95.9317 96.1293 Tableau 2-9 : Evolution de R²mult en fonction des paramètres m, n et q (système 4) - Méthode de Kortmann : yˆ ( k ) = 0.0173 + 0.99u ( k − 1) − 0.1045 u ( k − 3) + 0.1433 y ( k − 1) + 01715 y ( k − 2) − 0.0464 u 2 ( k − 2) + 0.0744 u ( k − 3) y ( k − 2) + 0.3391u ( k − 2) y ( k − 1) y ( k − 2) − 0.2011u ( k − 2) y 2 ( k − 1) + 0.1418 u ( k − 3) y 2 ( k − 1) + 0.1663u 2 ( k − 2) y ( k − 1) − 0.1141u ( k − 3) y ( k − 1) y ( k − 2) (2-63) − 0.0263u ( k − 4) y ( k − 3) y ( k − 4) + 0.6115 y ( k − 1) y ( k − 2) + 0.1598 y ( k − 1) y ( k − 3) − 0.2141 y 2 ( k − 1) y ( k − 2) 60 Chapitre 2 - Méthode AGB : yˆ ( k ) = 0 . 015 + 0 . 9928 u ( k − 1) − 0 . 606 u ( k − 2 ) + 0 . 0078 u ( k − 3 ) + 0 . 158 u ( k − 4 ) + 0 . 59 y ( k − 1) + 0 . 26 y ( k − 2 ) − 0 . 0574 y ( k − 4 ) + 0 . 0331 y ( k − 5 ) + 0 . 1371 u ( k − 2 ) u ( k − 3 ) + 0 . 0677 u ( k − 2 ) u ( k − 4 ) − 0 . 0713 u 2 ( k − 3 ) − 0 . 0062 u 3 ( k − 1) − 0 . 0592 u 2 ( k − 1) u ( k − 3 ) − 0 . 0434 u ( k − 1) u ( k − 2 ) u ( k − 4 ) − 0 . 021 u ( k − 1) u 2 ( k − 3 ) − 0 . 022 u ( k − 1) u ( k − 3 ) u ( k − 4 ) + 0 . 0354 u ( k − 1) u 2 ( k − 4 ) + 0 . 1418 u 2 ( k − 2 ) u ( k − 3 ) − 0 . 0624 u ( k − 2 ) u 2 ( k − 3 ) − 0 . 0438 u 3 ( k − 3 ) + 0 . 0294 u ( k − 2 ) u ( k − 3 ) u ( k − 4 ) − 0 . 1723 u ( k − 3 ) y ( k − 3 ) + 0 . 0254 u ( k − 1) y ( k − 1) + 0 . 3764 u ( k − 2 ) y ( k − 2 ) + 0 . 1217 u ( k − 3 ) y ( k − 2 ) (2-64) Les valeurs de R²mult ainsi que le temps demandé par chaque méthode sont donnés sur le tableau 2-10. D’après ce tableau, on constate que la plus grande valeur de R²mult est assurée par la méthode basée sur les réseaux de neurones artificiels. Cependant, cette méthode exige un temps très élevé pour l’estimation des paramètres du modèle NARMA. K AGB RN R²mult 0.9653 0.9116 0.9707 Temps (sec) 36.589 2198.2 2401.03 Tableau 2-10 : Les valeurs de R²mult et du temps de calcul (système 4) 2.5. Conclusion Dans ce chapitre, nous nous sommes intéressé à l’identification paramétrique et structurelle des modèles de type NARMA. Nous avons présenté en premier lieu l’algorithme de régression par pas échelonné modifié, qui est basé sur des critères statistiques. Ensuite, nous avons proposé deux méthodes heuristiques pour l’identification des modèles de type NARMA. La première méthode est basée sur les algorithmes génétiques binaires et la deuxième méthode constitue une combinaison entre le réseau de neurones artificiels à fonction d’activation polynomiale et l’algorithme génétique sous sa représentation réelle. Cette dernière permet de chercher les coefficients de la fonction polynomiale simultanément avec l’apprentissage du réseau de neurones. Ces algorithmes permettent la sélection des 61 Chapitre 2 termes ainsi que l’estimation des paramètres du modèle NARMA. Les résultats de simulation ont confirmé l’efficacité des algorithmes proposés. En effet, des modèles NARMA caractérisés par une précision acceptable et par une complexité raisonnable ont été déterminés. L'utilisation de telles méthodes est très recommandée quand nous n'avons aucune connaissance a priori de la structure du système. Cet avantage permet à cette approche de traiter les systèmes complexes. D’autre part, l’expression de la sortie du modèle NARMA en fonction des paramètres du réseau neuronal (méthode RN) peut être calculée par un programme approprié. Par conséquent, on évite la difficulté d’implémentation de l’algorithme de régression par pas échelonné modifié surtout dans le cas des systèmes multivariables. Dans le chapitre qui suit, nous allons exploiter les modèles identifiés pour la synthèse d’un correcteur prédictif sous contraintes qui ne nécessite pas le calcul de la dérivée. 62 Chapitre 3 Chapitre 3 La commande prédictive non linéaire Chapitre 3 La commande prédictive non linéaire 3.1. Introduction ............................................................................................................ 64 3.2. Calcul du prédicteur à j - pas ................................................................................ 65 3.3. Le critère de performance ..................................................................................... 66 3.4. Optimisation du critère........................................................................................... 67 3.4.1. Généralités sur les méthodes d’optimisation................................................ 67 3.4.2. La méthode de Nelder-Mead ....................................................................... 71 3.4.3. La méthode de Rosenbrock .......................................................................... 76 3.4.4. Améliorations des performances des méthodes d’optimisation................... 79 3.5. La loi de commande .............................................................................................. 86 3.6. Résultats de simulation .......................................................................................... 87 3.7. Conclusion............................................................................................................ 100 __________________________________________________________________________ 63 Chapitre 3 Chapitre 3 La commande prédictive non linéaire 3.1. Introduction Dans ce chapitre, on s’intéresse à la commande prédictive basée sur un modèle de type ARMA non linéaire. L’utilisation d’un modèle non linéaire conduit à un problème de commande non convexe. La minimisation du critère de performance est généralement effectuée par une approche classique tel que la méthode du gradient, qui nécessite le calcul de la dérivée de ce dernier. Cette méthode d’optimisation donne une loi de commande sous optimale à cause de la présence des minimums locaux [Mrabet, 2002], [Filali, 2002], [Kenneth, 1999], [Yonghong, 1996]. Dans notre travail, la loi de commande est calculée par deux méthodes différentes qui ne nécessitent pas le calcul de la dérivée. La première méthode utilise une version améliorée de l’algorithme de Nelder-Mead, appelé encore méthode du simplexe, et la deuxième exploite l’algorithme de Rosenbrock (méthode de recherche multidirectionnelle). Pour traiter correctement les contraintes et pour avoir plus de probabilité de convergence vers des optimums globaux, nous avons introduit des améliorations aux méthodes d’optimisation utilisées. Tous les développements dans ce chapitre se basent sur une approche discrète de la commande prédictive. On se focalisera essentiellement sur l’application à la Commande Prédictive à base de Modèle MPC. Le chapitre est organisé comme suit : Dans le deuxième paragraphe on présente le calcul du prédicteur. Le troisième paragraphe concerne le critère de performance. Le paragraphe quatre sera consacré au développement des méthodes d’optimisations utilisées ainsi que les améliorations introduites. Le calcul de la loi de commande est présenté dans le cinquième paragraphe. Les résultats de simulation seront donnés dans le paragraphe six et la conclusion dans le dernier paragraphe. 64 Chapitre 3 3.2. Calcul du prédicteur à j - pas : Le calcul de la loi de commande requiert la définition d’un prédicteur qui permet d’anticiper le comportement futur du processus sur un horizon fini. Le calcul du prédicteur yˆ ( k + Hp ) , c'est-à-dire la prédiction effectuée à l’instant k de la valeur de la sortie y à Hp pas d’échantillonnage en avance, nécessite l’exploitation de l’équation du modèle obtenu. La figure 3-1 présente le schéma de principe du prédicteur avec un horizon de prédiction égal à Hp, en utilisant un modèle de type NARMA. y(k-1) y(k-n) u(k-1) yˆ(k) Modèle NARMA u(k-m) Modèle NARMA yˆ(k +1) u(k-m+1) yˆ(k +1) Modèle NARMA yˆ(k + 2) yˆ(k +Hp−1) yˆ(k + Hp) Modèle NARMA u(k-m+2) u(k-m+Hp) Figure 3-1 : Schéma de principe du prédicteur L’équation du prédicteur sera déduite à partir de l’équation du modèle obtenu. La relation (3-1) donne l’expression du prédicteur, déduite à partir de l’équation (1-7) du modèle général de type NARMA. 65 Chapitre 3 yˆ ( k + l ) = y + ∑ bi u ( k − i + l ) + ∑ ∑ bij u ( k − i + l ) u ( k − j + l ) m m i =1 m i =1 j = i + ... + ∑ ... ∑ ∑ bi ... l u ( k − i + l )... u ( k − h + l ) m m i =1 m v= p h=v + ∑ ai yˆ ( k − i + l ) + ∑ ∑ aij yˆ ( k − i + l ) yˆ ( k − j + l ) n n i =1 n i =1 j = i + ... + ∑ ... ∑ ∑ ai ... l yˆ ( k − i + l )... yˆ ( k − h + l ) n n i =1 n v= p h=v + ∑ ∑ cij u ( k − i + l ) yˆ ( k − j + l ) m (3-1) n i =1 j =1 + ... + ∑ ...∑ ∑ ci ... l u ( k − i + l )... u ( k − v + l ) yˆ ( k − l ) m m i =1 n v =1 l =1 + ... + ∑ ∑ ...∑ cij ... l u ( k − i + l ) yˆ ( k − j + l )... yˆ ( k − l ). m n i =1 j =1 n l =1 3.3. Le critère de performance La détermination de la loi de commande prédictive est associée à la minimisation d’un critère quadratique, à horizon fini, qui fait intervenir deux termes quadratiques, le premier est lié aux erreurs futures entre la sortie prédite et la consigne et le deuxième correspond aux incréments futurs de commande. Le critère utilisé pour calculer la loi de commande est donné par la relation suivante : Hc ⎤ 1 ⎡ Hp J 2 (k ) = ⎢∑ (yˆ(k + i) -yc(k + i))2 + λ∑ ( Δu(k + i-1 )) 2 ⎥ 2 ⎣ i=1 i =1 ⎦ (3- 2) avec : Hp: horizon de prédiction sur la sortie, Hc: horizon de commande, λ: coefficient de pondération sur la commande qui permet de donner plus ou moins de poids à la minimisation de l’énergie (commande) par rapport à l’erreur de poursuite. 66 Chapitre 3 On note par ΔU le vecteur des incréments de commande, qui est donné par la relation suivante : ΔU = [Δu (k ) ... Δu (k + Hc − 1)]T (3-3) Puisque le modèle est non linéaire, le critère J2 est non convexe vis-à-vis du vecteur ΔU . Par conséquent, l’utilisation d’une méthode classique d’optimisation, telle que la méthode de descente du gradient, peut conduire à une solution locale. Dans ce qui suit, nous présentons deux méthodes d’optimisation utilisées pour l’obtention de la loi de commande, ces deux méthodes sont d’ordre zéro, i.e. ne nécessitent pas le calcul de la dérivée. 3.4. Optimisation du critère 3.4.1. Généralités sur les méthodes d’optimisation Un facteur très important qui n’est pas traité par la théorie de la commande prédictive mais qui a une importance pratique est celui de l’influence des temps de calculs. Si le temps de calcul de la loi de commande est supérieur au temps d’application (la période d’échantillonnage), on risque d’engendrer une instabilité dans le schéma de commande prédictive. On a intérêt, donc, à réduire le temps de calcul par rapport à la constante de temps du système à commander. Malgré l’arrivée de calculateurs qui sont de plus en plus puissants, cette contrainte reste importante dans le cas des systèmes rapides. Ceci explique le grand succès de la commande prédictive dans le domaine du génie des procédés par exemple et son succès plus que modeste dans d’autres domaines tel que l’aéronautique. En effet, le temps de calcul de la loi de commande dépend aussi bien de la complexité du système que de la méthode utilisée pour la résolution du problème. La formule (3-4) présente la formulation mathématique d’un problème d’optimisation de dimension n : 67 Chapitre 3 ⎧ Min f(x) ∈ ℜ n ⎪ ⎪ g i (x) ≤ 0 i = 1 ,...,p ⎨ ⎪ h j (x) = 0 j = 1 ,...,q ⎪x ⎩ l min ≤ x l ≤ x l max l = 1 ,...,n (3-4) où f(x) est le critère à minimiser appelé aussi fonction objectif x est un vecteur à n variables qui représentent les paramètres du problème à optimiser gi(x) et hj(x) représentent respectivement les contraintes d’inégalité et d’égalité, xl min et xl max désignent les contraintes de domaine. ℜ n est l’espace de recherche borné par les contraintes de domaine. La solution d’un problème d’optimisation est alors donnée par un ensemble de paramètres x* pour lesquels la fonction objectif présente une valeur minimale, en respectant les contraintes d’égalité, d’inégalité et de domaine. Selon les caractéristiques du problème d’optimisation, nous pouvons appliquer différentes méthodes de résolution pour identifier sa solution. Ces méthodes sont classées en deux grands groupes: les méthodes déterministes et les méthodes stochastiques. Dans les méthodes d’optimisation déterministes l’évolution vers la solution du problème est toujours la même pour un même contexte (point de départ). Ce sont en général des méthodes efficaces, peu coûteuses, mais qui nécessitent une configuration initiale (point de départ) pour résoudre le problème. Ce sont souvent des méthodes locales, c’est-à-dire qu’elles convergent vers l’optimum le plus proche du point de départ. Selon la dimension de la fonction objectif à optimiser, les méthodes déterministes peuvent être classifiées en unidimensionnelles ou multidimensionnelles. Parmi les méthodes déterministes unidimensionnelles, aussi appelées méthodes de Recherche Linéaire (Line Search Methods), on distingue la méthode de Dichotomie [Culioli, 1994], la méthode de la Section Dorée [Culioli, 1994] [Press, 1992] et la méthode de Brent [Brent, 1973] [Press, 1992]. La plupart de ces méthodes ne supposent pas que la fonction à minimiser soit différentiable, ni même continue, mais seulement uni modale [Culioli, 1994]. 68 Chapitre 3 Les méthodes déterministes multidimensionnelles, quant à elles, sont consacrées à l’optimisation de fonctions à un paramètre ou plus. Elles peuvent être classées selon l’utilisation de l’information des dérivées de la fonction objectif par rapport aux paramètres xl. Elles sont dites directes ou d’ordre 0 si elles n’utilisent que l’information de la valeur de la fonction elle-même. Dans le cas où elles nécessitent aussi le calcul du gradient de la fonction, elles sont dites indirectes ou d’ordre 1. Les méthodes d’ordre 0 offrent l’avantage de se passer du calcul du gradient, ce qui peut être intéressant lorsque la fonction n’est pas différentiable ou lorsque le calcul de son gradient représente un coût important. Les méthodes multidimensionnelles, peuvent être groupées en deux groupes: les méthodes analytiques ou de descente et les méthodes heuristiques ou géométriques. Les méthodes analytiques se basent sur la connaissance d’une direction de recherche souvent donnée par le gradient de la fonction. Les exemples les plus significatifs de méthodes analytiques sont la méthode de la Plus Grande Pente [Culioli, 1994], le Gradient Conjugué [Culioli, 1994] [Fletcher, 1987] [Press, 1992], la méthode de Powell [Powell, 1965] et les méthodes QuasiNewton [Culioli, 1994] [Fletcher, 1987] [Press, 1992]. Les méthodes heuristiques explorent l’espace par essais successifs en recherchant les directions les plus favorables. À l’opposé des méthodes analytiques, la plupart de ces méthodes sont d’ordre 0. Les implémentations de méthodes géométriques les plus souvent utilisées sont celles de la méthode du Simplexe [Nelder, 1965], la méthode de Rosenbrock [Rao, 1996] et la méthode de variations locales de Hooke et Jeeves [Cherruault, 1999]. Le deuxième groupe qui englobe les méthodes d’optimisation stochastiques s’appuie sur des mécanismes de transition probabilistes et aléatoires. Cette caractéristique indique que plusieurs exécutions successives de ces méthodes peuvent conduire à des résultats différents pour une même configuration initiale d’un problème d’optimisation. Malgré l’aspect du hasard, ces méthodes ont une grande capacité de trouver l’optimum global du problème, contrairement à la plupart des méthodes déterministes. Cependant, elles demandent un nombre important d’évaluations avant d’arriver à la solution du problème. Parmi les méthodes stochastiques les plus employées, nous distinguons le Recuit Simulé [Kirkpatrick, 1983], la Recherche Tabou [Glover, 1989] [Glover, 1990] [Hu, 1992] et les Méthodes Évolutionnistes tel que les Algorithmes Génétiques [Holland, 1975], [Goldberg, 1989]. 69 Chapitre 3 La figure 3-2 présente un organigramme de classification des méthodes d’optimisation les plus utilisées (stochastiques et déterministes). Cette classification n'est pas exhaustive, mais regroupe les principales stratégies d’optimisations employées à l'heure actuelle. Dans ce travail, on propose le calcul de la loi de commande par deux méthodes qui n’utilisent pas la dérivée. La première est basée sur l’algorithme du Simplexe (Nelder-Mead), et la deuxième exploite l’algorithme de Rosenbrock. Ces méthodes sont destinées à des problèmes d’optimisation sans contraintes. Dans nos travaux, nous avons associés à ces méthodes des éléments afin de traiter des fonctions non linéaires, non convexes et en présence des contraintes. Méthodes d’Optimisation Méthodes stochastiques Méthodes déterministes Méthodes Evolutionnistes Méthodes sans contraintes - Algorithmes Génétiques - Stratégies d’Evolution - Programmation Evolutionniste - Programmation Génétique Méthodes avec contraintes Méthodes de Transformation Méthodes Unidimensionnelles Recuit Simulé Recherche avec Tabou - Dichotomie Section Dorée Brent Fibonacci Armijo et Goldstein - Pénalité Lagrangien Augmenté Variables Mixtes Asymptotiques Mobiles Méthodes Multidimensionnelles Méthodes directes Méthodes Analytiques Méthodes Heuristiques - Gradient Conjugué Powell Plus grande pente Quasi Newton - Ellipsoïde Direction Admissibles Gradient Réduit Gradient Projeté - Simplex - Rosenbrock - Hooke et Jeeves Figure 3-2 : Les méthodes d’optimisation les plus importantes 70 Chapitre 3 3.4.2. La méthode de Nelder-Mead L’algorithme de Nelder-Mead, appelé encore algorithme du simplexe, permet de déterminer le minimum d’une fonction sans calcul de la dérivée. Dans notre cas, on a exploité cette propriété pour minimiser le critère de performance J2(k) donné par la relation (3-2). Avant de présenter le principe de l’algorithme de Nelder-Mead, il est nécessaire de donner quelques définitions. 3.4.2.1. • Définitions Un polytope, un polyèdre : est une figure géométrique de (n+1) points dans un espace à n dimensions. • Un simplexe : on appelle simplexe de ℜ n tout ensemble de points ( x1 , x 2 ,..., x n +1 ) de • Une réflexion : consiste à échanger le point de plus grande valeur avec son symétrique ℜ tel que f(x1 ) ≤ f(xi ) ∀ i ∈ [2, L , n + 1] . par rapport au barycentre des n points restants. • Une contraction : consiste à échanger le point de plus grande valeur avec son image par l’homothétie de centre (le barycentre) des n points restants et d’un facteur γ < 1 . On distingue deux cas : - soit une contraction externe, dans ce cas γ = 0.5 . • • - soit une contraction interne, γ = −0.5 . Une expansion : consiste à échanger le point de plus grande valeur avec son image par l’homothétie de centre (le barycentre) des n points restants et d’un rapport χ supérieur à 1 (généralement χ = 2 ). Le rétrécissement ou « multi-contraction » : consiste à échanger tous les points xi (i = 2, ..,n + 1) du polytope par (n) nouveaux points en utilisant la relation suivante : xi = (xi + x1) 2 , (i = 2, .., n + 1) 71 Chapitre 3 3.4.2.2. Méthode du polytope de Nelder-Mead C'est une méthode d'optimisation locale [Nelder, 1965] [Wright, 1996] qui est fréquemment utilisée. Cette méthode déterministe est dite "directe" : elle tente de résoudre le problème en utilisant directement la valeur de la fonction objectif, sans faire appel à ses dérivées. Cette méthode est surtout appréciée pour sa robustesse, sa simplicité de programmation, sa faible consommation de mémoire (peu de variables) et son faible temps de calcul. Comparé à d'autres méthodes de recherche directes, la méthode NM est relativement rapide. Le temps de calcul peut être amélioré suivant l'initialisation du simplexe. Une fois le nombre de points du simplexe est fixé, leurs valeurs influent la vitesse de l'algorithme. Dans [Yang, 2005] les auteurs ont exploité les Algorithmes Génétiques et la méthode NM pour l'optimisation globale, la population initiale de l'Algorithme Génétique devrait éviter le risque d'avoir un grand nombre d'individus dans la même région. D'autre part, puisque la méthode NM n'est pas une méthode d'optimisation globale, l’utilisation des points de départs multiples (multi initialisation) est avantageux et aide à vérifier qu'une valeur fortement optimale sera trouvée [Bortolot, 2005]. Deux autres propriétés de l'algorithme NM sont développées par M. A. Luersen et R. le Riche [Luersen, 2004]. Premièrement, l'algorithme NM peut converger à un optimum local, ce qui arrive en particulier quand le simplexe s'effondre dans un sousespace. Deuxièmement, la méthode peut échapper à une région qui serait un bassin d'attraction si le simplexe est assez large. D’autre part, si la taille du simplexe est petite, l'algorithme devient local. Cependant, cet algorithme est robuste car il est très tolérant aux bruits dans les valeurs de la fonction objectif. En conséquence, la fonction n'a pas besoin d'être calculée exactement et il est possible d'avoir recours à une approximation de la valeur de la fonction. Contrairement aux autres méthodes qui démarrent à partir d'un point initial, la méthode de Nelder-Mead utilise un "polytope" de départ. 3.4.2.3. Principe de la méthode de Nelder-Mead L’algorithme de Nelder-Mead est basé sur l’évaluation de la fonction (le critère J2) aux sommets du simplexe, sans calcul de dérivées [Kelley, 1999] [Wright, 1996]. Il est utilisé, dans notre cas, pour estimer itérativement le minimum du critère de performance donné par l’équation (3-2). 72 Chapitre 3 Pour cela on initialise un polyèdre sur une zone aléatoire, ce polyèdre de départ sera donc obtenu par le tirage aléatoire d’un point x1 dans l’espace solution. Les autres points xi (i = 2, ..,n + 1) sont choisis de manière à former une base, qui est généralement orthogonale. Dans chaque itération de l’algorithme de Nelder-Mead, un simplexe formé de (n+1) points est utilisé pour déterminer un pas d’essai. Chaque point du simplexe correspond à une valeur de la fonction. Les points xi (i = 1,L, n + 1 ) sont ordonnés, au début de l’itération, de manière à avoir: f(x 1 ) ≤ f(x 2 ) ≤ … ≤ f(x n +1 ). (3-5) Après avoir calculer un ou plusieurs points de test et évaluer la valeur de la fonction, de ces points, l’itération actuelle dégénère un ensemble de n+1 sommets qui définis un simplexe différent pour la prochaine itération. Les points de test sont obtenus moyennant quatre opérations possibles, à savoir : réflexion, expansion, contraction et rétrécissement. Chaque opération est associée à un paramètre. Les coefficients associés aux opérations de réflexion, expansion, contraction et rétrécissement sont respectivement notés par ρ , χ , γ et σ . Ces coefficients doivent remplir les conditions suivantes [Kelley, 1999] : ρ > 0, χ > 1, 0 < γ < 1, et 0 < σ < 1. (3-6) D’après C.T. Kelley [Kelley, 1999] et M.H. Wright [Wright, 1996], La séquence standard 1 1⎫ ⎧ associée aux coefficients {ρ , χ , γ , σ } est ⎨1, 2, , ⎬ 2 2⎭ ⎩ (3-7) Dans chaque itération de l’algorithme de Nelder-Mead deux transformations sont possibles : un seul point du simplexe actuel sera modifié. Le nouveau point de test accepté par l’algorithme remplace le plus mauvais point du simplexe actuel ( x n +1 ), et prépare le nouveau simplexe pour l’itération suivante, 73 Chapitre 3 ou bien un ensemble de n nouveaux points sera calculé et remplace les n actuels points xi (i = 2,L , n + 1 ) afin de former un nouveau simplexe pour l’itération suivante. Il s’agit de l’opération de rétrécissement. L’algorithme s’arrête lorsque la différence ( f ( x n +1 ) - f (x 1 )) , qui présente la distance entre les valeurs de la fonction du meilleur et du plus mauvais point, soit inférieure à un certain seuil de précision ε 0 . x3 x3 x1 x2 x x1 x2 x xr xr xe a- réflexion b- expansion x3 x1 x3 x3 xc′ x x2 x1 x x2 x2 x1 xc xr c- contraction externe d- contraction interne e- rétrécissement (Le simplexe original est présenté avec une ligne continue) Figure 3-3 : Les mouvements possibles dans la méthode du polytope de Nelder-Mead dans ℜ 2 La figure 3-3 montre l’effet des opérations de réflexion, de l’expansion, de la contraction (interne et externe) et du rétrécissement pour un simplexe à deux dimensions, soit un triangle 1 1⎫ ⎧ ( x1 , x 2 , x3 ) en utilisant la séquence standard des coefficients ⎨ρ = 1, χ = 2, γ = , et σ = ⎬ . 2 2⎭ ⎩ 74 Chapitre 3 3.4.2.4. L’algorithme de Nelder-Mead L’algorithme de Nelder-Mead, comme décrit par C. T. Kelley [Kelley, 1999] et M. H. Wright [Wright, 1998], [Wright, 1996], est résumé par les étapes suivantes : Etape 1 : Prendre les sommets ( x1 , x 2 ,..., x n +1 ) d’un polyèdre de départ vérifiant f(x 1 ) ≤ f(x 2 ) ≤ … ≤ f(x n +1 ). Etape 2 : Tant que la différence entre f(x n +1 ) et f(x 1 ) est supérieure à un certain seuil de précision ε 0 . 2.1 2.2 calculer x et f r = f(x( ρ )) x= et x( ρ ) = x + ρ ⋅ ( x − x n +1 ) réflexion : si f(x 1 ) ≤ f r < f(x n ). Alors 2.3 1 n ∑ xi , est le centroïde des n meilleurs sommets. n i =1 avec x n +1 = x( ρ ) , aller à 2.7. expansion : si f r < f(x 1 ) alors calculer f e = f(x( χ )) , avec x( χ ) = x + χ ⋅ ( x r − x ) si f e < f r alors x n +1 = x( χ ) et aller à 2.7 sinon x n +1 = x( ρ ) et aller à 2.7. 2.4 contraction externe : si f(x n ) ≤ f r < f(x n +1 ) alors calculer f o = f(x(γ )) si f o ≤ f r alors x n +1 = x(γ ) puis aller à 2.7 ,avec x(γ ) = x + γ ⋅ ( x r − x ) sinon aller à 2.6 2.5 contraction interne si f r ≥ f(x n +1 ) alors calculer f i = f(x(γ )) , avec x(γ ) = x − γ ⋅ ( x − x n +1 ) si f i < f(x n +1 ) alors x n +1 = x(γ ) puis aller à 2.7 sinon aller à 2.6. 2.6 rétrécissement pour i = 2, ..., n + 1 , calculer xi = 2.7 ( xi + x1 ) 2 , puis calculer f ( xi ) donner les sommets ( x1 , x 2 ,..., x n +1 ) tel que f(x 1 ) ≤ f(x 2 ) ≤ … ≤ f(x n +1 ). retour à l’étape 2. 75 Chapitre 3 3.4.3. La méthode de Rosenbrock Cette méthode d'optimisation a été publiée par Rosenbrock en 1960 [Rosenbrock, 1960], l’approche utilisée met en œuvre plusieurs stratégies d'optimisation différentes. Elle a été capable de trouver de bons compromis même avec des objectifs d'optimisation difficiles à atteindre. 3.4.3.1. Présentation de l'algorithme de Rosenbrock Dans le passé cet algorithme a déjà été utilisé avec succès par Rosenbrock pour l'optimisation de fonctions cibles jusqu'à 50 paramètres et plus (pour des filtres digitaux et analogiques ou la conception de circuit électrique). La stratégie est basée sur un algorithme de recherche très stable d'ordre 0 qui n'exige pas de dériver la fonction cible, bien qu'il s'approche de la méthode du gradient. Donc il combine les avantages des stratégies d'ordre 0 et 1. Il s'est avéré que cette approche simple est plus stable que beaucoup d'algorithmes sophistiqués et exige beaucoup moins de calculs de la fonction cible que des stratégies plus élaborées. À cause de cette stabilité inhérente et avec de bonnes heuristiques pour le calcul des longueurs de pas, cet algorithme est même approprié et a déjà prouvé sa valeur pour des problèmes d'optimisation impliquant des fonctions objectif fortement non linéaires et non monotones. La méthode de Rosenbrock (dite aussi ‘’rotation de coordonnées’’) permet de rendre l’exploration indépendante des directions initiales. C'est une méthode itérative qui se décompose en deux étapes répétitives. Une étape d’exploration qui réalise des améliorations successives via des directions privilégiées, et une étape globale qui permet de construire une nouvelle base dans la direction de deux sorties consécutives. La vitesse de convergence de la méthode dépend du choix de la base de départ. De plus, elle nécessite la reconstruction d’une nouvelle base à chaque itération ce qui rend la méthode coûteuse lorsque la dimension du vecteur des variables de décision est grande. Cependant, la simplicité de cette méthode offre une chance réelle à un utilisateur qui n'est pas un expert en optimisation de le comprendre et de fixer et ajuster ses paramètres correctement. Des procédures de contrôle heuristiques ont été mises en œuvre pour ajuster les 76 Chapitre 3 paramètres d'optimisation en cours d'exécution dès que l'optimisation ralentit ou se bloque. Cela permet de réussir des optimisations sans aucune intervention d'utilisateur. 3.4.3.2. Principe de fonctionnement La première itération de cet algorithme est une simple recherche d'ordre 0 dans les directions des vecteurs de base d'un système de coordonnées à n dimensions. En cas de succès, un essai fournissant une nouvelle valeur minimale de la fonction cible, la longueur de pas est augmentée, tandis que dans le cas d'un échec elle est diminuée et on essaye la direction opposée. Une fois qu'un succès a été trouvé et exploité dans chaque direction de la base, le système de coordonnées est réorienté pour que le premier vecteur de base soit dans la direction du gradient. Toutes les longueurs de pas sont réinitialisées et le processus est répété en utilisant le nouveau système de coordonnée. En initialisant les longueurs de pas à des valeurs plutôt grandes, on permet à la stratégie de laisser de côté des optimums locaux et de continuer la recherche de minimums plus globaux. Dans chaque itération de l’algorithme, une exploration itérative dans les n directions de recherche est effectuée. n étant le nombre de variables et les directions de recherche sont orthogonales et linéairement indépendantes. Quand une nouvelle solution est localisée à la fin de l’itération, un nouvel ensemble de directions de recherche d1 ,L , d n sera construit. La procédure est répétée jusqu’à ce que la valeur de la fonction objectif ne change que de moins d’une valeur donnée (précision). Construction des directions de recherche : Soient d1 ,..., d n , des vecteurs linéairement indépendants et dont leurs normes sont égaux à 1. En outre, supposant que les vecteurs sont mutuellement orthogonaux, c. à d., d it d j = 0 , pour i ≠ j . En commençant par le vecteur actuel x k , la fonction objectif est minimisée itérativement dans toutes les directions de recherche. La solution obtenue sera notée x k +1 . En particulier, x k +1 − x k = ∑ j =1 λ j d j , ou λ j est la distance déplacée le long de d j . n La nouvelle direction de recherche d1 ,..., d n est ainsi formée par la procédure de GramSchmidt (Gram-Schmidt orthogonalization), comme suit : 77 Chapitre 3 ⎧ si λ j = 0 ⎧d j ⎪ ⎪ n ⎪a j = ⎨ ⎪ ∑ λ i d i si λ j ≠ 0 ⎪ ⎩ i= j ⎪ ⎪ si j = 1 ⎧a j ⎪ ⎪ j −1 ⎨b j = ⎨ t a − ⎪ ⎪ j ∑ a j d i d i pour i =1 ⎩ ⎪ ⎪ bj ⎪d j = bj ⎪ ⎪ ⎩ ( d j ) j ≥ 2 = d j ; j = 1 ,..., n (3-8) Une particularité intéressante de cette méthode est qu’un minimum local est toujours encadré dès qu’un succès est suivi d’un échec. Quand cela arrive, le point du milieu des 3 derniers points est toujours plus bas que les 2 autres, ainsi on est dans une position favorable pour essayer une interpolation quadratique. C’est peut être la méthode la plus efficace à utiliser dans le cas d’une fonction unidimensionnelle générale. 3.4.3.3. L’algorithme de Rosenbrock L'algorithme de Rosenbrock est résumé dans les étapes suivantes [Bazaraa, 1993] : Étape d'initialisation Soient ε > 0 , un scalaire positif non nul (critère d’arrêt), d1 ,..., d n les directions de recherche et n le nombre de variables. Choisir un point initial x1, y1 = x1, Prendre k=j=1 ; x1 ∈ ℜ n et passer à l'étape principale. Étape principale 1. On cherche λ j la solution du problème suivant : min ( f ( y j + λ d j ) ) y j +1 = y j + λ j d j si j < n, j = j+1, retour à l’étape 1. Sinon aller à l’étape 2. 78 Chapitre 3 2. x k +1 = y n+1 , si x k +1 − x k < ε arrêter ; sinon y 1 = x k +1 , k=k+1, j=1, aller à l’étape 3 ; 3. Construire une nouvelle direction de recherche, ( d1 ,..., d n ) et retour à l'étape 1. Dans l’étape d’initialisation, les vecteurs directions de recherche ( d1 ,..., d n ) sont orthogonaux, linéairement indépendants et de norme égale à 1. Dans le cas où le nombre de variables n est supérieur ou égale à 2, les vecteurs sont initialisés de manière qu’ils forment une matrice identité de dimension n. [d1 d 2 ... d n ]T = In (3-9) Dans le cas monovariable (n = 1) , on aura un seul vecteur, donc une seule direction de recherche ; soit d 1 = [1 0 ] . Dans l’étape 3, les directions de recherche sont calculées suivant l’équation (3-8). 3.4.4. Améliorations des performances des méthodes d’optimisation Bien que les méthodes d’optimisation déterministes présentées précédemment possèdent des avantages, tel que la faible consommation du temps de calcul, la robustesse et elles ne nécessitent pas le calcul du gradient. Ces méthodes ne sont pas directement adaptables à notre problème. En effet, puisque ces méthodes sont locales (la solution n’est pas globale), elles dépendent de l’initialisation et ne tiennent pas compte des contraintes. D’autre part, les systèmes réels peuvent avoir des contraintes d’ordre pratique ou des limitations physiques qu’on doit prendre en considération. Pour exploiter ces méthodes dans le calcul de la loi de commande, on propose l’introduction de quatre améliorations aux algorithmes d’optimisation de Nelder-Mead et de Rosenbrock. Ces améliorations sont détaillées comme suit : 79 Chapitre 3 1) Application de la méthode de pénalité pour que la méthode d’optimisation traite les contraintes. 2) Application de l’approche CFON (Constrained First Objective Next) [Kurpati, 2002] dans l’initialisation du (des) point(s) de départ(s), en considérant un intervalle qui respecte les contraintes. Ce qui permet de réduire le temps de calcul. 3) Application de la notion de multi initialisation, pour avoir plus de probabilité de convergence vers un minimum global. 4) Application de la notion d’élitisme qui consiste à exploiter la dernière solution trouvée. Application de la méthode de pénalité Les méthodes d'optimisation proposées (Simplexe et Rosenbrock) sont capables de traiter des problèmes non linéaires. Cependant, ces méthodes ne sont pas conçues pour résoudre des problèmes avec des contraintes. L’emploi d’une fonction de pénalité exacte, avec les méthodes d’optimisation proposées, nous permet d’atteindre cet objectif en transformant un problème non contraint en un sous contraintes [Gesing, 1979] [Wright, 1998]. Dans ce cas, les contraintes sont placées dans la fonction objectif par l'intermédiaire d'un paramètre de pénalité de manière qu'il pénalise chaque violation des contraintes. En général, une fonction de pénalité appropriée doit avoir une pénalité positive pour les points infaisables et aucune pénalité pour les points faisables [Bazaraa, 1993]. Ainsi, une fonction de pénalité est généralement traitée comme suit : Si on considère le problème d’optimisation de l’équation (3-4) dans ℜ , alors la fonction de pénalité sera de la forme suivante : α (x) = avec ∑ p i =1 φ [g i ( x ) ] + ∑ ψ [h q i =1 i ( x )] (3-10) φ et ψ sont des fonctions continues qui vérifient les conditions suivantes : φ ( y) = 0 si y ≤ 0 et φ ( y) > 0 si y > 0 ψ ( y) = 0 si y = 0 et ψ ( y) > 0 si y ≠ 0 (3-11) 80 Chapitre 3 φ et ψ peuvent s’écrire comme suit : φ ( y ) = [Max{0, y}]n et ψ ( y ) = y n (3-12) ou n est un entier positif (généralement n=2). La forme de la fonction de pénalité sera alors : α (x) = ∑ [Max {0 , g p i=1 ( x ) }] + n i ∑ q i =1 hi ( x ) n (3-13) Le problème avec contraintes sera remplacé par un problème sans contraintes, en considérant la fonction auxiliaire f ( x ) + μ α ( x ) , ( μ > 0 ). Dans notre cas, la minimisation du critère J2 est faite en respectant les contraintes sur le signal de commande et sur ses incréments. Ces contraintes peuvent être écrites comme suit : u min ≤ u ( k + j ) ≤ u max , Δ u min ≤ Δ u ( k + j ) ≤ Δ u max , j = 0 , ... , Hc − 1 (3-14) j = 0 , ... , Hc − 1 (3-15) Le nombre de contraintes est égal à 4 Hc. Ces contraintes peuvent être réécrites de la manière suivante : Ω = {Δ U / f i ( Δ U ) ≤ 0 , i = 1, ... , 4 Hc } (3-16) Le but est de concevoir, à chaque période d’échantillonnage, une loi de commande qui minimise le critère de performance J2 (équation 3-2) en tenant compte des contraintes sur la commande et sur ses incréments (équation 3-16). Donc la fonction de pénalité sera comme suit : α (x) = ∑ [Max {0 , 4 Hc i =1 f i ( x ) }] 2 (3-17) 81 Chapitre 3 Algorithme de la méthode de pénalité Un sommaire de l'algorithme de la méthode de pénalité est donné comme suit [Bazaraa, 1993] : Étape d'initialisation Soit ε > 0 , un scalaire positif non nul (critère d’arrêt). Choisir un point initial x1, un paramètre de pénalité ν 1 > 0 et une grandeur scalaire β > 0 . Prendre k=1 et passer à l'étape principale. Étape principale 1. On cherche xk, solution du problème suivant : min ( f ( x) + ν k α ( x) ) tel que : x ∈ X . Prendre xk+1 une solution optimale et passer à l'étape 2. 2. Siν k α (x) < ε , arrêter ; sinon, prendreν k +1 = βν k , remplacer k par k+1 et retour à l'étape 1. Application de l’approche CFON L’approche CFON consiste à prendre le(s) point(s) d’initialisation de manière à satisfaire les contraintes dès le départ. Cette approche est employée dans [Kurpati, 2002], pour améliorer la manipulation des contraintes en utilisant les algorithmes génétiques multi objectif. Elle permet de réduire le temps de calcul, puisque tous les points d’initialisation respectent toutes les contraintes dès le début. Ainsi, l’emploi d’une fonction de pénalité avec cette approche évite l'apparition des solutions non faisables dans les boucles intérieures de l'algorithme d’optimisation. Multi-initialisation des algorithmes d’optimisation Cette étape, intégrée dans les algorithmes d’optimisation, permet d’effectuer l’optimisation en partant de différents points d’initialisation. Tous les points de départ respectent l’approche CFON. Après un certain nombre de répétitivité de l’optimisation, on garde le meilleur optimum. Cette étape permet d’augmenter la probabilité de convergence vers le minimum globale. Pour vérifier l’influence de cette approche, on a appliqué trois versions d’initialisation sur des fonctions benchmarks avec des minimums globaux connus. Les fonctions benchmarks (F1, F2, F3 et F4) sont données dans l’annexe 3 [Chelouah, 2000]. La première version proposée consiste à prendre une seule initialisation dans l’intervalle de recherche S (figure 34), la deuxième consiste à prendre plusieurs initialisations aléatoires dans le même intervalle S 82 Chapitre 3 (figure 3-4) et la troisième consiste à diviser l’intervalle S (espace de recherche) en plusieurs sous intervalles égaux (figure 3-5). La résolution du problème d’optimisation s’effectue chaque fois avec une initialisation appartenant à un sous intervalle, de façon à balayer tout l’intervalle de recherche S. A la fin, on garde la meilleure solution parmi les solutions partielles. Dans les figures 3-4 et 3-5, on considère un problème d’optimisation simple, avec deux contraintes c1 et c2. Dans ce cas l’intervalle de recherche S sera donc compris entre c1 et c2. S ℜ c1 c2 Figure 3-4 : Initialisation des algorithmes d’optimisation versions (NM-V1 et NM-V2 ou ROS-V1 et ROS-V2) S s1 c1 … sj … sn ℜ c2 Figure 3-5 : Initialisation des algorithmes d’optimisation versions (NM-V3 ou ROS-V3) Afin d’évaluer l'efficacité de ces méthodes, nous avons considéré des critères récapitulant les résultats de 100 minimisations par fonction. Les critères maintenus sont : le taux de minimisations réussies (%), la moyenne du temps de calcul dans chaque méthode d'optimisation et la moyenne de l'erreur quadratique (variance) entre la solution trouvée par l'optimisation et l’optimum global connu. L'optimisation est considérée réussie si l'erreur quadratique entre la solution obtenue et le minimum global est moins de 5 10-4. Les tableaux 3-1, 3-2 et 3-3 présentent respectivement les résultats obtenus avec les versions originales des algorithmes de Nelder-Mead et de Rosenbrock (NM-V0 et ROS-V0), les trois versions améliorées (NM-V1, NM-V2 et NM-V3) en exploitant l'algorithme de Nelder Mead ainsi que les trois versions (ROS-V1, ROS-V2 et ROS-V3) qui exploitent l'algorithme de Rosenbrock. 83 Chapitre 3 Fonctions Benchmarks (tests) F1 F2 F3 F4 Méthodes temps de taux de minimisations calcul moyen (s) réussies (%) NM-V0 0 ROS-V0 38 NM-V0 0 ROS-V0 50 NM-V0 0 ROS-V0 88 NM-V0 0 ROS-V0 100 0.0395 0.0039 0.2819 0.0044 0.0402 0.0052 0.0683 0.0250 erreur quadratique moyenne 136.4803 68.6517 5.8859 1.7130 0.3723 0.2985 0.0076 4.23 10-9 Tableau 3-1 : Résultats de l’optimisation avec versions originales (méthode NM et ROS) Fonctions Benchmarks (tests) F1 F2 F3 F4 Méthodes temps de taux de minimisations calcul moyen (s) réussies (%) erreur quadratique moyenne NM-V1 25 0.0547 0.0173 NM-V2 82 2.7533 6.252 10-4 NM-V3 100 4.7341 1.479 10-5 NM-V1 8 0.0421 0.161 NM-V2 94 3.0414 1.596 10-5 NM-V3 100 3.8663 5.04110-7 NM-V1 1 0.0020 0.148 NM-V2 27 0.0681 0.007 NM-V3 87 0.0742 2.68 10-5 NM-V1 100 0.0023 4.388 10-4 NM-V2 100 0.0253 4.388 10-5 NM-V3 100 0.0489 1.799 10-5 Tableau 3-2 : Résultats de l’optimisation avec réinitialisation (méthode NM) 84 Chapitre 3 Fonctions benchmarks F1 F2 F3 F4 temps de taux de Méthodes minimisations calcul moyen réussies (%) (s) erreur quadratique moyenne ROS -V1 44 0.0048 48.6253 ROS -V2 99 0.0025 1.3635 ROS -V3 100 0.0022 1.0396 10-8 ROS -V1 56 0.0042 1.5258 ROS -V2 100 0.0028 3.0969 10-9 ROS -V3 100 0.0016 3.4586 10-9 ROS -V1 74 0.0391 0.3563 ROS -V2 88 0.0358 0.2335 ROS -V3 100 0.0363 0 ROS -V1 100 0.0207 4.9297 10-9 ROS -V2 100 0.0222 4.6731 10-9 ROS -V3 100 0.0259 4.0415 10-9 Tableau 3-3 : Résultats de l’optimisation avec réinitialisation (méthode ROS) A partir des tableaux (3-1, 3-2 et 3-3), nous notons que la troisième version des algorithmes NM et ROS (NM-V3 et ROS-V3) présente les meilleurs taux de minimisation réussie, avec une moyenne de l'erreur quadratique la plus basse. On remarque aussi que le temps de calcul est relativement faible. Cependant, la méthode ROS est plus rapide que la méthode NM, ceci est dû au fait que la méthode ROS utilise un seul point de recherche, par contre la méthode NM utilise un ensemble de points qui forment un simplexe. Contrairement à la méthode ROS, les améliorations effectuées sont nettement visibles sur la méthode NM. D’après le tracé de la fonction F3 (annexe 3), cette fonction présente le cas le plus complexe avec presque 50 minimums locaux. Les résultats trouvés dans les tableaux 3-2 et 3-3 prouvent que les deux méthodes NM et ROS combinées avec la version V3 sont capables de trouver des bonnes performances, puisque l’erreur est très faible dans les deux cas. 85 Chapitre 3 3.5. La loi de commande : L’algorithme qui permet la détermination de la loi de commande, en utilisant la méthode de Nelder-Mead, est résumé par les étapes suivantes : 1. Calcul de la sortie du système y(k). 2. Calcul de la sortie yˆ ( k ) à l’aide de l’équation du modèle. 3. Calcul de la valeur uopt, qui permet la minimisation du critère de performance J(k), en utilisant l’algorithme de Nelder-Mead ou bien la méthode de Rosenbrock. 4. Poser u (k+i) = uopt, pour i = 0,…, Hp. 5. Incrémenter k, retour à l’étape 1. L’initialisation des algorithmes d’optimisation en exploitant l’approche CFON, impose que les valeurs initiales de la commande doivent respectées les contraintes données par les relations (3-14) et (3-15). La figure 3-6, décrit le principe général d’exploitation de cette approche pour le calcul de la loi de commande. Espace de recherche ℜ u(k-1) min[u(k-1) + ∆umin , umin] min[u(k-1) + ∆umax , umax] Figure 3-6 : Initialisation des algorithmes d’optimisation (approche CFON) dans ℜ Dans le cas des versions NM-V1 et ROS-V1, on considère une seule initialisation des variables dans le domaine de recherche S. Dans le cas des versions NM-V2 et ROS-V2, on prend plusieurs initialisations différentes de domaine de recherche S. Cependant, dans le cas des algorithmes NM-V3 et ROS-V3, on divise l’espace de recherche S en plusieurs sous 86 Chapitre 3 intervalles puis on considère une initialisation par sous intervalle. Les figures 3-7 et 3-8 présentent, respectivement, l’espace de recherche S considéré avec les versions V1, V2 et V3. S min[u(k-1) + ∆umin , umin] min[u(k-1) + ∆umax , umax] Figure 3-7 : Initialisation des algorithmes d’optimisation versions (NM-V1 et NM-V2 ou ROS-V1 et ROS-V2) S s1 … sj … min[u(k-1) + ∆umin , umin] sn min[u(k-1) + ∆umax , umax] Figure 3-8 : Initialisation des algorithmes d’optimisation versions (NM-V3 ou ROS-V3) 3.6. Résultats de simulation : Système 1: Nous considérons le système non-linéaire décrit par l'équation suivante [Narendra, 1990]: y(k) = y(k −1) + u2 (k −1) + d (k) 2 1+ y (k −1) (3-18) Le signal de perturbation d(k) est défini comme suit : ⎧1 d (k ) = ⎨ ⎩0 k ∈ [120 , ... ,150 ] ailleurs (3-19) 87 Chapitre 3 Pour (m=n=q=2), le modèle identifié en exploitant l'algorithme génétique binaire est comme suit : ym(k ) = 0.2632− 0.0121u(k −1) − 0.0165u(k − 2) + 0.3559 y(k −1) + 0.9991u(k −1) 2 − 0.2316u(k − 2) 2 + 0.0136u(k −1) y(k −1) − 0.0388 y(k −1) 2 (3-20) Les paramètres du correcteur prédictif considéré dans la simulation sont, Hc = 1, Hp = 4 et λ = 0.2. Les contraintes sur la commande et ses incréments sont comme suit : 0 ≤ u (k ) ≤ 3 (3-21) − 0.5 ≤ Δ u ( k ) ≤ 0.5 (3-22) Dans ces simulations, l'intervalle S de l'espace de recherches est donné par : S = [u ( k − 1) + Δ u min , u ( k − 1) + Δ u max ] (3-23) L’approche CFON est considérée. Si on considère la méthode NM (respectivement ROS), alors le simplexe (point) initial est ainsi choisi de manière à satisfaire les contraintes sur l'entrée et sur ses incréments. Nous avons considéré les trois versions (V1, V2 et V3) d'optimisation basées sur l'algorithme de NM (ROS) et la fonction de pénalité. Dans la première méthode, un seul simplexe (point) initial est choisi dans l'intervalle de recherches S (NM-V1, ROS-V1). Dans la deuxième version (NM-V2, ROS-V2), une initialisation multiple et aléatoire de simplexes (points) est employée. Dans ce cas chaque simplexe (point) se trouve sur l'intervalle S du domaine de recherches. Dans la troisième version (NM-V3, ROS-V3), l'intervalle S est divisé en 10 sous intervalles égaux. La résolution du problème d'optimisation est effectué en employant l'algorithme NM (ROS) avec la fonction de pénalité dans chaque sous intervalle. Quand tous les sous intervalles sont considérés, on aura un ensemble de solutions partielles. Nous récupérons, ainsi, la meilleure solution parmi les solutions partielles. Avant de donner les résultats obtenus par les trois versions des algorithmes NM et ROS, nous présentons les résultats de simulations obtenues avec la version originale de chacun des algorithmes NM et ROS (NM-V0 et ROS-V0). La figure 3-9 présente les résultats trouvés avec la méthode NM et la figure 3-10 présente les résultats trouvés avec la méthode ROS. 88 Chapitre 3 6 sortie y consigne yc perturbation d 4 2 0 0 20 40 60 80 100 120 140 160 180 200 140 160 180 200 140 160 180 200 commande u 3 2 1 0 -1 0 20 40 60 80 100 120 Incréments de commande 0.5 0 -0.5 0 20 40 60 80 100 120 Itérations k Figure 3-9 : Sortie du système, commande et incréments de la commande (NM-V0) 6 sortie y consigne yc perturbation d 4 2 0 0 20 40 60 80 100 120 140 160 180 200 140 160 180 200 140 160 180 200 commande u 5 0 -5 0 20 40 60 80 100 120 incréments de commande 1 0 -1 0 20 40 60 80 100 120 itérations k Figure 3-10 : Sortie du système, commande et incréments de la commande (ROS-V0) 89 Chapitre 3 Les évolutions de la consigne yc, de la sortie y, de la commande u et de ses incréments ∆u, obtenues par les trois versions des méthodes NM et ROS sont données dans les figures (3-11,…, 3-16). 6 sortie y consigne yc perturbation 4 2 0 0 20 40 60 80 100 120 140 160 180 200 120 140 160 180 200 140 160 180 200 commande u 3 2 1 0 -1 0 20 40 60 80 100 Incréments de commande 0.5 0 -0.5 0 20 40 60 80 100 Itérations k 120 Figure 3-11 : Sortie du système, commande et incréments de la commande (NM-V1) 6 sortie y consigne yc perturbation d 4 2 0 0 20 40 60 80 100 120 140 160 180 200 140 160 180 200 140 160 180 200 commande u 5 0 -5 0 20 40 60 80 100 120 incréments de commande 1 0 -1 0 20 40 60 80 100 120 itérations k Figure 3-12 : Sortie du système, commande et incréments de la commande (ROS-V1) 90 Chapitre 3 6 sortie y consigne yc perturbation 4 2 0 0 20 40 60 80 100 120 140 160 180 200 140 160 180 200 140 160 180 200 commande u 3 2 1 0 -1 0 20 40 60 80 100 120 Incréments de commande 0.5 0 -0.5 0 20 40 60 80 100 120 Iterations k Figure 3-13 : Sortie du système, commande et incréments de la commande (NM-V2) 6 sortie y consigne yc perturbation 4 2 0 0 20 40 60 80 100 120 140 160 180 200 140 160 180 200 140 160 180 200 commande u 3 2 1 0 -1 0 20 40 60 80 100 120 Incréments de commande 0.5 0 -0.5 0 20 40 60 80 100 120 Iterations k Figure 3-14 : Sortie du système, commande et incréments de la commande (ROS-V2) 91 Chapitre 3 6 sortie y consigne yc perturbation 4 2 0 0 20 40 60 80 100 120 140 160 180 200 140 160 180 200 140 160 180 200 commande u 3 2 1 0 -1 0 20 40 60 80 100 120 Incréments de commande 0.5 0 -0.5 0 20 40 60 80 100 120 Iterations k Figure 3-15 : Sortie du système, commande et incréments de la commande (NM-V3) 6 sortie y consigne yc perturbation d 4 2 0 0 20 40 60 80 100 120 140 160 180 200 140 160 180 200 140 160 180 200 commande u 5 0 -5 0 20 40 60 80 100 120 incréments de commande 1 0 -1 0 20 40 60 80 100 120 itérations k Figure 3-16 : Sortie du système, commande et incréments de la commande (ROS-V3) 92 Chapitre 3 Les figures 3-17 et 3-18 représentent respectivement les résultats obtenus avec la méthode NM-V3 et la méthode ROS-V3. Dans ce cas, les incréments de commande sont limités entre -0.1 et 0.1. 6 sortie y consigne yc perturbation 4 2 0 0 20 40 60 80 100 120 140 160 180 200 120 140 160 180 200 140 160 180 200 commande u 3 2 1 0 -1 0 20 40 60 80 100 Increments de commande 0.2 0 -0.2 0 20 40 60 80 100 Itérations 120 Figure 3-17 : Sortie système, commande et incréments de commande avec |∆u| ≤ 0.1 (NM-V3) 6 sortie y 4 consigne yc perturbation 2 0 0 20 40 60 80 100 120 140 160 180 200 120 140 160 180 200 140 160 180 200 commande u 3 2 1 0 -1 0 20 40 60 80 100 Incréments de commande 0.2 0.1 0 -0.1 -0.2 0 20 40 60 80 100 Iterations k 120 Figure 3-18 : Sortie système, commande et incréments de commande avec |∆u| ≤ 0.1 (ROS-V3) 93 Chapitre 3 Les résultats trouvés, en simulation, confirment les résultats obtenus avec les fonctions benchmark. En effet, avec la méthode NM : la version NM-V3 (figure 3-15) offre les meilleurs résultats trouvés. La version NM-V2 (figure 3-13) donne des résultats très proches de ceux trouvés avec NM-V3, offrant un signal de commande stable. L’effet des améliorations introduites est nettement visible sur la loi de commande et par conséquent sur la qualité du signal de sortie obtenue. La version originale NM-V0 (figure 3-9) est incapable de donner une loi de commande stable permettant de ramener le signal de sortie y(k) à la consigne désirée yc(k). La version NM-V1 (figure 3-11) permet de réaliser cet objectif, cependant, le signal de commande présente beaucoup de fluctuations ce qui entraîne une variance importante du signal de commande. La méthode ROS offre des résultats très intéressants. En effet, d’après les figures 3-10, 312, 3-14, 3-16 et 3-18, toutes les versions (ROS-V0, ROS-V1, ROS-V2 et ROS-V3) ont aboutis à des résultats comparables avec ceux obtenus avec la version NM-V3. Cependant, dans le cas des systèmes plus complexes, les deux dernières versions offrent plus de succès. Le traitement de la fonction benchmark (F3) de l’annexe 3, qui présente plusieurs minimums, montre que les versions ROS-V2 et ROS-V3 sont beaucoup plus performantes que la première version. Le tableau 3-4 donne les valeurs des critères de performances considérés (I1, I2 et I3 , Eq. 3-24 au 3-26) obtenues par les deux méthodes d'optimisation, avec les trois versions considérées. I1 = 1 N 2 ∑ (yc (k + i) - y(k + i)) N i =1 I2 = I3 = 1 N 2 ∑ u (i ) N i =1 1 N 2 ∑ Δu(i) N i =1 (3-24) (3-25) (3-26) 94 Chapitre 3 Méthode I1 I2 I3 0.3173 2.4795 0.0092 NM-V1 Nelder-Mead Rosenbrock NM-V2 0.2751 2.4979 0.0070 NM-V3 0.2745 2.4980 0.0070 ROS-V1 0.2625 2.4893 0.0066 ROS-V2 0.2625 2.4893 0.0066 ROS-V3 0.2625 2.4893 0.0066 Tableau 3-4 : Critères de performances A partir du tableau 3-4, on remarque que les trois versions ROS-V1, ROS-V2 et ROS-V3 de la méthode de Rosenbrock convergent vers la même valeur de la commande u. Cependant, les résultats des versions NM-V3 et NM-V2 sont comparables et leurs performances sont meilleures que la version NM-V1. Les résultats obtenues par la méthode de Rosenbrock sont légèrement meilleurs que ceux obtenus par la méthode de Nelder-Mead (NM-V2 et NM-V3). Dans le tableau 3-5, on présente le temps requis par chaque méthode pour calculer la loi de commande. Ces résultats sont obtenus en effectuant 20 simulations aléatoires avec un PC Pentium 4-2.4GHz. Méthode tmin Temps (s) tmoy tmax NM-V1 ROS-V1 0.0457 0.0784 0.0471 0.08 0.0488 0.084 NM-V2 ROS-V2 0.4482 0.0859 0.4546 0.1261 0.4582 0.1284 NM-V3 ROS-V3 0.4456 0.1282 0.4528 0.1307 0.4583 0.1311 Tableau 3-5 : Temps de simulation selon la méthode d'initialisation Le tableau 3-5, indique, respectivement, le temps minimum, le temps maximum et le temps moyen (tmin, tmax et tmoy) mis dans les 20 simulations aléatoires. On constate que la méthode de Rosenbrock est plus rapide que celle de Nelder-Mead. Ceci s’explique par le fait que la première méthode utilise un seul point de recherche, contrairement à la seconde qui exploite un ensemble de point (un simplexe). La taille du simplexe dans la deuxième méthode 95 Chapitre 3 a une influence sur le temps de calcul, puisque dans chaque itération il y a évaluation de tous les points du simplexe. Nous notons, aussi, que la première version V1 de chacune des deux méthodes est plus rapide que les deux autres, mais les résultats fournis en boucle fermée ne sont pas satisfaisants. La deuxième et la troisième version nécessitent presque le même temps de calcul. Cependant, la troisième version fournit toujours des performances plus satisfaisantes en boucle fermée. Système 2 : le deuxième système est décrit par [Chen, 1999] : x1 (k + 1) = 0.1x1 (k ) + 2 u (k ) + x 2 (k ) 1 + (u (k ) + x 2 (k ) ) 2 ⎛ ⎞ u ²(k ) ⎟⎟ x 2 (k + 1) = 0.1x 2 (k ) + u (k )⎜⎜1 + 1 ²( ) ²( ) x k x k + + 1 2 ⎝ ⎠ y (k ) = x1 (k ) + x 2 (k ). (3-27) Le modèle identifié est celui donné par l'algorithme génétique binaire (méthode AGB), puisqu’il présente la plus grande valeur de R²tot. Le modèle NARMA (m=3, n=2 et q=3) est donné par l’équation suivante : y ( k ) = −0.0444 + 2.1089 u ( k − 1) + 1.2887 u ( k − 2) + 0.2361 u ( k − 3) + 0.1337 u 2 ( k − 1) + 0.3172 u 3 ( k − 1) − 0.269 u 3 (k − 1) − 0.071 u 2 (k − 1) y (k − 2) − 0.715 u (k − 1) u (k − 2) y (k − 1) + 0.1465 u (k − 1) y 2 (k − 1) (3-28) Les paramètres du correcteur prédictif considéré dans la simulation sont, Hc = 1, Hp = 4 et λ = 0.2. Les contraintes sur la commande et ses incréments sont comme suit : − 0 .1 ≤ u ( k ) ≤ 0 .5 (3-29) − 0 . 05 ≤ Δ u ( k ) ≤ 0 . 05 (3-30) Les figures 3-21 et 3-22 comportent, respectivement, les évolutions de la sortie, du signal de commande ainsi que les incréments de commande obtenues par la méthode de NelderMead et la méthode de Rosenbrock. Dans cette simulation, nous avons considéré la version 3 des méthodes d’optimisation. 96 Chapitre 3 2 sortie y consigne yc perturbation 1.5 1 0.5 0 0 20 40 60 80 100 120 140 160 180 200 120 140 160 180 200 140 160 180 200 commande u 0.4 0.2 0 0 20 40 60 80 100 Incréments de commande 0.1 0.05 0 -0.05 -0.1 0 20 40 60 80 100 Itérations k 120 Figure 3-19 : Sortie système (S2), commande et incréments de commande (méthode NM-V3) 2 sortie y consigne yc perturbation 1.5 1 0.5 0 0 20 40 60 80 100 120 140 160 180 200 120 140 160 180 200 140 160 180 200 commande u 0.4 0.2 0 0 20 40 60 80 100 Incréments de commande 0.1 0.05 0 -0.05 -0.1 0 20 40 60 80 100 Iterations k 120 Figure 3-20 : Sortie système (S2), commande et incréments de commande (méthode ROS-V3) 97 Chapitre 3 A partir des figures 3-19 et 3-20, on remarque que les deux méthodes NM-V3 et ROS-V3 sont capables de fournir une loi de commande stable permettant de ramener la sortie du système à la consigne désirée. Cependant, la différence entre les deux méthodes se manifeste lors des changements de la valeur de la consigne ou bien en présence de perturbations. Ceci est remarquable dans la figure 3-19. En effet, entre les itérations 100 à 120 et 140 à 160, le signal des incréments de commande présente des fluctuations contrairement à la méthode ROS-V3 (figure 3-20) qui présente un signal plus stable. Le tableau 3-6 présente les critères de performances (I1, I2 et I3) ainsi que le temps de calcul moyen par itération des algorithmes NM-V3 et ROS-V3. Critères de performances Méthode I1 I2 I3 tmoy NM-V3 0.0129 0.0351 3.019 10-4 0.0161 ROS-V3 0.0115 0.0354 2.516 10-4 0.0175 Tableau 3-6 : Critères de performances (système S2) Le tableau 3-6 montre que les résultats obtenus avec les méthodes NM-V3 et ROS-V3 sont comparables de points de vue stabilité, énergie de commande et temps de calcul. La méthode ROS-V3 donne une valeur plus faible du critère de performance I3. Ceci s’explique par la présence des fluctuations au niveau du signal des incréments de commande avec la méthode NM-V3. Les figures 3-21 et 3-22 présentent le comportement du système en présence de bruit. Le bruit qui a été ajouté à la sortie, est un signal aléatoire d’amplitude égale à 0.3. 98 Chapitre 3 2 sortie y consigne yc perturbation 1.5 1 0.5 0 0 20 40 60 80 100 120 140 160 180 200 120 140 160 180 200 140 160 180 200 commande u 0.4 0.2 0 0 20 40 60 80 100 Incréments de commande 0.1 0.05 0 -0.05 -0.1 0 20 40 60 80 100 Itérations k 120 Figure 3-21 : Sortie système (S2), commande et incréments de commande (méthode NM-V3) avec bruit 2 sortie y consigne yc perturbation 1.5 1 0.5 0 0 20 40 60 80 100 120 140 160 180 200 120 140 160 180 200 140 160 180 200 commande u 0.4 0.2 0 0 20 40 60 80 100 Incréments de commande 0.1 0.05 0 -0.05 -0.1 0 20 40 60 80 100 Itérations k 120 Figure 3-22 : Sortie système (S2), commande et incréments de commande (méthode ROS-V3) avec bruit 99 Chapitre 3 D’après les figures 3-21 et 3-22, on constate que la commande calculée avec les deux méthodes NM-V3 et ROS-V3 est capable de rejeter le signal de perturbation d(k), dans un temps raisonnable. Cette commande permet de ramener la sortie à la consigne désirée, même en présence de bruit non négligeable. 3.7. Conclusion Dans ce chapitre, un contrôleur prédictif des systèmes non-linéaires sous contraintes a été présenté. Le modèle de type NARMA est employé pour caractériser le comportement du système. La loi de commande est obtenue en minimisant un critère quadratique non convexe. Le problème d'optimisation est résolu par deux méthodes utilisant respectivement l'algorithme de Nelder-Mead et l’algorithme de Rosenbrock. Ces méthodes, qui sont déterministes, locales et sans contraintes ont étés améliorées par l’introduction d’une fonction de pénalité, de l’approche CFON ainsi que l’utilisation de la notion de multi initialisation, ce qui favorise la convergence vers le minimum global. La qualité du signal de commande obtenu avec ces méthodes est très bonne, comparée aux résultats obtenus avec les versions originales. Le temps de calcul de la loi de commande qui est relativement faible ainsi que la bonne précision qu’offre la méthode de Rosenbrock améliorée nous encourage à traiter des systèmes plus complexes et même à dynamiques rapides. Comparée aux méthodes d’optimisation qui utilisent la dérivée, les méthodes que nous avons proposées ont l’avantage d’êtres exploitées facilement dans le cas des systèmes multivariables. Dans le chapitre suivant nous nous intéresserons à l’identification et à la commande des systèmes multivariables en exploitant les modèles de type NARMA. 100 Chapitre 4 Chapitre 4 Identification et commande prédictive des systèmes non linéaires multivariables Chapitre 4 Identification et commande prédictive des systèmes non linéaires multivariables 4.1. Introduction .......................................................................................................... 102 4.2. Le modèle NARMA (cas multivariable).............................................................. 103 4.3. Identification du modèle NARMA....................................................................... 105 4.4. Commande prédictive multivariable .................................................................... 109 4.4.1. Calcul du prédicteur .................................................................................. 109 4.4.2. Le critère de performance .......................................................................... 109 4.4.3. Optimisation du critère............................................................................... 110 4.5. Résultats de simulation......................................................................................... 113 4.6. Conclusion............................................................................................................ 124 ___________________________________________________________________ 101 Chapitre 4 Chapitre 4 Identification et commande prédictive des systèmes non linéaires multivariables 4.1. Introduction Dans ce chapitre, on s’intéressera à l’identification et la commande prédictive des systèmes multivariables en utilisant les modèles de type NARMA. Dans le deuxième chapitre, trois approches ont été développées pour la détermination des modèles de type NARMA. La première approche est basée sur l’algorithme de régression par pas échelonné [Kortmann, 1988], les deux autres méthodes exploitent respectivement les algorithmes AGB et RN. L’analyse des trois approches, présentées dans le deuxième chapitre, permet de conclure que la dernière méthode (RN) peut être développée pour traiter les systèmes multivariables. La méthode proposée dans ce chapitre constitue, donc, une combinaison du réseau de neurones artificiels à fonction d’activation polynomiale avec l’algorithme génétique sous sa représentation réelle. Dans ce cas, les individus de la population de l’algorithme génétique représentent les coefficients du polynôme de la fonction d’activation. Les paramètres du modèle NARMA sont alors estimés à partir des pondérations des connexions du réseau neuronal à la fin de la phase d’apprentissage. Le réseau de neurones utilisé dans le cas des systèmes monovariables nécessite uniquement la modification du nombre d’entrées et du nombre de neurones de la couche de sortie pour être exploité dans le cas des systèmes multivariables. L’algorithme génétique codé réel utilisé pour l’estimation des coefficients de la fonction d’activation polynomiale des neurones de la couche cachée reste pratiquement le même. Le seul changement dans cet algorithme concerne le critère de performance qui intervient dans la fonction d’évaluation (fitness). Pour le calcul de la loi de commande prédictive non linéaire sous contraintes, nous proposons l’extension des deux méthodes utilisées dans le chapitre 3 (NM et ROS). Ce chapitre est organisé comme suit : Dans le deuxième paragraphe, on présente la structure du modèle NARMA dans le cas multivariable. Le troisième paragraphe est réservé au développement de l’approche utilisée pour l’estimation structurelle et paramétrique du modèle. Le quatrième paragraphe traite la commande prédictive multivariable, en exploitant 102 Chapitre 4 les méthodes utilisées dans le chapitre 3, et le cinquième paragraphe présente les résultats de simulation. La conclusion est donnée dans le dernier paragraphe. 4.2. Le modèle NARMA (cas multivariable) On considère la classe des systèmes non linéaires qui peuvent être modélisés par la relation suivante : Y ( k ) = G (φ ( k )) + E ( k ) avec φ ( k ) = [Y ( k − 1) ... Y ( k − n ) U ( k − 1) ... U ( k − m ) ] (4-1) G est une fonction non linéaire supposée inconnue. Les paramètres m et n représentent, respectivement, l’ordre de la régression sur les entrées et l’ordre de la régression sur les sorties. E(k) est le vecteur de bruit. U (k ) ∈ ℜ ( m1 ,1) et Y (k ) ∈ ℜ ( n1 ,1) désignent, respectivement, le vecteur des entrées et le vecteur des sorties. Pour un système multivariable ayant m1 entrées et n1 sorties, chaque sortie peut être exprimée en fonction du vecteur d’observation par la relation suivante : y i (k ) = g i (φ (k )), i = 1, ..., n1 (4-2) avec gi est une fonction non linéaire et φ (k ) est formé par les anciennes mesures de toutes les sorties et de toutes les entrées : φ ( k ) = [u 1 ( k − 1) ... u 1 ( k − m ) ... u m 1 ( k − 1) ... u m 1 ( k − m ) y1 ( k − 1) ... y1 ( k − n ) ... y n1 ( k − 1) ... y n1 ( k − n ) ] (4-3) La structure dynamique de la sortie du modèle général de type NARMA, opérant dans un environnement déterministe, peut être régie par une équation mathématique discrète. Pour des raisons de simplification nous proposons une présentation matricielle, donnée par l’équation (4-4). Dans cette présentation nous avons intégré la valeur moyenne de chaque sortie du système dans la matrice des paramètres, puisque leurs valeurs seront identifiées avec l’ensemble des paramètres du modèle. 103 Chapitre 4 [ Yˆ ( k ) = yˆ 1 ( k ) yˆ 2 ( k ) ⎡ y1 ⎢b ⎢ 111 ⎢ M ⎢ ⎢ b11m ⎢ b121 ⎢ ⎢ M ⎢b 12 m =⎢ ⎢ M ⎢b ⎢ 1m1m ⎢ a111 ⎢ ⎢ M ⎢ a1n n ⎢ 1 ⎢ M ⎢ c1r ⎣ yi ... yˆ n1 ( k ) y2 b211 ... ... b21m b221 ... ... b22 m ... b2 m1m ... a 211 ... a 2 n1n ... c2 r ... ] (4-4) T 1 y n1 ⎤ ⎡ ⎤ ⎥ ⎢ ⎥ bn111 u1 (k − 1) ⎥ ⎥ ⎢ ⎥ M ⎥ ⎢ M ⎥ ⎥ ⎢ bn11m ⎥ ⎢ u1 (k − m) ⎥ ⎥ bn1 21 ⎥ ⎢ u 2 (k − 1) ⎥ ⎥ ⎢ M ⎥ ⎢ M ⎥ ⎥ ⎢ ⎥ bn1 2 m u 2 (k − m) ⎥ ⎥ ⋅⎢ M ⎥ ⎢ M ⎥ ⎥ ⎢ ⎥ bn1m1m u m1 (k − m) ⎥ ⎥ ⎢ ⎥ a n111 ⎥ ⎢ y1 (k − 1) ⎥ ⎥ ⎢ M ⎥ ⎢ M ⎥ ⎥ ⎢ ⎥ a n1n1n y n1 (k − n) ⎥ ⎥ ⎢ M ⎥ ⎢ M ⎥ ⎢ ⎥ cn1r ⎦ ⎣ y n1 (k − n) ... y n1 (k − n)⎥⎦ T représente la valeur moyenne de la sortie i. (i = 1, ..., n1 ) du modèle peut être écrite sous la forme suivante: Chaque sortie yˆ i (k ) ˆ T Φ(k ) yˆ i (k ) = Θ i (4-5) avec Θ̂ i T et Φ (k ) désignent, respectivement, le vecteur de paramètres et le vecteur d'observation. ˆ = [y Θ i i Φ (k ) = [1 où θ i1 θ i 2 ... θ ir ]T v1 v2 ... vr ]T (4-6) (4-7) v 1 ( k ) = u 1 ( k − 1), ... , v m ( k ) = u 1 ( k − m ), v m + 1 ( k ) = u 2 ( k − 1 ), ... , v 2 m ( k ) = u 2 ( k − m ), v 2 m + 1 ( k ) = u 3 ( k − 1 ), ... , v m 1 m ( k ) = u m 1 ( k − m ), v m 1 m + 1 ( k ) = y 1 ( k − 1), ... , v n1 n + m 1 m ( k ) = y n1 ( k − n ),..., (4-8) M v n1 n + m 1 m + 1 ( k ) = u 1 ( k − 1 ) 2 , ... , v r ( k ) = y n1 ( k − n ) q 104 Chapitre 4 Le nombre de termes possibles qui peuvent intervenir dans l’expression de chaque sortie du modèle yˆ i (k ) est donné par la relation suivante : r= (n1 n + m1 m + q )! −1 q! (n1 n + m1 m)! (4-9) où le paramètre q représente le degré de non linéarité du modèle. Par conséquent, le nombre de modèles possibles qu’on peut avoir pour un système qui comporte n1 sorties est : w = ( 2 r − 1) n1 (4-10) Pour m, n, q, n1 et m1 donnés, il existe un nombre très élevé de combinaisons possibles. Par exemple, si le système comporte deux entrées (m1=2) et deux sorties (n1=2), pour m=n=q=2, on aura 3.0949 1026 combinaisons différentes. Il est donc nécessaire d’utiliser une procédure permettant la sélection d’un modèle parmi l’ensemble des modèles possibles qui représente le comportement du système réel. On présente par la suite la méthode utilisée pour l’identification du modèle multivariable de type NARMA. 4.3. Identification du modèle NARMA Notre méthode exploite hors ligne un fichier formé par les mesures de l’entrée et de la sortie du système : (Y(k),U(k), k=1,…,N). En effet, les paramètres du modèle NARMA peuvent être estimés en utilisant un réseau de neurones à fonction d’activation polynomiale, comme décrit dans le premier chapitre. La structure du réseau de neurones est de type «feedforward» avec une seule couche cachée. La typologie de ce réseau dans le cas multivariable est représentée dans la figure 4-1. Le réseau neuronal reçoit en entrée les éléments du vecteur d’observation φ(k ) et délivre à sa sortie une estimation des sorties actuelles. Les sorties de la première couche sont données par : S (k ) = f (W 1Φ T (k )) = [s1 (k ) ... s M (k )]T (4-11) 105 Chapitre 4 W1 est la matrice des pondérations entre la couche d’entrée et la couche cachée. f étant une fonction d’activation polynomiale et M désigne le nombre de neurones de la couche cachée. La relation (4-12) présente une fonction d’activation polynomiale d’ordre q : f ( x i ) = f 0 + f 1 x i + f 2 x i2 + ... + f q x iq (4-12) En tenant compte des poids des connexions entre la couche d’entrée et la couche cachée, la quantité xi sera exprimée par : xi = ∑ + ... + ∑ m j =1 ∑w m1m 1 ji u m1 ( k j = m ( m1 −1) +1 m1 m + n j = m1 m +1 + ... + w1ji u1 ( k − j ) + w1ji y1 ( k ∑ ∑w 2m j = m +1 1 ji u 2 ( k − j + m) − j + m ( m1 − 1)) + − j + m1m ) + m1m + n1n w1ji y n1 ( k j = m1 m + n ( n1 −1) +1 ∑w m1 m + 2 n 1 ji y 2 ( k j = m1 m + n +1 − j + m1m + n ) − j + m1m + n ( n1 − 1)). i ∈ [1, M ] (4-13) ∑/ f y1(k-1) y1(k-n) ∑ yn1(k-1) yˆ1( k ) yn1(k-n) ∑ u1(k-1) yˆ n1(k ) u1(k-m) um1(k-1) um1(k-m) Couche d'entrée ∑/ f W1 W2 Couche cachée Couche de sortie Figure 4-1 : Réseau de neurones à une seule couche cachée (cas d’un système multivariable) 106 Chapitre 4 Comme dans le cas des systèmes monovariables, le choix de M dépend aussi de la complexité du système. En effet, un nombre faible de neurones cachés peut conduire à un modèle non précis alors que le choix d’un nombre élevé de neurones cachés peut rendre le modèle beaucoup plus compliqué, puisqu’il s’agit des combinaisons non linéaires de toutes les entrées et de toutes les sorties. Ainsi, le nombre de neurones cachés est choisi de façon à réaliser un compromis entre la précision et la complexité du modèle estimé. L’équation de la sortie du réseau sera alors: Yˆ (k ) = [ yˆ1 (k ) ... yˆ n1 (k )]T = W 2 S (k ) où yˆ j ( k ) = w ij2 f ( x i ). ∑ i =1 M (4-14) j = [1,.., n1 ] W2 désigne la matrice des pondérations entre la couche cachée et la couche de sortie. En remplaçant la quantité f ( xi ) de l’équation (4-12) par sa valeur dans (4-14), on ( obtient : ( yˆ j ( k ) = w12j f 0 + f 1 x1 + ... + f q x1q ( ) + w 22 j f 0 + f 1 x 2 + ... + f q x 2q ) ) q 2 . + ... + w Mj f 0 + f 1 x M + ... + f q x M (4-15) En remplaçant (xi) par son expression donnée par (4-13), l’expression de la sortie du réseau sera comme suit : yˆ j ( k ) = Θ TF Φ ( k ) où ( Θ F = h f 0 , f 1 ..., f q , W 1 , W 2 ), (4-16) h: fonction non linéaire, dim Θ F = [r + 1, 1]. Les pondérations (W1 et W2) des connexions du réseau de neurones sont déterminées avec l’algorithme de rétro propagation du gradient. La déduction des paramètres du modèle NARMA se fait alors à partir des pondérations (W1 et W2) et des paramètres du polynôme considéré (f0, f1,…,fq). En effet, par identification de l’équation (4-16) avec l’équation (4-4), on trouve les coefficients des termes du modèle NARMA. 107 Chapitre 4 Puisque la qualité du modèle obtenu, avec cette méthode, dépend du choix adéquat des coefficients (f0, f1,…,fq) du polynôme. Nous proposons l’incorporation de la même procédure utilisée dans le chapitre 2, c’est à dire la combinaison de l’algorithme génétique réel avec l’algorithme de rétro propagation du gradient pour déterminer ces coefficients. La figure 4-2 présente la structure du superviseur utilisé pour la détermination du modèle NARMA. U(k) Y(k) Procédé - Modèle RNA + Ŷ(k) Algorithme d’adaptation W1, W2 Algorithme Génétique f0, f1,…,fq Figure 4-2 : Structure du superviseur basé sur l'algorithme génétique Le critère J3 qui intervient dans la fonction d’évaluation (fitness) de l’algorithme génétique codé réel s’exprime comme suit : J3 = ∑ (( y ( k ) − yˆ ( k ) ) N k =1 2 1 1 + ... + ( y n1 ( k ) − yˆ n1( k ) ) 2 ) (4-17) Avec cette méthode, le calcul des paramètres du modèle NARMA se fait de la même façon que dans le cas monovariable, il faut tenir compte seulement du nombre d’entrées et de sorties. Par ailleurs, le rôle de l’algorithme génétique est de trouver les valeurs des coefficients de la fonction d’activation des neurones alors que l’algorithme de propagation du gradient s’occupe de l’ajustement des poids des connexions. 108 Chapitre 4 4.4. Commande prédictive multivariable Le principe de la commande prédictive présenté dans la cas monovariable reste le même dans le cas multivariable. Cependant, pour traiter les systèmes multivariables, on doit tenir compte de certaines spécificités tel que : le modèle utilisé, le calcul du prédicteur, le critère de performance et la méthode d’optimisation adoptée. - 4.4.1. Calcul du prédicteur Le calcul du prédicteur multivariable yˆ (k + j ) j = 1, ... , Hp , nécessite l’exploitation de l’équation du modèle. A partir de l’équation (4-4), qui décrit le modèle général de type NARMA, l’équation du prédicteur sera donnée comme suit : [ Yˆ ( k + l ) = yˆ 1 ( k + l ) yˆ 2 ( k + l ) y2 ⎡ y1 ⎢b ⎢ 111 b211 ⎢ M ⎢ ⎢ b11m b21m ⎢ b121 b221 ⎢ ⎢ M ⎢b b22m 12m =⎢ ⎢ M ⎢b b ⎢ 1m1m 2m1m ⎢ a111 a211 ⎢ ⎢ M ⎢ a1n n a2n n 1 ⎢ 1 M ⎢ ⎢c c2r ⎣ 1r 4.4.2. ... ... ... ... ... ... ... ... ... ] T yˆ n1 ( k + l ) ; l = 1, ... , Hp ... (4-18) 1 yn1 ⎤ ⎡ ⎤ ⎥ ⎢ ⎥ bn111 ⎥ ⎢ u1 (k − 1 + l ) ⎥ ⎥ M M ⎥ ⎢ ⎥ ⎢ ⎥ u1 (k − m + l ) bn11m ⎥ ⎢ ⎥ ⎥ ⎢ ⎥ bn1 21 u 2 (k − 1 + l ) ⎥ ⎢ ⎥ M M ⎥ ⎢ ⎥ ⎥ ⎢ ⎥ bn1 2m u 2 (k − m + l ) ⎥ ⋅⎢ ⎥ M M ⎥ ⎢ ⎥ ⎥ ⎢ ⎥ bn1m1m um1 (k − m + l ) ⎥ ⎢ ⎥ ⎥ an111 ⎥ ⎢ y1 (k − 1 + l ) ⎥ ⎢ ⎥ M M ⎥ ⎢ ⎥ ⎥ an1n1n ⎥ ⎢ yn1 (k − n + l ) ⎥ ⎢ ⎥ M M ⎥ ⎢ ⎥ ⎥ ⎢ cn1r ⎦ ⎣ y n1 (k − n + l ) ... y n1 (k − n + l )⎥⎦ T Le critère de performance Le critère de performance utilisé contient (n1 + m1 ) termes quadratiques. Les (n1 ) premiers termes sont en relation avec les erreurs futures entre les sorties prédites et les 109 Chapitre 4 consignes désirées, tandis que, les m1 autres termes correspondent aux incréments futurs de commandes liées aux m1 entrées du système. Le critère utilisé est ainsi donné par la relation suivante : 1 ⎡ n1 j J 4 (k ) = ⎢∑∑ (yˆ j (k + i) -y j c(k + i))2 + 2 ⎣ j =1 i =1 Hp ⎤ 2 λ ( Δu (k + i1 ) ) ⎥ ∑∑ j j j =1 i =1 ⎦ m1 Hc j (4-19) avec : m1 : nombre des signaux de commande du système, n1 : nombre de sorties du système, Hpj: horizon de prédiction sur la sortie j, Hcj : horizon de commande sur le signal de commande j, λj: coefficient de pondération sur la commande (pour l’entrée j). Dans le cas où les valeurs de l’horizon de prédiction pour toutes les sorties sont identiques, les valeurs de l’horizon de commande pour toutes les entrées sont égales et les coefficients de pondération sur la commande sont identiques, l’équation du critère sera alors comme suit : m1 Hc ⎤ 1 ⎡ n1 Hp 2 J (k ) = ⎢∑∑ (yˆ j (k + i) -y j c(k + i)) + λ∑∑ ( Δu j (k + i-1 )) 2 ⎥ 2 ⎣ j =1 i =1 j =1 i =1 ⎦ 4.4.3. (4-20) Optimisation du critère Pour le calcul de la loi de commande, nous proposons l’extension des méthodes d’optimisation présentées dans le chapitre 3, à savoir la méthode du simplexe et celle de Rosenbrock, pour le cas multivariable. Afin de traiter les contraintes et pour avoir des bonnes performances, nous avons intégré les mêmes améliorations introduites dans le cas monovariable, c’est-à-dire, la méthode de pénalité, l’approche CFON, la notion de multi initialisation ainsi que la notion d’élitisme. Au niveau de l’initialisation des algorithmes d’optimisations on considère une matrice de dimension (m, n+1) dans le cas de la méthode de Nelder-Mead et m directions de recherche 110 Chapitre 4 pour la méthode de Rosenbrock. m désigne le nombre de variables à optimiser et n+1 représente la taille de chaque simplexe. La fonction de pénalité est définie de manière à respecter, pour chaque entrée du système, les contraintes sur le signal de commande et sur ses incréments. Les contraintes peuvent être écrites comme suit : u imin ≤ u i ( k + j ) ≤ u imax , Δ u imin ≤ Δ u i ( k + j ) ≤ Δ u imax , j = 0 , ... , Hc − 1 , i = 1, ... , m (4-21) j = 0 , ... , Hc − 1 , i = 1, ... , m (4-22) Hc : Horizon de commande. En partant l’équation (3-17) du chapitre 3, nous pouvons écrire la fonction de pénalité dans le cas multivariable comme suit : α (x) = ∑ ∑ [Max {0 , m i =1 4 Hc j =1 ] f ij ( x ) } 2 (4-23) L’application de l’approche CFON et de la notion multi initialisation, en utilisant les versions d’optimisation NM-V3 et ROS-V3, nécessite l’intégration d’une procédure qui comporte deux étapes : Etape 1 : Fixer les intervalles espaces de recherche pour chaque variable (Sj, j=1,…,m) en respectant globalement toutes les contraintes. Etape 2 : Diviser chaque espace de recherche Sj en n sous intervalles égaux si (i = 1,…, n) et refaire à chaque fois l’initialisation des différentes variables de manière à balayer toutes les combinaisons possibles des sous intervalles si. Le choix du nombre de divisions de l’espace de recherche n dépend de la complexité du problème. En effet, pour un problème d’optimisation simple on prend n faible et pour un cas complexe on choisi un nombre plus important. Dans le chapitre 3 nous avons présenté l’algorithme de Rosenbrock dans le cas général (mono et multivariable). 111 Chapitre 4 L’algorithme de Nelder-Mead (cas multivariable) L’algorithme de Nelder-Mead, appliqué aux problèmes d’optimisation multivariables, est résumé par les étapes suivantes : Etape 1 : Prendre m polyèdres de départ (m : nombre de variables) Prendre les sommets ( x1 , x 2 ,..., x n +1 ) de chaque polyèdre vérifiant f(x 1 ) ≤ f(x 2 ) ≤ … ≤ f(x n +1 ). Etape 2 : Tant que la différence entre f(x n +1 ) et f(x 1 ) est supérieure à un certain seuil de précision ε 0 . 2.1 calculer x et f r = f(x( ρ )) avec x= pour chaque polyèdre 1 n ∑ xi , est le centroïde des n meilleurs sommets, de n i =1 chaque polyèdre, et 2.2 x( ρ ) = x + ρ ⋅ ( x − x n +1 ) réflexion : si f(x 1 ) ≤ f r < f(x n ). Alors x n +1 = x( ρ ) pour chaque polyèdre, aller à 2.7. 2.3 expansion : si f r < f(x 1 ) alors calculer f e = f(x( χ )) , avec x( χ ) = x + χ ⋅ ( x r − x ) si f e < f r alors, pour chaque polyèdre prendre x n +1 = x( χ ) et aller à 2.7 sinon x n +1 = x( ρ ) et aller à 2.7. 2.4 contraction externe : si f(x n ) ≤ f r < f(x n +1 ) alors calculer f o = f(x(γ )) ,avec x(γ ) = x + γ ⋅ ( x r − x ) si f o ≤ f r alors pour chaque polyèdre prendre x n +1 = x(γ ) puis aller à 2.7 sinon aller à 2.6 2.5 contraction interne si f r ≥ f(x n +1 ) alors calculer f i = f(x(γ )) , avec x(γ ) = x − γ ⋅ ( x − x n +1 ) si f i < f(x n +1 ) , pour chaque polyèdre prendre x n +1 = x(γ ) puis aller à 2.7 sinon aller à 2.6. 112 Chapitre 4 2.6 rétrécissement pour i = 2, ..., n + 1 , calculer pour chaque polyèdre prendre xi = (xi + x1 ) 2 , puis calculer f ( xi ) 2.7 donner les sommets ( x1 , x 2 ,..., x n +1 ) de chaque polyèdre prendre tel que f(x 1 ) ≤ f(x 2 ) ≤ … ≤ f(x n +1 ). retour à l’étape 2. 4.5. Résultats de simulation A. Identification des systèmes multivariables Au cours de la phase d’identification, le signal d’entrée utilisé est un bruit blanc de moyenne nulle et de variance égale à un. Les coefficients de corrélation totaux R tot2 (i=1, …, n1) i sont calculés afin de décider si le modèle obtenu est acceptable ou non. En effet, le modèle est acceptable si ces critères sont proches de 1 [Haber, 1990]. ∑ ( y ( k ) − yˆ ( k )) N 2 R tot =1− i k =1 i ∑ ( y ( k )) N k =1 2 i , i=1, …, n1 (4-24) 2 i Les paramètres de réglage utilisés dans les simulations sont résumés dans le tableau 4-1. Algorithme génétique réel Réseau de neurones artificiels Nombre de générations 40 - Nombre d'individus 100 - Probabilité de croisement 0.8 - Probabilité de mutation 0.1 - Nombre d'époques - 5000 Coef. d’apprentissage - 0.005 Tableau 4-1 : Paramètres de réglage des algorithmes Trois systèmes multivariables non linéaires, avec deux entrées et deux sorties, sont considérés. Le premier est présenté par une structure de type NARMA. Le deuxième exemple est un benchmark qui présente un système plus compliqué. Cependant, le troisième exemple 113 Chapitre 4 concerne une application pratique, il s’agit d’un procédé de régulation de niveau d’eau caractérisé par trois réservoirs interconnectés. Le choix (m1=n1= m=n=q=2), permet d’avoir 44 termes possibles dans l’expression générale de chacune des sorties du modèle NARMA. Système 1 : Le premier système est décrit par une équation de type NARMA donnée comme suit : ⎧ y 1 ( k ) = 0 . 8 y 1 ( k − 1) + u 1 ( k − 2 ) + 0 . 4 ( u 1 ( k − 2 )) 2 ⎪ − 1 . 2 u 1 ( k − 1) u 2 ( k − 2 ) − 0 . 1 y 2 ( k − 1). ⎪ ⎨ 2 ⎪ y 2 ( k ) = 0 . 5 y 2 ( k − 1) + u 2 ( k − 2 ) + ( u 1 ( k − 1)) ⎪ + 0 . 5 y 2 ( k − 2 ) u 2 ( k − 1). ⎩ (4-25) Système 2 : Le deuxième système multivariable est caractérisé par l’équation 4- 26. Il présente un cas plus complexe [Song, 2006], [Petlenkov, 2007]. a1 y1 (k − 1) y1 (k − 2) ⎧ y k + a4 u1 (k − 2) + a5 u1 (k − 1) + a6 u2 (k − 2) = ( ) 1 ⎪ 1 + a2 y12 (k − 1) + a3 y22 (k − 2) ⎪⎪ ⎨ (4-26) ⎪ b y (k − 1) sin( y2 (k − 2)) ⎪ y 2 (k ) = 1 2 2 + b4 u 2 (k − 2) + b5 u2 (k − 1) + b6 u1 (k − 2) 1 + b2 y2 (k − 1) + b3 y12 (k − 2) ⎪⎩ Les paramètres du système sont : a1=0.7 ; a2=1 ; a3=1 ; a4=0.3 ; a5=1 et a6=0.2 ; b1=0.5 ; b2=1 ; b3=1 ; b4=0.5 ; b5=1 et b6=0.2 . Système 3 : Le troisième système multivariable est schématisé par la figure 4-3(b). Il présente un procédé de régulation de niveau d’eau. Ce procédé de laboratoire est constitué de trois réservoirs cylindriques connectés entre eux par des valves et qui possèdent une connexion avec un bassin de niveau plus bas (figure 4-3(a)). Le niveau d’eau dans chaque réservoir 114 Chapitre 4 dépend de la quantité d’eau qui entre et qui sort du réservoir. Cette quantité dépend du débit donné par les vannes V1 et V2. Le paramètre hi (i=1,2,3) désigne le niveau de l’eau dans le réservoir i, q1 et q2 sont les débit dans les réservoirs 1 et 2, Sj représente la section de la colonne j, qij représente le débit de l’eau qui passe du réservoir i au réservoir j et qi0 représente le débit sortant du réservoir i. Dans notre cas, les sorties du modèle NARMA sont les niveaux dans les réservoirs 1 et 2 (y1=h1, y2=h2) et les entrées sont les débits d’entrées (u1=q1, u2=q2). Un fichier de 1600 mesures (entrées/sorties) est considéré pour la modélisation du système, les 1100 premières mesures sont utilisées pour l’identification du modèle et les 500 dernières mesures sont utilisées pour la validation de ce dernier. Figure 4-3(a) : Procédé de régulation du niveau d’eau à trois réservoirs interconnectés 115 Chapitre 4 Le procédé de la figure 4-3(a) peut être schématisé par la figure 4-3(b), comme suit : V1 V2 q2 q1 R3 S3, h3 S1, h1 R1 q13 q30 q10 S2, h2 R2 q23 q20 Figure 4-3(b) : Schéma du système à trois réservoirs (système multivariable 3) Les coefficients de la fonction d’activation trouvés par l’algorithme génétique réel, pour les trois systèmes sont donnés dans le tableau 4-2. Coefficients de la fonction d’activation polynomiale f0 f1 f2 Système 1 15 -4.16 13.5 Système 2 2.5 14.2 25 Système 3 5 26.3 35 Tableau 4-2 : Coefficients de la fonction d’activation polynomiale pour les trois systèmes identifiés Le tableau 4-3 présente les termes sélectionnés ainsi que les paramètres estimés des modèles NARMA correspondants aux trois systèmes identifiés. 116 Chapitre 4 N Termes Possibles 0 y 1 u1(k-1) u1(k-2) u2(k-1) u2(k-2) y1(k-1) y1(k-2) y2(k-1) y2(k-2) u1²(k-1) u1(k-1) u1(k-2) u1(k-1) u2(k-1) u1(k-1) u2(k-2) u1²(k-2) u1(k-2) u2(k-1) u1(k-2) u2(k-2) u2²(k-1) u2(k-1) u2(k-2) u2²(k-2) u1(k-1) y1(k-1) u1(k-1) y1(k-2) u1(k-1) y2(k-1) u1(k-1) y2(k-2) u1(k-2) y1(k-1) u1(k-2) y1(k-2) u1(k-2) y2(k-1) u1(k-2) y2(k-2) u2(k-1) y1(k-1) u2(k-1) y1(k-2) u2(k-1) y2(k-1) u2(k-1) y2(k-2) u2(k-2) y1(k-1) u2(k-2) y1(k-2) u2(k-2) y2(k-1) u2(k-2) y2(k-2) y1²(k-1) y1(k-1) y1(k-2) y1(k-1) y2(k-1) y1(k-1) y2(k-2) y1²(k-2) y1(k-2) y2(k-1) y1(k-2) y2(k-2) y2²(k-1) y2(k-1) y2(k-2) y2²(k-2) 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 Système 1 Système 2 Système 3 Système 3 (rectifié) y1(k) y2(k) y1(k) y2(k) y1(k) y2(k) y1(k) y2(k) 0 0 0.0465 -0.0142 0.0113 0.0005 0.0113 0 0.0037 -0.0028 0.0037 -0.0028 0 0 0.9918 0.0001 1.0001 0 0.4076 0.1826 0.0027 0.0034 0.0027 0 0 0.0012 0.9988 -0.0023 0.0032 -0.0023 0.0032 0.0001 1 0.1880 0.5282 0.0040 -0.0074 0.0040 -0.0074 0.7999 0 -0.0749 0.0104 0.4862 0.1909 0.4862 0.1909 -0.0553 0.3090 -0.0553 0.0034 0.0001 0 0.0610 0.0015 0.3090 -0.100 0.5 0.0174 -0.0318 0.0342 0.4367 0.0342 0.4367 0.0500 0.4573 0 0 0.0024 0.0040 0.0500 0.4573 0 1 0.0148 -0.0120 0.0000 0.0000 0 0 0 0 0.0297 -0.0240 0.0000 0.0001 0 0 0 0 0.0028 -0.0017 -0.0003 0.0006 0 0 -1.200 0 0.0103 0.0104 -0.0000 0.0010 0 0 0.0000 0.0000 0 0 0.3998 0 0.0043 -0.0128 0 0 -0.0258 -0.0235 0.0002 -0.0002 0 0 -0.0000 -0.0002 0 0 0 0 -0.0062 -0.0179 0 0 -0.0178 0.0065 0.0012 -0.0005 0 0 0 -0.0035 0.0246 -0.0006 0.0033 0 0.0033 0.0090 0 0 0 0.0090 0.0409 -0.0006 0.0090 0 0 0 -0.0533 0.0357 -0.0020 -0.0013 0 0 0 0 0.0525 -0.0186 -0.0010 -0.0002 0 0 -0.0157 -0.0005 -0.0019 0 0 0.0224 -0.0004 -0.0018 0 0.0006 0 0 0 0 0.0126 0.0193 0 0 0 -0.0181 0.0192 0.0028 0 0 0.0196 -0.0079 0.0018 -0.0002 0 0 0 0 0.0026 -0.0245 0.0005 0.0018 0 0 0 0 0.0087 -0.0058 0.0002 0.0019 0 0 0 0 0.0291 0.0257 0.0016 -0.0012 0 0 0 0 -0.0148 0.0089 0.0006 0.0005 0 0 0 0 -0.0112 -0.0489 0.0017 -0.0039 0 0.0240 -0.0006 -0.0043 0 -0.0043 0.0212 0.0003 -0.0044 0 -0.0044 0 0.5 0.0213 0.0028 0 -0.0039 0 0 0.0083 0 0 0.0006 0.0031 0.0025 -0.0045 0.0025 0 0 -0.0065 -0.0147 -0.0053 0.0019 -0.0053 0 0 0.0149 0.0455 -0.0048 -0.0043 -0.0048 -0.0043 0 0 0.0042 -0.0036 0.0732 -0.0336 0.0732 -0.0336 0 0 -0.0042 0.0032 0.0476 0.0019 0.0476 0 -0.0389 0.0488 -0.0389 -0.0045 0 0 0 0.0050 0.0017 0.0488 0 0 0.0070 0.0024 0.0175 -0.0516 0.0175 -0.0516 0.0180 0.0215 0.0180 0 0 0.0042 -0.0017 0.0215 0 0 -0.0029 -0.0001 0.0529 -0.0440 0.0529 -0.0440 0 0 -0.0035 0.0036 0.0336 -0.0410 0.0336 -0.0410 0.0002 0 -0.0100 0.0190 -0.0438 0.0829 -0.0438 0.0829 0 0 -0.0088 0.0305 -0.0466 0.0588 -0.0466 0.0588 0.0001 0 -0.0028 0.0611 -0.0489 0.0345 -0.0489 0.0345 Tableau 4-3 : Termes et paramètres estimés pour les trois systèmes multivariables 117 Chapitre 4 Dans le tableau 4-3, la colonne notée (système 3 (rectifié)), comporte seulement les termes dont les paramètres sont inférieurs ou égales à 2 10-3. Le tableau 4-4 présente la somme des erreurs quadratiques normalisées des deux sorties du modèle du système 3 et du modèle rectifié correspondant. Erreurs quadratiques normalisées (modèle et modèle rectifié) du système 3 Erreur par rapport à la sortie y1 -5 1.5446 10 Erreur par rapport à la sortie y2 -5 1.3019 10 Tableau 4-4 : Valeurs de la somme des erreurs quadratiques normalisées entre le modèle trouvé et le modèle rectifié (système 3) Le tableau 4-5 comporte les valeurs des coefficients de corrélation totaux (R²tot ) et du temps de calcul des simulations effectuées avec un PC Pentium 4, 2 GHz. Système 1 R²tot Temps (sec) Système 2 y1(k) y2(k) y1(k) 0.9953 0.9998 0.9926 2295.120 y2(k) 0.9911 2331.56 Système 3 y1(k) y2(k) 0.9936 0.9931 2531.56 Tableau 4-5 : Les valeurs de R²tot et du temps de calcul (système multivariable) Le nombre de neurones de la couche caché est : M = 5. Commentaires : On remarque, d’après le tableau 4-3, que les termes et les paramètres du système 1 sont identifiés avec succès. Le bon choix des coefficients de la fonction d’activation a permis d’avoir une bonne précision. On constate, d’après le tableau 4-5 que les valeurs de R²tot sont très proche de 1 ce qui prouve que les modèles trouvés sont précis. Les performances des méthodes d’identification peuvent être améliorées en augmentant le nombre d’itérations, puisque l’identification se fait hors ligne. Afin de valider le modèle correspondant au système 3, on a comparé les évolutions des sorties du système avec celles du modèle identifié. Les séquences d’entrées sont choisies de type échelon de différents niveaux. La figure 4-4 présente respectivement les séquences d’entrées (u1, u2), les sorties système/modèles et l’erreur de modélisation pour la première 118 Chapitre 4 sortie y1 et les sorties système/modèles et l’erreur de modélisation pour la deuxième sortie y2. On remarque que les erreurs de modélisation sont très faibles, ce qui permet de conclure que le modèle obtenu caractérise convenablement la dynamique du système considéré. entrées système (u1, u2) 3 2 u1 1 u2 0 0 200 400 600 800 1000 1200 sorties 1 (système, modèle et modèle rectifié) et erreur de modélisation 1400 1600 0.6 0.4 y1 0.2 y1m erreur1 y1m(rectifié) 0 -0.2 0 200 400 600 800 1000 1200 1400 1600 sorties 2 (système, modèle et modèle rectifié) et erreur de modélisation 0.6 0.4 y2 y2m erreur2 0.2 y2m(rectifié) 0 -0.2 0 200 400 600 800 itérations k 1000 1200 1400 1600 Figure 4-4 : Evolutions des séquences d’entrées (u1, u2), des sorties (système et modèle) et des erreurs de modélisation Les sorties du système, du modèle et du modèle rectifié sont pratiquement confondues. Le tableau 4-4 montre aussi que l’erreur entre le modèle et le modèle rectifié est très faible. Ceci nous permet de considérer le modèle rectifié, qui présente moins de termes avec une précision comparable au modèle obtenu initialement. 119 Chapitre 4 B. Commande des systèmes multivariables Pour appliquer nos algorithmes de commande, nous avons considéré le système multivariable non linéaire de l’équation (4-26). Dans une première étape, nous avons introduit des perturbations (d1 et d2) respectivement aux deux sorties (y1 et y2) du système considéré. Ensuite, nous avons supposé des variations au niveau des paramètres du système de l'ordre de ± 10 % et ± 20 %. La loi de commande a été calculée en utilisant le modèle nominal. Ainsi, le système objet de l’étude, avec incertitudes sur les paramètres et en présence de perturbations, est donné par l’équation 4-27. a1 y1 (k −1) y1 (k − 2) ⎧ ⎪y1 (k) = 1+ a y2 (k −1) + a y2 (k − 2) + a4 u1 (k − 2) + a5 u1 (k −1) + a6 u2 (k − 2) + d1 (k) 2 1 3 2 ⎪⎪ ⎨ ⎪ b y (k −1) sin(y2 (k − 2)) ⎪y2 (k) = 1 2 2 + b4 u2 (k − 2) + b5 u2 (k −1) + b6 u1 (k − 2) + d2 (k) 1+ b2 y2 (k −1) + b3 y12 (k − 2) ⎪⎩ (4-27) Les caractéristiques des signaux de perturbation d1(k) et d2(k), ajoutés respectivement aux sorties du système y1(k) et y2(k), sont définies par l’équation (4-28). d1 (k ) = 0.2 d 2 (k ) = 0.2 k ∈ [120, ... ,170] ; 0 ailleurs k ∈ [80, ... ,120] ; 0 ailleurs (4-28) Dans les algorithmes de commande, nous avons considéré les versions d'optimisation (NM-V3 et ROS-V3) associés avec l’approche CFON ainsi que la fonction de pénalité. Les paramètres du correcteur prédictif considéré dans la simulation sont, Hc = 1, Hp = 4 et λ = 0.2. Les contraintes sur les entrées et leurs incréments sont comme suit : − 0 .1 ≤ u1 ( k ) ≤ 1; − 0 .1 ≤ u 2 ( k ) ≤ 1 − 0 .1 ≤ Δ u1 ( k ) ≤ 0 .1; − 0 .1 ≤ Δ u 2 ( k ) ≤ 0 .1 (4-29) Le tableau 4-6 présente les paramètres ai et bi (i=1,…,6) du système nominal et du système incertain. 120 Chapitre 4 Paramètres du système a1 a2 a3 a4 a5 a6 b1 b2 b3 b4 b5 b6 Système nominal 0.7 1 1 0.3 1 0.2 0.5 1 1 0.5 1 0.2 Système incertain ( ± 10 %) 0.77 0.9 1.1 0.33 0.9 0.22 0.55 1.1 0.9 0.55 1.1 0.18 Système incertain ( ± 20 %) 0.56 1.2 0.8 0.36 1.2 0.16 0.4 1.2 1.2 0.4 0.8 0.24 Tableau 4-6 : Paramètres des systèmes nominal et incertains Dans ces simulations, on a considéré les paramètres du système nominal pour les itérations 1 jusqu'à 200. Entre les itérations 201 et 400, nous avons utilisé les mêmes séquences des signaux des consignes et des perturbations pour le système avec incertitude de (+/- 10 %) sur les paramètres et entre les itérations 401 et 600 une incertitude de (+/- 20 %) sur les paramètres. Les figures 4-5 et 4-6 présentent les résultats trouvés avec la méthode NM et les figures 4-7 et 4-8 présentent les résultats trouvés avec la méthode ROS. 2 s ortie y1 cons igne yc1 1.5 perturbation d1 1 0.5 0 0 100 200 300 400 500 600 400 500 600 400 500 600 commande u1 1 0.5 0 0 100 200 300 Incréments de commande 0.1 0.05 0 -0.05 -0.1 0 100 200 300 itérations k Figure 4-5 : Sortie y1, commande et incréments de commande (NM-V3) 121 Chapitre 4 s ortie y2 cons igne yc2 1.5 perturbation d2 1 0.5 0 0 100 200 300 400 500 600 400 500 600 400 500 600 commande u2 1 0.5 0 0 100 200 300 Incréments de commande 0.1 0.05 0 -0.05 -0.1 0 100 200 300 itérations k Figure 4-6 : Sortie y2, commande et incréments de commande (NM-V3) 2 sortie y1 consigne yc1 perturbation d1 1.5 1 0.5 0 0 100 200 300 400 500 600 400 500 600 400 500 600 commande u1 1 0.5 0 0 100 200 300 incréments de commande 0.1 0.05 0 -0.05 -0.1 0 100 200 300 itérations k Figure 4-7 : Sorties y1, commande et incréments de commande (ROS-V3) 122 Chapitre 4 2 sortie y2 1.5 consigne yc2 perturbation d2 1 0.5 0 0 100 200 300 400 500 600 400 500 600 400 500 600 commande u2 1 0.5 0 0 100 200 300 incréments de commande 0.1 0.05 0 -0.05 -0.1 0 100 200 300 itérations k Figure 4-8 : Sortie y2, commande et incréments de commande (ROS-V3) Le tableau 4-7 présente les valeurs des critères de performances I1, I2 et I3, donnés par les équations (3-24 - 3-26) ainsi que le temps moyen par itération, obtenues par les deux méthodes d'optimisation avec le système nominal. Méthode Nelder-Mead Sortie y1 Sortie y2 Rosenbrock Sortie y1 Sortie y2 I1 0.0041 0.0049 0.0042 0.0049 I2 0.1803 0.2285 0.1802 0.2286 I3 7.0547 10-4 5.4768 10-4 7.0238 10-4 5.4668 10-4 tmoy (sec) 1.3128 0.2767 Tableau 4-7 : Valeurs des critères de performances et temps de calcul moyen (système multivariable) On remarque, d'après le tableau 4-7, que les performances des deux méthodes d'optimisation sont très proches, assurant un signal de commande de faible variance et une erreur de poursuite quadratique minimale. Cependant, la méthode de Rosenbrock est pratiquement cinq fois plus rapide que celle de Nelder-Mead. 123 Chapitre 4 4.6. Conclusion Le travail présenté, dans ce chapitre, concerne la modélisation et la commande prédictive des systèmes non linéaires multivariables en exploitant les modèles NARMA. Nous avons proposé l’extension de la méthode qui combine le réseau de neurones artificiels à fonction d’activation polynomiale et l’algorithme génétique réel, présentée dans le chapitre 2, pour l’identification des modèles multivariables de type NARMA. Cette méthode permet la sélection des termes ainsi que l’estimation des paramètres du modèle NARMA. Le modèle trouvé sera par la suite exploité pour la synthèse de la commande prédictive. Trois exemples ont été présentés pour valider notre méthode d’identification. Le premier exemple est un cas simple, le deuxième représente un cas plus compliqué et le troisième est une application pratique qui présente un procédé physique de laboratoire. Les résultats trouvés ont confirmé l’efficacité de la méthode proposée. En effet, les modèles NARMA déterminés caractérisent avec une précision acceptable et avec une complexité raisonnable le comportement des systèmes étudiés. Dans une deuxième étape, nous avons appliqué les algorithmes de commande (NM-V3 et ROS-V3) pour commander un système qui possède deux entrées et deux sorties. Les performances trouvées montrent la robustesse des algorithmes de commande proposés. En effet, les signaux de commande respectent d'une part les contraintes imposées et permettent d'autre part aux sorties du système d'atteindre les consignes désirées. De plus, le correcteur prédictif permet d'éliminer les perturbations et traite convenablement les incertitudes sur les paramètres du système. 124 Conclusion générale Conclusion générale Dans ce mémoire, nous nous sommes intéressé à l’identification et la commande prédictive des systèmes non linéaires mono et multivariables, en exploitant les modèles NARMA. Pour l’identification des modèles de type NARMA, nous avons présenté en premier lieu l’algorithme de régression par pas échelonné modifié. Dans l'implémentation de cet algorithme nous avons généralisé son fonctionnement de manière à considérer des ordres élevés, à savoir les ordres de régressions sur l'entrée et la sortie ainsi que le degré de non linéarité du modèle. Parmi les difficultés rencontrées dans la programmation de cet algorithme, on distingue le calcul des coefficients de corrélation partiels, des tests statistiques partiels et des critères d'informations partiels, puisqu'il s'agit de la permutation des termes afin de tester l'influence de chacun de ceux qui forment la structure du modèle courant. Ensuite, nous avons proposé deux autres méthodes heuristiques. La première méthode est basée sur les algorithmes génétiques binaires et la deuxième méthode constitue une combinaison entre le réseau de neurones artificiels à fonction d’activation polynomiale et l’algorithme génétique sous sa représentation réelle. Cette dernière permet de chercher les coefficients de la fonction polynomiale simultanément avec l’apprentissage du réseau de neurones. Les méthodes que nous avons proposées sont moins compliquées que la méthode de Kortmann et offrent des résultats similaires. D’autre part ces méthodes, contrairement à la méthode de Kortmann, sont capables de traiter les systèmes multivariables. Cet avantage nous a permis de faire l’extension de la deuxième méthode proposée pour la modélisation des systèmes multivariables. Les méthodes proposées sont également efficaces quand nous n'avons aucune connaissance a priori de la structure du système. Cet avantage permet à ces approches de traiter les systèmes complexes. Pour valider ces méthodes, nous avons considéré quatre exemples monovariables et trois multivariables, dont une application pratique. Les résultats trouvés, pratiques ou de simulations, ont confirmé l’efficacité et la robustesse des méthodes proposées. En effet, les modèles NARMA déterminés caractérisent avec une précision acceptable et avec une complexité raisonnable le comportement des systèmes étudiés. Par la suite nous avons proposé un contrôleur prédictif des systèmes non-linéaires sous contraintes, qui exploite les modèles de type NARMA. La loi de commande est obtenue en minimisant un critère quadratique non convexe. Le problème d'optimisation est résolu par deux méthodes utilisant l'algorithme de Nelder-Mead et l’algorithme de Rosenbrock qui 125 Conclusion générale n’utilisent pas la dérivée. Ces méthodes, qui sont à l'origine déterministes, locales et sans contraintes, combinées avec la fonction de pénalité, l’approche CFON ainsi que l’utilisation de la notion de multi initialisation permettent une meilleure convergence vers le minimum global. Les méthodes d’optimisation proposées présentent deux avantages essentiels, d’une part elles ne nécessitent pas le calcul de la dérivée et d’autre part elles sont facilement exploitables dans le cas des systèmes multivariables. Dans ce cadre nous avons appliqué les algorithmes de commande avec les versions (NM-V3 et ROS-V3) pour commander un système qui possède deux entrées et deux sorties. Les performances trouvées montrent la robustesse des algorithmes de commande proposés. En effet, la commande calculée respecte les contraintes imposées et permet de ramener les sorties du système aux consignes désirées même en présence de perturbation ou voir même dans le cas d’un système incertain. D’après les résultats de simulations, on constate que la méthode ROS requière peu de temps de calcul, par rapport à celle de NM et offre des bonnes performances. Ceci nous permet de recommander cette approche pour traiter les systèmes rapides. Dans la méthode d’identification basée sur la combinaison du réseau de neurones artificiels à fonction d’activation polynomiale et l’algorithme génétique sous sa représentation réelle, l’expression de la sortie du modèle NARMA peut être calculée par un programme informatique approprié qui prend en considération l'ensemble des paramètres du réseau de neurones. Par conséquent, on peut généraliser le fonctionnement de cet algorithme pour des degrés de non linéarités plus importants. Suite à ces travaux, plusieurs perspectives peuvent êtres proposées. Parmi lesquelles on peut citer: - Simplification, d'avantage, la complexité du modèle NARMA, en éliminant les termes ayant une faible contribution dans la sortie du modèle. - Optimisation de la structure et des paramètres du réseau de neurones permettant l'identification du modèle NARMA. - Amélioration des performances de l'algorithme NM par optimisation de ses paramètres, surtout pour la définition de la taille du simplexe. - Mise en œuvre pratique des approches proposées sur des procédés réels. 126 Bibliographie [Adetona, 2004] O. Adetona, S. Sathananthan, and L. H. Keel, 2004. Approximations of the NARMA model of Non-affine Plants, Proceeding of the 2004 American Control Conference, Boston, Massachusetts, June 30-July 2, 2004. [Akaike, 1974] H. Akaike, 1974. A new look at the statistical model – validation, IEEE. Transactions on Automatic Control, Vol. AC 19, N° 6, pp. 716-723. [Bai, 1990] E.W. Bai and D. Li, 1990. Convergence of the Iterative Hammerstein System Identification Algorithm, Automatica, Vol. 26, N° 4, pp. 651677. [Bai, 2002] E. W. Bai, 2002. A blind approach to the Hammerstein-Wiener model identification, Automatica, Vol. 38, N° 6, pp. 967-979. [Bai, 2003] E. W. Bai, 2003. Frequency domain identification of Hammerstein models, IEEE Transactions on Automatic Control, Vol. 48, pp. 530-542. [Bai, 2008] E. W. Bai, 2008. Identification of an additive nonlinear system and its applications in generalized Hammerstein model, Automatica, Vol. 44, pp. 430-436. [Bariani, 1988] J. P. Bariani, 1988, Conception et réalisation d’un logiciel de CAO en automatique : identification structurelle et commande PID, Thèse de 3eme cycle, université de NICE. [Bauer, 2002] D. Bauer and B. Ninness, 2002. Asymptotic properties of leastsquares estimates of Hammerstein-Wiener models, International Journal of Control, Vol. 75, N° 1, pp. 34-51. [Bazaraa, 1993] M. S. Bazaraa, H. D. Sherali and C. M. Shetty, 1993. Nonlinear Programming theory and algorithms – second edition, John Wiley and Sons, Inc. [Ben Abdennour, 2001] R. Ben Abdennour, P. Borne, M. Ksouri et F. M’Sahli, 2001, Identification et commande numérique des procédés industriels, Edition Technip, Paris. [Billings, 1981] S. A. Billings and I. J. Leontaritis, 1981. Identification of nonlinear systems using parameter estimation techniques, in Proc. IEE Conference on Control Application, Warwick, U.K, pp. 183-187. 127 [Borne, 1997] P. Borne, 1997. Analyse et régulation des processus industriels, Eyrolles. [Bortolot, 2005] Z. J. Bortolot and R. H. Wynne, 2005. Estimating forest biomass using small foot print LiDAR data: An individual tree-based approach that incorporates training data, ISPRS Journal of Photogrammetry and Remote Sensing, N° 59, pp. 342-360. [Boucher, 1996] P. Boucher et D. Dumur, 1996. La commande prédictive, Edition Technip. [Boujmil, 1991] M. H. Boujmil, 1991, Identification paramétrique et structurelle des systèmes non linéaires discrets, déterministes, stochastiques et monovariables, DEA, ENSET de Tunis. [Boyd, 1985] S. Boyd, L. Chua, 1985. Fading Memory and the Problem of Approximating Nonlinear Operators with Volterra Series, IEEE Transactions on Circuits and Systems, Vol. CAS-32, N° 11. [Brent, 1973] R. P. Brent, 1973. Algorithms for Minimization without Derivatives, Prentice-Hall. [Camacho, 1999] E. F. Camacho et C. Bordons, 1999. Model Predictive Control, Springer-Verlag, London, 1999. ISBN 3-540-76241-8. [Chelouah, 2000] R. Chelouah and P. Siarry, 2000. A continuous genetic algorithm designed for the global optimization of multimodal functions, Journal of Heuristics, N° 6, pp. 191-213. [Chen, 1999] L. Chen and K. S. Narendra, 1999. On a new Approach to the design of tracking controllers for nonlinear dynamical systems, Proceedings of the American Control Conference Sans Diego, California, pp. 3534-3538. [Cherruault, 1999] Y. Cherruault, 1999. Optimisation: Méthodes Locales et Globales”, Presses Universitaires de France, ISBN 2-130-49910-4. [Clarke, 1979] D. W. Clarke, P. J. Gawthop. 1979. Self-tuning control, Proceedings IEEE, Vol. 126, pp. 633-640. [Clarke, 1987] D. W. Clarke, C. Mohtadi et P.S. Tuffs. 1987. Generalized predictive control- part I et II. Automatica, Vol. 23, N° 2, pp. 137-160. [Clarke, 1988] D. W. Clarke, 1988. Application of Generalized Predictive Control to Industrial Processes, IEEE control systems Magazine, Vol. 8, pp. 49-55. 128 [Clarke, 1989] D. W. Clarke and T. Mohtadi, 1989. Properties of Generalized Predictive Control, Automatica, Vol. 25, N° 6, pp. 859-875. [Culioli, 1994] J. C. Culioli, 1994. Introduction à l’Optimisation, Ellipses, ISBN 2729-89428-4. [Cutler, 1980] C. R. Cutler, B. C. Ramaker, 1980. Dynamic matrix control - a computer control algorithm, Automatic Control Conference. San Fransisco, Californie. [Filali, 2002] S. Filali, K. Kemih and A. Kias, 2002. Constrained predictive control using gradient method, Advances in Modelling Analysis -C- AMSE Journal, Vol. 57, N° 3, pp. 35-46. [Fletcher, 1987] R. Fletcher, 1987. Practical Methods of Optimization, John Wiley & Sons, ISBN 0-471-49463-1. [Fnaich, 2000] Fnaich F., 2000. Algorithmes Génétiques Application à l’Automatique, CA2I’2000, 19,20 et 21 Février, Hammamet, Tunisie. [Gesing, 1979] W. Gesing and E. J. Davison, 1979. An Exact Penalty Function Algorithm for Solving General Constrained Parameter Optimization Problems. Automation, Vol. 15, pp. 175-188. Pergamon Press Ltd. [Giannakis, 2001] G. B. Giannakis and E. Serpedin, 2001. A bibliography on nonlinear system identification, Signal Processing, Vol. 81, pp. 533-580. [Glover, 1989] F. Glover, 1989. Tabu Search - Part I, ORSA Journal on Computing, Vol. 1, N° 3, pp. 190-206. [Glover, 1990] F. Glover, 1990. Tabu Search - Part II, ORSA Journal on Computing, Vol. 2, N° 1, pp. 4-32. [Goldberg, 1989] D. E. Goldberg, 1989. Genetic Algorithms in Search, Optimization and Machine Learning, Addison Wesley, ISBN 0-201-15767-5. [Greblicki, 2000] W. Greblicki, 2000. Continuous time Hammerstein system identification, IEEE Transactions on Automatic Control, Vol. 45, pp. 1232-1236. [Haber, 1990] R. Haber, and H. Unbehauen, 1990. Structure identification of Non linear systems–A survey on Input/output Approaches, Automatica, Vol. 26, N° 4, pp. 651-677. [Hazem, 2004] M. A. Hazem and M. M. Bayoumi, 2004. Volterra system identification using adaptive genetic algorithms, Applied Soft Computing, N° 5, pp. 75-86. 129 [Hazem, 2005] M. A. Hazem and M. M. Bayoumi, 2005. An adaptive evolutionary algorithm for Volterra system identification, Pattern Recognition Letters N° 26, pp. 109-119. [Hernando, 1988] D. Hernando and A. A. Desrochers, 1988. Modelling of Nonlinear Discrete-time Systems from Input-Output Data, Automatica, Vol. 24, N° 5, pp. 629-641. [Holland, 1975] J. H. Holland, 1975. Adaptation in Natural and Artificial System, The University of Michigan Press, ISBN 0-472-08460-7. [Hu, 1992] N. Hu, 1992. Tabu Search Method with random moves for globally optimal design, International Journal for Numerical Methods in Engineering, Vol. 35, N° 5, pp. 1055-1070. [Hunter, 1986] I. W. Hunter and M. J. Korenberg, 1986. The identification of nonlinear biological systems: Wiener and Hammerstein cascade models, Biolog. Cybernet., Vol. 55, pp. 135-144. [Kelley, 1999] C. T. Kelley, 1999. Detection and Remediation of Stagnation in the Nelder-Mead. Algorithm using a sufficient decrease condition, SIAM journal on optimization. Vol. 1, pp 43-55, http : // epubs.siam.org/sam-bin/dbq/article/31520. [Kenneth, 1999] E. Kenneth, J. Jer-Nan and S. Richard, 1999. Direct adaptive predictive control using gradient descent, The Journal of the Acoustical Society of America, Vol. 105, Issue 2, February 1999, p.1300. [Ki, 1997] H. Ki and J. Cohen, 1997. Linear and Nonlinear ARMA Model Parameter Estimation Using an Artificial Neural Network, IEEE Trans. on Biomedical Engineering, Vol. 44, N° 3, pp. 168-174. [Kirkpatrick, 1983] S. Kirkpatrick, C. D. Gelatt, M. P. Vecchi, 1983. Optimization by Simulated Annealing, Science Magazine, Vol. 220, N° 4598, pp. 671680. [Kortmann, 1988] M. Kortmann and H. Unbehauen, 1988. Structure detection in the identification of non linear systems, APII, N° 22, pp. 5-25. [Kurpati, 2002] A. Kurpati, S. Azarm and J. Wu, 2002. Constraint handling improvements for multiobjective genetic algorithms, Structural and Multidisciplinary Optimization Revue, Springer-Verlag, Vol. 23, pp. 204-213. [Labrousse, 1977] C. Labrousse, 1977. Statistique,Tome 3, Edition Dunod, Paris. 130 [Landau, 1988] I. D. Landau, 1988. Identification et commande des systèmes, édition Hermès. [Leontaritis, 1985a] I. J. Leontaritis and S. A. Billings, 1985. Input–output parametric models for nonlinear systems—Part I: Deterministic nonlinear systems, International Journal of Control, Vol. 41, N° 2, pp. 303-328. [Leontaritis, 1985b] I. J. Leontaritis and S. A. Billings, 1985. Input-Output Parametric Models for Nonlinear Systems—Part II: Stochastic Nonlinear Systems, International Journal of Control, Vol. 41, N° 2, pp. 329-344. [Li, 1993] C. J. Li and Y. C. Jeon, 1993. Genetic Algorithm in Identifying Nonlinear Autoregressive with Exogenous Input Models for Nonlinear Systems, Proceeding of American Control Conference, San Francisco, CA, June 2-4, IEEE, Piscataway, NJ, pp. 2305-2309. [Li, 1994] C. J. Li and Y. C. Jeon, 1994. A Learning Controller Based on NonLinear ARX Inverse Model Identified by a Genetic Algorithm, Proceeding of Symposium on Intelligent Process Control in ASME International Mechanical Engineering Congress and Exposition, Chicago, IL, Nov. 6-11, Vol. 55, N° 1, pp 447-458. [Liu, 2003] W. Liu, L. Zhengpei, L. Fu and W. Yaqi, 2003. A systematic method to identify nonlinear dynamics of BWR by using the reactor noise, Progress in Nuclear Energy, Vol. 43, N° 1-4, pp. 209-216. [Luersen, 2004] M. A. Luersen and R. Le Riche, 2004. Globalized Nelder-Mead method for engineering optimization, Computers and Structures, N° 82, pp. 2251-2260. [Luh, 1998] G. C. Luh and G. Rizzoni, 1998. Nonlinear System Identification Using Genetic Algorithms with Application to Feedforward Control Design, Proceedings of the American Control Conference, Philadelphia, Pennsylvania, June1998, pp. 2371-2375. [Madár, 2005] J. Madár, J. Abonyi and F. Szeifert, 2005. Genetic programming for the identification of nonlinear input-output models, Industrial and Engineering Chemistry Research, Vol. 44, N° 9, pp. 3178-3186. [Man, 1996] K. F. Man , K.S. Tang and S. Kwong, 1996. Genetic algorithms: Concepts and Applications in engineering design, IEEE Transactions on Industrial Electronics, Vol. 43, N° 5, pp. 519-533. 131 [Mrabet, 2002] M. Mrabet, F. Fnaiech, A. Chaari and K. Al-Haddad, 2002. Nonlinear Predictive Control Based on NARX Models with Structure Identification, IECON 02, 28th Annual Conference of the Industrial Electronics Society, IEEE 2002. Vol. 3, pp. 1757-1762. [Morari, 1994] M. Morari, 1994. Model Predictive Control: Multivariable Control Technique of Choice in the 1990s?, Advances in Model-Based Predictive Control. Oxford University Press. [Narendre, 1966] K. S. Narenda and P. Gallman, 1966. An iterative method for the identification of nonlinear systems using a Hammerstein model, IEEE Transactions on Automatic Control, Vol. AC-11, pp. 546-550. [Narendra, 1990] K. S. Narendra and K. Parthasarathy, 1990. Identification and control of dynamical systems using neural networks, IEEE Transactions on Neural Networks, Vol. 1, N° 1, pp. 4-27. [Narendra, 1997] K. S. Narendra, and S. M. Mukhopadhyay, 1997. Adaptive Control Using Neural Networks and Approximate Models, IEEE Transactions on Neural Networks, Vol. 8, N° 3, pp. 475-485. [Nelder, 1965] J. A. Nelder, R. Mead, 1965. A Simplex Method for Function Minimization, Computer Journal, Vol. 7, pp. 308-312. [Ong, 2005] C. S.Ong, J. J. Huang and G. H. Tzeng, 2005. Model identification of ARIMA family using genetic algorithms, Applied Mathematics and Computation, Vol. 164, Issue 3, pp. 885-912. [Petlenkov, 2007] E. Petlenkov, 2007. NN-ANARX Structure Based Dynamic Output Feedback Linearization for Control of Nonlinear MIMO Systems. Proceedings of the 15th Mediterranean Conference on Control and Automation, Athens-Greece. [Powell, 1965] M. J. D. Powell, 1965. An efficient method for finding the minimum of a function of several variables without calculating derivations, Computer Journal, Vol. 7, pp. 155-162. [Press, 1992] W. H. Press, 1992. Numerical Recipes in C: The Art of Scientific Computing, Cambridge University Press, ISBN 0-521-43108-5. [Propoi, 1963] A. I. Propoi., 1963. Use of linear programming methods for synthesizing sampled-data automatic systems. Automation and Remote Control, Vol. 24, pp. 837-844. [Pukkila, 1988] T. M. Pukkila and P. R. Krishnaiah, 1988. On the use of autoregressive order determination criteria in Univariate white noise 132 tests, IEEE Transactions on Acoustic Speech and Signal Processing, Vol. ASSP-36, pp. 764-774. [Qin, 1997] S. J. Qin, and T.J. Badgwell, 1997. An overview of industrial model predictive control technology, Chemical Process Control-V, edited by Jeffrey C. Kantor, Carlos E. Garcia and Brice Carnahan, pp.232-256, Tahoe, California. [Rao, 1996] S. S. Rao, 1996. Engineering Optimization: Theory and Practice, John Wiley & Sons, ISBN 0-471-55034-5. [Richalet, 1976] J. Richalet, A. Rault, J.L. Testud and J. Papon, 1976. Algorithmic control of industrial processes, 4th IFAC Symposium on Identification and System Parameter Estimation. Tribilsi, URSS. [Richalet, 1978] J. Richalet, A. Rault, J.L. Testud and J. Papon, 1978. Model Predictive Heuristic Control : Application to Industrial processes, Automatica, Vol. 14, N° 2, pp. 413-428. [Rosenbrock, 1960] H. H. Rosenbrock, 1960. An automatic method for finding the greatest or least value of a function, Computer Journal, N° 3, pp. 175-184. [Ruano, 2003] A. E.Ruano, P. J. Fleming, C. Teixeira, K. Rodriguez-Vazquez and C. M. Fonseca, 2003. Nonlinear identification of aircraft gas-turbine dynamics, Neurocomputing Journal, Vol. 55, pp. 551-579. [Sheng, 2001] L. Sheng, K. H. Ju and K. H. Chon, 2001. A new Algorithm for linear and Nonlinear ARMA Model parameter estimation using Affine Geometry, IEEE Transactions on Biomedical Engineering, Vol. 48, N° 10, pp. 1116-1124. [Sheng, 2003] L. Sheng and K. H. Chon, 2003. Nonlinear Autoregressive and Nonlinear Autoregressive Moving Average Model parameter estimation by Minimising Hyper surface Distance, IEEE Transactions on Signal processing, Vol. 51, N° 12, pp. 3020-3026. [Song, 2006] F. Song and P. Li, 2006. MIMO Decoupling Control Based on Support Vector Machines αth-order Inversion, Proceeding of the 6th World Congress on Intelligent Control and Automation. Dalian, China, pp. 1002-1006. [Wen, 2004] Y. Wen and L. Xiaoou, 2004. Fuzzy identification Using Fuzzy Neural Networks with stable learning algorithms. IEEE Transactions on Fuzzy Systems, Vol. 12, N° 3, pp. 411-420. 133 [Westwick, 1996] D. Westwick, M. Verhaegen, 1996. Identifying MIMO Wiener systems using subspace model identification method, Signal Processing Journal, Vol. 52, pp. 235-258. [Wigren,1994] T. Wigren, 1994. Convergence analysis of recursive identification algorithm based on the nonlinear Wiener model, IEEE Transactions on Automatic Control, Vol. 39, pp. 2191-2206. [Wright, 1996] M. H. Wright, 1996. Direct Search Methods : Once Scorned, Now Respectable”. Numerical Analysis 1995, proceedings of the 1995 Dundee Biennial Conference in Numerical Analysis, Edition D.F. Griffiths and G.A. Watson, Vol. 344, pp. 191-208. [Wright, 1998] M. H. Wright, 1998. Optimization methods for base station placement in wireless systems, Proceedings of the IEEE VTC’98 Conference, pp. 387-391. [Yang, 2005] X. H. Yang, Z. F. Yang, G. H. Lu and J. Q. Li, 2005. A grey-encoded, hybrid-accelerated, genetic algorithm for global optimizations in dynamical systems, Communications in Nonlinear science and Numerical Simulation, N° 10, pp. 355-363. [Yang, 2000] Z. J. Yang, T. Fujimoto and K. Kumamaru, 2000. A Genetic Algorithm Approach to Identification of Nonlinear Polynomial Models, IFAC System Identification, Santa Barbara, California, USA, 2000. [Yonghong, 1996] T. Yonghong, A. R. Van Cauwenberghe, 1996. Optimization techniques for the design of a neural predictive controller, Journal Neurocomputing, Vol. 10, Issue 1, pp. 83-96. [Zhu, 1999] Y. Zhu, 1999. Distillation column identification for control using Wiener model, Proceeding of American Control Conference 1999, pp. 3462-3466. [Zhu, 2002] Y. Zhu, 2002. Estimation of an N-L-N Hammerstein-Wiener model, Automatica, Vol. 38, N° 9, pp. 1607-1614. 134 Liste des publications Conférences internationales : • Brahim Tlili, Faouzi Bouani, ‘’Commande prédictive non linéaire’’, Conférence Internationale Francophone d’Automatique (CIFA’2004), Novembre 2004, DOUZ, Tunisie, Art. N°285. • Brahim Tlili, Faouzi Bouani, Mekki Ksouri, ‘’Estimation des paramètres du modèle NARMA à l’aide des réseaux de neurones et des algorithmes génétiques’’, Actes de la 6° conférence francophone de MOdélisation et SIMulation. Avril 2006, Rabat, Maroc. Lavoisier, Doc et Tech. Vol.1, pp.47-53, Mars 2006. http://www.isima.fr/mosim06/actes/articles/1-Metaheuristiques/89.pdf. • Brahim Tlili, Faouzi Bouani, Mekki Ksouri,‘’A derivative-free constrained predictive controller’’, Proceedings of the 10th WSEAS International Conference on Systems, Vouliagmeni, Athens, Greece, July 10-12, 2006 (pp:358-363). • Brahim Tlili, Faouzi Bouani, Mekki Ksouri, ‘’ Identification of multivariable NARMA models using Artificial Neural Networks’’, Industrial Simulation Conference ISC 2008, Lyon, France, June 9-11, 2008. pp. 40-44. Revues scientifiques : • Brahim Tlili, Faouzi Bouani, Mekki Ksouri, ‘’Constrained model predictive control using a derivative-free optimization method’’, WSEAS Transactions on Systems, Issue 10, Vol. 5, October 2006 (pp:2307-2313). http://www.worldses.org/journals/systems/systems-october2006.doc • Brahim Tlili, Faouzi Bouani, Mekki Ksouri, ‘’ Identification des systèmes non linéaires par des modèles de type NARMA’’, à paraître au Journal Européen des Systèmes Automatisés JESA. Vol. 9, Novembre 2008. 135 ANNEXES 136 ANNEXE 1 1) Le F_Test ou le test de Fisher : Le F_Test est utilisé pour comparer deux variances observées afin de savoir si elles diffèrent entre elles de manière significative ou non. La loi de Fisher-Snedecor à n1 et n2 degrés de liberté [Labrousse, 1997] a) Définition Soient S21 et S22 deux variances de deux variables aléatoires indépendantes X et Y, avec X = (x1, x2, …, xi, …, xn1) et Y = (y1, y2, …, yi, …, yn2) La quantité : S 21 / n 1 F= 2 S 2 / n2 (A1.1) est appelée variable aléatoire de Fisher à n1 et n2 degrés de liberté. Les tables permettent de déterminer par simple inspection, pour α donné, et pour n1 et b) Tables de la Loi de Fisher n2 degrés de liberté la valeur X α telle que la probabilité de F supérieure à X α est égale à α (Pr ( F > X α ) = α ), avec α un seuil de signification. En général, on choisit α = 0.05. Ces tables sont données dans le paragraphe 2 de cet annexe. c) Principe du test Il s’agit d’un test statistique obtenu en utilisant un raisonnement par l’absurde : - On suppose que S21 et S22 ne diffèrent pas entre elles de manière significative. Sous cette hypothèse (appelée hypothèse nulle H 0 ) on peut déterminer, au seuil α , un intervalle I de confiance pour F . - Si l’hypothèse était juste, la valeur de F a une probabilité α d’être à l’extérieur de l’intervalle. Si l’on constate que F ∈ I cela est compatible avec l’hypothèse donc on l’accepte. Si par contre l’on constate que F ∉ I alors l’hypothèse faite n’a que α chance sur 100 d’être vraie. Si α est suffisamment petit, on peut, avec un risque faible, rejeter l’hypothèse nulle ( H 0 ). 137 d) Différentes étapes du F_Test On va appliquer les différentes étapes du F_Test sur un critère utilisé pour estimer l’ordre exact n d’un système. Ce critère est défini comme suit [Bariani, 1988], [Boujmil, 1991]: Soient V (n) : la variance du résidu pour l’ordre n, V ( n̂1 ) : la variance du résidu pour l’ordre n̂1 , V ( n̂ 2 ) : la variance du résidu pour l’ordre n̂ 2 . Le critère est basé sur l’indépendance statistique des quantités V ( n̂1 ) et [V ( n̂1 ) - V ( n̂ 2 )]. Pour n̂ 2 > n̂1 ≥ n, la quantité : F= V(nˆ1 ) - V(nˆ 2 ) N − nˆ 2 V(nˆ 2 ) nˆ 2 - nˆ1 (A1.2) a une distribution de Fisher à ( n̂ 2 – n̂1 ) et (N- n̂ 2 ) degrés de liberté avec N est le nombre de mesures et n1 est le nombre de paramètres du modèle. Nous pouvons donc établir six étapes du F_Test utilisées pour estimer l’ordre n : Etape 1 Définir les hypothèses H 0 est l’hypothèse nulle : n1 étant l’ordre exact du modèle (c'est-à-dire il n’y a pas de différence significative entre les deux variances). Etape 2 Définir un paramètre qui sous l’hypothèse nulle obéit à une loi connue. F= V N (nˆ1 ) - V N (n̂ 2 ) N − nˆ 2 ⋅ nˆ 2 - n̂ 1 V N (n̂ 2 ) Etape 3 Définir un seuil α que l’on appellera seuil de signification, soit α = 0.05. Etape 4 Définir une région critique associée à ce seuil : Par exemple pour α = 0.05 on lit sur les abaques de Fisher à (n1 = 1) et (n2 > 120) degrés de liberté une valeur Xα = 3.84. Etape 5 Calculer la valeur du paramètre en fonction des données du problème, soit F cette valeur. Etape 6 Décider (A1.3) Si F ≤ 3.84 on accepte H 0 avec un risque d’erreur inférieur ou égal à 5% (c'est-à-dire que n1 = n, l’ordre exact recherché). Les écarts entre les deux variances sont dus au hasard. 138 2) Table de Fisher-Snedecor 139 ANNEXE 2 42 0.1008 0.8103 -0.3830 -0.0246 -0.1490 0.5467 0.1016 0.1382 -0.0765 -0.0194 0.0094 0.0755 0.0926 0.0058 -0.1685 0.1899 0.0631 -0.0204 -0.0221 -0.0147 0.0130 -0.0117 0.0009 0.0193 0.1151 0.0052 0.0627 0.0090 -0.0119 0.0317 -0.1598 0.0738 -0.0305 0.0187 0.0654 0.0178 -0.0035 -0.0085 0.0317 -0.0000 -0.0542 -0.1153 -0.0047 43 0.0296 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 86 0.0443 0.0911 0.1078 -0.1242 -0.0498 -0.0085 -0.1170 0.1985 0.0948 -0.0352 0.0376 -0.0791 0.0995 -0.1126 0 -0.0132 0.0206 -0.0043 -0.0078 0.0028 -0.0034 -0.2038 0.0521 -0.0282 0.0118 0.0355 0.0162 -0.0031 0.0107 -0.0087 0.0019 -0.0550 0.0386 -0.0493 0.0203 -0.0030 0.4138 -0.0321 -0.0524 -0.0316 -0.0845 -0.0427 0.0109 87 -0.0589 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 130 0.0129 0.0129 0.0675 -0.0002 -0.0411 -0.0088 -0.0013 -0.0343 0.0426 0.0206 -0.0015 0.0030 -0.0509 0.0722 -0.0354 0.0150 -0.0020 0.0115 -0.0040 -0.0012 0.0003 0.0031 0.1239 -0.1606 -0.0945 -0.0096 -0.0137 0.1104 0.1694 0.1526 0.0387 -0.0424 0.1047 -0.0103 -0.0135 -0.0444 0.0346 0.0033 -0.0500 0.0680 -0.0047 -0.0262 0.0130 131 0.0016 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 174 -0.0694 -0.0524 -0.0444 -0.0094 -0.1216 0.0021 0 0.0043 -0.0080 -0.0303 -0.0055 -0.0260 -0.0076 -0.0039 0.0257 0.0243 0.0231 -0.0097 -0.0299 0.0276 -0.0399 0.0096 0.0132 -0.0096 0.0121 -0.0211 0.0070 -0.0004 0 -0.0058 -0.0035 0.0137 -0.0075 0.0189 0.0086 -0.0122 0 -0.0100 -0.0508 0.3694 0.2845 0.0607 -0.0308 175 -0.1426 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 218 -0.2917 -0.0061 -0.0540 0.0820 -0.1573 0.0199 0.0079 -0.0125 -0.0161 -0.1048 0.1340 0.0593 0.0979 0.0056 -0.1437 -0.2177 -0.0714 0.0318 0.0149 -0.0377 0.0449 0.0019 -0.0358 0.0077 0.0436 0.0893 0.0190 -0.0023 -0.0018 0.0680 -0.0199 -0.0035 0.0164 -0.0094 0.0086 -0.0171 0.0142 0.0066 -0.0222 0.0085 -0.0019 0.0033 -0.0011 219 0.0016 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 Tableau A2 : Les paramètres du modèle NARMA estimés par la méthode (RN), (système 4. Chapitre I). 140 ANNEXE 3 Liste des fonctions benchmarks F1 (1 variable): F1(x) = sin(x) + sin(2x/3) Domaine de recherche: 3.1<x<20.4. 3 minima locaux minimum global: x* = 17.0393; F1(x*) = -1.90596. F2 (1 variable): F2(x)=sin(x) + sin(10x/3) + Ln(x) – 0.84x Domaine de recherche: 2.7<x<7.5. 3 minima locaux 1 minimum global: x* = 5.19 ; F2(x*) = -3.8717. F3 (2 variables): F3(x1,x2) = x1² + x2² - cos(18x1) – cos(18x2). Domaine de recherche: (-1,-1) < (x1,x2) < (1,1) A peut près 50 minima locaux 1 minimum global: (x1*,x2*) = (0,0); F3(x1*,x2*) = -2. F4 (2 variables): F4(x1,x2) = (x1² +1)² + (x2²+1)² - 2 (x1 +x2+1)². Domaine de recherche: (-3,-3) < (x1,x2) < (3,3) 1 minimum global: (x1*,x2*) = (1.3247,1.3247); F4(x1*,x2*) = -11.45851. 141 Tracé de la fonction F1 2 1.5 1 F1 (x) 0.5 0 -0.5 -1 -1.5 -2 4 6 8 10 12 14 16 18 20 x Tracé de la fonction F2 0 -0.5 -1 F2 (x) -1.5 -2 -2.5 -3 -3.5 -4 3 3.5 4 4.5 5 5.5 6 6.5 7 7.5 x 142 Tracé de la fonction F3 0 F3 (x1, x2) -0.5 -1 -1.5 -2 1 0.5 1 0.5 0 0 -0.5 -0.5 -1 x2 -1 x1 Tracé de la fonction F4 -6 F4 (x1, x2) -7 -8 -9 -10 -11 -12 2 0 -2 x2 -3 -2 0 -1 1 2 3 x1 143