Implementation of a Scheduling and Pricing Model for Natural Gas

Since 1999, the Australian state of Victoria has operated a natural gas spot market to both determine daily prices for natural gas and develop an optimal schedule for the market based on an LP (Linear Programming) approximation to the underlying inter-temporal nonlinear aspects of the gas flow optimization problem. This market employs a dispatch optimization model and a related market clearing model. Here we present the model employed for both the operational scheduling and price determination. The basic dispatch optimization formulation covers the key physical relationships between pressure, flow, storage, with flow controlled by valves, and assisted by compressors, where flow and storage are measured with respect to energy rather than in terms of mass. But we also discuss a range of sophisticated mathematical techniques which have had to be employed to create a practical dispatch tool, including iterating between piecewise and successive linearization; iterating between barrier and simplex algorithms to manage numerical accuracy and solution speed issues, and special methods developed to deal with scheduling flexibility. The market clearing model is a variation on the dispatch optimization model which replaces the gas network with an infinite storage tank with unlimited transport capacity. We address the performance of the model including accuracy and run time.

1. 1.

   The use of common carriage is due to the similarities with electricity markets. Victoria actually use the term “market carriage” to represent its pool based market, contrasting with “contract carriage” which relates to more traditional scheduling under bilateral contracts for access to each individual transmission pipelines.

2. 2.

   As of the Victorian gas market recently coming under Australia’ National Gas Rules the PTS is now called the DTS, or Determined Transmission System.

3. 3.

   Based on information supplied by AEMO for winters 2007 to 2009 inclusive

4. 4.

   Storage of LNG is limited, and it is possible to use a significant proportion of the storage quite quickly if LNG is used too freely. Peak LNG vaporization for 1 h can take about 1 day to replace at the (emergency) maximum rate of LNG production, though more typically takes several days to replace.

5. 5.

   In this formulation, variables are generally upper case, and constants lower case, whereas Read et al. [1] use lower case for variables, and upper case for constants. Upper and lower limits are still represented by over and under bars, respectively. Unless otherwise stated, all variables in this formulation are positive.

6. 6.

   slack is a variable, yet it is also a vector, and as such it is stated in lower case.

7. 7.

   Although this terminology is common in such models these are not really “slack” variables in the traditional sense. They do indicate how far the final solution point is from the constraint, but it lies outside the feasible region, not inside. Thus they have the form of generalized slack variables, which are positive when the constraint is slack, but negative when it is violated.

8. 8.

   Market participants submit supply offers to the market, these ‘bids’ are composed of up to 10 price and quantity tranches, where we call each individual tranche a ‘step’.

9. 9.

   The “slack” variables applied to such constraints, have very small penalties, so as to encourage, but not force, tie-breaking.

10. 10.

    Parameters associated with a number of these constraints limits/bounds may also be scaled in a pre-processing step so as to resolve conflicts between quantities of gas previously scheduled and quntities actually observed, for example.

11. 11.

    PS indicates the number of “pressure states”, although actually there are PS + 1, including state 0.

12. 12.

    For a 2 dimensional representation like this, every point can actually be defined as a weighted sum of three particular grid points. But, in order to define a convex feasible region, we need to leave the model to determine which grid points it prefers to use.

13. 13.

    Together, these two constraints actually mean all weights must lie in the range [0,1].

14. 14.

    Another potential advantage of the convex combination approach to piece-wise linearization is that, if and when convexity issues arise, they could be dealt with by employing an integer type formulation, employing Special Ordered Sets (SOS2) as in Martin et al. [14], to force the model to apply weights only to adjacent grid points. A similar strategy is mentioned by DeWolf and Smeers (2000a), but dismissed on computational efficiency grounds. Still, that approach has been applied successfully in modeling HVDC link losses in the New Zealand electricity market model.

15. 15.

    Early experience with the model showed that most cases solved within 14 iterations. Also, cases which did not solve showed negligible improvement in convergence after 14 successive approximations.

16. 16.

    In practice, it is usual for the axes to be reversed to match compressor manufactures data relating to efficiency curves, however we present them in this orientation for consistency

17. 17.

    There are actually two null operational states one for the minimum inlet pressure, and the other for the maximum outlet pressure, although they appear as one point in Fig. 6

18. 18.

    In theory, operating states between zero and the minimum operating level could actually be achieved, on average, by operating efficiently for only part of the time. If so, the convex approximation to performance in that region could actually be valid. Hydro generators, for example, can be validly operated, and represented, in that way. In this case, though, the savings in operational efficiency would probably be outweighed by increased startup/shutdown costs.

19. 19.

    This lower limit here defines the minimum running boundary of the feasible region if we force the compressor to be “on”, thus eliminating the null states. As shown, this boundary is actually non-convex, but it can reasonably be approximated by a straight line, at the cost of eliminating a relatively small set of operating states that are seldom utilized, in practice.

20. 20.

    In other words, the pressure can decline significantly without implying the increase in flow rate one would otherwise expect for such a decline.

21. 21.

    This being important because it determines who faces penalties and/or receives constrained-on/off payments when the actual dispatch schedule differs from the market trading schedule.

22. 22.

    Energy = Work = Force × Displacement and Pressure = Force÷Area; so combining these two results in Pressure × Length = Force × Displacement÷Area which then equals the flow of energy past a point, which is Work÷Area = Energy÷Area

Authors and Affiliations

1. ICF International, Fairfax, Virginia, USA

   W. Pepper

2. Market Reform, Melbourne, Australia

   B. J. Ring

3. University of Canterbury, Christchurch, New Zealand

   E. G. Read & S. R. Starkey

Corresponding author

Correspondence to E. G. Read .

Editors and Affiliations

1. , Industrial and Systems Engineering, University of Florida, Weil Hall 303, Gainesville, 32611, Florida, USA

   Alexey Sorokin

2. , Division of Economics and Business, Colorado School of Mines, 15th Street 816, Golden, 80401, Colorado, USA

   Steffen Rebennack

3. , Department of Industrial & Systems Engin, University of Florida, Weil Hall 401, Gainesville, 32611-6595, Florida, USA

   Panos M. Pardalos

4. Research & Development, EnerCoRD - Energy Consulting,, Plastira Street 4, Athens, 171 21, Greece

   Niko A. Iliadis

5. Rio Praia de Botafogo, Centro Empresarial, -A-Botafogo 2281701, Rio de Janeiro, 22250-040, Rio de Janeiro, Brazil

   Mario V. F. Pereira

Appendix: Detailed Pressure Flow Equations for Flow Rates and Linepack

In this Appendix we present the key equations for deriving natural gas flows and pressures in the Victorian pipe network. The derivation is based on six initial equations described in Eqs. 16, 17, 18, 19, 20, 21) below.

The ideal gas law equation 1 describes pressure, P _l at a point l along a pipeline of length 0 < l < L. R is the ideal gas constant in units of kPa*m³/(°K × kg), T is the temperature of the pipeline in °K, and z is the supercompressibility of gas, while ρ _l is the density in units of kg/m³.

The gas flow rate q _l at a point l, in kg/s, is described in (17). Where A is the pipe cross-sectional area in units of m², ν _l is the velocity of gas in m/s in the pipe at point l, at the point l measured at sea level.

The Fanning equation 18 describes the rate of change of pressure at position l along a pipeline. D is the pipe diameter in units of m while f is the Fanning friction factor.

The Reynolds Number, Re, is defined by (19). In this equation μ is the viscosity of gas measured in kg/ms.

The Fanning friction factor formulation in (20) utilizes the Blasius formulation [18], where η is the pipeline efficiency (a fraction). It describes friction as gas flows along a pipe. This makes use of a turbulent friction factor, fturb.

Smooth pipelines have no turbulence while rough pipelines have more turbulence. Turbulence becomes more important as the Reynolds Number increases. Given e, a measure of the roughness of the pipeline in mm, then the relative roughness of a pipeline can be defined with respect to its diameter as (e/D). Using empirical data the form of the turbulent friction factor shown in (21) was developed by gas system engineers working on the development of the MCE. Now, using Eq. 17 to substitute for ν _l in (18) we get:

Equation 22 can be further refined by using Eq. 16:

Where:

By integrating (23) with respect to l and observing that for l = 0, P _l = P _o, the origin pressure of the pipeline gives:

Using (25) to define the value of P _l where l = L, i.e. the destination pressure, and assuming a constant flow (q _l = q), friction factor, and supercompressibility, then the steady state flow rate can be derived as:

Here Eq. 26 is closely related to the Weymouth panhandle equation referred to by Zheng et al. [8] and Midthun et al. [13]. It is used later to define friction factors as a function of flow rate. However, we can derive another flow rate equation by assuming no friction arises from turbulence (fturb = 0), and substituting the Reynolds Number from (20) into (19) and the resultant equation for f into (22) to give:

Where:

Integrating (27) with the assumption that the flow is constant along the pipe segment and the requirement that for l = 0, P _l = P _o (the origin pressure) we get:

Using (29) to define the value of P _l where l = L, i.e. destination pressure, then for a non-constant friction factor and constant supercompressibility the steady state flow rate can be derived as:

Equations 6 and 30 describe flow rates under different assumptions about friction factors. These are used later to define more general flow rate equations. However, before exploring that, it is necessary to consider the linepack equations. The linepack, I, in a pipe segment can be calculated by integrating the volume of gas in each slice of pipeline along the length of the pipeline:

Assuming a constant friction factor this can be rewritten as:

We also have^{Footnote 22}:

Substituting this into (32), we derive:

Substituting the expression for P _l from (25), which was based on a constant fraction factor, into (33), then integrating over l, we get, for a non-zero flow rate:

Given the rule (f(x))ⁿ differentiated by x gives n × df/dx × f(x)^n-1 then we can integrate dF(l) by the reverse transformation to give:

Evaluation of (36) for F(0) and for F(L), and substituting these into (34) gives:

Further, substituting for q from (26) gives linepack for non-zero flow of:

Where the flow rate is zero, then P _l = P ₀ for all l, and (16) implies

Equations 38 and 39 can be represented in terms of an average pressure on the pipeline P _a:

Hence

Or:

Given user defined values of typical low and high pressures in the system, P _low and P _high, and corresponding supercompressibility values z _low and z _high it is possible to define an average supercompressability z _a as:

To take advantage of the linear relationship between linepack and pressure in (40), the MCE formulation uses this equation to compute linepack for all cases, including the case when the pipeline inlet and outlet pressure values are different. The term (P ₀ /z) is replaced by the average value of the supercompressibility-adjusted-pressures at the inlet and outlet nodes. Further adjustments are made to these equations to allow for altitude. The model also combines the results of (26) and (30) to determine flows which address the impact of a varying friction factor, non-zero fturb, and varying supercompressibility along the pipeline, for given pressures at the origin and destination of the pipe. All of these values are refined using an iterative approach that converges quickly and tests have demonstrated that a further iteration past the current stopping point would typically impact final flows by less than 0.2%.

Reprints and permissions

© 2012 Springer-Verlag Berlin Heidelberg

Policies and ethics