Dispatch, Pricing and Power Flows

Disclaimer Reasonable care has been taken to ensure that the information in this paper is up to date at the time of issue. Potential users of EMarket should, however, ensure that they evaluate EMarket and this paper through an appropriate evaluation process in consultation with Energy Link. The authors are also reliant on certain information external to EMarket and Energy Link, for which no liability or responsibility can be accepted.

Introduction

This technical bulletin is intended to provide users and interested parties with a detailed explanation of how EMarket’s IQP dispatch and pricing module functions, including the DC power flow calculations embedded within it.

Dispatch

Dispatch and pricing in EMarket and EMO are primarily determined by the methodology specified in Schedule 1 of Part 2 of the New Zealand Electricity Governance Rules (EGRs). The market implementation of these methodologies is currently Trans Power's SPD model which is a large scale LP.

Reserves modelling is required within EMarket, primarily due to their impact on generation. The reserve constraints are straight forward to model in principle, however, a major concern with reserves is the large number of reserve constraints which combine to significantly increase the processing time, often for little or no additional accuracy.

RiskOffsets are also required for the modelling of reserves. These are calculated for SPD using TransPower's RMT tool which is a complex non linear dynamic model of the generation and load on the Grid. RMT is too complex to implement within EMarket.

In EMarket and EMO dispatch, pricing and power flow are achieved by Energy Link’s IQP software module.

Optimum Dispatch

The dispatch objective function from Schedule G6 of Part G of the EGRs is to maximise the net purchaser surplus.

EGR Dispatch Objective

where the D_i are the demand bids and the BP_i the bid prices, G_j the dispatched generation with offers OP_j, R_k the dispatched reserves with offers OR_k.

This objective function is used to form the forecast prices and the pre-dispatch schedule using the demand bids BP_i, generator offers OP_j and reserve offers OR_k. But provisional and final prices are calculated using metered demand, so in effect it is performed by minimising the Cost function shown below. Note that the offers Ox have units of $/MWh whereas G and R have units of MW, so Cost has units of $/hour.

Pricing Objective Function

Both objective functions are subject to the constraints listed in the following table (This is not a complete list of all of the constraints within SPD). Note that j indexes generators, i indexes either lines or island, depending on context, n indexes nodes and c indexes class of reserve.

Constraint	Formulation
Generator output greater than zero	G_j ≥ 0 ∀j
Generator output less than maximum output	G_j - G_j,max ≤ 0 ∀j
Generation plus reserve less than generator's maximum output	G_j +R_j - G_{j, max} ≤ 0 ∀j
Reserve less than a specified input, rslope_j, proportion of generation	R_j - rslope_jG_j ≤ 0 ∀j
Nodal Energy Balance Constraints: Energy conserved at each and every node. Note that the Z_in and Z_out include any HVDC flows in or out of node
AC line flow, Z_i, within thermal limits which are usually the same in both directions	Z_i - Z_{i, max} ≤ 0 ∀j Z_{i, min} - Z_i ≤ 0 ∀j
HVDC MW within limits but note maximum transfer different in each direction	As for AC lines
Losses match power flows	Handled by DC power flow formulation
Island Reserve Constraints: For each island and class of reserve (6s and 60s) total reserves greater than or equal to island risk S_i
For each island and class of reserve (6s and 60s) island risk greater than or equal to island risk adjustment factor times maximum generator unit output less RiskOffset	S_{i, c} ≥ IRAF_{i, c} max (G_j - RiskOffset)
For each island and class of reserve (6s and 60s) island risk greater than or equal to island risk adjustment factor times HVDC MW into island less HVDC overload factor Z_DCOL for the island	S_{i, c} ≥ IRAF_{i, c} (Z_DC - Z_DCOL)
For each island and class of reserve (6s and 60s) island risk greater than or equal to island risk adjustment factor times island minimum risk	S_{i, c} ≥ IRAF_{i, c} S_{min, i, c}
Security constraints - we understand these are not usually in use	n/a
Ramping – ability of generators to achieve new output within a half hour	n/a

The calculation of the island reserve risk in EMarket includes the minimum island risk entered by the user, so hydro risk in the South Island is covered.

The IQP Dispatch Algorithm

IQP, our iterated quadratic program, was developed solely by Energy Link to find optimal dispatches of generation and reserve. Quadratic programming is a widely used optimisation technique and avoids many of the pitfalls that can occur in more general problems, particularly if the objective function is convex.

Formulation of the Optimisation Problem

IQP is used to find an optimum solution for dispatching generation on the grid based on a set of generation and reserve offers, fixed nodal demand, and a large set of constraints including line limits and reserves constraints.

The value of the objective function will henceforth be referred to as the Cost of a dispatch, or Cost(d) and the object of a dispatch to minimise Cost(d).

The generation must equal the demand across the grid including losses. This condition can be expressed as the Global Conservation of Energy Constraint.

Global Conservation of Energy Constraint

where D_Tot is the figure for total demand and L is the total losses on the grid (L includes the losses on the Grid between GIPs and GXPs. It does not include distribution losses that are incurred downstream of each GXP).

Losses can be calculated from net nodal injection, where net injection at a node is the difference between injection and off-take at the node. IQP uses the standard DC power flow formulation where one node is designated the swing node or swing bus. Net nodal generation at the swing node balances the total net generation of all other nodes to satisfy the global conservation of energy constraint.

Losses in a DC Power Flow

: i,j = 0 to N - 1

where N is the total number of nodes and the y_i are the net nodal injections, i.e. y_i is the total generation at node i minus total demand at node i minus distributed losses. B is the susceptance matrix for the Grid, where the susceptance is approximately the reciprocal of the line reactance.

Since y_i is a linear function of G_i, the global conservation constraint can be transformed into a general quadratic form.

General Form of the Global Conservation of Energy Constraint

where α_ij is a matrix of constant coefficients, β_i a column vector of constant coefficients and δ is the total demand D_Tot. The first term is quadratic in G_j and is the losses term. The second term is linear in G_j and is the generation term.

The other constraints on the power flows are the line limits. Power flows through transmission lines are given as part of the DC power flow equations and these are constrained to lie within maximum and minimum flow limits.

Line Constraints

H is called the transfer admittance matrix for the Grid.

Total reserve constraints for each region are:

Reserves Constraints

: i ∈ I

where I is the set of reserve offers in an island and G_i is a significant generation dispatch in that region.

PLSR Reserve Constraints

R_i < rslope_iG_i , R_i < G_{max, i} - G_i

where PLSR reserve offer i is associated with generation offer j and rslope_j is the maximum reserve that can be provided as a proportion of generation G_j.

Apart from the general form of the global conservation of energy constraint (see section 3 2.1.3 on page 4), all constraints are linear. While the problem is in this form it can only be described as a general optimisation problem. It is unclear even whether the problem is convex or if it has a single local minimum. It can be demonstrated however that it has a single local minimum under certain conditions (An important condition is that negative offers are not on the margin. Negative offers are not currently permitted under the EGRs) and also that these conditions can be used to formulate our particular problem of optimising the dispatch of generation on the grid.

One difficulty to overcome is the presence of zero price offers: if a combination of these can satisfy all conditions then the solution is degenerate (Degenerate simply means that there is more than one dispatch that satisfies the conditions that the Cost objective function is minimised. In this case the minimum Cost is zero). It is necessary to make an arbitrary decision in this case and IQP actually assigns a price of $0.01/MWh to all zero offers. This is also consistent with the EGRs which prevent degeneracy through the Must-run Dispatch Auction in which generators bid for the right to offer zero, but only up to a major portion of total demand.

Now we can make some assumptions about the energy available to be dispatched.

Available Energy

EA is a variable representing available energy. This somewhat arbitrary value is the total demand that could be satisfied by a dispatch of plant if net injection at the swing bus was adjusted to make EA equal to total demand.

Since losses are only a small fraction of generation for all realistic situations we can make an assumption for the partial derivative of EA.

Partial Derivative of the Available Energy

Given that no negative price offers are included, we can replace the pricing objective function (on page 2) with an inequality constraint.

Available Energy Constraint

EA ≥ D_Tot

The partial derivative shows that EA can always be reduced by reducing generation, which in turn can be achieved at a lower or equal Cost. So for any dispatch satisfying the available energy constraint there can always be found one with less or equal cost that satisfies the global conservation of energy constraint (on page 3). This means a minimal solution to the problem with the global conservation of energy constraint, will also be a minimal solution to the problem with the available energy constraint.

If Γ(E) is defined as the set of dispatches where EA(d) ≥ E and which also satisfy all the linear constraints, then the optimisation process can now be stated as:

Problem 1: minimise (Cost(d)) where d ∈ Γ(D_Tot).

A set π(P) can be defined as the set of dispatches that satisfy all the linear constraints and Cost(d) ≤ P. Since Cost is a linear function of the dispatch variables, this set is entirely defined by linear constraints. The maximisation problem:

Problem 2: maximise (EA(d)) where d ∈ π(P)

Is a member of the quadratic programming family since it involves only linear constraints and has a quadratic objective function. Furthermore EA(d) is a convex function so any local maxima of this problem will also be a global maxima.

It can be shown that if d1_Opt is a solution to Problem 1 and d2_Opt is a solution to:

Problem 2a: maximise(EA(d)) where d ∈ π(P_s) and P_s = Cost(d1_Opt)

then d2_Opt is also a solution to Problem 1.

This can be shown:

Firstly, it is true that d1_Opt ∈ π(P_s) since Cost(d1_Opt) = P_s.

Therefore EA(d2_Opt) >= EA(d1_Opt) and this means that d2_Opt ∈ Γ(D_Tot)

Since Cost(d2_Opt) ≤ P_s, the cost of an optimal solution to Problem 1, d2_Opt must also be an optimal solution to Problem 1.

The Iteration Process

We found above we can solve Problem 1 by:

finding the cost of the optimal dispatch, P_s
then solving the quadratic programming problem, Problem 2.

This is done in the IQP algorithm by first making an educated guess for a value for P, then solving Problem 2 and using the EA_max of the solution to re-estimate P, repeating until a solution is found with an EA_max = D_Tot. It usually takes only two or three iterations in order to find a very good estimate of P_s.

Figure 1 - Available Energy as a Function of Cost

One by-product of the optimisation of Problem 2 is the Lagrange multiplier for the cost constraint which is the inverse of . This allows a new figure of P to be estimated by the Newton method which gives convergence on a solution for P_s very quickly.

Figure 1 illustrates a typical curve. In this case there is an infeasible region to the left caused by the northward constraint on the HVDC link (The DC equations balance the demand in the North Island with generation from the swing node, which in this example is in the South Island, causing a high northward flow when dispatched generation is low).

The point where EA_max meets the Total Demand line, D_Tot, is the solution point. Since the curve has a smoothly decreasing slope toward higher values, the Newton method works well from a first estimate of P that is not too much greater than P_s.

To deal with overshoots that are a feature of using the Newton method the algorithm keeps a track of a low P value and a high P value and uses the bisection method if these limits are exceeded by the Newton method. If a solution is first found infeasible then P is increased until either a solution is found or the whole problem is declared infeasible.

Convergence and Accuracy of IQP

The iteration process continues until convergence is obtained, i.e. until the available energy, equal to total generation less the total losses, is within 0.0001 MW of the required demand, which typically occurs within 3 to 5 iterations.

If the IQP algorithm fails to converge quickly to this value then it switches to a modified bisection method which generally converges. If a dispatch is ultimately infeasible then IQP returns an error code 1, indicating a no solution error.

IQP can correctly produce negative prices, but will not accept negative offers. Zero offers are set to $0.01/MWh so if there is enough low priced generation to satisfy demand plus losses, then the losses determine the dispatch.

Inclusion of Losses

Figure 2 - Modelling of Losses

A DC power flow (on page 10) does not include losses in the flows calculated in each line of the Grid, but the losses can easily be calculated separately. Losses are therefore included in dispatch for each line by calculating the line loss on the line and adding these to the demand at the node at the receiving end of the line, as occurs in SPD.

In the figure to the right the dark arrow at the receiving node shows the actual demand. The lighter arrow shows the demand including actual demand plus the line losses on this line.

After IQP is run once with nil distributed losses, the line losses are recalculated and the dispatch is recalculated.

Nodal Pricing

Nodal energy prices are the shadow prices of the nodal energy balance constraints applied in SPD and the island reserve prices are the shadow prices of the island reserve constraints applied in SPD.

The shadow price of a constraint is the change in the optimum Cost when the constraint is relaxed by one unit. To calculate the shadow price re-optimisation is required when the constraint is relaxed, so shadow prices must not be confused with the change in Cost from a simple change in generation by 1 MW at a node.

In IQP, the shadow prices are available as Lagrange multipliers on the appropriate constraints.

Note that the EGRs also specify that a DC power flow (on page 10) is to be used when dispatching and calculating prices. This is a power flow that makes simplifying assumptions about the relationship of the reactance and resistance of each line and also doesn't calculate the losses on each line.

SPD is a large scale LP and its DC power flow is incorporated directly into the LP formulation. The more traditional approach is to use a separate DC power flow and a non-linear program, as occurs within IQP.

A standard DC power flow calculates very quickly and implicitly achieves global energy balance at all nodes except perhaps at the "swing bus." IQP therefore does without the many nodal energy balance constraints within SPD and uses only the one global conservation of energy constraint:

where L is the total losses on the Grid.

With only one energy constraint however, IQP requires formulae that calculate the desired nodal energy prices from the Lagrange multiplier of the global conservation of energy constraint and other Lagrange multipliers that are available from IQP. The following derivation shows how this is achieved.

In addition to aggregating the nodal energy balance constraints into a global conservation of energy constraint, combine the reserve and risk constraints as follows:

The IRAFs are now almost always equal to 1. Form the Lagrangian function where

is the cost component due to generation and

is the cost component due to reserves.

Lagrangian for the Cost Objective Function

Now differentiate with respect to nodal generation G_j to get the following equations which must hold at minimum points of the Cost function.

Nodal Pricing Equations

where the P_n are the nodal prices and the P_G,n can be thought of as the energy component of the nodal price at node n. The ρ_x are the various Lagrange multipliers or shadow prices from IQP.

Note the shadow price terms for reserve apply only if generation at the node is setting the reserve risk either directly or via the HVDC link.

The above equation has been verified in extensive testing during the investigation leading to the first prototypes of EMarket in 1997.

The location factor is defined as , where P_r is the appropriate reference nodal price, which then gives the full equation for location factors.

Marginal Location Factors

Note that in practice the location factors are typically calculated directly from the appropriate nodal prices. The equation shows that in the absence of active line constraints the major determinant of location factors is the relative effect of marginal losses through the terms and . In this case the location factors can be calculated directly from the power flow input data noting that

where B is the susceptance matrix described in section 4 and losses are calculated according to the result shown in section 4.2. The B_nj are calculated directly from the transfer admittance matrix described in the next section.

DC Power Flow

IQP incorporates a more or less standard DC power flow based on three assumptions:

per unit voltages are all equal to 1;
phase angle differences between nodes are all small;
line reactance exceeds line resistance.

SPD works with resistance, R_i, reactance, X_i, and with susceptance of line i, B_i expressed as (In some texts the sign of B is positive):

for AC lines and optimises the HVDC link separately based on its resistive losses. IQP is entirely based on a DC power flow and treats the HVDC link as an AC line in the power flow. Since the HVDC link is not within a loop in the grid, this approach is very accurate for total HVDC transfer.

In a grid with n nodes and m lines we will refer to the n-1 × n-1 diagonal matrix with the B_i as elements as the Ω matrix.

The DC power flow is formulated as a matrix equation - the ‘A’ matrix has a 1 in a column if the line in that row starts at the node in that column, -1 if the node is the ending node (Choice of starting and ending node is arbitrary and only affects the sign of the flows in the lines). Everywhere else it has zeros. In other words the A matrix specifies how lines are connected between nodes.

The m×m R matrix contains the line per unit resistances on the diagonal, zeros elsewhere. The susceptance matrix (For a full derivation of the DC power flow equations in the nodal pricing context see Appendix D of Spot pricing of electricity, Schweppe, Caramanis, Tabors and Bohn, Kluwer Academic Publishers, 1997), B, is given by:

The transfer admittance matrix, H, is an m×n-1 matrix given by:

which also leads to

Transfer Admittance Matrix

Each element of the H (transfer admittance) matrix H_ij is the rate that flow on line i varies with net injection at node j, i.e.

where Z_i is the flow in line i and y_j the net injection at node j. The net injection is just generation less demand y_i = G_i - D_i. Where demand is not able to be dispatched on price, as occurs in the pricing runs of SPD, it is taken as fixed so that for the purpose of pricing

where G_j is the total injection at node j.

The DC power flow matrix equation for the line flows is then

Z = Hy

where Z is the vector of all line flows and y the vector of net injections at all but the swing node. The swing node may be chosen arbitrarily but in EMarket and EMO it is chosen outside of IQP as a well connected node to maximise the speed at which the full power flow with losses converges.

Losses

The assumption of X_i > R_i leads to losses in line i given by

L_i = R_i Z_i² and total losses are given by

L = y^T By.

The per unit system used in electrical engineering requires a base power to which all parameters are referenced, including the line resistance, reactance and susceptance. In the case of IQP we use the same as in SPD, ie. 100 MW. Note that the equations above require scaling in some instances to get sensible units. For example, Z = 100Hy and give units of MW. The net injection y_i also needs to be scaled down by 100.

A DC power flow does not include losses in the line flows, even though the losses can be calculated. The typical approach to including losses is either to ignore losses (which is surprisingly accurate for pricing purposes) but creates the problem of a non-physical dispatch for the System Operator to deal with, or to add losses to the demand at the node at the downstream end of each line (this is also what happens in SPD). The latter requires iterations to get convergence but convergence is fast and stable in IQP.

Note that any imbalance between losses, generation and demand will always appear at the swing node. Since IQP performs a full market dispatch which balances losses, demand and generation, the swing node will always be in balance unless an infeasibility occurs.

Differences between EMarket and EMO

Both EMarket and EMO include a wide range of inputs for the grid so both explicitly take account of line limits, equation constraints and losses.

EMarket and EMO are designed for different purposes so there are some differences in the way the power flow works in the two models, as shown below:

	EMO	EMarket
Grid data and detail	Grid data obtained from em6 or similar data source, hence every node, line and related constraint modelled, including equation constraints.	Capable of solving the full grid but not generally done due to the reduction in performance – typically uses an aggregated grid with approx 180 nodes and 200 lines. Transformer losses between modelled nodes and lower voltage nodes are modelled at each node by entering transformer resistance values. Default line limits can be entered for summer and winter in line with Transpower time zones, or re-rated on the fly during simulations runs from the EMarket Schedule.
Line outages	The power flow is re-solved at every period modelled so full or partial outages are catered for. This slows the power flow solve down considerably but techniques for solving large, sparse matrices are employed without loss of accuracy.	The power flow equations are not re-solved each modelled period in EMarket so full outages during a simulation run can not be modelled. Partial outages (re-ratings) are modelled simply by scheduling the re-rating. Re-rating to zero is not allowed so 1 MW is the minimum line limit.

Aggregated Lines in EMarket

SPD models hundreds of nodes and all individual lines and transformers. EMarket can model the full grid but to do so would slow EMarket down sufficiently to be counterproductive in most applications, with little or no gain in power flow accuracy.

Most users model the grid with about 180 nodes and 200 aggregated lines which gives high accuracy, excellent performance and has enough nodes that aggregating lines is straightforward – in most cases lines are either entered directly or aggregated with one or more lines in parallel. Using 180 nodes more or less eliminates the need to aggregate loops within the grid.

For n lines in series the per unit impedances (The use of Z to denote impedance should not be confused with its earlier use to denote power flow) of the lines are simply summed:

Z = Z₁ + ... + Z_n

where Z = R + jX using complex notation.

For the more common case of lines in parallel the impedances of the aggregated line is the reciprocal of the sum of the reciprocals of the individual line impedances. For two lines parallel:

The line limits also need to be adjusted. Lines in parallel typically have similar impedances and similar limits applied to them in SPD, in which case the limit of the aggregated line is the sum of the individual limits.

Where parallel lines have significantly different impedance then the lines with lower resistance will carry more power than the lines with higher per unit resistance. Most of the time in EMarket simulations are undertaken with a limited number of lines with limits enabled (The user would typically enable constraints on key lines and monitor all other lines to check if any are exceeding their limits to an extent that would warrant re-running with their constraints enabled), so the effort required to set limits on aggregated lines is typically minimal.

As the New Zealand grid is relatively static in its configuration, at least for the purposes for which EMarket is used, once a grid is aggregated it can be maintained with ease. Users will often liaise with Energy Link to obtain a grid suitable for immediate use.

Transformer and Fixed Loss Modelling in EMarket

As EMarket is typically run with an aggregated grid with 3-letter nodes representing all nodes at a location, transformer losses would go unaccounted for without explicit provision being made. EMarket’s grid configuration includes, for each node, an island for reserves purposes, and an associated resistance value. The resistance value is chosen by the user to represent the combined resistance of the transformers at each modelled node. EMarket models demand at the nodes as passing through a line with this resistance, and thus accounts for transformer losses on the grid with the associated direct or indirect impact on total losses, power flows, dispatch and nodal prices.

Fixed losses on the grid arising from transformers and the HVDC link amount to around 30 MW and are typically added to demand at strategically selected nodes, using functions in the EMarket Schedule.

DC Power Flow Accuracy

In a theoretical sense, the DC power flows are constructed in the same way as the DC power flow in SPD and thus by construction the power flow modelling and agreement with SPD is excellent. SPD is, of course, a linear model, but the impact on the accuracy of its power flows, as compared to the full quadratic loss power flows implemented in IQP, is small in most lines on the grid.

Back to Technical Bulletins