# 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.

Regarding comments and improvements, please contact Dr. Harley Mackenzie [HARD software] at hjm@hardsoftware.com. I welcome your input.

## Contents

## 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 HARD software hjm@hardsoftware.com 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 HARD software hjm@hardsoftware.com 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.

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 HARD software hjm@hardsoftware.com 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]

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