GLPK/Electricity markets

This tutorial contains a set of examples demonstrating the features of competitive electricity markets used for market-based dispatch. Although the Australian National Electricity Market (NEM) is used as the reference for these examples, this market shares many characteristics with other competitive electricity markets and can offer useful insights into marked-based dispatch more generally.

Single-period electricity dispatch may be framed as a flow network problem — with each region (and later, each transmission line) represented by a node and each loss-free link represented by a bi-directional arc. The examples presented here are instead constructed as linear programs. The two formulations are mathematically equivalent (although the preferred algorithms normally differ). Nonetheless, it is useful to retain a network view in mind as you study these examples.

Simple single unit dispatch edit

The purpose of this example is to demonstrate the dispatch of a single generation unit for a defined regional price over a series of time periods, taking account of fuel costs and ramp rates, and maximising the total profit for the unit.

/*

 Simple single unit dispatch

 Dr. H J Mackenzie
 2010-03-24

*/

# set of dispatch intervals

set I;

# dispatch price is $

param regionalprice {I};

# unit characteristics

param unit_max_capacity >= 0;
param fuel_cost >= 0;
param max_ramp_rate >= 0;
param start_dispatch >= 0;

# dispatch variables

var dispatch {I} >= 0;
var ramp {I}, >= - max_ramp_rate, <= max_ramp_rate;
var profit {I};

# objective function

maximize totalprofit: sum {i in I} profit[i];

# constraints

s.t. initial_dispatch: dispatch[0] = start_dispatch;
s.t. dispatch_profit {i in I}: profit[i] = 
 dispatch[i] * (regionalprice[i] - fuel_cost);
s.t. dispatch_ramp {i in I}: ramp[i] = 
 dispatch[i] - dispatch[if i > 0 then i-1 else 0];
s.t. unit_capacity {i in I}: dispatch[i] <= 
 unit_max_capacity;

# solve the problem

solve;

# output input and determined values

printf {i in I} "%d regionalprice = %.1f; dispatch = %.1f; ramp = %.1f; profit = %.1f\n", 
i, regionalprice[i], dispatch[i], ramp[i], profit[i];

data;

param unit_max_capacity := 100; /* MW */
param fuel_cost := 30;      /* $/MWh */
param max_ramp_rate := 20;    /* max MW up or down */
param start_dispatch := 0;    /* present dispatch point MW */

param : I :  regionalprice :=
     0   15
     1   17
     2   18
     3   22
     4   55
     5   40
     6   65
     7   10
     8   12
     9   4
;

end;

The solution is found in about 10 iterations and shows the dispatch of the unit prior to profitable running to maximise the total profit of $2680 for the entire set of time periods. Unfortunately in the real world, the regional price is not known prior to the dispatch interval and so the approach is only useful for the analysis of optimal dispatch after the dispatch has occurred or, in the case of pre-dispatch pricing, with the implicit assumption that the dispatch of the unit is insignificant with respect to the regional price.

*  10: obj = 2.680000000e+003 infeas = 0.000e+000 (0)
OPTIMAL SOLUTION FOUND
Time used:  0.0 secs
Memory used: 0.1 Mb (135015 bytes)
0 regionalprice = 15.0; dispatch = 0.0; ramp = 0.0; profit = 0.0
1 regionalprice = 17.0; dispatch = 0.0; ramp = 0.0; profit = 0.0
2 regionalprice = 18.0; dispatch = 20.0; ramp = 20.0; profit = -240.0
3 regionalprice = 22.0; dispatch = 40.0; ramp = 20.0; profit = -320.0
4 regionalprice = 55.0; dispatch = 60.0; ramp = 20.0; profit = 1500.0
5 regionalprice = 40.0; dispatch = 80.0; ramp = 20.0; profit = 800.0
6 regionalprice = 65.0; dispatch = 60.0; ramp = -20.0; profit = 2100.0
7 regionalprice = 10.0; dispatch = 40.0; ramp = -20.0; profit = -800.0
8 regionalprice = 12.0; dispatch = 20.0; ramp = -20.0; profit = -360.0
9 regionalprice = 4.0; dispatch = 0.0; ramp = -20.0; profit = 0.0

Single region energy dispatch with regulation edit

This example is similar to the previous example but introduces the complication of a market dispatched ancillary service for system regulation. In the real Australian NEM, there are six market dispatched services for raises of 6 seconds (R6SEC), 60 seconds (R60SEC), 5 minutes (R5MIN) and regulation (RREG) and for lowers of 6 seconds (L6SEC), 60 seconds (L60SEC), 5 minutes (L5MIN) and regulation (LREG). In this example, only the RREG service is considered for a single region where the market cost for the region is minimised. The example does not consider real world constraints regarding the dispatch of ancillary services where varying amounts of service are available at different output levels, implemented in the Australian market as a defined trapezium for each ancillary service and the physical characteristics of each unit.

/*

 Least cost unit dispatch for energy with regulation

 Dr. H J Mackenzie
 2010-03-24

*/

# set of units

set U;

# maximum capacity of the units in MW

param capacity {U};

# offer price of the energy in $/MWh

param offerprice {U};

# maximum regulation capacity of the units in MW

param regcapacity {U};

# regulation offer price of the regulation energy in $/MWh

param regofferprice {U};

# market parameters

param dispatch_demand >= 0;
param regdispatch_demand >= 0;

# dispatch variables

var dispatch {U} >= 0;
var cost {U};
var regdispatch {U} >= 0;
var regcost {U};

# objective function

minimize market_cost: sum {u in U} (cost[u] + regcost[u]);

# constraints

s.t. dispatch_cost {u in U}: cost[u] = dispatch[u] * offerprice[u];
s.t. regdispatch_cost {u in U}: regcost[u] = regdispatch[u] * regofferprice[u];
s.t. feasible_dispatch {u in U}: dispatch[u] + regdispatch[u] <= capacity[u];
s.t. feasible_regdispatch {u in U}: regdispatch[u] <= regcapacity[u];
s.t. dispatch_matches_demand: sum {u in U} dispatch[u] = dispatch_demand;
s.t. regdispatch_matches_demand: sum {u in U} regdispatch[u] = regdispatch_demand;

# solve the problem

solve;

# output input and determined values

printf {u in U} "%d capacity = %.1f; offerprice = %.1f; regcapacity = %.1f; regofferprice = %.1f\n", 
 u, capacity[u], offerprice[u], regcapacity[u], regofferprice[u];
printf {u in U} "%d dispatch = %.1f; cost = %.1f; regdispatch = %.1f; regcost = %.1f;\n", 
 u, dispatch[u], cost[u], regdispatch[u], regcost[u];

data;

param dispatch_demand := 45;   /* MW */
param regdispatch_demand := 10; /* MW */

param : U :  capacity  offerprice  regcapacity  regofferprice:=
     1   10     10      5       1
     2   10     15      5       7
     3   10     20      5      13
     4   10     25      5      19
     5   10     30      5      25
     6   10     35      5      31
     7   10     40      5      37
     8   10     45      5      43
     9   10     50      5      49
    10   10     55      5      55
;

end;

The results show that the (electrical) energy dispatch for the region is influenced by the provision of the ancillary service therefore leads to a higher total market cost than if just energy were considered. Note that the additional of the regulation service has lead to a dispatch that is not obvious if only examining the energy dispatch.

*  13: obj = 1.090000000e+003 infeas = 0.000e+000 (0)
OPTIMAL SOLUTION FOUND
Time used:  0.0 secs
Memory used: 0.1 Mb (137757 bytes)
1 capacity = 10.0; offerprice = 10.0; regcapacity = 5.0; regofferprice = 1.0
2 capacity = 10.0; offerprice = 15.0; regcapacity = 5.0; regofferprice = 7.0
3 capacity = 10.0; offerprice = 20.0; regcapacity = 5.0; regofferprice = 13.0
4 capacity = 10.0; offerprice = 25.0; regcapacity = 5.0; regofferprice = 19.0
5 capacity = 10.0; offerprice = 30.0; regcapacity = 5.0; regofferprice = 25.0
6 capacity = 10.0; offerprice = 35.0; regcapacity = 5.0; regofferprice = 31.0
7 capacity = 10.0; offerprice = 40.0; regcapacity = 5.0; regofferprice = 37.0
8 capacity = 10.0; offerprice = 45.0; regcapacity = 5.0; regofferprice = 43.0
9 capacity = 10.0; offerprice = 50.0; regcapacity = 5.0; regofferprice = 49.0
10 capacity = 10.0; offerprice = 55.0; regcapacity = 5.0; regofferprice = 55.0
1 dispatch = 5.0; cost = 50.0; regdispatch = 5.0; regcost = 5.0;
2 dispatch = 5.0; cost = 75.0; regdispatch = 5.0; regcost = 35.0;
3 dispatch = 10.0; cost = 200.0; regdispatch = 0.0; regcost = 0.0;
4 dispatch = 10.0; cost = 250.0; regdispatch = 0.0; regcost = 0.0;
5 dispatch = 10.0; cost = 300.0; regdispatch = 0.0; regcost = 0.0;
6 dispatch = 5.0; cost = 175.0; regdispatch = 0.0; regcost = 0.0;
7 dispatch = 0.0; cost = 0.0; regdispatch = 0.0; regcost = 0.0;
8 dispatch = 0.0; cost = 0.0; regdispatch = 0.0; regcost = 0.0;
9 dispatch = 0.0; cost = 0.0; regdispatch = 0.0; regcost = 0.0;
10 dispatch = 0.0; cost = 0.0; regdispatch = 0.0; regcost = 0.0;

Two region energy dispatch edit

In this example we look at the dispatch of two regions given varying capacity at different price bands in each region and different regional demand for a single time period. The objective function seeks to find the lowest cost solution for both regions. For this example it is assumed that there are no link losses between each region. Energy in this context refers to electrical energy.

/*

 Least cost unit dispatch for energy with two regions
 
 Dr. H J Mackenzie
 HARD software
 hjm@hardsoftware.com

 2010-03-25

*/

# set of units

set U;

# maximum capacity of the units in MW for each region a and b

param regiona_capacity {U};
param regionb_capacity {U};

# offer price of the energy in $/MWh

param offerprice {U};

# market parameters

param regiona_dispatch_demand >= 0;
param regionb_dispatch_demand >= 0;
param ab_link_capacity >= 0;

# dispatch variables

var regiona_dispatch {U} >= 0;
var regionb_dispatch {U} >= 0;
var unit_dispatch {U} >= 0;
var ab_link_dispatch; /* positive flow is defined as flow from region a to region b */
var cost {U};

# objective function

minimize market_cost: sum {u in U} cost[u];

# constraints

s.t. dispatch_cost {u in U}: cost[u] = ((regiona_dispatch[u] + ab_link_dispatch) 
 * offerprice[u]) + ((regionb_dispatch[u] - ab_link_dispatch) * offerprice[u]);
s.t. feasible_regiona_dispatch {u in U}: regiona_dispatch[u] <= regiona_capacity[u];
s.t. feasible_regionb_dispatch {u in U}: regionb_dispatch[u] <= regionb_capacity[u];
s.t. total_unit_dispatch {u in U}: unit_dispatch[u] = regiona_dispatch[u] + regionb_dispatch[u];
s.t. feasible_ab_link_dispatch: ab_link_dispatch <= ab_link_capacity;
s.t. feasible_ba_link_dispatch: -ab_link_dispatch <= ab_link_capacity;
s.t. regiona_dispatch_matches_demand: sum {u in U} regiona_dispatch[u] 
 - ab_link_dispatch = regiona_dispatch_demand;
s.t. regionb_dispatch_matches_demand: sum {u in U} regionb_dispatch[u] 
 + ab_link_dispatch = regionb_dispatch_demand;

# solve the problem

solve;

# output input and determined values

printf {u in U} "%d regiona_capacity = %.1f; regionb_capacity = %.1f; offerprice = %.1f\n", 
 u, regiona_capacity[u], regionb_capacity[u], offerprice[u];


printf {u in U} "%d dispatch = %.1f; ab link dispatch = %.1f; cost = %.1f\n", 
 u, unit_dispatch[u], ab_link_dispatch, cost[u];

data;

param ab_link_capacity := 10;     /* MW */

param regiona_dispatch_demand := 9;  /* MW */
param regionb_dispatch_demand := 22; /* MW */

param : U :  regiona_capacity  regionb_capacity  offerprice:=
     1   10         0         10
     2   10         0         20
     3   10         0         30
     4   0        10         15
     5   0        10         25
     6   0        10         35
;

end;

The solution shows the dispatch in each region, with the link running at full capacity due to the lower demand and lower cost of generation in region A.

*   5: obj = 4.800000000e+002 infeas = 0.000e+000 (0)
OPTIMAL SOLUTION FOUND
Time used:  0.0 secs
Memory used: 0.1 Mb (140601 bytes)
1 regiona_capacity = 10.0; regionb_capacity = 0.0; offerprice = 10.0
2 regiona_capacity = 10.0; regionb_capacity = 0.0; offerprice = 20.0
3 regiona_capacity = 10.0; regionb_capacity = 0.0; offerprice = 30.0
4 regiona_capacity = 0.0; regionb_capacity = 10.0; offerprice = 15.0
5 regiona_capacity = 0.0; regionb_capacity = 10.0; offerprice = 25.0
6 regiona_capacity = 0.0; regionb_capacity = 10.0; offerprice = 35.0
1 dispatch = 10.0; ab link dispatch = 10.0; cost = 100.0
2 dispatch = 9.0; ab link dispatch = 10.0; cost = 180.0
3 dispatch = 0.0; ab link dispatch = 10.0; cost = 0.0
4 dispatch = 10.0; ab link dispatch = 10.0; cost = 150.0
5 dispatch = 2.0; ab link dispatch = 10.0; cost = 50.0
6 dispatch = 0.0; ab link dispatch = 10.0; cost = 0.0

Modeling of quadratic line losses edit

Losses in transmission lines are quadratic of the form            and therefore need to be approximated with a linear form for solving with an LP formulation. The way to approach the problem is to approximate the function            with the formulation:   subject to   This is an example of separable programming being used to approximate a convex non-linear function (see HP Williams ‘Model Building in Mathematical Programming’[1] – Chapter on non-linear models). Separable functions are functions which can be broken down into single variable functions, and these in turn can be approximated by piecewise linear functions.


 

To do:
• upload graph image file of losses and linear approximation


For the GLPK formulation, I have modeled both positive and negative flows over the link to check that it works correctly for flow in both directions. It is very easy to confuse the signs and to get this wrong, so checking both flows is a good check for the technique. This technique will be required for all of the inter-connector flows between regions and so will need to be correct for later more sophisticated models.

/*

 Linearised loss model

 Dr. H J Mackenzie
 2010-04-13
 
 Problem description
 
  Losses on the link are loss = 0.01 x flow^2 modeled 
  as a 10 step linear approximation [-10, 10]
  
  Modeling of loss for loss = k x flow^2
  
  k =     0.01
  steps =   10
  min flow = -10
  max flow = 10
  
  interval  flow  losses  slope  midpoint  mp losses
        10   1  
  1      -8   0.64   -0.18  -9     0.81
  2      -6   0.36   -0.14  -7     0.49
  3      -4   0.16   -0.1   -5     0.25
  4      -2   0.04   -0.06  -3     0.09
  5      0   0     -0.02  -1     0.01
  6      2   0.04    0.02   1     0.01
  7      4   0.16    0.06   3     0.09
  8      6   0.36    0.1   5     0.25
  9      8   0.64    0.14   7     0.49
  10     10   1     0.18   9     0.81
*/

# set of regions

set REGIONS;

# set of transmission lines

set LINES;

param line_start {LINES} symbolic in REGIONS;
param line_end {LINES} symbolic in REGIONS;
param max_flow {LINES} >= 0;
param min_flow {LINES} <= 0;

param no_loss_segments;
set LOSS_SEGMENTS := 1.. no_loss_segments;

param loss_segment_length {LINES, LOSS_SEGMENTS} >= 0;
param loss_at_minimum_flow {LINES} >= 0;
param loss_coefficient {LINES, LOSS_SEGMENTS}; # can be positive or negative

# dispatch variables

param midline_flow {l in LINES} >= min_flow[l], <= max_flow[l];

var line_loss {LINES};
var loss_segment_dispatch {l in LINES, s in LOSS_SEGMENTS} >= 0 , <= loss_segment_length[l, s];
var line_start_flows {LINES};
var line_end_flows {LINES};

# objective function

minimize line_losses: sum {l in LINES} line_loss[l];

# constraints

s.t. transmission_line_flow {l in LINES}: sum {s in LOSS_SEGMENTS} loss_segment_dispatch[l, s] 
 + min_flow[l] = midline_flow[l];
s.t. transmission_line_loss {l in LINES}: sum {s in LOSS_SEGMENTS} loss_segment_dispatch[l, s] 
 * loss_coefficient[l, s] + loss_at_minimum_flow[l] = line_loss[l];
s.t. transmission_line_start_flows {l in LINES}: midline_flow[l] 
 + 0.5 * line_loss[l] = line_start_flows[l];
s.t. transmission_line_end_flows {l in LINES}: midline_flow[l] 
 - 0.5 * line_loss[l] = line_end_flows[l];

# solve the problem

solve;

# output input and determined values

printf "\n\nLine losses\n\n%-10s %-10s %10s %10s %10s\n", 'Line', 'Segment', 'Length', 'Coeff', 'Dispatch';
printf {l in LINES, s in LOSS_SEGMENTS} "%-10s %-10d %10d %10.2f %10.2f\n", 
 l, s, loss_segment_length[l, s], loss_coefficient[l, s], loss_segment_dispatch[l, s];

printf "\n\nLine dispatch\n\n%-10s %-10s %-10s %10s %10s %12s %10s %10s %10s %10s\n", 
 'Line', 'Line start', 'Line end', 'Max flow', 'Min flow', 'Minflow loss', 'Line flow', 
 'Line start', 'Line end', 'Line loss';
printf {l in LINES} "%-10s %-10s %-10s %10.2f %10.2f %12.2f %10.2f %10.2f %10.2f %10.2f\n", 
 l, line_start[l], line_end[l], max_flow[l], min_flow[l], loss_at_minimum_flow[l], 
 midline_flow[l], line_start_flows[l], line_end_flows[l], line_loss[l];

printf "\n";

data;

set REGIONS := 'REGION_A', 'REGION_B';

param : LINES : line_start line_end  max_flow min_flow loss_at_minimum_flow midline_flow :=
'AtoBneg'    'REGION_A' 'REGION_B' 10    -10    1.0          -6
'AtoBpos'    'REGION_A' 'REGION_B' 10    -10    1.0          6
;

param no_loss_segments := 10;

param   : loss_segment_length loss_coefficient :=
'AtoBneg' 1   2          -0.18
'AtoBneg' 2   2          -0.14
'AtoBneg' 3   2          -0.1
'AtoBneg' 4   2          -0.06
'AtoBneg' 5   2          -0.02
'AtoBneg' 6   2           0.02
'AtoBneg' 7   2           0.06
'AtoBneg' 8   2           0.1
'AtoBneg' 9   2           0.14
'AtoBneg' 10  2           0.18
'AtoBpos' 1   2          -0.18
'AtoBpos' 2   2          -0.14
'AtoBpos' 3   2          -0.1
'AtoBpos' 4   2          -0.06
'AtoBpos' 5   2          -0.02
'AtoBpos' 6   2           0.02
'AtoBpos' 7   2           0.06
'AtoBpos' 8   2           0.1
'AtoBpos' 9   2           0.14
'AtoBpos' 10  2           0.18
;
end;

The solution works consistently for both flow directions and gives reasonable results for a nominal flow measured at the midpoint of the line.

*  10: obj = 7.200000000e-001 infeas = 0.000e+000 (0)
OPTIMAL SOLUTION FOUND
Time used:  0.0 secs
Memory used: 0.1 Mb (141776 bytes)


Line losses

Line    Segment    Length   Coeff  Dispatch
AtoBneg  1          2   -0.18    2.00
AtoBneg  2          2   -0.14    2.00
AtoBneg  3          2   -0.10    0.00
AtoBneg  4          2   -0.06    0.00
AtoBneg  5          2   -0.02    0.00
AtoBneg  6          2    0.02    0.00
AtoBneg  7          2    0.06    0.00
AtoBneg  8          2    0.10    0.00
AtoBneg  9          2    0.14    0.00
AtoBneg  10         2    0.18    0.00
AtoBpos  1          2   -0.18    2.00
AtoBpos  2          2   -0.14    2.00
AtoBpos  3          2   -0.10    2.00
AtoBpos  4          2   -0.06    2.00
AtoBpos  5          2   -0.02    2.00
AtoBpos  6          2    0.02    2.00
AtoBpos  7          2    0.06    2.00
AtoBpos  8          2    0.10    2.00
AtoBpos  9          2    0.14    0.00
AtoBpos  10         2    0.18    0.00


Line dispatch

Line    Line start Line end   Max flow  Min flow Minflow loss Line flow Line start  Line end Line loss
AtoBneg  REGION_A  REGION_B    10.00   -10.00     1.00   -6.00   -5.82   -6.18    0.36
AtoBpos  REGION_A  REGION_B    10.00   -10.00     1.00    6.00    6.18    5.82    0.36

References edit

  1. Williams, H. Paul (1999). Model Building in Mathematical Programming. Wiley. ISBN 978-0-471-99788-7.