FR2900753A1 - Procede d'optimisation automatique d'un reseau de transport de gaz naturel - Google Patents

Abstract

Le procédé d'optimisation automatique est appliqué à un réseau de transport de gaz naturel en régime permanent comprenant à la fois un ensemble d'ouvrages passifs tels que des canalisations (101 à 112) ou des résistances, et un ensemble d'ouvrages actifs comprenant des vannes de régulation (31, 32) , des vannes d'isolement (51), des stations de compression (41), des dispositifs de stockage ou d'alimentation (21, 22), des dispositifs de consommation (61 à 65), des éléments (41A) de dérivation des stations de compression (41) et des éléments (31A, 32A) de dérivation des vannes de régulation (31, 32), les ouvrages passifs et les ouvrages actifs étant reliés entre eux par des jonctions (1.1 à 1.13). Le procédé d'optimisation comprend la détermination de valeurs pour des variables continues. On choisit comme état initial de l'optimisation des intervalles de valeurs pour les variables continues et des ensembles de valeurs pour les variables discrètes. On explore les possibilités de valeurs pour les variables en construisant au fur et à mesure un arbre avec des branches reliées à des noeuds décrivant les combinaisons de valeurs envisagées en utilisant une technique de séparation des variables et d'évaluation, les valeurs des grandeurs recherchées étant considérées comme optimales lorsque des contraintes prédéterminées ne sont plus transgressées ou sont transgressées au minimum et un critère économique prédéterminé constituant une fonction objectif est minimisé.

Description

La présente invention a pour objet un procédé d'optimisation automatique

d'un réseau de transport de gaz naturel en régime permanent, le réseau de transport de gaz naturel comprenant à la fois un ensemble d'ouvrages passifs tels que des canalisations ou des résistances, et un ensemble d'ouvrages actifs comprenant des vannes de régulation, des vannes d'isolement, des stations de compression avec chacune au moins un compresseur, des dispositifs de stockage ou d'alimentation, des dispositifs de consommation, des éléments de dérivation des stations de compression et des éléments de dérivation des vannes de régulation, les ouvrages passifs et les ouvrages actifs étant reliés entre eux par des jonctions, le procédé d'optimisation comprenant la détermination de valeurs pour des variables continues telles que la pression et le débit du gaz naturel en tout point du réseau de transport, et la détermination de valeurs pour des variables discrètes telles que l'état de démarrage des compresseurs, l'état d'ouverture des stations de compression, l'état d'ouverture des vannes de régulation, l'état des éléments de dérivation des stations de compression, l'état des éléments de dérivation des vannes de régulation, l'orientation des stations de compression et l'orientation des vannes de régulation.

La présente invention vise à permettre de déterminer notamment les valeurs optimales de pression et de débit en tout point d'un réseau de transport de gaz naturel en régime permanent. L'invention vise également à permettre de déterminer de façon optimale et automatique non seulement des variables continues, telles que le débit, qui peuvent prendre toutes les valeurs comprises dans un intervalle, mais également des variables discrètes qui ne peuvent prendre qu'un nombre fini de valeurs.

Adtre d'exemple, l'ouverture d'une vanne est une variable discrète, dès lors que cette vanne ne peut être qu'ouverte (ce qui peut être représenté par exemple par un 1) ou fermée (ce qui peut alors être représenté par un 0).

Le procédé selon l'invention vise ainsi à permettre de déterminer de façon automatique et de manière optimale notamment des facteurs tels que l'ouverture des vannes, le démarrage des compresseurs, l'orientation des ouvrages actifs (station de compression et vannes de régulation), l'état des éléments de dérivation (by-pass) de ces ouvrages actifs, voire l'adaptation série ou parallèle de certains compresseurs.

POL.r déterminer par le t.eu! les oaraotërütguss d'un réseau Je transport gaz, q~~!!e ~ue soi! !a s~Uon ~hysi~~ on x~.`~/ en / aditionn:5.ément iv m* b~ 's es (enr. iëea loi i'}nohhoff S'cii.L=, d'un graphe cOmp"^^ (somne 51 une re/ztionodentée entre deux noeuds. L.= o.5a poosederd un a5h~u1 "ETAT" qui indique si l'arc est activé ou désactivé. 2 Selon la loi des noeuds, il y a pour tous les noeuds du réseau égalité entre la quantité de gaz entrant dans un noeud et la quantité de gaz sortant de ce noeud et globalement tout ce qui entre dans le réseau doit en sortir. De manière u , selon la loi des noeuds, on obtient le système d'équations linéaires suivant :

Mnmu= Enonsmmmations + Eressources + Cntatmns avec B : matrice d'incidence du réseau traduisant la correspondance entre les arcs et les nœuds du réseau, s==: vecteur des quantités s'écoulant en chaque arc, iU E , : vecteur des quantités livrées aux consommations, E" "=: vouteu, dee quanhtéo émisae uu injeotéeu au, wunumna (dispooitifs d omokage ou d'alimentation), Cmm* : vecteur des quantités de gaz carburant consommées par les stations de compression. 15 La loi des noeuds permet de définir ainsi un système d'équations linéaires. Toutes les quantités entrant ou sortant sont algébriques et leur signe est défini par le choix d'une convention. On peut considérer que ce qui entre dans un noeud est positif tandis que ce qui en sort est négatif.

Selon la loi des mailles, la somme , le long d'une maille, des 70 différences de pression du gaz entre deux noeuds consécutifs est nulle. La loi des mailles permet ainsi de définir un système d'équations : =8 maiile avec 4P: différence des pressions entre deux noeuds consécutifs d'une maille. 25 Parce que l'on connaît les formules de perte de charge dans les canalisations selon la forme suivante : P,2 ûP22 =ox{]xkd. on peut exprimer également de façon équivalente la loi des maiUasà l'aide de différences de pression au carré : =O maille avec /P2 : différence des pne/ns .~~ues entre deux nœuds consécutifs d'une maille. La loi des mailles vué/me d'équaiior^ aires. 30 con L' )1-oc uu~ !e nox'x ~ Si l'on considère un réseau de N noeuds et M mailles, on en déduit le nombre d'arcs égal à N + M - 1, auxquels correspondent autant d'équations indépendantes, à savoir N-1 équations selon la loi des noeuds et M équations selon la loi des mailles. Les deux lois de Kirchhoff permettent de déterminer des débits (dans la mesure 5 où la loi des mailles remplace les différences de pressions quadratiques par leur expression équivalente fonction du débü, en Qénénsl de la forme sPz= oxQo où o est considéré constant.

Lorsque le système d'équations de ces deux lois est résolu, les débits sont partout connus et l'imposition d'une pression particulière en un noeud quelconque du 10 réseau permet de connaître les pressions en tous les noeuds.

De manière traditionnelle, les procédés de simulation visant à déterminer les variables continues en tout point d'un réseau comprennent une première phase de résolution des deux lois de Kirchhoff et d'obtention des débits partout et une deuxième phase d'imposition d'une pression en un noeud particulier et d'obtention des pressions 15 partout.

Généralement, le processus itère plusieurs fois entre la phase n 1 et la phase n 2 car les coefficients a intervenant dans les relations de la loi des mailles ne sont pas parfaitement constants et dépendent très légèrement des pressions et des débits.

Cette démarche impose deux fortes restrictions. La première restriction est 20 qu'elle ne s'applique que sur des réseaux ne comportant que des canalisations ou, de façon plus générale, des ouvrages passifs. En effet, les ouvrages passifs présentent une relation entre la différence des pressions en amont et en aval de l'ouvrage et son débit. Cette relation est l'équation de perte de charge proprement dite. Disposant de cette relation, il est toujours possible de remplacer les différences de pressions par leur 25 expression fonction du débit. En nevenohe, un ouvrage aotif, tel qu'une vanne de régulation ou une station de compression, ne présente pas nécessairement une telle relation ou du moins, si cette équation existe, elle contient au moins une inconnue supplémentaire.

Les ouvrages actifs constituent des organes de commande du réseau en 30 introduisant des inconnues supplémentaires telles quo, par exemple, le degré d'ouverture d'une vanne de régulation. Connaissant le degré d'ouverture et considérant un certain nombre de coefficients caractéristiques fournis par le constructeur, on relie les pressions en amont, en aval et le débit àoepo' çeniaged,ouverturo.

Pour les stations de compre. onnue introdL est la puiss.1ce motrice 3~ impression puissance dépensée compression, ,uisqI le nette dernière est ro au débit et au taux de oompress~.,rapport entre sa po en avs- ion en am(' procédés d'ouvrages estent à l'utilisateur de lui-même la valeur cie ces inconnues. 40 Implicitement, !esovruOes actifs ne le sont alors plus car ils présent en~ ule véritable 4 équation de perte de charge (ou de gain s'il s'agit de compression). Typiquement, le contournement proposé par ces procédés consiste à demander à l'utilisateur d'imposer soit la puissance de compression s'il s'agit d'une station de compression, soit le degré d'ouverture de la vanne s'il s'agit d'une détente, etc. L'imposition de ces grandeurs établit 5 un lien entre le débit de l'ouvrage et ses pressions amont et aval. Disposant alors d'une telle relation, il est donc possible de résoudre la seconde loi de Kirchhoff. Toute la difficulté consiste à déterminer quelle puissance des stations de compression ou quel degré d'ouverture des vannes de régulation imposer. Il n'est pas toujours possible, du moins en un temps raisonnable, de trouver manuellement selon une 10 approche essais/erreurs un jeu de valeurs convenables notamment pour un réseau complexe où les mailles sont interconnectées les unes avec les autres. La seconde restriction est la nécessité d'imposer une pression en un noeud particulier de la phase n" 2. Du fait de la première restriction, le réseau est supposé composé uniquement d'ouvrages passifs. En imposant cette pression particulière et 15 après avoir résolu les deux lois de Kirchhoff, les pressions peuvent être connues partout.

Si le réseau ne comporte qu'une unique source, il paraît naturel d'imposer la pression au noeud particulier qu'est le noeud de cette source. En général, on impose la pression la plus haute possible en ce point et l'ensemble des pressions en tous les noeuds constitue alors le régime de pression maximal. Une autre approche est de choisir 20 au noeud source une pression la plus basse possible dans la limite où les pressions en tous les nœuds ne soient pas inférieures à un seuil fixé. L'ensemble des pressions en tous les noeuds constitue alors le régime de pression minimal.

Si le régime de pression minimal présente des pressions supérieures au régime de pression maximal, cela signifie que l'on ne peut trouver une pression au nœud source 25 qui, à la fois :

soit inférieure à la pression maximale de ce noeud, ù soit supérieure à une valeur limite qui permet de satisfaire tous les seuils de pression minimale en tous les noeuds. Le réseau est dit saturé.

30 Dans le cas d'un réseau ne comportant qu'une seule source, le débit de celle-ci injecté sur 'e réseau est partestement déterminé par la loi des noeuds. Ce n'est plus le cas si une seconde acp nlud,^ correspo lébit de cette seconde s( par le nombre de ne.s-us ri Lst s une inconnue additionnelle 1troduite dans le problème. est p/éseniea/' -Lradhionne sieur, que !es ùeu~ soc/cas s~:~ ~ p@rumarana|isa1ionde langueur n0xset de dame:re très grand. L'introduction de cette nouvelle maille poae, dans le système d'équations, l'équation suppLimentaire manquante. L'équiiibre entre le nombre d'inconnues du problème et le nombre d'équations est alors rétabli. En général, la canalisation fictive est construite de telle façon qu'elle présente une loi de perte de charge constante (P` = Cste). Avec cette canalisation fictive, l'imposition d'une pression en une seule des deux sources du réseau permet de connaître les pressions en tous les noeuds du réseau.

Cette méthode présente toutefois l'inconvénient que si la constante de la loi de perte de charge de la canalisation fictive est trop importante, alors la résolution de la seconde loi de Kirchhoff conduit à trouver un débit qui sort du réseau pour la seconde source, ce qui peut ne pas être souhaitable puisqu'il s'agit d'une source, autrement dit d'une entrée de gaz sur le réseau.

L'approche du calcul de réseau dans toute sa généralité par simulation n'est donc pas satisfaisante puisque la recherche des valeurs optimales des puissances et des pressions et des débits à imposer doit être prise en charge manuellement. Pour remédier à ces inconvénients, on a déjà proposé de disposer d'un nombre d'inconnues supérieur au nombre d'équations, de sorte qu'il existe plusieurs solutions au problème posé et que le choix d'une solution particulière va s'opérer selon un critère donné, qui détermine une optimisation.

Certains procédés connus sont toutefois conçus pour le calcul de réseaux en régime dynamique plutôt qu'en régime permanent.

D'autres procédés d'optimisation pour le calcul de réseaux en régime permanent ou dynamique imposent des conditions et contraintes particulières qui rendent ces procédés incomplets ou peu souples.

La présente invention vise à remédier aux inconvénients précités et à permettre de déterminer automatiquement de façon optimale tous les degrés de liberté d'un réseau de transport de gaz en régime permanent, avec minimisation d'un critère économique et non transgression des contraintes, ou transgression des contraintes a minima.

L'invention vise plus particulièrement à effectuer une hybridation d'une méthode d'optimisation combinatoire et continue afin de déterminer les valeurs de l'ensemble des variables discrètes et continues, de façon entièrement automatique. Ces buts sont atteints, conformément à l'invention, grâce à un procédé d'optimisation automatique d'un réseau de transport de gaz naturel en régime permanent, le réseau de transport de gaz naturel cornp = ant à 'a fois ensemi d'ouvrages passifs 'mis r e de cana i_ tEo d~S ~3I ~rh 1 i ~~ ...tifs yuiüüvi i, co vct ico ü Vwi' I ic compres_._e avec chacune au moins un compresseur, res disnc stuc rie etoekeria d'alin latic des dispositifs de consommation, des été ~ts de c n us passi s e? . rue, _.iies entre eux ~~s cuoos, le procédé d'opti ri sxt c mprenant la détermination de valeurs pour des varàb'es continues telles que la pression et le débit du gaz naturel en tout point du réseau de transport, et la détermination de valeurs pour des variables discrètes telles que l'état de démarrage des compresseurs, l'état d'ouverture des stations de compression, l'état d'ouverture des vannes de régulation, l'état des éléments de dérivation des stations de compression, l'état des éléments de dérivation des vannes de régulation, l'orientation des stations de compression et l'orientation des vannes de régulation, caractérisé en ce que l'on choisit comme état initial de l'optimisation des intervalles de valeurs pour les variables continues et des ensembles de valeurs pour les variables discrètes, en ce que l'on explore les possibilités de valeurs pour les variables en construisant au fur et à mesure un arbre avec des branches reliées à des noeuds décrivant les combinaisons de valeurs envisagées en utilisant une technique de séparation des variables, c'est-à-dire de découpage conduisant à la génération de nouveaux noeuds dans l'arbre, et d'évaluation, c'est-à-dire de détermination avec une probabilité forte des branches de l'arbre qui peuvent conduire à des feuilles constituant une solution finale optimisée, de manière à parcourir en priorité ces branches à plus forte probabilité de réussite, les valeurs des grandeurs recherchées étant considérées comme optimales lorsque des contraintes prédéterminées ne sont plus transgressées ou sont transgressées au minimum et un critère économique prédéterminé constituant une fonction objectif est minimisé, cette fonction objectif étant de la forme

g = a x Régime + R x Energie + y x Cible avec : a, l et y sont des coefficients de pondération.

Régime représente un facteur de minimisation ou de maximisation de la pression en des points donnés du réseau tels que tout point aval d'un dispositif de stockage ou d'alimentation, tout point amont et tout point aval d'une station de compression ou d'une vanne de régulation, et tout point amont d'un dispositif de consommation, Energie représente un facteur de minimisation de la consommation d'énergie de compression,

Cible représente un facteur de maximisation ou de minimisation du débit d'un tronçon du réseau situé entre deux jonctions ou de la pression d'une jonction particulière, et lesdites contraintes prédéterminées comprenant d'une part des contraintes d'égalité comprenant la loi de perte de charge dans les canalisations et la loi des noeuds régissant le calcul des réseaux, et d'autre part des contraintes d'inégalités comprenant des contraintes de débit minimal et maximal, des contraintes de pression minimale et maximale des ouvrages actifs ou passifs, des contraintes de puissance de compression des stations de compression.

Plus généralement, le problème de configuration optimale des ouvrages actifs est modélisé sous la forme d'un programme P, d'optimisation se présentant sous la forme suivante : 7 2 ImIn(x,s,e) f (x,$) -=g(x)+uxI "~/.r~./ x)=sa E R".a/ e Rp'oe E Rq.ee avoo; x est l'ensemble des variables des des débits 0 et des pressions P, ,_fonction objectif constituant _critère économique --,__iù, C(x) est l'ensemble des p contraintes d'inégalité linéaires et non linéaires sur les ouvrages actifs, 5 f est une matrice dont les coefficients sont nuls ou égaux aux valeurs maximales des contraintes, e est le vecteur des variables binaires, de dimension p pour que l'équation le faisant intervenir soit cohérente mais le nombre de variables binaires est réellement : 3 x le nombre d'ouvrages actifs, CE(x) est l'ensemble des q contraintes d'égalité linéaires ou non linéaires, s est une variable d'écart qui, lorsqu'elle est non nulle, représente la transgression d'une contrainte, 10 a est un coefficient représentant le degré d'autorisation de transgression de contraintes. Selon un mode particulier de réalisation, les variables sont représentées par des intervalles, la technique de séparation des variables est appliquée aux seules variables discrètes et des bornes de la fonction objectif sont calculées en utilisant l'arithmétique 15 des intervalles.

Selon un autre mode particulier de réalisation, les variables sont représentées par des intervalles, la technique de séparation des variables est appliquée à la fois aux variables discrètes et aux variables continues, la séparation comprenant le découpage de l'espace de définition des variables continues, l'exploration s'effectuant séparément sur 20 des parties de l'ensemble réalisable et l'intervalle de variation de la fonction objectif étant évalué sur chacune de ces parties.

Dans ce cas, avantageusement, lors de l'exploration des possibilités de valeurs pour les variables avec une technique de séparation des variables et d'évaluation, on établit d'abord une liste de noeuds à explorer triés selon un critère de mérite M calculé 25 pour chaque noeud, tant que la liste de noeuds à explorer n'est pas vide, pour chaque noeud courant, on évalue si ce noeud courant peut contenir une solution, si c'est le cas, on découpe l'intervalle correspondant à la variable considérée selon une loi de séparation pour établir une liste de noeuds fils, pour chaque noeud fils on évalue des bornes minimale et maximale de la fonction objectif et on évalue si le noeud fils peut améliorer la 30 situation courante, si c'est le cas on effectue une propagation de la contrainte sur ses variables, si la propagation ne conduit pas à des intervalles vides, on évalue des bornes minimale et maximale de la fonction objectif et on vérifie qu'}l n'y a pas d'impossibilité à ce que le noeud fils contienne au moins une solution faisable, on effectue un test pour déterminer a'il y a encore des valeurs diaonbh o non inmhunciéea, c'est-à-ùira des 35 va, ladies D our L uelles une valeur , se ..- n'a pu être e, à le mérite du noeud ae mente

ùinn L. ..ee,ier en oeuvre les .nAya:

Le critère de mérite M est tel qu'un noeud est exploré en priorité lorsqu'il présente la plus faible borne minimale de la fonction objectif. jour la meilleure soiu. rer dans la liste d urante s'il y a !i( lies, triée seloi 8 Lors des tests d'élimination des noeuds ne pouvant pas contenir l'optimum, on met en oeuvre !'une des méthodes consistant à utiliser la monotonie de la fonction objectif, à oti!iaer un test de contraintes transgressées ou à utiliser un test de valeur d'objectif moins bonne que la valeur courante. 5 Lors de la séparation d'un nœud courant en noeuds fils, on divise le domaine de variation d'une ou plusieurs variables choisies selon des critères basés sur le diamètre d'intervalles rattachés aux variables. Le procédé comprend en outre un critère d'arrêt basé sur le temps d'exécution ou sur l'évaluation de certains diamètres d'intervalles. 10 En complément de la propagation des contraintes, on procède à une actualisation de la borne maximale de l'optimum de la fonction objectif en utilisant les conditions d'optimalité du problème d'optimisation dites de Fritz-John. Lorsqu'à un noeud du procédé de séparation et d'évaluation, toutes les variables discrètes ont été instanciées, on met en oeuvre en outre un processus d'optimisation non 15 linéaire basé sur une méthode de points intérieurs. De façon alternative, à chaque noeud du procédé de séparation et d'évaluation, on met en outre en oeuvre un processus d'optimisation non linéaire basé sur une méthode de points intérieurs.

D'autres caractéristiques et avantages de l'invention ressortiront de la description 20 suivante de modes particuliers de réalisation, donnés à titre d'exemples, en référence aux dessins annexés, sur lesquels :

û la Figure 1 est un schéma-bloc montrant les modules principaux d'un système d'optimisation automatique d'un réseau de transport de gaz selon l'invention ; 25 û la Figure 2 est une vue schématique d'un exemple de partie de réseau de transport de gaz ; la Figure 3 est une vue schématique d'un exemple de configuration d'une station de compression se situant sur un point d'interconnexion d'un réseau de transport de gaz ;

30 û !2Figure4es1un u800hëma\iciL [ a sàon ! ~ez ~u8deséps

xartie de réseau à

__ 6 usaax~n~ d _~eada pression d'initialisation pour différents nceu Js de la partie réseau de la figure 5 ; 9 ù la Figure 7 est un tableau donnant des exemples d'intervalles de débit d'initialisation pour différents arcs de la partie de réseau de la Figure 5 ; ù la Figure 8 est un tableau donnant les résultats de tests effectués sur la partie de réseau de la figure 5; 5 ù la Figure 9 est un tableau donnant des résultats des intervalles de pression pour les différents noeuds de la partie du réseau de la Figure 5 dans les cas du tableau de la Figure 8 où la propagation n'est pas stoppée ; la Figure 10 est un tableau donnant les résultats des intervalles de débit pour les différents arcs de la partie du réseau de la Figure 5 dans les cas du lO tableau de la Figure 8 où la propagation n'est pas stoppée ; la Figure Il est un organigramme illustrant un exemple de mise en oeuvre du procédé d'optimisation selon l'invention ; ù la Figure 12 est un diagramme montrant un arbre de calcul qui représente la propagation/rétro-propagation de contraintes ; et 15 _ la Figure 13 est une vue schématique d'un exemple de réseau de transport de gaz naturel auquel est applicable l'invention. La présente invention s'applique d'une manière générale à tous les réseaux de transport de gaz, notamment de gaz naturel, même si ces réseaux sont très étendus, à l'échelle d'un pays ou d'une région. De tels réseaux peuvent comprendre plusieurs 20 milliers de canalisations, plusieurs centaines de vannes de régulation, plusieurs dizaines stations de compression, plusieurs centaines de ressources (points d'entrée du gaz sur le réseau) et plusieurs milliers de consommations (points de sortie du gaz sur le réseau).

Le procédé selon l'invention vise à déterminer automatiquement tous les degrés de liberté d'un réseau en régime permanent et ceci de manière optimale.

25 Les valeurs sont optimales dans le sens où les contraintes ne sont pas transgressées et un critère économique est minimisé ou, si cela n'est pas possible, les contraintes sont transgressées au minimum.

Les degrés de liberté sont les pressions, les débits, les démarrages de compresseurs, les états ouvert/fermé, en ligne/en dérivation (by-pass) et les orientations 30 directe ou inverse des ouvrages actifs.

Pou n réseau ree~ o!ue entaines de variables à valeur entière xarexempi urouvert [em,~ sus des quelques milliers de varit rnno

=' rend possible le lancement de cobo!en série, c'est- 35 à-dire sans intervention hL. .a. Ce caractère autonome du calcul est d'un intérêt majeur dans un contexte de réseaux pouvant donner lieu à une mutpcté de scénarios d'acheminement. 10 La Figure 1 est un schéma-bloc illustrant les principaux modules mis en oeuvre dans le cadre de la définition d'un réseau de transport de gaz. Le module 5 constitue un modeleur qui est un ensemble permettant la modélisation du réseau. On entend par là sa description physique via ses ouvrages et sa structure (sous-réseaux connexes, blocs de pression,...). Ce modeleur inclut également de préférence des moyens de simulation (ou d'équilibrage) du réseau en débits et pressions, Le module 8 constitue par sa part un coeur de calcul permettant d'assurer une optimisation du réseau. 10 Le module d'optimisation 8 comprend essentiellement un solveur 6 dont les fonctions (notamment mise en oeuvre de la technique de séparation des variables et d'évaluation) seront explicitées plus loin et un solveur non linéaire convexe 7 qui peut agir en complément du solveur 6.

La Figure 2 montre de façon schématique une parlie de réseau de transport de 15 gaz comprenant différents points de prélèvements de gaz pour des consommations locales C. A chaque point de prélèvement est associée une contrainte de pression dépendant des besoins de consommation.

La partie du réseau de transport comprend également des points d'approvisionnement en gaz pour fournir au réseau du gaz à partir de ressources locales 20 R qui peuvent être par exemple des réserves de gaz stockées dans des cavités souterraines.

La capacité du tronçon de réseau dépend à ta fois du niveau des consommations C et des mouvements d'approvisionnement à partir des ressources R.

Dans un réseau de transport de gaz, la pression du gaz décroît progressivement 25 au cours du transit. Afin que le gaz soit acheminé avec un respect de la contrainte de pression admissible pour le consommateur, le niveau de pression doit être remonté régulièrement à l'aide de stations de compression réparties sur le réseau.

Chaque station de compression comprend au moins un compresseur et compte généralement de 2 à 12 compresseurs, la puissance totale des machines installées 30 pouvant être comprise entre environl MW et 50 MW. ne doit pas dépasser la pression =iguna 3 U!ust,e un exemple de contigu l'une o!`sa' / oompxeoso~. lc/ : n un point d'inAs:zo.~nexi~o 1~ ~.~ /pas n.ëre vision int point 1.r, Une J'api 1--., i sur quelle est placee vanne de JaL régu!adonde pression 30 est ' let lient raccordée au point d'inten:onne^k` 1.0. Un ou plusieurs compresseurs 40 sont disposés sur une troisième canalisation qui prend naissance au point d'interconnexion ou jonction 1.0. 11 Selon un exemple de réalisation typique, on peut avoir une pression de 51 bar sur la première canalisation 100, une pression de 59 bar sur la deuxième canalisation en amont de la vanne de régulation 30, une pression de 51 bar sur la deuxième canalisation en aval de la vanne de régulation 30 et une pression de 67 bar sur la troisième 5 canalisation en aval des compresseurs 40.

La présente invention vise à optimiser automatiquement les mouvements de gaz sur des réseaux complexes, le procédé offnan1à la fois une grande robustesse et une grande précision.

Dans la suite de la description, on considérera que l'expression "ouvrege actif " 10 englobe les vannes de régulation et les statïons de compression ainsi que les vannes d'isolement, les ressources et les stockages.

L'expression ouvrage passif recouvre les canalisations et les résistances.

Le procédé selon l'invention a pour but de rechercher les réglages adéquats pour les ouvrages actifs et d'établir une carte de débits et de pressions du réseau, pour 15 optimiser un critère économique.

Le critère économique est composé des trois termes différents :

ù le régime de presmion: minimise ou maximise les pressions en aval desstockages et ressources, en amont et aval des stations de compression et des vannes de régulation et en amont des consommations,

20 ù l'économie : minimise la consommation d'énergie de compression,

la cible : maximise ou minimise le débit d'un arc ou la pression d'un noeud particulier.

Dans le problème mathématique d'optimisation, ce critère est appelé fonction objectif. Dans cette fonction, chaque terme est pondéré par un coefficient (a, 13 et y) qui 25 lui donne plus ou moins d'importance :

g =oxRégime +BxEne i +yxcible Les degrés de liberté sont :

ù les pressions en chaque lise_

ù les ,its en chaque

` nui peuvent prendre toutes les valeurs co" Les ,e s de liberté sont : l'ouverture/la fermeture des ouvrages actifs, dan 12 ù le by-pass des stations de compression et des vannes de régulation, ù l'orientation des stations de compression et des vannes de régulation, ù le démarrage des compresseurs, pour les paramètres discrets ou variables discrètes, qui ne peuvent prendre qu'un 5 nombre fini de valeurs. Le but est de trouver les valeurs des variables qui minimisent le critère économique. La recherche des valeurs des variables est soumise à des contraintes de différents types : ù des contraintes d'égalité : loi de perte de charge dans les canalisations, loi des noeuds. Ces contraintes sont intrinsèques au réseau, on ne peut donc les transgresser ;

ù des contraintes d'inégalités : contraintes de débit minimal et maximal, de pression minimale et maximale des ouvrages, contraintes de puissance de compression des stations, contraintes de vitesse minimale et maximale du gaz à chaque noeud, contraintes déprimogènes pour les vannes de régulation et pour les stations de compression, contraintes de pompage et de gavage des turbocompresseurs, contraintes de pressions de refoulement minimale et maximale des compresseurs, contraintes d'énergies journalières minimale et maximale des consommations, etc. Ces contraintes sont inhérentes aux ouvrages du réseau ou liées aux contraintes contractuelles du réseau (exemple : pression minimale pour un client); elles donnent des limites à ne pas dépasser, mais certaines peuvent être transgressées. Mathématiquement, ces contraintes sont de deux types : linéaires ou non linéaires.

25 Pour effectuer une modélisation d'un réseau de transport de gaz dans sa totalité, on peut considérer qu'à chaque état d'un ouvrage actif correspond une variable binaire e (qui prend la valeur 1 lorsque l'état est actif ou 0 dans le cas contraire, par exemple 1 pour ouvert et O pour fermé). On peut ainsi modéliser le choix entre chacun des états uniquement avec des contraintes linéaires. Le principe est illustré ci-dessous dans le cas 30 d'une station de compression. Exemple pour une station de compression : 10 15 20 Soit pressions Ple triplet variables continues débi ?t de tien de cor ession. 35 es 4 va~/a~!es bha/ s e'- ; lUX 4 ë1s~s ~|~erno~~S ' inverse - /x* poevan am Vui/e S. /muaRé/nen 1. . Cf (x), . C.(x)` (,` ` e~ 4 contraintes pour ces 4 états disjomo".a. Par exon pie. pour l'état direct, Cu,x) est le vecteur des contraintes de débits minimal et maximal, de taux de compression minimal et maximal et de puissances minimale et maximale. 13 Soient Cf rnax, Cbmax, Cbmax'max une estimation des valeurs maximales de ces contraintes quel que soit x. Dans l'exemple de l'état direct, Cd ,ax est le vecteur des débits mininnal et nnuxima!, des taux de compression noinimal et maxinnal et des puissances minimale et maximale.

5 Les contraintes linéaires vont donc s'écrire sous la forme :

ù CJx\~(1 -ef).{f m^"

ù C~hd~(1'eb).Cualma

ù C~(X) 5- (1 -eb).C*ma"

ù Ci(x)~(1 'ej.C, "am

10 ef + eb + ed + e, = 1 afin d'assurer le choix d'un et d'un seul des 4 états.

En partant de ce principe, on peut également effectuer une modélisation en ne conservant que les 3 variables eb, ed, ei, ce qui réduit la combinatoire.

Ces variables vont s'intégrer dans les contraintes de la manière suivante : ù Cf(x) ~(ou+*u+ei).{fmuu

15 ù C^hd≤. (1 eb).Cb maya ù Ci(x) ≤ (1 'eJ.Ci mu, eu+eu+e,~1 afin d'assurer le choix entre un des 4 états, l'état fermé correspondant aux 3 variables nulles.

20 Ainsi le problème de configuration optimale des ouvrages actifs est modélisé sous la forme d'un programme d'optimisation mixte (associant variables continues et variables binaires) non linéaire (car une partie des contraintes Cf(x).Cb(x)^Cd(x)'C,ùdcst non linéaire).

Le programme général peut donc s'écrire sous la forme suivante : 1_1 utitère (exemple : le coût de :az carburant les stations de 25 14 C,(x), l'ensemble des contraintes linéaires (contraintes de bornes) et non linéaires sur les ouvrages actifs ; ces contraintes sont des contraintes d'inégalité au nombre de p, f3, un vecteur dont les coefficIents sont nuls ou égaux aux valeurs maximales des contraintes, e le vecteur des variables binaires, . de dimension p pour que l'équation le faisant intervenir soit chre'e mais le nombre de varia ss L rs uréellement : 3 x le nombre d'ouvrages actifs, ,. ensemble des contraintes à ,éai s (exemple : loi des nœuds), et non linéaire ~xample: équations de perte de charg des canalisations). Elles sont au nombre de q. Le procédé selon l'invention vise à fournir une réponse quel que soit l'état de 10 saturation du réseau. C'est-à-dire qu'il est demandé au procédé d'autoriser, s'il ne peut faire autrement, de transgresser certaines contraintes afin de proposer un résultat, même en cas de saturation. L'autorisation de violation des contraintes est tempérée puisqu'on cherchenaà la minimiser et qu'elle conduira de toute façon à un message de saturation. Prenant en compte cette demande, on modifie légèrement l'écriture du problème en 15 introduisant les variables s qui, si elles sont non nulles, représentent la transgression des contraintes. f (ua)=g(x)+ux C(x\+8e~s, p . ' Cs( ) =ss xsR".s,e~n.aceRn,ee 60,1 aveo: x est l'ensemble des variables des débits Ci et des pressions P, 20 C g(x) est la fonction objectif constituant le critère économique d'optimisation, i(x) est l'ensemble des p contraintes d'inégalité linéaires et non linéaires sur les ouvrages actifs, p est un vecteur dont les coefficients sont nuls ou égaux aux valeurs maximales des contraintes, e est le vecteur des variables binaires, de dimension p pour que l'équation le faisant intervenir soit cohérente mais le nombre de variables binaires est réellement : 3 x le nombre d'ouvrages actifs,, CE(x) est l'ensemble des q contraintes d'égalité linéaires ou non linéaires, 25 oest une variable d'écart qui, lorsqu'elle est non nulle, représente la transgression d'une contrainte, a est un coefficient représentant le degré d'autorisation de transgression de contraintes. Notons que, à variables binaires fixées, le programme PI, qui n'est pas rigoureusement équivalent à Po, a une solution proche de Po si le coefficient a est choisi 30 suffisamment grand puisqu'alors les variables d'écart si et SE sont recherchées très fortement proches de O. C'est un problème combinatoire de taille importante puisqu'il compte plusieurs centaines de variables à valeur entière en sus des quelques milliers de variables continues. 35 Cette mixité du type dos var i' ocnÜnueCes~ ~ounqmoi e deux types d'optimisation sont de préférence au final une solution exacte. Le pur._ [e:ve en pr m~. ^ chDi de 40 os~~rotion daa v~~~~~~^ =. , Dn, par la suite . caet-F:'ttliia "na.~i~ x" L, "m~ . Cette techniqu~ une ukaase de méthodes d'optimisatio~ capables detraiter ééupmbiUmOS OU interviennent des variables discrètes. Le caractère discret d'une variable s'oppose au caractère cofltinu: +d8 ' obmisa5nn cornb/natoi.a 6t tiqu es > de ta Dride afin d obtenir 15 û une variable continue peut prendre n'importe quelle valeur dans un intervalle donné. Dans le cadre du makzul de réseau, ce sera le cas des pressions exprimées en bar, par exemple : P E [40,80], une variable discrète ne peut prendre qu'un certain nombre de valeurs. H 5 s'agit souvent de variables binaires qui représentent par exemple le sens de fonctionnement dune station de compression par exemple x = O (sens direct) oux=1 (sens indirect). La méthode B&B est une méthode arborescente et consiste, au fur et à mesure qu'on construit l'arbre, à réduire le domaine de variation des variables. Cette méthode est 10 couramment utilisée pour obtenir le minimum g/obal d'un problème d'optimisation mettant en jeu des variables binaires. Pour résoudre avec la méthode B&B un problème mixte qui est un problème traitant à la fois des variables discrètes et continues, plusieurs variantes peuvent être envisagées :

B8B1 : la méthode B&B ne sépare que sur les variables binaires. Les variables sont représentées par des intervalles. On pourra alors calculer les bornes de la fonction objectif en utilisant l'arithmétique des intervalles. B&B2: la méthode B&B sépare à la fois sur les variables binaires et continues ; ceci implique une représentation par intervalles. Dans ce cas, le principe de séparation (branch) va consister à découper l'espace de définition des variables continues au lieu de fixer les variables discrètes à l'une de leur valeur. Ainsi, on va explorer séparément des parties de l'ensemble réalisable et borner (bounding) l'intervalle de variation de la fonction objectif sur ces sous-parties.

25 La mise en place d'une méthode B&B de séparation des variables et d'évaluation nécessite donc le choix de stratégies concernant :

la sélection du noeud à examiner : en fonction de la date d'arrivée des noeuds dans la pile, de leur positionnement ou de la valeur d'une fonction de mérite calculée en chaque noeud candidat,

30 l'évaluation des bornes de la solution courante qui permet d'avancer dans la méthode B&B,

l'élimination des pouvant pas contenir l'optimum (test de contraintes transgressé. de valeur d'objectif moins hennie que la valeur courante. utilisation de la monotonie de la fonction objete -iine de varie en d'une ou plu el anables (chdsi - des c'itëree basés sur le diamètre d'intervalles rattachés à la (aux variable(s), le diamètre ou la largeur d'un intervalle correspondant à la différence entre sa borne maximale et sa borne minimale), 15 20 16 le critère d'arrêt basé sur le temps d'exécution ou sur l'évaluation de certains diamètres.

Pour le problème de la configuration optimale des ouvrages actifs, les méthodes B&B consistent à fixer progressivement l'état des ouvrages actifs, et àévaluer àchaque étape, parmi ces combinaisons partieUes, celles qui sont susceptibles de mener à la combinaison globale la plus favorable.

Un exemple va être décrit en référence à la Figure 4.

[)n considère un réseau de gaz sur !mquol on aurait plusieurs stations de compression. On cherche, par exemple, à minimiser le gaz carburant sur le réseau, Si 10 l'on choisit au début de l'arbre de B&B, la station de compression n 1 et que l'on teste la variable binaire associée à son état (e =1). est la borne minimale de la fonction objectif calculée au noeud i, sachant l'ensemble des décisions qui ont déjà été prises. ~ est la borne maximale de la fonction objectif associée à la meilleure 15 combinaison d'états connus lors de l'exploration du noeud i. Si f! > fil:a), (avec fr =0 ), alors il est sûr que la station 1 orientée en sens indirect ( = 0 ) ne conduit pas à la solution optimale. En revanohe, si /', ~ /'l l'exploration continue en fixant une autre variable binaire. Toutes les variables binaires sont ainsi fixées progressivement. Si on n'effectue 20 aucune coupe dans une branche, on obtient une configuration réalisable, c'est-à-dire que l'on a fixé l'ensemble des variables binaires et que l'on respecte l'ensemble des contraintes.

Différentes techniques peuvent être associées à la technique de séparation des variables et d'évaluation.

25 On peut en particulier utiliser la propagation de contraintes qui permet d'exploiter les informations de l'équation ou de l'inéquation pour diminuer les intervalles des variables de cette équation.

Oh ' àre que l'équation non linéaire C(x) et, de manière g/Th cherche avec : 17 La propagation de contraintes peut être basée sur la construction d'un arbre de calcul qui représente C(x). Dans un premier temps, on calcule à partir des feuilles de l'arbre, qui sont les variables et les constantes, la valeur des noeuds intermédiaires et du noeud racine oornespondan1é la valeur de la contrainte (ce qui équivaut à appliquer les 5 règles de l'arithmétique d'intanxsUe), puis on propage la valeur de l'intervalle de la uontoainte, depuis la racine de l'arbre vers les feuilles pour réduire les espaces de définition des variables, L'algorithme de propagation d'une contrainte sur ses variables est le suivant : Etape 1, pn: construction de l'arbre de calcul de la contrainte C, les feuilles sont les variables intervalles x, ou des constantes réelles, dans chaque noeud est stocké le résultat de l'opération partielle et unitaire qu'il représente, par exemple Xa + Xb, le dernier calcul est effectué à la racine. û Etape2. rétro-propagation : descente de l'arbre de la racine vers les feuilles. Aohaque noeud, on tente de réduire le résultat partiel calculé en 1. Par exemple : xa +xu=[a,b] xa :=([a.b]-xb) nxa et xb :=([a.b]'o) nzu La Figure 12 illustre la propagation/rétro-propagation des contraintes pour l'équation suivante donnée à titre d'exemple :

2x3x2 +xi =3avec xi =[1'3]`pour ie (1'2,31 20 La première étape de l'algorithme est présentée dans la partie gauche de la Figure 12 : en partant des valeurs des variables et des constantes, on effectue chaque opération unitaire constituant l'expression jusqu'à obtenir la valeur du membre de gauche de l'expression en haut de l'arbre ; ce noeud est le noeud racine. La seconde étape de l'algorithme est explicitée par la partie droite de la Figure 25 12 : on veut le membre de gauche égal à une valeur précise, on redescend donc l'arbre à partir de la racine, grâce aux opérations inverses à celles utilisées dans la première partie, on cherche à réduire les intervalles de chaque noeud et particulièrement celui des variables. Dans l'exemple, la propagation a permis de réduire chaque intervalle des variables de [1.3] à [1 ,1], c'est-à-dire que les variables ont été instanciées à 1, 30 uniquement grâce à la propagation.

L'algorithme de propagation sur l'ensemble des contraintes d'un problème s'effectue de la fan,û, .sä 'aree bbn, dans une file tric oe[on un ohbèr de o 2. !Bouc e s_,r !a file des contraintes Ta file n'est pas vide 10 15 35 18 Extraction de la "meilleure" contrainte C (pour le critère MI Propagation de C. Si la propagation a conduit â un Intervalle vide pour au moins une variable i Sortie de la boucle . il n'y a pas de solution au problèTe Sinon { Pour chaque variable modifiée par la propagation de C Pour chaque contrainte où intervient cette variable ( Si la contrainte n'est pas déjà résolue, ajout dans la file 15 Selon un exemple de réalisation, seul l'âge de la contrainte intervient dans le critère de mérite M, i.e. la file est équivalente dunepile BFO. Cependant, on peut utiliser un critère plus complexe. Par exemple, une variable très réduite par la propagation d'une contrainte pourrait conduire à insérer les contraintes dans lesquelles elle intervient dans 20 la file avec un mérite élevé. On notera qu'une contrainte est dite résolue lorsqu'elle ast déjà vérifiée quelles que soient les valeurs que les variables prennent dans leurs intervalles (autrement dit si l'intervalle résultant de la propagation sur la contrainte ne contient que des valeurs 25 Pour une contrainte C de fonction d'inclusion C(X) = C~~~~C C est résolue si: s!: ù C est une contrainte d'égalité et C(X) = 0, ù C est une contrainte d'inégalité et C(X). 0. ù C est une contrainte d'inégalité négative et C(X) ~O. 30 Lorsqu'une contrainte est résolue, sa propagation ne conduira plus à aucune réduction des intervalles de ses variables. La technique de propagation de contraintes peut être utilisée par exemple pour déterminer l'orientation des ouvrages actifs des réseaux de transport de gaz. Les ouvrages actifs peuvent être simplement considérés orientés en sens direct lorsque le une r ou4variables t ires comme indiq Jlus ha _ivre de la technique de la propagation des contraintes peut êt ffectué Je d'un ibliothèque d'arithmétique d( propagat. .. traiter 40 dis. méfhodaode propagation des contraintes peuvent d'une part servir à réduire ùaocnbi. i dans des temps réduits, lors d'une première étape pouvant précéder un processus d'optimisation, exact ou approché, et d'autre part être intégrées aux méthodes 10 19 B&B pour permettre en chaque noeud un meilleur calcul des bornes de la fonction objectif et éventuellement des coupes supplémentaires.

En particulier, dans ce dernier cas où la propagation des contraintes s'effectue au sein d'un noeud de l'arbre de recherche et est utilisée pour élaguer les noeuds qu'elle 5 permet de déclarer infaisables, et pour diminuer le diamètre des intervalles des variables, on considère dans la file initiale des contraintes à propager, les contraintes où interviennent la ou les variables dont la séparation a conduit à la création du noeud en cours d'évaluation. Si ce noeud est la racine de Y rbre, alors on place toutes les contraintes dans la file.

10 A titre d'exemple de mise en oeuvre d'une technique de propagation des contraintes, on se référera aux Figures 5 à Io.

Si l'on se reporte à la Figure 5, on voit un réseau simple de transport du gaz comprenant une ressource R, une consommation C, un premier compresseur CP1 et un deuxième oompreooeurCP2. Le réseau comprend des noeuds No à N4 (jonctions ou 15 points d'interconnexion) et des arcs ! à vil (canalisations ou tronçons comportant les compresseurs CPi, CP2, la ressource R et la consommation C).

Le réseau définit cinq variables de pression aux noeudoNo à N4 et sept variables de débit dans les arcs 1 à vll.

La Figure Sdonne un exemple d'intervalles de pression d'initialisation (en bar) 20 aux différents noeuds No à N4.

La ressoumeA a une pression de consigne de 40 bar. C'est pourquoi son intervalle d'initialisation est un intervalle de largeur nulle.

Le noeud N4 de consommation a une pression minimale de livraison de 42 bar, d'où l'initialisation dans l'intervalle [40, 60].

25 La Figure 7 donne un exemple d'intervalles de débit d'initialisation (en m3/h) dans les arcs 1 à vil.

La ressource R et la consommation C correspondant aux arcs / et Nl ont des débits imposés à 800 000 "`3/h. Leurs intervalles sont donc initialisés à des intervalles de largeur nulle.

30 Les arcs Ill ~~7t '-~~aen~ar e IV et VI oone ,ndp~~ sama!iaetior

P!~a

A. 'L.t toutes les des 35 compresseurs CP1. CP2 (tests Al à A4). 20 B. On laisse libre l'orientation du compresseur CP1 et on fixe celle du compresseur CP2 (tests B1 et B2). C. On laisse libres les orientations des deux compresseurs CPi, CP2 (test C). Les résultats de ces tests Ai àA4. B1. B2 et C sont présentés sur le tableau de 5 la Figure 8. Dans les trois cas où la propagation n'est pas stoppée (tests Al, B1 et C), on obtient les résultats identiques présentés sur les tableaux des Figures 9 et 10. La Figure 9 indique les intervalles de pression résultats (en bar) aux différents noeuds N3 à Na. 10 La Figure 10 indique les intervalles de débit résultats (en m3/h) pour les différents arcs 1 à vil.

Dans ces exemples, on voit que les informations contenues dans les contraintes sont utilisées pour réduire les intervalles des variables et permettent également de fixer la valeur de certaines variables discrètes (ici l'orientation de chaque compresseur). En 15 particulier, on voit que si l'orientation d'un ou des deux compresseurs est laissée libre, en appliquant la seule méthode de propagation des contraintes, on peut conclure que le compresseur libre doit s'orienter dans le sens direct. La méthode de propagation des contraintes ainsi que la méthode de séparation des variables et d'évaluation (B&B) fait appel au calcul par intervalles dont les 20 caractéristiques principales seront rappelées ci-dessous.

En arithmétique par irtenvallem, on ne manipule plus des nombres, qui approchent plus ou moins fidèlement une valeur, mais des intervalles contenant cette valeur. Par exemple, on peut tenir compte d'une erreur de mesure en remplaçant une valeur mesurée x avec une incertitude e par l'intervalle [x-e.x+s]. On peut également 25 remplacer une valeur par sa plage de validité telle qu'une pression P d'une ressource représentée par un intervalle [4, 68] bar. Enfin, si l'on désire obtenir un résultat valide pour tout un ensemble de valeurs, on utilise un intervalle contenant ces valeurs. En effet, l'objectif de l'arithmétique par intervalles est de fournir des résultats qui contiennent à coup sûr la valeur ou l'ensemble cherché. On parle alors de résultats garantis, validés ou 30 encore certifiés. ~rar!~ / ~cbeur inbama!de 8P. <~~ st /R3 à un lent de droite, à un rectangle et à un parallélépipède. Un vecteur intervalle est ace- ,n hyper-parallélépipède. Par la suite, on Comme aë ennentpas ua ^ ~.ou ensemble des interv vecteur intervalle ur !A 35 une in erres- -D corre ind respectives- el pësenr ~ss enxa!iec qui ne en p.ùieurs dimensions un ;antes sont des intervalles et '-rand adma "oun-e..` nput !pn 21 utilisera indifféremment les termes de vecteur d'intervalles, de pavé, de boîte ou même On note les objets intervalles par des caractères gras : x. On note x le minimum _ n 5 X ≤ Y<=> xi ≤ y, pour i 1..n. On note w(x) la largeur de x, (avec w pour width) ou encore son diamètre: vv(x)=x _x Le centre midbdet son rayon rad(x) sont définis par: mid(x)=~-~û2 oû~ w(x) rmju)=---=---2 2 Une fonction F : lRn û* IR est une fonction d'inclusion do/our X e !R Si Xe X alors j(X) F(X). On désigne par l'adjectif << ponctuel un objet numérique usuel (c'est-à-dire un nombre réel, ou un vecteur, une matrice de nombres réels) et on le confond avec 15 l'intervalle de diamètre nul. Le résultat d'une opération entre deux intervalles x et y est le plus petit intervalle (au sens de l'inclusion) contenant tous les résultats de l'opération appliquée entre tous les éléments x de x et tous les éléments y de y, c'est-à-dire contenant l'ensemble : 20 (xVy; x e o, ye y) De même, le résultat d'une fonction F(z) est le plus petit intervalle contenant l'ensemble : h(z) ze o~ Si l'on oon~ mrS~'a~i~~oore~~ ' " x ., / [Yi~=I.A+~~+y] [Y, YI_ b~ ] x[y, y]=[nin(Icxy.x xy.)s x y.xx y).nlaxblv y.x x!~~x ~~x xyD ~ _ _ _ 10 22 ùc

e )max(~~ x ~sd8e [~x [~x] ù~u ' \ lKrnex(x^.x )]sinon _ / _ 1 ~x~~=[ninM/~_).max~ /~)] oiOe[x, n _ [x,njx(1 /[y, }]) siOe [y, y] [~[Ic,ù~ ~f~ ]=~~~,v~]si0~x Les propriétés algébriques traditionnelles (c'est-à-dire en arithmétique ponctuelle) telles que la réciprocité entre l'addition et la soustraction ou la distributivité de la multiplication par rapport à l'addition ne sont plus vérifiées : ù la soustraction n'est plus la réciproque de l'addition. En effet : 10 _ la division n'est plus non plus la réciproque de la multiplication, par le même raisonnement que ci-dessus, on obtient : _ la multiplication d'un intervalle par lui-même n'est pas égale à l'élévation au carré. Prenons l'exemple où x = [-3,2] 15 x xx =[-SQ] 2 ù la multiplication n'est pas distributive par rapport à l'addition. Prenons x=[-2.3], y=[1.4] et z =[-2.1] : x x (y + z) = [-1 0,15] 20 x x y + x xz =[14'16] ù La multiplication est en fait sous-distributive par rapport à l'addition, c'est-àdire que : x x (y+o) cxxy+xxz peut aussi finir des 1ctions élémentaires telles que sinus. 25 !'exponem/eUe. etc. prena in arguù es intervalles. Pour cela, on utilise la définition abstraite oi~'~~~ une fonction mot de~u~~ [zoU8nnent les formules permettant de es ca' 23 Par contre, on ne sait définir les fonctions élémentaires que sur des intervalles contenus dans leur domaine de définition : par exemple, le logarithme ne sera défini que pour des intervalles strictement positifs. L'arithmétique par intervalles permet de calculer avec des ensembles et d'obtenir 5 des informations générales et précieuses pour l'optimisation globale d'une fonction. Pour éviter qu'il y ait une surestimation des résultats, il est préférable d'utiliser pour la fonction à prendre en compte une expression dans laquelle chaque variable n'apparaît qu'une fois. Différentes méthodes de sépavariables et d'évaluation (13&13) utilisant 10 l'arithmétique des intervalles seront décrites ci-dessous. Une méthode B&B peut être caractérisée en 5 étapes : 1. sélection : choix du noeud à examiner, 2. évaluation des bornes (bounding), 3. élimination : destruction des noeuds ne pouvant pas contenir l'optimum, 15 4. séparation : construction de 2 noeuds fils en divisant le domaine de variation d'une variable, 5. critère d'arrêt. Différentes solutions peuvent être choisies pour ces 5 étapes afin d'améliorer la qualité du procédé. 20 Soit le problème d'optimisation miny,xif(X). Le vecteur d'intervalles de dimension n, X IRn, est la zone de recherche. La fonction f : Rn-->R est la fonction objectif. {]n note j* le minimum global du prob!ëme. X^ un point optimal tel que KXi =/*, et l'ensemble de ces points X* : minx xl(X) et X* =[K e X f(X) = 25 {]n note les objets intervalles par des caractères gras : x. On note x le minimum _ n de x et x son maximum. On a alors x =- [ x , x ], et on considère l'ordre partiel sur xi ≤ y, pm/r!=1.n. On note x) ' w(x)=xùx Le centre mid(x) et son rayon rad(x) sont définis par : 24 mk(o)=~-~-2 x _~ w(x ~dx)=---'=--- 2 2 Une fonction F : /R"-+ !H est une fonction d'inclusion de/'aur Xe !An Si Xe X alorsf(X) F(X). 5 Voici différentes règles de sélection du noeud à examiner parmi la liste de noeuds en attente. Bien sûr, ces stratégies peuvent être combinées : par exemple la stratégie "Best-fint est souvent combinée avec la stratégie "D!dastfint comme second critère a'il yades ex-oequn. 1. Plus ancien en premier ( Oldest-first ) 10 Cette stratégie consiste à examiner en premier le noeud qui a été créé le plus tôt. 2. Plus profond en premier ( Depth-First ) Cette stratégie consiste à examiner en premier le noeud qui se trouve au niveau le plus profond de l'arbre, Le. le noeud qui a le plus d'ascendants. 3. Meilleur en premier ( Best-First ) [Règle de Moore-Skelboe] 15 Cette stratégie consiste à privilégier !enoeudquioornaspnndaup!uspeth F(( ' Le. celui qui a la plus faible borne inférieure de l'optimum. 4. Index de rejet ( Reject Index ) a. Optimum connu Définissons pour chaque noeud correspondant au vecteur d'intervalles X le 20 paramètre : ~ &/^~)=~---~~~ ~' w(F(X)) Remarquons que si w(F(X)) est nul, alors il ny a pas lieu dévaluer t puisque le noeud ne sera pas découpé. Le noeu électionné est alc ~e!ui qui correspond à la plus forte valeur de pf''i pf. tre nécessite cependant de connaître optimum à ' avance, ce qui luoi s"o du"F ë\~déf" oo~ . b. Optimum estimé 25 La variante du paramètre t quand on ne connaît pas l'optimum à ravance peut s'écrire : pf (f,,X)= w(F(X)) où k est l'indice de l'itération considérée. L'indice k correspond globalement au 5 nombre de nœuds examinés et fk est une approximation del* à l'itération k. Remarquons que la règle Best First n'est donc jamais qu'un cas particulier de êf pour!equedfk =~~ En efhet, si Yoea1 l'intervalle du noeud présentant la plus faible borne inférieure de F ( meilleur noeud ), on a o!orepf(Yo) =Ost pf négatif pour tous les autres noeuds. 10 D'autres possibilités pour fi, peuvent être : ou encore : c. Avec contraintes 15 Pour un problème contraint de la forme :

(min f (X) 0,i =1..p !Xe Rn

les Reject Index définis ci-dessus ne tiennent pas compte du tout descontraintes et risquent de sélectionner des nœuds qui présentent de bonnes valeurs de p/ mais conduiront à des noeuds infaisables. 20 Certains auteurs proposent donc de définir un index de faisabilité construit de la façon suivante.

Pour une contrainte Ci et pour un noeud correspondant à un domaine de variation 25 Dans !e cas où v,=D, ksksidade la contrainte 1 peut être décidée directement, et pao(X) peut ë're fixé à 1 si }{ satisfait CÏ, -1 sinon. Remarquons que si puo(X) < O. alors X. de manière certaine, ne satisfait pas C/ car C/(X 0 . Au uon1naine, si ' 26 puci(X) =1. alors C(X)≤ 0 et donc X satisfait Ci de manière certaine. Dans tous les autres cas, l'état de transgression de C, est indéterminé. Pour les X qui ne sont pas certainement infaisables , c'est-à-dire pour lesquels i=1..p,yoc(X). 0, définissons un index de faisabilité global pour l'ensemble des p pu(X ) = pom(X) Ainsi construit, cet index global possède 2 propriétés:

ù pu(X) = 1 X est certainement faisable , Le. faisable de manière certaine

ù po(X) [0,1] <=:i X est indéterminé.

10 Cela permet alors de définir un index de rejet modifié intégrant l'index de faisabilité :

ounKfk,XO=naQQxnKfk,XO Sino(X) =1,ie. si X est certainement faisable , alors on se ramène au simple Reject Index . En revanche, si X est indéterminé, ce nouvel index prend en compte le 15 degré de faisabilité de X. Cela permet de définir une nouvelle règle de sélection de nooud: on sélectionne le nœud avec la plus grande valeur de pupf. Un dernier critère permet d'hybrider le critère pupf avec le critère "best-first " classique basé sur la valeur de F(X) : pgpf6~({,x) pupf (fk,X) ≠ 0 20 avec M une valeur très grande fixée préalablement.

En effet si =0, alors soit 8f(fk'X) = 0, ce qui signifie ù dans le cas où J =f- que l'on est sûr de ne pas améliorerf; ooüyu(fk.X) =Omsquisignifie qu'il existe au moins une contra nte telle que C,(X) = 0 . De tels X ne semblent pas très prometteurs. Uoio n ès grande, 25 sidérer / intenam retape d'évaluatio ~~erAu r!eS_- celles das oon~cs~n~eo s'U y a. Pour les nlélhcdes de 8&b u![!can( !nr]:m~Jque d'inten/aUos, les fonctions d'inclusion sont généralement obtenues par extension naturelle des fonctions usuelles. Exemple : Si/': x-*,o'o`sù s=[-5,2j. alors F : x -+"u-e est une fonction d'inclusion de/' sur x avec : i n ( ~-o u-e `x ).mex(u 'u )j ai8ebx,x -n~ex(xe.-=~s~on et e` elX`! [eX, ` = =e ] Pour l'étape d'élimination, plusieurs méthodes sont possibles. 1. Test de faisabilité Si le problème est un problème soumis à p contraintes d'inégalité Ci : min f(X) i=1..z 10 Soit Ci une fonction d'inclusion de la contrainte Ci. A chaque examen d'un noeud correspondant au domaine de variation de X, les p contraintes Ci(X) sont évaluées. Si 3 ie(1'p)/[-=,0nC/(X)=0.alors Ueet sûr que le noeud ne peut contenir de solution faisable. H peut donc être élagué. 2. Test de coupure << Cut Off Test>>) 15 C'est le critère le plus simple et le plus connu d'élimination : il s'agit de rejeter tous les noeuds pour lesquels /~/< FQQ, où f est la borne supérieure courante de 3. Test de point milieu (' Middle Point Test ) Certaines publications ne font pas de distinction entre le "Cu1 Off Test et le 20 Middle Point Test (MPT). Le MPT ne serait en fait qu'un moyen supplémentaire de calculer une borne supérieure de/ Le "Cut Off Test " consiste à prendre initialement comme borne supérieure F(X) puis de l'actualiser àchaque division d'intervalle. Pour un p/5b!ènne contraint. ( aSotiC n'/at posslb!e que /c ,u'on sait que con5ent au m n poit f~ Da~s [ on n5 qui est me 25 supérieure de !optimum. S'il s'agit d'un problème contraint. 3t toutefois nécessaire de s'ass e midtX) est un point faisable. Pour ur -e [-on contraint, si la fonc ion objectif est strictement monotone par rapport à la composante xi d'un vecteur d'intervalle X. alors l'optimum ne peut se 30 trouver à l'intérieur de xi, Pour déterminer si f est strictement monotone par rapport aux 27 c -z 5 28 composantes de X, on évalue les n composantes de la fonction d'inclusion du gradient de _f sur X. Si pour i, l'intervalle rësu:rni ne contient pas la valeur 0, alors f est strictement monotone par rapport Au/. Dans ce cas, on peut réduire la composante xi à un réel : xi se rédubà x si la 5 1èrne composante de la fonction d'inclusion du gradient as ,u~ n intervalle qui a une borne supérieure strictement négative, et xi se réduit à x, si la composante de la fonction d'inclusion du gradient est un intervalle qui a une borne inférieure strictement positive. Pour l'étapde séparationméthodes sont également envisageables : 1. Bissection sur une variable 10 Dans toutes les règles suivantes, on sélectionne la variable j qui maximise une fonction de mérite D. On sépare donc sur la variable j telle que j=arg(maxi=I,nD(i)). a. Plus grand diamètre La fonction de mérite est ici simplement le diamètre de la variable: D(i)= oi). La difficulté d'utiliser cette fonction de mérite est liée à la nécessité de s'affranchir des 15 facteurs d'échelles. Par exemple, si l'on traite un problème de calcul de réseau, il faudra bien échelonner les variables pour pouvoir comparer les diamètres des pressions avec ceux des variables binaires. Pour pouvoir contourner cet obstacle, une règle proche de celle-ci et ne faisant pas non plus intervenir d'information sur les dérivées peut être définie : (ojsjOex/ 20 Di w( ,) ,~"~~ ~xvu x/ avec mig(X) =min, m!xiDnpournaüud!iser!anlagnbude: mag(X) =rnaxx m!xi Cette variante permet ainsi de normaliser le diamètre des intervalles considérés. b. Règle de Hansen b3. , ~mw ~Fest la Doai /o/ux/u/ Jsio ent idë, de sépa ~ur!avr~~/ nao au/ 25 29 D(i)=w[(xi ùmkj(xi)x\7F,(X)] L'idée sous-jacente est ici de réduire le diamètre de w(F(X)) qui, après calcul, se ramène à la somme sur toutes les directions du terme D(i). d. Règle de Ratz bis 5 L'idée sous-jacente est la même, mais on va jusqu'au second ordre : DU)= x/ùmkKojx Vf(nid(xj)+ où FG est l'élément de coordonnées (i,k) de la matrice des dérivées secondes (hessien) clef. Pour des méthodes qui calculent de toute façon le gradient et le hessien, par 10 différentiation automatique, cette règle n'est pas beaucoup plus coûteuse que les autres.

2. Multi-section a. Multi-section statique

Jusqu'ici, nous avons considéré qu'a partir d'un noeud, 2 noeuds fils étaient créés grâce à la bissection du pavé X e !Rn dans une direction unique. Cependant, il peut être 15 pertinent de retenir plusieurs directions de séparation. Par exemple, on peut couper l'intervalle de variation de chaque variable en 2, 2n noeuds fils sont alors créés. On peut aussi couper l'intervalle d'une direction en 3 parties, créant ainsi 3 noeuds fils, ou encore les intervalles de 2 variables en 3, créant 32 fils, etc.

b. Multi-section adaptative

20 Notons (a) la règle de plus grand diamètre présentée en 1.a. (b) la règle qui sépare les intervalles de toutes les variables en 2, (c) la règle qui sépare les intervalles de toutes les variables en 3.

Une règle hybride (adaptative) utilisera 3 paramètres Pl, P2 et pf pour déterminer quelle règle utiliser.

25 Les paramètres pl et pz sont deux ses qui seront à ajuster. pf est le Reject Index défini plus haut, et est une fonction -tu :Donsidéré. Les noeuds qui aL~~ Reject < p- seront séparés suivant a régie eux tels que i ron és suivant la règle (b) et ceux tels qui

30 L = ua|!= .' (e peu, tout à fait être définie à partir de variantes de pf, telles que ae:.ni plus haut par exemple.

Différents critères d'arrêt peuvent être utilisés. /ùmid(xj)

30 1. Diamètre de la zone de recherche

Un critère d'arrêt peut être l'examen d'un noeud N tel que w(X) s c où X est l'intervalle de variations des variables pour N. Bien sûr, cela suppose un bon échelonnement des variables. 2. Diamètre de la fonction obiectif

Un critère d'arrêt peut être l'examen d'un noeud N tel que w(F(X)) s c où X est l'intervalle de variations des variables pour N.

3. Temps maximal d'exécution

Un critère d'arrêt complémentaire peut être un temps maximal d'exécution au-10 delà duquel on arrête l'algorithme, quels que soient les résultats obtenus. Un critère d'arrêt de ce type est nécessaire en complément éventuellement d'un autre pour éviter des explorations trop longues.

On décrira maintenant en référence à la Figure 11 un exemple d'organigramme illustrant la méthode de B&B (séparation des variables et évaluation) et de propagation 15 des contraintes appliquée dans un solveur pour rechercher une solution optimale et exacte dans le cadre de la configuration d'un réseau de transport de gaz.

Pour mettre en oeuvre cette technique, une bibliothèque d'intervalles est mise en place pour permettre la gestion des variables exprimées sous forme de nombres ou d'intervalles.

20 Par ailleurs des procédures de différentiation automatique basées sur des arbres de calcul permettent de calculer les valeurs des dérivées premières et secondes à partir d'une expression mathématique.

Des moyens sont également mis en oeuvre pour effectuer le calcul de développements de Taylor aux ordres 1 et 2.

25 Dans l'organigramme de la Figure 11, les étapes 201, 202 et 203 correspondent à des étapes globales du procédé de B&B, tandis que les étapes 204, 206, 208, 211, 212, 214 sont appliquées à chaque pas du procédé de B&B. Les références 205, 207, 209, 210 correspondent à des tests aboutissant à une réponse oui ou non qui permet le choix de la pr-=mua, à suivre.

30 C particul corresp..d au choix de l'arbre à explorer. Létape 202 consi; en une séparation en ne is fils. L'étape comprend une série d'opérations effectuées ur chaque noeud fils.

ainsi ce 203, ilyad'abor_ à u e t ce phis 205 est ensL e pue. Si la r (Il, ù aur à l'éta;~e 203 pou: traüer un autre noeua ti s. Si la réponse au test 205 est Il y a passage à une etape 206 de propagation, étro-propagation telle que celle proposée par exomple par F. Messine. 31 Après l'étape 206, on effectue un nouveau test d'élagage 207. Si la réponse est oui, il y a retour à l'étape 203, si en revanche la réponse est non, on peut passer directement à un autre test 210, mais selon un mode de réalisation préférentiel, on procède d'abonjà l'étape 208 à une résolution du système d'optimalité de Fritz-John, qui sera décrite plus en détail plus loin. En sortie de l'étape 208, un nouveau test d'élagage 209 permet de revenir à l'étape 203 si la réponse est oui et de passer au test 210 si la réponse est non (absence d'élagage).

Le test 210 permet d'examiner si les variables discrètes sont toutes instanciées Si toutes les variables discrètes ne sont pas toutes instanciées, on passe à une étape 211 de mise àjnur éventuelle de la meilleure solution, puis à une étape 212 de oa!oul du mérite du noeud pour insertion dans la file des feuilles et on revi ndè l'étape 203 de calcul pour un autre noeud fils.

Si le test 210 permet de déterminer que toutes les variables discrètes sont inetamoiéem, on peut passer à une étape 214 de mise à jour éventuelle de la meilleure solution et on revient à l'étape 203 de calcul pour un autre noeud fils, sans calcul de mérite ni sous-arbre.

A titre de variante, si le test 210 permet de déterminer que toutes les variables discrètes sont inatanciéee, on peut d'abord peSs rà une étape 213 de mise en oeuvre d'un solveur non linéaire qui permet d'effectuer une optimisation non linéaire basée par exemple sur une méthode des points intérieurs.

Après l'étape 213'Uya passage à l'étape 214 précédemment décrite.

L'exemple de la Figure 11, sans les étapes 208, 209 et 213 est à nouveau explicité ci-dessous.

On part d'une liste triée de noeuds à explorer (étape 201). Le tri est effectué selon un mérite calculé pour chaque noeud. On peut par exemple effectuer une exploration selon la méthode du meilleur en premier ( Best FirsX"\ qui a été mentionnée plus haut. Dans ce cas, un noeud est exploré en priorité lorsqu'il présente la plus faible borne min de la fonction objectif.

Plusieurs fois au cours du procédë, on effectue un test d'élagage (étapes 205` 207). Si le noeud ne peut pas améliorer !a solution courante, il ne sera pas exploré plus avant.

Le principe dnn`oédéB partager un noeud ei 1s fils (étape ~quin , borne eure e~ !a 30' [IUT., deux /ra~~|esO) [ - .?, aloe ces cieux nouveaux noeuis ans une liste de noeuds fils du noeud courant. Puis pour cholua noeud fils /é~ape . on évalue la fonction objectif, c'est-à-dire qu'on évalue les bornes de la fonction objectif à partir des intervalles des variables de ce noeud (étape 204). L'algorithme résultant peut être par exemple le suivant : Tant que la liste L des nœuds à explorer n'est pas vide NoeudCourant e L.PremierElement ; eoeudCourant.TestElagage false nœud courant_ peut contenir r liteen NoeudCourant.Separe ; / - _ : . , a l'intelvalle eelon -ne loi de erati(e Pour e 0 à NoeudCourant.LesteNoeudeflils.taille //peur chae_ue = NoeudFils = NoeudCourant.ListeNoeudeFila[i) ; NoeudFils.Evalue8ornes ; //évaluation des bernes vin et max ue la fonction objectif NoeudFils.TestElagage faire Ras = NoeudFils.Propage ; /:propagation Res 1= 0 fila propagation ne conduit pas a des intervalles vides NoeudFils.EvalueBornes ; //évaluation des bornes min max de la fonction objectif Si NoeudFila.TestElagage e false Si NoeudFils.Faisable true lion vèrefie que le noeud fils contient au moine une solution faisable TestMiseAJourSolution ; //met; à jour la meilleure solution courante s'il y a lieu Si NeeudFila.lnutancie - Talas /ri] y a encore des variables discrètes non inetenclees NoeudFils.CalculMerite ; L.Insert(NoeudFils) ; Fin Si

pin si Fin si Fin Si Fin Pour FinSi Fin Tant que A titre de variante, un noeud pourrait être séparé en plus de deux noeuds fils {mufti-section, par exemple quadri-section). On indiquera ci-dessous quelques compléments concernant l'étape 208 de résolution du système d'optimalité de Fritz-John qui peut apporter une réponse au problème d'actualisation de la borne max de l'optimum en permettant de statuer sur la faisabilité d'un noeud. Considérons le problème d'optimisation suivant : rnin f {X)

C (X) 0,i 'l ..p ≤ 0,i nrhg let' ira 3oudrp D cc Tucker c il ,,ns d'optirrFli -sente '' produire i système d'équations aegénéré dès lors que certaines conti t' tes sont linéairement dépendantes en la solution. Pour obtenir une approche plus robuste, on utilise les conditions d'optimalité de Fritz-John présentées ci-dessous. 33 Les conditions de Fritz-John énoncent qu'il existe X. et ph,... pnqui vérifient le système d'optimalité suivant : X\+ 3.~C|(X)+ K~C/(X]=0 Remarqules muWplicateumpl peuvent être positifs ou négatifs tandis 5 que les multiplicateurs 2,.i sont exclusivement positifs. Une première différence entre les conditions de KKT et celles de Fritz-John réside dans ce que ces dernières introduisent le multiplicateur de Lagrange Xo 1.

Une deuxième différence concernant toujours les multiplicateurs de Lagrange est que, pour les conditions de Fritz-John, les multiplicateurs et pi et peuvent être 10 initi !isào, respectivement, avec les intervalles [0`1] et [-1,1] tandis que, pour les conditions KKT, les multiplicateurs X, et pi sont initialisés, respectivement, avec les

Les conditions d'optimalité de Fritz-John n'incluent pas, à l'origine, de condition de normalisation. Dans ce o a, on peut remarquer qu'il y a (n + p + q + 1) variables et 15 (n + p + q) équations, donc plus de variables que d'équations. Aussi, peut-on considérer la condition de normalisation suivante :

++~~ph+.+enpiv=1 oùej =[1.1 +go[j 1.q (CN1) où eo est le plus petit nombre tel que, selon la précision machine, 1+eo soit strictement supérieur à 1. ou: iaoaaoun on en ida'va!~e mN ) )0=t.p ~O. C'est un Programme Sous Contraintes en Intervalles (ICSP).

Notons alors :

ft(A,M) = Xo .+X + + + eqga -1 et R.2(A,M) = + 2,p + 1.2 + (.1.2 où A= (u,G-.-gq)T (CN1) s'écrit alors :

Ri(A,M) =0

et (CN2) : R,(A,M) = 0

Pour résoudre le système des conditions d'optimalité de Fritz-John, on pose : t = (X,A,M)T et : R, (t) 2,,0Vf(X 2,VC;(X) g, X) i=i j=1 X.,C(X) (D(t) = où k = 1 ou 2 2LpC(X) CE(X) CMX) Notons tl une boîte de dimension N, où N = n + p + q + 1, contenant t. Soit J la jacobienne de d). Pour i, j = 1..N : 34 Ji} Les j p eitlien is ir ei En utilisa a normalis, réai 1cot sont des ra intervenir les lles. Ainsi ers. tige nul muitiplicatel. iman , lut t L _(grange ,i dans la jacobienne c intervalles et augmer.._ les risque_ c 'obtenir une matrce sir gulière. Une méthode de Newton peut alors soit échouer, soit être peu efficace. Dans ce cas, il faut 35 envisager de découper les intenxa!ea. Mais partager les intervalles des multiplicateurs implique, a priori, énormément de calculs supplémentaires. D'où la préconisation d'utUiseryCN1> et l'ordre des variables de t comme indiqué ci-dessus. D'autant que (CN1) présente un caractère linéaire favorable. 5 En utilisant (CN1), certaines méthodes de Newton ne nécessitent pas d'initialiser un intervalle pour les multiplicateurs de Lagrange. Il peut pourtant être intéressant d'en disposer dans certains cas. En particuUer, on peut avoir besoin d'une estimation des valeurs des multiplicateurs, ce qui est le cas dans le problème de calcul de réseau. Une telle estimation pour un multiplicateur peut être obtenue en retenant le milieu de son 10 intervalle ; il faut donc disposer d'un encadrement. Pour le déterminer, on peut utiliser la méthode suivante : On pose :

1 ù 1 e, en A(X) ' !Vy(X) \;7C(X) ù VCF(X) VC(X) VCM X)

Si on résout : 15 A(X) on obtiendra l'encadrement voulu pour les multiplicateurs de Lagrange.

L'utilisation des conditions d'optimalités de Fritz-John au sein du solveur peut être de deux utilités. La première est qu'elles peuvent réduire davantage l'espace des solutions en complément ou en remplacement de la propagation des contraintes à partir d'un certain niveau de l'arbre de la méthode B&B. La seconde vient du fait que la résolution des conditions d'optimalité de Fritz-John est un opérateur de Newton. Peut alors s'appliquer le théorème de Moore-Nickel qui énonce que si un opérateur de Newton permet de réduire un intervalle de définition d'une variable au moins, alors l'espace des solutions courant contient forcément un optimum. Ainsi la résolution de ces conditions d'opLirnaU1é peut aussi être un critère de mise à jour de la borne max de l'optimum de la fonction objecf''. ! ~!c!or du système t méthode itérative d On se ^ conditions d'optimci-d--- (SL' ~~u. itve effectuée. par 36 A est une matrice m x n de réels ou d'intervalles, X est le vecteur des variables de dimension n, B est un vecteur de dimension m de réels ou d'intervalles.

La méthode Gauas-Sgkjel est une méthode itérative faisant suite à une amélioration de la méthode de Jacobi.

5 Une méthode itérative pour résoudre un système système linéaire tdque(S-consiste à construire une suite de vecteurs Xk qui converge vers la solution X". Dans la pratique, les méthodes itératives sont rarement utilisées pour résoudre les systèmes linéaires de petites dimensions uar, dans ce cas, elles sont généralement plus coûteuses que les méthodes directes. Toutefois, ces méthodes s'avèrent efficaces (en termes de coût) dans 10 les cas où le système linéaire (SL) est de grande dimension et contient un grand nombre de coefficients nuls. On dit alors que la ma1rioeA est creuse; c'est le cas lors d'un calcul de réseau.

La méthode itérative de Jacobi consiste à résoudre lieme équation en fonction de X, pour obtenir : 15 ;≠i On construit le terme Xk à partir des composants de Xx-' Or, .. lors du calcul de a'' les composantes X^/ pour / '< i sont connues. La méthode Gauss-Seidel substitue X^ par pour j <i 20 Dans le problème du calcul de réseau, les éléments de A, X et B sont des intervalles. L'algorithme est donc le suivant : k=O

SE=0

25 ipéraUc les é'éments diagonaux de A necVntenar /=1à Am

A,, et no, , ,2est-à-dire nor, ;e,5 un point, Alors Fin Si Fin Pour // Calcul des composantes de x Tant Que SE ≠ 0 et k < nombre itérations maximal k = k + 1 e = SEO) SE = SE - {SEO i e.ligne Bi û 1 tmp=1x A.N 15 20 L H Test de fin 10 XX = X, n tmp Si XX c Xi Alors // inclusion stricte x i = XX Pour j = 1 à A.N, j Si Ai,i SE Alors SE = SE + Fin Si Fin Pour Fin Si Fin Tant Que méthode LU décompose la matrice A du systèr elon le produit A= LU où L est une matrice Mar inférieure à diagonale unité : 38 Lm ( et Uest une matrice triangulaire supérieure : Wx ù 01,, u= : ! O W' ~./ Le système devient donc : LW.X = B (SL) que l'on peut décomposer en deux systèmes : =B La résolution de (SUI) puis de (SL2) est grandement facilitée par la forme triangulaire de L et U.

La Figure 13 montre un exemple de réseau auquel est applicable le procédé d'optimisation automatique selon l'invention. Ce réseau comprend un ensemble de points d'interconnexion (jonctions ou nœuds) 1.1 à 1.13 qui permettent de relier entre elles des canalisations passives 101 à 112 ou des tronçons de canalisation comportant des ouvrages actifs tels que des vannes de régulation 31, 32, une station de compression 41, une vanne d'isolement 51, des consommations 61 à 65 ou des ressources 21, 22. OeSconduits d8by-pass 31A,32A.41Asont associés auxvannes de régulation 31, 32 et à la station de compression 41.

Claims (12)

REVENDICATIONS

1. Procédé d'optimisation automatique d'un réseau de transport de gaz naturel en régime permanent, le réseau de transport de gaz naturel comprenant à la fois un ensemble d'ouvrages passifs tels que des canalisations (101 à 112) ou des résistances, et un ensemble d'ouvrages actifs comprenant des vannes de régulation (31, 32) , des vannes d'isolement (51), des stations de compression (41) avec chacune au moins un compresseur, des dispositifs de stockage ou d'alimentation (21, 22), des dispositifs de consommation (61 à 65), des éléments (41A) de dérivation des stations de compression (41) et des éléments (31A, 32A) de dérivation des vannes de régulation (31, 32), les ouvrages passifs et les ouvrages actifs étant reliés entre eux par des jonctions (1.1 à 1.13), le procédé d'optimisation comprenant la détermination de valeurs pour des variables continues telles que la pression et le débit du gaz naturel en tout point du réseau de transport, et la détermination de valeurs pour des variables discrètes telles que l'état de démarrage des compresseurs, l'état d'ouverture des stations de compression, l'état d'ouverture des vannes de régulation, l'état des éléments de dérivation des stations de compression, l'état des éléments de dérivation des vannes de régulation, l'orientation des stations de compression et l'orientation des vannes de régulation, caractérisé en ce que l'on choisit comme état initial de l'optimisation des intervalles de valeurs pour les variables continues et des ensembles de valeurs pour les variables discrètes, en ce que l'on explore les possibilités de valeurs pour les variables en construisant au fur et à mesure un arbre avec des branches reliées à des noeuds décrivant les combinaisons de valeurs envisagées en utilisant une technique de séparation des variables, c'est-à-dire de découpage conduisant à la génération de nouveaux noeuds dans l'arbre, et d'évaluation, c'est-à-dire de détermination avec une probabilité forte des branches de l'arbre qui peuvent conduire à des feuilles constituant une solution finale optimisée, de manière à parcourir en priorité ces branches à plus forte probabilité de réussite, les valeurs des grandeurs recherchées étant considérées comme optimales lorsque des contraintes prédéterminées ne sont plus transgressées ou sont transgressées au minimum et une fonction objectif est minimisée, cette fonction objectif étant de la forme g = a x Régime + x Energie + y x Cible avec : a, 5 et y sont des coefficients de pondération Régime représente un facteur de minimisation ou de maximisation de la pression en des points donnés du réseau tels que tout point aval d'un dispositif de stockage ou d'alimentation, tout point amont et tout point aval d'une station de compression ou d'une vanne de régulation, et tout point amont d'un dispositif de consommation, Energie représente un facteur de minimisation de la consommation d'énergie de 40 compression, 39Cible représente un facteur de maximisation ou de minimisation du débit d'un tronçon du réseau situé entre deux jonctions ou de la pression d'une jonction particulière, et lesdites contraintes prédéterminées comprenant d'une part des contraintes d'égalité comprenant la loi de perte de charge dans les canalisations et la loi des noeuds régissant le calcul des réseaux, et d'autre part des contraintes d'inégalités comprenant des contraintes de débit minimal et maximal, des contraintes de pression minimale et maximale des ouvrages actifs ou passifs, des contraintes de puissance de compression des stations de compression.

2. Procédé selon la revendication 1, caractérisé en ce que le problème de configuration optimale des ouvrages actifs est modélisé sous la forme d'un programme P, d'optimisation se présentant sous la forme suivante : min1X,s,el f(x,$) = g(x)+ax(fs'12 C,(x)+(3.e <s, CE(x) = SE xE R ,s, E RP,sE E Rq,ee {0,1r avec : x est l'ensemble des variables des débits Q et des pressions P, g(x) est la fonction objectif constituant le critère économique d'optimisation, Ci(x) est l'ensemble des p contraintes d'inégalité linéaires et non linéaires sur les ouvrages actifs, 13 est un vecteur dont les coefficients sont nuls ou égaux aux valeurs maximales des contraintes, e est le vecteur des variables binaires, CE(x) est l'ensemble des q contraintes d'égalité linéaires ou non linéaires, s est une variable d'écart qui, lorsqu'elle est non nulle, représente la transgression d'une contrainte, a est un coefficient représentant le degré d'autorisation de transgression de contraintes.

3. Procédé selon la revendication 1 ou la revendication 2, caractérisé en ce que les variables sont représentées par des intervalles, en ce que la technique de séparation des variables est appliquée aux seules variables discrètes et en ce que des bornes de la fonction objectif sont calculées en utilisant l'arithmétique des intervalles.

4. Procédé selon la revendication 1 ou la revendication 2, caractérisé en ce que les variables sont représentées par des intervalles, en ce que la technique de séparation des variables est appliquée à la fois aux variables discrètes et aux variables continues, la séparation comprenant le découpage de l'espace de définition des variables continues, l'exploration s'effectuant séparément sur des parties de l'ensemble réalisable et l'intervalle de variation de la fonction objectif étant évalué sur chacune de ces parties.

5. Procédé selon la revendication 4, caractérisé en ce que lors de l'exploration des possibilités de valeurs pour les variables avec une technique de séparation des variables et d'évaluation, on établit d'abord une liste de noeuds à explorer triés selon un critère de mérite M calculé pour chaque noeud, tant que la liste de noeuds à explorer n'est P,pas vide, pour chaque noeud courant, on évalue si ce noeud courant peut contenir une solution, si c'est le cas, on découpe l'intervalle correspondant à la variable considérée selon une loi de séparation pour établir une liste de noeuds fils, pour chaque noeud fils on évalue des bornes minimale et maximale de la fonction objectif et on évalue si le noeud fils peut améliorer la situation courante, si c'est le cas on effectue une propagation de la contrainte sur ses variables, si la propagation ne conduit pas à des intervalles vides, on évalue des bornes minimale et maximale de la fonction objectif et on vérifie qu'il n'y a pas d'impossibilité à ce que le noeud fils contienne au moins une solution faisable, on effectue un test pour déterminer s'il y a encore des valeurs discrètes non instanciées, c'est-à-dire des variables pour lesquelles aucune valeur précise et définitive n'a pu être décidée, on met à jour la meilleure solution courante s'il y a lieu et on calcule le mérite du noeud pour l'insérer dans la liste des feuilles, tirée selon ce critère de mérite.

6. Procédé selon la revendication 5, caractérisé en ce que le critère de mérite M est tel qu'un noeud est exploré en priorité lorsqu'il présente la plus faible borne minimale de la fonction objectif.

7. Procédé selon la revendication 5 ou la revendication 6, caractérisé en ce que lors des tests d'élimination des noeuds ne pouvant pas contenir l'optimum, on met en oeuvre l'une des méthodes consistant à utiliser la monotonie de la fonction objectif, à utiliser un test de contraintes transgressées ou à utiliser un test de valeur d'objectif moins bonne que la valeur courante.

8. Procédé selon l'une quelconque des revendications 5 à 7, caractérisé en ce que lors de la séparation d'un noeud courant en noeuds fils, on divise le domaine de variation d'une ou plusieurs variables choisies selon des critères basés sur le diamètre d'intervalles rattachés aux variables.

9. Procédé selon l'une quelconque des revendications 5 à 8, caractérisé en ce qu'il comprend, en outre, un critère d'arrêt basé sur le temps d'exécution ou sur l'évaluation de certains diamètres d'intervalles.

10. Procédé selon l'une quelconque des revendications 5 à 10, caractérisé en ce qu'en complément de la propagation des contraintes, on procède à une actualisation de la borne maximale de l'optimum de la fonction objectif en utilisant les conditions d'optimalité du problème d'optimisation dites de Fritz-John. 5

11. Procédé selon l'une quelconque des revendications 5 à 10, caractérisé en ce que lorsqu'à un noeud du procédé de séparation et d'évaluation, toutes les variables discrètes ont été instanciées, on met en oeuvre en outre un processus d'optimisation non linéaire basé sur une méthode de points intérieurs.

12. Procédé selon l'une quelconque des revendications 5 à 9, caractérisé en ce qu'à chaque noeud du procédé de séparation et d'évaluation, on met en outre en oeuvre un processus d'optimisation non linéaire basé sur une méthode de points intérieurs.