Optimisation sans dérivées pour le
problème de calage d’historique

Un nouvel algorithme d’optimisation sans dérivées adapté aux fonctions partiellement séparables (DFO-PSOF pour “Derivative Free Optimization for Partially Separable Objective Functions”) est présenté dans ce chapitre. Il s’agit d’un algorithme à région de confiance utilisant des modèles d’interpolation quadratique de la fonction objectif. On rappelle qu’une fonction partiellement séparable est une fonction qui peut s’écrire sous la forme : F(x1, . . . , xp) = Xn i=1 fi(xi,1, . . . , xi,pi ) (2.1) avec pour tout i ∈ {1, …, n}, pi < p. L’adaptation à la séparabilité partielle de la fonction objectif se fait au niveau de la construction des modèles, similairement à l’adaptation aux fonctions de type moindres carrés [113, 112] : au lieu de construire un modèle pour la fonction objectif globale, la méthode DFO-PSOF construit un modèle pour chaque sous-fonction objectif. Chacune de ces fonctions ne dépendant que d’un nombre très restreint de paramètres, les modèles d’interpolation peuvent en particulier être construits avec seulement un petit nombre de points sans détériorer leur précision. On présente en section 2.1 un état de l’art des méthodes d’optimisation locales utilisées en ingénierie pétrolière. En section 2.2 sont rappelés les résultats nécessaires à l’introduction du nouvel algorithme. Ce dernier est ensuite décrit en section 2.3 avec la preuve de sa convergence en section 2.4. Enfin, on montre en section 2.5 les premiers résultats numériques obtenus à partir de son implémentation en C++. 

État de l’art

 De nombreux tests ont été réalisés dans la littérature et à IFPEN pour déterminer des méthodes d’optimisation efficaces dans le domaine pétrolier afin d’apporter des solutions “optimales” à des problèmes variés tels que le placement des puits [5, 7, 13], la maximisation de la production [73, 31, 114, 109] ou le calage d’historique [11, 18, 36]. Pour simplifier l’écriture, on considère que ces problèmes s’écrivent sous la forme : Trouver (x ∗ 1 , …, x∗ p ) ∈ [b1, B1] × … × [bp, Bp] tel que (x ∗ 1 , …, x∗ p ) = argminF(x1, . . . , xp) (2.2) où F(x1, . . . , xp) = Pn i=1 fi(x1, . . . , xp) est supposée continûment différentiable avec ∇F continu et Lipschitzien sur le domaine vérifiant les contraintes.Par exemple, le calage d’historique est un problème de minimisation de l’erreur de production sur chaque puits et peut donc s’écrire sous la forme d’une minimisation d’une somme de fonctions. Les contraintes de bornes sont présentes pour assurer que chacun des paramètres ne varie pas au-delà de tout sens physique ou pour éviter d’ajouter des minimums locaux, la fonctions objectif étant périodique en certains paramètres (les paramètres de déformation graduelle par exemple). La fonction objectif étant par exemple périodique pour les paramètres de déformation graduelle, il est naturel de contraindre ces derniers à rester dans la limite d’une période. Il existe une très grande diversité de méthodes d’optimisation dans la littérature dont bon nombre ont été appliquées à des problèmes pétroliers. On peut trouver par exemple des méthodes stochastiques [5, 43] ou déterministes [73, 36], des méthodes globales [13, 48, 79] ou locales [11, 70], etc. Les méthodes globales, bien que permettant généralement un meilleur calage, demandent la plupart du temps un nombre d’évaluations beaucoup plus important que les autres méthodes et sont difficilement utilisables pour des cas non synthétiques de réservoir (dans la plupart des cas, elles sont mêmes trop coûteuses pour les cas synthétiques). On se concentrera ainsi dans ce manuscrit sur les méthodes locales. Plus précisément, on décrira plus en détail les méthodes de descente et les méthodes à région de confiance avant de proposer une nouvelle méthode adaptée aux caractéristiques de la fonction objectif du problème de calage d’historique.

Méthodes de descente 

Les méthodes présentées ici utilisent une information locale sur la fonction objectif qui garantit la plupart du temps, sous des hypothèses de régularité, une convergence de l’algorithme mais n’offre aucune assurance quant au caractère global de l’optimum trouvé. Les hypothèses de régularité peuvent être difficile à vérifier dans le cas du problème de calage d’historique. Par exemple, la condition de différentiabilité de la fonction objectif peut être très contraignante et entraîner des problèmes de convergence dans des cas réels. Les algorithmes décrits dans cette section consistent à donner une estimation initiale de l’optimum puis à la mettre à jour en se déplaçant dans une direction de descente, comme illustré ci-dessous :

1. Choisir une estimation initiale de l’optimum x0 et une tolérance .
2. À l’itération k ≥ 0 : (a) Trouver une direction de descente dk. (b) Mettre à jour l’estimation de l’optimum, xk+1 = xk+αkdk avec un pas αk bien choisi. 
3. Conditions d’arrêt (par exemple) : kxk+1−xkk kx0k < . La direction de descente dk est en général calculée à partir des dérivées au premier et second ordre de la fonction objectif. Un exemple simple est de prendre dk = −∇F(xk), auquel cas dk est appelée la direction de plus forte pente (ou steepest descent). Les méthodes de type quasi Newton s’appuient quant à elles sur une estimation de la Hessienne de la fonction objectif en prenant dk = −H˜ −1 k ∇F(xk). Il est donc important de calculer des approximations fiables des dérivées aux ordres un et deux de la fonction objectif. On présente ci-dessous cinq approches possibles de ce problème, pour la plupart très coûteuses en terme de nombre de simulations.

Approximation du gradient par méthode adjointe

 Une méthode courante en problèmes inverses pour calculer le gradient de la fonction objectif consiste à construire un problème adjoint au problème direct et, à partir de la résolution de ces deux problèmes d’obtenir une évaluation du gradient. Les domaines d’application peuvent être très vastes : Oberail et al. [76] utilisent par exemple ce type de méthode dans le milieu de l’imagerie médicale, Jarny et al. [57] l’utilisent pour résoudre des problèmes de conduction tandis que Fisher et al. [40] en présentent une application dans le domaine de l’optimisation de forme. Une telle méthode s’avère très efficace puisqu’une simple simulation du système et une résolution du problème adjoint permettent d’obtenir une bonne estimation du gradient. Cependant, il n’existe pas de méthode systématique pour construire le problème adjoint ; la construction d’un tel problème peut ainsi se révéler très délicate et nécessite une connaissance approfondie du problème direct. Dans les problèmes inverses issus de l’ingénierie pétrolière, le problème direct correspond à la simulation des écoulements dans le réservoir [78]. Lorsqu’un problème adjoint est disponible, cette méthode offre d’excellents résultats [73, 51, 66, 96]. Cependant, les simulateurs industriels étant très complexes, il est impossible de produire un code adjoint dans la plupart des cas. Il est donc nécessaire de disposer d’une méthode alternative pour calculer le gradient pour le calage d’historique.  

Approximation du gradient par différences finies 

On peut estimer le gradient de F en approximant chacune de ses composantes par différences finies décentrées [75] : ∂F(x) ∂xi ≈ F(x + hiei) − F(x) hi où ei est le i ème vecteur de la base canonique de R p . Cette méthode d’approximation du gradient a notamment été utilisée pour des problèmes pétroliers dans [11, 49, 105]. Elle requiert par contre (p + 1) évaluations de la fonction objectif, ce qui devient très coûteux lorsque la dimension du problème augmente, phénomène mis en évidence dans [5] par exemple. De plus, un problème récurrent avec de telles approximations consiste à choisir correctement le paramètre de pas hi . En effet, le développement de Taylor de F d’où est issue la formule de différences finies nous incite à choisir hi le plus petit possible. Toutefois, dans le cas où une erreur e est commise dans l’évaluation de F (ce qui est le cas dans un problème réel de simulation de réservoir), cette erreur se répercute en e hi dans le calcul de ∂F(x) ∂xi . Une petite erreur peut alors avoir un très grand impact sur l’évaluation de ∇F. Le calibrage du paramètre hi doit alors se faire à partir d’une estimation de l’erreur commise lors de l’évaluation de la fonction objectif, ce qui peut s’avérer très délicat dans la pratique lorsque les problèmes sont bruités.

Gradient stochastique – SPSA 

Une méthode proposée par Spall dans  permet de calculer une approximation du gradient d’une fonction à l’aide de perturbations aléatoires. Le gradient stochastique (ou SPSA) basique s’écrit : ∇F(x) ≈ F(x + δ) − F(x) δ où δ = (δ1, . . . , δp) représente le vecteur de perturbation et  > 0 est l’amplitude de la perturbation. La notation avec δ en dénominateur correspond ici à une division terme à terme de 37 vecteurs, c’est à dire que la coordonnées i du gradient estimé de F est : (∇F(x))i ≈ F(x + δ) − F(x) δi La i ème coordonnée de la perturbation δ est choisie dans la méthode de Spall en faisant un tirage aléatoire à partir d’une loi de Bernouilli symétrique ±1, mais d’autres résultats [65] ont montré qu’il peut être avantageux de choisir une autre loi. A noter que l’écriture avec δi en dénominateur interdit de choisir une loi qui donne la possibilité de tirer une des coordonnées de δ à 0, l’une des hypothèses nécessaires étant même de choisir 1 δi uniformément borné. Deux propriétés importantes sont par ailleurs à signaler pour cette estimation du gradient : — elle propose toujours une direction de descente, en l’occurrence la direction δ, — son espérance est égale au gradient réel de la fonction objectif. la première propriété permet d’assurer la convergence d’un algorithme de descente utilisant ce gradient. La seconde permet d’affiner grandement l’approximation du gradient en moyennant plusieurs estimations calculées à partir de différents tirages de δ. Le principal intérêt de cette méthode est le faible coût pour avoir une estimation du gradient puisque deux évaluations suffisent à l’obtenir. Elle a d’ailleurs été utilisée dans pour des problèmes de calage d’historique, dans pour des problèmes de maximisation de la production et dans  pour l’optimisation du placement de puits. Malheureusement, le manque de précision dans ce calcul ralentit souvent beaucoup la vitesse de convergence des algorithmes d’optimisation, particulièrement dans le das où les problèmes traités sont de surcroit bruités.