Procédé d'optimisation automatique d'un réseau de transport de gaz naturel

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 (41 A) 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

[0001] 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.

[0002] 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.

[0003] A titre 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).

[0004] 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.

[0005] Pour déterminer par le calcul les caractéristiques d'un réseau de transport de gaz, quelle que soit la modélisation physique retenue, on prend en compte traditionnellement la loi des noeuds et la loi des mailles (encore dénommées lois de Kirchhoff du fait de leur emprunt à la théorie des circuits électriques).

[0006] Un réseau de transport de gaz peut être représenté sous la forme d'un graphe composé de noeuds (sommets) et d'arcs qui établissent une relation orientée entre deux noeuds. Les arcs possèdent un attribut « ETAT » qui indique si l'arc est activé ou désactivé.

[0007] 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.

[0008] De manière synthétique, selon la loi des noeuds, on obtient le système d'équations linéaires suivant :

avec

B : matrice d'incidence du réseau traduisant la correspondance entre les arcs et les noeuds du réseau,

E_arcs : vecteur des quantités s'écoulant en chaque arc,

E_{consommations} : vecteur des quantités livrées aux consommations,

E_ressources : vecteur des quantités émises ou injectées aux ressources (dispositifs de stockage ou d'alimentation),

C_station : vecteur des quantités de gaz carburant consommées par les stations de compression.

[0009] La loi des noeuds permet de définir ainsi un système d'équations linéaires.

[0010] 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.

[0011] Selon la loi des mailles, la somme algébrique, le long d'une maille, des 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 :

avec ΔP : différence des pressions entre deux noeuds consécutifs d'une maille.

[0012] Parce que l'on connaît les formules de perte de charge dans les canalisations selon la forme suivante :

on peut exprimer également de façon équivalente la loi des mailles à l'aide de différences de pression au carré :

avec ΔP² : différence des pressions quadratiques entre deux noeuds consécutifs d'une maille.

[0013] La loi des mailles permet ainsi de définir un système d'équations non linéaires.

[0014] On connaît déjà des procédés de calcul de réseau qui abordent le problème en supposant celui-ci parfaitement déterminé, c'est-à-dire en supposant que le nombre d'inconnues est égal au nombre d'équations.

[0015] 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.

[0016] Les deux lois de Kirchhoff permettent de déterminer des débits (dans la mesure où la loi des mailles remplace les différences de pressions quadratiques par leur expression équivalente fonction du débit, en général de la forme ΔP² = α × Q² où α est considéré constant.

[0017] 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 réseau permet de connaître les pressions en tous les noeuds.

[0018] 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 partout.

[0019] Généralement, le processus itère plusieurs fois entre la phase n° 1 et la phase n°2 car les coefficients α 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.

[0020] Cette démarche impose deux fortes restrictions. La première restriction est 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 expression fonction du débit. En revanche, un ouvrage actif, 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.

[0021] Les ouvrages actifs constituent des organes de commande du réseau en introduisant des inconnues supplémentaires telles que, 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 à ce pourcentage d'ouverture.

[0022] Pour les stations de compression, l'inconnue introduite est la puissance motrice de compression (puissance dépensée pour la compression) puisque cette dernière est reliée au débit et au taux de compression (rapport entre sa pression en aval et sa pression en amont).

[0023] Généralement, les procédés de calcul de réseau permettant la simulation d'ouvrages actifs imposent à l'utilisateur de fixer lui-même la valeur de ces inconnues. Implicitement, les ouvrages actifs ne le sont alors plus car ils présentent une véritable é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 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.

[0024] 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 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.

[0025] 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 après avoir résolu les deux lois de Kirchhoff, les pressions peuvent être connues partout.

[0026] 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 au noeud source une pression la plus basse possible dans la limite où les pressions en tous les noeuds ne soient pas inférieures à un seuil fixé. L'ensemble des pressions en tous les noeuds constitue alors le régime de pression minimal.

[0027] 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 noeud source 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.

[0028] Le réseau est dit saturé.

[0029] Dans le cas d'un réseau ne comportant qu'une seule source, le débit de celle-ci injecté sur le réseau est parfaitement déterminé par la loi des noeuds. Ce n'est plus le cas si une seconde source est présente sur le réseau car le nombre de noeuds n'est pas modifié (donc pas d'équation supplémentaire) et une inconnue additionnelle correspondant au débit de cette seconde source est introduite dans le problème.

[0030] Les procédés traditionnels de calcul de réseau contournent ce cas en créant une maille fictive entre les deux sources. Cette maille est qualifiée de fictive car l'on suppose que les deux sources sont reliées par une canalisation de longueur nulle et de diamètre très grand. L'introduction de cette nouvelle maille pose, dans le système d'équations, l'équation supplémentaire manquante. L'équilibre 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² = C^ste). 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.

[0031] 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.

[0032] 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.

[0033] 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.

[0034] 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.

[0035] 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.

[0036] 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.

[0037] 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.

[0038] 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 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,
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

avec : α, β et γ sont des coefficients de pondération.

[0039] 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,

[0040] Energie représente un facteur de minimisation de la consommation d'énergie de compression,

[0041] 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.

[0042] 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 :

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,

C_l(x) est l'ensemble des p contraintes d'inégalité linéaires et non linéaires sur les ouvrages actifs,

β est une matrice dont les coefficients sont nuls ou égaux aux valeurs maximales des contraintes,

cohérente mais le nombre de variables binaires est réellement : 3 x le nombre d'ouvrages actifs,

C_E(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,

α est un coefficient représentant le degré d'autorisation de transgression de contraintes.

[0043] 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 des intervalles.

[0044] 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 des parties de l'ensemble réalisable et l'intervalle de variation de la fonction objectif étant évalué sur chacune de ces parties.

[0045] 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é 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 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, triée selon ce critère de mérite.

[0046] Le procédé selon l'invention peut en particulier mettre en oeuvre les caractéristiques avantageuses suivantes :

[0047] 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.

[0048] 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.

[0049] 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.

[0050] 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.

[0051] 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.

[0052] 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.

[0053] 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.

[0054] D'autres caractéristiques et avantages de l'invention ressortiront de la description 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 ;
• 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 ;
• la Figure 4 est une vue schématique montrant le processus d'exploration d'un arbre selon la technique de séparation des variables et d'évaluation ;
• la Figure 5 est une vue schématique d'un exemple de partie de réseau à laquelle est appliqué le procédé d'optimisation selon l'invention,
• la Figure 6 est un tableau donnant des exemples d'intervalles de pression d'initialisation pour différents noeuds de la partie de réseau de la figure 5 ;
• 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 ;
• 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 tableau de la Figure 8 où la propagation n'est pas stoppée ;
• la Figure 11 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
• la Figure 13 est une vue schématique d'un exemple de réseau de transport de gaz naturel auquel est applicable l'invention.

[0055] 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 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).

[0056] 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.

[0057] 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.

[0058] 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 directe ou inverse des ouvrages actifs.

[0059] Pour un réseau réel, il existe plusieurs centaines de variables à valeur entière (par exemple 1 pour ouvert et 0 pour fermé) en sus des quelques milliers de variables continues (pressions et débits).

[0060] Le procédé selon l'invention rend possible le lancement de calcul en série, c'est-à-dire sans intervention humaine. Ce caractère autonome du calcul est d'un intérêt majeur dans un contexte de réseaux pouvant donner lieu à une multiplicité de scénarios d'acheminement.

[0061] 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.

[0062] 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.

[0063] Le module 8 constitue par sa part un coeur de calcul permettant d'assurer une optimisation du réseau.

[0064] 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.

[0065] La Figure 2 montre de façon schématique une partie de réseau de transport de 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.

[0066] 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 R qui peuvent être par exemple des réserves de gaz stockées dans des cavités souterraines.

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

[0068] Dans un réseau de transport de gaz, la pression du gaz décroît progressivement 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.

[0069] 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 pouvant être comprise entre environ1 MW et 50 MW.

[0070] La pression au refoulement des compresseurs ne doit pas dépasser la pression maximale de service (PMS) de la canalisation.

[0071] La Figure 3 illustre un exemple de configuration d'une station de compression qui se situe en même temps en un point d'interconnexion 1.0 du réseau. Une première canalisation 100 d'approvisionnement est raccordée au point d'interconnexion 1.0. Une deuxième canalisation d'approvisionnement sur laquelle est placée une vanne de régulation de pression 30 est également raccordée au point d'interconnexion 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.

[0072] 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 canalisation en aval des compresseurs 40.

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

[0074] Dans la suite de la description, on considèrera que l'expression « ouvrage actif » englobe les vannes de régulation et les stations de compression ainsi que les vannes d'isolement, les ressources et les stockages.

[0075] L'expression « ouvrage passif » recouvre les canalisations et les résistances.

[0076] 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 optimiser un critère.

[0077] Le critère est composé des trois termes différents :

• le régime de pression : minimise ou maximise les pressions en aval des stockages et ressources, en amont et aval des stations de compression et des vannes de régulation et en amont des consommations,
• l'énergie : 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.

[0078] 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 (α, β et γ) qui lui donne plus ou moins d'importance :

[0079] Les degrés de liberté sont :

• les pressions en chaque noeud,
• les débits en chaque arc,

pour les variables continues, qui peuvent prendre toutes les valeurs comprises dans un intervalle.

[0080] Les degrés de liberté sont :

• l'ouverture/la fermeture des ouvrages actifs,
• 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 nombre fini de valeurs.

[0081] 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.

[0082] Mathématiquement, ces contraintes sont de deux types : linéaires ou non linéaires.

[0083] 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 0 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 d'une station de compression.

[0084] Exemple pour une station de compression :

Soit x = (Q,P_amont,P_aval) le triplet des variables continues de débits Q et de pressions P_amont et P_aval de la station de compression.

Soient e_f, e_b, e_d, e_i les 4 variables binaires associées aux 4 états alternatifs - fermé, bipassé, direct et inverse - ne pouvant se produire simultanément. Soient C_f (x), C_b(x), C_d(x), C_i(x), les 4 contraintes pour ces 4 états disjonctifs. Par exemple, pour l'état direct, C_d(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.

Soient C_{f max}, C_{b max}, C_{d max}, C_{i max} une estimation des valeurs maximales de ces contraintes quel que soit x. Dans l'exemple de l'état direct, C_{d max} est le vecteur des débits minimal et maximal, des taux de compression minimal et maximal et des puissances minimale et maximale.

[0085] Les contraintes linéaires vont donc s'écrire sous la forme :

• C_f(x) ≤ (1 - e_f).C_{f max},
• C_b(x) ≤ (1 - e_b).C_{b max},
• C_d(x) ≤ (1 - e_d).C_{d max},
• C_i(x) ≤ (1 - e_i).C_{i max},
• e_f + e_b + e_d + e_i = 1 afin d'assurer le choix d'un et d'un seul des 4 états.

[0086] En partant de ce principe, on peut également effectuer une modélisation en ne conservant que les 3 variables e_b, e_d, e_i, ce qui réduit la combinatoire.

[0087] Ces variables vont s'intégrer dans les contraintes de la manière suivante :

• C_f(x) ≤ (e_b + e_d + e_i).C_{f max},
• C_b(x) ≤ (1 - e_b).C_{b max},
• C_d(x) ≤ (1 - e_d).C_{d max},
• C_i(x) ≤ (1 - e_i).C_{i max},
• e_b + e_d + e_i ≤ 1 afin d'assurer le choix entre un des 4 états, l'état fermé correspondant aux 3 variables nulles.

[0088] 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 C_f(x), C_b(x), C_d(x), C_i(x) est non linéaire).

[0089] Le programme général peut donc s'écrire sous la forme suivante :

avec :

x, l'ensemble des variables des débits et de pressions (Q, P),

g(x), une fonction objectif a priori non linéaire. Il s'agit du critère économique (exemple : le coût de fonctionnement des ouvrages actifs tel que le gaz carburant consommé par les stations de compression),

C_i(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,

β, 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 cohérente mais le nombre de variables binaires est réellement : 3 x le nombre d'ouvrages actifs,

C_E(x), l'ensemble des contraintes d'égalité linéaires (exemple : loi des noeuds), et non linéaires (exemple: équations de perte de charge des canalisations). Elles sont au nombre de q.

[0090] Le procédé selon l'invention vise à fournir une réponse quel que soit l'état de 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 cherchera à 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 introduisant les variables s qui, si elles sont non nulles, représentent la transgression des contraintes.

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,

C_i(x) est l'ensemble des p contraintes d'inégalité linéaires et non linéaires sur les ouvrages actifs,

β 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,,

C_E(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.

[0091] Notons que, à variables binaires fixées, le programme P₁, qui n'est pas rigoureusement équivalent à P₀, a une solution proche de P₀ si le coefficient α est choisi suffisamment grand puisqu'alors les variables d'écart s_I et s_E sont recherchées très fortement proches de 0.

[0092] 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.

[0093] Cette mixité du type des variables oblige à faire de l'optimisation combinatoire et continue. C'est pourquoi plusieurs méthodes mathématiques pouvant s'adapter à ces deux types d'optimisation sont de préférence combinées de façon hybride afin d'obtenir au final une solution exacte.

[0094] Le procédé selon l'invention met en oeuvre en premier lieu une technique de séparation des variables et d'évaluation, dite « Branch & Bound » (notée par la suite B&B), c'est-à-dire « ramifier » et « borner ». Cette technique couvre une classe de méthodes d'optimisation capables de traiter des problèmes où interviennent des variables discrètes. Le caractère discret d'une variable s'oppose au caractère continu :

• une variable continue peut prendre n'importe quelle valeur dans un intervalle donné. Dans le cadre du calcul de réseau, ce sera le cas des pressions exprimées en bar, par exemple : P ∈ [40,80],
• une variable discrète ne peut prendre qu'un certain nombre de valeurs. Il s'agit souvent de variables binaires qui représentent par exemple le sens de fonctionnement d'une station de compression par exemple x = 0 (sens direct) ou x = 1 (sens indirect).

[0095] 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 couramment utilisée pour obtenir le minimum global d'un problème d'optimisation mettant en jeu des variables binaires.

[0096] 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 :

• B&B₁: 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&B₂: 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.

[0097] 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,
• l'évaluation des bornes de la solution courante qui permet d'avancer dans la méthode B&B,
• l'élimination des noeuds ne pouvant pas contenir l'optimum (test de contraintes transgressées, de valeur d'objectif moins bonne que la valeur courante, utilisation de la monotonie de la fonction objectif),
• la séparation du noeud courant en noeuds fils (2 ou plus) en divisant le domaine de variation d'une ou plusieurs variables (choisie(s) selon des critères 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),
• le critère d'arrêt basé sur le temps d'exécution ou sur l'évaluation de certains diamètres.

[0098] 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 partielles, celles qui sont susceptibles de mener à la combinaison globale la plus favorable.

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

[0100] On considère un réseau de gaz sur lequel on aurait plusieurs stations de compression. On cherche, par exemple, à minimiser le gaz carburant sur le réseau. Si 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

[0101] 

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.

[0102] 

est la borne maximale de la fonction objectif associée à la meilleure combinaison d'états connus lors de l'exploration du noeud i.

[0103] Si

alors il est sûr que la station 1 orientée en sens indirect

ne conduit pas à la solution optimale.

[0104] En revanche, si

l'exploration continue en fixant une autre variable binaire. Toutes les variables binaires sont ainsi fixées progressivement. Si on n'effectue 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.

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

[0106] 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.

[0107] On ne considère que l'équation non linéaire C(x) et, de manière générale, on cherche à résoudre :

avec :

IR est l'ensemble des intervalles,

X est un vecteur d'intervalles de dimension n.

[0108] 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 correspondant à la valeur de la contrainte (ce qui équivaut à appliquer les règles de l'arithmétique d'intervalle), puis on propage la valeur de l'intervalle de la contrainte, depuis la racine de l'arbre vers les feuilles pour réduire les espaces de définition des variables.

[0109] L'algorithme de propagation d'une contrainte sur ses variables est le suivant :

• Etape 1, propagation : construction de l'arbre de calcul de la contrainte C, les feuilles sont les variables intervalles x_i 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 x_a + x_b,
  le dernier calcul est effectué à la racine.
• Etape 2, rétro-propagation : descente de l'arbre de la racine vers les feuilles.
  A chaque noeud, on tente de réduire le résultat partiel calculé en 1.
  Par exemple : x_a + x_b = [a,b] ⇒ x_a := ([a,b]-x_b) ∩ x_a et x_b := ([a,b]- x_a) ∩ x_b

[0110] La Figure 12 illustre la propagation/rétro-propagation des contraintes pour l'équation suivante donnée à titre d'exemple :

[0111] 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.

[0112] La seconde étape de l'algorithme est explicitée par la partie droite de la Figure 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, uniquement grâce à la propagation.

[0113] L'algorithme de propagation sur l'ensemble des contraintes d'un problème s'effectue de la façon suivante :

1. Initialisation de la file des contraintes à propager

[0114] Pour cela, on insère, sans doublon, dans une file triée selon un critère de mérite M toutes les contraintes.

2. Boucle sur la file des contraintes

[0115]  Tant que la file n'est pas vide { Extraction de la "meilleure" contrainte c (pour le critère M) Propagation de C Si la propagation a conduit à un intervalle vide pour au moins une variable { Sortie de la boucle : il n'y a pas de solution au problème } 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 } } } }

[0116] 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 à une pile FIFO. 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 la file avec un mérite élevé.

[0117] On notera qu'une contrainte est dite résolue lorsqu'elle est 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 acceptables).

[0118] Pour une contrainte C de fonction d'inclusion C(X) = └C(X),C(X)┘, C est résolue si :

• C est une contrainte d'égalité et C(X) = 0,
• C est une contrainte d'inégalité positive et C(X) ≥ 0,
• C est une contrainte d'inégalité négative et C(X) ≤ 0.

[0119] Lorsqu'une contrainte est résolue, sa propagation ne conduira plus à aucune réduction des intervalles de ses variables.

[0120] 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 débit est positif et en sens indirect lorsque le débit est négatif. On peut aussi effectuer une modélisation complète de la configuration des ouvrages actifs en faisant intervenir 3 ou 4 variables binaires comme indiqué plus haut. La mise en oeuvre de la technique de la propagation des contraintes peut être effectuée à l'aide d'une bibliothèque d'arithmétique des intervalles et de propagation des contraintes capable de traiter des variables discrètes.

[0121] Les méthodes de propagation des contraintes peuvent d'une part servir à réduire la combinatoire 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 B&B pour permettre en chaque noeud un meilleur calcul des bornes de la fonction objectif et éventuellement des coupes supplémentaires.

[0122] 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 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 l'arbre, alors on place toutes les contraintes dans la file.

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

[0124] 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 compresseur CP2. Le réseau comprend des noeuds N₀ à N₄ (jonctions ou points d'interconnexion) et des arcs l à VII (canalisations ou tronçons comportant les compresseurs CP1, CP2, la ressource R et la consommation C).

[0125] Le réseau définit cinq variables de pression aux noeuds N₀ à N₄ et sept variables de débit dans les arcs I à VII.

[0126] La Figure 6 donne un exemple d'intervalles de pression d'initialisation (en bar) aux différents noeuds N₀ à N₄.

[0127] La ressource A a une pression de consigne de 40 bar. C'est pourquoi son intervalle d'initialisation est un intervalle de largeur nulle.

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

[0129] La Figure 7 donne un exemple d'intervalles de débit d'initialisation (en m³/h) dans les arcs I à VII.

[0130] La ressource R et la consommation C correspondant aux arcs I et VII ont des débits imposés à 800 000 m³/h. Leurs intervalles sont donc initialisés à des intervalles de largeur nulle.

[0131] Les arcs III et V sur lesquels se trouvent les compresseurs CP1 et CP2 respectivement présentent des intervalles de débit plus réduits que les arcs II, IV et VI correspondant à de simples canalisations.

[0132] Plusieurs tests sont effectués :

A. On teste tout d'abord toutes les combinaisons d'orientation des compresseurs CP1, CP2 (tests A1 à A4).

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 CP1, CP2 (test C).

[0133] Les résultats de ces tests A1 à A4, B1, B2 et C sont présentés sur le tableau de la Figure 8.

[0134] Dans les trois cas où la propagation n'est pas stoppée (tests A1, B1 et C), on obtient les résultats identiques présentés sur les tableaux des Figures 9 et 10.

[0135] La Figure 9 indique les intervalles de pression résultats (en bar) aux différents noeuds N₀ à N₄.

[0136] La Figure 10 indique les intervalles de débit résultats (en m³/h) pour les différents arcs I à VII.

[0137] 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 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.

[0138] 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 caractéristiques principales seront rappelées ci-dessous.

[0139] En arithmétique par intervalles, 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 ε par l'intervalle [x - ε,x + ε]. On peut également 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 encore certifiés.

[0140] Comme cela a été implicitement admis jusqu'à présent, les intervalles qui ne contiennent pas de « trou », sont des sous-ensembles fermés connexes de R. On notera IR l'ensemble des intervalles. On peut les généraliser en plusieurs dimensions : un vecteur intervalle x ∈ IRⁿ est un vecteur dont les n composantes sont des intervalles et une matrice intervalle A ∈ IR^mxn est une matrice dont les composantes sont des intervalles. Une représentation graphique d'un vecteur intervalle de IR, IR² et IR³ correspond respectivement à un segment de droite, à un rectangle et à un parallélépipède. Un vecteur intervalle est donc un hyper-parallélépipède. Par la suite, on utilisera indifféremment les termes de vecteur d'intervalles, de pavé, de boîte ou même d'intervalle.

[0141] On note les objets intervalles par des caractères gras : x. On note x le minimum de x et x son maximum. On a alors x = [x, x ], et on considère l'ordre partiel sur IRⁿ :

[0142] On note w(x) la largeur de x, (avec w pour width) ou encore son diamètre :

[0143] Le centre mid(x) et son rayon rad(x) sont définis par :

[0144] Une fonction F : lRⁿ → IR est une fonction d'inclusion de f sur X ∈ IRⁿ. Si X ∈ X alors f(X) ∈ F(X).

[0145] 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 l'intervalle de diamètre nul.

[0146] 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 :

[0147] De même, le résultat d'une fonction F(z) est le plus petit intervalle contenant l'ensemble :

[0148] Si l'on considère les opérateurs traditionnels +, -, x, ², / ou √, on peut définir les formules suivantes, plus utilisables en pratique que la définition théorique ci-dessus :

[0149] 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 :
  
  
• 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] :
  
  
  
  
• la multiplication n'est pas distributive par rapport à l'addition. Prenons x = [-2,3], y = [1,4] et z = [-2,1] :
  
  
  
  
• La multiplication est en fait sous-distributive par rapport à l'addition, c'est-à-dire que :
  
  

[0150] On peut aussi définir des fonctions élémentaires telles que le sinus, l'exponentielle, etc. prenant en argument des intervalles. Pour cela, on utilise la définition abstraite ci-dessus.

[0151] Si on s'intéresse à une fonction monotone, on déduit facilement les formules permettant de les calculer.

[0152] 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.

[0153] L'arithmétique par intervalles permet de calculer avec des ensembles et d'obtenir des informations générales et précieuses pour l'optimisation globale d'une fonction.

[0154] 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.

[0155] Différentes méthodes de séparation des variables et d'évaluation (B&B) utilisant l'arithmétique des intervalles seront décrites ci-dessous.

[0156] 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,

4. séparation : construction de 2 noeuds fils en divisant le domaine de variation d'une variable,

5. critère d'arrêt.

[0157] Différentes solutions peuvent être choisies pour ces 5 étapes afin d'améliorer la qualité du procédé.

[0158] Soit le problème d'optimisation min_x∈xf(X). Le vecteur d'intervalles de dimension n, X ∈ IRⁿ, est la zone de recherche. La fonction f : Rⁿ→R est la fonction objectif.

[0159] On note f* le minimum global du problème, X* un point optimal tel que f(X*) = f*, et l'ensemble de ces points X* :

[0160] On note les objets intervalles par des caractères gras : x. On note x le minimum de x et x son maximum. On a alors x = [x , x], et on considère l'ordre partiel sur IRⁿ :

[0161] On note w(x) la largeur de x, (avec w pour width) ou encore son diamètre :

[0162] Le centre mid(x) et son rayon rad(x) sont définis par :

[0163] Une fonction F : IRⁿ → IR est une fonction d'inclusion de ƒ sur X ∈ lRⁿ. Si X ∈ X alors ƒ(X) ∈ F(X).

[0164] 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-first » est souvent combinée avec la stratégie « Oldest-first » comme second critère s'il y a des ex-æquo.

1. Plus ancien en premier (« Oldest-first »)

[0165] Cette stratégie consiste à examiner en premier le noeud qui a été créé le plus tôt.

2. Plus profond en premier (« Depth-First »)

[0166] Cette stratégie consiste à examiner en premier le noeud qui se trouve au niveau le plus profond de l'arbre, i.e. le noeud qui a le plus d'ascendants.

3. Meilleur en premier (« Best-First ») [Règle de Moore-Skelboe]

[0167] Cette stratégie consiste à privilégier le noeud qui correspond au plus petit F(X), i.e. celui qui a la plus faible borne inférieure de l'optimum.

4. Index de rejet (« Reject Index »)

a. Optimum connu

[0168] Définissons pour chaque noeud correspondant au vecteur d'intervalles X le paramètre :

[0169] Remarquons que si w(F(X)) est nul, alors il n'y a pas lieu d'évaluer pƒ* puisque le noeud ne sera pas découpé.

[0170] Le noeud sélectionné est alors celui qui correspond à la plus forte valeur de pf*. Le calcul de ce paramètre nécessite cependant de connaître l'optimum à l'avance, ce qui n'est pas toujours le cas. C'est pourquoi des variantes du « Reject Index » s'appuyant sur des estimations de l'optimum ont été développées.

b. Optimum estimé

[0171] La variante du paramètre pƒ* quand on ne connaît pas l'optimum à l'avance peut s'écrire :


où k est l'indice de l'itération considérée. L'indice k correspond globalement au nombre de noeuds examinés et ƒ_k est une approximation de ƒ* à l'itération k.

[0172] Remarquons que la règle « Best First » n'est donc jamais qu'un cas particulier de pƒ pour lequel ƒ_k = ƒ_k. En effet, si Y₀ est l'intervalle du noeud présentant la plus faible borne inférieure de F (« meilleur noeud »), on a alors pƒ(Y₀) = 0 et pƒ négatif pour tous les autres noeuds.

[0173] D'autres possibilités pour f_k peuvent être :


ou encore :

c. Avec contraintes

[0174] Pour un problème contraint de la forme :

[0175] les « Reject Index » définis ci-dessus ne tiennent pas compte du tout des contraintes et risquent de sélectionner des noeuds qui présentent de bonnes valeurs de pf mais conduiront à des noeuds infaisables.

[0176] Certains auteurs proposent donc de définir un index de faisabilité construit de la façon suivante.

[0177] Pour une contrainte C_i et pour un noeud correspondant à un domaine de variation X, définissons :

[0178] Dans le cas où w(C_i(X)) = 0, la faisabilité de la contrainte i peut être décidée directement, et pu_Ci(X) peut être fixé à 1 si X satisfait C_i, -1 sinon. Remarquons que si pu_ci(X) < 0, alors X, de manière certaine, ne satisfait pas C_i car C_i(X) > 0. Au contraire, si pu_Ci(X) = 1, alors C_i(X) ≤ 0 et donc X satisfait C_i de manière certaine. Dans tous les autres cas, l'état de transgression de C_i est indéterminé.

[0179] Pour les X qui ne sont pas « certainement infaisables », c'est-à-dire pour lesquels ∀ i= 1..p, pu_ci(X) ≥ 0, définissons un index de faisabilité global pour l'ensemble des p contraintes :

[0180] Ainsi construit, cet index global possède 2 propriétés :

• pu(X) = 1 ⇔ X est « certainement faisable », i.e. faisable de manière certaine
• pu(X) ∈ [0,1] ⇔ X est indéterminé.

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

[0182] Si pu(X) = 1, i.e. 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 degré de faisabilité de X. Cela permet de définir une nouvelle règle de sélection de noeud : on sélectionne le noeud avec la plus grande valeur de pupf.

[0183] 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) :

avec M une valeur très grande fixée préalablement.

[0184] En effet si pupf(f_k,X) = 0, alors soit pf(f_k,X) = 0, ce qui signifie - dans le cas où f_k =f- que l'on est sûr de ne pas améliorer f; soit pu(f_k,X) = 0 ce qui signifie qu'il existe au moins une contrainte telle que C_i(X) = 0. De tels X ne semblent pas très prometteurs. C'est pourquoi on fixe M à une valeur très grande.

[0185] On considérera maintenant l'étape d'évaluation.

[0186] Dans cette étape, il s'agit d'évaluer les bornes de la fonction objectif, mais aussi celles des contraintes s'il y en a. Pour les méthodes de B&B utilisant l'arithmétique d'intervalles, les fonctions d'inclusion sont généralement obtenues par extension « naturelle » des fonctions usuelles.

Exemple :

[0187] Si f : x → x² - e^x et x = [-5,2], alors F : x → x² - e^x est une fonction d'inctusion de f sur x avec :

et

[0188] Pour l'étape d'élimination, plusieurs méthodes sont possibles.

1. Test de faisabilité

[0189] Si le problème est un problème soumis à p contraintes d'inégalité C_i :

[0190] Soit C_i une fonction d'inclusion de la contrainte C_i. A chaque examen d'un noeud correspondant au domaine de variation de X, les p contraintes C_i(X) sont évaluées. Si 3 i ∈ {1,p}/[-∞,0] ∩ C_i(X) = ∅, alors il est sûr que le noeud ne peut contenir de solution faisable. Il peut donc être élagué.

2. Test de coupure (« Cut Off Test »)

[0191] C'est le critère le plus simple et le plus connu d'élimination : il s'agit de rejeter tous les noeuds pour lesquels f^* ≤ f < F(X), où f est la borne supérieure courante de l'optimum.

3. Test de point milieu (« Middle Point Test »)

[0192] Certaines publications ne font pas de distinction entre le « Cut Off Test » et le « Middle Point Test » (MPT). Le MPT ne serait en fait qu'un moyen supplémentaire de calculer une borne supérieure de f*. Le « Cut Off Test » consiste à prendre initialement comme borne supérieure F(X) puis de l'actualiser à chaque division d'intervalle. Pour un problème contraint, l'actualisation n'est possible que lorsqu'on sait que X contient au moins un point faisable. Dans le MPT, on prend f(mid(X)) qui est aussi une borne supérieure de l'optimum. S'il s'agit d'un problème contraint, il est toutefois nécessaire de s'assurer que mid(X) est un point faisable.

4. Test de monotonie

[0193] Pour un problème non contraint, si la fonction objectif est strictement monotone par rapport à la composante x_i d'un vecteur d'intervalle X, alors l'optimum ne peut se trouver à l'intérieur de x_i. Pour déterminer si f est strictement monotone par rapport aux 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ésultant ne contient pas la valeur 0, alors f est strictement monotone par rapport à x_i.

[0194] Dans ce cas, on peut réduire la composante x_i à un réel : x_i se réduit à x_i si la i^ème composante de la fonction d'inclusion du gradient est un intervalle qui a une borne supérieure strictement négative, et x_i se réduit à x_i si la i^ème composante de la fonction d'inclusion du gradient est un intervalle qui a une borne inférieure strictement positive.

[0195] Pour l'étape de séparation, plusieurs méthodes sont également envisageables :

1. Bissection sur une variable

[0196] 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(max_i=1.nD(i)).

a. Plus grand diamètre

[0197] La fonction de mérite est ici simplement le diamètre de la variable : D(i)=ω(x_i). La difficulté d'utiliser cette fonction de mérite est liée à la nécessité de s'affranchir des 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.

[0198] 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 :

avec mig(X) = min_x∈Xi|X|. On pourrait utiliser la magnitude : mag(X) = max_x∈Xi|X|.

[0199] Cette variante permet ainsi de normaliser le diamètre des intervalles considérés.

b. Règle de Hansen

[0200] Ici,


où ∇F_i est la i^ème composante de la fonction d'inclusion du gradient de f. L'idée est de séparer sur la variable qui a le plus d'impact sur f.

c. Règle de Ratz

[0201] Ici,

[0202] 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

[0203] L'idée sous-jacente est la même, mais on va jusqu'au second ordre :


où H_ik est l'élément de coordonnées (i,k) de la matrice des dérivées secondes (hessien) de f.

[0204] Pour des méthodes qui calculent de toute façon le gradient et le hessien, par différentiation automatique, cette règle n'est pas beaucoup plus coûteuse que les autres.

2. Multi-section

a. Multi-section statique

[0205] Jusqu'ici, nous avons considéré qu'à partir d'un noeud, 2 noeuds fils étaient créés grâce à la bissection du pavé X ∈ IRⁿ dans une direction unique. Cependant, il peut être pertinent de retenir plusieurs directions de séparation. Par exemple, on peut couper l'intervalle de variation de chaque variable en 2, 2ⁿ 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 3² fils, etc.

b. Multi-section adaptative

[0206] 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.

[0207] Une règle hybride (adaptative) utilisera 3 paramètres P₁, P₂ et pf pour déterminer quelle règle utiliser.

[0208] Les paramètres p₁ et p₂ sont deux seuils qui seront à ajuster. pf est le « Reject index » défini plus haut, et est une fonction du noeud considéré.

[0209] Les noeuds qui auront un « Reject Index » pf < p₁ seront séparés suivant la règle (a), ceux tels que p₁ <pf< p₂ seront séparés suivant la règle (b) et ceux tels que pf > p₂ seront séparés selon la règle (c).

[0210] Une telle règle peut tout à fait être définie à partir de variantes de pf, telles que pupf défini plus haut par exemple.

[0211] Différents critères d'arrêt peuvent être utilisés.

1. Diamètre de la zone de recherche

[0212] Un critère d'arrêt peut être l'examen d'un noeud N tel que w(X) ≤ ε 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 objectif

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

3. Temps maximal d'exécution

[0214] Un critère d'arrêt complémentaire peut être un temps maximal d'exécution au-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.

[0215] 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 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.

[0216] 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.

[0217] 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.

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

[0219] 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 procédure à suivre.

[0220] De façon plus particulière, l'étape 201 correspond au choix de la meilleure feuille de l'arbre à explorer. L'étape 202 consiste en une séparation en noeuds fils. L'étape 203 comprend une série d'opérations effectuées pour chaque noeud fils.

[0221] Ainsi, à l'étape 203, il y a d'abord passage à une étape 204 de calcul des bornes, puis un test d'élagage 205 est ensuite effectué. Si la réponse est oui, il y a retour à l'étape 203 pour traiter un autre noeud fils. Si la réponse au test 205 est non, il y a passage à une étape 206 de propagation/rétro-propagation telle que celle proposée par exemple par F. Messine.

[0222] 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'abord à 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).

[0223] Le test 210 permet d'examiner si les variables discrètes sont toutes instanciées ou non.

[0224] Si toutes les variables discrètes ne sont pas toutes instanciées, on passe à une étape 211 de mise à jour éventuelle de la meilleure solution, puis à une étape 212 de calcul du mérite du noeud pour insertion dans la file des feuilles et on revient à l'étape 203 de calcul pour un autre noeud fils.

[0225] Si le test 210 permet de déterminer que toutes les variables discrètes sont instanciées, 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.

[0226] A titre de variante, si le test 210 permet de déterminer que toutes les variables discrètes sont instanciées, on peut d'abord passer à 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.

[0227] Après l'étape 213, il y a passage à l'étape 214 précédemment décrite.

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

[0229] 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 First ») 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.

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

[0231] Le principe du procédé B&B est de partager un noeud en noeuds fils (étape 202). A titre d'exemple, on choisit la loi de séparation suivante : on sépare l'intervalle de la variable du noeud courant qui a le plus grand diamètre (la plus grande différence entre la borne supérieure et la borne inférieure de son intervalle) en deux intervalles. On range alors ces deux nouveaux noeuds dans une liste de noeuds fils du noeud courant. Puis pour chaque noeud fils (étape 203), 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).

[0232] L'algorithme résultant peut être par exemple le suivant

[0233] A titre de variante, un noeud pourrait être séparé en plus de deux noeuds fils (multi-section, par exemple quadri-section).

[0234] 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.

[0235] Considérons le problème d'optimisation suivant :

[0236] L'approche la plus naturelle pour résoudre ce problème d'optimisation est de considérer le système d'équations issu des conditions d'optimalité de Karush-Kuhn-Tucker (KKT). Cependant, ces conditions d'optimalité présentent l'inconvénient de produire un système d'équations dégénéré dès lors que certaines contraintes 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.

[0237] Les conditions de Fritz-John énoncent qu'ii existe λ₀,..., λ_p et µ₁,... µ_q qui vérifient le système d'optimalité suivant :

[0238] Remarquons que les multiplicateurs µ_i peuvent être positifs ou négatifs tandis que les multiplicateurs λ_i sont exclusivement positifs.

[0239] 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 λ₀ ≠ 1.

[0240] Une deuxième différence concernant toujours les multiplicateurs de Lagrange est que, pour les conditions de Fritz-John, les multiplicateurs λ_i et µ_j et peuvent être initialisés, respectivement, avec les intervalles [0,1] et [-1,1] tandis que, pour les conditions KKT, les multiplicateurs λ_i et µ_j sont initialisés, respectivement, avec les intervalles [0,+∝] et [-∝,+∝].

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


où ε₀ est le plus petit nombre tel que, selon la précision machine, 1 + ε₀ soit strictement supérieur à 1.
ou :

[0242] Dans le cas d'un problème d'optimisation en intervalles :

[0243] C'est un Programme Sous Contraintes en Intervalles (ICSP).

[0244] Notons alors :

et

où

[0245] (CN1) s'écrit alors:

et (CN2) :

[0246] Pour résoudre le système des conditions d'optimalité de Fritz-John, on pose :

et :

[0247] Notons t^l une boîte de dimension N, où N = n + p + q + 1, contenant t. Soit J la jacobienne de Φ. Pour i, j = 1..N :

[0248] Les j premiers arguments de J_ij(t,t') sont des intervalles, les suivants sont des réels. En utilisant la normalisation linéaire (CN1), la jacobienne de Φ ne fera intervenir les multiplicateurs de Lagrange que sous forme de réels et non d'intervalles. Ainsi pour résoudre Φ(t) = 0, il n'est nul besoin d'initialiser un intervalle pour les multiplicateurs.

[0249] Utiliser (CN2) implique que les multiplicateurs de Lagrange apparaissent dans la jacobienne en intervalles et augmente les risques d'obtenir une matrice singulière. Une méthode de Newton peut alors soit échouer, soit être peu efficace. Dans ce cas, il faut envisager de découper les intervalles. Mais partager les intervalles des multiplicateurs implique, a priori, énormément de calculs supplémentaires.

[0250] D'où la préconisation d'utiliser (CN1) et l'ordre des variables de t comme indiqué ci-dessus. D'autant que (CN1) présente un caractère linéaire favorable.

[0251] 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 particulier, 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 intervalle ; il faut donc disposer d'un encadrement. Pour le déterminer, on peut utiliser la méthode suivante :

[0252] On pose :

[0253] Si on résout :

on obtiendra l'encadrement voulu pour les multiplicateurs de Lagrange.

[0254] 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'optimalité peut aussi être un critère de mise à jour de la borne max de l'optimum de la fonction objectif.

[0255] La résolution du système linéaire ci-dessus (SL) peut être effectuée, par exemple, par la méthode itérative de Gauss-Seidel (ou de propagation des contraintes) ou par la méthode LU.

[0256] On se donne un système linéaire tel que celui posé par la linéarisation des conditions d'optimalité d'un problème d'optimisation, de la forme :

[0257] A est une matrice m × 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.

[0258] La méthode Gauss-Seidel est une méthode itérative faisant suite à une amélioration de la méthode de Jacobi.

[0259] Une méthode itérative pour résoudre un système linéaire tel que (SL) consiste à construire une suite de vecteurs X_k 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 car, 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 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 matrice A est creuse ; c'est le cas lors d'un calcul de réseau.

[0260] La méthode itérative de Jacobi consiste à résoudre l'i^ème équation en fonction de X_i pour obtenir :

[0261] On construit le terme X^k à partir des composants de X^k-1 :

[0262] Or, lors du calcul de X^k, les composantes

pour j < i sont connues. La méthode Gauss-Seidel substitue

par

pour j < i.

[0263] 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 :

[0264] La méthode LU décompose la matrice A du système (SL) selon le produit suivant :


où L est une matrice triangulaire inférieure à diagonale unité :

et U est une matrice triangulaire supérieure :

[0265] Le système devient donc :

que l'on peut décomposer en deux systèmes :

[0266] La résolution de (SL1) puis de (SL2) est grandement facilitée par la forme triangulaire de L et U.

[0267] La Figure 13 montre un exemple de réseau auquel est applicable le procédé d'optimisation automatique selon l'invention.

[0268] Ce réseau comprend un ensemble de points d'interconnexion (jonctions ou noeuds) 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.

[0269] Des conduits de by-pass 31A, 32A, 41A sont associés aux vannes de régulation 31, 32 et à la station de compression 41.

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

avec :

α, β et γ 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.

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 :

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,

C₁(x) est l'ensemble des p contraintes d'inégalité linéaires et non linéaires sur les ouvrages actifs,

β est un vecteur dont les coefficients sont nuls ou égaux aux valeurs maximales des contraintes,

e est le vecteur des variables binaires,

C_E(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,

α 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 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.

11. Procédé selon i'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.

Dessins