FERC Technical Report on Loss Estimation

Marginal Loss Calculations for the DCOPF

Brent Eldridge∗1,2, Richard P. O’Neill†1, and Anya Castillo‡3

1Federal Energy Regulatory Commission§

, Washington, DC, USA

2Johns Hopkins University, Baltimore, MD, USA

3Sandia National Laboratories¶, Albuquerque, NM, USA

January 24, 2017

Abstract

The purpose of this paper is to explain some aspects of including a marginal line loss approximation in the DC optimal power flow (DCOPF). The DCOPF optimizes electric generator

dispatch using simplified power flow physics. Since the standard assumptions in the DCOPF

include a lossless network, a number of modifications have to be added to the model. Calculating

marginal losses allows the DCOPF to optimize power generation, so that generators that are

closer to demand centers are relatively cheaper than generators that are far away. The problem

formulations discussed in this paper will simplify many aspects of practical electric dispatch

implementations in use today, but will include sufficient detail to demonstrate a few points with

regard to the handling of losses.

First, we examine marginal line loss approximations in the DCOPF and how different methods

effect LMP pricing. The methodology explained in this paper begins with a feasible AC power

flow solution, called the base point or operating point. This base point includes information

about network power flows and bus voltages that affect the calculation for marginal line losses.

We show that when these aspects are ignored, prices no longer reflect the network’s physics.

Various DCOPF model formulations also affect the accuracy of the optimal solution’s physics.

Selecting a reference bus simplifies calculations in the DCOPF. We compare a few common

formulations of the DCOPF and show that one of these formulations can distort flows and

results in a Kirchoff’s Current Law violation at the reference bus. Correcting this formulation

results in a model with optimal solutions that are independent of the reference bus.

Additionally, we propose a novel method for updating the loss approximation without solving

for a new base point. If the update procedure converges, then it gives a solution to a nonlinear

problem. Results show rapid convergence properties on all networks tested.

∗B. Eldridge can be contacted at b.eldridge@jhu.edu

†R. P. O’Neill can be contacted at richard.oneill@ferc.gov

‡A. Castillo can be contacted at arcasti@sandia.gov

§The views presented are the personal views of the authors and not the Federal Energy Regulatory Commission or

any of its Commissioners.

¶Sandia National Laboratories is a multi-program laboratory managed and operated by Sandia Corporation, a

wholly owned subsidiary of Lockheed Martin Corporation, for the U.S. Department of Energys National Nuclear

Security Administration under Contract DE-AC04-94AL85000.

1
Contents

1 Introduction 3

1.1 Current Practices . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4

1.2 Literature Review . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5

1.3 Notational Conventions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7

2 Power Flow Derivations 9

2.1 DC Power Flow . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10

2.2 Marginal Line Losses . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11

2.3 Alternative Line Loss Derivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14

3 Model Formulation 14

3.1 Loss Distribution Factors and Kirchoff’s Current Law . . . . . . . . . . . . . . . . . 15

4 Model Comparison 18

4.1 LMP Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19

5 Quadratic Update Procedure 20

5.1 Motivating Example . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20

5.2 Algorithm Description . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21

5.3 Algorithm Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24

6 Conclusion 25

List of Figures

4.1 LMP Comparison of Linear Models . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19

5.1 Convergence Results for Algorithm 2. . . . . . . . . . . . . . . . . . . . . . . . . . . 25

List of Tables

1.1 ISO Loss Factor Methodologies. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6

3.1 Distorted Power Flow without Loss Distribution Factor . . . . . . . . . . . . . . . . 17

3.2 Consistent Power Flow with Loss Distribution Factor . . . . . . . . . . . . . . . . . . 17

3.3 KCL Violations without Loss Distribution Factor. . . . . . . . . . . . . . . . . . . . 18

4.1 IEEE 300-bus Test Case Solution Statistics . . . . . . . . . . . . . . . . . . . . . . . 20

5.1 Two Node Example . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21

5.2 Solutions for Initial and Final Bids . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21

List of Algorithms

1 Zero-Centered Quadratic Update . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22

2 Generic Quadratic Update . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24

2
1 Introduction

All independent system operators (ISOs) in the US implement locational marginal cost pricing

[1, 2, 3, 4, 5, 6, 7] in which market participants pay or receive the cost of delivering the next

unit of power at their node in the network. The marginal cost pricing approach is economically

efficient in a competitive market because the price signal to each node reflects the increase in

system cost required to serve the next unit of demand. Marginal loss prices are a component of

marginal pricing and reflect the portion of the change in cost that is due to a change in system

line losses. The locational marginal price (LMP) is the primary economic signal in ISO markets

and decomposes into the marginal loss component, marginal congestion component, and marginal

energy component.

However, approximations within the market dispatch model differ from the network physics. The

approximations result in prices that do not reflect physical measurements, and this causes problems

in the market since prices do not reflect actual marginal costs [8, 9, 10, 11]. This paper focuses

on making the modeling approximation as close as possible to the actual physics because this will

ensure that prices accurately reflect the marginal cost of electricity.

The magnitude of loss payments also justifies a closer look at current practices. In PJM in 2014,

total marginal loss costs were $1.5 billion, compared to $1.9 billion in total congestion costs [12].

It is important that this money is charged accurately since prices that accurately reflect locational

prices are a cornerstone of ISO market design.

Since losses are approximately quadratic, marginal losses are about twice the average losses. As a

result, ISOs will collect more revenue for line losses than what is paid to generators, and this money

is returned to demand based on load ratios in both the DAM and RTM. The over-collection must be

repaid with a rebate to market participants [8, 9]. The methodolgy to determine the rebate can have

significant effects. For example, a 2010 CAISO study showed that two alternative loss allocation

methodologies would change regional allocations by $18.8 million and $13.8 million compared to

the filed methodology [13]. Similarly, marginal loss rebate policies can also be exploited by market

participants [10, 14]. In addition to accurate prices, allocation methodologies are also potentially

important but will not be discussed further in this paper.

It is important to precisely study the loss approximation so that the market optimization will

model the actual network physics as closely as possible. For example, poor loss modeling can be

exploited by financial market participants who will place bids to correct for a poor loss approximation. Poor loss estimation can be caused from load forecast bias or can be inherent to the power

flow and loss estimate methodology, and a consistent over or under-estimation between day ahead

and real time marginal loss components is all that is needed for financial market participants to

place bids based on the mis-estimation. These bids can correct the mis-estimation by aiding price

convergence, but this would be unnecessary if the market cleared using a better loss approximation. A better solution may be for the market software to have a good loss approximation from the

outset. MISO changed its loss modeling to limit such behavior, as prompted by its market monitor

[11].

The general name for the dispatch problem that ISOs solve is the optimal power flow (OPF). The

AC Optimal Power Flow (ACOPF) accounts for alternating current’s mathematical complexities.

However, the ACOPF is a large scale, nonlinear, non-convex optimization problem and requires

more time to solve using existing methods than current practice allows [15]. Dispatch models must

solve quickly in order to be practical in day-ahead and real-time markets (DAM and RTM), which

is why today’s ISOs solve linear programming models. The DCOPF is named a bit awkwardly

because it is not modeling “direct current” power, but is really a linearization of the ACOPF [16].

3
1.1 Current Practices

ISOs typically implement the DCOPF with linear power transfer or generation shift factors which

relate power generation and demand to power transfers across transmission lines in the network.

We call this implementation the distribution factor model. The sensitivities in the model can be

linearized inputs from a feasible AC power flow solution. A loss approximation is also incorporated

into the market software’s economic dispatch. Typically, the approximation is based on historical

ratios or an AC power flow solution that predict load flow in the time period being dispatched.

The distribution factor model requires the selection of a reference bus which is assumed to be

the marginal source (or sink) of any changes in power consumed (or produced). Power flows to and

from the reference bus are “summed” using the superposition principle, and therefore the effect of

the reference bus gets canceled out in a lossless model. Although the reference bus simplifies the

mathematics, its inclusion in the model can distort power flows when line losses are considered.

A common alternative to the distribution model approach is called the “BΘ” model and also

results in a linear model. However, the BΘ model takes a few orders of magnitude longer to solve

and therefore is not used to clear markets. Therefore, this paper will focus on the distribution

factor model implementation of the DCOPF.

Marginal loss calculations in the DCOPF are sensitive to many things, including the input data,

the approximation approach, and the selection of a reference bus or slack bus. We use the terms

reference bus and slack bus interchangeably. Input data may include physical properties of the

transmission network and perhaps voltage angles and magnitudes. For example, a feasible AC

power flow can be used to supply the input data that defines the base point from which the

dispatch model optimizes.

Loss factors define the sensitivity of system losses to power injections or withdrawals at a specified

bus on the network. They can be positive or negative. When the loss factor at a bus is positive,

a small injection at that node will result in a small increase in system losses. Loads pay a lower

price because system losses are decreased by a small increase in demand. Generators receive a

lower price because system losses are increased by a small increase in production. On the other

hand, when the loss factor at a bus is negative, a small injection will result in a small decrease in

system losses. Loads pay a higher price because system losses are increased by a small increase in

demand. Generators receive a higher price because system losses are decreased by a small increase

in production. The change in system losses is relative to delivering power to the reference bus. By

definition, the loss factor at the reference bus is zero.

Typical dispatch models assume that losses are linearized around a base point solution. This base

point may be taken from the state estimator, AC power flow analysis, or the results of a dispatch

optimization. Each ISO’s processes for estimating losses are described below, and a summary is

located in Table 1.1.

California ISO (CAISO) determines marginal loss factors by linearizing around an AC power

flow solution base point [17]. The AC power flow is calculated throughout iterations of the Security

Constrained Unit Commitment (SCUC) process. System losses are calculated at after each AC

power flow. In addition, loss sensitivities and shift factors are calculated from linearizing around

the AC power flow solution then fed into the SCUC and SCED optimization models. In the

Integrated Forward Market (IFM) piece of CAISO’s DAM, SCUC uses generation and demand

bids to determine power flows. The IFM is followed by Reliability Unit Commitment (RUC) which

uses generation bids and a load forecast.

Although Electric Reliability Council of Texas (ERCOT) uses LMPs, the LMPs only include

components for energy and congestion [18]. In the absence of transmission congestion, LMPs are

uniform. Losses are added during the settlement process and are based on linear interpolation or

4
extrapolation of forecasted on-peak and off-peak transmission loss factors [19].

ISO-New England (ISO-NE) uses loss factors to calculate LMPs every hour in the DAM and every

five minutes in the RTM. The state estimator provides information regarding transmission losses

for both Day-Ahead and Real-Time LMPs. ISO-NE bases its DAM on the expected transmission

configuration and the bids and offers from market participants. The RTM clears off information

from the state estimator [20].

In Midcontinent ISO (MISO), the Energy Management System (EMS) state estimator calculates

the total system losses using a combination of an AC power flow and a statistical model based

on system measurements [21]. In real time, loss factors can be calculated directly from the state

estimator, and MISO monitors the calculated real-time loss factors to make sure these are adequate

for settlement purposes. For day ahead, MISO uses recent solutions from the state estimator with

similar load and wind characteristics as day ahead interval.

The New York ISO (NYISO) uses marginal loss factors that reflect expected scheduled and

unscheduled power flows on the network [22]. In the DAM, expected unscheduled power flows are

generally determined using a 30-day moving average of on and off-peak flows. Unscheduled power

flows in the RTM are based on current power flows. NYISO calculates LMPs multiple times for

each time period. The first set of LMPs settle the DAM and are taken from unit commitment and

dispatch optimization. The final LMPs settle the RTM and are taken from a dispatch optimization.

Prior to implementing marginal loss pricing, PJM used generic on-peak and off-peak loss factors,

adding 3% and 2.5% to on and off peak demands, respectively. Because loss factors were not part

of the economic dispatch, the result was less than optimal [23]. PJM now calculates loss factors

in the DAM and RTM based on transmission characteristics, generation levels and load levels, and

state estimator data [6].

The Southwest Power Pool (SPP) operates its DAM and RTM with marginal loss pricing. In

the RTM, losses are estimated using the current state estimator solution. SPP estimates future

operating conditions and performs a power flow study to calculate marginal losses for the DAM [7].

Table 1.1 provides a brief summary of the processes used by each ISO.

1.2 Literature Review

The DCOPF has long been of interest to academic research, and the following section will review a

small sample of academic work on the subject. DCOPF is a subset of the more generic optimal power

flow (OPF) problem, which is a large-scale, nonlinear, non-convex problem that is exceptionally

difficult to solve. This problem was first formalized as an optimization problem by Carpentier in

1962 [24]. This sparked interest in formulations for electricity markets [25].

Many surveys give a more comprehensive review of the various methods for solving OPF problems

than we provide here [26, 27, 28, 29, 30, 31]. The surveys serve a dual purpose to anyone interested in mathematical optimization because of their close tracking of mathematical programming

developments. In particular, linear programming, quadratic programming, generalized reduced gradient, Newton’s method, and conjugate gradient were all common approaches through the 1980s

and 1990s. More recently, semidefinite programming [32, 33, 34, 35] and second order conic programming [36] have shown promising results. Many of the recent advances focus on solving the

ACOPF problem [15, 37], but the linear DCOPF problem remains the standard problem for electric

dispatch applications [38, 39].

Computational performance has always been the main advantage of using linear OPF models, as

well as its easy integration with standard economic theory [25]. It was first formulated and solved

by Wells in 1968 [40], which reports solution times of a few minutes on power networks consisting

of 100 nodes. This work led to interest in more efficient DCOPF formulations and loss sensitivity

5
Table 1.1: ISO Loss Factor Methodologies.

ISO Used in

Dispatch

Base Point (DAM) Base Point (RTM) Update Frequency

CAISO Yes SCUC AC power flow with

generation and demand

bids or load forecast

Same as DAM Every hour in DAM

and every fifteen minutes in RTM

ERCOT No Interpolation during settlement processSame as DAM N/A

ISO-NE Yes State estimator solution

with estimated future operating conditions

Current state estimator solutionEvery hour in DAM

and every five minutes

in RTM

MISO Yes Recent state estimator solutions with similar demand and wind characteristicsCurrent state estimator solutionMonitored in real time,

with updates possible

up to every minute

NYISO Yes SCUC AC power flow using

a 30-day moving average of

on and off peak power flows

SCUC AC power

flow

During intermediate

dispatch runs between

DAM clearing and

RTM clearing

PJM Yes State estimator solution

with estimated future operating conditions

Current state estimator solutionEvery hour in DAM

and every five minutes

in RTM

SPP Yes AC power flow after RUC

with estimated future operating conditions

Current state estimator solutionEvery hour in DAM

and every five minutes

in RTM

calculations [16, 41, 42, 43, 44]. In addition to computational advantages, the current formulations

of the DCOPF show impressive–but not perfect–accuracy relative to AC power flows [38].

Various new approaches to the DCOPF remain an active area of research. Iterative approaches

to the DCOPF in [39, 45, 46, 47] have shown some success at improving the physical accuracy of

the model. LMP pricing is one of the most important DCOPF applications [47, 48, 49, 50, 51, 52],

but the DCOPF is certainly not restricted to real time dispatch. It is also an important aspect in

transmission expansion planning [53, 54], renewable energy and storage integration [55], and other

applications that are not enumerated here.

However, many previously proposed methodologies either require nonlinear or piecewise linear

constraints [45, 54, 56] or additional solutions to an AC power flow [47]. Some use dispatch formulations that are not scalable for large networks, such as the BΘ power flow approximation [49]

or current-based models [50]. Some use a modeling construct called fictitious nodal demand which

embeds line losses in demand [39, 45, 46, 51] or exclude losses all together [52].

The key advantage of the proposed methodology is that it is based on the same DCOPF formulation used in current market software and can be easily integrated into ISO operations.

The rest of the paper is organized as follows. Section 2 derives power flow and line loss sensitivities. Section 3 formulates the distribution factor model. Section 4 compares LMPs from

three different loss approximations. Section 5 describes an iterative method to improve the loss

6
approximation. Section 6 concludes the paper and is followed by references.

1.3 Notational Conventions

The following sections will present various equations and optimization problems concerning optimal

power flow. To make this as readable as possible, we will make an honest effort to adhere to the

following conventions. Unfortunately there are some places where historical precedent supersedes

our best efforts.

Scalars, vectors, and matrices will all follow a general convention. Scalars will be lowercase (x),

vectors will be uppercase (X), and matrices will be bolded uppercase (X). Elements of vectors

will be noted with a subscript (xn), and we may also superscript some scalars to denote separate

values (x

p

n, x

g

n). Vectors are generally column vectors, and may be differentiated using subscripts

in some cases (Xp, Xg), as will matrices (Xd). [X]ij denotes the element of X in the i

th row and

j

th column. Matrix transpose is denoted by the symbol >.

Parameters and variables are more difficult to separate by a convention. Instead, this difference

will be noted immediately after the parameter or variable is introduced. They are also listed as

parameters or variables in the nomenclature list at the end of this section for easy reference. Fixed

variables are denoted by an overline (x) and optimal solutions are denoted by an asterisk (x

∗

).

Dual variables will use the Greek alphabet, although one exception to this rule is θ and φ, which

are used for voltage angles in common practice.

Sets names are in calligraphic font (N ). We have simplified the notation a bit by making no

distinction between generator and node indices, i.e., the power injection at each node is equal to

the sum of injections from all generators at that node. Therefore, the formulations in this paper

may require some additional bookkeeping to be implemented in an actual optimization program.

Nomenclature

Dual Variables

αmin, αmax Dual variables to the minimum and maximum generation output constraints.

λ Dual variable to the power balance constraint.

μ Dual variable to the transmission constraint.

σ Dual variable to the loss function constraint.

Λ Vector of locational marginal prices (LMPs).

Functions

◦ Hadamard product.

> Matrix or vector transpose.

c(·) Linear or convex cost function.

Indices

h Iteration index.

i, j, n Nodes or bus indices, i, j, n ∈ N .

7
k Transmission line index, k ∈ K.

Parameters

aijk Transformer tap turns ratio from node i to j on branch k.

bijk Branch susceptance from node i to j on branch k.

gijk Branch conductance from node i to j on branch k.

`

0 Line loss function constant.

rk Resistance on branch k.

ui(n) Unit injection equal to 1 if i = n or zero otherwise.

xk Reactance on branch k.

γk Loss approximation second-order coefficient.

ηk Loss approximation range translation.

ξk Loss approximation domain translation.

ω Damping parameter.

1 Vector of ones.

LF Loss factor vector with elements `fi.

Pd Power demand vector with elements p

d

i

.

Pmin, Pmax Minimum and maximum generator output vectors.

Tmax Transmission branch power flow limit vector.

W Weighting vector with elements wi.

Φ Transformer phase angle change vector with elements φijk for branch k from node i to

j.

A Network incidence matrix.

Bd Diagonal branch susceptance matrix with elements bijk for branch k from node i to j.

B Nodal susceptance matrix.

Gs Diagonal shunt conductance matrix with elements g

sh

i

for bus i.

I Identity matrix.

L Marginal loss matrix.

S Marginal branch flow matrix.

T Transmission sensitivity matrix with elements t(ijk,n)(PTDFs).

Sets

8
K Set of transmission lines, {1, . . . , K}.

N Set of nodes or buses, {1, . . . , N}.

Variables

` Total network real power losses.

`k Line losses on branch k.

y Slack variable for a reference bus withdrawal.

P Net real power injection vector with elements pi for node i.

Pg Power generation injection vector with elements p

g

i

for generators for node i.

Pt Power flow vector with elements pijk for branch k from node i to j.

V Voltage magnitude vector with elements vi for node i.

Θ Voltage angle vector with elements θi for node i.

∆P Marginal real power injection vector.

∆V Vector of marginal change in voltage magnitude.

∆L Marginal line loss vector.

∆Θ Marginal voltage angle vector.

∆Θ Voltage angle sensitivity matrix.

2 Power Flow Derivations

The purpose of this section is to derive the linear DC power flow equations that are used by ISOs.

Three standard assumptions are made to derive the DC power flow approximation from nonlinear

AC power flow equations. These assumptions result in no losses, so line loss approximations are

calculated based on an AC-feasible operating point or base point. Two approximations are presented

for line losses, one which uses voltage information from the base point solution and one which uses

some assumptions to make losses a quadratic function of real power flows.

We begin with the AC real power flow equations for a single branch and network balancing

constraints for each node. Notationally, a branch k ∈ K connects nodes i and j ∈ N , and power

flows pijk from i to j and pjik from j to i are given by,

pijk = gijk

v

2

i

a

2

ijk

−

vivj

aijk

(gijk cos(θi − θj − φijk) + bijk sin(θi − θj − φijk)), (1)

pjik = gijkv

2

j −

vivj

aijk

(gijk cos(θj − θi + φijk) + bijk sin(θj − θi + φijk)), (2)

where the parameters are the branch conductance gijk, branch susceptance bijk, tap transformer

turns ratio aijk, transformer phase angle change φijk, and the variables are the voltage magnitude

vi and voltage angle θi.

9
According to Kirchoff’s Current Law (KCL), power flows balance at each node or bus. The

amount of power generated minus the amount consumed at a node must be equal to the amount

flowing in or out. Losses in the shunt conductance, g

sh

i

, are also accounted for in the nodal balance.

For now, we simplify power generation (an injection) and consumption (a withdrawal) for now

using the net injection pi at node i ∈ N , which by convention is positive for a net injection and

negative for a net withdrawal. For real power, the network balance equations are:

pi =

X

k

pijk + v

2

i

g

sh

i

, ∀i ∈ N . (3)

Or, in matrix form:

P = APt + Gs(V ◦ V ), (4)

where A is an (N × K) network incidence equal to 1 for branch a k assumed to flow into node i,

−1 if the branch is assumed to flow out of node i, and 0 if branch k is not connected to node i, Pt

is a vector of transmission flows, Gs is an N × N diagonal matrix of shunt conductances, V is a

vector of nodal voltage magnitudes, and ◦ is the element-by-element (or Hadamard) product.

2.1 DC Power Flow

Many industry applications rely on DC power flow approximations. DC power flow equations are

preferable in many instances because they are linear and can be solved quickly. Conversely, AC

power flow equations model the system more accurately, but are non-convex. It can even be difficult

to find a feasible solution to AC power flow equations in a large scale system such as one of the

main US power grids. Therefore, the common DC power flow approximation makes three main

assumptions:

• Voltage is close to one per unit (p.u.) at all buses,

• Voltage angle differences are small, i.e., sin θ ≈ 0 and cos θ ≈ 1,

• Line resistance is negligible compared to reactance, i.e., rk  xk and therefore gijk  bijk.

Line conductance and susceptance are physical properties of the lines, and are respectively gijk =

rk

r

2

k+x

2

k

and bijk =

−xk

r

2

k+x

2

k

. Under the assumptions above (and also assuming that the tap turns ratio

aijk = 1), (1) and (2) reduce to

pijk = gijk

v

2

i

a

2

ijk

−

vivj

aijk

(gijk cos(θi − θj − φijk) + bijk sin(θi − θj − φijk))

≈ gijk − (gijk cos(θi − θj − φijk) + bijk sin(θi − θj − φijk))

≈ gijk − (gijk + bijk(θi − θj − φijk))

≈

−1

xijk

(θi − θj − φijk).

Or, in matrix form:

Pt = Bd



A>Θ + Φ, (5)

where Bd is a (K × K) diagonal matrix with values −1

xijk

, A is the network incidence matrix, Θ is

an (N × 1) vector of nodal voltage angles, and Φ is a K × 1 of transformer phase angle changes.

Similar derivations can be followed in reference [57].

10
To reduce solution time in practice, the transmission constraint (5) can be simplified using power

transfer distribution factors (PTDFs), also called shift factors. The PTDF is the sensitivity of real

power flow to the change in power injected at a particular bus. The injection (or withdrawal) is

assumed to be withdrawn (or injected) at the reference bus. The PTDFs form a K × N sensitivity

matrix T with elements t(ijk,n).

t(ijk,n) =

dpijk

dpi

≈

d

dpi

(

−1

xijk

(θi − θj − φijk))

=

−1

xijk 

dθi

dpi

−

dθj

dpi



.

It is common to define an N × N nodal susceptance matrix, B.

B = ABdA>

Now we introduce an N × N matrix [I − W1

>], whose nth column pairs a unit injection at node

n with a reference bus withdrawal. I is the identity matrix. Reference bus withdrawals are defined

by a weighting vector W that sums to one. The n

th column of matrix ∆Θ represents the marginal

change in bus voltage angles corresponding to the injections and withdrawals in the n

th column of

[I − W1

>].

[I − W1

>] = B∆Θ.

Then we can construct the PTDF matrix:

T = −BdA>∆Θ = −BdA>B

−1

[I − W1

>]

The matrix B is singular, so the row and column corresponding to the reference bus are removed

to make the matrix invertible. The PTDF at the reference bus is zero by definition.

This is equivalent to solving the following set of N equations N times:

ui(n) − wi =

X

k

−1

xijk

(θi − θj ), ∀i ∈ N , (6)

with n = {1, . . . , N}, where ui(n) is equal to 1 if i = n or zero otherwise, and wiis the i

th element

of the weight vector W.

2.2 Marginal Line Losses

Line losses `k across branch k are equal to the difference between the amount leaving node i and

flowing into node j. Since real power flow pijk is assumed to be the power flowing out of node i

from line k, and similarly −pjik is the amount flowing into node j from line k, losses are simply the

summation of (1) and (2):

`k = pijk + pjik = gijk 

v

2

i

a

2

ijk

+ v

2

j − 2

vivj

aijk

cos(θi − θj − φijk)

!

(7)

11
To linearize (7), let `k = `k(x), where x = [θi

, θj , vi, vj ]. Assuming a small distance ∆x from some

base point x, the first order Taylor series of `(x) at x = x+∆x = [θi

, θj , vi, vj ]+[∆θi, ∆θj , ∆vi, ∆vj ]

is:

`k(x + ∆x) = `k(x) + ∇`k(x)∆x

= `k(x) + ∂`k

∂θi

∆θi +

∂`k

∂θj

∆θj +

∂`k

∂vi

∆vi +

∂`k

∂vj

∆vj . (8)

Instead of vi and θi, it would be helpful to have `k as a function of the control variables pn. A

physical intuition of AC power systems is that voltage angles are tightly coupled with real power

and voltage magnitudes with reactive power, and therefore we will assume that voltages do not

change as a function of pn. Accordingly, we have ∆θi =

dθi

dpn

∆pn and ∆vi =

dvi

dpn

∆pn = 0. Then,

`k(P + ∆P) = `k(P) +X

n

∂`k

∂θi

dθi

dpn

+

∂`k

∂θj

dθj

dpn



∆pn +



∂`k

∂vi

dvi

dpn

+

∂`k

∂vj

dvj

dpn



∆pn



= `k(P) +X

n

∂`k

∂pn

∆pn. (9)

The partial derivatives introduced in (8) can be calculated directly from (7):

∂`k

∂θi

= 2gijk

vivj

aijk

sin(θi − θj − φijk) (10a)

∂`k

∂θj

= −2gijk

vivj

aijk

sin(θi − θj − φijk) (10b)

This creates a K × N matrix L of the sensitivity of line losses with respect to voltage angles:

[L]kn =







2gijk

vivj

aijk

sin(θi − θj − φijk), if n = i,

−2gijk

vivj

aijk

sin(θi − θj − φijk), if n = j,

0, otherwise.

(11)

Calculating the total derivatives dθi

dpn

and dθj

dpn

will require solutions to network equations. By

linearizing (1) and (2), we can use the network equations (3) with a real power injection ∆pn at

some node n ∈ N and solve for marginal changes in θi, and θj at each branch. The relevant partial

derivatives for changes in θ are:

∂pijk

∂θi

=

−vivj

aijk



bijk cos(θi − θj − φijk) − gijk sin(θi − θj − φijk)



, (12a)

∂pijk

∂θj

=

−vivj

aijk



−bijk cos(θi − θj − φijk) + gijk sin(θi − θj − φijk)



. (12b)

We construct a K ×N matrix S of the sensitivity of the branch power flows to changes in voltage

angles:

[S]kn =







−vivj

aijk



bijk cos(θi − θj − φijk) − gijk sin(θi − θj − φijk)



, if n = i,

−vivj

aijk



− bijk cos(θi − θj − φijk) + gijk sin(θi − θj − φijk)



, if n = j,

0, otherwise.

(13)

12
Since the network equations are already satisfied at the base point, only the marginal changes

need to be balanced. The injection ∆pn is paired with a withdrawal w = (w1, ..., wi, ..., wN ) at

an arbitrarily chosen reference bus to keep the equations feasible. The reference bus is defined by

assigning weights wi to each bus i ∈ N such that P

i wi = 1. The weights are often assigned by

the proportion of total load at each bus. The network equations (3) can be used to simultaneously

solve for the marginal changes to voltage angle at each node given a marginal injection at node n

and withdrawal at the reference bus. The marginal network equations are:

∆pi − wiy =

X

k

−vivj

aijk



bijk cos(θi − θj − φijk) − gijk sin(θi − θj − φijk)



(∆θi − ∆θj ), ∀i ∈ N .

(14)

The marginal injection ∆pi

is fixed to 1 at bus n and 0 otherwise. The same amount as the

injection may not be feasible to withdraw from the reference bus due to network losses, so we add

a slack variable y for the amount withdrawn. Assuming it exists, the solution to this system of

equations is {∆θ

∗

i

: i ∈ N }, and the derivative relating θi and pn is dθi

dpn

=

∆θ

∗

i

∆pn

.

Then we define the loss factor `fn as follows:

`fn :=X

k



∂`k

∂θi

dθi

dpn

+

∂`k

∂θj

dθj

dpn



=

X

k

2gijk

vivj

aijk



sin(θi − θj − φijk)





∆θ

∗

i

∆pn

−

∆θ

∗

j

∆pn



(15)

The result is that each loss factor is calculated using an AC power flow solution, {vi

, θi: i ∈ N }

and solving a square N × N system of equations. Calculating all loss factors amounts to inverting

the matrix (A>S), which can be done using direct methods in O(N4). Iterative methods such as

Jacobi or Gauss-Seidel or descent methods like conjugate gradient will likely be much faster.

This derivation can be succinctly described in vector and matrix form:

∆L = L∆Θ,

∆P = A>S∆Θ,

1

>∆L = 1>L



A>S

−1

∆P,

where L is a K ×N marginal loss matrix defined by (11), A is a K ×N network-incidence matrix, S

is a K × N marginal branch flow matrix defined by (13), and 1 is a vector of K ones. ∆L, ∆P and

∆Θ are respectively vectors of marginal changes in line losses (K × 1), nodal real power injections

(N × 1), and marginal changes in voltage angle (N × 1).

The loss factors can then be calculated as follows:

LF > = 1

>L



A>S

−1

[I − W1

>]. (16)

Then, we can calculate a constant `

0

to compensate the marginal term such that the following

equation gives a linear approximation of line losses ` that is exact at the base point solution,

` = `

0 + LF > (Pg − Pd), (17)

where ` is the total real power losses in the network, `

0

is a constant, LF is the vector of loss factors

`fi, and Pg is the vector of net real power generation, and Pd is the vector of real power demands.

Because (17) is a linear equation, it can be easily integrated into electricity market optimization

software.

13
2.3 Alternative Line Loss Derivation

Alternatively, a set of loss factors can be derived using a common method ([44, 45, 46, 51, 53, 54, 58])

which assumes that all voltages are equal to 1 and approximates,

cos(θi − θj ) ≈ 1 −

(θi − θj )

2

2

. (18)

After some substitution will yield a losses as a quadratic function of pijk times line resistance rk,

` =

X

k

rkp

2

ijk. (19)

This approximation is in fact the second order Taylor series expansion of the cosine function at

θi − θj = 0. Section 5 expands on this by generalizing the base point to values other than zero.

In the DC approximation, line flows pijk are determined using shift factors t(ijk,n)that relate

nodal injections to real power flow. This equation takes the form

pijk =

X

n

t(ijk,n)pn (20)

Let pn be the net injections at the base point solution, and pijk be calculated by substituting pn

into (20). First, we take a first order Taylor series of (19) at pijk:

` =

X

k



2rkpijkpijk − rijkp

2

ijk

. (21)

Then substitute (20):

` =

X

k

 

2rkpijkX

n

t(ijk,n)pn − rijkp

2

ijk!

=

X

n

X

k



2rkpijkt(ijk,n)



pn −

X

k

rijkp

2

ijk (22)

This gives the loss function in the same form as (17), so the loss factors can now be defined as:

`fn := X

k



2rkpijkt(ijk,n)



(23)

The constant `

0

can also be defined as the second term in (22) would imply, but a more accurate

method is to equate (17) to the actual losses at the base point, if known.

This approximation loses some fidelity compared to (16) due to the assumption on voltages and

the cosine approximation (18).

3 Model Formulation

To formulate the model, we start from the model from Litvinov [48], used by ISO-NE. Dual variables

are indicated by [·].

min c(Pg) (24a)

s.t. 1

>(Pg − Pd) = ` [λ] (24b)

` = `

0 + LF >(Pg − Pd) [σ] (24c)

T(Pg − Pd − D`) ≤ Tmax [μ] (24d)

Pmin ≤ Pg ≤ Pmax [αmin, αmax] (24e)

14
where c(·) is a linear or convex cost function; the decision variables are a vector of power generation

injections Pg and total system losses `; parameters are a vector of power demand withdrawals Pd,

loss function constant `

0

, loss factor vector LF, loss distribution factor vector D, PTDF matrix

T, transmission limit vector Tmax, generation output limit vectors Pmin and Pmax, and the dual

variables are for the power balance constraint λ, the loss function constraint σ, the transmission

constraint μ, and the generation output limit constraints αmin and αmax.

The only part of the model that has not been discussed at this point is the loss distribution factor

D, an (N × 1) vector that allocates line losses into nodal withdrawals. D is normalized to one, i.e.,

1

>D = 1. In later sections, we choose to set the elements of D for each bus to be proportional to

the line losses in the branches connected to that bus.

Excluding the term D` from (24d) causes the optimal solution to (24) to be dependent on an

arbitrarily chosen reference bus. PTDFs require the selection of a reference bus which is assumed to

be the source (or sink) of all power consumed (or produced). Power flows to and from the reference

bus are “summed” in constraint (24d), and the effect of the reference bus gets canceled out in a

lossless model. When losses are included, total injections are greater than total withdrawals, so

there will be more power flowing into the reference bus than out. Loss distribution factors fix this by

adding additional power withdrawals in the approximate location of the line losses. Consequently,

total injections will equal total withdrawals, so the solution becomes independent of the reference

bus. See reference [48] for more detail.

Excluding the term D` also causes a violation of Kirchoff’s Current Law. An example demonstrating this is in the next section.

LMPs are obtained by considering the dual problem of (24), which is given below.

max λ1

>Pd + σ



`

0 − LF >Pd



+ μ

> (Tmax + TPd) + αmin>Pmin − αmax>Pmax (25a)

s.t. λ1 + σLF + μ

>T + αmin − αmax = c [Pg] (25b)

λ + σ − μ

>TD = 0 [`] (25c)

μ, αmin, αmax ≥ 0 (25d)

Constraint (25b) forms the basis for LMPs, with its terms commonly decomposed into three

components:

λ

E := λ1, (26a)

λ

L

:= σLF, (26b)

λ

C := μ>T, (26c)

Λ := λ

E + λL + λC. (26d)

where λ

E is the marginal cost of energy at the reference bus, λL is the marginal cost of losses, and

λ

C is the marginal cost of congestion. Each component is calculated with respect to the reference

bus. Many formulations assume that λ = −σ, but in fact the economic interpretation of (25c)

shows that losses become more expensive in a congested network.

If the vector Λ gives the LMPs used to remunerate generators, then αmin and αmax are vectors

of each generator’s losses and profits, respectively.

3.1 Loss Distribution Factors and Kirchoff’s Current Law

As stated in the above, loss distribution factors are an important aspect of the DCOPF with

losses. This brief section provides an example which shows that the solution will change based on

15
which reference bus is selected when loss distribution factors are removed from the formulation.

Additionally, excluding the loss factors causes the solution to violate KCL at the reference bus.

A detailed look at how losses can be modeled in the DCOPF framework can be found in [16]. In

this example, we will use relatively intuitive values for D, setting each element proportional to the

losses on adjacent lines. The notation P

k(i)

is used to indicate a sum over the subset of branches

k ∈ K that are connected to node i.

di =

1

2

×

P

k(i)

`k

P

k

`k

New model parameters need to be calculated to perform this analysis. As mentioned previously,

the reference bus can be defined by a weighting vector W which sums to one. For example, to

select bus 1 as the reference bus, then the first element of W is one and the rest are zero. To select

a “load-weighted” reference bus, let each element of W be proportional to the load at each bus.

Reference [48] provides a simple way to update parameter values from one reference bus defined by

W to another defined by Wc.

Tb = T − TWc1

>

LFd =



LF − Wc>LF1



/



1 − Wc>LF

b`

0 = `

0

/



1 − Wc>LF

We will compare results of the model formulation (24) with results from a traditional formulation

(27), below. The later formulation is equivalent to the former with D set to zero. We solve both

formulations on the the 6-bus network from Wood and Wollenberg, available in MATPOWER [59].

The analysis was implemented in GAMS based on code available from [60]. We selecting in sequence

each of the 6 buses to be the reference bus and then finally selecting a load-weighted reference bus

(denoted ‘LW’).

min c(Pg) (27a)

s.t. 1

>(Pg − Pd) = ` [λ] (27b)

` = `

0 + LF >(Pg − Pd) [σ] (27c)

T(Pg − Pd) ≤ Tmax [μ] (27d)

Pmin ≤ Pg ≤ Pmax [αmin, αmax] (27e)

First, we look at the resulting power flow values from the transmission constraints (24d) and

(27d). As shown in Tables 3.1 and 3.2, the traditional model without loss distribution factors

distorts power flow on the network. For example, the line with the highest power flow can range

from 48.3 MW to 53.1 MW, which represents around 10% of its value. In contrast, the formulation

with loss distribution factors results in consistent power flow, i.e., the transmission constraint is

unaffected by changing the reference bus.

Secondly, we look at the solution’s adherence to KCL. KCL simply states that the current flowing

into a node is equal to the flow out. Since voltage is assumed constant in the DC power flow, here

KCL can be applied to power instead of current. KCL at each node can be checked the following

way. First, the power on each branch k is calculated.

pijk =

X

n

t(ijk,n)(p

g

n − p

d

n − dn`)

16
Table 3.1: Distorted Power Flow without Loss Distribution Factor

Branch power flow (MW), Model (27)

Branch Reference Bus

(to.from) 1 2 3 4 5 6 LW Range

1.2 -0.2 2.9 2.5 1.9 1.9 2.5 2.1 3.2

1.4 24.0 26.2 26.0 27.4 25.9 26.0 26.4 3.4

1.5 19.5 20.9 21.5 20.7 22.2 21.5 21.4 2.7

2.3 -1.7 -2.1 0.6 -1.8 -1.0 -0.4 -1.1 2.7

2.4 48.5 46.4 47.1 51.1 47.8 47.0 48.7 4.6

2.5 19.6 18.9 19.8 19.4 20.9 19.8 20.0 2.0

2.6 23.3 22.9 24.9 23.2 24.1 26.0 24.4 3.2

3.5 24.3 23.8 22.3 24.1 25.1 23.2 24.1 2.7

3.6 50.8 50.9 48.3 50.8 50.7 53.1 51.5 4.8

4.5 2.6 2.6 3.1 1.8 3.7 3.1 2.9 1.9

5.6 -4.1 -3.7 -3.2 -4.0 -4.8 -2.4 -3.7 2.4

Table 3.2: Consistent Power Flow with Loss Distribution Factor

Branch power flow (MW), Model (24)

Branch Reference Bus

(to.from) 1 2 3 4 5 6 LW Range

1.2 2.1 2.1 2.1 2.1 2.1 2.1 2.1 0.0

1.4 26.1 26.1 26.1 26.1 26.1 26.1 26.1 0.0

1.5 21.1 21.1 21.1 21.1 21.1 21.1 21.1 0.0

2.3 -1.1 -1.1 -1.1 -1.1 -1.1 -1.1 -1.1 0.0

2.4 48.0 48.0 48.0 48.0 48.0 48.0 48.0 0.0

2.5 19.7 19.7 19.7 19.7 19.7 19.7 19.7 0.0

2.6 23.9 23.9 23.9 23.9 23.9 23.9 23.9 0.0

3.5 23.8 23.8 23.8 23.8 23.8 23.8 23.8 0.0

3.6 50.6 50.6 50.6 50.6 50.6 50.6 50.6 0.0

4.5 2.8 2.8 2.8 2.8 2.8 2.8 2.8 0.0

5.6 -3.8 -3.8 -3.8 -3.8 -3.8 -3.8 -3.8 0.0

Then we calculate the “mismatch” between power flow into and out of each bus n.

KCLn = p

g

n − p

d

n −

X

k(n)

(pnjk − pink) − dn`

Since dn = 0 for all n in the traditional model, this calculation is general for both models. In the

preferred model, 24, KCL is satisfied at all nodes. In the traditional model without loss distribution

factors, 27, there is a mismatch at the reference bus which is equal to the losses approximated by

the model. The mismatch resulting from each reference bus selection (including the load-weighted

reference bus, denoted LW) is shown in Table 3.3.

In conclusion, the DCOPF with losses can give a solution that significantly distorts physical laws

if it is not formulated properly. One of these distortions is a large amount of uncertainty in the

accuracy of power flow on the transmission lines. The other related distortion is that the solution

does not satisfy KCL. When such inaccuracies are present in the dispatch model, operators are

17
Table 3.3: KCL Violations without Loss Distribution Factor.

Reference Bus

Bus 1 2 3 4 5 6 LW

1 6.7 0

2 6.7 0

3 6.7 0

4 6.7 2.2

5 6.7 2.2

6 6.7 2.2

Total 6.7 6.7 6.7 6.7 6.7 6.7 6.7

forced to operate the grid more conservatively. We believe that more accurate dispatch models will

allow operators to operate the grid closer to its physical limits and therefore more efficiently.

4 Model Comparison

This section will demonstrate the importance of initializing the OPF model with a good base point

solution when clearing an electricity market. We test three initializations of (24c) in particular:

ignoring losses, using (23), or using (16) to calculate marginal losses. Each initialization uses

progressively more information from the base point solution, a locally optimal solution to the

ACOPF. We demonstrate this on the IEEE 300-bus network from the University of Washington

test case archive [61], available in MATPOWER [59]. The analysis was implemented in GAMS

based on code available from [60].

First, we look at the canonical version of the DCOPF without an approximation for line losses.

With the distribution factor formulation, this is the same as parameterizing the model with LF =

`

0 = 0. The model therefore reduces to being “lossless”, so we compensate this inherent inaccuracy

by proportionally increasing demand to account for line losses. That is, if there are `

∗

losses in

the base point solution, then we create a new parameter for demand, Ped := Pd



1 + `

∗/1>Pd



.

This initialization will be labeled simply ‘DCOPF’ because it has been reduced to the form of the

standard ‘lossless’ DCOPF model.

Second, we initialize (24c) using the quadratic approximation (23). In this method, the dispatch

at the base point solution is used to calculate line flows, and then a quadratic approximation for

line losses is used to calculate marginal line losses. Total line losses are calculated to be equal to

actual line losses at the base point solution’s dispatch levels. This initialization will be labeled

‘DCOPF-Q’ for its use of the quadratic approximation.

Lastly, we initialize (24c) using the AC linearization (16). This parameterization uses the voltage

angles and magnitudes in the base point solution to accurately calculate marginal losses as a function

of nodal injections and withdrawals. In brief, it uses the most information from the base point

solution of the three initializations presented here. This initialization will be labeled ‘DCOPF-L’

for DCOPF with losses.

Each of the three model initializations uses the same PTDFs. Since the goal is to analyze the

model initializations effect on line loss estimation accuracy, we prefer to hold all other aspects of

the model constant. The network is not congested since we wish to discern aspects of the loss

approximation independent of network congestion.

The LP model solutions are compared to an ACOPF solution to the problem. The ACOPF uses

an exact representation of power flow in the OPF problem, and therefore is generally considered a

18
36 38 40 42 44 46 48

AC LMP($)

34

36

38

40

42

44

46

48

DC LMP($)

DCOPF (lossless)

DCOPF-Q

DCOPF-L

Figure 4.1: LMP Comparison of Linear Models

better solution. However, the ACOPF problem is non-convex and therefore is too difficult to solve

in practice and may not always find the globally optimal solution. Nonetheless, it is used here as

a benchmark for our DCOPF results.

Similar prices in the two models suggests that the linear model is a good approximation of

marginal losses, while similar objective function values suggests that the linear model gives an

efficient dispatch solution. ACOPF LMPs are the dual variable of the real power balance constraint

in an ACOPF solution [52], while DCOPF LMPs are calculated from (26).

4.1 LMP Results

Figure 4.1 shows that locational information can be an important factor in pricing. In the IEEE

300-bus example problem, we see first that prices from the ACOPF range from $37.19/MWh to

$46.76/MWh. Considering that transmission losses are only 1.2% of total demand in this example,

the total price spread may be surprising.

The most simplistic model, ‘DCOPF’, produces only a single price for each node in the system, at

$41.19/MWh. This can create inefficiencies because some generators with costs under $41.19/MWh

but a large effect on system losses will be selected to produce ahead of generators that are apparently

more expensive, but are located such that their marginal effect on losses is very low or negative.

The ‘DCOPF-Q’ model does a better job of differentiating locations based on their marginal

19
Table 4.1: IEEE 300-bus Test Case Solution Statistics

Model LMP MAPE Cost Deviation Time* (s)

DCOPF 3.77% -0.172% 0.061

DCOPF-Q 1.54% -0.114% 0.085

DCOPF-L 0.24% 0.005% 0.070

ACOPF - - 0.772

*average of ten trials

affect on losses, but it also mis-estimates the marginal affect by a large amount at some buses.

For example, the largest overestimate is at bus 7049, where the ACOPF LMP is $41.45 but the

DCOPF-Q LMP is $45.96, or about 10% higher.

The ’DCOPF-L’ model performs the best of all three linear models, producing prices that are

very similar to the ACOPF LMPs. It’s worst mis-estimation is at bus 250, where it overestimates

the price by only 3.8%.

Summary comparisons of the three models are given in Table 4.1. Solution time was measured

on a laptop computer with a 2.30 GHz processor and 8GB of RAM. Results are summarized with

two statistics defined by,

LMP MAPE = 1

N

X

i

|λ

∗

i − λ

AC

i

|

λ

AC

I

× 100%,

Cost Deviation =

c(P

∗

g

) − c(P

AC

g

)

c(P AC

g

)

× 100%.

The relative performance of the three models will obviously vary depending on the network being

studied, but our general belief is that in most cases, the three models will either perform similarly

or the DCOPF-L will perform significantly better because it can be tuned to the current operating

conditions of the network.

It should be emphasized that loss function accuracy leads to better prices and dispatch without

any additional computational cost. On the IEEE 300-bus test case, the improved methodology

decreases LMP MAPE more than sixfold compared to the quadratic approximation and more than

fifteen-fold compared to the lossless model. A similar comparison can be made with regard to

how closely the model arrives at the actual cost of dispatch. In summary, the DCOPF-Q model

performed similarly to the lossless DCOPF while the DCOPF-L model generated results more

similar to the ACOPF solution.

5 Quadratic Update Procedure

5.1 Motivating Example

The base point in the previous section was the ACOPF solution, which is optimistic to say the least.

In practice, such a good base point is not possible. The following motivating example changes the

problem so that the base point is not the optimal solution in order to show the importance of the

DCOPF’s marginal loss approximation.

Consider the two node problem described in Table 5.1. There are three generators the initially

have identical costs, and the line connecting the nodes has a specified resistance. For simplicity we

assume the voltage at both nodes is 1, so line losses are equal to p

2

12rk.

20
Table 5.1: Two Node Example

Generators

Bus Initial Bid ($) Final Bid ($) Capacity (MW)

A 1 30.00 29.50 10

B 1 30.00 29.75 100

C 2 30.00 30.00 100

Transmission Load

From To rk (Ω) Bus Amount (MW)

1 2 0.0005 2 90

Table 5.2: Solutions for Initial and Final Bids

Solution

Dispatch 1 2 3

Gen A 10 MW 10 MW 0 MW

Gen B 84.46 MW 0 MW 0 MW

Gen C 0 MW 80.05 MW 90 MW

Flow 94.46 MW 10 MW 0 MW

Losses 4.46 MW 0.05 MW 0 MW

Initial Cost $2,833.84 $2,701.50 $2,700.00

Final Cost $2,807.73 $2,696.50 $2,700.00

A few potential solutions are given in Table 5.2. Solution 3 is clearly optimal for the initial bids,

and the dispatch cost is given on the ‘Initial Cost’ line. Suppose that in the next time period,

generators A and B reduce their bids after purchasing new gas contracts on the spot market.

Instead of $30, the new bids are $29.50 for generator A and $29.75 for generator B. The new costs

are shown on the ‘Final Cost’ line of Table 5.2, and Solution 2 is optimal.

However, there is a key point that current practices miss in this scenario! Suppose that Solution

3 is used as a base point to calculate loss factors. There are no losses in the network since p12 = 0,

so loss factors at both nodes are zero. Therefore the dispatch model would select the cheapest

generators, A and B, corresponding to Solution 1. If the model were given an accurate set of loss

factors then it would have selected Solution 2. Otherwise, the actual cost of dispatch will cost

almost 4% more than the optimal solution.

5.2 Algorithm Description

We propose a sequential linear programming (SLP) method to update loss factors in such a case.

This results in a more accurate representation of marginal losses, which results in more accurate

prices and more efficient dispatch.

The core idea in the SLP methodology comes from Section 2.3 and the fact that ` =

P

k

rkp

2

ijk

gives a decent approximation for line losses. When linearized, this function splits into linear terms

(2rkpijkt(ijk,n)pn) analogous to `fn and some constant terms (rijkp

2

ijk) analogous to `

0

. Therefore,

line losses can be updated with new values p

∗

ijk each time the model is solved. If p

∗

ijk = pijk, then

the optimal solution is the same as the base point solution and the model has a good representation

of marginal line losses.

However, we have also seen that by itself, a linearization of the quadratic approximation can

result in significant pricing errors (described in Section 4). We therefore wish to combine the

21
quadratic approach with the more accurate loss factor parameterization in Section 2.2. We will

start with the more accurate parameterization and then update it as if the `fn and `

0

terms were

functions of pijk.

An update procedure assuming (22) is given in Algorithm 1. However this is just the same

quadratic loss approximation as before, which was shown to be less accurate than our proposed

method. Because it is based on a second order Taylor series expansion of the cosine terms with a

“zero” base point, this algorithm will be called the “Zero-Centered” Quadratic Update.

To find the a better quadratic approximation than this zero-centered approach, we find a second

order Taylor series expansion around the current operating point. After solving the LP, the second

order Taylor series expansion is used to update parameters in the linear loss function approximation.

The n

th order Taylor series approximation of a function f : Rn → R at a base point x ∈ R

n

can

be written,

f(x) = Xn

k=0

∇kf(x)

k!

(x − x)

k

.

Applying this to line losses as a function of power flows, `(Pt), we have the first order approximation,

`(Pt) ≈ `(Pt) + ∇`(Pt)(Pt − Pt)

= `

0 + ∇`(Pt)T(Pg − Pd).

This is a rephrasing of the methodology in Section 2.2, and clearly we have LF = ∇`(Pt)T. The

following would need to be calculated to extend this to a second order approximation,

`(Pt) ≈ `(Pt) + ∇`(Pt)(Pt − Pt) + 1

2

(Pt − Pt)

>∇2

`(Pt)(Pt − Pt).

The losses on a particular branch do not depend on the power flow across other branches, so ∇2`(Pt)

is a diagonal matrix. For each branch, we only need to compute `

00

k

(pijk).

Next, we will assume that the function takes a quadratic form. It should be similar to (19), so

we try the following,

`k = γk(pijk + ξk)

2 + ηk. (28)

Algorithm 1 Zero-Centered Quadratic Update

Require: t(ijk,n), dn, rk, p

g

n, pd

n

, `

1: pijk ←

P

n

t(ijk,n)(p

g

n − p

d

n − dn`) . ∀k ∈ K

2: ` ←

P

k

rijkp

2

ijk

3: `fn ←

P

k

2rkpijkt(ijk,n) . ∀n ∈ N

4: `

0 ← ` −

P

n

`fn(p

g

n − p

d

n

)

5: solve (24), h = 1

6: while P

i

(p

g∗

n − p

g

n)

2 ≤ tol and h ≤ hmax do

7: p

g

n ← p

g∗

n , pijk ← p

∗

ijk . ∀n ∈ N , ∀k ∈ K

8: ` ←

P

k

rijk(pijk)

2

9: `fn ←

P

k

2rkpijkt(ijk,n) . ∀n ∈ N

10: `

0 ← ` −

P

n

`fn(p

g

n − p

d

n

)

11: solve (24), h ← h + 1

12: end while

22
Any quadratic function can be defined by changing the values of γ, ξk and ηk, so this form can be

assumed without loss of generality. Unfortunately, the initial loss function does not provide enough

information to calculate all three of these coefficients. Instead, the following is derived from (7),

d`k

dpn

=

d`k

dθ

dθ

dpn

= 2gijk

vivj

aijk

sin(θi − θj − φijk)

dθ

dpn

= 2rk

vivj

aijk

1

xk

pijk

dθ

dpn

= 2rk

vivj

aijk

pijkt(ijk,n).

This gives good reason to believe that γk = rk

vivj

aijk

is a good guess, and it also happens to give

the same quadratic approximation as before if voltages are equal to their nominal values and the

turns ratio is 1. Now we proceed to find ξk and ηk. The first order linearization of (28) simplifies,

`k ≈ 2rk

vivj

aijk

(pijk + ξk)pijk + rk

vivj

aijk

(ξ

2

k − p

2

ijk) + ηk.

To put this in the same terms as (17), define `fkn, `

0

k

, `fn and `

0

,

`fkn := 2rk

vivj

aijk

(pijk + ξk)t(ijk,n)

,

`

0

k

:= rk

vivj

aijk

(ξ

2

k − p

2

ijk) + ηk,

`fn := X

k

`fkn,

`

0

:= X

k

`

0

k

.

The initial `fkn and `

0

k

are known, so solve for ξk and ηk,

ξk =

`fknaijk

2rkvivj t(ijk,n)

− pijk,

ηk = `

0

k − rk

vivj

aijk

(ξ

2

k − p

2

ijk).

Algorithm 2 implements the basic SLP, with the following few numerical side notes.

The assignment of ξk requires an arbitrary selection for the index n for `fkn and t(ijk,n). This

can be a source of numerical errors, so we choose n = arg maxm(|t(ijk,m)| : m ∈ {i, j}) to minimize

these errors.

Another numerical issue can occur when calculating ξk if γk is very small or zero due to very low

resistance on the line. In this case, we set a tolerance value εγ and let ξk = γk = 0 if γk < εγ.

In Line 6 of the algorithm in Fig. 2, any convergence criterion can be set that fits operational

needs. We’ve shown a criteria with p

g

n, but there are many other options.

Lastly, we define a damping parameter ω that changes Line 7 of the algorithm to the following

base point update rule,

p

g,h+1

n = ωp

g,h

n + (1 − ω)p

∗

g

(29)

p

h+1

ijk = ωp

h

ijk + (1 − ω)p

∗

ijk (30)

23
Each iteration in this SLP solves an approximation of a nonlinear program (NLP). This NLP

is the same formulation as (24) except that the constraint (24c) is replaced with (28). Since this

constraint is an equality instead of a greater-than-or-equal-to constraint, the problem is non-convex,

so a locally optimal solution is not guaranteed to be the globally optimal solution. However, we

would like to know if the SLP is converging to the globally optimal solution, so we solve a problem

with the following relaxation of (28),

` ≥

X

k



γk(pijk + ξk)

2 + ηk



. (31)

This relaxation makes the problem convex, and therefore any locally optimal solution is also a

globally optimal solution. Furthermore, if this constraint holds at equality in the optimal solution,

then it is also the solution to the unrelaxed problem, and this was true for all cases solved. We will

refer to this problem as the ‘DCOPF-QCP’.

5.3 Algorithm Results

We provide results for implementing Algorithm 2 on a selection of test cases from the University

of Washington test case archive [61] as well as few other that are available in MATPOWER [59].

The analysis was implemented in GAMS based on code available from [60].

Including the damping parameter ω improved the convergence speed of all test cases, and the

118- and 300-bus cases did not converge unless the damping parameter was used. After some trial

and error, ω = 0.25 for the smaller cases (<100 buses) and ω = 0.75 for the larger cases (118-

and 300-bus) showed good results. Generally, setting ω too large can slow down convergence, but

setting it too small may lead it to diverge.

The results in Fig. 5.1 were obtained by uniformly increasing demand parameters by 5% compared to the base point solution and randomizing generator costs by multiplying by a normal

random variable, N(1, 0.02). The randomization step was necessary because many of the generators have identical cost functions in the original data sets.

Convergence was measured with the standard L2 norm, defined as the square root of the sum of

squared differences. We compared P

h

g

, Ph

t and Λh with values from the previous iteration, where

the letter h denotes the h

th iterative solution to (24). Results were similar for each of Ph

g

, Ph

t and

Algorithm 2 Generic Quadratic Update

Require: t(ijk,n), dn, rk, aijk, `fi, `0, p

g

n, pd

n

, `, vi

1: pijk ←

P

n

t(ijk,n)(p

g

n − p

d

n − dn`) . ∀k ∈ K

2: γk ← rkvivj/aijk . ∀k ∈ K

3: ξk ← `fkn/2γkt(ijk,n) − pijk . ∀k ∈ K , n = arg maxm(|t(ijk,m)| : m ∈ {i, j})

4: ηk ← `

0

k − γk(ξ

2

k − p

2

ijk) . ∀k ∈ K

5: solve (24), h = 1

6: while P

i

(p

g∗

n − p

g

n)

2 ≤ tol and h ≤ hmax do

7: p

g

n ← p

g∗

n , pijk ← p

∗

ijk . ∀n ∈ N , ∀k ∈ K

8: ` ←

P

k

γk(pijk + ξk)

2 + ηk

9: `fn ← 2

P

k

γk(pijk + ξk)t(ijk,n) . ∀n ∈ N

10: `

0 ← ` −

P

n

`fn(p

g

n − p

d

n

)

11: solve (24), h ← h + 1

12: end while

24
0 1 2 3 4 5

iteration h

10-6

10-3

100

103

L2 norm

P

h

t − P

h−1

t

0 1 2 3 4 5

iteration h

10-6

10-3

100

103

L2 norm

P

h

g − P

h−1

g

0 1 2 3 4 5

iteration h

10-6

10-3

100

103

L2 norm

Λ

h − Λh−1

0 1 2 3 4 5

iteration h

 1e-06%

0.0001%

 0.01%

 1%

 100%

Absolute % Di

fference

c(P

h

g

) − c(P

QCP

g

)

case6ww

case9

case14

case24

case30

case39

case57

case118

case300

Figure 5.1: Convergence Results for Algorithm 2.

Λ

h

, so only P

h

g

is shown to save space. Fig. 5.1 also shows convergence with respect to the objective

function of the DCOPF-QCP.

In our selection of test cases, we see that all appear to converge within 0.01% of the DCOPF-QCP

solution after 20 iterations. Convergence is rapid as shown by the approximately linear results in

Figure 5.1. Comparisons between networks are difficult to interpret since the networks are different

sizes and this is not discounted in the L2 norm. However, all are converging to zero. Although

there is no guarantee for convergence, it was fairly easy to achieve the results using a very simple

damping method. Step size constraints may also be useful in larger or more complex networks, but

their use was unnecessary in our selection of test cases.

6 Conclusion

The DCOPF is at the core of many applications in today’s electricity markets. The fact that it

can be solved as an LP makes this problem formulation computationally advantageous, but this

comes at the expense of approximating the physics of power flow. Therefore, we believe that it

is important for DCOPF implementations to use the best approximations that the linear problem

formulation will allow.

Accuracy of the loss approximation should not be ignored. A feasible base point can provide

25
information about voltage angles and voltage magnitudes that are omitted from traditional DCOPF

formulations, and this additional information about the current operating point can change the

calculation for marginal line losses. In turn, changes in the calculation for marginal line losses can

have a significant effect on prices.

In addition to an accurate calculation for marginal line losses, we have discussed two other

important aspects of model accuracy. First is the inclusion of loss distribution factors in the

transmission constraint of the model formulation. Excluding loss distribution factors from the

model distorts power flows and leads to a violation of KCL in the solution. Despite this problem,

loss distribution factors are ignored in many DCOPF formulations that include losses.

Lastly, we provided an algorithm that can be used to improve the accuracy of the loss function

when the solution to the DCOPF-L is significantly different than the base point solution. This

update procedure assumes a quadratic loss function that is exact at the base point solution. The

DCOPF-L’s loss function can then be continuously updated until it converges to a solution. This

may lead to significant differences in prices when the base point solution materially differs from

the optimal solution to DCOPF-L. If this were not a common occurrence, then there would be no

reason to optimize dispatch.

The analysis presented in this paper can be of use to many researchers and practitioners interested

in modeling electricity markets. Inaccuracy of the dispatch model’s marginal terms can have a

significant effect on how much each resource is dispatched and how much they are remunerated, so

it is important to limit this inaccuracy. However, even though some inaccuracy will be present in

any model that linearizes an inherently nonlinear process, the methods explained in this paper can

be used to lessen the effect of these inaccuracies.

26
References

[1] California ISO, Business Practice Manual for Market Operations, October 2016.

https://bpmcm.caiso.com/BPM%20Document%20Library/Market%20Operations/BPM_

for_Market%20Operations_V50_clean.doc.

[2] Electric Reliability Council of Texas, ERCOT and QSE operations practices during the operating hour, June 2016. http://www.ercot.com/content/wcm/key_documents_lists/89328/

ERCOT_And_QSE_Operations_Practices_During_The_Operating_Hour_062116.doc.

[3] ISO New England, ISO New England Manual for Market Operations, Manual

M-11, May 2016. http://www.iso-ne.com/static-assets/documents/2016/05/m11_

market-operations_rev52_20160525.doc.

[4] Midcontinent ISO, Energy and Operating Reserve Markets Business Practice Manual,

March 2016. https://www.misoenergy.org/Library/BusinessPracticesManuals/Pages/

BusinessPracticesManuals.aspx.

[5] New York ISO, Day-Ahead Scheduling Manual, April 2016. http://www.nyiso.com/public/

webdocs/markets_operations/documents/Manuals_and_Guides/Manuals/Operations/

dayahd_schd_mnl.pdf.

[6] PJM, Manual 11: Energy & Ancillary Services Market Operations, July 2016. http://www.

pjm.com/documents/manuals/manuals-archive.aspx.

[7] Southwest Power Pool, Market Protocols SPP Integrated Marketplace, July 2016. https:

//www.spp.org/documents/43795/integrated%20marketplace%20protocols%2040.pdf.

[8] Calpine Corporation, Calpine & FPL Energy Request for Further Information

Prior to Considering Changes in the Allocation of Loss Over-collections, September

2003. www.iso-ne.com/committees/comm_wkgrps/mrkts_comm/mrkts/mtrls/2003/

sep16172003/A10_Calpine%20&%20FPL%20Energy%20Request_09-17-03.doc.

[9] Calpine Corporation, Current Form of Marginal Loss Pricing is Inequitable, Sends

the Wrong Market Signals and May Lead to Inefficient Outcomes, September 2003.

www.iso-ne.com/committees/comm_wkgrps/mrkts_comm/mrkts/mtrls/2003/oct222003/

A2_FPL%2520and%2520Calpine%2520Marginal%2520Loss%2520Memo_10-21-03.doc.

[10] Federal Energy Regulatory Commission, Order Assessing Civil Penalties IN15-3-000, May

2015. www.ferc.gov/CalendarFiles/20150529114118-IN15-3-000.pdf.

[11] Potomac Economics, 2011 State of the Market Report for the MISO Electricity Markets, June

2012. www.potomaceconomics.com/uploads/midwest_reports/2011_SOM_Report.pdf.

[12] Monitoring Analytics, LLC, State of the Market Report for PJM: Section 11 Congestion and

Marginal Losses, March 2015. www.monitoringanalytics.com/reports/PJM_State_of_the_

Market/2014/2014-som-pjm-volume2-sec11.pdf.

[13] H. Liu, “Regional marginal losses surplus allocation impact study,” tech. rep., California ISO,

October 2010. http://www.caiso.com/2828/2828977521d30.pdf.

[14] W. W. Hogan, Electricity Market Design Flaws and Market Manipulation, 2014. https:

//www.hks.harvard.edu/fs/whogan/Hogan_MDFMM_02_03_14.pdf.

27
[15] M. B. Cain, R. P. O’Neill, and A. Castillo, “Optimal power flow papers and formulations,”

tech. rep., Federal Energy Regulatory Commission, July 2015. https://www.ferc.gov/

industries/electric/indus-act/market-planning/opf-papers.asp.

[16] B. Stott, J. Jardim, and O. Alsa ̧c, “DC power flow revisited,” IEEE Trans. on Power Systems,

vol. 24, no. 3, pp. 1290–1300, 2009.

[17] CAISO Market Operations, Market Optimization Details, November 2009. http://www.

caiso.com/23cf/23cfe2c91d880.pdf.

[18] Electric Reliability Council of Texas, ERCOT Market Education: Transmission 101, March

2013. http://www.ercot.com/content/wcm/training_courses/101/trn101_m1_120313.

pdf.

[19] Electric Reliability Council of Texas, Transmission and Distribution Losses, December 2013.

www.ercot.com/content/mktrules/protocols/current/13-010109.doc.

[20] E. Litvinov, Locational Marginal Pricing. ISO-New England, November 2011. www.iso-ne.

com/support/training/courses/wem301/wem301_lmp.pdf.

[21] D. Chatterjee and R. Sutton, Modeling and Pricing of Losses in Market Applications. Midcontinent Independent System Operator, March 2015. www.misoenergy.org/_layouts/MISO/

ECM/Redirect.aspx?ID=195505.

[22] New York System Operator, Inc., Market Administration and Control Area Services Tariff

(MST) 17, Attachment B, July 2013. http://www.nyiso.com/public/markets_operations/

documents/tariffviewer/index.jsp.

[23] PJM Interconnection, Marginal Losses Implementation Training, November 2008. http://www.pjm.com/~/media/training/new-initiatives/ip-ml/

marginal-losses-implementation-training.ashx.

[24] J. Carpentier, “Contribution to the economic dispatch problem,” Bulletin de la Societe Francoise des Electriciens, vol. 3, no. 8, pp. 431–447, 1962.

[25] M. C. Caramanis, R. E. Bohn, and F. C. Schweppe, “Optimal spot pricing: practice and

theory,” IEEE Trans. on Power Apparatus and Systems, no. 9, pp. 3234–3245, 1982.

[26] B. H. Chowdhury and S. Rahman, “A review of recent advances in economic dispatch,” 1990.

[27] M. El-Hawary, “Optimal economic operation of large scale electric power systems: a review,”

in Athens Power Tech, 1993. APT 93. Proceedings. Joint International Power Conference,

vol. 1, pp. 206–210, IEEE, 1993.

[28] M. Huneault, F. Galiana, and Q. St Bruno, “A survey of the optrial power flow literature,”

IEEE Trans. on Power Systems, vol. 6, no. 2, p. 18, 1991.

[29] J. A. Momoh, M. El-Hawary, and R. Adapa, “A review of selected optimal power flow literature

to 1993. part i: Nonlinear and quadratic programming approaches,” IEEE Trans. on Power

Systems, vol. 14, no. 1, pp. 96–104, 1999.

[30] J. A. Momoh, M. El-Hawary, and R. Adapa, “A review of selected optimal power flow literature

to 1993. part ii: Newton, linear programming and interior point methods,” IEEE Trans. on

Power Systems, vol. 14, no. 1, pp. 105–111, 1999.

28
[31] Z. Qiu, G. Deconinck, and R. Belmans, “A literature survey of optimal power flow problems in

the electricity market context,” in Power Systems Conference and Exposition, 2009. PSCE’09.

IEEE/PES, pp. 1–6, IEEE, 2009.

[32] X. Bai, H. Wei, K. Fujisawa, and Y. Wang, “Semidefinite programming for optimal power

flow problems,” International Journal of Electrical Power & Energy Systems, vol. 30, no. 6,

pp. 383–392, 2008.

[33] B. C. Lesieutre, D. K. Molzahn, A. R. Borden, and C. L. DeMarco, “Examining the limits

of the application of semidefinite programming to power flow problems,” in Proc. of the 49th

Allerton Conf. on Communication, Control and Computing, 2011.

[34] J. Lavaei and S. H. Low, “Zero duality gap in optimal power flow problem,” IEEE Trans. on

Power Systems, vol. 27, no. 1, pp. 92–107, 2012.

[35] D. K. Molzahn, J. T. Holzer, B. C. Lesieutre, and C. L. DeMarco, “Implementation of a largescale optimal power flow solver based on semidefinite programming,” IEEE Trans. on Power

Systems, vol. 28, no. 4, pp. 3987–3998, 2013.

[36] R. A. Jabr, “Radial distribution load flow using conic programming,” IEEE Trans. on Power

Systems, vol. 21, no. 3, pp. 1458–1459, 2006.

[37] A. Castillo and R. P. O’Neill, “Survey of approaches to solving the ACOPF,” tech. rep., Federal

Energy Regulatory Commission, 2013.

[38] T. J. Overbye, X. Cheng, and Y. Sun, “A comparison of the ac and dc power flow models

for lmp calculations,” in Proceedings of the 37th Annual Hawaii International Conference on

System Sciences, 2004, pp. 9–pp, IEEE, 2004.

[39] F. Li and R. Bo, “Dcopf-based lmp simulation: algorithm, comparison with acopf, and sensitivity,” IEEE Trans. on Power Systems, vol. 22, no. 4, pp. 1475–1485, 2007.

[40] D. Wells, “Method for economic secure loading of a power system,” Proceedings of the Institution of Electrical Engineers, vol. 115, no. 8, pp. 1190–1194, 1968.

[41] B. Stott and J. Marinho, “Linear programming for power-system network security applications,” IEEE Trans. on Power Apparatus and Systems, vol. 3, no. PAS-98, pp. 837–848, 1979.

[42] R. R. Shoults, W. M. Grady, and S. Helmick, “An efficient method for computing loss formula

coefficients based upon the method of least squares,” IEEE Trans. on Power Apparatus and

Systems, no. 6, pp. 2144–2152, 1979.

[43] O. Alsac, J. Bright, M. Prais, and B. Stott, “Further developments in LP-based optimal power

flow,” IEEE Trans. on Power Systems, vol. 5, no. 3, pp. 697–711, 1990.

[44] Y.-C. Chang, W.-T. Yang, and C.-C. Liu, “A new method for calculating loss coefficients [of

power systems],” IEEE Trans. on Power Systems, vol. 9, no. 3, pp. 1665–1671, 1994.

[45] T. N. D. Santos and A. L. Diniz, “A dynamic piecewise linear model for dc transmission losses

in optimal scheduling problems,” IEEE Trans. on Power Systems, vol. 26, no. 2, pp. 508–519,

2011.

[46] V. N. Bharatwaj, A. Abhyankar, and P. Bijwe, “Iterative dcopf model using distributed slack

bus,” in 2012 IEEE Power and Energy Society General Meeting, pp. 1–7, IEEE, 2012.

29
[47] Z. Hu, H. Cheng, Z. Yan, and F. Li, “An iterative lmp calculation method considering loss

distributions,” IEEE Trans. on Power Systems, vol. 25, no. 3, pp. 1469–1477, 2010.

[48] E. Litvinov, T. Zheng, G. Rosenwald, and P. Shamsollahi, “Marginal loss modeling in lmp

calculation,” IEEE Trans. on Power Systems, vol. 19, no. 2, pp. 880–888, 2004.

[49] T. Orfanogianni and G. Gross, “A general formulation for lmp evaluation,” IEEE Trans. on

Power Systems, vol. 22, no. 3, pp. 1163–1173, 2007.

[50] J.-C. Peng, H. Jiang, G. Xu, A. Luo, and C. Huang, “Independent marginal losses with application to locational marginal price calculation,” IET generation, transmission & distribution,

vol. 3, no. 7, pp. 679–689, 2009.

[51] F. Li, “Fully reference-independent LMP decomposition using reference-independent loss factors,” Electric Power Systems Research, vol. 81, no. 11, pp. 1995–2004, 2011.

[52] H. Liu, L. Tesfatsion, and A. Chowdhury, “Locational marginal pricing basics for restructured

wholesale power markets,” in 2009 IEEE Power & Energy Society General Meeting, pp. 1–8,

IEEE, 2009.

[53] B. F. Hobbs, G. Drayton, E. B. Fisher, and W. Lise, “Improved transmission representations

in oligopolistic market models: quadratic losses, phase shifters, and DC lines,” IEEE Trans.

on Power Systems, vol. 23, no. 3, pp. 1018–1029, 2008.

[54] D. Z. Fitiwi, L. Olmos, M. Rivier, F. de Cuadra, and I. P ́erez-Arriaga, “Finding a representative

network losses model for large-scale transmission expansion planning with renewable energy

sources,” Energy, vol. 101, pp. 343–358, 2016.

[55] A. Castillo, X. Jiang, and D. Gayme, “Lossy dcopf for optimizing congested grids with renewable energy and storage,” in American Control Conference (ACC), 2014, pp. 4342–4347,

IEEE, 2014.

[56] R. Jabr et al., “Modeling network losses using quadratic cones,” IEEE Trans. on Power Systems, vol. 20, no. 1, pp. 505–506, 2005.

[57] H. Ronellenfitsch, D. Manik, J. H ̈orsch, T. Brown, and D. Witthaut, “Dual theory of transmission line outages,” arXiv preprint arXiv:1606.07276, 2016.

[58] F. C. Schweppe, M. C. Caramanis, R. D. Tabors, and R. E. Bohn, Spot pricing of electricity.

Kluwer Academic Publishers, 1988.

[59] R. D. Zimmerman, C. E. Murillo-S ́anchez, and R. J. Thomas, “MATPOWER: Steady-state

operations, planning and analysis tools for power systems research and education,” IEEE

Trans. on Power Systems, vol. 26, no. 1, pp. 12–19, 2011.

[60] L. Tang and M. Ferris, “Collection of power flow models: Mathematical formulations,” 2015.

http://www.neos-guide.org/sites/default/files/math_formulation.pdf.

[61] “Power systems test case archive,” tech. rep., University of Washington, 1999. http://www.

ee.washington.edu/research/pstca.

30