A stochastic optimization model for natural gas sale companies

1. A stochastic optimization model for natural gas sale companies M.T. Vespucci (*), F. Maggioni (*), E. Allevi (#), M. Bertocchi (*), M. Innorta (*) (#) University of Brescia (*) University of Bergamo

2. Structure of presentation • description of some principles of gas market liberalization • details of problem to be solved • deterministic model • stochastic model

3. Basic principles of liberalized gas market • unbundling of production and transportation activities both at national and local levels • national level: • shippers: production, import, re-gasification, wholesale commercialisation • national distributor: transport on national network and storage • local level: • local distributors: transport on local networks • gas sellers: purchase gas from shippers and sell it to final consumers • protection of so called “small final consumers” realized by Regulatory Authority by setting a maximum price they may be required to pay

4. Small final consumers Consumers whose annual consumption does not exceed 8 · 106 MJoule domestic customers (cooking, cooking/heating) commercial activities, crafts and small industries medium and large industries classes 1- 6 : high consumption proportion is for heating (depends on weather conditions) industrial customers: consumption for production (independent on weather conditions)

5. Maximum price for classes 1 - 6 set by Regulatory Authority Gas maximum price set (and periodically revised) by Regulatory Authority on the basis of following splitting of cost QE + QVI + QL + QT + QS + TD + QF + QVD where QE : raw material cost QVI : wholesale commercialization cost QL : costs of rigassification of liquid gas QT : trasportation cost QS : storage cost TD : distribution cost QF : fixed retail commercialization cost QVD : variable retail commercialization cost shipper costs national distributor costs local distributor cost gas seller costs

6. Shipper – gas seller interaction • Gas seller purchases gas from shipper on the basis of a contract: • one contract for each citygate operated by gas seller • for each thermal year (July 1st – June 30th) • In the contract are indicated • gas volume required by gas seller for next thermal year (Va_citygate) • gas volume required in particular in winter months (Vw_citygate) • maximum daily consumption (capacity) requested by gas seller (Cg_citygate) • purchase price fixed by shipper (P) • In the contract it is also specified how to compute penalties, • which are due by gas seller if daily consumption exceeds daily capacity.

7. Gas seller problem • Gas seller needs a model for determining optimal decisions about • 1. number of final customers to supply in each consumption class • (so called citygate customer portfolio) • 2. sell prices to apply to each consumption class • Different customer portfolios determine different citygate consumption patterns. • In a citygate where mainly industrial customers are served • citygate consumption tends to be constant along the year, • since it is used more for production than for heating • it is easier for gas seller to determine the citygate daily capacity to require • Shippers prefer a more constant citygate consumption along the year, •  a lower purchase price for gas seller is set by shipper

8. In a citygate where mainly “small” customers are served • (consumption strongly dependent on weather conditions) • citygate consumption may have big fluctuations in the year • daily consumption may exceed daily capacity • a higher purchase price for gas seller is set by shipper • penalties are likely to be paid by gas seller • Revenue side: • gas seller can fix “freely” only prices for industrial consumers: • these prices have to be set so as to maximize gas seller profits • at the same time the price proposed to the industrial customer • does not have to result in the customer buying gas from another gas seller

9. Purchase price and citygate consumption profile • How does shipper fix price P paid by gas seller ? • a mathematical relation is not known • historical data show that better citygate consumption profiles • correspond to lower prices P • Consumption profile indicators: • ratio among winter consumption and annual consumption constant consumption in all months • ratio of average daily gas consumption of gas seller and daily • capacity Cg_citygate (citygate loading factor)

10. Shipper consumption preferences • citygate consumption as constant as possible along the year more than two thirds of gas consumption concentrated in winter months • average daily use of “virtual pipeline” as high as possible less than half of daily capacity used on average Linear regression model of purchase price P onto LF_citygate: P = QT + QS +intercept + slopeP · LF_citygate Model based on both LF_citygate and _citygate not significant, _citygate and LF_citygate being highly correlated.

11. Computation of penalties If daily consumption in month i exceeds daily capacity, penalties are applied by shipper to gas seller. Percentages ik and unitary penalties ik (Euro/m3) for computing penalties are set in the contract. Example: citygate Sotto il Monte (thermal year 2003 - ’04) intervals with different unitary penalties are numbered from 0 to K

12. Gas seller commercial policies: sell prices • protected customers (classes 1 – 6): maximum price, qvj , fixed by Regulatory Authority with possible discount, sj, 0 sj < 1, fixed by gas seller • industrial customers (classes 7 – 10) Pj” : fixed by gas seller, who relates selling price for class j to • gas purchase price, P, paid to shipper • consumption profile of customers of class j

13. Prices for industrial customers: customer consumption profiles • On the basis of historical data of average monthly consumption • per customer of class j in month i (Vm_customerij) • 2 indicators of consumption profiles of customer j are computed • average monthly consumption of customer j | | | | | | gas consumption from November to March gas consumption of thermal year 2) average daily consumption of customer j • > 1coefficient set by Regulatory Authority | | | | | n° of days in month i 

14. Example: citygate Sotto il Monte - thermal year 2003-’04 Vm_customerij : average monthly consumption per customer of class j in month i

15. Example:  ratio and loading factor per consumer class in citygate Sotto il Monte (thermal year 2003-04) Citygate gas demand depends on gas demand of each consumption class weigthed by the number of customers in each class (ncj) • First 6 consumption classes: • their consumption behaviours have a negative impact on citygate • consumption profile • therefore gas seller would tend to apply high price to counterbalance • the high purchase price due to bad citygate profile • but regulations protect these customers by setting a price cap

16. Last 4 consumption classes: • their consumption behaviours have a positive impact on citygate consumption profile • gas seller should determine a sell price Pj” that attracts these customers, • in order to pay a lesser purchase price P to shipper • and therefore increase profits, in particular from first 6 classes positive term subtracted if LF_customerj > LF_citygate (better than) positive term subtracted if _customerj < _citygate (better than) e.g.: coeff1 =2.126, coeff2 =2.554, rechargej = 1.1

17. Data of deterministic model • Vm_customerij : average monthly consumption per customer of class j • in month i of previous thermal year (historical data) • : maximum number of customers available belonging to class j • QT, QS, , qvj for 1 j 6 : set by Regulatory Authority • intercept, slopeP : determined by regression on historical data • sj for 1 j 6, rechargej for 7 j 10, coeff1 and coeff2 : set by gas seller • ik and ik (widths of intervals from 0 to K–1 for computing penalties)

18. Consumption estimates • Assumption: • average monthly consumption of each class may be estimated • by previous year consumption • Therefore, following estimates are used for next thermal year: • annual consumption of customer of class j • peak consumption of customer of class j (for penalty computation)

19. The deterministic optimization model • Find values of • number of customers per class • ncj , 1  j  10 (integer) • contract parameters • Vm_citygatei , 1  i  12 • Va_citygate,Vw_citygateand Cg_citygate • variables related to penalties computation • surplusik , 5  i  9 and 0  k  2 • that maximize gas seller profits

20. Objective function: gas seller profits revenues from classes 1 – 6 revenues from classes 7 – 10 costs penalties

21. subject to • constraints defining monthly, winter and annual citygate consumption (citygate demand in month i) (citygate annual demand) (citygate winter demand) • constraints for computing penalties • lower and upper bounds on number of customers per class, restricted to be integer

22. The gas seller citygate model is a nonlinear mixed integer model with linear constraints nonlinearities, coming from definitions of LF_citygate and _citygate, appear only in the objective function. Simulation framework: based on ACCESS 97, for database management MATLAB, release 12, for data visualization GAMS, release 21.5, for optimisation Optimization solver: DICOPT (in GAMS framework) solves a sequence of NLP subproblems, by CONOPT2, and MIP subproblems, by CPLEX

23. Model validation • analyse the impact on overall citygate management of industrial consumer • (_customer7 = 0.463 and LF_customer7= 0.628) marginal profits give indications about possible further reduction of price for first 6 classes through parameter sj

24. Model validation • investigate commercial policies towards customers of consumption class 6 • (whose annual consumption is 73% of total annual consumption of first 6 classes)

25. Dependence on temperature of consumption of classes 1 – 6 • Gas consumption of first 6 classes strongly depends • on temperature variations along months • Investigate impact of weather conditions on optimal gas seller decisions • build scenarios of future temperatures • stochastic version of the model • numerical experiments

26. Building scenarios of future temperatures Data: minimum and maximum daily temperature (degree Celsius) measured in Bergamo from 1.1.1994 to 30.11.2005 Tt Observing historical data: temperature is a mean reverting process, reverting to some cyclical function t

27. Histogram of daily temperature differences • temperature differences • between 2 subsequent days • approximate a normal distribution • temperature process to be • modeled as a Brownian Motion dWt : Wiener process at : speed of mean reversion mean value, which the process reverts to t : process volatility Tt: process to be modelled

28. Deterministic model of temperature • (phase angle): max and min temperatures do not necessarily occur at January 1st and July 1st it models cyclic behaviour in the year global warming trend assumed to be linear By using addition formulas for sin function, a linear model in unknown parameters a1, a2, a3 and a4 is obtained

29. Values of a1, a2, a3 and a4 that correspond to following parameter values in model of temperature A = 13.33 B = 6.8891 · 10-5 C = 10.366  = – 1.7302

30. Estimation of volatility t We only need a value of volatility for each month  t is taken as a piece-wise function,constant during each month initi : number of the day in the year at which month i begins mean of squared differences between temperature values of two subsequent days

31. Estimation of speed of reversion at Efficient estimator of at is (Bibby and Sorensen, 1995)

32. Process simulation From the differential equation we obtain the approximation scheme t, 1  t  365, independent standard normally distributed random variables by which nscen temperature scenarios for 365 days ahead are built

33. Representation of scenarios Scenarios are then expressed in Heating Degree Days Tts max {18° – Tts , 0} mean temperature in monthi, as available consumption data refer to months expected value of random variable Tmis over scenarios deviation of mean temperature in month i from mean value over scenarios

34. Monthly consumption, dependent on temperature, of first 6 classes random variable slopeSij : consumption variation for unitary temperature variation in month i for class j • Representation of scenarios: for each month i • interval between is divided into nint sub-intervals • with middle point Dir being the representative point of sub-interval • (r : sub-interval index, 1 r nint) • to temperature differenceDir we associate • – monthly gas consumption • – probability PRir based on frequency • (n° of scenarios of month i belonging to sub-interval r)

35. Objective function of stochastic model Expected value of profits, given by revenues from classes 1 – 6 revenues from classes 7 – 10 costs penalties

36. Expected value of profits • By computing the expected value of profits we obtain that • terms representing expected values of revenues from classes 1 – 6 and 7 – 10: • coincide with corresponding terms of deterministic model • terms representing expected values of costs and of penalties: • computed under the hypothesis that • random variables Dis and Djs , ij , are independent C P where

37. Model validation • sensitivity of solution to generation of different numbers of scenarios • (a) 1000 scenarios • (b) 10000 scenarios • given a number of scenarios, how solution changes as representations • of a given number of scenarios become more and more refined • i.e. number of sub-intervals increases: nint = 5, 10, 15, …… • Note: • both in (a) and in (b), • as scenario representation becomes more and more refined (nint increases), • optimal profit converges to a value between 151 700 and 151 710, • Such value is lower than the value (154 265) obtained in the deterministic case, • as expected, because cold scenarios have been taken into account

38. Solution of stochastic model: 1000 scenarios  5

39. Solution of stochastic model: 10000 scenarios  5

40. Computed solutions Case study: citygate Sotto il Monte (thermal year 2004-05) optimal values of deterministic model optimal values of stochastic model • in the deterministic solution no penalties are paid • solution of stochastic model: by requiring in the contract 4 484 725 m3 gas • (expected value of annual volume), gas seller has an expected value of profit • lower than in the deterministic case, since the expected value of penalties is • positive. This solution, though, allows gas seller to have the same purchase price • of the deterministic case and therefore the same selling price for the industrial • customer.

41. Future work • As number of scenarios and number of sub-intervals increase, • complexity of the problem also increases (numerical difficulties) •  devise a new algorithm that decouples computation of Cg_citygate • from all other decision variables, so that reduced problem is linear • Scenarios not represented on a monthly basis, but treated as vectors of 12 • realizations  scenario reduction techniques • There exists a relation between purchase price P and international price indeces, • since gas seller must choose the index of reference among a certain number of • admitted choises (e. g. oil price index, etc.) • Study the influence on P of future variations of these indices to help gas seller • in taking his decision.