Electrical Power and Energy Systems 32 (2010) 329–332

Contents lists available at ScienceDirect

Electrical Power and Energy Systems

journal homepage: www.elsevier.com/locate/ijepes

Short Communication

Tracing the flows of electricity

Paulo Barcia a,*, Rui Pestana b
a
Universidade Nova de Lisboa, Faculdade de Economia, Campus de Campolide, 1099-032 Lisboa, Portugal
b
REN, Rede Eléctrica Nacional, S.A. Rua Cidade de Goa 4, 2685-039 Sacavem, Portugal

a r t i c l e i n f o a b s t r a c t

Article history: Consider an electrical power system and let G denote the set of nodes with net generation, L the set of
Received 3 October 2008 nodes with a net load and assume the DC model of the network. Our aim is to determine how much
Received in revised form 27 July 2009 power is each node of G sending to each node of L in such a way that the flow decomposition arising from
Accepted 25 September 2009
this assignment is such that each branch flow is broken in parcels of the (as much as possible) same ori-
entation in order to provide a sensible tool for transmission cost allocation.
Ó 2009 Elsevier Ltd. All rights reserved.
Keywords:
Power systems
DC model
Transmission cost allocation
Linear programming
Transportation problem

1. Introduction very large networks, such as the full European grid, we give a fast
heuristic capable of finding good quality solutions.
Tracing the flows of electricity is of paramount importance for In what follows the DC model of the network will be assumed
transmission cost allocation. To this date European UCTE members throughout.
still debate how transmission costs should be properly allocated to
its different partners. This paper is a contribution to this ongoing
discussion. 2. The model
Several approaches to this problem have been attempted in the
literature: Bialek [1] uses proportional sharing, Usaola [2] suggests The key ingredient for determining a realistic transmission cost
the differential method, Visakha et al. [3] employ a special type of allocation in a power grid is to establish how much power is each
electrical distances while Oloomi-Buygi et al. [4] consider system generator sending to each load in the network. In fact, once that
non-linearity. Camacho et al. [5] provide an interesting assessment information is made available, the linearity of the DC model allows
of the current inter-TSO compensation algorithms in the internal do decompose every flow on the network in ‘‘partial flows” accord-
electricity market of the European Union. ing to each pair generation to load and therefore to produce a sen-
Here we assume, as an ‘‘economy principle” of nature, that elec- sible transmission cost allocation.
tricity flows through paths that minimize the total MW km cov- Consider the example in Fig. 1 of a simple 4-bus network with
ered in the entire system in order to obtain a linear two generators and two loads. Fig. 1 shows the generators, loads
programming model that may be useful for large networks. We and flows observed as well as the lengths of the branches of
introduce the notion of ‘‘power distances” which, if used within network.
the linear programming model, provides an allocation of genera- If one asks how much power is each generator sending to each
tion to loads that induces a flow decomposition such that each load several answers are compatible with the flows that are ob-
branch flow is broken in parcels of the (as much as possible) same served on the network.
orientation. This feature and the fact that such a model runs on Let pij denote the amount of power sent from bus i to bus j.
moderate computing time for large networks such as the UCTE sys- A first natural solution is as follows: p21 = 1000, originating the
tem makes it a sensible tool for transmission cost allocation. For ‘‘partial flows” shown in Fig. 2; p23 = 0, p41 = 100 originating the
‘‘partial flows” in Fig. 3; and p43 = 1000 producing the ‘‘partial
flows” shown in Fig. 4.
Adding the diagrams corresponding to Fig. 2, Fig. 3 and Fig. 4
* Corresponding author. Tel.: +351 213 801 600; fax: +351 213 870 933. one obtains the original flows generations and loads. Note also that
E-mail addresses: barcia@fe.unl.pt (P. Barcia), rui.pestana@ren.pt (R. Pestana). all the ‘‘partial flows” thus obtained are flowing in the same sense

0142-0615/$ - see front matter Ó 2009 Elsevier Ltd. All rights reserved.
doi:10.1016/j.ijepes.2009.09.009
330 P. Barcia, R. Pestana / Electrical Power and Energy Systems 32 (2010) 329–332

1000 MW 1100 MW 100 MW

1100 MW
~ 100 MW ~ 100 MW
~ 100 MW ~
1100 MW 100 MW
BUS 1 BUS 2 BUS 3 BUS 4 BUS 1 BUS 2 BUS 3 BUS 4
1 km 1000 km 1 km 1 km 1000 km 1 km

1100 MW 1000 MW 100 MW

Fig. 1. Fig. 3.

1000 MW
1000 MW

1000 MW
~ ~
~ ~
1000 MW
BUS 1 BUS 2 BUS 3 BUS 4
1 km 1000 km 1 km BUS 1 BUS 2 BUS 3 BUS 4
1 km 1000 km 1 km

1000 MW
1000 MW
Fig. 2.
Fig. 4.

of the original flows suggesting that this solution is a natural way

of assigning generation to loads for this network and therefore
decomposing each branch flow according to each pair generation
to load allowing for transmission cost allocation. 0.4 MW
0.4 MW 0.4 MW
There are, of course, infinitely many solutions to this problem. A 1 km
1 km 1 km
second possible solution is p21 = 0, p23 = 1000, p41 = 1100 and 1 MW i j 1 MW
p43 = 0. Remark that, for this second solution, the flow of 0.6 MW 0.6 MW
100 MW from bus 3 to bus 2 is obtained adding 1100 MW to
2 km 2 km
1000 MW flowing in the opposite direction. We would like to avoid
this type of unnatural solutions. To do so we shall assume, as a first Fig. 5.
step, that generation is assigned to loads in such a way that the to-
tal MW km covered in the network is minimized.
Consider an arbitrary network with a set G of generators and a
set L of loads. Let gi denote the generation at node i of G, lj the load ing to determine a ‘‘proper” distance dij from bus i to bus j. A first
at node j of L and dij a distance (this distance can be a geographical obvious approach would be to take dij = 3 km, the shortest distance
distance or something else as we shall see later on) from node i of G from bus i to bus j. However electric power has the bad habit of not
to node j of L. taking the shortest path, instead it flows through paths that are
Our first step will be to assume, as an ‘‘economy principle” of determined by the network branch reactances. Assume we are
nature, that generations are assigned to loads according to the fol- sending 1 MW from bus i to bus j and that the reactances are such
lowing model: that 0.4 MW follow the shortest path and 0.6 MW the longest one.
X In this case 0.4  1 + 0.4  1 + 0.4  1 + 0.6  2 + 0.6  2 =
Min dij pij 3.6 MW km will be covered to send 1 MW from bus i to bus j.
ij
X Therefore it is reasonable to state that the ‘‘power distance” from
pij ¼ lj for each node j of L i to j is dij = 3.6 ‘‘power kilometres”.
X
i We now give a precise definition of the power distance dij. Let
pij ¼ g i for each node i of G dFkt/dpij be the flow on branch kt originated by sending 1 MW from
j bus i to bus j. The value of dFkt/dpij can be easily computed from the
pij P 0 for every i of G and every j of L DC equations of the network (see Appendix). If lkt denotes the
length of branch kt the power distance dij is defined as
This is a classical transportation model that can be successfully X
solved for very large networks. Solving this model using geograph- dij ¼ jdF kt =dpij jlkt
kt
ical distances for the previous example one would, of course, get the
first solution as an optimal one. where the summation is taken over all the branches kt of the
network.
We are now ready to discuss two important properties of power
3. Power distances
distances.
In this section we define distances suitable for our purpose. Proposition 3.1. The power distance dij is non-negative, symmetrical
Consider the example shown in Fig. 5 and suppose that we are try- and the triangle inequality also holds.
 P. Barcia, R. Pestana / Electrical Power and Energy Systems 32 (2010) 329–332 331

Proposition 3.1 states that dij is a distance in the mathematical The quality of this heuristic was accessed by solving the full
sense. This is of paramount importance if dij is to be used in algo- network transportation model and the heuristic whenever both
rithms that have such a requirement. solutions were available. The heuristic solutions obtained in the
It is natural to assume that generations are assigned to loads situations studied so far were proven to be at most 2% from
minimizing the total MW km covered in the entire network. How- optimality.
ever the following proposition gives a critical insight on the utiliza- Once an assignment of generations to loads is available solving
tion of power distances in our model: DC power flows for each pair generation to load will do the tracing
of electricity flows therefore allowing for transmission cost
Proposition 3.2. Using power distances in the transportation model
allocation.
of the previous section is equivalent to minimize the total amount of
A code implementing this algorithm can be obtained from the
MW km covered by the ‘‘partial flows” flowing opposite to each branch
authors.
flow.
Proof of Proposition 3.1. It follows directly from the definition of
Proposition 3.2 implies that using power distances in the trans- ‘‘power distances” that it is obviously positive (dij P 0) and
portation model will produce an allocation of generation to loads symmetrical (dij = dji). It is easy to verify that the triangle inequality
that most naturally fits the branch flows observed on the network (dij 5 dil + dlj) also holds. The formula for dFkt/dpij implies:
and therefore allowing to decompose each flow in fragments corre-
sponding to generation to load pairs (the ‘‘partial flows”) having (as jdF kt =dpij j ¼ jAkt;i  Akt;j j ¼ jðAkt;i  Akt;l Þ þ ðAkt;l  Akt;j Þj
much as possible) the same orientation. ¼ jdF kt =dpil þ dpkt =dplj j

4. Solving the real thing Therefore we have jdFkt/dpijj 5 jdFkt/dpilj + jdFkt/dpljj and the triangle
inequality holds. h
We implemented our model successfully, using the routine RE-
LAXT in [6], to solve real life situations for the UCTE network of Proof of Proposition 3.2. The proof of Proposition 3.2 relies on a
2005 and 2006 with moderate computing times of about 10 min- simple technical lemma that we state and prove now.
utes on an Intel Core Duo laptop.
When addressing the full European network, of about 10,000 Lemma. Let ai, i 2 I, be a finite set of real numbers and denote by
nodes and 15,000 branches, a heuristic may be needed in order O  I the index set of the numbers ai with a sign (zero is taken as posi-
P
to obtain a feasible solution for the model. tive) opposite to the sign of their sum i2Iai. The following equality
This heuristic solves the problem for each European country holds:
 
separately and manages to glue together the partial solutions in or- X  X X
 
der to obtain a feasible solution for the full network. This is done as  a¼ jai j  2 jai j
 i2I i  i2I i2O
follows.
In the ‘‘middle” of each interconnection branch we place an ex-
tra node called an x-node. For each country each x-node is viewed Proof. Let P = {i 2 I: ai P 0} and N = {i 2 I: ai < 0}, i.e., P is the index
as a load (generation) if the corresponding branch is exporting set of the positive numbers and N the index set of the negative ones.
(importing) with a value corresponding to the interconnection P
If i2Iai P 0 then O = N and we have
branch flow.  
X  X X X X X
For each country we solve a ‘‘pure transportation” model using  
 a¼ ai ¼ ai þ ai ¼ jai j  jai j
RELAXT in [6]. Once this is done loop flows between the x-nodes in  i2I i  i2I i2P i2N i2P i2N
the full network, if any, can be removed. This makes the x-graph X X X X X
¼ jai j þ jai j  2 jai j ¼ jai j  2 jai j
(the digraph of flows amongst x-nodes) cycle free, therefore it will
i2P i2N i2N i2I i2N
have (at least) one x-node with in-degree zero. This means that X X
this x-node will have no generator x-node assigned to it but it will ¼ jai j  2 jai j:
i2I i2O
be assigned to some loads, maybe other x-nodes, as suggested in
P
Fig. 6 below. If i2Iai < 0 then the lemma follows from the previous reasoning
We can then solve a transportation problem for that x-node applied to the numbers ai. h
assigning generators to loads (possibly other x-nodes) between We are now ready to prove Proposition 3.2.
the two different countries. After that is done the previous x- Recall that the objective function of the model in Section 2 is
P
node with zero in-degree can be removed from the x-graph W = ijdijpij.
and we still will have a cycle free digraph. Therefore we can find If power distances are used we have:
another x-node with in-degree zero in the remaining x-graph ! !
X X  X X  
and iterate the procedure until all x-nodes are solved thus W¼ pij dF kt =dpij lkt ¼ lkt 
pij dF kt =dpij 
obtaining an assignment of generations to loads for the full ij kt kt ij
network.
where lkt denotes the length of branch kt.
P
Recall that Fkt = ijpijdFkt/dpij is the flow on branch kt and
remark that the total MW km covered in the entire network
P
X = ktlktjFktj is a constant because generations gi(i2G) and loads lj
(j2L) are kept fixed in our model.
Denote by Okt the index set of the pairs ij (generation to load) for
which the corresponding ‘‘partial flow” pijdFkt/dpij on branch kt has a
sign opposite to the sign of total flow Fkt on that branch. The equality
 
X  X   X  
 
jF kt j ¼  pij dF kt =dpij  ¼ pij dF kt =dpij   2 pij dF kt =dpij :
 ij  ij ij2Okt
Fig. 6.
332 P. Barcia, R. Pestana / Electrical Power and Energy Systems 32 (2010) 329–332

must hold since it is the result of having pij P 0 and of applying the where we set href = 0)h = Z(g  l) where Z = Y1 denotes the network
lemma to the numbers pijdFkt/dpij. impedance matrix.
We have It is well known from the DC model theory (see for instance
Kirschen et al. [7]) that the sensitivity coefficient Akt,i measuring
X X X X X the impact on the flow of branch kt of a power injection at node
X¼ lkt jF kt j ¼ lkt pij jdF kt =dpij j  2 lkt pij jdF kt =dpij j
i is given by Akt, i = (Zki  Zti)/xkt where xkt denotes the reactance
kt kt ij kt ij2Okt
of the branch kt. In this expression Zrs denotes the element of line
or r and column s of Z with the caveat that if either r or s is the index
X X ref of the reference node we set Zrs = 0.
X¼W2 lkt pij jdF kt =dpij j Now suppose that we are ‘‘sending” an extra 1 MW from gener-
kt ij2Okt ator i of G to load j of L and we ask what happens to the flow Fkt on
P P each branch kt, i.e., what is the value of the derivative dFkt/dpij,
Therefore minimizing W = X + 2 kt lkt ij2Okt pijjdFkt/dpijj is
P P where pij denotes the amount of power sent from node i of G to
equivalent to minimize ktlkt ij2OktpijjdFkt/dpijj the total amount
of MW km covered by the ‘‘partial flows” flowing opposite to each node j of L. Clearly we must have
branch flow. h dF kt =dpij ¼ Akt;i  Akt;j ¼ ½ðZ ki  Z ti Þ  ðZ kj  Z tj Þ=xkt

Acknowledgements an expression which is easily computable from the network branch

reactances and where the previous caveat also applies.
We are grateful to the editor and two anonymous referees for
many comments and suggestions that helped to improve this pa- References
per. We are also thankful to Prof. Teresa Correia de Barros (Lisbon
[1] Bialek J. Tracing the flows of electricity. Proc Gener Transm Distrib
Technical University) for suggesting the term ‘‘power distance”. 1996;143(4):313–20.
[2] Usaola J. A transaction-based method for allocation of transmission grid cost
and losses. Electr Power Syst Res 2006;76:395–403.
Appendix. Computing dFkt/dpij from the DC model [3] Visakha K, Thukaram D, Jenkins L. Transmission charges of power contracts
based on relative electrical distances in open access. Electr Power Syst Res
We take the DC equations of the network Yh = g  l, where Y de- 2004;70:153–61.
[4] Oloomi-Buygi M, Salehizadeh M. Considering system non-linearity in
notes the network admittance matrix, h the vector of the node volt-
transmission pricing. Int J Electr Power Energy Syst 2008;30:455–61.
age angles and g  l the net power injection vector, i.e., generation [5] Camacho L, Perez-Arriaga I. An assessment of inter-TSO compensation
minus load at each node. Of course we must have total generation algorithms in the internal electricity market of the European Union. Int J
Electr Power Energy Syst 2007;29:699–712.
equal to total load since the DC model allows for no losses.
[6] Bertsekas D, Tseng P. The relax codes for minimum cost network flow. Ann
The network equations Yh = g  l imply (after removing, as Operat Res 1988;13(1–4):125–92.
usual, the line and column corresponding to a reference node ref [7] Kirschen D, Strbac G. Fundamentals of power systems economics. Wiley; 2004.