Zhihao Cen
email: zhihao.cen@polytechnique.edu

J Frédéric Bonnans
email: frederic.bonnans@inria.fr

Christel Thibault
email: thibault.christel@total.com

Energy contracts management by stochastic programming techniques

Keywords:

e onsider the prolem of optiml mngement of energy onE trtsD with ounds on the lol @time stepA mounts nd glol @whole peE riodA mounts to e trdedD integer onstrint on the deision vriles nd unertinty on pries onlyF efter uilding (nite stte wrkov hin y using vetoril quntiztion tree methodD we rely on the stohsti dul dynmi proE grmming @hhA method to solve the ontinuous relxtion of this stohsti optimiztion prolemF en heuristi for omputing su optiml solutions to the integer optimiztion prolemD sed on the fellmn vlues of the ontinuous relxtionD is providedF gomining the previous tehniquesD we re le to del with highEdimension stte vriles prolemsF xumeril tests pplied to relE isti energy mrkets prolems hve een performedF

Introduction

he physil energy mrkets re rtiulted round long term supply nd deE mnd ontrtsF heir purpose is to ommodte oth produers nd onE sumers ording to their needsF yn one sideD produer will e keen on selling energy on long term sis for sh )ow ertinties llowing him to shedule other investmentsF yn the other sideD it ensures the onsumer of ertinty on volumes nd ts s hedge ginst short term unertintiesF vong term energy ontrts typilly hve the following feturesF pirst of llD they re ssoited with grnulrity de(ning the periods t whih the energy is exhngedF eondlyD the mounts of trded energy re ontrtully upper nd lower ounded during eh period s well s on the overll ontrtF he motivtion ehind these ounds is generlly to llow the onsumer to modE ulte his energy purhseF pinllyD the prie t whih the energy is ought or sold t given period is funtion of the index spot prie whih is stohsti vrileF hereforeD the optiml mngement of long term ontrts implies for the option9s holder to suessively hoose the est possile deision over the quntity to e ought or sold ording to the prie informtion ville t tht momentF ih deision tken will impt the reminder of the ontrtF sn dditionD we might onsider the sitution where the deisions to e mde re disreteF sndeedD unlike pipeline nturl gs where it is generlly possile to ssume tht ny quntity of gs n e exhnged s long s it ful(lls the ontrtul oundriesD this is no longer true for se shipped ommodities suh s vique(ed xturl qs @vxqAF sn ftD the quntity of energy to e lloted is represented y the numer of vxq vessels to llote to destintionF he omintion of stohsti nd disrete deisions mkes the prolem very di0E ultD nd it is nturl to onsider ontinuous relxtion of the prolemD t lest s (rst stepF por surveys on stohsti progrmmingD we refer to SD IVD PWF nless the unE derlying rndom spe hs (nite supportD our model is n optimiztion prolem over in(niteEdimensionl funtion spesD whih is di0ult to solveF enlytil solutions re not ville exept for extremely simple nd unrelisti sesF xumeril pproximtion methods hve een lrgely studied in the litertureF wost pproximtion shemes re sed on disretiztion of the underlying rndom spe y senrio treeF e survey nd evlution of populr senrio genertion tehniques is provided in IWF emong themD let us mention smplingE sed methods PWD ghpter QD moment mthing methods ITD tree redution methods VD IQD reomining tree PHF yne the underlying proility spe eomes disreteD (nite dimensionl optimiztion lgorithms n e pplied without di0ultyF roweverD otining high preision for the optiml vlue nd getting the optiml solution is still ig hllengeF enother min pproximtion method is to projet the ontrol funtion or ostEtoEgo funtion into (nite dimension funtionl spe RD QHD QIF e hve hene to solve (nite dimensionl prolemD whose dimension is the rdinlity of the funtionl speF he dvntge is to void disretizing the underlying proility speF xumeril simultion n e done sed on wonteEgrlo methodsD nd thus the onvergene results re given y lssil wonteEgrlo methods properties s well s rilert spe theoryF his method hs een suessfully pplied PUF roweverD in prtieD the result depends strongly on the sis funtions hosenF purthermoreD when the stte vrile @nd rndom RR n°7289 vrileA dimension inresesD the dimension of sis funtion spe usully grows exponentillyF eprt from these two min pproximtion methodsD few tehniques hve een proposed for ertin spei( prolemsF por exmpleD ereir nd into PQ suggested stohsti dul dynmi progrmming @hhA lgorithm to solve the prolem where rndomness only ppers on rightEhndEside of onstrint nd is independent from one stge to notherF hilpott nd qun PS hve studied the onvergene of this methodF hpiro hs reently nlysed the sttistil properties nd rte of onvergene of this method PVF his method n lso e pplied to ewe proesses IUF pinllyD we note tht wo et lF PP use hh pproh omined with (nite stte wrkov hin representtion of the evolution of the prie @whih is slrAF e spei( feture of our prolem is tht the rndom vriles re funtions of some frownin motionD nd therefore we hve to design n pproxE imtion proedure in order to redue the prolem to one in whih the rndom vriles hve wrkov hin representtionF he min ontriution of this rtile is to use the vetoril quntiztion tree method introdued in flly nd gès PD ID whih tkes dvntge of the ft tht the rndom vrile is qussinF e pply the stohsti dul deomposition method to the ontinuous relxtion of the quntized prolemF he numeril test results present good performnes when the two methods re ominedF his rtile is orgnised s followsX in setion PD we present the model nd its propertiesF sn setion QD we introdue the quntiztion tree methodsD nd then in setion R the lgorithmF pinllyD numeril tests re performed over relisti energy mngement prolems in setion SF e heuristi method for omputing integer suoptiml solutions sed on the fellmn vlues otined y ontinuous relxtion is lso providedF 2 Problem formulations 2.1 Framework e now formulte the lss of prolems of onern in the following settingF pirst of llD we hve disrete time wrkov proess (ξ t ), 0 ≤ t ≤ T in the proility spe L 2 (Ω, (F t ), P; R d )D where F t denotes the noni (ltrtion ssoited with (ξ t ) X F t := σ (ξ s , 0 ≤ s ≤ t)F e my write its dynmis s

ξ t+1 = f (W t , ξ t , α t ) t = 0, . . . , T -1, @IA where ξ 0 is deterministiD (α t ) t∈[0,T -1] is sequene of deterministi prmeE tersD nd (W t ) t∈[0,T -1] is sequene of independent rndom vriles in R d D independent of F t Y (nlly f is mesurle mppingF he deision prolem hs the following expressionX inf (ut) E T -1 t=0 c t (ξ t ) • u t + g(ξ T , x T ) sujet to u t ∈ U t , x t+1 = x t + A t u t ,
x 0 = 0, x T ∈ X T lmost surely;

F Lemma 2.1. If T -1 t=0 A t U t ∩ X T +1 = ∅,
problem @PA has a nonempty set of solutions.

Proof. he set of fesile poliies is losed nd onvexD nd heneD wekly losedF he (nl ost E[g(ξ T , x T )] is lFsFF onvex funtion over L 2 (Ω, F T , P; R m )D nd therefore is wekly lFsFF in this speY so is the ll ost funtion @c t (ξ t ) • u t AF ine the minimizing sequene u k is ounded in L 2 (Ω, (F t ), P; R n )D extrting if neessry susequeneD we my ssume tht it wekly onverges in this speF ssing to the limit in oth the stte eqution nd in the onstrint on the deision nd using the wek lFsFF of the ost funtionD we otin tht the wek limit of u k is solution of @PAF he ondition T -1 t=0 A t u t ∩ X T = ∅ is ssumed to e lwys stis(ed in this rtile in order to mke sure of the existene of solutionF Remark PFP. he funtion to e minimized is not stritly onvexF herefore the minimum is not neessrily uniqueF 2.2 Dynamic programming principle he prolem under study hs the desired struture for the dynmi progrmE ming priniple to hold @see PWD ghpter QAF vet us de(ne the reliztion up to time t of the rndom proess s ξ [t] = (ξ 0 , . . . , ξ t )D nd of the stte vrile vlue s x [t] = (x 0 , . . . , x t )F he fellmn vlue t stge tD denoted y Q(t, x [t] , ξ [t] )D is the vlue of the vrint of the originl prolem otined when strting t time t with stte vrile @storgeA x t X

Q(t, x [t] , ξ [t] ) := inf E T -1 s=t c s (ξ s ) • u s + g(ξ T , x T ) sujet to u s ∈ U s , F s -mesurle, x s+1 = x s + A s u s , x T ∈ X T lmost surely.
@QA sn view of the wrkov property of (ξ t )D (ξ s ) t≤s≤T only depends on ξ t insted of

ξ [t] F end x t is F t-1 mesurle euse u t is F t mesurleD then Q(t, x [t] , ξ [t]
) is tully funtion of (t, x t , ξ t )F he optiml vlue of the originl prolem is Q(0, 0, ξ 0 )F he following dynamic programming principle is esily estlishedX

Q(t, x t , ξ t ) = inf c t (ξ t ) • u t + Q(t + 1, x t+1 , ξ t ) sujet to u t ∈ U t , x t+1 = x t + A t u t ; @RA where Q(•, •, •) is the onditionl expettion of the vlue t time t + 1 otined in stte x t+1 X Q(t + 1, x t+1 , ξ t ) := E Q(t + 1, x t+1 , ξ t+1 )ṽF t = E Q(t + 1, x t+1 , ξ t+1 )ṽx t+1 , ξ t @SA
where the seond equlity omes from the ft tht x t+1 is F t mesurleF pinlly the fellmn vlue funtion t (nl stge isX

Q(T, x T , ξ T ) = g(ξ T , x T ) if x T ∈ X T , +∞ otherwiseF @TA
2.3 Feasibility he property of relatively complete recourse does not hold generlly in this settingF sndeedD it my hppen tht for given time t nd stte vrile x t D where the ssoited fellmn vlue is (nite @implying the existene of sequene of fesile deisions with the orresponding (nl time elongs to X T AD some u t ∈ U t my result in non fesile stte vrile x t+1 := x t + A t u t D whose ssoited fellmn vlue will therefore e in(niteF sn other wordsD the (nl stte onstrint indues impliit onstrints on the deision u t t time tF he hh pproh inludes mehnism for deling with this situtionD y generting fesiility uts @similr to fender9s fesiility uts for two stge prolem ut omined with kwrd reursionAF roweverD generting fesiE ility uts my e long proess nd we prefer to del with them y diret pprohD y tking dvntge of the spei( property tht no unertinty oE urs in the stte dynmisF sn order to gurntee tht the (nl onstrints re ful(lledD let us de(ne the sets of fesile ontrolX

U ad := (u t ) T -1 t=0 A t u t ∈ X T , u t ∈ U t , t = 0, . . . , T -1 . @UA nd U ad = L 2 (Ω, (F t ), P; U ad )F et stge t ∈ [0, T -1]
nd for given stte x t is ssoited the @possily emptyA lol set of fesile ontrol

U ad t (x t ) := u t x t + T -1 s=t A s u s ∈ X T , u s ∈ U s , s = t, . . . , T -1 . @VA nd U ad t = L 2
(Ω, F t , P; U ad t )F reneD we n lso de(ne X t = xṽU ad t (x) = ∅ D suh tht only for x t ∈ X t D there exists fesile ontrol leding to the (nl stte fesile set X T +1 X and Q(t, x t , ξ t-1 ) are convex and lower semi-continuous functions of x t .

X t := x t ∃u s ∈ U s s ≥ t, s.t. x t + T -1 s=t A s u s ∈ X T = X T - T -1 s=t A s U s . @
Proof. e prove the onvexity nd lower semiEontinuous property y kE wrd indutionF sn view of @TAD oth properties hold for funtion Q(T, x T , ξ T )F he onditionl expettion @SA keeps lerly the onvexity on x t sine the onditionl expettion opertion is linerF es to lower semiEontinuous propE erty of Q(t + 1, x t+1 , ξ t )D we pply the ptou lemm on

Q(t + 1, x t+1 , ξ t+1 ) ≤ lim inf x k t+1 →xt+1 Q(t + 1, x k t+1 , ξ t+1 )X Q(t + 1, x t+1 , ξ t ) = E[Q(t + 1, x t+1 , ξ t+1 )ṽF t ] ≤ lim inf x k t+1 →xt+1 E[Q(t + 1, x k t+1 , ξ t+1 )ṽF t ] = lim inf x k t+1 →xt+1 Q(t + 1, x k t+1 , ξ t+1 ).
@IHA he onvexity follows y indution over tD from T -1 to 0D using @RAF he lower semiEontinuous propertyD in view of theorem PTD heorem IFIU on prmetri minimiztionD holds s the expression on right hnd side of @RA is proper @whih mens tht it is lwys greter thn -∞D nd not lwys equl to +∞AD lower semiEontinuous nd levelEounded in u t lolly uniformly in x t PTD he(nition IFITF fy indutionD oth properties follow for ll t = 0, . . . , T F Remark PFR. Q(t, x t , ξ t ) nd Q(t+1, x t+1 , ξ t ) re not in generl onvex funtions of ξ t F e denote y f * the penhel onjugte of n extended rel vlued funtion f on n iuliden spe XD identi(ed with its dulD de(ned for eh x * ∈ X y

f * (x * ) := sup x∈X (x * • x -f (x)) . @IIA vet f * * denote the ionjugte of f D de(ned y f * * (x) := sup x * ∈X * (x • x * -f * (x * )) .
fy the woreu penhel theorem PTD hm IPFID f = f * * whenever f is onE vexD lower semiEontinuous nd proper F epplying this theorem to Q(t, x t , ξ t ) llows to write this funtion s the supremum of fmily of 0ne funtions of xD depending on (t, ξ t ) X

Q(t, x t , ξ t ) = Q * * (t, x t , ξ t ) = sup x * t ∈R m x * t • x t -Q * (t, x * t , ξ t ), @IPA
where Q * denotes the onjugte of Q with respet to x onlyF RR n°7289 husD we n use the stohsti dul deomposition pproh on the stte vrile x t SD PWD ut not on the rndom vrile ξ t F he orresponding formultion of the dynmi progrmming priniple is

Q(t, x t , ξ t ) = inf c t (ξ t ) • u t + Q(t + 1, x t+1 , ξ t ) @IQA = inf c t (ξ t ) • u t + Q * * (t + 1, x t+1 , ξ t ) @IRA sujet to x t+1 = x t + A t u t , (fesiility) u t ∈ U ad t (x t ), (optimlity ut) Q * * (t + 1, x t+1 , ξ t ) ≥ E Q * * (t + 1, x t+1 , ξ t+1 )ṽF t , Q * * (t + 1, x t+1 , ξ t+1 ) ≥ x * • x t+1 -Q * (t + 1, x * , ξ t+1 ), ∀x * .
his is the formultion we will pply in our lgorithmF roweverD in this formultionD we still need to lulte onditionl expetE tionsF wny numeril methods hve een proposedD suh s hi or wonteE grlo methodsF rereD we hoose to disretize the rndom speF gonsidering the size of the portfolio tht we re interested in @see seond exmple in seE tion SAD stndrd senrio trees will use n inrese of the dimension of the prolem suh tht it will low up numerillyF o void this pointD we rely on the vetoril quntiztion tree whih uilds wrkov hin with (nitely mny sttesF 2.5 Heuristic method for integer solution sn this pperD our min motivtion onsists of solving the vxq vessels routing prolemF he deision vriles desriing how mny rgos should e 'eted to route re integersF he ontinuous prolem presented in @PA n e onsidE ered s ontinuous relxtion of the following integer prolem desriing the vxq vessel routing prolemX

inf E T -1 t=0 c t (ξ t ) • u t + g(ξ T , x T ) sujet to u t ∈ L 2 (Ω, (F t ), P; Z n ∩ U t ), x t+1 = x t + A t u t ,
x 0 = 0, x T ∈ X T lmost surely; @ISA sn order to ssure the fesiility of the integer prolemD U ad ∩ Z n×T = ∅ is lso ssumed to e stis(edF es in setion PFQD we will de(ne

U int,ad := (u t ) T -1 t=0 A t u t ∈ X T , u t ∈ U t ∩ Z n , t = 0, . . . , T -1 . @ITA nd U int,ad t (x t ) := u t x t + T -1 s=t A s u s ∈ X T , u s ∈ U s ∩ Z n , s = t, . . . , T -1 . @IUA INRIA
Energy contracts management by stochastic programming techniques W yur heuristi method onsists in (nding the optiml integer ontrol using the fellmn vlue otined from the ontinuous @relxtionA prolemX

Q(t, x t , ξ t ) = inf c t (ξ t ) • u t + Q * * (t + 1, x t+1 , ξ t ) sujet to x t+1 = x t + A t u t , (fesiility) u t ∈ U int,ad t (x t ), (optimlity ut) Q * * (t + 1, x t+1 , ξ t ) ≥ E Q * * (t + 1, x t+1 , ξ t+1 )ṽF t , Q * * (t + 1, x t+1 , ξ t+1 ) ≥ x * • x t+1 -Q * (t + 1, x * , ξ t+1 ), ∀x * .
@IVA where Q * * (t + 1, x t+1 , ξ t ) is identil to its equivlent in @IRA orresponding to the ontinuous relxtion formultionF 3

Vectorial Quantization Tree he quntiztion methods were introdued in the erly IWSH9s nd hve een pplied in informtion siene nd signl proessingF eentlyD flly nd gès used this tehnique in order to disretize the proility spe nd lter pplied it to (nnil engineering ID PF sn this setionD the min ides nd some importnt results of vetoril quntiztion re presentedF 3.1 Optimal quantization yptiml quntiztion of rndom vetors onsists in (nding the est possiE le pproximtion @in L p norm senseA of R d -vlued rndom vetor ξ ∈ L p (Ω, F, P; R d ) y mesurle funtion φ(ξ) tking vlues in (nite set Γ of rdinlity t most N vluesF reneD we wish to solve the minimiztion prolemX

min ξ -φ(ξ) p φ : R d → R d , mesurle, #(φ(R d )) ≤ N @IWA
nd let us denote y Γ := φ(R d ) the set of quntized pointsF vet us denote the set of quntized points y ξ := φ(ξ)F he ove distne @IWA is de(ned s the miniml distortion D µ,p N = ξξ p p of quntiztionD where µ is the proility mesure of ξF st is ovious to get tht φ = proj Γ D where proj Γ denotes the projetion over the quntized points set Γ = ξ 1 , . . . , ξ n following the nerest neighor rule of

L p distneF e sy tht C(ξ i ) is Voronoi tessellation of N Etuple Γ if for ll ξ i ∈ ΓD C(ξ i ) is forel set stisfyingX C(ξ i ) = proj -1 Γ (ξ i ) ⊂ ξ |ξ -ξ i | p = min ξ j ∈Γ |ξ -ξ j | p . @PHA
nd we de(ne the proility ssoited to ξ i X

p i = P(ξ ∈ C(ξ i )) @PIA
pollowing the de(nitionD N Etuple hs in(nitely mny oronoi tesselltionsF roweverD eh tesselltion lwys hs the sme losure nd the sme oundryF sf the distriution of ξ is strongly ontinuousD iFeFD the mesure µ of ny hyperE plne is 0D we will not other with oundries nd my think s if there ws only one tesselltionF his hypothesis is ssumed to e stis(ed in our pplitionF hereforeD the prolem of (nding the optiml quntiztion eomes glol optimiztion prolemD whih is in generl di0ult to solveF roweverD we n esily derive the following properties @eFgF PD property IEQAX he miniml distortion is deresing with respet to N F lim

N →∞ min ξ∈Γ,#(Γ)≤N ξ -ξ p = 0F
yne importnt property is the rte of the onvergene to zeroD whih is referred s dor9s theoremF Theorem 3.1. [12,Theorem 6.2

] If E[|ξ| p+η ] < ∞ for some η > 0, then, lim N N p/d min #(Γ)≤N ξ -ξ p p = J p,d |g| d/(d+p) (u)du 1+p/d where P(du) = g(u)λ d (du) + ν, ν ⊥ λ d (λ d Lebesgue measure on R d ). The constant J p,d corresponds to the case of the uniform distribution on [0, 1] d .
e n otin without di0ulty the following orollryX Corollary 3.2.

min ξ∈Γ,#(Γ)≤N ξ -ξ p = O(N 1/d ).
3.2 Quantization (recombining) tree he quntiztion tree ( ξt ) is (nite stte wrkov hinD suh tht t eh stge tD ξt is n optiml quntiztion of the originl rndom vrile ξ t in the disrete time proessF vet us denote y N t the numer of quntiztion points t stge tD nd set N := T t=0 N t suh tht N 0 := 1F e (rst hoose the numer of totl quntiztion points N D whih depends on the lultion pity of the mhine nd on the required preision of the resultF he next step is to determine the numer N t of quntiztion points tken t eh stge tF sf (ξ t ) ∈ L 2 (Ω, (F t ), P; R d ) D following PD proposition V in order to minimize the L 2 distortion of the whole proessD we get tht

N t = a t T τ =0 a τ N where a t := ξ t d/2(1+d) d/(d+2) t = 0, . . . , T.
@PPA e still need to evlute the proility trnsitions for the wrkov hinF vet p ij t e the trnsition proility from ξ i t to ξ j t+1 D de(ned s

p ij t := P proj Γt+1 (ξ t+1 ) = ξ j t+1 ṽproj Γt (ξ t ) = ξ i t @PQA
sn prtieD n nlytil omputtion of the proility trnsitions is out of reh for generl distriutionsF yne e0ient wy is to ompute the trnsition y wonteEgrlo methodD iFeF p ij t is estimted y pij

t pij t := Mp m=1 1 {proj Γ t (ξ m t )=ξ i t } 1 {proj Γ t+1 (ξ m t+1 )=ξ j t+1 } M m=1 1 {proj Γ t (ξ m t )=ξ i t } @PRA INRIA
Energy contracts management by stochastic programming techniques II where M p is the numer of trjetories in wonteEgrlo simultionF por the detil of quntiztion tree uildingD prtiulrly the extended omE petitive lerning vetor quntiztion @gvA lgorithmD we refer to PF Remark QFQ. )ug PR onsidered the prolem of generting smll size senrio tree @in order to e le to solve in n e'etive wy stohsti progrmming prolemsA nd set this prolem in the frmework of minimizing the distnes etween proility lwsF sn this pper we use quntiztion pproximtion tht lso minimizes distnes etween proility lwsD ut in di'erent setting sine we look for wrkovEhin pproximtionsF sn PRD )ug studied similr lgorithm to disretize ontinuousEstte proE esses into nonEreomining treeF fut the quntiztion theory nd gv lgorithm over wht he hs proposedF 3.3 Dual dynamic programming formulation on quantization tree he lgorithm presented in this setion for quntized prolem is similr to the one in wo et lF PPF sn the sequel we denote the quntized proess y ( ξt )D the ssoited (ltrtion y ( Ft )F end we represent the ssoited fellmn vlue y Q(t + 1, x t+1 , ξt ) to the prolem similr to @PA in whih rndom vriles (ξ t ) re repled y ( ξt )F he orresponding dynmi progrmming priniple reds

Q(t, x t , ξt ) := inf c t ( ξt ) • u t + Q(t + 1, x t+1 , ξt ) sujet to u t ∈ U ad t (x t ), x t+1 = x t + A t u t ; @PSA where Q(•, •, •) is the onditionl expettion of the vlue t time t + 1 otined in stte x t+1 X Q(t + 1, x t+1 , ξt ) := E Q(t + 1, x t+1 , ξt+1 )ṽ Ft @PTA nd Q(T, x T , ξT ) = g( ξT , x T ) if x T ∈ X T , +∞ otherwiseF @PUA
sn the lgorithm to e studied elowD only (nite numer of 0ne lower ounds of the fellmn vlues Q(t, •, ξt ) re villeF henote y M i t the orresponding set for the reliztion of ξ i t F e see tht the mximum of these 0ne lower ounds is itself onvex lower ound of the fellmn vlues

Q(t, x t , ξ t )D denoted y X θ(t, x, ξ i t , M i t ) := max x * • x -e; (x * , e) ∈ M i t .
@PVA end the ssoited onditionl expettion ppering in the dynmi progrmE ming eqution is then lower ound of Q(t + 1, x t+1 , ξ t )D when the fellmn vlues Q(t+1, x t+1 , ξ t+1 ) re repled y the lower ounds θ(t+1, x,

ξ i t+1 , M i t+1 )X ϑ(t+1, •, ξ i t , M t+1 ) := E[θ(t+1, •, ξ t+1 , M j t+1 )ṽ Ft ] = ξ j t+1 ∈Γt+1 pij t θ(t+1, •, ξ j t+1 , M j t+1 ) @PWA RR n°7289
where M t+1 := M j t+1 : ξ j t+1 ∈ Γ t+1 F e my then write the resulting approximate dynamic programming equa-tionD in whih the fellmn vlue is repled y its lower estimteD in the following formX

inf c t (ξ i t ) • u t + ϑ(t + 1, x t+1 , ξ i t , M t+1 ) sujet to u t ∈ U ad t (x t ), x t+1 = x t + A t u t , ϑ(t + 1, x t+1 , ξ i t , M t+1 ) = ξ j t+1 ∈Γt+1 pij t θ(t + 1, x t+1 , ξ j t+1 , M j t+1 ), θ(t + 1, x t+1 , ξ j t+1 , M j t+1 ) ≥ x * • x t+1 -e, ∀(x * , e) ∈ M j t+1 .
@QHA trting from the initil stte x 0 = 0D we egin y generting trjetories for (ξ t )D nd for eh of these trjetoriesD we solve the ove liner progrm t eh time stepF pinllyD we otin suoptiml poliy @in the spirit of pproxiE mte dynmi progrmmingAF e will ll this lultion forwrd opertionF imultneouslyD eh time we solve @QHAD from the dulity theory for liner proE grmsD we otin n 0ne lower ound of θ(t, •, ξt ) @nd hene lso of Q(t, •, ξt )AF his is wht we will ll kwrd elementry stepF e perform kwrd indution y propgting this informtion over prtiulr trjetoryD from the (nl time to the initil oneF sn the ove disussion we do not use nymore the originl @non quntizedA rndom vrileF roweverD it is soundD in order to otin n upper ound of the optiml vlueD to use the originl proility lw for omputing trjetoriesF et the sme timeD for omputing the deision vrile t eh stepD we solve the pproximte dynmi progrmming prolem where ξt is tken s the projetion of ξ t over the quntized gridF he prolem to e solved then is

inf c t (ξ t ) • u t + ϑ(t + 1, x t+1 , proj Γt (ξ t ) = ξ i t , M t+1 ) sujet to u t ∈ U ad t (x t ), x t+1 = x t + A t u t , ϑ(t + 1, x t+1 , ξ i t , M t+1 ) = ξ j t+1 ∈Γt+1 pij t θ(t + 1, x t+1 , ξ j t+1 , M j t+1 ) θ(t + 1, x t+1 , ξ j t+1 , M j t+1 ) ≥ x * • x t+1 -e, ∀(x * , e) ∈ M j t+1 . @QIA
Remark QFR. his formultion @QIA orresponds to the multiut version in liner reourse prolemD see SD setion SFQF 3.4 Aggregated formulation roweverD during the numeril testsD we me to relise tht this formultion @QIA ws still too omputtionlly intensiveF e therefore provide n ggregted version of the ove prolemX

inf c t (ξ t ) • u t + ϑ(t + 1, x t+1 , proj Γt (ξ t ) = ξ i t , M i t+1 ) sujet to u t ∈ U ad t (x t ), x t+1 = x t + A t u t , ϑ(t + 1, x t+1 , ξ i t , M i t+1 ) ≥ x * • x t+1 -e; ∀(x * , e) ∈ M i t+1 , @QPA INRIA
Energy contracts management by stochastic programming techniques IQ where the set M i t+1 of 0ne lower ounds is otined y ggregting the new lower ounds for eh ξ j t+1 ∈ Γ t+1 @eh time they re genertedA into single one tking into ount the proility trnsitionsX

M i t+1 :=      (x * , e) = ξ j t+1 ∈Γt+1 pij t ((x * ) j , e j ), ((x * ) j , e j ) ∈ M j t+1     
. @QQA hereforeD this expression still provides lower oundsF st undoutedly loses some informtion provided y optimlity uts in @QIAF reneD more itertions will e required to onverge @s proved in our numeril testsAF fut this formuE ltion simpli(es suprolem @QIAD nd numeril tests indite tht it gretly redues the g timeF Remark QFS. his formultion @QPA orresponds to the vEhpe method in liner reourse prolemD see SD setion SFIF he version of heuristi method @IVA for the integer solutions follows the sme ideX

inf c t (ξ t ) • u t + ϑ(t + 1, x t+1 , proj Γt (ξ t ) = ξ i t , M i t+1 ) sujet to u t ∈ U int,ad t (x t ), x t+1 = x t + A t u t , ϑ(t + 1, x t+1 , ξ i t , M i t+1 ) ≥ x * • x t+1 -e; ∀(x * , e) ∈ M i t+1 .
@QRA roweverD it is only vlid for the forwrd pss sine we nnot otin the dul informtion from mixed integer progrmming prolemF 4 Algorithm efter hving disretized the proility speD it is nturl to use the stohsti dul dynmi progrmming method on the resulting (nite stte wrkov hinF ine the fellmn vlue @nd its pproximtion y olletion of optimlity utsA is nonEonvex with respet to the underlying rndom vrileD we need to lulte the uts on every vertex of the quntiztion treeF he lgorithm itertes over forwrd pss nd kwrd pssF he forwrd pss genertes set of rndom trjetories following the dynmi of wrkov hin @IA nd runs the optimiztion on the )y on these trjetoriesF st provides n upper ound of the optiml vlueF he kwrd pss updtes the optimlity uts on eh vertex nd omputes lower ound t the (rst stgeF 4.1 Forward pass he forwrd pss genertes trjetories ξ m , m = 1, . . . , M f following the dyE nmi of wrkov hin @IAD iFeF out of smple for quntiztion treeF por eh trjetoryD strting from the (rst stge @x 0 = 0AD the urrent deision u t nd the stte vrile x t+1 for the next stge re omputed y solving @QPAF efter hving omputed ll the trjetoriesD n ojetive vlue v m for eh trjetory m = 1, . . . , M f is produedF husD the forwrd sttisti n e oE tined s followsX

v = 1 M f M f m=1 v m nd s = 1 M f M f m=1 (v m -v) 2 .
Remark RFI. elso notie tht the ontrol solution (u t ) omputed in this lgoE rithm is suEoptiml sine we pproximte Q(t, x t , ξ t ) y θ(t, x t , ξt , M t ) @respF Q(t + 1, x t+1 , ξ t ) y ϑ(t + 1, x t+1 , ξt , M t+1 )AF roweverD (u t ) elongs to the fesile spe (U ad t ) @it neessrily elongs to L 2 sine U t is omptAD hene is fesile for the originl prolemF es it is minimiztion prolemD eh v m is suEoptimlD whih implies tht v is n upper ound in the sttistil sense for the originl prolemD the 95% symptoti on(dene intervl is given y [v -1.96s, v + 1.96s]F Remark RFP. vike ny wonteEgrlo type methodsD the onvergene rte is 1/ M f D whih is slowF sn dditionD the numer of suprolems to solve in the forwrd pss is proportionl to the numer of forwrd pss trjetories M f F hereforeD some vrine redution tehniques re very useful hereF sn prtieD we use quite e0ient ontrol vrite tehniqueD see setion SF Remark RFQ. ine the lultion of eh trjetory is independentD it is possile to use prllel omputing in order to speed up the omputtionF 4.2 Backward pass he gol of the kwrd pss is to updte the optimlity uts used in @QPAF st follows the quntiztion treeD iFeF in smple lultionF e still use the formuE ltion of suEprolem @QPAD where ξt ∈ Γ t nd x t hve een stored during the previous forwrd pssF yne the suEprolem @QPA hs een solvedD the optiE mlity ut is omputed from the dul vlue ssoited to the dynmi eqution x t+1 = x t + A t u t F e then dd this ut into the optimlity uts olletion of vertex M i t F

Remark RFR. yptimlity uts hve to e updted for eh vertex of the qunE tiztion tree due to the nonEonvexity of Q(t, x t , ξ t ) @respF Q(t, x t , ξt )) with respet to ξ t @respF ξt AF he numer of suEprolems to solve in the kwrd pss is proportionl to the quntiztion tree size N times the numer of trjeE tories denoted M b F ine the stte vrile vlues re stored during the forwrd pssD M b is lwys smller thn M f F fut if we re using ll M f trjetories for the kwrd pss @M f is reltively lrge owing to remrk RFPAD the numer of suEprolems would e N × M f whih is omputtionlly unmngeleF oD one ide is to selet only smll susmple of M b smples to e omputed in the kwrd pss in order to llow the progrm to (nish in resonle timeF enother reson is tht in prtieD some stte vlues x t re very lose to eh otherD using the suEprolems to return lmost the sme dul informtion whih will not help us muh pproximte well Q(t, x t , ξt )F roweverD the proess of seleting M b trjetories mong the M f ones deE pends on the nture of the prolem on whih the lgorithm is ppliedF sn the numeril tests in setion SD the seletion follows the rejetion method whih onsists in heking whether the distne etween the new x t nd the existing ones is elow prede(ned threshold or notF his is rndom wy of stisfying the requirement tht they re not too lose to eh otherF INRIA Remark RFS. imilrly to the forwrd pssD the kwrd pss lgorithm my lso use prllel omputtionD euse t stge tD the suEprolems @QIA using di'erent (x t , ξ i t ) re independent from eh otherF efter solving the suEprolem t the (rst stgeD n optiml vlue v is oE tinedD whih is lower ound for the quntized prolemD sine this vlue is lulted y fenders deompositionF roweverD v is lulted ording to the quntiztion treeD nd would di'er if lulted on nother treeF he erE ror used y this pproximtion n e oundedD s disussed in the next susetionF Remark RFT. he stopping ondition is to hek whether v ∈ [v-1.96s, v+1.96s]F roweverD the on(dene intervl is sometimes too lrge mking this ondition too strightforwrd to e stis(edF o void this prolemD we dd nother stopping ondition stting tht the kwrd pss vlue v should e stleX |v it -v it-1 | ≤ |v it |D where it ∈ Z ++ is the numer of itertion of the lgorithmF eentlyD hpiro hs exposed in PV tht the ove stopping ondition is reltively wek if the on(dene intervl is lrgeF re proposed tht the upper ound is v + 1.96s nd the following stopping onditionX

v ∈ [v + 1.96s -ε, v + 1.96s + ε] @QSA
where ε is prmeterF sn this rtileD we tke the stopping onditions

v ∈ [v -s, v + s] @QTA
where is prmeter whih is striter thn the one proposed y ereir nd into PQF

Convergence result and error analysis on quantization

sn this setionD we estlish the onvergene result of the ove lgorithm nd provide n error nlysis of disretiztion y quntiztionF e (rst show the onvergene of the lgorithm when the forwrd pss follows the disretized distriution ! quntiztion treeF es we hve notied in remrk QFRD the lgorithm is equivlent to vEhpe methodF he onvergene result on quntized forwrd pss follows P rgumentsF he onvergene result of vEhpe method on (nite proility distriuE tion is given in firge nd vouveux SD hpterUD theorem ID where the min rgument is (niteness of optimlity uts nd fesiility utsF husD fter generting ll fesile uts nd optimlity utsD the kwrd pss vlue does not inrese nymoreF ine now the forwrd pss follows disretized distriutionD the forwrd vlue of wonte grlo method is itself rndom vrile whose men is the upper ound of vEhpe method whih onverges to lower ound y (nite onvergene of vEhpe methodF Remark RFU. sn prtieD we n (rst simulte the forwrd pss following the disretized distriution nd let the forwrd vlue nd kwrd vlue onvergeF hen we n hnge the forwrd simultion to following originl distriution to ompute the (nl resultF xowD we study the onvergene of the upper ound nd the lower ound one we let forwrd pss follow ontinuous distriutionF sn order to estimte the error used y quntiztionD we refer to the result in reitsh nd ömish IRD ISF Theorem 4.8. [14, Theorem 2.1] Assume that the solution set U ad is nonempty, and bounded. Then, there exists positive constants L and δ such that the estimate

|v(ξ) -v( ξ)| ≤ L( ξ -ξ 2 + D f (ξ, ξ)) @QUA
holds for all random element ξ ∈ L 2 (Ω, (F t ), P; R d ) with ξξ 2 ≤ δ, where

D f (ξ t , ξt ) = inf u∈S(ξ),û∈S( ξ) T -1 t=2 max u t -E[u t ṽF t ( ξ)] 2 , ût -E[û t ṽF t (ξ)] 2 , @QVA
where S(ξ) (resp. S( ξ)) is the optimal solution of with random variable ξ(resp. ξ), and v(ξ) (resp. v( ξ)) are the related optimal values. ine ( ξt ) is the quntiztion of

(ξ t )D then F t ( ξ) ⊂ F t (ξ)F he optiml ontrol ût is F t ( ξ) mesurleD hene it is F t (ξ) mesurleF henD we hve for @QVAX D f (ξ t , ξt ) = inf u∈S(ξ),û∈S( ξ) T -1 t=2 u t -E[u t ṽF t ( ξ)] 2 @QWA
yur frmework is little di'erent from the prolem studied in IRD ISD sine our prolem possesses (nl ost funtion g(ξ T , x T +1 )F roweverD their result is still suitle in our frmework if the vipshitz property of this funtion holdsF e will omit the proofD whih is similr to tht in IRF his error ound onsists in P termsD the (rst one of the right hnd side @QUA is the l 2 Edistne etween the originl rndom proess nd the pproximted oneD the seond one is the distne etween the P (ltrtions over the optiml ontrolF he (rst term n e mesured y theorem QFID whih provides the onvergene rte to HF roweverD it is very di0ult to provide onvergene rte to @QWAF Remark RFW. here re other referenes whih estlish stility results of opE timl vlue of ontrol wrkov proessF vet us mention qordienko et lF II on verge ost wrkov ontrol proess sed on untorovih distne nd the referenes thereinF roweverD their prolem is to minimize n verge or disE ounted ost for n in(nite horizon prolemF hereforeD pplying their nlyst in our lgorithm is not oviousF 4.4 Algorithm he whole lgorithm is summrized in (gure IF

5

Numerical tests sn this setionD two numeril tests re performedF he (rst one is swing optionD whih is n exmple where oth stte vrile nd rndom vrile re slrF he seond one onerns vxq portfolio optimiztionD whih is fr more omplexF he lgorithm hs een suessfully pplied on oth exmplesF

INRIA

Step 0 Construction quantization tree (possible o-line)

uild one quntiztion tree nd ompute trnsition proilityF

Step 1 Initialisation let u 0 ∈ U ad 0 the urrent optiml solution nd ompute stte vrile vlue

x 1 Y let v 0 = -∞F
Step 2 Forward pass rete M f di'usion proesses following dynmi @IAY for t = 1, . . . , T -1 doX for m = 1, . . . , M f doX lulte ξm t = proj Γt (ξ m t )D solve the suprolem @QPAY ompute the forwrd sttisti v nd sF

Step 3 Backward pass for t = T -1, . . . , 1 doX for m = 1, . . . , M b doX solve the suprolem @QPA for ll verties in Γ t Y lulte new optimlity ut nd dd it into M t Y t = 0@(rst stgeA solve the (rst stge prolem nd lulte the kwrd vlue v it F

Step 4 Check stopping condition

if v ∈ [v -1.96s, v + 1.96s] nd |v it -v it-1 | ≤ |v it-1 | yY else go to tep PF pigure IX elgorithm
5.1 Price model sn this pperD we ssume tht the rndom proessD tking prt in the modeling of futures ommodity prie evolutionD follows the elerted flk model TF elthough simplistiD this model is equivlent to the one ftor qussin rtw model @see UA whih is well dpted for forwrd urve modelingF his model is disussed in the prtieEsed litertureD iydelnd nd olynie WD ghpter SD nd vi et lF IHF ine we re using the vetoril quntiztion method to disretize the underlying rndom proessD we re not limited to tht modelF yther (nnil prie models n lso e ppliedF he prie model tkes s prmeters the mrket voltilities σ s well s the forwrd urves F 0 F he flk model sttes tht future pries of ommodities re mrtingles whose expettions orrespond to the originl forwrd urve level F 0 due to lk of ritrge opportunitiesF he forwrd ontrt prie F t s t time sD with mturity @lso lled tenorA tD follows under the equivlent mrtingle mesure @iwwA P the following dynmiX

ln (F s t+1 ) i (F s t ) i = σ i W i t - 1 2 (σ i ) 2 i = 1, . . . , d; @RHA
where W i t is stndrd norml distriution N (0, 1)D with orreltion corr(W i t , W j t ) = ρ ij F henD the spot pries t time t re simply otined y tking the tEexpiring forwrd ontrts t time tX ξ t = F t t D

ξ i t := (F t t ) i = (F t 0 ) i exp t-1 s=0 σ i W i s - 1 2 (σ i ) 2 t i = 1, . . . , d; @RIA
where (F t 0 ) is the originl forwrd urve oservedF etting

S i t := exp t-1 s=0 σ i W i s - 1 2 (σ i ) 2 t ,
we then hve

ξ i t = (F t 0 ) i S i t F

Swing option

wing options re firly ommon in the energy mrketF st llows its holder to purhse totl mount of ommodity during (xed time periodF e si swing option ontrt n e formulted sX

max E T -1 t=0 (ξ t -K t ) • u t sujet to u t ∈ L 2 (Ω, (F t ), P; [0, 1]) T -1 t=0 u t ∈ [L, U ] lmost surely @RPA
where K t is the strike prieD ssumed to e deterministiD L nd U re lower nd upper ounds of the (nl onditionF vet us (rst )ip the sign of the ojetive funtion to (K t -ξ t ) in order to swith to minimiztion prolemF sn this testD we onsider tht the strike prie equls to the originl forwrd prie oserved nowdys K t = F t 0 without loss of generlityF sn the se where (L, U ) ∈ Z 2 D frdou et lF Q show tht the optiml ontrol hs ngEng propertyD iFeF u t ∈ {0, 1}F husD dynmi progrmming n e diretly pplied in this frmework y disretizing the stte spe x t F frdou et lF hve pplied this ide to disretize the stte vrile nd use quntiztion tree to pproximte the underlying rndom speF roweverD in the presene of more omplex onditions on the ontrol vrile suh s u t ∈ [l(x t ), u(x t )] @like in most of storge pities tehnil feturesAD the ngEng property will not hold nymoreD nd furthermore the disretiztion of the stte spe x t my eome intrtleF 5.2.1 Control variate ine the forwrd pss produes lrge on(dene intervl with smll M f D the following ontrol vrite tehnique hs een introdued to redue itF vet us onsider the following vrileX

v cv (ξ) = T -1 t=0 (K t -ξ t ). @RQA yserve tht E[v cv (ξ)] = 0 sine E[ξ t ] = F t 0 = K t F INRIA 5.2.2 Numerical result
sn the following testD without losing ny generlityD we tke F (0, t) = 1, ∀t = 0, . . . , T -1 in the prie model @RHAF he results re ompred to nother method using ngEng strtegy nd the vongst'Ehwrtz method PIF huring the simultionD the quntiztion tree hs 3000 points in totlF e set prmeters M f = 1000 nd M b = 10 in wonteEgrlo simultionF e use 1000 senrios in the vongst'Ehwrtz methodF [4.04, 4.40] [5.91, 6.51] [7.86, 8.62] [1.77, 1.97] [2.66, 3.00] [3.68, 4.16] ffCvwF FiF [4.19,4.47] [6.23,6.71] [8.22,8.94] [1.81,2.17] [2.69,3.31] [3.57,4.42] le IX untiztion tree C dul dynmi progrmming versus ngEng strtegy C vongst' hwrtz method for priing swing optionF pigure PX hh itertion for priing swing option using quntiztion methodF [L, U ] = [20, 30]D σ = 0.5/ √ 50 pigure P shows the onvergene etween the upper nd lower ounds for the swing option priing prolemF st only tkes S itertions for the lgorithm to onverge whih represents less thn I minute of g timeF xotie lso tht the stndrd devition ssoited to the upper ound is less thn I7 of the optiml vlueF gompred to other numeril methodsD we see tht there is no rel dvntge in using this method in terms of omputtion timeF roweverD the exmple in the next setion will evidene its superiority over the other lgorithms whih fil to provide urte results in resonle timeF 5.3 Dynamic portfolio optimization e now onsider gs trding portfolioF e trding ompny purhses lique(ed nturl gs @vxqFA from set of produing ountries indexed t prie formul nd sells it to onsuming ountries t nother other prie formul @see (gure Q for the min mrket nd tle P for the prie formuleAF ennul quntity nd prie formule hve een greed ontrtullyD the ltter re funtions of the future pries of mjor energy mrkets ξ t D suh s rude oil prie @ysvAD north emerin nturl gs prie @xe xqAD nd iurope nturl gs prie @i xqAF he gol is to (nd strtegy tht optimizes the portfolio vlueF he deision is mde over the quntity of vxq trnsported on eh routeF At that point, we relax the discrete nature of the decision variables. he (rst di0ulty lies in the numer of indies we hve to onsider when we model the portfolioX ξ t ∈ R 3 F he seond di0ulty is linked to the dimension of the stte vrile x t ∈ R m D orresponding to the totl numer of prodution nd delivery ountriesD whih is lso in high dimensionD m = 6 in the simultionF pigure QX e (tive supply nd demnd portfolioD s well s the possile routesF X produing ountryY •X onsuming ountryF

T = 50 [L, U ] = [0, 25] [L, U ] = [
E T -1 t=0 ct (ξ t ) • u t sujet to u t ∈ U t @monthly onstrintA T -1 t=0
A t u t ∈ X T @nnul onstrintA @RRA vet us )ip the sign of the ojetive funtionX c t = -c t to swith to minimizE tion prolemF he forwrd prie F t 0 in @RHA is red from the energy mrketF rere is the vlue used in numeril test @see (gure R E SAF pigure RX ysv prie @6 GrrelA pigure SX i xq prie @6 G wwftuA pigure TX xe xq prie @6 G wwftuA he quntiztion tree ws proessed with N = 36000 R 3 Evlued elementry quntizers dispthed on T = 12 stgesF he estimtion of trnsition proility ws rried out y wonteEgrlo method using 10 9 smplesF 5.3.1 Control variate imilrly to the previous exmpleD we keep using ontrol vrite to redue the lrge on(dene intervlF vet us onsider the following vrileX

v cv (ξ) = T -1 t=0 c t (ξ t )u * t @RSA
where u * t is one of the optiml solutions for prolem

inf T -1 t=0 E [c t (ξ t )] u t sujet to u t ∈ U t T -1 t=0 A t u t ∈ X T @RTA
vet v * cv e the optiml vlueF his is deterministi liner progrmD whih is strightforwrd to solveF he only tehnil point is to evlute E[c t (ξ t )]F yne possile method onsists in using quntiztionX

E[c t (ξ t )] ≈ ξ i t ∈Γt p i t c t (ξ i t ) glerly we hve E[v cv (ξ)] = v * cv F
Remark SFI. he @RTA is nothing else ut the deterministi ounterprt formuE ltionD whose vlue is lso known in (nnil engineering s the intrinsi vlue of the optionF hus the optiml vlue v * cv is onsidered s good referene to ompre the optiml vlue otined y stohsti progrmmingF

Results

he prmeters driving the stohsti proesses in the test re set s followsX voltility σ 1 = σ 2 = σ 3 = 40%D nd orreltion ρ 12 = 0.7 ρ 13 = 0.2 ρ 23 = 0.4F end the prmeters driving the hh prt of the lgorithm re forwrd pss smples M f = 3000 to 6000D kwrd pss smple M b = 8 to 12F he progrm is written in gCC with gplex IHFID nd the tests were performed on g with PFP qrz hul gore ewh g nd ITqfyte min memoryF he numeril results re shown in pigure UF xotie tht the vlue of the kwrd pss is lmost stle from the 6th itertionD nd enters the on(dene intervl of the upper ound severl timesF he stopping ondition tken here is @QTA proposed y hpiro in PV where we tke = 0.5F purthermoreD the on(dene intervl of the upper ound generted with ontrol vrite represents ir 4% of the forwrd pss vlueF his shows tht the ontrol vrite proves to e dequte in this exmpleF woreoverD the lgorithm onverges in 17 itertionsD whih is e0ientF roweverD if we test this lgorithm with lrger voltilitiesD or with more produingGonsuming ountriesD more itertions will e required to onvergeF felow re the grphs representing the g time for eh itertion @(gure V nd (gure WAX his progrm does not onsume lot of time during eh itertionF rowever the lultion time inreses shrply when the itertion goes onF e know tht

  0, 20.0] [0.0, 4.0] 0.08OIL -0.8 le PX gonstrints nd prie formule * The quantity unit is in TBtu. MMBtu stands for a million British thermal unit, a TBtu is a trillion British thermal unit thus equivalent to a Million MMBtu. ** The price unit is $/MMBtu. oD the prolem is sup

  INRIAPQ pigure UX hh itertion result for portfolio optimiztionF pigure VX forwrd g time pigure WX kwrd g time the forwrd omputtion time is proportionl to the numer of forwrd smples M f D nd tht the kwrd omputing time is proportionl to the numer of verties in the quntiztion tree N nd to the numer of kwrd smples M b F reneD we n modify these prmeters to redue the lultion time when it eomes too longF fut if we do soD we will lose the preision y enlrging the on(dene intervlFInteger solution es explined in setion PFSD we n ompute one su optiml integer solution y omining the fellmn vlues of the ontinuous relxtion with the integer onstrints on the deision vriles during the forwrd stepF he results re presented in tle QF he ontrol vrite vlue of EVFQWRT is otined when solving @RTA tking ount of the integer onstrints on the onE trol vrileF he heuristi gives n expeted vlue of EIQFRIRPD with stndrd devition of HFWQQIUD for SHHH forwrd smplesF es mentioned in remrk SFID the deterministi ounterprt vlue gives good insight on the gin otined y using stohsti progrmmingF sn the ontinE upper ound stndrd devition of uFF lower ound intrinsi of optiml vlue of ontinuous prolem nd integer proE lem uous relxtion seD the optiml vlue got through stohsti progrmming is out S times etter thn the deterministi ounterprt vlueF gonerning the integer vrile prolemD lthough only n upper ound ould e genertedD it is nevertheless out SH7 etter thn the one of the deterministi ounterE prtF herefore we onlude tht the optiml vlue of the stohsti integer progrm lies somewhere etween -41.00 @ontinuous relxtion optiml vlueA nd -13.41 @integer suEoptiml vlueAF eduing this intervl is the ojet of our future workF Acknowledgments he uthors thnk to two nonymous referees for their useful remrksD in prtiulr y ringing our ttention to ppers wo et lF PP nd qordienko et lF IIF References I F flly nd qF gèsF irror nlysis of the quntiztion lgorithm for ostle prolemsF Stochastic Processes & Their ApplicationsD IHT@IAXI!RHD PHHQF P F flly nd qF gèsF e quntiztion lgorithm for solving disrete time multiEdimensionl optiml stopping prolemsF BernoulliD W@TAXIHHQ!IHRWD PHHQF Q yF frdouD F fouthemyD nd qF gèsF yptiml quntiztion for the priing of swing optionsF Applied Mathematical FinanceD IT@PAXIVQ!PIUD PHHWF R hFF fertseks nd tFxF sitsiklisF Neuro-Dynamic ProgrammingF ethen ienti(D IWWTF S tFF firge nd pF vouveuxF Introduction to Stochastic ProgrammingF pringerD xew orkD PHHHF T pF flkF he priing of ommodity ontrtsF Journal of Financial Econ-omyD QXITU!IUWD IWUTFU FeF trrow hFg reth nd eF wortonF fond priing nd the term struE ture of interest rtes X new methodology for ontingent lim vlutionF EconometricaD THXUU!IHSD IWWPF V tF hupováD xF qröweEuuskD nd F ömishF enrio redution in stohsti progrmming!n pproh using proility metrisF Math. Program. Ser. AD WSXRWQ!SIID PHHQF INRIA Energy contracts management by stochastic programming techniques PS W eF iydelnd nd uF olynieF Energy and Power Risk Management: New Developments in Modeling, Pricing and HedgingF tohn iley 8 onsD snFD rookenD xtD PHHQF IH qFviD pF wrgotD nd xF eomndiF en pproximte dynmi progrmE ming pproh to enhmrk prtieEsed heuristis for nturl gs storE ge vlutionF Oper. Res.D SVXSTR!SVPD wy PHIHF II ivgueni qordienkoD inrique vemusEodrÃguezD nd úl wontesEde yF everge ost mrkov ontrol proessesX stility with respet to the kntorovih metriF Mathematical Methods of Operations ResearchD UHXIQ! QQD PHHWF IHFIHHUGsHHIVTEHHVEHPPWETF IP F qrf nd rF vushgyF Foundations of quantization for probability distri-butionsF pringerEerlgD ferlinD PHHHF IQ rF reitsh nd F ömishF enrio redution lgorithms in stohsti progrmmingF Computational Opti. and Appli.D PRXIVU!PHTD PHHQF IR rF reitshD F ömishD nd gF trugrekF tility of multistge stohsE ti progrmsF SIAM J. on OptimizationD IUXSII!SPSD PHHTF IS rolger reitsh nd erner ömishF tility nd senrio trees for mulE tistge stohsti progrmsF sn qerd snfngerD editorD Stochastic Program-mingD volume ISH of International Series in Operations Research & Management ScienceD pges IQW!ITRF pringer xew orkD PHIIF IT uF roylndD wF uutD nd FF lleF e heuristi for momentEmthing senrio genertionF Computational Optimization and ApplicationsD PRXITW! IVSD PHHQF IU qF snfnger nd hFF wortonF gut shring for multistge stohsti liner progrms with interstge dependenyF Mathematical ProgrammingD USXPRI! PSTD IWWTF IV F ull nd FF lleF Stochastic ProgrammingF tohn iley 8 onsD ghihesterD IWWRF IW wF uut nd FF lleF ivlution of senrioEgenertion methods for stohsti progrmmingF Pacic Journal of OptimizationD Q@PAXPSU ! PUID PHHUF PH gF uühler nd F igerskeF heomposition of multistge stohsti proE grms with reomining senrio treesF PI pF vongst' nd iF hwrtzF luing merin options y simultion X simple lest squres pprohF Review of Financial StudiesD IRXIIQ!IRVD PHHIF PP fF woD eF qjelsvikD nd eF qrundtF sntegrted risk mngement of hydro power sheduling nd ontrt mngementF IEEE Transactions on Power SystemsD IT@PAXPIT!PPID PHHIF PQ wFFpF ereir nd vFwFFqF intoF wultiEstge stohsti optimiztion pplied to energy plnningF Mathematical ProgrammingD SPXQSW!QUSD IWWIF RR n°7289

  t is the ontrol vrileD ssumed to e squred integrleD tking vlue in ompt nonempty onvex polyhedrl set U t ⊂ R n D U t := L 2 (Ω, F t , P; U t ) represents the fesile speD x t ∈ R m is the stte vrile representing the inventory levelsD whih is oviously F t-1 mesurle y forwrd indutionD X T is the set of dmissile (nl stteD ssumed to e nonempty polyhedrl suset of R m D tking ple in the (nl stge onstrint over stte vrile x T D A t ∈ R m×n desries how the deisions impt the inventoriesD c t (•) : R d → R n is running ost unitD ssumed to e vipshitzD nd (nlly g(ξ, x) : R d ×R m → R is onvex lower semiEontinuous @lFsFFA funtion of xD uniformly lower ounded nd n often e interpreted s penlty funtion for (nl stte x T

@PA

INRIA where u

  For any t = 0 to T , the Bellman value functions Q(t, x t , ξ t )

	Energy contracts management by stochastic programming techniques	U
	ine X T nd U t re polyhedrl setsD so is X t F sf we reple the deision u t y
	sequene (u t , . . . , u T -1 ) suh tht x t ∈ X t D we get the property of reltively
	omplete reourseF he ounterprt is n inrese of dimension of the deision
	vrileD whih in our vxq pplition @see seond exmple in setion SA is
	eptleF
	2.4 Stochastic dual dynamic programming (SDDP) for-
	mulation
	e next tke dvntge of prtil onvexity property of the fellmn funtionsF
	Proposition 2.3.
		WA
		INRIA

RR n°7289

http://www.inria.fr ISSN 0249-6399


Zhihao Cen
email: zhihao.cen@polytechnique.edu

J Frédéric Bonnans
email: frederic.bonnans@inria.fr

Christel Thibault
email: thibault.christel@total.com

Sensitivity analysis of energy contracts management problem by stochastic programming techniques

Keywords: sensitivity, stochastic programming, Danskin theorem, dual dynamic programming

We consider a model of medium-term commodity contracts management. Randomness takes place only in the prices on which the commodities are exchanged whilst state variable is multi-dimensional. In [9], we proposed an algorithm to deal with such problem, based on quantization of random process and a dual dynamic programming type approach. We obtained accurate estimates of the optimal value and a suboptimal strategy from this algorithm.

 and other classical methods are provided and evidence the accuracy of the estimate of marginal prices.

Analyse de sensibilité de problémes de gestion de contrats d'énergie par des techniques de programmation stochastique

Résumé : Nous considérons un modèle de gestion de contrats de produit de base en moyen terme. Les aléas ne portent que sur le dynamique de prix. Dans [START_REF] Bonnans | Energy contracts management by stochastic programming techniques[END_REF], nous avons proposé un algorithme pour ce type de problème, basé sur la méthode de quantication de processus aléatoire, et l'approche de programmation dynamique duale. Nous avons obtenu un estimateur précis de la valeur optimale et une stratégie sous-optimale par cet algorithme. Dans cet article, nous analysons la sensibilité par rapport aux paramètres dans le processus aléatoire (modèle de prix). Nous discutons un estimateur de coût marginal basé sur le théorème de Danskin. A la n, des tests numériques appliqués à des problèmes réels dans le marché l'énergie ont été réalisés. Les comparaisons du résultat par la méthode dans [START_REF] Bonnans | Energy contracts management by stochastic programming techniques[END_REF] et celui par des autres méthodes, donnent des estimateurs précis sur le coût marginal.

Mots-clés :

sensibilité, programmation stochastique, théorème Danskin, programmation dynamique duale

We consider a model of medium-term commodity contracts management.

Randomness takes place only in the prices on which the commodities are exchanged whilst state variable is multi-dimensional. In [START_REF] Bonnans | Energy contracts management by stochastic programming techniques[END_REF], we proposed an algorithm to deal with such a problem, based on quantization of the random process and a dual dynamic programming type approach. We obtained accurate estimates of the optimal value and a suboptimal strategy from this algorithm. In this paper, we analyse the sensitivity with respect to parameters driving the price model. We discuss the estimate of marginal price based on the Danskin's theorem. Finally, some numerical results applied to realistic energy market problems have been performed. Comparisons between results obtained in [START_REF] Bonnans | Energy contracts management by stochastic programming techniques[END_REF] and other classical methods are provided and give evidence of the good accuracy of the estimate of marginal prices.

1 Motivation

We study a class of problems in the following settings. First of all, we consider a uniform discretization of time horizon [0, T ], 0 = t 0 < t 1 < • • • < t K = T , and a discrete time Markov process (ξ k ) := (ξ t k ) (in the following, we will replace in index t k by k for sake of clarity) in the probability space L 2 (Ω, (F k ), P; R d ).

The canonic ltration associated with (ξ k ) is denoted by F k := σ(ξ s , 0 ≤ s ≤ k).

We may write its dynamics as

ξ k+1 = f (W k , ξ k , α k ) 0 ≤ k ≤ K -2, (1)
where ξ 0 is deterministic, (α k ) is a sequence of deterministic parameters, and (W k ) is a sequence of independent squared integrable random variables, independent of F k , nally f is a measurable mapping.

The stochastic dynamic decision problem has the following expression:

inf

(u k )∈U E K-1 k=0 c k (ξ k ) • u k subject to u k ∈ U k , almost surely, x k+1 = x k + A k u k , x 0 = 0,
x T ∈ X T almost surely;

(

where u k ∈ R n is the control variable, U := L 2 (Ω, (F k ), P; R n ) is the functional space of the control variable u k , U k is a compact nonempty polyhedral set in R n , x k ∈ R m represents the state variable, X T is the set of admissible nal states assumed to be a nonempty polyhedral set in R m taking place in the nal stage constraint over state variable x T , A k ∈ R m×n is the dynamic matrix, c k (•) : R d → R n is the running cost per unit, assumed to be Lipschitz.

Such a formulation (2) is motivated by the fact that most medium and long-term energy contracts can be modelled as multi stage stochastic dynamic programs. Indeed, most of the contracts dene the terms on how the commodity will be exchanged throughout the duration of the contract. These terms generally state that the quantity to be taken at time t k (modelled by the control u k ) is bounded by minimal and maximal amounts. They also specify the price at which the commodity will be exchanged at time t k ; this price is generally a function of the spot price (modelled by c k (ξ k )) observed at that period. Since the spot price is unknown at the beginning, the problem becomes stochastic.

The aim of the holder of such a contract is to optimize its decisions with respect to the price information available to him at that time. Here, we choose to maximize the expectation of the revenue. In order to be compliant with the formulation of (2), we will ip the sign of the real cost-revenue functions to transform a maximizing problem into a standard minimization problem.

Most of the articles in the literature focus on the procedures on how to solve problem (2) and hence how to get the contract's price (see Carmona and Touzi [START_REF] Carmona | Optimal multiple stopping and valuation of swing options[END_REF], Carmona and Ludkovski [START_REF] Carmona | Valuation of energy storage: an optimal switching approach[END_REF], and Barrera-Esteve et al [START_REF] Barrera-Esteve | Numerical methods for the pricing of swing options: a stochastic control approach[END_REF]). Although essential, this information needs to be coupled with sensitivity information in a trading perspective. Indeed, the values of sensitivities (called the greeks) allow traders to replicate the price variations originating from the contracts they have to manage on a regular basis. See [START_REF] Hull | Options, Futures and other derivatives[END_REF]Chapter 17] for further discussions on hedging strategies using the greeks applied to vanilla options.

A lot of eorts have been made in order to assess proper sensitivities from discontinuous payos such as digital or barrier options. However, only a handful of papers focus on optimal exercise type options and most of them concern Bermudan options which are common in the banking industry. Recently, Cont and Fournié have developed functional Itô formulae in order to study the sensitivities of path dependent options without control [START_REF] Cont | Functional Ito calculus and stochastic integral representation of martingales[END_REF][START_REF] Cont | A functional extension of the ito formula[END_REF]. Our purpose in this paper is to provide an alternative method to estimate sensitivities for problem (2) based on Danskin's theorem which is well known among optimization practitioners. The advantage of Danskin's theorem is that it can be associated to any discretization type method such as partial dierential equations (PDE), Monte Carlo or optimal quantization. Also, it does not require to reestimate the conditional expectations once they have been computed during the price estimation procedure.

This paper is organised as follows: section 2 presents the resolution framework of problem (2) articulated around optimal quantization and stochastic dual dynamic programming (SDDP). Section 3 introduces the model driving the prices and section 4 focuses on Danskin's theorem through the computation of sensitivities. Convergence results are also provided in section 4. Finally, numerical tests are performed in section 5.

A review of quantization discretization and stochastic dual dynamic programming approach

In this section, we briey present the algorithm which was introduced in our previous article [START_REF] Bonnans | Energy contracts management by stochastic programming techniques[END_REF].

Discretization

In many real life problems, the discrete time random process (ξ k ) may possibly take values within an innite set. In order to make the problem numerically tractable, we start by discretizing the random process by a tree. Because of the Markov property of the random process, a recombining tree is the method of choice for reducing complexity. Bally et al [START_REF] Bally | A quantization algorithm for solving discrete time multi-dimensional optimal stopping problems[END_REF][START_REF] Bally | Error analysis of the quantization algorithm for obstacle problems[END_REF] introduced the vectorial quantization tree method as a new discretization scheme. The main idea is to replace ξ k by the approximation ξk whose support is a nite set

Γ k = ξ i k , 1 ≤ i ≤ N k .

INRIA

The quantization ξ is said to be optimal if it minimizes the error associated with the L 2 norm:

ξ ∈ argmin ξ -ξ 2 := K-1 k=0 ξ k -ξk 2 2 1/2 : ξk ∈ Γ k , K-1 k=0 #(Γ k ) = K-1 k=0 N k ≤ N ( 3
)
where N is the number of quantized points in the tree. Then, we dene a

Voronoi tessellation associated with the quantized points ξ i k ∈ Γ k by

C i k := ξ k : |ξ k -ξ i k | ≤ min j =i |ξ k -ξ j k | .
The nal step consists in building the recombining tree by computing probability transitions:

p ij k := P ξk+1 = ξ j k+1 | ξk = ξ i k = P ξ k+1 ∈ C j k+1 | ξ k ∈ C i k
based on the law of (ξ k ). Bally and Pagès also proposed in [START_REF] Bally | A quantization algorithm for solving discrete time multi-dimensional optimal stopping problems[END_REF] the competitive learning vectorial quantization (CLVQ) algorithm aimed at building a vectorial quantization tree based on the stochastic gradient method.

Stochastic dual dynamic programming algorithm

After having discretized the random process (ξ k ) with the help of vectorial quantization, we now focus on the resolution of problem [START_REF] Bally | Error analysis of the quantization algorithm for obstacle problems[END_REF]. This problem can be naturally decomposed in time by the dynamic programming principle. We denote by Q(k, x k , ξ k ) the Bellman value:

Q(k, x k , ξ k ) = inf (us)∈U E K-1 s=k c s (ξ s ) • u s subject to u s ∈ U s almost surely, x s+1 = x s + A s u s , x 0 = 0, x T ∈ X T almost surely. (4)
Then, the dynamic programming principle can be written as

Q(k, x k , ξ k ) = inf u k ∈U k c k (ξ k ) • u k + Q(k + 1, x k+1 , ξ k ) (5)
where

Q(k + 1, x k+1 , ξ k ) = E[Q(k + 1, x k+1 , ξ k+1 ) | F k ]
and nal time condition:

Q(T, x T , ξ T -1 ) = χ X T (x T ) = 0 if x T ∈ X T ∞ otherwise
In that perspective, we have opted for the stochastic dual dynamic programming (SDDP) algorithm which is well suited to our problem. The stochastic dual dynamic programming method was rst introduced by Pereira and Pinto [START_REF] Pereira | Multi-stage stochastic optimization applied to energy planning[END_REF].

Recently, Philpott and Guan [START_REF] Philpott | On the convergence of stochastic dual dynamic programming and related methods[END_REF] and Shapiro [START_REF] Shapiro | Analysis of stochastic dual dynamic programming method[END_REF] have analysed the convergence property of this method. The main idea of SDDP algorithm is to approximate the Bellman value function Q(k, x k , ξ k-1 ) by the supremum of a family of ane function of x k , depending on t k and ξ k , denoted by ϑ(k + 1, x k+1 , ξk , O ξk k+1 ).

The method is based on successive iterations of forward and backward passes.

The forward pass generates a set of random trajectories following the dynamic of (ξ k ) and runs the optimization on the y over these trajectories using the following approximation of dynamic programming:

inf u k ∈U k c k (ξ k ) • u k + ϑ(k + 1, x k+1 , ξk , O ξk k+1 ) subject to x k+1 = x k + A k u k , ϑ(k + 1, x k+1 , ξ i k , O ξk k+1 ) ≥ x * • x k+1 -e; ∀(x * , e) ∈ O ξk k+1 , (6)
where O ξk k+1 is the set of optimality cuts cumulated at vertex ξk ∈ Γ k . The statistical average performed over the value associated to each forward pass provides an upper bound of the optimal value of (2).

The backward pass updates the optimality cuts O ξk k+1 on each vertex in the tree using:

inf u k ∈U k c k ( ξk ) • u k + ϑ(k + 1, x k+1 , ξk , O ξk k+1 ) subject to x k+1 = x k + A k u k , ϑ(k + 1, x k+1 , ξ i k , O ξk k+1 ) ≥ x * • x k+1 -e; ∀(x * , e) ∈ O ξk k+1 . (7)
Therefore the backward pass generates a lower bound of the optimal value at the rst stage.

The stopping criteria has been discussed in [START_REF] Shapiro | Analysis of stochastic dual dynamic programming method[END_REF]. The criteria taken in this article is that

|v -v| ≤ ρσ(v) (8)
where v is the forward pass value v is the backward pass value σ(v) is the standard deviation of the forward pass value and ρ is a parameter.

Remark 2.1. Note that we are using the original continuous random process (ξ k ) in the forward pass, whereas the discretized random process ( ξk ) is used in the backward pass. The convergence property follows 3 steps.

The convergence result of L-Shape algorithm on nite probability distribution is given in Birge and Louveaux [7, chapter 7, theorem 1], where the main argument is the niteness of optimality cuts and feasibility cuts. The convergence of classical SDDP method is studied in Philpott and Guan [START_REF] Philpott | On the convergence of stochastic dual dynamic programming and related methods[END_REF] in nite probability framework, and Shapiro [START_REF] Shapiro | Analysis of stochastic dual dynamic programming method[END_REF] in general probability framework.

During the numerical test, we can rst simulate the forward pass following the discretized distribution ( ξk ) and let the forward value and the backward value converge. This convergence follows that the forward value of Monte Carlo method is itself a random variable whose mean is the upper bound of L-Shape method which converges to lower bound by nite convergence of L-Shape method. Remark 2.2. The method presented here deviates from standard stochastic dual dynamic programming in [START_REF] Pereira | Multi-stage stochastic optimization applied to energy planning[END_REF] because the random variable intervenes only in the objective function and is not present in the right hand side of state dynamics.

This causes the partial convexity of Q(t, x t , ξ t ), which is non convex with respect to ξ t , forcing us to proceed to discretization of ξ t . As a result, the optimality cuts cannot be shared among discretized points ξ i t ∈ Γ t . Therefore, this method is clearly outperformed by others [START_REF] Barrera-Esteve | Numerical methods for the pricing of swing options: a stochastic control approach[END_REF][START_REF] Carmona | Valuation of energy storage: an optimal switching approach[END_REF] in terms of CPU time compared when state variable is in low dimension, as it will be shown in section 5.3. However, the main advantage resides in its ability to deal with high dimension (state variable) problems.

Price model

In this paper, we assume that the random process, taking part in the modeling of futures commodity price dynamics, follows the celebrated Black model [START_REF] Black | The pricing of commodity contracts[END_REF].

Although simplistic, this model is equivalent to the one factor Gaussian HJM model (see [START_REF] Jarrow | Bond pricing and the term structure of interest rates : a new methodology for contingent claim valuation[END_REF]) which is well adapted for forward curve modeling. This model is discussed in the practice-based literature, Eydeland and Wolyniec[18, Chapter 5], and Lai et al [START_REF] Lai | An approximate dynamic programming approach to benchmark practice-based heuristics for natural gas storage valuation[END_REF]. Since we are using the vectorial quantization method to discretize the underlying random process, we are not limited to that model.

Other nancial price models can also be applied.

The price model takes as parameters the market volatilities σ, correlations ρ as well as the forward curves F 0 . The Black model states that future prices of commodities are martingales under the risk neutral probability P whose expectations correspond to the original forward curve level F 0 due to lack of arbitrage opportunities. The forward contract price F t s at time s, with maturity (also called tenor) t, follows the dynamic:

ln (F t k+1 ) i (F t k ) i = σ i W i k - 1 2 (σ i ) 2 i = 1, . . . , d; (9)
where

W i k is a standard normal distribution N (0, 1), with correlation corr(W i s , W j s ) = ρ ij .
Then, the spot prices at time t k are simply obtained by taking the t k -expiring forward contracts at time t k :

ξ k = F k k , ξ i k := (F k k ) i = (F k 0 ) i exp k-1 s=0 σ i W i s - 1 2 (σ i ) 2 k i = 1, . . . , d; (10)
where (F k 0 ) is the original forward curve observed. Setting

S i k := exp k-1 s=0 σ i W i s -(σ i ) 2 k/2 , we then have ξ i k = (F k 0 ) i S i k .
Finally, a simple computation on the price model [START_REF] Bonnans | Perturbation analysis of optimization problems[END_REF] shows that:

∂ξ i k ∂(F k 0 ) i = S i k , (11)
∂ξ i k ∂σ i = ξ i k k-1 s=0 W i s -σ i k . ( 12
)
Remark 3.1. In the context of our problem, it is reasonable to consider that

(F, σ) ∈ R d×K + × R d + , where R + = {x > 0} is positive value.
Under this condition, both ξ t , D F,σ ξ t have probability measure absolutely continuous with respect to Lebesgue measure, i.e. the value set does not reduce to singleton. Furthermore, the density function of log normal is bounded. This point is essential to continue analysis below since this property provides continuity even if the integrand is not continuous.

Sensitivity analysis

The principal contribution of this paper lies in the application of Danskin's theorem for computing certain sensitivities of problem [START_REF] Bally | Error analysis of the quantization algorithm for obstacle problems[END_REF]. In this section, we begin by stating this theorem and apply it to our sensitivity analysis. Then, we provide a result on the convergence properties of the sensitivity estimate of our approximation algorithm in section 2 to the actual sensitivities of the original problem (2).

In the subsequent analysis, given two Banach spaces X and Y , we say that φ : X → Y is directionally dierentiable at a point x ∈ X in the direction h ∈ X if the following limit exists:

φ (x; h) := lim t↓0 φ(x + th) -φ(x) t . (13)

Danskin's theorem and its applications

Let us denote by v(F, σ) (resp. U * (F, σ)) the optimal value function (resp.

the optimal solution set) of the original problem [START_REF] Bally | Error analysis of the quantization algorithm for obstacle problems[END_REF]. In this section, we show

that v(F, σ) is directionally dierentiable at every F ∈ R d×K + and σ ∈ R d + , and
Fréchet dierentiable almost everywhere.

We rst study the objective function of problem (2) which can be written as:

f (F, σ, u) = E K-1 k=0 c k (ξ k (F k 0 , σ)) • u k . ( 14
)
f (F, σ, u) is a continuous function on R d×K + ×R d + ×U according to the continuity of c k and the continuity of (ξ k ) in the price model [START_REF] Bonnans | Perturbation analysis of optimization problems[END_REF]. Then, we show the Fréchet dierentiability by the following lemma.

D F,σ f (F, σ, u) = E K-1 k=0 D F,σ c k • u k . ( 15
)
INRIA Proof. In view of the Lipschitz continuity of c k , Rademacher's theorem [START_REF] Rockafellar | Variational analysis[END_REF]Theorem 9.60] states that it is almost everywhere Fréchet dierentiable, and its derivative is bounded. Note that in the derivative of f (F, σ, u) will appear the derivative of c k as an integrand against measure of ξ k . Furthermore, at time t k , the random variable ξ k follows a log-normal distribution according to our price model, whose probability measure is absolutely continuous with respect to the Lebesgue measure. Thus, using chain rule of derivative on the composite mapping as well as the dominated convergence theorem ends the proof.

Following this lemma, we immediately deduce another property of f (F, σ, u).

Proposition 4.2. The derivative

D F,σ f is continuous with respect to (F, σ) ∈ R d×K + × R d
+ , and it is weakly continuous with respect to u. Proof. In view of (15), D F,σ f (F, σ, u) is weakly continuous with respect to u.

Now, let us consider a sequence (F n , σ n ) → (F * , σ * ) in R d×K + × R d + , and denote by ξ n k = ξ k (F n , σ n ) and ξ * k = ξ k (F * , σ * ). We have to show that D F,σ c k (ξ n k ) L 1 --→ D F,σ c k (ξ * k )
. By chain rule of derivative, at the points where c k is dierentiable, we have

D F,σ c k (ξ k ) = D ξ k c k (ξ k )D F,σ ξ k (F, σ). (16)
The term D F,σ ξ k (F, σ) is computed in [START_REF] Brezis | Functional Analysis, Sobolev Spaces and Partial Dierential Equations[END_REF] and [START_REF] Carmona | Valuation of energy storage: an optimal switching approach[END_REF]. As mentioned in remark

c k : R d → R n is Lipschitz, Dc k ∈ L ∞ (Ω, F, P) ⊂ L q (Ω, F, P), 1 ≤ q < ∞.
Since the space of continuous functions with compact support C c (Ω) is dense in L q (Ω, F, P), we have ∀ > 0, ∃g, such that g -

Dc k q ≤ . g(ξ n k ) -Dc k (ξ n k ) q = ξ |g(ξ) -Dc k (ξ)| q pdf n k (ξ)dξ 1/q ≤ pdf n k 1/q ∞ g -Dc k q ≤ pdf n k 1/q ∞ We have same result for (F * , σ * ) g(ξ * k ) -Dc k (ξ * k ) q ≤ pdf * k 1/q ∞ Moreover, ξ n k a.s. --→ ξ * k and D F,σ ξ n k L p --→ D F,σ ξ * k , 1 ≤ p < ∞.
By continuity of g and the dominated convergence theorem, ∃N p , we have that ∀n ≥ N p

g(ξ n k ) -g(ξ * k ) q ≤ . Thus, for n ≥ N q Dc k (ξ n k ) -Dc k (ξ * k ) q ≤ Dc k (ξ n k ) -g(ξ n k ) q + g(ξ n k ) -g(ξ * k ) q + g(ξ * k ) -Dc k (ξ * k ) q ≤ ( pdf n k 1/q ∞ + pdf * k 1/q ∞ + 1) (17) Let → 0, we obtain Dc k (ξ n k ) L q --→ Dc k (ξ * k ). Finally, taking p -1 + q -1 = 1, D F,σ c k (ξ n k ) L 1 --→ D F,σ c k (ξ * k ) follows by combi- nation of convergence of D ξ k c k (ξ k ) in L p and convergence of D F,σ ξ k in L q .
Next, we recall Danskin's theorem: Theorem 4.3. [10, Theorem 4.13] Let V be a nonempty, compact topological space and ψ : R n × V → R. Let φ(x) be the optimal value function and V * (x) be the optimal solution set:

φ(x) := min v∈V ψ(x, v) V * (x) := argmin v∈V ψ(x, v) If ψ(•, v) is dierentiable for every v ∈ V and D x ψ(x, v) is continuous on R n ×V , then φ(x)
is locally Lipschitz continuous, directionally dierentiable, and

φ (x; h) = inf v∈V * (x) D x ψ(x, v)h. (18)
In particular, if for some x ∈ R n , the set V * (x) = {v * } is a singleton, then the min-function is dierentiable at x and

φ (x) = D x ψ(x, v * ). (19)
In order to show the directional dierentiability of v(F, σ), we only need to verify the compactness assumption of the feasible solution space in Danskin's (

) 20
Proof. The almost every Fréchet dierentiability is a direct consequence of Rademacher's theorem [START_REF] Rockafellar | Variational analysis[END_REF]Theorem 9.60] and the fact that v(F, σ) is locally Lipschitz by Danskin's theorem. And if u 1 and u 2 are two optimal solutions such that D F,σ f (F, σ, u 1 ) = D F,σ f (F, σ, u 2 ), then [START_REF] Eydeland | Energy and Power Risk Management: New Developments in Modeling, Pricing and Hedging[END_REF] states that f is not Fréchet dierentiable at (F, σ). Then, we have a contradiction.

The direct application of Danskin's theorem leads to one method to compute the value of sensitivities. We will focus on the sensitivities δ F (dF ) and δ σ (dσ) and give an expression the directional derivative of the optimal value v(F, σ) with respect to the forward curve F and volatility vector σ in the directions dF and dσ:

δ F (dF ) := v (F, σ; dF, 0) = inf (u * k )∈U * E K-1 k=0 D F c k dF k 0 • u * k (21) δ σ (dσ) := v (F, σ; 0, dσ) = inf (u * k )∈U * E K-1 k=0 D σ c k dσ • u * k (22)
Danskin's theorem requires all the optimal solutions which is generally difcult to obtain in stochastic programs which are not strongly convex. In most cases, we can only get one optimal policy (u * k ), and the whole optimal policy INRIA set is hard to compute. Therefore, we can simply obtain an upper bound of the sensitivity value:

δ F (dF ) ≤ δF (dF, u * ) := E K-1 k=0 D F c k dF k 0 • u * k (23) δ σ (dσ) ≤ δσ (dσ, u * ) := E K-1 k=0 D σ c k dσ • u * k (24)
However, by corollary 4.4, v(F, σ) is Fréchet dierentiable almost everywhere. Thus, the inequalities in ( 23) and ( 24) become equality at almost every F and σ. So, knowing one optimal policy u * ∈ U * , allows us to get the derivative value almost everywhere.

Furthermore, we can study the derivatives associated to the standard base of F and σ. At almost every F and σ, we have:

(δ F ) i k := δ F (e i k ) = E ∂c k (ξ k ) ∂(F k 0 ) i • u * k i = 1, . . . , d (25)
(δ σ ) i := δ σ (e i ) = E K-1 k=0 ∂c(ξ k ) ∂σ i • u * k i = 1, . . . , d (26)
Then, at almost every F and σ:

δ F (dF ) = K-1 k=0 d i=1 (dF ) i k (δ F ) i k δ σ (dσ) = d i=1 (dσ) i (δ σ ) i . ( 27
)

Convergence of sensitivity estimate

In this section, we start by presenting two discrete approximations of problem

(2) helping us describe the forward and backward schemes made in our algorithm in section 2. Then we show some properties of the solutions of the discrete approximation problems and present their convergence to the solution of problem [START_REF] Bally | Error analysis of the quantization algorithm for obstacle problems[END_REF]. Eventually, we deduce the convergence of the sensitivity estimates of the discrete approximation problems to the sensitivities of continuous problem.

We consider a sequence of optimal vectorial quantization trees ξm ∈ Γ m , such that the total number of quantized points is N m and lim m→∞ N m = ∞. Denote by ( Fm k ) the canonic ltration associated to the discretized random process ( ξm

k ): Fm k := σ( ξm s , 0 ≤ s ≤ k).
Let us consider two discrete approximations of problem [START_REF] Bally | Error analysis of the quantization algorithm for obstacle problems[END_REF]. The rst one enforces a constant strategy in each Voronoi cell. We call it as the quantization approximation of problem (2):

v m Q (F, σ) = inf (u m k )∈U E K-1 k=0 c k (ξ k ) • u m k subject to u m k ∈ U k , almost surely u m k (ξ k , x m k ) = ξ i,m k ∈Γ m k u i,m k (x m k )1 C i,m k (ξ k ) x m k+1 = x m k + A k u m k , x m 0 = 0 x m T ∈ X m T almost surely (28) RR n°7574
By forward induction, we can deduce that u k is Fm k adapted and that x k is Fm k-1 measurable. This quantization approximation describes the formulation made in the forward pass of our algorithm. Since u k is Fm k adapted, we may write the cost function as a function of ( 28)

as E K-1 k=0 c k (ξ k ) • u m k = E K-1 k=0 E c k (ξ k ) • u m k | Fm k = E K-1 k=0 E c k (ξ k )| Fm k • u m k . ( 29
)
In this formulation the decision u i,m k to be taken when the state belongs to the cell C i k at step t k actually depends on the previous cells along the realization of ξ. We may give another approximation formulation by replacing in each Voronoi

cell E c k (ξ k )| Fm k by c k (E[ξ k | ξ k ∈ C i k ]) in
v m S (F, σ) = inf (u m k )∈U (ξ i,m k ) k p i,m K-1 k=0 c k (ξ i,m k ) • u i,m k subject to u i,m k ∈ U k (non-anticipativity) u i,m k = u j,m k , if ξ i,m [k] = ξ j,m [k] x i,m k+1 = x i,m k + A k u i,m k , x i,m 0 = 0 x i,m K ∈ X i,m K (30)
where ξ i,m In the same way as ( 25) and ( 26), we obtain the sensitivity estimates corresponding to the above discretized approximation problems by Danskin's theo-INRIA rem: The quantization approximation problem (28) only diers from (30) by re-

(δ m Q,F ) i k := (v m Q ) (F, σ; e i k , 0) ≤ ( δm Q,F ) i k := E ∂c k (ξ k ) ∂(F k 0 ) i • (u * ,m Q ) k i = 1, . . . , d (31)
(δ m Q,σ ) i := (v m Q ) (F, σ; 0, e i ) ≤ ( δm Q,σ ) i := E K-1 k=0 ∂c(ξ k ) ∂σ i • (u * ,m Q ) k i = 1, . . . , d (32)
(δ m S,F ) i k := (v m S ) (F, σ; e i k , 0) ≤ ( δm S,F ) i k := E ∂c k ( ξk ) ∂(F k 0 ) i • (u * ,m S ) k i = 1, . . . , d (33)
(δ m S,F ) i := (v m S ) (F, σ; e i , 0) ≤ ( δm S,σ ) i := E K-1 k=0 ∂c( ξk ) ∂σ i • (u * ,m S ) k i = 1, . . . , d (34)
placing the coecient in objective function E[c(ξ t ) | ξ t ∈ C i t ] by c(E[ξ t | ξ t ∈ C i t ]
). Thus, concavity of value function v m Q with respect to the coecient in the objective function still holds. With the same arguments, we get the almost everywhere Fréchet dierentiability for the optimal value function of the quantization approximation problem (28).

Therefore, we are able to obtain exact sensitivity values for both discrete approximation problems following Danskin's theorem almost everywhere, i.e.

the inequalities (31)-(34) become equalities almost everywhere.

Lemma 4.6. A subsequence of the optimal strategy u * ,m Q of the quantization approximation problem (28) converges weakly in L 2 (Ω, (F t ), P; R n×K ) to an optimal strategy u * of continuous problem [START_REF] Bally | Error analysis of the quantization algorithm for obstacle problems[END_REF].

Proof. We denote by v * Q the limit inferior of the associated value function

v * Q := lim inf m→∞ v * ,m
Q , and the subsequence associated to this limit inferior u * ,mn1 Q . The feasible set of ( 28) is a subset of U taking value in U which is non empty, convex and bounded. The optimal solution u * ,mn1 Q is therefore a bounded sequence in a Hilbert space U . Therefore, we deduce the existence of a weakly convergent subsequence u * ,mn2 In order to obtain the reverse inequality, let us build a sequence which is feasible for the quantization approximation problem (28) and which converges to an optimal solution of problem [START_REF] Bally | Error analysis of the quantization algorithm for obstacle problems[END_REF]. The feasible solution set U is included in L 2 , where C ∞ is dense. Then for one optimal solution (u * k ) ∈ U * of the continuous problem ( 2), there exists one sequence

(u n,c k ) ∈ C ∞ converging to (u * k ) such that K-1 k=0 u n,c k -u * ≤ 1/n. We build u n,m k (ξ [k] ) = E u n,c k (ξ [k] ) | Fm k ,
which is one feasible solution of the quantization approximation problem (28).

The convergence of u n,m k to u n,c k in L 2 is due to the convergence ξ-ξm → 0 and the dominated convergence theorem. We denote by v n,m Q the value associated

with strategy u n,m k . Therefore, v m Q ≤ v n,m Q . Let us consider the sequence v n,n Q ,
and let n → ∞, we obtain the second inequality v * Q ≤ v. Thus, we get v * Q = v and hence prove the lemma.

One immediate corollary is Corollary 4.7. The sensitivity estimates of the quantization approximation problem δm Q,F and δm Q,σ dened in (31) and (32) converge to the sensitivity values of continuous problem δ F and δ σ (dened in [START_REF] Hull | Options, Futures and other derivatives[END_REF] and ( 22)) at almost every

F ∈ R d×K + and σ ∈ R d + .
Lemma 4.8. A subsequence of the optimal strategy u * S of scenario approximation problem (30) converges weakly in L 2 (Ω, (F t ), P; R n×K ) to an optimal strategy u * of continuous problem [START_REF] Bally | Error analysis of the quantization algorithm for obstacle problems[END_REF].

Proof. In view of lemma 4.6, we just need to compare the optimal value function v m S and v m Q . Both discretized approximation problems have the same feasible solution space L 2 (Ω, ( Fm k ), P) → (U k ). By a straightforward computation, we can obtain that

v m Q -v m S = K-1 k=0 E[c k (ξ k )(u * ,m Q ) k ] - K-1 k=0 E[c k ( ξk )(u * ,m S ) k ] ≤ K-1 k=0 E[c k (ξ k )(u * ,m S ) k ] - K-1 k=0 E[c k ( ξk )(u * ,m S ) k ] ≤ K-1 k=0 [c k ] lip U k ∞ ξ t -ξt 2 ≤ C QS K-1 k=0 ξ k -ξk 2 ≤ 2C QS ξ -ξ 2 INRIA
where U ∞ = sup {|u| ∈ U}, and C QS is a constant which only depends on [c k ] lip and U ∞ . This error bound can be obtained on (v m S -v m Q ) by using same technique. Thus we have

|v m Q -v m S | ≤ 2C QS ξ -ξ 2
Following Zador's theorem (see below theorem 4.10), we get that lim m→∞ v m Q = lim m→∞ v m S . Then, the lemma can be proved by following the same argument as in the proof of 4.6.

Then, we have a similar corollary as corollary 4.7 for the scenario approximation problem: Corollary 4.9. The sensitivity estimates of the scenario approximation problem δm S,F and δm S,σ dened in (33) (34) converge to the sensitivity values of continuous problem δ F and δ σ dened in [START_REF] Hull | Options, Futures and other derivatives[END_REF] and [START_REF] Pereira | Multi-stage stochastic optimization applied to energy planning[END_REF], at almost every

F ∈ R d×K + and σ ∈ R d + .
Proof. This corollary is obtained following the fact that ( 5 Algorithm and numerical tests

Algorithm

After getting one optimal ( -optimal) strategy (u * ,tr-sddp k ) thanks to our discretization and SDDP algorithm described in section 2, we use the Monte Carlo method to compute the expectations in (31)-(32).

We simulate a sample ξ s , s = 1, . . . , S of trajectories and compute the value of sensitivities for each trajectory:

s (δ F ) i k = ∂c k (ξ s k ) ∂(F k 0 ) i • s u * ,tr-sddp k , i = 1, . . . , d; (35)
s (δ σ ) i = K-1 k=0 ∂c k (ξ s k ) ∂σ i • s u * ,tr-sddp k , i = 1, . . . , d. (36)
Therefore we obtain the mean and standard deviation of this sample:

µ tr-sddp (δ F ) i k = 1 S S s=1 s (δ F ) i k , σ tr-sddp (δ F ) i k = 1 S -1 S s=1 ( s (δ F ) i k -µ tr-sddp (δ F ) i k ) 2 1/2 ; (37)
µ tr-sddp

(δσ) i = 1 S S s=1 s (δ σ ) i , σ tr-sddp (δσ) i = 1 S -1 S s=1 ( s (δ σ ) i -µ tr-sddp (δσ) i ) 2 1/2 . (38)

Comparison of methods

We shall now address the question of the accuracy of the results obtained from the application of Danskin's theorem. In that perspective, we benchmark them against the sensitivities obtained by the mean of other algorithms.

The rst comparison lies in the discretization of state variable (x k ). Indeed, many popular algorithms impose the discretization of (x k ) as a prerequisite whilst our SDDP based algorithm removes this constraint. Although more ecient than SDDP in low dimensions, the (x k ) discrete algorithms become numerically intractable when the dimension of (x k ) increases.

The second source of comparison concerns the discretization of the random process (ξ k ). In our framework, we have chosen the optimal quantization tree method as a way to discretize the randomness. As long as (ξ k ) remains in low dimension, we may also implement a PDE based method.

Finally, the last and probably most important criteria of comparison resides in the way the sensitivity values are computed. We challenge the numerical results obtained by the application of the Danskin's theorem against the nite dierence method which is as simple as it is popular.

Danskin + quantization tree + state space discretization

The rst comparison method is the algorithm proposed in [START_REF] Bardou | Optimal quantization for the pricing of swing options[END_REF] that combines the quantization tree method and the dynamic programming principle on discretized state space (x k ). Let us stress again that this method proves to be very ecient as long as the state variable (x k ) remains in low dimension. Let v tr-disc stands for the Bellman value and u * ,tr-disc for the optimal strategy:

v tr-disc (k, xk , ξk ) = min u k ∈U k c k ( ξk )u k + E[v tr-disc (k + 1, xk + u k , ξk+1 ) | ξk ] ; (39) u * ,tr-disc (k, xk , ξk ) = argmin u k ∈U k c k ( ξk )u k + E[v tr-disc (k + 1, xk + u k , ξk+1 ) | ξk ] ; (40)
where

E[v tr-disc (k + 1, x, ξk+1 ) | ξk = ξ i k ] = ξ j k+1 ∈Γ k+1 p ij k v tr-disc (k + 1, x, ξ j k+1 ).
We can reuse the same formula as ( 35)-(38) with the optimal policy u * ,tr-disc to compute the sensitivity values that will be denoted by (δ F ) tr-disc and (δ σ ) tr-disc .

INRIA

Thanks to this procedure, we are able to compare two methods to compute the optimal solutions obtained on the same quantization tree as the one used in our SDDP based algorithm.

Danskin + PDE + state space discretization

If the random process (ξ k ) is in low dimension, one natural idea is to use a PDE method to compute conditional expectations. Therefore, the second comparison method presented in this paragraph combines a PDE method with the dynamic programming principle on a discretized state space (x k ). In order to use a PDE method, we need to build one continuous time process from the discrete time random process. To be consistent with our price model ( 9), we will use the continuous version of the Black model:

dF • t = σF • t dW t ; (41)
and make an interpolation of F t 0 between t k and t k+1 :

F t 0 = t -t k t k+1 -t k F k+1 0 + t -t k+1 t k -t k+1 F k 0 , t k ≤ t ≤ t k+1 . (42)
Let v pde-disc (k, x k , ξ k ) stand for the Bellman value, and u * ,pde-disc (k, x k , ξ k )

for the optimal strategy:

v pde-disc (k, xk , ξ k ) = min u k ∈U k c k (ξ k )u k + E[v pde-disc (k + 1, xk + u k , ξ k+1 ) | ξ k ] ; (43)
u * ,pde-disc (k, xk , ξ k ) = argmin

u∈U k c k (ξ k )u k + E[v pde-disc (k + 1, xk + u k , ξ k+1 ) | ξ k ] . (44)
Let ṽ(t, x,

ξ t ) := E[v pde-disc (k + 1, x, ξ k+1 ) | ξ t ], t k ≤ t ≤ t k+1 . Combining
the Feynman-Kac formula and Itô's lemma on the continuous extension of our price model (41) (42), we compute the conditional expectations by solving the following PDE backwards from t k+1 to t k :

   ∂ṽ ∂t + 1 2 d i=1 (σ i ) 2 (ξ i ) 2 ∂ 2 ṽ ∂(ξ i ) 2 + i =j ρ ij σ i σ j ξ i ξ j ∂ 2 ṽ ∂ξ i ∂ξ j = 0, t k ≤ t ≤ t k+1 ; ṽ(t k+1 , •, ξ t k+1 ) = v pde-disc (t k+1 , •, ξ t k+1 ). ( 45
)
This parabolic partial dierential equation is straightforward to solve numerically in dimension 1. For example, we may solve it by the use of a nite dierence method with an implicit scheme, which is introduced in appendix A.

Once we have obtained the optimal solution u * ,pde-disc , we may still apply Danskin's theorem by using the same formula (35)-(38). We denote the sensitivity result by (δ F ) pde-disc and (δ σ ) pde-disc . Note that this method does not depend on the quantization tree. Hence, the results obtained from the method presented in the current paragraph give evidence of the accuracy of discretizing (ξ k ) by the quantization tree method.

Finite dierence (fd.) + PDE + state space discretizaton

The last method used to compare the accuracy of the sensitivity estimates is probably the most used in practice. Once an ecient numerical method to compute the optimal value v(F, σ) is available, the nite dierence method may be used to approximate the derivative value:

(δ f d F ) i k = v(F + e i k , σ) -v(F -e i k , σ ) 2 ;
(46)

(δ f d σ ) i = v(F, σ + e i ) -v(F, σ -e i ) 2 .
(47)

In the following numerical tests, we use PDE and the state variable discretization to compute the optimal value v(F, σ). We denote by (δ F ) f d-pde and (δ σ ) f d-pde the sensitivity estimates. This method allows us to compare the relevance of the values obtained from Danskin's theorem.

Swing option

The rst test is performed on swing options which are fairly common in the energy markets. A swing option allows its holder to purchase a total amount of commodity during a predetermined period of time. A typical swing option contract has the following parameters: a feasible control set U k = [0, 1] for every time step t k , a nal state set X K = [L, U ], and a time horizon K = 50 in our example. The running cost function c(ξ k ) is equal to ξ k -K k , where the strike

K k = F k 0 .
Under certain circumstances, swing option featuring one dimension control possesses the bang-bang property. Bardou et al [START_REF] Bardou | Optimal quantization for the pricing of swing options[END_REF][START_REF] Bardou | When are swing options bang-bang and how to use it?[END_REF] show that there always exists a bang-bang optimal strategy when (L, U ) ∈ Z 2 . Therefore, the discretization on x t takes values in Z only.

In the current test, we take the lower bound of the exercise right L equal to 20, and upper bound of the exercise right U equal to 30. We set F k 0 in price model curve to 1.0 + 0.2 sin(2πt k /T ), and volatility σ to 30%/ √ K. The quantization tree built contains N = 5000 vertices dispatched on K = 50 stages.

The estimation of transition probability is carried out by Monte-Carlo method using 10 7 samples. The SDDP algorithm generates 5000 scenarios in forward pass and uses 20 scenarios in backward pass. The stopping criteria (8) takes ρ = 1. Concerning the parameters intervening in the PDE resolution (see section A), we set the lower bound ξ to 0.1, the upper bound ξ to 5.0, the number of discretization space steps J to 200, and nally the number of time discretization L to 5. Finally, we take = 0.01 in the nite dierence estimation. The program is written in C++ with GNU Linear Programming Kit(GLPK) 4.44, and the test was performed on a PC with a 3.07GHz Intel Xeon CPU and 12GByte main memory. The results obtained for optimal values are given in table

INRIA

The forward value (by Monte Carlo) is a random variable whose expectation is a value associated to a suboptimal policy, i.e. an upper bound of optimal value. Furthermore, our backward value is a lower bound of discretized problem, which is not the lower bound of the original problem. Therefore, it is possible that the forward value is smaller than the backward value.

The sensitivity values with respect to each contract of the forward curve are displayed in gure 1. Finally, the sensitivity values with respect to σ are given in table .

Small commodity portfolio case study

We now consider a gas trading portfolio. A trading company purchases natural gas from a set of producing countries indexed at a price formula and sells it to consuming countries at another other price formula (see gure 2 for the main market and table 3 for the price formulae). Annual quantity and price formulae have been agreed contractually, the latter are functions of the future prices of major energy markets ξ t , such as crude oil price (OIL), north American natural gas price (NA NG), and Europe natural gas price (EU NG). The objective is to analyse the sensitivity of the optimal value of the portfolio with respect to the forward price 4 and volatility of the random process. * ξ 1 = OIL, ξ 2 = EU NG, and ξ 3 = NA NG.

The quantization tree is processed with N = 36000 R 3 -valued elementary quantizers dispatched on T = 12 stages. The estimation of transition probability is carried out by Monte-Carlo method using 10 9 samples. In this numerical test, the other comparison methods are not adapted because of the large dimension of the state variable m = 6. Discretizing both state variable and the random space will make the problem numerically intractable. Therefore, we can only provide results based on Danskin's theorem combined with SDDP algorithm in section 2. The optimal value is presented in table 5 A Appendix Implicit scheme of nite dierence for 1 dimension PDE Generally, the parabolic partial dierential equation is very dicult to solve numerically in high dimension. We refer to Allaire [START_REF] Allaire | Numerical analysis and optimization: An introduction to mathematical modelling and numerical simulation[END_REF] and Ern and Guermand [START_REF] Ern | Theory and practice of nite elements[END_REF] for more details.

In the case that ξ t is in dimension one, the equation ( 45) is written as    ∂ṽ ∂t + 1 2 σ 2 ξ 2 ∂ 2 ṽ ∂ξ 2 = 0, t k ≤ t ≤ t k+1 , ξ ∈ [ξ, ξ], ṽ(t k+1 , •, ξ t k+1 ) = g(t k+1 , •, ξ t k+1 ).

(48)

To apply the nite dierence scheme, we need to discretize both time t and space ξ. Let the time step be h t = (t k+1 -t k )/L, and the space step be h ξ = (ξ -ξ)/J, where L, J are positive integers. We denote by ṽl j an approximation of ṽ(t k + lh t , ξ + jh ξ ) where 0 ≤ l ≤ L and 0 ≤ j ≤ J. (51)

INRIAFinally,

  we can change the forward simulation to following continuous distribution. Using continuous distribution in forward sampling gives an indication on the error induced by the discretization of the continuous distribution.

Lemma 4 . 1 .

 41 The objective function f (F, σ, u) is Fréchet dierentiable with respect to F ∈ R d×K + and σ ∈ R d , and its derivative is

  3.1, ξ * k and its derivative D F,σ ξ * k (resp. ξ n k and D F,σ ξ n k ) both have bounded density function pdf * k (resp. pdf n k ) with respect to Lebesgue measure. Since

theorem 4 . 3 .

 43 In our framework, the feasible solution set U taking value in (U k ) is a bounded closed subset of L 2 ∩ L ∞ . From the Banach Alaoglu Bourbaki theorem[11, theorem 3.16], it is a weak compact subset in L 2 . Corollary 4.4. The optimal function v(F, σ) is Fréchet dierentiable at almost every F ∈ R d×K + and σ ∈ R d + , and at the point where v(F, σ) is dierentiable, we have Dv(F, σ) = D F,σ f (F, σ, u), for all u ∈ U * .

  [k] := (ξ i0,m 0 , . . . , ξ i k ,m k ) ∈ (Γ m 0 , . . . , Γ m k ) is a sub-trajectory in the quantization tree until time t k , p i,m := K-2 k=0 p i k i k+1 ,m k is the probability associated to one trajectory ξ i,m [K-1] . This scenario approximation describes the formulation made in the backward pass in our algorithm. Denote by v m Q (F, σ) and u * ,m Q (F, σ) (resp. v m S (F, σ) and u * ,m S (F, σ)) the optimal value function and one optimal solution of the quantization approximation problem (28) (resp. of the scenario approximation problem (30)).

Lemma 4 . 5 .

 45 If the distribution of ξ k is absolutely continuous with respect to the Lebesgue measure, then the optimal value functions of the two discrete discrete approximation problems v m S (F, σ) and v m Q (F, σ) are Fréchet dierentiable at almost every F ∈ R d×K and σ ∈ R d + . Proof. The scenario approximation problem (30) can be viewed as a large linear program. According to the linear programming theory, v m S is piecewise linear, concave with respect to the coecient of the objective function dened in whole space, so v m S is locally Lipschitz and then almost everywhere Fréchet dieren- tiable by Rademacher's theorem [24, Theorem 9.60]. Under the assumption on the distribution of ξ k , and the fact that (F, σ) → ξm (F, σ) is injective and Fréchet dierentiable, the chain rule of Fréchet derivative of composite mapping holds almost everywhere and we deduce the lemma for scenario approximation problem (30).

Q

  that converges in the weak topology to u * Q following [11, Corollary 3.30]. Since the feasible set is closed and convex, it is weakly closed and hence u * Q is feasible. Since u * ,m Q is always feasible for the original problem (2), we get v ≤ v * ,m Q for every m. Let m go to ∞ and we obtain the rst inequality v ≤ v * Q .

  ξm k ) converges to (ξ k ) and that u * ,m S converges weakly in U to u * and the [11, Proposition 3.13]. Finally, we recall Zador's theorem: Theorem 4.10. [20, Theorem 6.2] If E[|ξ| p+η ] < ∞ for some η > 0, then, lim N N p/d min |Γ|≤N ξξ p p = J p,d |g| d/(d+p) (u)du 1+p/d where P(du) = g(u)λ d (du) + ν, ν ⊥ λ d (λ d Lebesgue measure on R d ). The constant J p,d corresponds to the case of the uniform distribution on [0, 1] d . This theorem implies that min |Γ|≤N ξξ p = O(N 1/d ).

Figure 1 :

 1 Figure 1: Sensitivity values with respect to F obtained by the 4 methods. Since the standard deviations of the rst 3 methods represent less than 1% of the sensitivity values, we will not plot them for sake of clarity.

Figure 2 :

 2 Figure 2: A ctive supply and demand portfolio, as well as the possible routes. : producing country; •: consuming country.

  The parameters driving the stochastic processes in the test are set as follows: volatility σ 1 = σ 2 = σ 3 = 40%, and correlations ρ 12 = 0.7 ρ 13 = 0.2 ρ 23 = 0.4. And the parameters driving the SDDP part of the algorithm are forward pass samples M f = 3000, backward pass sample M b = 10. The stopping criteria (8) takes ρ = 1. The program is written in C++ with Cplex 10.1, and the tests were performed on a PC with a 2.2 GHz Dual Core AMD CPU and 16GByte main memory.

  ) = g(t k+1 , •, ξ + jh ξ ).

∂ξ 2

 2 need the boundary condition at bounds of ξ. By the arguments of the problem context, we take the following boundary conditions: = 0, ξ = ξ, ∂ 2 ṽ ∂ξ 2 = 0, ξ = ξ.

Table 1 :

 1 Optimal values for the 3 methods.

	1.

Table 2 :

 2 2. Sensitivity value with respect to σ obtained by the 4 methods.

	tree+sddp.	tree+discret.	pde.+discret.	fd.+ pde.
	mean	std. dev.	mean	std. dev.	mean	std. dev.
	δ σ -32.181	0.695	-32.198	0.698	-32.207	0.700	-35.558
	We remark that the sensitivity values obtained by all 4 methods are close to
	each other. This result proves the accuracy of the sensitivity values obtained by
	our original algorithm (Danskin + quantization tree + SDDP). Furthermore,
	we also observe the closeness of the deltas obtained by nite dierences to those
	obtained by Danskin for the case where the discretization methods are the same
	(i.e. PDE+discretization).
	Compared to other methods, our algorithm does not procure any advantage
	in CPU time for such a low dimension problem. In order to show the advantage

Table 3 :

 3 Constraints and price formulaeThe forward price F t 0 in (9) is read from the energy market. We take the following values in the numerical test (see in table4).

	Time	OIL *	EU NG *	NA NG *
	0	73.0	5.2	5.15
	1	73.5	5.3	5.30
	2	74.0	5.2	5.22
	3	74.5	4.5	4.51
	4	75.0	4.4	4.35
	5	75.5	4.3	4.33
	6	76.0	4.2	4.18
	7	76.5	4.3	4.32
	8	77.0	4.4	4.37
	9	77.5	5.2	5.21
	10	78.0	5.3	5.25
	11	78.5	5.4	5.42

* The quantity unit is in TBtu. MMBtu stands for a million British thermal unit, a TBtu is a trillion British thermal unit thus equivalent to a Million MMBtu. ** The price unit is $/MMBtu.

INRIA

Table 4 :

 4 Forward price

Table 6 :

 6 , and sensitivity values are in tables 6 and 7. sensitivity value with respect to forward price * σ here means standard derivation.

		upper bound	std. dev. of upper bound	lower bound
	optimal value	-39.867	2.101	-41.900
	Table 5: optimal value by quantization tree + SDDP method

Table 7 :

 7 sensitivity values with respect to volatility

RR n°7574

http://www.inria.fr ISSN 0249-6399