See discussions, stats, and author profiles for this publication at: https://www.researchgate.

net/publication/315705432

Solution of DC Railway Traction Power Flow Systems Including Limited Network

Receptivity

Article in Power Systems, IEEE Transactions on · March 2017

DOI: 10.1109/TPWRS.2017.2688338

CITATIONS READS
55 1,069

2 authors, including:

Izudin Džafić
University of Sarajevo
127 PUBLICATIONS 2,051 CITATIONS

SEE PROFILE

All content following this page was uploaded by Izudin Džafić on 15 December 2019.

The user has requested enhancement of the downloaded file.

 1

This is a preprint of an article whose final and definitive form was published in IEEE Transactions on Power Systems:
R. A. Jabr and I. Džafić, “Solution of DC railway traction power flow systems including limited network receptivity,” IEEE
Transactions on Power Systems, vol. 33, no. 1, pp. 962– 969, Jan 2018.
Link: https://ieeexplore.ieee.org/document/7888589
DOI: 10.1109/TPWRS.2017.2688338
 2

Solution of DC Railway Traction Power Flow

Systems including Limited Network Receptivity
Rabih A. Jabr, Fellow, IEEE, and Izudin Džafić, Senior Member, IEEE

Abstract
The voltage solution of DC railway traction power networks is classically obtained via the current injection (CI) method,
which is based on solving a sequence of nodal voltage equations. Specialized techniques, which build on the CI method, have
been proposed for simulating limited network receptivity due to voltage rise constraints and nonreversible substations. These
techniques may require a multitude of power reduction steps for modeling the local controller operation of trains in regenerative
braking mode, and they consequently lead to a large computational effort. This paper proposes a sensitivity based approach for
computing the regenerative train power that can be returned to the network without causing over voltage. In case of non-receptive
substations, each regenerating train is switched to a voltage-current source model and the CI method is used to further adjust the
power that can be received by the network; a two-phase approach is used to compute the regenerative train resistance without
recourse to iterations. The proposed method is tested on network models with branched lines, detailed return circuits, and having
up to 144 trains. The computational performance comparisons show that the proposed method for simulating local controllers can
be significantly faster than the classical power reduction method.

Index Terms
Load flow analysis, load flow control, power control, traction power supplies, voltage control.

N OMENCLATURE
cj Index of catenary node with regenerating train connection and voltage higher than the limit.
Ggk Element in row g and column k of the modified conductance matrix.
[H] Sensitivity matrix.
Igk Current directed from node g to node k.
Ik Current injection at node k.
m Iteration counter in the CI method.
N Total number of nodes.
PT k Train power defined according to the active sign convention.
Rk Resistance of the regenerative train in the voltage-current source model at catenary node k.
Rgk Element in row g and column k of the modified resistance matrix.
Vk Voltage at node k.
Vjmax Maximum permissible voltage at node j.
Vs Substation voltage.
∆PT j Computed change in power injection at catenary node j.
∆Vj Deviation between the maximum voltage magnitude and its computed value at catenary node j.

I. I NTRODUCTION
ESIGNING a DC rail traction power system or modifying an existing one is achieved by computer-aided techniques based
D on simulation. Traction simulators are employed to establish equipment ratings, to assess train regeneration effects, and to
verify that catenary voltages remain within permissible limits and that rail potential and earth leakage currents are kept according
to specifications for safety and interference mitigation [1], [2]. Traction simulators are based on power flow calculations carried
out at different time instants sampled from the study horizon; the locations and power demands of the vehicles vary over the
study period, and each instant has the trains represented using the stationary equivalents for moving vehicular loads [3]. A
typical simulation study may require hundreds of sampled instants and consequently a significant computational effort; this
effort increases with the sampling rate employed over the study horizon [4]. The computational effort is further exacerbated
with the simulation of regenerative train breaking and nonreversible substations, as this requires additional power flow iterations
for modeling the behavior of the train local controllers [5]. This paper proposes the use of a sensitivity approach as a means
to enhance computational performance when modeling local controllers in DC traction power flow.
R. A. Jabr is with the Department of Electrical & Computer Engineering, American University of Beirut, P.O. Box 11-0236, Riad El- Solh / Beirut 1107
2020, Lebanon (email: rabih.jabr@aub.edu.lb).
I. Džafić is with the International University of Sarajevo, Hrasnička cesta 15, 71210 Sarajevo, Bosnia (email: idzafic@ieee.org).
 3

A. Related Work
Transport via railways is increasing, thus creating a need for techniques to study the expansion planning and investment in
railway power systems [6], [7]. DC traction power flow is employed in studies for supplying railway traction power equipment
or for supporting tendering and contracting activities [8]. The simulators used by the industrial companies commonly employ
the conductance matrix iterative approach with the trains represented either as conductances or as current sources [1], [2],
[8]; a comparison reported in [5] showed that the current source train model is superior in terms of convergence speed. The
conductance matrix iterative approach using the current source train model is known as the current injection (CI) power flow.
The conductance matrix that arises in the CI method is a sparse representation of the nodal voltage equations for a ladder
type network; this matrix is large because it includes rail nodes that turn out in finite cell modeling of return circuits [9], but
typically more than 80 percent of its elements are zeros [10]. The solution of the nodal equations at each iteration of the CI
method can be carried out using direct matrix factorization methods [10], or indirect methods such as incomplete Cholesky
conjugate gradient [11]. Other techniques for traction power flow include Newton-Raphson type methods [4], optimal power
flow solutions [12], and the forward-backward sweep power flow [13]. The DC traction power flow has been also extended to
model the AC distribution network [12], [14], [15], and onboard train storage [16], [17]. The simulation of trains that utilize
combined regenerative and rheostatic braking requires the modeling of local controllers that limit the maximum pantograph
voltage; a survey of the related literature [1], [2], [4], [5], [10]–[17] shows that this type of control has been reported in [5],
[12], [13], [16]. The local voltage regulation employed during braking so that the catenary voltage does not exceed a prescribed
limit is referred to as squeeze control [13]; squeeze control reduces the regenerated power that is returned to the network by
diverting part of it to the rheostatic braking system, with the aim of bringing the catenary voltage within its preset limit.
However, limited network receptivity is not only attributed to voltage rise limitations, but also to nonreversible substations; in
this respect, only [5] simulates non-receptive substations while refs. [12], [13], [16] all consider DC transit systems that are
fitted with inverting substations, as described in [18]–[20]. Both the CI method [5] and the forward-backward sweep technique
[13] have been also investigated for the power flow solution of AC distribution networks [21], [22], and the superiority of the
CI method for handling meshes is well understood.

B. Technical Contribution
This paper builds on the CI method and the associated techniques [5] to simulate the local control of the regenerative chopper
in blended regenerative/rheostatic braking so as to handle voltage rise constraints and non-receptive supplies. It presents two
contributions that aim to improve the computational performance of the power flow as compared to [5]. These contributions
are:
• a sensitivity method to handle local controllers that is competitive with the power reduction and balance method in [5].
A closed-form solution for the elements of the sensitivity matrix is derived.
• a hybrid voltage-current source model to reset the pantograph voltage of regenerating trains when substations are reverse
biased. As compared to [5], a two-phase simulation approach is proposed for computing the regenerative train resistance
without iteration; this extends the voltage source model in [5] to effectively handle non-receptive substations in networks
with detailed models of return circuits.
This study also brings attention to the analysis of this special type of power networks with moving vehicular loads, at the time
the electrification of railway systems is becoming increasingly important [6], [7].

C. Paper Structure
The rest of this paper is organized as follows. Section II presents the formulation of the DC traction power flow problem
and its solution using the CI method. The main contributions of the paper, i.e. the simulation of local controllers and the two
phase approach for handling non-receptive substations, are described in Section III. Section IV includes numerical results on
four DC traction networks with branched lines, and a comparative analysis with the approach in [5]. The paper is concluded
in Section V.

II. DC R AILWAY T RACTION P OWER F LOW S OLUTION

Consider a DC railway traction network having N nodes, NS substations, and NT trains. The electrical performance of the
network can be obtained from a sequence of power flow solutions representing different time instants over the study horizon,
with the trains frozen in time for each instant; this approach is referred to as the stationary equivalent method for moving loads
[3]. Fig. 1 shows a portion of a DC railway traction system forming a ladder type network. Each train is connected electrically
between the catenary (or third rail) and the track, whose running rails are used for traction current return. The earth leakage
currents are modeled to flow through a ladder network composed of finite cells, each represented by an equivalent π-circuit
[9]. The substation is connected between the positive busbar p and the negative busbar n, which itself is grounded and also
connected to the track for current return. Consider a train connected between a catenary node k and track node g; the train
power is related to the node voltages and the current Igk in the train via:
 4

PT k = (Vk − Vg )Igk (1)

where PT k is defined according to the active sign convention, i.e. PT k < 0 indicates that the train is absorbing power in the
motoring mode and PT k > 0 is for the train supplying power in the regenerative braking mode. Fig. 2 shows that the train
can be represented via two current source injections, the first at the catenary node k and the second at the track node g:
PT k
Ik = Igk = (2)
Vk − Vg
−PT k
Ig = −Igk = (3)
Vk − Vg
A substation connected between the positive busbar p and the negative busbar n stipulates two equations:

Vp − Vn = Vs (4)

Ip + In = 0 (5)
where (5) is effectively a super-node equation in which the sum of current injections into nodes p and n is set to zero; this is
depicted in Fig. 3. The network nodal voltages are related to the nodal injection currents through the modified conductance
matrix [23]:
    
1 −1 0 ··· ··· 0 V1 Vs
 .. . . .. .. .. ..   ..   .. 
 .
 . . . . .  .  
   . 
 Gk1 · · · Gkk Gkg ··· GkN   Vk   Ik 
  Vg  =  (6)
    
 Gg1 · · · Ggk Ggg ··· GgN Ig 
    
 . .. .. .. .. ..  .   .. 
 .. . . . . .   ..   . 
GN 1 ··· GN k GN g ··· GN N VN 0
The network voltage solution is obtained by iterating the current injection values (2)-(3) with voltage solution from the
network (6), until the scheduled train power values are satisfied within tolerance. The iterative procedure is known as the
current injection power flow; it consists of the following steps:
Step 1: (Initialization) Form the [G] matrix on the left hand side of (6) and factorize it using optimally-ordered sparse LU
decomposition. Initialize all catenary node voltages to Vs , all track node voltages to zero, and the iteration counter
 m = 0.
Step 2: (Current Update) Increment the iteration counter m by 1 and compute the node current injection vector I (m) :

(m) PT k
Ik = −Ig(m) = (m−1) (m−1)
(7)
Vk − Vg
 
Step 3: (Voltage Update) Solve for V (m) in (6) using forward/backward substitution.
Step 4: (Convergence Test) Check if the power mismatch values (8), evaluated using the most recently available voltage and
current, are within tolerance for all trains:
(m) (m)
(Vk − Vg(m) )Ik − PT k ≤ ε (8)

If the convergence test is satisfied, stop and print the solution, otherwise go to Step 2.
A flowchart of the CI power flow is given in Fig. 4.

III. R EGENERATIVE T RAINS AND N ONREVERSIBLE S UBSTATIONS

Voltage rise in DC traction networks is associated with the operation of regenerative braking systems and gets more
pronounced when the supply is non-receptive, i.e when the supply is incapable of absorbing the regenerated power and
presents a high resistance to the train [24]. In practice, the rise in voltage is brought under control via a chopper-controlled
rheostatic brake. This section presents a two-phase approach to simulate the local control of voltage in the CI method (Section
II). Phase-1 (subsection III-A) is a sensitivity method to efficiently simulate the adjustment of regenerative power that is fed
to the network, when interleaved with the CI power flow iterations. Phase-2 (subsection III-B) employs the phase-1 solution
to compute the resistance of each regenerating train; it then proceeds with the solution of the network having non-receptive
substations via the CI method but with the regenerating trains represented using a voltage-current source model.
 5

Fig. 1. DC traction power network.

Fig. 2. Current source train model.

A. Regenerative Trains (phase-1)

The trains operating in the regenerative braking mode increase the catenary voltage while supplying their braking power
to the network. If the catenary voltage of any regenerating train exceeds the prescribed limit, the regenerated power that is
returned to the network must be reduced to bring the voltage within its maximum limit via squeeze control [5], [13]; this
is achieved by means of a local controller that dissipates part of the regenerated power in braking resistors [24]. Cai et al.
model the behavior of the local controller in power flow through successive power reduction, i.e. the power supplied by the
regenerating train is gradually reduced until the right power balance is found and the train voltage is at its limit [5]. Although
this approach is effective in controlling voltage, it may require a large number of CI iterations to achieve power balance.
The alternative technique presented herein is a sensitivity matrix approach for simultaneously adjusting the power injected
into the network from all regenerating trains with over-voltage. Consider for illustration that the trains connected at catenary
nodes c1, . . . , j, k, . . . , cn are operating in regenerative braking mode and have higher than the prescribed voltage limit (Vjmax );

Vp Vn Vs
Ip In

Fig. 3. Super-node substation model.

 6

Fig. 4. Flowchart of the CI power flow.

the corresponding sensitivity matrix equation is given by:

[H] [∆PT ] = [∆V ] (9)

 ∂Vc1 ∂Vc1 ∂Vc1 ∂Vc1

∂PT c1 ··· ∂PT j ∂PT k · · · ∂P T cn
 .. .. .. .. .. .. 

 . . . . . .


∂Vj ∂Vj ∂Vj ∂Vj
··· · · · ∂PT cn
 
 ∂PT c1 ∂PT j ∂PT k

[H] =  ∂Vk ∂Vk ∂Vk ∂Vk
 (10)

 ∂PT c1 ··· ∂PT j ∂PT k · · · ∂P T cn


 .. .. .. .. .. .. 

 . . . . . .


∂Vcn ∂Vcn ∂Vcn ∂Vcn
∂PT c1 ··· ∂PT j ∂PT k · · · ∂PT cn
   
∆PT c1 ∆Vc1
 ..   .. 

 . 


 . 

 ∆PT j   ∆Vj 
[∆PT ] = 
 , [∆V ] = 
  (11)
 ∆PT k   ∆Vk 
 
 ..   .. 
 .   . 
∆PT cn ∆Vcn
The elements of the voltage correction vector [∆V ] are given by
(m)
∆Vj = Vjmax − Vj , j = c1, . . . , cn (12)
(m)
where Vj is the voltage solution from the CI algorithm. The vector [∆PT ] is obtained by solving the linear system in (9),
and is used to update the regenerative power that gets injected into the network according to:
(m)
PT j = PT j + ∆PT j , j = c1, . . . , cn (13)
(m)
where PT j denotes the train power injection used in the last CI solution; note that if PT j falls outside its limit values, then
it is clamped at the violated limit. In summary, when a stable voltage solution is obtained from the CI method, the catenary
voltages of the regenerative trains are first checked to see if they violate their maximum prescribed limits. If there are violated
voltages, the sensitivity matrix equation (9) is solved to update the train powers employed in the next iteration of the CI
algorithm, which is initialized using the last voltage solution. The iterations between the CI voltage update and the sensitivity
computation/current update are continued until a stable voltage solution is found that (i) satisfies the train power mismatch
(8) and (ii) does not violate any voltage limit on regenerating trains. The implementation of the sensitivity based CI method
requires the following change to steps 2-4 in Section II, starting from a CI solution:
 7

Fig. 5. Flowchart of the phase-1 solution.

Step 2: (Sensitivity Computation/Current Update) For regenerating trains that exceed the prescribed voltage limit, adjust the
train
 (m)injection
 powers using (9)-(13). Increment the iteration counter m by 1 and compute the node current injection vector
I using (7).  
Step 3: (Voltage Update) Solve for V (m) in (6) using forward/backward substitution.
Step 4: (Convergence Test) Check if all regenerative train voltages do not exceed their prescribed limits by εV , and if the
power mismatch values (8) are within tolerance ε for all trains. If the convergence test is satisfied, stop and print the solution,
otherwise go to Step 2.
A similar procedure is implemented to restrict the power consumed by motoring trains whenever their catenary voltage drops
below a prescribed minimum value. Fig. 5 shows a flowchart of the phase-1 solution.
1) Sensitivity Matrix: The sensitivity analysis requires the computation of the partial derivatives in (10). Towards this end,
the modified conductance matrix in (6) is inverted to express the nodal voltages in function of the injection currents via a
modified resistance matrix:
    
V1 R11 · · · R1k R1g · · · R1N Vs
 ..   .. .. .. .. .. ..   .. 
 .   .
   . . . . .  . 
 
 Vk   Rk1 · · · Rkk Rkg · · · RkN   Ik 
 Vg  =  Rg1 · · · Rgk Rgg · · · RgN   Ig  (14)
    
    
 .   . .. .. .. .. .. .
 ..   .. .   .. 
  
. . . .
VN RN 1 ··· RN k RN g ··· RN N 0
∂Vk
To compute a sensitivity element on the diagonal of [H], for instance ∂PT k , consider only the injection currents that are
related to PT k as being nonzero:
PT k −PT k
Ik = , Ig = (15)
Vk − Vg Vk − Vg
The voltages Vk and Vg at the catenary and track connection nodes become:
PT k
Vk = (Rk1 + · · · ) Vs + (Rkk − Rkg ) (16)
Vk − Vg
PT k
Vg = (Rg1 + · · · ) Vs + (Rgk − Rgg ) (17)
Vk − Vg
Taking the derivative of (16)-(17) with respect to PT k , and making use of the chain-rule, gives:
  " ∂Vk # " Rkk −Rkg #
A B ∂PT k Vk −Vg
∂Vg = Rgk −Rgg (18)
C D ∂P V −V
Tk k g

where
Rkk − Rkg Rkg − Rkk
A=1+ 2 PT k , B= 2 PT k (19)
(Vk − Vg ) (Vk − Vg )
 8

Rgk − Rgg Rgk − Rgg

C= 2 PT k , D =1− 2 PT k (20)
(Vk − Vg ) (Vk − Vg )
The solution to (18) gives the sensitivity coefficients:
∂Vk (Rkk − Rkg ) (Vk − Vg )
= 2 (21)
∂PT k (Vk − Vg ) + PT k (Rkk + Rgg − Rkg − Rgk )
∂Vg (Rgk − Rgg ) (Vk − Vg )
= 2 (22)
∂PT k (Vk − Vg ) + PT k (Rkk + Rgg − Rkg − Rgk )
∂Vj
To compute an off-diagonal sensitivity element in [H], for instance ∂PT k , consider:
PT k
Vj = (Rj1 + · · · ) Vs + (Rjk − Rjg ) (23)
Vk − Vg
The derivate of Vj with respect to PT k , is obtained from (23) using the chain-rule:
 
∂Vj Rjk − Rjg Rjk − Rjg ∂Vk ∂Vg
= − 2 PT k − (24)
∂PT k Vk − Vg (Vk − Vg ) ∂PT k ∂PT k
Substituting (21) and (22) in (24) gives:
∂Vj (Rjk − Rjg ) (Vk − Vg )
= 2 (25)
∂PT k (Vk − Vg ) + PT k (Rkk + Rgg − Rkg − Rgk )

B. Nonreversible Substations (phase-2)

Even if all the regenerative train power injections are reduced to prevent any voltages exceeding their limits, the regenerated
power returned to the network may still be greater than the total motoring power. This will cause power to flow back into the
supply through some or all of the substations; while this may be possible with inverting substations [18]–[20], it is not realistic
with nonreversible substations as the rectifiers become reverse biased. Nonreversible substations are commonly available in
practical DC traction systems, and this makes it necessary to account for them in the power flow solution.
One approach is to check the substation power after a voltage solution is obtained. If reverse power is found to be flowing
through a substation, the substation is turned off and a new stable voltage solution is sought via the CI method. A problem
arises when all the substations are turned off; in this case the power flow solution with all trains modeled as current sources
cannot be reliably computed, due to the absence of any slack power source. The regenerative train model in Fig. 6 is used to
circumvent this problem: whenever one or more substations are turned off, the regenerating trains are converted into a voltage
source in series with a resistance connected to node k, and a current source Ik (equal to the current supplied by the regenerative
train) connected to node g. The voltage source is set to the maximum allowed catenary voltage Vkmax and the resistance Rk
is chosen according to:
PT k V max − Vk
Ik = = k (26)
Vk − Vg Rk
(Vkmax − Vk ) (Vk − Vg )
⇒ Rk = (27)
PT k
The choice of Rk according to (27) is such that when none of the substations are turned off, then the CI voltage solution
with all trains modeled as current sources (Fig. 2) will also be a solution to the network with the regenerative trains represented
as in Fig. 6. This is evident because for the same voltage solution, (26) stipulates that the regenerative train injection currents
are identical when using the models from Fig. 2 and Fig. 6. Consequently when any of the substations having reverse power
flow is turned off, the regenerative trains will have to supply less current to the network. This results in less voltage drop over
Rk and therefore a higher catenary voltage; this voltage cannot however exceed the maximum prescribed value Vkmax , which
is chosen as the voltage source value in Fig. 6. A flowchart of the phase-2 solution is given in Fig. 7.

IV. N UMERICAL R ESULTS

The CI power flow method (Section II) and the two-phase approach (Section III) were programmed in Matlab. The programs
were run on an iMAC having a 2.7 GHz quad-core Intel Core i5 processor with 4 MB L3 cache and 8 GB of RAM. The
parameters employed in the stopping criteria are ε = 0.0001 W in the CI method and εV = 0.0001 V for checking the
maximum permissible catenary voltage violations.
 9

Fig. 6. Voltage-current source regenerative train model.

Fig. 7. Flowchart of the phase-2 solution.

A. Benchmark Test Network with 18 Trains

The benchmark test network consists of a double track 750 V DC traction system whose diagram and data are available in
[1], [2], [25]; the upper permissible catenary voltage is set to 787.5 V. The system has 3 substations and 18 trains, and an
upper track and a lower track each having a branch. The ground network is represented using the finite cell model with a cell
defined for each 100 m of the track length. This gives a DC power flow problem with 779 nodes. Nineteen simulation cases
were studied, with case #0 having all trains in motoring mode, case #1 having the first train in regenerative braking mode and
the remaining ones in motoring mode, case #2 having the first 2 trains in regenerative braking mode and the remaining ones
in motoring mode, etc. Table I shows for each case the CI iterations and the execution time using the proposed sensitivity
based CI method (phase-1). The proposed method was compared with the power reduction and balance method in [5], and the
computational performance results are given in columns 4 and 5 in Table I. It is evident that the proposed method and that in
[5] are identical up to case #7, as there is no requirement to reduce the regenerated power that is fed in the network due to
over-voltage. The situation changes in cases #8 to #18, where it is seen that [5] requires a large number of power reduction
steps with the increasing number of regenerating trains; on the other hand, the sensitivity method maintains roughly the same
number of CI iterations and computing time. The speed-up factor (SUF = time of the method in [5]/time of the sensitivity
method), which is given in column 6 of Table I, increases to around 7 as the number of regenerating trains increases. Notably,
the proposed approach maintains convergence within 26 CI iterations even for the extreme case #18 in which all trains are in
regenerative braking mode.
Tables II and III show representative comparative results of the train power and catenary voltages for case #9. The results
in column 2 of Tables II and III reveal that the sensitivity method reduced the power of the regenerating trains #1, #3, and
#9 from 150 kW to satisfy the permissible pantograph voltage limit. The computation of the substations’ power, as given in
Table IV, shows that the substation at node #25 is conducting power to the AC network side. Given that such a situation is
not permissible for nonreversible substations, the computation proceeds to phase-2 in which the substation is blocked and the
regenerating trains are switched to the voltage-current source model in Fig. 6; the resistance in Fig. 6 is held constant at the
value given by (27). Column 4 in Table II gives the train power from the CI method with the regenerating trains switched to
their voltage-current source model (phase-2), and Table V shows that the substation at node #25 is now no longer conducting
power. Columns 3 and 5 in Tables II and III present the power flow results using the method in [5]; a power reduction step of
1 kW is employed to get results comparable in accuracy to those obtained from the sensitivity method. The results show that
the proposed method gives results almost identical to [5], with the minor differences attributed to using a fixed power reduction
step. Note that none of the catenary voltages in Table III exceed the set limit of 787.5 V, and the voltages corresponding to
 10

Table I
C OMPARATIVE P ERFORMANCE R ESULTS FOR THE N ETWORK WITH 18 T RAINS

Proposed PF Ref. [5]

Case # SUF
iter. time [ms] iter. time [ms]
0 43 148 43 148 1.0
1 38 147 38 147 1.0
2 32 148 32 148 1.0
3 25 150 25 145 1.0
4 17 145 17 145 1.0
5 16 145 16 145 1.0
6 16 145 16 145 1.0
7 15 146 15 146 1.0
8 26 163 52 165 1.0
9 28 164 1275 244 1.5
10 24 165 1777 291 1.8
11 26 166 2632 382 2.3
12 26 166 4239 528 3.2
13 28 166 5956 687 4.1
14 27 166 6275 729 4.4
15 27 166 7436 840 5.1
16 27 167 8816 960 5.8
17 25 166 9744 1049 6.3
18 26 166 11299 1196 7.2

Table II
T RAIN P OWER FOR CASE # 9 ( IN K W - ACTIVE SIGN CONVENTION )

Reversible SS Nonreversible SS
Train #
Proposed Ref. [5] Proposed Ref. [5]
1 75.75 76 47.84 48.07
2 150 150 95.39 95.56
3 126.18 126 105.49 105.38
4 150 150 147.40 147.40
5 150 150 147.80 147.80
6 150 150 149.25 149.25
7 150 150 148.32 148.32
8 150 150 149.16 149.16
9 71.76 71 45.32 44.91
10 -150 -150 -150 -150
11 -150 -150 -150 -150
12 -150 -150 -150 -150
13 -150 -150 -150 -150
14 -150 -150 -150 -150
15 -150 -150 -150 -150
16 -150 -150 -150 -150
17 -150 -150 -150 -150
18 -150 -150 -150 -150

motoring trains are substantially lower but still above the 650 V limit stipulated in [13]; while it is possible to reduce the
power drawn by motoring trains to raise the lower voltage values in a scheme quite similar to squeeze control [13], networks
such as the Swedish railway locomotive system will still remain in operation for voltage drops of around 40% [6].

B. Larger Test Networks

The proposed method was also contrasted with [5] on larger test instances with 36, 72, and 144 trains. These test instances
were formed starting from the typical parameter values and structure of the traction network in [1], [2]. Table VI shows the
number of trains, the number of substations, and the total number of nodes (including ground circuit modeling) in the three
networks; the electronic companion [25] includes schematics of the network interconnections and the complete data sets. The
traction systems were simulated under three scenarios involving 50%, around 75%, and 100% of the trains in regenerative
braking mode. The required CI iterations and computing time are presented in Table VII. It is noticed that the computational
performance of the competing power reduction and balance method [5] deteriorates at a fast pace, with both the increasing
number of trains and their proportion in regenerative braking mode; this is contrary to the performance of the proposed method
where the computing time remains within 0.5 s even for the network with 144 regenerating trains. The speed-up factor (SUF
in the last column of Table VII) shows that the computational advantage of the proposed technique relative to [5] increases
with system size and with the number of regenerating trains to reach around 98.
 11

Table III
C ATENARY VOLTAGES FOR CASE # 9 ( IN V )

Reversible SS Nonreversible SS
Train #
Proposed Ref. [5] Proposed Ref. [5]
1 787.50 787.40 787.50 787.44
2 773.95 773.91 778.94 778.90
3 787.50 787.44 787.50 787.45
4 746.95 746.93 747.68 747.66
5 742.86 742.85 743.55 743.53
6 746.02 746.02 746.23 746.23
7 774.79 774.78 774.93 774.92
8 758.70 758.70 758.86 758.85
9 787.50 787.32 787.50 787.39
10 687.94 687.91 696.97 696.93
11 676.93 676.91 679.74 679.71
12 668.64 668.62 669.93 669.91
13 664.23 664.22 665.27 665.25
14 718.90 718.89 719.21 719.21
15 702.50 702.50 702.79 702.78
16 659.86 659.85 662.83 662.81
17 697.54 697.54 698.83 698.82
18 713.04 713.03 713.31 713.30

Table IV
S UBSTATION VOLTAGE , C URRENT, AND O UTPUT P OWER FOR CASE # 9 (P HASE -1: R EVERSIBLE S UBSTATIONS )

Node # Positive [V] Negative [V] Current [A] Power [kW]

25 750.55 0.55 -176.93 -132.70
26 741.62 -8.38 400.07 300.05
27 744.27 -5.73 229.40 172.05

Table V
S UBSTATION VOLTAGE , C URRENT, AND O UTPUT P OWER FOR CASE # 9 (P HASE -2: N ONREVERSIBLE S UBSTATIONS )

Node # Positive [V] Negative [V] Current [A] Power [kW]

25 766.40 -4.27 0 0
26 741.91 -8.09 398.17 298.63
27 744.61 -5.39 224.25 168.18

Table VI
M AIN C HARACTERISTICS OF THE L ARGER T EST N ETWORKS

Trains (NT ) Substations (NS ) Nodes (N )

36 6 1558
72 12 3116
144 24 6304

Table VII
C OMPARATIVE P ERFORMANCE R ESULTS FOR THE L ARGER T EST N ETWORKS

Reg. Proposed PF Ref. [5]

SUF
Trains Trains iter. time [ms] iter. time [ms]
36 18 32 170 2420 466 2.7
36 28 29 174 13733 2065 11.8
36 36 26 177 22216 3301 18.7
72 36 32 188 4090 1244 6.6
72 56 29 207 27070 7529 36.5
72 72 26 214 44228 12288 57.4
144 72 32 346 8104 4001 11.6
144 112 29 410 53557 25679 62.7
144 144 26 430 87718 42116 98.0
 12

V. C ONCLUSION
This paper presented a two-phase approach for simulating local controllers in DC traction power flow applications that model
limited network receptivity. The approach builds on the commonly adopted CI power flow technique. Phase-1 is a sensitivity
method for adjusting the regenerated power that is fed back to the network in blended regenerative/rheostatic braking; the
power adjustment is computed to keep the catenary voltages within their permissible limits, but assumes that the stations are
reversible. Phase-2 accounts for limited network receptivity due to nonreversible substations; this is achieved in the CI method
by switching each of the regenerating trains to a voltage-current source model whose resistance value is computed from the
voltage solution in phase-1 without iteration, in contrast to [5]. Numerical results show that the sensitivity based CI method
obviates the need for many of the power reduction steps in [5], and gives a speed-up factor that significantly increases with
the network size and the number of regenerating trains.

R EFERENCES
[1] R. J. Hill and Y. Cai, “Simulation of rail voltage and earth current in a PC-based traction power simulator,” WIT Trans. on Engineering Sciences, vol. 3,
pp. 1–8, 1993.
[2] ——, “An efficient computational technique for a DC rail traction power flow simulator incorporating track branches,” WIT Trans. on Modelling and
Simulation, vol. 5, pp. 1–12, 1993.
[3] S. N. Talukdar and R. L. Koo, “The analysis of electrified ground transportation networks,” IEEE Trans. Power Appar. Syst., vol. 96, no. 1, pp. 240–247,
Jan. 1977.
[4] T. Kulworawanichpong, “Multi-train modeling and simulation integrated with traction power supply solver using simplified Newton-Raphson method,”
Journal of Modern Transportation, vol. 23, no. 4, pp. 241–251, 2015.
[5] Y. Cai, M. R. Irving, and S. H. Case, “Iterative techniques for the solution of complex dc-rail-traction systems including regenerative braking,” IEE
Proc. - Gener. Transm. Distrib., vol. 142, no. 5, pp. 445–452, Sep. 1995.
[6] L. Abrahamsson and L. Söder, “Fast estimation of relations between aggregated train power system data and traffic performance,” IEEE Trans. Vehicular
Technology, vol. 60, no. 1, pp. 16–29, Jan. 2011.
[7] L. Abrahamsson, T. Kjellqvist, and S. Östlund, “High-voltage DC-feeder solution for electric railways,” IET Power Electron., vol. 5, no. 9, pp. 1776–1784,
Nov. 2012.
[8] R. J. Hill, Y. Cai, S. H. Case, and M. R. Irving, “Discussion of: Iterative techniques for the solution of complex DC-rail-traction systems including
regenerative braking,” IEE Proc. - Gener. Transm. Distrib., vol. 143, no. 6, pp. 613–615, Nov. 1996.
[9] J. G. Yu and C. J. Goodman, “Modelling of rail potential rise and leakage current in DC rail transit systems,” in IEE Colloquium on Stray Current
Effects of DC Railways and Tramways, Oct. 1990, pp. 2/2/1–2/2/6.
[10] Y. Cai, M. R. Irving, and S. H. Case, “Modelling and numerical solution of multibranched DC rail traction power systems,” IEE Proc. - Electr. Power
Appl., vol. 142, no. 5, pp. 323–328, Sep. 1995.
[11] C. L. Pires, S. I. Nabeta, and J. R. Cardoso, “ICCG method applied to solve DC traction load flow including earthing models,” IET Electr. Power Appl.,
vol. 1, no. 2, pp. 193–198, Mar. 2007.
[12] M. Coto, P. Arboleya, and C. Gonzalez-Moran, “Optimization approach to unified AC/DC power flow applied to traction systems with catenary voltage
constraints,” Int. J. Elec. Power Energy Syst., vol. 53, pp. 434 – 441, 2013.
[13] P. Arboleya, B. Mohamed, C. González-Morán, and I. El-Sayed, “BFS algorithm for voltage-constrained meshed DC traction networks with nonsmooth
voltage-dependent loads and generators,” IEEE Trans. Power Syst., vol. 31, no. 2, pp. 1526–1536, Mar. 2016.
[14] C. L. Pires, S. I. Nabeta, and J. R. Cardoso, “DC traction load flow including AC distribution network,” IET Electr. Power Appl., vol. 3, no. 4, pp.
289–297, Jul. 2009.
[15] P. Arboleya, G. Diaz, and M. Coto, “Unified AC/DC power flow for traction systems: A new concept,” IEEE Trans. Vehicular Technology, vol. 61, no. 6,
pp. 2421–2430, Jul. 2012.
[16] P. Arboleya, M. Coto, C. González-Morán, and R. Arregui, “On board accumulator model for power flow studies in DC traction networks,” Elec. Power
Syst. Res., vol. 116, pp. 266 – 275, Nov. 2014.
[17] Z. Jia, Z. Yang, F. Lin, and X. Fang, “Dynamic simulation of the DC traction power system considering energy storage devices,” in IEEE Conference
and Expo on Transportation Electrification (ITEC Asia-Pacific), Aug. 2014, pp. 1–6.
[18] T. Suzuki, “DC power-supply system with inverting substations for traction systems using regenerative brakes,” IEE Proc. B - Electr. Power Appl., vol.
129, no. 1, pp. 18–26, Jan. 1982.
[19] B. Mellitt, Z. S. Mouneimne, and C. J. Goodman, “Simulation study of DC transit systems with inverting substations,” IEE Proc. B - Electr. Power
Appl., vol. 131, no. 2, pp. 38–50, Mar. 1984.
[20] D. Cornic, “Efficient recovery of braking energy through a reversible DC substation,” in Electrical Systems for Aircraft, Railway and Ship Propulsion,
Oct. 2010, pp. 1–9.
[21] I. Džafić, H.-T. Neisius, M. Gilles, S. Henselmeyer, and V. Landerberger, “Three-phase power flow in distribution networks using Fortescue
transformation,” IEEE Trans. Power Syst., vol. 28, no. 2, pp. 1027–1034, May 2013.
[22] I. Džafić, B. C. Pal, M. Gilles, and S. Henselmeyer, “Generalized π Fortescue equivalent admittance matrix approach to power flow solution,” IEEE
Trans. Power Syst., vol. 29, no. 1, pp. 193–202, Jan. 2014.
[23] J. D. Irwin and R. M. Nelms, Basic Engineering Circuit Analysis - 9th Ed. Hoboken, NJ: Wiley, 2008.
[24] R. G. Fletcher, “Regenerative equipment for railway rolling stock,” Power Engineering Journal, vol. 5, no. 3, pp. 105–114, May 1991.
[25] “DC railway traction network data sets,” https://www.dropbox.com/s/b9e69tfg0z8e65b/DC_Railway_Net.zip?dl=0, accessed: 2017-2-15.

Rabih Jabr (M’02-SM’09-F’16) was born in Lebanon. He received the B.E. degree in electrical engineering (with high distinction) from the American
University of Beirut, Beirut, Lebanon, in 1997 and the Ph.D. degree in electrical engineering from Imperial College London, London, U.K., in 2000.
Currently, he is a Professor in the Department of Electrical and Computer Engineering at the American University of Beirut. His research interests are in
mathematical optimization techniques and power system analysis and computing.
 13

Izudin Džafić (M’05-SM’13) received his Ph.D. degree from University of Zagreb, Croatia in 2002. He is currently an Associate Professor in the Department
of Electrical Engineering at the International University of Sarajevo, Bosnia. From 2002 to 2014, he was with Siemens AG, Nuremberg, Germany, where
he held the position of the Head of the Department and Chief Product Owner (CPO) for Distribution Network Analysis (DNA) R&D. His research interests
include power system modeling, development and application of fast computing to power systems simulations. Dr. Džafić is a member of the IEEE Power
and Energy Society and the IEEE Computer Society.

View publication stats