Non-linear Switched Model for Accurate Voltage Estimation and Power Flow Analysis of DC Railway Systems

Simulation tools are crucial to efficiently design the infrastructures and operations of DC electrical railway systems, including potential innovative technologies such as reversible traction power substations and energy storage systems. For this purpose, it is essential to accurately estimate the evolution of the voltage and power flows within the DC traction network, with fast computation time. This paper therefore proposes a new simulation approach for fast and accurate voltage estimation and power flow analysis of DC railway systems. It is based on the use of non-linear switched models for traction power substations and trains. The Modified Nodal Analysis is extended to consider such models, including the voltage drop control of the different subsystems, avoiding the necessity to use complex numerical iterative solvers. This new approach is validated and compared to an existing dynamical model and a conventional static model. The comparisons prove the relevance of the new approach, which provides validated and accurate results (less than 2 % error compared to the validated dynamical model) with fast computation time (speed up of 500 compared to the dynamical model). It can therefore be used to study, design, size, and optimize DC traction systems with new technologies aimed at saving braking energy.


Introduction
The use of simulation techniques is essential to design future infrastructures and plan the operations of efficient transport systems. Given the advantages of high transport capacity, punctuality, short isolation distance and environmental aspects, electric DC railway transportation is a preferred solution in urban area. However, global warming and the increase in the cost of energy lead to the need to further improve the overall efficiency [1]. Indeed, cost and energy optimizations are more and more expected in the transport sector to reduce its ecological footprint and make its economical aspect more resilient [2]. To this end, several approaches are investigated to improve overall efficiency: improving the design of trains (reduction in mass, aerodynamics, etc.), improving operations (optimization of timetables, eco-driving, etc.) [3], [4] and improving the infrastructure. Different solutions are therefore envisaged such as reversible Traction Power Substations (TPS) [5]- [8], Energy Storage Systems (ESS) [2], [8]- [16], or renewable energies [16]. Optimization techniques are also often used to improve the operation of specific subsystems (i.e. control of reversible TPS or ESS) [5], [8], [10] or of the entire railway system. Indeed, the choice of a specific solution, its design, its placement and sizing, and the evaluation of its economic interest are complex tasks. For example [2] optimizes the overall economic impact induced by the peak power demands of DC railway systems by using ESS. [7] proposes a new DC voltage control to optimize voltage regulation and power sharing using reversible TPS. [16] optimizes the operation of electric railways with renewable energies and ESS by considering a large-scale nonlinear optimization problem.
Several approaches and simulation tools have therefore been proposed to estimate the potential energy savings of DC railway systems. The results obtained with these tools can be used as inputs variables for an optimization problem to improve the design and operation of new railway systems, including potential innovative solutions [17], [18]. To this end, it is important to maintain good accuracy on the simulation results and a low computational load despite the increase in the complexity of the studied railway systems.
Conventional TPS, which connect the AC Distribution Network (ACDN) to the DC Traction Network (DCTN), and the over-voltage protection of trains play an important role in the solving procedure of the simulation tools. Indeed, they represent non-linear behaviours, which require complex models that are generally solved using numerical iterative techniques [19], [20]. Most existing and commercial tools consider the train as a static current source and use the conductance matrix iterative approach known as the Current Injection (CI) method [12]- [15], [18]- [20]. These currents are then injected into a Modified Nodal Analysis (MNA) method to solve and analyse the DCTN in a specific configuration according to the injected currents [21], [22]. Algebraic equations with iterative methods, such as Newton-Raphson, Point-Jacobi, or Zollenkopf's bifactorisation and incomplete Cholesky conjugate gradient methodologies are then necessary to make the problem converge according to the non-linear behaviours of the different subsystems [19], [20]. The objective of these iterative solvers is therefore mainly to accommodate the over-voltage protections of regenerative trains supplied by a DCTN with a limited receptivity [17], [19]. In this case, part of the braking energy must be dissipated in the rheostatic brakes [12]. The receptivity is limited by the voltage increase due to the resistive part of the DCTN, but also by the non-reversible TPS which act as equivalent capacitors when they are blocked [23]. However, most simulation tools simplify this second phenomenon and Auto-generated PDF by ReView IET Electrical Systems in Transportation IET-DCRailwayModels2. 1.docxM ainDocument only consider the increase in voltage due to the resistive part of the DCTN [19]. This problem is generally solved by considering reversible TPS or ESS [5]- [16], which contribute to improving receptivity. Other works aiming at estimating the energy consumption neglect this phenomenon because it occurs when all the TPS are blocked, which does not lead to a significant error on the global energy consumption [24]. This widely accepted simplification leads to simulation errors on power flows and voltage estimations when dealing with DC railway systems with low receptivity and regenerative trains [19]. This can have an impact on the design and sizing of potential innovative solutions, which are particularly interesting in these specific scenarios (i.e. low receptivity, low possibility of energy recovery, etc.). The worst scenario is when all TPS are blocked at the same time, leading to totally disconnect the DCTN from the ACDN [24].
The design, sizing, and optimization of railway systems with innovative solutions require dedicated models capable to accurately estimate the voltages and the power flows within the DCTN, with a fast computation time. Such models must be able to deal with all possible studied cases. To this end, a dynamical model has been proposed to highlight the physical interactions between the different subsystems, even when all TPS are blocked [23], [24]. This model was developed using the Energetic Macroscopic Representation (EMR) [25], a graphical description tool which highlights the energy properties of complex systems while respecting physical causality [26]. This model has proven its accuracy and has been experimentally validated [23]. It requires a small time-step to compute the dynamics of the system, which leads to an important computation time. This model, despite its good accuracy, cannot therefore be used for optimization purposes. A first comparison between this dynamical model and the static model used in [12], which uses the conventional CI method, was carried out. It has been shown that the static model is not able to deal with the specific scenario where all the TPS are disconnected. Voltage and energy errors have been estimated [24]. This paper proposes a new simulation approach for DC railway systems. It aims to obtain a model with a similar accuracy in terms of power flow and voltage estimations than that of the dynamical model, but with faster computation time. The novelty consists in the use of EMR to systematically respect the physical causality. Switched models are therefore introduced to consider the non-linear behaviours of the different subsystems. The conventional MNA is extended to consider such models. Furthermore, respecting physical causality and the limitations of subsystems directly in switched models avoids the need to use complex numerical iterative solvers. The dynamical, static, and switched models are compared and evaluated according to their accuracy.
In section 2, the general description of studied systems is given. Section 3 describes the models of the basic subsystems of DC railway systems. The existing dynamical model and the conventional static model are explained in section 4. Then, section 5 develops the new simulation approach. Finally, simulation results, comparisons, analysis, and discussions are provided in sections 6 and 7.

Studied DC Railway Systems
In this paper, lowercases are used for time-dependent variables, while uppercases are used for variables without time dependence, and underlined variables are for vectors or matrices, while non-underlined variables are for scalars.

Studied conventional DC Railway Systems
A typical DC railway system is shown in Fig. 1. It consists in 3 basic subsystems: TPS, Trains (T), and DCTN. The ACDN supplies the different TPS, which convert and transfer energy from the ACDN to the DCTN. The most typical DC voltages are 750 V, 1500 V and 3000 V. Conventional TPS consist in transformers and diode rectifiers. The energy transfer is unidirectional, from the ACDN to the DCTN. The DCTN allows energy to be distributed to the different trains. It consists in a catenary or a third rail in both directions of the railway track. Parallel electrical connections between the two directions are often used to reduce the overall impedance of the DCTN (see Fig. 1) [20].
The electrical connection between the train and the DCTN is made through the pantograph or the third rail (Fig. 2). The equipment is supplied through one or several input filter(s) composed of smoothing inductor(s) and DC bus(es). The auxiliaries (Aux.) (lights, compressors, air conditioning, etc.), rheostatic brake, and traction subsystems are all connected to the DC bus(es). The traction subsystem consists in inverters, machines, gearboxes, and wheels. The performances of the train in traction phases can be affected by the current or torque limitations as well as by the "overcurrent protection" (Fig. 3.a) [12], [25]. In such cases, the traction torque is linearly decreased (coefficient koc in Fig. 3.a) when the DC bus voltage uf is lower than a defined limit (high traffic can induce important voltage drops and high currents). Consequently, the possible acceleration is limited, which may induce a delay in the scheduled timetable. When koc = 1 (uf > Uoc-max), the full torque is available; when koc = 0 ibk Aux.

Filter
Brake Aux.   [11]. But it is not always possible, and it depends on the traffic conditions. In this case, the rheostatic brake must dissipate all or part of the braking energy when the DC voltage reaches a defined limit. It is activated by the "over-voltage protection" (Fig. 3.b) [15], [18]. This does not cause any delay, but this reduces the overall efficiency. Indeed, when this protection is activated (coefficient kov Fig. 3.b), the energy recovery is limited or cancelled. The current ibk of the rheostatic brake is expressed by (1), with itot the total current (auxiliary plus traction sub-systems). The relationship between the voltage thresholds and the average no-load DC voltage Ess-0 of the TPS is given by (2).
in traction mode in braking mode

General Problematics and Innovations
Conventional DC railway systems have proven their robustness and effectiveness. Nevertheless, in recent years, new improvements have been proposed to increase the efficiency. Additionally, new solutions enable smooth power peaks, stabilized voltage, or off-grid traction possibility.
The optimization of the timetable and drive cycles allows the reduction of the energy consumption without the need to modify the existing infrastructure [3], [4]. But this is more effective in cases of high traffic. Furthermore, this is very sensitive to the small disturbances (ex: delays) that can occur during operations. In addition, this does not allow voltage stabilization or off-grid traction possibility. Another solution is to install inverters in parallel with the diode rectifier of the TPS to achieve reversibility while keeping the robustness of the diode rectifier ( Fig. 4.a) [5], [8]. This solution improves the efficiency of the entire system by recovering the braking energy of the trains on the ACDN. Nevertheless, it does not allow the reduction of the power peaks, the stabilization of the DC voltage, or the off-grid traction possibility. The voltage drop control of the inverter is presented on Fig. 4.a, with pinv the actual power absorbed by the inverter, Pinv-max its maximum power capability, and uinv the actual DC voltage of the catenary [12]. In this case, the inverter is only used to allow power flow from the DCTN to the ACDN. Another possibility is to install wayside ESS ( Fig. 4.b). In this case the braking energy can be stored and reused later during the traction phases to support the DCTN. This improves overall efficiency as well as voltage regulation and reduces power peaks at the TPS. The corresponding voltage drop control is shown on Fig. 4.b, with pess the power absorbed or injected by the ESS, Pess-max its maximum power capability, and uess the actual DC voltage at the connection with the DCTN [11], [12]. Other possibility is to consider onboard ESS. In such a case, on-board ESS can offer off-grid traction possibilities. More details about such technologies and other examples are in [5]- [16].
Assessment or optimization studies are required before implementing such solutions. The cost and energy profitability need to be estimated, especially to choose the most appropriate solution, and design and size the selected technology. Simulation tools are widely used to carry out such studies. However, existing tools need to be improved to be more accurate and to consider new components and controls. For instance, conventional railway systems were designed based on maximum and minimum values of the main variables (i.e. max. power at TPS, min. voltage of the DCTN, etc.) without the need to know precisely the power flows and voltage variations in the DCTN. However, new systems with innovative technologies require more advanced studies with more accurate models and simulation tools. Indeed, as previously presented, the controls for innovative solutions as well as for braking systems are mainly based on the DC voltage level [11], [12], [25]. These new simulation tools must accurately predict the evolution of voltages and power flows within the DCTN, with a fast computation time, and with sufficiently simple parameters for the models [24].

Models of the Basic Subsystems
The models of the different subsystems (TPS, trains, and DCTN) of DC railway systems will be the common and basic models used by the different simulation approaches (see Sections 4 and 5). The Energetic Macroscopic Representation (EMR) formalism is used as a common and unified representation tool to organize models [26]. EMR is a graphical description that highlights the energy properties of complex systems. It organizes the system into interconnected basic elements: source (green oval), accumulation (orange crossed rectangle), mono-physical conversions (orange square) and distribution (orange double square) of energy. All elements are connected according to the action and reaction principle. The product of the action variables with the reaction variables leads to the power exchanged [26]. In addition, all the components are described respecting the physical causality (i.e. the integral causality) [27].

Model of the Conventional TPS
By neglecting the harmonics, the equivalent electrical circuit of the TPS is shown in Fig. 5.a. Its EMR is composed of different elements that convert energy from the ACDN to the DCTN (Fig. 5.b). The ACDN is a source of energy which imposes its AC voltage vg-3φ on the transformer. A coupling element represents the transformer including its iron losses

Model of the Train
The train model takes into account its different components (see § 2.1). It determines the current itot, absorbed or injected, on the DC bus as a function of the power paux of the auxiliaries, the power ptr of the traction subsystem, and the voltage uf of the DC bus (3). The power paux is assumed to be known and is simulated by a time-dependent power profile. The power ptr is estimated as a function of the line profile (slope, wind, etc.), the mechanical dynamics of the train (mass), the resistive efforts (friction, aerodynamic, etc.), the components losses, the electro-mechanical limitations (current, torque, voltage, acceleration, etc.) and the speed of the train. The speed can be predetermined based on the scheduled drive profiles or can be generated by an adaptive velocity generator to adapt the electrokinematical behaviour of the train according to its different limitations and protections. More details are available in [24] and [25].
The equivalent electrical circuit of the train consists in a current source itot (3) (aux. + traction subsystems) connected in parallel with the rheostatic brake on the DC bus Cf (Fig. 5). The brake is considered as a current source that imposes the current ibk (1). Finally, the resistance Rf represents the losses in the input filter [23]- [25]. The EMR of the train will be described in the following sections because it depends on the implementation of the DC bus and braking subsystem models.

Model of the DCTN
The model of the catenary (or third rail) and return rail (or ground) segments of the DCTN assumes linear resistance distribution along the line. The resistances rc-ik and rr-ik are respectively the resistances of the catenary and return rail of the segment between two subsystems i and k. An example is given in Fig. 6.a, where TPS1, T1, T2, TPS2 and the parallel electrical connection between the two directions are respectively the subsystems numbers 1, 2, 3, 4 and 5. The parallel connection is considered as a subsystem, but without any source. These resistances are expressed by (4), with Δx-ik the distance between both subsystems (length of the segment), and Rl-c and Rl-r the linear resistances (in Ω/m) of the catenary and return rail, respectively. These resistances depend on the positions of the subsystems. Parallel electrical connections between both directions are often used to reduce the overall impedance of the DCTN and can be taken into account in the model. Nevertheless, simplifications leading to an equivalent catenary system for the two directions are often considered ( Fig. 6.b) [5]- [20]. This simplified model assumes that all the subsystems are supplied by a common catenary, thus reducing the number of nodes to consider. In this case, the equivalent resistances of each segment are expressed by (5)    where Kp-c is a weighting coefficient representing the reduction of the equivalent resistance as a function of the number of parallel connections. The overall DCTN impedance can be divided by 2 with regular parallel connections (Kp-c = 0.5), while it is not divided without connection (Kp-c = 1).

Existing Simulation Approaches
Two existing approaches for the simulation of DC railway systems are presented. They will be compared to the new simulation approach in Section 7. Both existing approaches use the TPS and DCTN models previously presented. The simplified DCTN model is considered instead of the complete model for simplicity reason. The first existing approach considers a dynamical model [23]- [25] while the second uses a static model [12]. Both approaches use the conventional MNA [21], [22], which is briefly explained in the following subsection. Both approaches are then presented.

Conventional Modified Nodal Analysis method
The MNA objective is to solve complex DC networks composed of several subsystems, which can be current or voltage sources. The conventional MNA is expressed by (6), where G is the conductance matrix (with geq-ik = req-ik -1 ), R is the matrix that considers the input resistances of each voltage source, and B is a matrix corresponding to the Kirchhoff's currents equations and is composed of ±1 for the voltage source nodes whose branch relations are introduced (it is zero elsewhere) [21], [22]. The vectors J and E include the known variables of the subsystems, respectively for the current sources and voltage sources. The vectors V and I correspond to the unknown variables, respectively the voltages at each node (connection of the DCTN with the subsystems) and the currents of each voltage source. These unknown variables are obtained by inverting the matrix relation (6).
The sizes of the matrices and vectors depend on the number of nodes nbi and voltage sources nbm. These sizes are given in Table 1 for the simplified DCTN model (see Fig. 6.b). In this case, nbi is equal to the number of subsystems (4 in Fig. 6.b) and V only represents the catenary voltage because the return rail is assumed at to be 0 V everywhere (ground connection). However, the MNA can be adapted to distinguish the catenary and return rail voltages (each subsystem thus induces 2 nodes instead of 1; 1 on the catenary and 1 on the return rail) and to consider the complete DCTN model (each parallel connection is considered as a subsystem inducing 2 supplementary nodes). In such a case, the sizes of the matrices/vectors can thus drastically increase.
Considering the simplified equivalent DCTN model, G is determined by (7), where i and j are the node numbers corresponding respectively to the row and column in G. For a given node i (source node), which is connected to one or more other nodes k(i) (destination nodes); Gi,j (with j = i) is the negative value of the sum of all conductances geq-ik(i) connected to the node i; Gi,j (with j = k(i) as destination node) is the positive value of geq-ik(i); it is zero elsewhere. B is determined by (8) if node is a current source if node is a voltage source Finally, V corresponds to the voltage value of the catenary at each node i, and I corresponds to the current of each voltage source m. The formulation of the MNA is therefore strongly dependent on the models used for each subsystem (voltage or current source). Two different simulation approaches are hereafter presented.

Dynamical Model-based Simulation Approach
The first simulation approach is based on a dynamical model (DM) of the train DC bus (see Fig. 5) [23]- [25]. The respect of the integral causality leads to represent the train (from the point of view of the DCTN) as a voltage source uf with a series resistance Rf. The DC bus voltage uf is expressed by (12), where uf-init is the initial voltage and Cf is the value of the capacitor. When a train (ex: T1) is connected to a node i on the DCTN and corresponds to a voltage source m, then it imposes on the MNA the known input variables Ji = 0 (not a current source) and Em = uf-1 (voltage source). The unknown variables are Vi = ut-1 (catenary voltage) and Ii = it-1 (train current). The other elements are Bi,m = 1 and Rm,m = Rf-1. With the DM, a train is always considered as a voltage source. According to the TPS model, the ON-state model (ex: TPS2) is a voltage source. When it is connected to another node i and corresponds to another voltage source m, it imposes on the MNA the known input variables Ji = 0 (not a current source) and Em = ess0-2 (voltage source). The other elements are Bi,m = 1 and Rm,m = Rss-2. The output variables are Vi = ess-2 (catenary voltage) and Ii = -iss-2 (TPS current). However, the TPS is a current source in OFF-state. It is connected to the node i but it does not correspond to a voltage source m. It thus imposes on the MNA the known input variable Ji = -iss-2 = 0 (current source equal to zero). The unknown variable is Vi = ess-2 (catenary voltage). The corresponding element in B is equal to zero.
The complete EMR of the railway system composed of 2 TPS and 2 trains (see Fig. 6.b) is given in Fig. 7. The train (ex: T1) is represented by a current source itot-1, a coupling element representing the parallel connection on the DC bus with the rheostatic brake (13), and an accumulation element representing the DC bus, which imposes the voltage uf-1 on the DCTN. The brake current ibk-1 is imposed by the braking management strategy (1). The interactions with the DCTN are the voltages uf of the DC buses and the currents it of the trains. These variables are merged into 2 vectors to be considered as inputs to the MNA, as explained above. The output variables of the TPS (yss) are also merged to provide inputs to the MNA, which returns the variables zss to the TPS models with respect to the interaction principle. The DCTN is therefore represented by a central conversion element, which corresponds to the MNA, and 2 coupling elements that merge the TPS and trains variables, respectively. The MNA also provides variables, such as the catenary voltages ut, which are not highlighted in the EMR. Finally, all TPS are coupled to the ACDN to get the overall energy consumption.
Algorithm 1 shows the global solving procedure with this simulation approach. The initial variables are the TPS states, the DC buses voltages, the powers of auxiliaries and traction subsystems, and the positions of the subsystems; step 1 consists in calculating the currents if (13) absorbed on the DC buses as a function of the DC buses voltages (see Fig. 3), the powers (3), and the braking energy management strategy (1); step 2 identifies the current/voltage sources according to the TPS states (the initial states are those from the previous time-step; the trains are voltage sources) (see Fig.  5.c); step 3 formulates and solves the MNA (6); step 4 updates the TPS states according to the MNA results (see Fig.  5.c); step 5 saves the new states; in step 6 if a least one state is different from its initial state, a new calculation is required (GOTO step 2), otherwise go to step 7; step 7 updates the voltages of the DC buses (12) for the next time-step as a function of the currents of the trains it and the currents of the DC buses if, which were respectively determined from the MNA results (vector I) (step 3) and the train models (step 1); Finally step 8 ends the procedure. The next time-step can start.
The respect of the physical causality and interactions (imposed by the EMR) leads a physical description which is perfectly able to predict the behaviours of DC railway systems. This model has been validated using experimental measurements on a real subway line [23]. It can therefore be used as reference in this paper.

Static Model-based Simulation Approach
The second existing approach is derived from the conventional approach used in [12]. The model has been adapted for use with the EMR methodology, but it provides results strictly identical to those obtained with the CI method used in [12]. The train uses a static model (SM), which neglects the DC bus capacitor. Consequently, the physical interactions are simplified in certain cases, in particular in the vg-3φ ess-1  (14). The corresponding element in B is zero. With this approach, the DC buses of the trains are not represented. The differences with Fig. 7 are that the trains directly impose the current it = if on the DCTN, which then returns to the trains the voltages uf (14). The rest of the model is identical to the dynamical model-based approach (DM).

= − . (14)
The main drawback of this approach, contrary to the dynamical model-based approach, is the inability to solve the system when the DCTN is totally disconnected from the ACDN, which means that all TPS are in the OFF-state [24]. This can occur frequently when the net power of the line is null (low receptivity line with regenerative trains). In such a case, there is no voltage reference to solve the DCTN. This situation has no impact on the overall energy consumption. However, this approach assume that no energy is exchanged in the DCTN, that the DC voltage is equal to the rated value, and that the current is equal to zero for each subsystem.
Algorithm 2 shows the solving procedure. The initial variables and steps 1, 3, 4, 6 and 8 are the same as with Algorithm 1; step 2 identifies the current and voltage sources based on the TPS states (see Fig. 5.c) (the initial states are those from the previous time-step; the trains are current sources); in step 5 if all TPS are in the OFF-state, thus the system cannot be solved (all the voltages are considered equal to the rated voltage and all the currents are zero) and Algorithm 1 directly stops and goes to step 8, but if at least one TPS is in the ON-state, thus the procedure continues; step 7 calculates the DC buses voltages (14); Finally step 8 ends the procedure. The next time-step can start.

Discussion
A first comparisons between both existing approaches have been recently performed to highlight the advantages and drawbacks of each approach [24]. The main advantage of the dynamical model-based approach (DM) is its ability to simulate the entire system and to provide accurate simulation results in all the studied cases. Conversely, the static modelbased approach (SM) is not able to provide the correct solution when all the TPS are blocked. More generally the DM allows better estimations of power flows and voltages within the DCTN, which is important for studying innovative supply structures as already discussed in this paper.
However, regarding the computation time the SM is clearly faster than the DM. This point is a great advantage for carrying out intensive simulations in order to optimize and design the DC railway system and its subsystems and components. Indeed, the DM needs to compute the high dynamics of the DC buses, which requires a small time-step and lead to important computation times. This is not necessary with the SM.

New Simulation Approach
The objective of the new simulation approach is to have similar accuracy than the DM but with faster computation time. It uses non-linear switched models (SWM) and does not need to compute the DC buses capacitors.

Principle of the New Simulation Approach
The principle of the non-linear SWM is to consider a current/voltage switched model for the TPS (as in the previous sections) but also for the trains. With this approach, the DC bus capacitor of the train is not taken into account, but instead a specific SWM is proposed. Based on the overvoltage protection characteristics, the behaviour of the rheostatic brake can be divided into two states. Firstly, it is a current source ibk = 0 when it is not activated. This is the case when the train is in traction phase (itot ≥ 0) OR when the train is in braking phase without activation of the over-voltage protection (uf < Uov-min) (Fig. 8.a). Secondly, the brake is a voltage source uf with a drop control when the train brakes (itot < 0) AND the over-voltage protection is activated (uf ≥ Uov-min) (Fig. 8.b). In this case, uf is expressed by (15) as a function of the threshold voltage Uov-max, the over-voltage coefficient kov (see Fig. 3.b), and Δov (16). Furthermore, based on (1) and (13) and with if = it, kov can be expressed as a function of the currents it and itot (17). Finally, the voltage uf is determined by (18) after injecting (17) in (15). Based on this approach, the characteristic of the brake (current ibk versus voltage uf) is presented on Fig. 8.c. Two equivalent models are thus used for the train (including the brake): a current source itot when the brake is in the OFF-state (with it = if = itot) and a voltage source uf with a drop control (over-voltage protection) when the brake is in the ON-state (18). Simple switching conditions are defined (see Fig. 8.c).

Extended Modified Nodal Analysis Method
The idea behind this non-linear SWM approach is to extend the conventional MNA. The equivalent voltage source uf (with the drop control) can be directly included in the MNA when the brake is activated, thus avoiding the need to use specific and complex iterative solvers. This means that the over-voltage protection (see Fig. 3.b), or any other subsystems with voltage drop controls, can be directly included in the extended MNA. The extended MNA is expressed by (19) and is solved by inverting the relationship.  The sizes of the new matrices and vectors are defined in Table 2, with nbi the number of nodes, nbm the number of voltage sources, and nbw the number of voltage sources with drop controls. First, Z is a zeros matrix. Then, F is determined by (20), with w the voltage source with the drop control and m the voltage source (with or without drop control). Fm,w is equal to 1 if the voltage source m is a source with drop control w, it is zero elsewhere. Contrary to the conventional MNA, E is now only constructed by the constant voltage sources without drop control. Em is therefore zero if the voltage source m uses a drop control (21). Finally, O is an identity matrix and D is a vector containing the slopes of the different drop controls. In the case of a train (ex: T1), where the brake is activated (ON-state), the element Dw is given by (22). The vector W corresponds to the threshold voltages of the voltage drop controls (i.e. Uov-max (18) in this paper). The vector U corresponds to the unknow variables, which are the voltages uf of the DC buses of each train where the overvoltage protection is activated. This extended MNA solves complex DC networks composed of several subsystems, which can be current or voltage sources, with or without voltage drop control. In addition, it is possible to switch between a system with one or more voltage source(s) with drop control(s), and a system without any voltage drop control. In such a case, the additional matrices and vectors must be added dynamically according to the studied case (i.e. all the new elements are empty if there is no drop control).

Non-linear Switched Model
The new simulation approach is thus based on a nonlinear SWM for the TPS and the trains. The TPS model is the same as in Section 4. However, a SWM is also implemented for the train depending on the state of its rheostatic brake (see Fig. 8). The extended MNA is dynamically adapted to take into account the over-voltage protection. A train (ex: T1), which is connected to a node i and where the brake is in the OFF-state, is therefore a current source and imposes the current itot as input of the MNA with Ji = itot-1 = it-1. The outputs of the MNA are Vi = ut-1 (catenary voltage) and uf-1 (14). The corresponding elements in B are zeros. However, when the brake is in the ON-state, the train is considered as a voltage source m with a drop control w. In this case it imposes as inputs of the extended MNA the known variables Ji = 0 (not a current source), Em = 0 and Ww = Uov-max_1 (voltage source with drop control). The other elements are Bi,m = 1, Rm,m = Rf-1, Fm,w = 1, Ow,w = 1, and Dw (22). The outputs of the extended MNA are Vi = ut-1 (catenary voltage), Ii = it-1 (train current), and Uw = uf-1 (DC bus voltage). Contrary to the two previous approaches, a train can now be modelled as a current source or as a voltage source with a drop control depending on the behaviour of its rheostatic brake (see Fig. 8).
Algorithm 3 shows the solving procedure. The TPS and trains states, DC buses voltages, auxiliary and traction powers, and trains positions are the initial variables; step 1 defines the currents itot (3); step 2 identifies the current-  voltage sources based on the TPS (see Fig. 5.c) and trains (see Fig. 8.c) states; step 3 formulates and solves the extended MNA (19); step 4 updates the TPS and trains states according to the MNA results from step 3; step 5 saves the new states; in step 6 if a least one state is different from its initial state a new calculation is required (GOTO step 2), if not continue to step 7; step 7 calculates the actual DC buses voltages of the trains in current sources (14) (these voltages are directly provided in U for the trains with the brake in the ON-state). The procedure stops at step 8. The next time-step can start.
With this approach, the rheostatic brake and overvoltage protection are merged into the extended MNA (Fig. 9). The switched model of the rheostatic brake is therefore not represented on the EMR because it is intrinsically taken into account during the formulation of the extended MNA. Certain variables are not directly available but can be obtained based on intermediate results of the extended MNA, such as the train currents it (vector I), the catenary voltages ut (vector V) and the brake current ibk (23).

Discussion
The new proposed approach does not simulate the high dynamics of the DC buses, which should allow to increase the simulation time-step and thus substantially reduce the computation time compared to the DM, even if the size of the extended MNA is more important than the size of the conventional MNA due to the added elements. In addition, the new approach must lead to high accuracy in the estimations of power flows and voltages due to the implementation of the real over-voltage protection of the train in the non-linear switched model. This new SWM will be compared with the two existing approaches in terms of simulation results accuracy in the next sections.

Simulation Scenarios and Analysis
The test line is derived from an existing conventional subway line described in [25]. But only the first 6.9 km of the line are taken into account in this paper (8 passenger stations and 3 TPS), for simplicity reason. The main data of the line is provided in Table 3 [25]. The velocity and power profile (traction plus auxiliaries) have been defined to respect the kinematical limitations (torque, acceleration, etc.) [25].
Different scenarios, with different headways (from 1 to 20 min), are defined to compare the different approaches (Fig. 10). According to the headway, the line requires between 1 and 18 trains (Fig. 10.a), which induces different distances travelled (Fig. 10.b) and capacities (in kp.km/h for a thousand persons transported over 1 km in 1 h) (Fig. 10.c). This also leads to blocking all the TPS at the same time more or less frequently (in % of the simulated period) (Fig. 10.d).
First simulation results are obtained with the validated DM (Fig. 11). These results will be used as references for the comparisons in the following section. The higher the traffic is, the lower the average voltage of the DCTN is ( Fig. 11.a). This is due to the increase in the absorbed power, which induces greater voltage drops in the catenary. Note that at low traffic, this average voltage may be higher than the no-load voltage of the TPS. This is explained because the DC voltage increases when the trains brake, especially when the DCTN  Fig. 9. EMR of DC railway system with the non-linear SWM-based simulation approach  is disconnected from the ACDN (blocked periods). The increase in the traffic also improves the reuse of braking energy by other trains in traction phases ( Fig. 11.b). This reuse of the braking energy is 0 % when a single train is running (headway of 20 min) and can reach 30 % with high traffic. That means that up 30 % of the energy absorbed by a train can be recovered during the braking phases. Since each train transports an average of 150 passengers (light rail vehicle), the energy absorbed on the DCTN (consumed by trains only) to transport 1 passenger over 1 km is presented in Fig. 11.c. Since reuse of braking energy improves with increasing traffic, energy consumption per passenger and per km is lower with high traffic than with low traffic. Finally, Fig. 11.d shows the total energy consumption absorbed on the ACDN, including the DCTN and TPS losses. In this case, it is observable that the optimal energy consumption is for an headway of about 3 min: lower traffic results in lower reuse of braking energy, but higher traffic induce more losses in the system, which are not compensated by the better energy reuse.

Comparisons of the Different Approaches
The comparison of the overall energy consumption is not performed in this paper but can be found in [24]. The conclusion of this paper was that the differences are not in terms of overall energy consumption, but in terms of local power flows (energy absorbed, and energy recovered) and voltage estimations within the DCTN. The SM and SWM are therefore compared to the reference results obtained with the validated DM (see Section 6) on these aspects.

Comparisons of the Simulation Results
First, the comparisons between the different approaches are achieved with the same simulation time-step of 1 ms, which is necessary to compute the DC bus dynamics with the DM approach. Comparisons are made only on the accuracy of the results. Indeed, the 3 approaches are adapted to respect the physical causality imposed by the EMR and do not require any iterative numerical solvers (no convergence problem). As the DM has already been experimentally validated, it is thus used as the reference. The study quantifies the errors (in %) on the energy absorbed and the energy recovered on the DCTN by all trains, for the different scenarios (Fig. 12). The black bars correspond to the SM, while the red bars correspond to the SWM. In addition, the errors are determined over the entire simulation period ("all") but also over the "blocked" periods only (all TPS blocked). The errors on the entire simulation period with the SM can reach up to 6.6 % (6 min headway) and 30.9 % (9 min headway), respectively on the energy absorbed ( Fig. 12.a) and the energy recovered (Fig. 12.c). With the SWM, the maximum errors are only 0.04 % on the energy absorbed (9 min headway) and 0.24 % on the energy recovered (7 min headway). By focusing on the "blocked" periods ( Fig. 12.b and Fig. 12.d), these errors are 100 % for the energy absorbed and 100 % for the energy recovered with the SM (not able to deal with this simulation case), while they only reach 0.45 % (7 min headway) and 0.7 % (9 min headway), respectively, with the SWM. These errors depend on the scenarios. They can be important for certain scenarios with the SM, while the results are accurate in all cases with the SWM.
A similar comparison on the estimated average DCTN voltages, which supply the different trains, is presented on Fig. 13. The differences (in Volts) between the SM (black bars) and the reference (DM) reach up to -25 V (20 min headway) and -200 V (20 min headway) on average, respectively for the entire simulation period (Fig. 13.a) and for the blocked periods only (Fig. 13.b). For the same periods, the differences between the SWM (red bars) and the reference reach up to ±0.4 V (5 min headway) and ±3.4 V (2 min headway), respectively. The SWM is therefore also more accurate than the SM for estimating the voltages.
The catenary voltage ut (Fig. 14.b) and the current it (Fig. 14.c) of the train 1 (T1) are displayed to better highlight the differences between the 3 approaches. The same power (auxiliary plus traction) (Fig. 14.a) (ON-state) to the ACDN. The accuracy of the conventional SM there fore strongly depends on the ratio between the total simulation period and the blocked periods (see Fig. 10.d) as previously presented.

Analysis and Discussion
The comparisons prove the relevance of the SWM, which provides equivalent results to those of the validated DM, even during the blocked periods. The SWM is capable of locally and accurately estimating catenary voltages as well as power flows within the DCTN. As previously discussed, this is particularly important for studying, designing, and sizing new railway systems, including new technologies aimed at saving braking energy. In the studied cases, an error of up to 30.9 % on the recovered energy can be made with the SM. In addition, unlike the DM, the SWM is not limited by the simulation time-step required to compute the DC bus. It is thus possible to use larger time-steps with the SWM. The computational load can therefore be considerably reduced, which is crucial for carrying out optimization studies.
The study presents the impacts of the simulation timestep on the computation time and results (with the 5 min headway scenario). The DM requires a time-step of 1 ms, which leads to a computation time of 612 s for a simulation period of 1 h. This time-step cannot be increased for the DM because it is limited by the dynamics of the DC bus. Table 4 summarizes the speed-up allowed by the SWM compared to the DM. For example, the computation time is divided by 515 with the SWM (with a time-step of 1 s) compared to the DM (with a time-step of 1 ms). Table 4 also highlights the impact of the time-step on the accuracy of the voltage estimation over the entire simulation period ("All.") and over the "blocked" periods only. The results present the differences between the average voltages estimated with the SWM (with different time-steps) and the DM (with a time-step of 1 ms). Table 5 summarizes the impact of the time-step on the error of estimation of the power flow (i.e. energy absorbed, and energy recovered) for the wh ole simulation period ("All.") and for the "blocked" periods only. These two tables show that larger the time-ste p of the SWM is, the more the errors (compared to the DM) increase. However, these errors are reasonable, and the results obtained with the SWM remain accurate by using a time-step to 500 ms (even up to 1 s). Then, the accuracy declines exponentially for larger time-steps. The SWM can therefore allow an acceleration of the computation time by around 500 (with a time-step of 1 s) compared to the validated DM, while keeping a good accuracy.

Conclusion
A new simulation approach is proposed for DC railway systems. It is based on the use of switched models for a better analysis of power flows and a better estimation of DC voltages. The switched models represent the non-linear behaviours of the diode-based Traction Power Substations (TPS) and the rheostatic brakes of the trains. The Modified-Nodal Analysis is extended to consider such models and include the voltage drop controls of the different subsystems. The comparisons are made with a complex dynamical model, which has been experimentally validated, and an equivalent static model, which uses conventional simulation approach. The new approach gives results similar to those of the dynamical model, but with a faster computation time (speedup of about 500). Furthermore, it allows better accuracy compared to the conventional static model and can simulate DC railway systems in all cases, even in the specific case where all TPS are blocked at the same time. Indeed, the error on the recovered energy can reach up to 30 % with the static model in the worst scenario, while it is less than 2 % error with the proposed switched model. More comparisons with other recent power flow solvers will be performed in a future work in terms of accuracy and computation time.
The proposed simulation approach is then used to analyse power flows and accurately estimate the voltage of a DC traction network under different scenarios, with fast computing capabilities. This is crucial for studying new technologies such as energy storage systems or reversible TPS, whose controls and managements are based on the DC voltage levels. This approach can be considered for