Dispositif De Stockage

Paul Charton
email: paul.charton@univ-lorraine.fr

GESTION OPTIMALE D'UNE FERME ÉOLIENNE COUPL ÉE À UN

Les producteurs d'électricité doivent planifier leur production et s'y tenir à sous peine de devoir payer des pénalités financières. Du fait du caractère aléatoire de leur production les producteurs utilisant des sources d'énergie renouvelable sont particulièrement sujets à ces pénalités. Nous étudions ici l'utilisation d'un dispositif de stockage permettant de stocker et restituer de l'énergie. Dans un premier temps nous décrivons notre modèle de production et de stockage en temps continu. Nous en dérivons un problème de contrôle stochastique optimal. Puis nous construisons un schéma numérique permettant d'approcher le gain du producteur et la stratégie de stockage/déstockage optimale.

Abstract. Electricity producers have to plan their delivery. If their effective delivery does not match the planed one, they have to pay financial penalties. Due to the stochastic nature of their production, renewable energy producers are more exposed to those penalties. We study here the use of a storage device. This device allows to store and release energy. We begin by introducing our production and storage model in continuous time. With this modelization, we obtain an optimal stochastic control problem. Then we describe a numerical method to approximate the maximal expected gain of the producer and the corresponding storage/release strategy.

Introduction

Du fait de la quantité limitée d'énergies fossiles et nucléaires disponibles et de leur coût écologique, la tendance actuelle est au développement des énergies renouvelables. En particulier les énergies solaire et éolienne sont en plein essor. L'inconvénient de ces moyens de production réside dans leur caractère aléatoire. En effet, sur certains marchés, les producteurs d'électricité renouvelable sont soumis aux mêmes règles que les autres producteurs. Ils doivent donc décider la veille (J -1) quelles puissances ils s'engagent à fournir sur les différentes plages horaire (0h-1h, 1h-2h,..., 23h-24h) du jour J. Tout manquement à l'engagement entraîne des pénalités financières. Dans ce type de problème on rencontre principalement deux difficultés :

(1) La puissance fournie, sur chaque plage horaire, par un tel moyen de production n'est pas forcément constante.

(2) La différence entre la prédiction de production au jour J -1 et la production effective du jour J peut s'avérer importante. Pour limiter ces pénalités, une première approche, étudiée dans [START_REF] Hedman | Comparing hedging methods for wind power : Using pumped storage hydro units vs. options purchasing[END_REF], est de mettre en place une stratégie de couverture financière. Cette approche a deux principaux défauts. Premièrement, le marché de l'électricité n'est pas complet, il n'y a donc a priori pas de raison de pouvoir mettre en place une couverture exacte. Deuxièmement, cette méthode permet de se couvrir financièrement contre les pénalités de sous-production et de sur-production mais ne garantit pas l'équilibre du réseau, ce qui peut être dangereux (dommage électrique ou coupure pour certains usagers).

Nous nous intéressons ici à une autre méthode : coupler l'unité de production à un dispositif de stockage. Ce dispositif permet de stocker une partie de l'énergie produite sous une autre forme. On peut par exemple penser à des roues d'inertie ou des turbines à air comprimé. Cette méthode a été étudiée dans le cas discret par exemple dans [START_REF] Castronuovo | On the optimization of the daily operation of a wind-hydro power plant[END_REF], [START_REF] Leonhard | Sustainable electrical energy supply with wind and pumped storage -a realistic long-term strategy or utopia[END_REF] et [START_REF] Schainker | Executive overview : energy storage options for a sustainable energy future[END_REF]. Dans ces travaux les auteurs considèrent l'énergie fournie heure par heure. Lorsque cette quantité est supérieure à ce que le producteur s'était engagé à fournir, alors il paie des pénalités de surproduction proportionnelles à la différence. Et lorsqu'elle est inférieure, il paie des pénalités de sous-production. On peut en fait considérer dans leur modèle que la puissance fournie est constante sur chacune des plages horaires. Cependant, comme l'illustre la figure 1, ce n'est pas le cas en pratique. Dans ce cas le producteur paie à la fois des pénalités de sous-production et des pénalités de sur-production, ce qui n'est pas possible dans les modèles précédents. Un autre élément à prendre en compte est le coût d'utilisation du dispositif de stockage. Ce coût n'est pas forcément une fonction linéaire de l'énergie chargée ou déchargée. Par exemple pour certains types de batterie, le nombre de cycles charge/décharge supportés diminue significativement si les décharges sont trop profondes, c'est-à-dire si le niveau d'énergie stockée est trop bas. Il faut alors changer plus régulièrement les batteries. Puiser une certaine quantité d'énergie dans ces batteries n'aura donc pas le même coût selon leur niveau de charge. Si on ne considère que la variation d'énergie stockée sur la période, il est alors difficile de modéliser ce type de coût. C'est pourquoi nous proposons un modèle en temps continu. Dans cet article nous considérons une ferme éolienne. Nous modélisons la puissance fournie par la ferme comme étant un processus stochastique à temps continu, solution d'une EDS. Pour cette première approche, nous supposons que les pénalités et les prix de vente sont connus, nous discuterons dans la suite de la signification et de la pertinence de ce choix. À chaque instant le producteur doit choisir si il charge ou décharge le dispositif de stockage. Il cherche alors une stratégie pour maximiser le gain moyen sur la période. Il s'agit donc d'un problème de contrôle stochastique optimal.

Nous utilisons une approche orientée équation aux dérivées partielles (EDP). Comme montré dans [START_REF] Yong | Stochastic Controls : Hamiltonian Systems and HJB Equations[END_REF], notre problème de contrôle se ramène à la résolution d'une équation de Hamilton-Jacobi-Bellman (HJB). Il existe des méthodes génériques de différences finies (voir par exemple [START_REF] Bonnans | Consistency of generalized finite difference schemes for the stochastic HJB equation[END_REF]), pour approcher numériquement la solution d'une équation de HJB dans le cas général. Ce genre de résolution nécessite de résoudre un système non-linéaire de grande taille. Pour résoudre ce système les méthode directes sont, si elles existent, trop coûteuses en temps de calcul. Dans [START_REF] Forsyth | Numerical methods for controlled Hamilton-Jacobi-Bellman PDEs in finance[END_REF] les auteurs proposent une méthode itérative pour résoudre ce système. Cette méthode consiste à calculer itérativement et alternativement la solution discrétisée pour un contrôle donné en tout point de la grille, puis le contrôle qui maximise le hamiltonien pour la solution ainsi construite, jusqu'à convergence de l'algorithme. Cette méthode a l'avantage d'être applicable à toute équation de HJB, mais le désavantage d'être lente.

La méthode précédente peut être optimisée pour les problèmes de type bang-bang, c'est-à-dire ceux où l'ensemble des valeurs que peut prendre le contrôle est soit fini, soit peut être réduit à un ensemble fini. C'est notamment le cas lorsque l'équation de HJB considérée est convexe par rapport au choix du contrôle. La recherche du maximum consiste alors à comparer un nombre fini de valeurs ce qui diminue la complexité de l'algorithme. Cependant, pour un certain nombre de problèmes, il existe des méthodes plus efficace qui s'appuient sur les particularités de l'équation.

Notre contrôle n'a d'influence que sur le niveau du stock qui est un processus à variation finie. En ce sens notre problème est donc similaire à celui étudié dans [START_REF] Chen | A semi-Lagrangian approach for natural gas storage valuation and optimal operation[END_REF]. Les auteurs de ces travaux utilisent une approche semi-lagrangienne pour obtenir un nouveau schéma de différence finie. Le schéma numérique obtenu demande de résoudre un problème d'optimisation en chaque point du maillage, puis de résoudre un système linéaire. La résolution du problème d'optimisation est équivalent en terme de temps de calcul à la recherche de maximum de l'algorithme précédent. Cependant, il n'y a plus d'itération. La complexité numérique est donc beaucoup moins importante. De plus, dans le cas ou le problème de contrôle est de type bang-bang, la complexité de ce schéma numérique est du même ordre que celle d'un schéma classique de différences finies. C'est cette méthode que nous avons adaptée pour résoudre notre problème de contrôle.

Nous commenons par présenter notre modèle pour la production et le stockage. Ce modèle respecte la capacité de production maximale de la ferme, et la capacité de stockage maximale du dispositif. Les processus obtenus sont alors positifs et bornés de manière déterministe. Ceci nous permet de formuler ensuite un problème de contrôle stochastique optimal. Dans ces conditions, notre problème est équivalent à la résolution d'une équation de HJB modifiée. Cette équation est définie sur l'ensemble des points respectant les bornes du stock et de la production. En particulier cette équation est vérifiée sur les bords de ce domaine. Ceci évite d'avoir à introduire des conditions au bord artificielles. Puis, à l'aide d'un raisonnement heuristique nous introduisons notre schéma numérique de résolution adapté de [START_REF] Chen | A semi-Lagrangian approach for natural gas storage valuation and optimal operation[END_REF]. Enfin, nous présentons quelques résultats numériques obtenus.

Le modèle et le problème de contrôle stochastique optimal associé

Fixation des prix et des pénalités

Une journée J est divisée en 24 périodes de production d'une heure [0, 1], [START_REF] Barles | Convergence of approximation schemes for fully nonlinear second order equations[END_REF][START_REF] Bonnans | Consistency of generalized finite difference schemes for the stochastic HJB equation[END_REF], ... Un producteur d'énergie désirant vendre de l'électricité lors de la période [i-1, i] doit choisir la puissance constante (W E i ) qu'il s'engage à délivrer durant cette période. Il doit ensuite transmettre ses différents engagements de production à l'organisme régulateur du réseau avant l'instant t s le jour J -1. L'organisme régulateur du réseau communique ensuite à l'instant t prix le jour J -1, les prix P i e(kWh) -1 auxquels l'énergie délivrée dans les limites de l'engagement est achetée sur les différentes périodes.

Lors de la journée J, le producteur délivre au réseau à chaque instant une puissance de W del t kW. Tous les producteurs ne tiennent pas forcément leurs engagements ce qui induit un déséquilibre énergétique. L'organisme régulateur utilise alors diverses méthodes pour assurer l'équilibre énergétique, ce qui a un coût. Ce coût est répercuté sous forme de pénalités financières payées par les producteurs n'ayant pas tenu leurs engagements. Ces pénalités sont proportionnelles à la contribution au déséquilibre. Elles dépendent de la période ainsi que des manques ou excès de livraison par rapport à l'engagement pris. Ces pénalités sont communiquées après la période de production à l'instant t pen . On peut alors attribuer à chacune de ces périodes respectivement C + i et C - i représentant le prix de l'énergie sur-produite est sous-produite. Le gain effectué par le producteur sur la journée se calcule de la manière suivante :

24 i=1 i i-1 min(W del t , W E i )P i -C + i (W del t -W E i ) + -C - i (W E i -W del t ) + dt , (1)
où pour tout réel x, x + est la partie positive de x définie par x + = max(x, 0). Dans la suite on se placera du point de vue d'un producteur d'électricité éolienne. On supposera que la capacité de production d'électricité éolienne est faible par rapport à l'énergie nécessaire au réseau. On pourra donc supposer que les pénalités sont indépendantes de la production de la ferme éolienne. Cette hypothèse permet une première approche simplifiée de ce problème. Il serait intéressant d'introduire une dépendance dans les travaux ultérieurs. Nous cherchons à maximiser l'espérance de gain, on peut donc supposer que les pénalités sont constantes et égales à leur moyenne. On se placera après l'instant t prix , les prix et engagements seront alors connus. La figure 2 récapitule la chronologie des différents événements.

temps (t) 0h 1h 2h 23h 24h J J -1 t s W E 1 W E 2 W E 24 Fixés à t s t prix P 1 P 2 P 24 Fixés à t prix t pen C + 1 C + 2 C + 24 C - 1 C - 2 C - 24 Fixés à t pen

Production

Nous allons maintenant considérer un cadre probabiliste. Soit (Ω, F = (F t ) t≥0 , P) un espace de probabilité filtré, dans lequel la filtration F représente l'information disponible au cours du temps.

On note W = (W t ) t∈[0,24] le processus représentant la puissance fournie sur l'intervalle [0, 24] et on le définit comme étant solution de l'EDS suivante :

W t = W 0 + t 0 σ(s, W s )dB s + t 0 b(s, W s )ds, (2)
où W 0 est une variable aléatoire L 2 F 0 -mesurable, (B t ) t≥0 est un mouvement brownien pour la filtration

F et σ, b : [0, T ] × R + → R sont deux fonctions lipschitziennes.
La ferme éolienne a une capacité de production maximale de M W kW. Le processus W représentant la production, doit donc être à valeurs dans [0, M W ]. Pour s'en assurer on fait les hypothèses suivantes :

∀t ∈ [0, 24], b(t, M W ) ≤ 0, σ(t, M W ) = 0 et P (W 0 ≤ M W ) = 1, (H 1 ) ∀t ∈ [0, 24], b(t, 0) ≥ 0, σ(t, 0) = 0 et P (W 0 ≥ 0) = 1. (H 2 )
Nous avons alors le résultat suivant :

Proposition 1.1.

(1) Si l'hypothèse (H 1 ) est vérifiée alors P (∀t ∈ [0, 24]; W t ≤ M W ) = 1.

(2) Si l'hypothèse (H 2 ) est vérifiée alors P (∀t ∈ [0, 24]; W t ≥ 0) = 1.

Stockage et contrôle

On associe à la ferme éolienne un dispositif permettant de stocker et de restituer de l'énergie. Ce dispositif a une capacité de stockage maximale de M Q kWh. À chaque instant le producteur doit donc choisir s'il emmagasine de l'énergie, restitue de l'énergie ou n'utilise pas son dispositif. Nous représentons cela de la manière suivante. Soient f stock , f del : [0, M W ] × [0, M Q ] × [0, 1] → R deux fonctions lipschitziennes telles que :

f del (w, q, 0) = 0, ∀(w, q) ∈ [0, M W ] × [0, M Q ], (H 3 ) ∀(w, q) ∈ [0, M W ] × [0, M Q ], la fonction u → f stock (w, q, u) est croissante, (H 4 ) ∀(w, q) ∈ [0, M W ] × [0, M Q ], la fonction u → f del (w, q, u) est décroissante, ( H 5 )
f stock (w, M Q , u) = 0, ∀(w, u) ∈ [0, M W ] × [0, 1], ( H 6 )
f del (w, 0, u) = 0, ∀(w, u) ∈ [0, M W ] × [-1, 0], ( H 7 )
f del (w, q, u)

+ w ≥ 0, ∀(w, q, u) ∈ [0, M W ] × [0, M Q ] × [0, 1], ( H 8 )
f stock (w, q, u) ≤ -f del (w, q, u), ∀(w, q, u)

∈ [0, M W ] × [0, M Q ] × [-1, 1]. (H 9 )
Supposons qu'à l'instant t ∈ [0, T ] Q t kWh sont stockés, on choisit en fonction des informations disponibles

u t ∈ [-1, 1]. La variation d'énergie du stock à l'instant t est alors de f stock (W t , Q t , u t )dt kWh et on délivre W del t = W t + f del (W t , Q t , u t
) kW au réseau. D'après les hypothèses précédentes, le choix de u t s'interprète de la manière suivante. Si u t = 0 on n'utilise pas le dispositif et toute l'énergie produite est délivrée sur le réseau (H 3 ). Si u t > 0 alors on stocke une partie de la production (H 4 ) on délivre moins d'énergie au réseau (H 5 ). Plus u t est proche de 1 plus on stocke vite et moins on délivre d'énergie au réseau. Si u t < 0 alors cette fois on puise dans le stock pour délivrer une puissance supérieure.

Les autres hypothèses servent à vérifier des propriétés physiques. Si le stock est vide alors on ne peut plus y puiser d'énergie (H 7 ) et s'il est plein on ne peut plus stocker (H 6 ). On ne peut stocker plus vite que l'on ne produit (H 8 ) et il ne peut pas y avoir de création d'énergie lors du stockage ou de la restitution (H 9 ).

Dans la suite nous cherchons à maximiser le gain en utilisant au mieux ce dispositif, en gardant à l'esprit que les choix sont fait uniquement en fonction des informations disponibles et donc doivent être adaptés à la filtration F. On définit les stratégies de la manière suivante : Définition 1.2. Soit A l'ensemble des processus F-adaptés à valeurs dans [-1, 1], on appelle A l'ensemble des stratégies admissibles.

À toute stratégie admissible u et à tout niveau initial de stock Q 0 , on associe le processus Q u solution de :

Q u t = Q 0 + t 0 f stock (W s , Q u s , U s )ds. ( 3
)
Le processus Q u représente l'évolution du stock au cours du temps, associé à la stratégie u. On représente alors la puissance délivrée par le processus W del,u = (W del,u t ) t≥0 défini par :

W del,u t = W t + f del (W t , Q u t , u t ). ( 4
)
Le processus Q u est défini de manière unique et reste compris entre les bornes du stock. Comme nous l'avons précisé dans l'introduction, l'utilisation de ce dispositif peut avoir un coût. Nous représentons ce coût à l'aide d'une fonction c : [0,

M W ] × [0, M Q ] × [-1, 1] → R lipschitzienne. A chaque instant t, l'utilisation du dispositif a un coût de c(W t , Q u t , u t ) eh -1 .

Gain

Dans notre modèle et en particulier d'après (1) et (4) le gain de la journée, associé à la stratégie u ∈ A, s'exprime de la manière suivante

G u = 24 i=1 i i-1 min(W del,u t , W E i )P i -C + i (W del,u t -W E i ) + -C - i (W E i -W del,u t ) + -c(W t , Q u t , u t ) dt := 24 0 g(t, W t , Q u t , u t )dt,
où g est définie de manière évidente sur chacun des intervalles [i -1, i].

On cherche à maximiser l'espérance de gain du producteur dans son ensemble et pas uniquement sur cette période de production. Or la stratégie choisie détermine le niveau du stock à la fin de la journée et ce niveau influe directement sur le gain de la journée suivante. Par exemple, un niveau de stock trop bas est désavantageux si la production est inférieure aux engagements. De même, si le niveau de stock est trop important, on pourra difficilement absorber la surproduction. Comme nous sommes dans un modèle markovien nous allons représenter les gains sur les périodes à l'aide d'un gain moyen potentiel dépendant uniquement de la production et du niveau du stock à la fin de la période (t = 24). Soit h

: [0, M W ] × [0, M Q ] → R une fonction lipschitzienne, h(W 24 , Q u 24
) représente le gain moyen potentiel sur les périodes suivantes. C'est-à-dire le gain moyen lorsqu'on utilise au mieux le dispositif de stockage.

Nous obtenons alors un problème de contrôle stochastique optimal. Un des outils classique de ce domaine est l'équation de HJB. Pour l'utiliser on introduit la fonction v :

[0, 24] × [0, M W ] × [0, M Q ] → R définie par : v(t, w, q) = sup u∈A E 24 t g(W s , Q u s , u s )ds + h(W 24 , Q u 24 )|W t = w, Q u t = q .
Cette fonction associe à un instant t, une puissance produite w et un niveau de stock q, le gain moyen réalisé à partir de l'instant t lorsque le dispositif de stockage est utilisé au mieux. L'objectif devient alors de trouver une stratégie u ε ∈ A telle que :

v(0, W 0 , Q 0 ) = E 24 t g(W s , Q u ε s , u ε s )ds + h(W 24 , Q u ε 24 ) + ε,
où ε est un réel positif fixé. Idéalement, il faudrait trouver u 0 , mais ce n'est pas toujours possible car le supremum d'un ensemble n'est pas toujours atteint. On cherche alors a être au plus près de ce supremum en fixant ε petit. On parle dans ce cas de stratégie ε-optimale. Nous allons dans la suite expliquer comment à la fois approcher la fonction v et construire une stratégie approchant la stratégie optimale par une approche EDP.

Équation de HJB et approximation numérique

Équation de HJB

Comme nous le démontrons dans dans [START_REF] Charton | Optimal Operation of a Wind Farm equipped with a Storage Unit[END_REF], la fonction v est solution au sens de viscosité de l'équation de HJB suivante :

   ∂ t v + Lv + sup u∈[-1,1] {f stock (w, q, u)∂ q v + g(t, w, q, u)} = 0, ∀(t, w, q) ∈ [0, 24[×[0, M W ] × [0, M Q ], v(24, w, q) = h(w, q), ∀(w, q) ∈ [0, M W ] × [0, M Q ],
(HJB) où L est le générateur infinitésimal de (W t ) t∈[0,24] . Cette équation de HJB est particulière au niveau de ses conditions au bord. En effet on fixe uniquement la valeur terminale. Comme le processus (W t , Q u t ) t∈[0,24] est à valeurs dans [0, M W ] × [0, M Q ], les points de la frontière spatiale se comportent comme des points intérieurs pour le processus (W t , Q u t ) t∈[0,24] , on peut donc y étendre l'équation de HJB. Dans ce même article nous avons montré que l'équation (HJB) vérifie un théorème de comparaison : Théorème 2.1. Si u et u sont respectivement sur-solution et sous-solution de l'équation (HJB) alors :

u ≥ u.
Ce résultat implique l'unicité et la continuité de la solution de l'équation (HJB). Il intervient aussi dans la preuve de la convergence de l'algorithme que nous allons présenter dans la prochaine partie.

Schéma quasi-variationnel et résultats numériques

La non linéarité de l'équation (HJB) n'intervient que sur un terme d'ordre 1 ou terme de transport. En cela notre problème est similaire à celui discuté dans [START_REF] Chen | A semi-Lagrangian approach for natural gas storage valuation and optimal operation[END_REF]. Dans cet article les auteurs obtiennent un schéma numérique par un raisonnement heuristique. Ils utilisent ensuite les arguments développés par Guy Barles dans [START_REF] Barles | Convergence of approximation schemes for fully nonlinear second order equations[END_REF] pour prouver que ce schéma est convergent. C'est cette approche que nous avons décidé d'adapter à notre problème.

Nous allons ici juste donner les grandes lignes permettant d'obtenir un schéma convergent. Le détail du schéma et la preuve de sa convergence sont disponibles dans [START_REF] Charton | Optimal Operation of a Wind Farm equipped with a Storage Unit[END_REF]. Les calculs que nous allons effectuer sont purement formels. Ainsi on supposera par exemple que la fonction v est suffisamment régulière pour que toutes les dérivées utilisées aient un sens. En général ce n'est pas le cas. Cependant les résultats que nous utilisons pour prouver la convergence nécessitent uniquement de pouvoir effectuer ces opérations avec des fonctions suffisamment régulières. D'après l'équation (HJB), pour tout τ 1 < τ 2 et pour toute fonction

Q; [τ 1 , τ 2 ] → [0, M Q ] on a : τ2 τ1 sup u∈[-1,1] {∂ t v(t, w, Q(t)) + L(t, w)v(t, w, Q(t)) + f stock (w, Q(t), u)∂ q v(t, w, Q(t)) + g(w, Q(t), u)}dt = 0.
En étudiant le cas des stratégies constantes par morceaux et en utilisant la densité de ces stratégies dans l'ensemble des stratégies admissibles on peut montrer que l'intégrale et le sup commutent :

sup u∈U τ2 τ1 ∂ t v(t, w, Q(t)) + L(t, w)v(t, w, Q(t)) + f stock (w, Q(t), u t )∂ q v(t, w, Q(t)) + g(w, Q(t), u t )dt = 0,
où U est un certain ensemble de fonctions à valeur dans [-1, 1]. Alors en particulier pour la fonction Q u solution de :

Q u (t) = q + t τ1 f stock (w, Q u (s), u s )ds, (5)
on obtient :

sup u∈U τ2 τ1 ∂ t v(t, w, Q u (t)) + L(t, w)v(t, w, Q u (t)) + ∂ t Q u (t)∂ q v(t, w, Q u (t)) + g(w, Q u (t), u t )dt = 0,
cela revient à dire qu'à tout instant on choisit la stratégie réalisant le supremum. En utilisant une intégration par parties il suit :

v(τ 1 , w, q) = sup u∈U v(τ 2 , w, Q u (τ 2 )) + τ2 τ1
L(t, w)v(t, w, Q u (t)) + g(w, Q u (t), u t )dt .

Rappelons que nous connaissons la valeur terminale de v, on cherche donc à résoudre cette équation de manière rétrograde. On utilise la méthode des rectangles à gauche pour approcher les diverses intégrales (ce qui correspond à un schéma implicite). On obtient alors :

Q u (τ 2 ) = q + (τ 2τ 1)f stock (w, q, u τ1 ) := Q(u τ1 ), et v(τ 1 , w, q) -(τ 2τ 1 )L(τ 1 , w)v(τ 1 , w, q) = sup

uτ 1 ∈[-1,1]
v(τ 2 , w, Q(u τ1 )) + (τ 2τ 1 )g(w, q, u τ1 ) .

On utilise ensuite les opérateurs de différences finies classiques pour approcher l'opérateur L et on approche v par interpolation linéaire lorsque Q(u τ1 ) n'est pas sur la grille de discrétisation. Notre schéma se résout donc en deux temps. Premièrement, en chaque point de la grille on résout un problème d'optimisation sur des quantités connues (on procède de manière rétrograde donc toutes les approximations au temps τ 2 sont connues). Ensuite, on résout un système linéaire pour obtenir les approximations au temps τ 1 . On garde en mémoire la solution du problème de maximisation car elle donne une approximation de la stratégie optimale au temps τ 1 sur les points de la grille. Les figures 3 et 4 représentent les résultats obtenus par cette méthode à un instant t fixé.

Figure 1 .

 1 Figure 1. Exemple de production sur une plage horaire

Figure 2 .

 2 Figure 2. Chronologie de la fixation des prix de l'électricité et des pénalités

Proposition 1 . 3 .

 13 Si les hypothèses (H 3 )-(H 9 ) sont vérifiées alors, pour tout Q 0 ∈ [0, M Q ], et tout u ∈ A : (i) Il y a existence et unicité de la solution de l'équation (3). (ii) Le processus Q u , solution de (3), est à valeurs dans [0, M Q ].

Figure 3 .Figure 4 .

 34 Figure 3. Approximation de v à un temps t


Yanne K Chembo

Curtis R Menyuk

Spatiotemporal Lugiato-Lefever Formalism for Kerr Comb Generation in Whispering Gallery Mode Resonators

Keywords: numbers: 42.62.Eh, 42.65.Hw, 42.65.Sf, 42.65.Tg

We demonstrate that frequency (Kerr) comb generation in whispering gallery mode resonators can be modeled by a variant of the Lugiato-Lefever equation that includes higher-order dispersion and nonlinearity. This spatio-temporal model allows us to explore pulse formation in which a large number of modes interact cooperatively. Pulse formation is shown to play a critical role in comb generation, and we find conditions under which single pulses (dissipative solitons) and multiple pulses (rolls) form. We show that a broadband comb is the spectral signature of a dissipative soliton, and we also show that these solitons can be obtained by using a weak anomalous dispersion and sub-critical pumping.

The development of frequency combs -equidistant frequency lines from a short-pulse laser -revolutionized the measurement of frequencies [1] and has opened up a host of potential applications in fundamental and applied physics, including the measurement of physical constants, the detection of earth-like planets, chemical sensing, the generation, measurement, and distribution of highly accurate time, and the generation of low-phase-noise microwave radiation [2]. Ti:sapphire lasers were used as the original source of frequency combs, but in the past seven years, alternative fiber laser sources have developed [2]. Recently, Del'Haye, et al. [3] demonstrated that it is possible to use the whispering gallery modes in microresonators in combination with the Kerr effect to generate an equidistant frequency comb that is also referred to as a Kerr comb. Since many applications in both fundamental and applied science would benefit from the small size, simplicity, robustness, and low power consumption of these whispering gallery mode sources, a considerable worldwide effort has gone into understanding and controlling them [4], and there have been several efforts to develop mathematical models of these sources [5][6][7].

A complete modal expansion has been derived to describe the growth of the Kerr combs from noise [5,6]. This model predicts a cascaded growth in which a primary comb is first generated, which then generates a secondary comb, and later higher order combs. This model is in complete agreement with experiments [5]. However, it is difficult and computationally expensive to use this mode expansion beyond the primary comb generation because calculation time grows as the third power of the number of modes [6]. Moreover, it is difficult to * E-mail: yanne.chembo@femto-st.fr study pulse formation in this model, since a large number of modes interact cooperatively.

Pulse formation plays a critical role in comb generation. It is desirable to find a spatio-temporal model that would be analogous to the Haus modelocking equation (a variant of the nonlinear Schrödinger equation) that has been used with great success to study modelocked lasers [8]. Early research had investigated the modulational instability in a synchronously pumped nonlinear dispersive ring cavity, and predicted the formation of stable temporal dissipative structures for both the normal and the anomalous dispersion regime of the fiber [9,10]. In recent work, Matsko, et al. [11] have used a two-timescales approach similar to that of Haus [8], and they have built a spatiotemporal model that successfully predicted comb generation through mode-locking and pulse formation.

In this communication, we show that a variant of the Lugiato-Lefever equation (LLE) [12] in which the evolution is in time and the conjugate (transverse) variable is the resonator's azimuthal angle is the appropriate spatiotemporal model for Kerr comb generation in whispering gallery mode resonators. We demonstrate that the previously-derived modal expansion is equivalent to a variant of the LLE that includes higher-order dispersion and nonlinearity, and the standard LLE closely approximates the modal expansion when the bandwidth is small compared to the carrier frequency and does not approach an octave. We then use this model to study how multiple pulses (rolls) and single pulses (dissipative solitons) form in the cavity.

The LLE can be considered a variant of the nonlinear Schrödinger equation (NLSE) that includes damping, driving, and detuning. It has been the subject of intensive mathematical study [13][14][15]. Unexpected phenomena such has convective instabilities [16] and excitabil-ity [17] are also known to occur in this system. The identification of the LLE as the fundamental equation governing the field evolution in whispering gallery mode resonators allows us to take advantage of this prior work, while at the same time suggesting efficient computational methods for studying the field evolution, including Kerr comb generation.

Our starting point is the modal equations that were derived in ref. [6]. The modal structure of the whispering gallery mode resonators is well-understood. It has been shown that these resonators can sustain several families of eigenmodes that trap light inside the cavity by total internal reflection. We are only interested in the fundamental family of modes that can be characterized by their toroidal structure. Provided that the polarization is fixed, the members of this family are unambiguously defined by an integer wavenumber ℓ that characterizes each member's angular momentum and can be interpreted as the total number of reflections that a photon makes during one round trip in the cavity.

We denote the eigennumber of the pumped mode as ℓ 0 . If we restrict ourselves to the spectral neighborhood of ℓ 0 , the eigenfrequencies can be expanded in a Taylor series, whose first N elements are

ω ℓ = ω ℓ0 + N n=1 ζ n n! (ℓ -ℓ 0 ) n , (1)
where ω ℓ0 is the eigenfrequency at ℓ = ℓ 0 , ζ 1 = dω/dℓ| ℓ=ℓ0 = ∆ω FSR is the free-spectral range of the resonator or the intermodal angular frequency, and

ζ 2 = d 2 ω/dℓ 2 ℓ=ℓ0
is the second-order dispersion coefficient. In effect, we are making a polynomial fit of N -th order to the frequency ω ℓ , which is a discrete function of ℓ. Since there are many thousands of frequencies and the ω ℓ are smoothly varying, we may treat the ω ℓ as a continuous function of ℓ and set the ζ n equal to the Taylor expansion coefficients.

The quantity ζ 2 corresponds to ζ in refs. [5,6]. In the case of a disk resonator with main radius a, we find ζ 1 = c/n 0 a, where c is the velocity of light and n 0 is the index of refraction at ω ℓ0 . This intermodal angular frequency is linked to the round-trip period of a photon through the resonator as T = 2π/ζ 1 . The second-order dispersion ζ 2 denotes the lowest-order deviation from frequency equistance of the modes. When ζ 2 = 0, the eigenfrequencies are equidistant to lowest order and are separated by ζ 1 . The dispersion is normal when ζ 2 < 0 and anomalous when ζ 2 > 0. In fact, ζ 2 is the sum of two contributions -the geometrical dispersion (generally normal) and the material dispersion (which can be normal or anomalous). Explicit expressions are given in refs. [5,6] for a spherical resonator.

We now consider the elements of the modal expansion in a range of ℓ-values in which the expansion of Eq. ( 1) is valid. The slowly varying envelopes A ℓ obey the equa-tions [5,6]

dA ℓ dt = - 1 2 ∆ω ℓ A ℓ + 1 2 ∆ω ℓ F ℓ e i(Ω0-ω ℓ )t δ(ℓ -ℓ 0 ) (2)
-ig 0 ℓm,ℓn,ℓp

A ℓm A * ℓn A ℓp e [i(ω ℓm -ω ℓn +ω ℓp -ω ℓ )t] ×Λ ℓmℓnℓp ℓ δ(ℓ m -ℓ n + ℓ p -ℓ) ,
where δ(x) is the Krönecker delta-function that equals 1 when x = 0 and equals zero otherwise. The mode fields have been normalized so that |A ℓ | 2 corresponds to the photon number in the mode ℓ. The mode bandwidth ∆ω ℓ = ω ℓ /Q 0 is inversely proportional to the loaded quality factor Q 0 and to the phonon lifetime τ ph,ℓ = 1/∆ω ℓ . The four-wave mixing gain is g 0 = n 2 c ω 2 ℓ0 /n 2 0 V 0 , where is Planck's constant, n 2 is the Kerr coefficient at ℓ = ℓ 0 , and V 0 is the effective mode volume. The variation of the coefficient Λ ℓmℓnℓp ℓ is due to variations in the mode overlap and the Kerr coefficient with mode number. The parameter F 0 denotes the amplitude of the external excitation, while Ω 0 is the angular frequency of the pump laser, and it is assumed to be close to ω ℓ0 . Typically, resonant pumping only occurs when |Ω 0 -ω ℓ0 | ∆ω ℓ0 .

When the spectrum is restricted to the neighborhood of ℓ 0 , then Λ ℓmℓnℓp ℓ ≃ 1, but when the spectrum begins to approach an octave of bandwidth, higher-order corrections, analogous to the higher order terms in Eq. ( 1) must be kept.

At lowest order, we find

Λ ℓmℓnℓp ℓ = 1 + η ℓ (ℓ -ℓ 0 ) + η ℓm (ℓ m -ℓ 0 ) + η ℓn (ℓ n -ℓ 0 ) + η ℓp (ℓ p -ℓ 0 ) , (3)
where η ℓ = ∂Λ/∂ℓ is evaluated when ℓ, ℓ m , ℓ n and ℓ p are all equal to ℓ 0 . The coefficients η ℓm , η ℓn , and η ℓp are all analogously defined. From the symmetry of Eq. (3) under interchange of A ℓm and A ℓp , it follows that η ℓm = η ℓp . As is the case for ω ℓ , the quantity Λ ℓmℓnℓp ℓ varies smoothly as a function of its arguments and may be treated as a continuous function, so that the coefficients in Eq. (3) correspond to the terms in a Taylor expansion. Since the mode overlap is a maximum when ℓ = ℓ m = ℓ n = ℓ p = ℓ 0 , we expect that these first-order contributions will be dominated by variations of the Kerr coefficient due to changes in the mode volume and the contribution from the Raman effect. The second-order contributions from the decrease in the mode overlap will also become important as the bandwidth approaches an octave. A detailed investigation of the competing nonlinear contributions and their relative importance remains to be carried out [18].

The spatiotemporal slowly varying envelope of the total field A(θ, t) may now be written

A(θ, t) = ℓ A ℓ (t) exp[i(ω ℓ -ω ℓ0 )t -i(ℓ -ℓ 0 )θ] , (4)
where θ ∈ [-π, π] is the azimuthal angle along the cir- cumference. From Eq. ( 4), it follows that

∂A ∂t = ℓ dA ℓ dt + i(ω ℓ -ω ℓ0 )A ℓ e [i(ω ℓ -ω ℓ 0 )t-i(ℓ-ℓ0)θ] .( 5
)
Equation ( 4) independently yields

i n ∂ n A ∂θ n = ℓ (ℓ -ℓ 0 ) n A ℓ e [i(ω ℓ -ω ℓ 0 )t-i(ℓ-ℓ0)θ] , (6)
so that if we set ∆ω ℓ = ∆ω ℓ0 , the evolution equation ( 5) can be rewritten as

∂A ∂t = - 1 2 ∆ω ℓ0 A -ig 0 |A| 2 A + 1 2 ∆ω ℓ0 F 0 e iσt (7)
+ N n=1 i n+1 ζ n n! ∂ n A ∂θ n +g 0 η ℓ ∂ ∂θ (|A| 2 A) + 2η ℓm |A| 2 ∂A ∂θ -η ℓn A 2 ∂A * ∂θ ,
where σ = Ω 0 -ω ℓ0 is the detuning between the laser and cavity resonance frequencies. Note that while Eq. ( 7) could be used to account for a frequency dependent nonlinearity, it still neglects the frequency dependence of the absorption. It is useful to translate the frequency of the carrier envelope to remove the explicit time-dependence of the driving term by making the transformation A → A exp(iσt). It is also useful to transform the θ-coordinate to remove the group velocity motion by making the transformation θ → θ -ζ 1 t mod [2π]. When we make these transformations, and the spectrum is restricted to the neighborhood of ℓ = ℓ 0 , so that ζ n , n ≥ 3 can be neglected and so that Λ ℓmℓnℓp ℓ ≃ 1, then Eq. ( 7) becomes This equation can be rewritten in the form of the normalized Lugiato-Lefever equation

∂A ∂t = - 1 2 ∆ω ℓ0 A -iσA + 1 2 ∆ω ℓ0 F 0 -ig 0 |A| 2 A -i ζ 2 2 ∂ 2 A ∂θ 2 . ( 8
)
∂ψ ∂τ = -(1 + iα)ψ + i|ψ| 2 ψ -i β 2 ∂ 2 ψ ∂θ 2 + F , (9)
where the field envelope has been rescaled so that ψ = (2g 0 /∆ω ℓ0 ) 1/2 A * and the time has been rescaled so that τ = ∆ω ℓ0 t/2. The dimensionless parameters of this normalized equation are the frequency detuning α = -2σ/∆ω ℓ0 , the dispersion β = -2ζ 2 /∆ω ℓ0 , and the external pump F = (2g 0 /∆ω ℓ0 ) 1/2 F * 0 . This LLE with periodic boundary conditions is the exact counterpart of the modal expansion as long as higher-order dispersion and the variation of ∆ω ℓ and Λ ℓmℓnℓp ℓ in Eq. ( 3) can be neglected. Since the higher-order corrections can be calculated, it is always possible to check the validity of the LLE. As a consequence, all the conclusions that were previously obtained using the modal expansion and confirmed experimentally can also be obtained by solving the LLE and its extensions when the bandwidth becomes large. These two twin-models are however useful in different and complementary ways. On the one hand, the modal expansion must be used to determine threshold phenomena when a small number of modes are involved. On the other hand, the LLE is appropriate to use when many hundreds or thousands of modes interact since it does not refer to the individual modes. In particular, as we will show later, it is useful in the study of Kerr combs or pulse (solitons or Turing rolls) growth and propagation in which a large number of modes interact cooperatively.

We now discuss the meaning and order of magnitude of the variables and parameters that appear in the dimensionless LLE. The detailed analysis of ref. [6] showed that the threshold value for the Kerr comb is given by |A| 2 th = ∆ω ℓ0 /2g 0 , which exactly corresponds to |ψ| 2 th = 1. Hence, we would expect ψ to be on the order of 1 in our theoretical analysis and numerical simulations. Note that the stationary intra-cavity field ψ 0 before comb generation can be obtained by setting all the derivatives to zero in Eq. ( 9). The dimensionless time can be rewritten as t/2τ ph , where we recall that τ ph is the photon lifetime, which is typically a few µs in ultrahigh-Q whispering gallery mode resonators. Using the external field's phase as our reference, the amplitude F is real and positive and F 2 is proportional to the external pump power. The real parameter α equals the ratio of the detuning to half the linewidth, so that we expect to be off resonance when |α| > 1. This parameter is easily tuned in experiments. The parameter β equals the ratio of the walkoff due to second-order dispersion to the half-linewidth. This parameter is negative for anomalous dispersion and positive for normal dispersion. It is interesting to note that the LLE formalism sheds new light on the role of dispersion on Kerr comb formation: in Eq. ( 9), no comb generation can occur when the dispersion parameter β is null, since ψ becomes independent of the azimuthal angle θ ("flat" solutions). The case of the generalized LLE in Eq. ( 7) is however less trivial to analyze in this regard. The overall dispersion has to be a relatively small parameter in all cases. In refs. [5,6] for example, it was found that β ∼ -0.01 was sufficient to generate combs. In practice, the dispersion magnitude |β| should preferably not be too large (∼ 1), as it would mean that the deviation from equidistance of the mode frequencies would be too strong for the nonlinearity to compensate, and the emergence of a wide-span Kerr comb would be suppressed [19,20].

Numerical simulations were performed using the splitstep Fourier algorithm, which is commonly used in simulations of the 1-D generalized NLSE. We note that this simulation method inherently assumes periodic boundary conditions because it is based on the fast Fourier transform (FFT). For studies of pulses in optical fibers -in which the pulse duration can often be many orders of magnitude smaller than the separation between pulses -this feature is often an annoyance that leads to boundary-induced artifacts. In our case, the physical boundary conditions are periodic. The split-step Fourier algorithm is therefore a remarkably cost-effective computational tool, enabling the simulation of Kerr comb dynamics with a laptop computer in a few minutes regardless of its spectral span, as opposed to a few days with the modal expansion for wide-span combs.

For the numerical simulation of the LLE, Eq. ( 9), we have considered a calcium fluoride resonator with a radius a = 2.5 mm. The polar eigennumber of the TEpolarized pump mode is ℓ 0 = 14350, corresponding to a wavelength of 1560.5 nm and ω ℓ0 = 2π × 192.24 THz in vacuum. The index of refraction is n 0 = 1.43. The resonator is critically coupled with a loaded quality factor Q 0 = 3×10 9 , corresponding to a central mode bandwidth of ∆ω ℓ0 = ω ℓ0 /Q ≃ 2π × 64 kHz. The free spectral range is equal to ∆ω FSR = ζ 1 = 2π × 13.36 GHz.

In Fig. 1, we present simulation results that correspond to those that have already been obtained from the modal expansion, namely the multiple-FSR solution. It can be seen that the solution corresponds to the formation of "rolls" in the time domain. We take advantage of the LLE formulation to show particular solutions that would be difficult to observe in a strictly modal study. Figure 2a shows the formation of a single-peaked cavity soliton. This dissipative localized structure is obtained after sub-critical pumping. The figure corresponds to a single pulse of width ∆T ≃ 500 fs, circulating inside the cavity with a round-trip time T = 2π/∆ω FSR ≃ 75 ps. In the frequency domain, this pulse consists of several tens of modelocked whispering gallery modes, whose frequencies have been nonlinearly shifted so that they are equidistant, as shown in Fig. 2b. In order to obtain a single soliton, it is necessary to lower the power so that it is below threshold and to reduce the magnitude of the anomalous dispersion so that it is closer to zero. Because this cavity soliton is subcritical, it emerges abruptly and at a low pump power. The solitonic Kerr combs do not grow as the pump power increases; instead, they are destroyed. By contrast, super-critical combs like the one shown in Fig. 1 have spectral components that grow in number and power as the external pump power grows. However, the dissipative cavity solitons are robust, and their full width at half maximum decreases with the dispersion parameter |β|.

In conclusion, we have demonstrated that the Kerr comb evolution in whispering gallery mode resonators can be modeled using the Lugiato-Lefever equation and its extensions, where the evolution variable is time, and the conjugate (transverse) variable is the azimuthal angle. We have shown that low-threshold, wide-span combs can emerge as dissipative cavity solitons. With different parameters, rolls appear, corresponding to narrowerband combs. So, when the goal is to optimize the bandwidth a correct parameter choice is critical. Additionally, we demonstrated that the LLE can be extended to incorporate higher-order dispersion and nonlinearity, and we expect that it can be extended without too much difficulty to include Rayleigh and Brillouin scattering [START_REF] Boyd | Nonlinear Optics[END_REF]. A key advantage of the spatio-temporal model that we have developed is that the LLE has already been the subject of extensive mathematical study, and it should be possible to take advantage of this earlier work to shed additional light on the dynamical properties of Kerr combs. We expect that a better understanding of Kerr comb generation in whispering gallery resonators will be the result, leading to new resonator designs that can produce octave-spanning combs [START_REF] Chembo | [END_REF]23].

Y. K. C acknowledges financial support from the European Research Council through the project NextPhase (ERC StG 278616), and insightful discussions with M. San Miguel about the LLE. C. R. M acknowledges hospitality and support from the FEMTO-ST Institute within the framework of the LabEx Action, and also acknowledges useful discussions with T. Sylvestre.

FIG. 1 :

 1 FIG. 1: Multiple-FSR comb with α = 0, β = -0.0125 and F = 1.71. The intracavity pump is |ψ0| 2 = 1.2, so that it is above threshold |ψ| 2 th = 1. The parameters correspond to Fig. 10b in ref. [6]. (a) Time-domain dynamics that consists of 19 smooth pulses ("roll" Turing pattern solution). (b) FFT of the roll solution, corresponding to an order 19 multiple-FSR comb.

FIG. 2 :

 2 FIG. 2: Temporal cavity soliton with α = 1.7, β = -0.002 and F = 1.22. The intracavity pump is |ψ0| 2 = 0.95, so that it is below threshold. (a) Cavity soliton formation. (b) FFT of the soliton, corresponding to a wide-span, low-threshold Kerr frequency comb.