Modeling, Simulationand Analysis of On-Board Hybrid Energy Storage Systems for Railway Applications

In this paper, a decoupled model of a train including an on-board hybrid accumulation system is presented to be used in DC traction networks. The train and the accumulation system behavior are modeled separately, and the results are then combined in order to study the effect of the whole system on the traction electrical network. The model is designed specifically to be used with power flow solvers for planning purposes. The validation has been carried out comparing the results with other methods previously developed and also with experimental measurements. A detailed description of the power flow solver is beyond the scope of this work, but it must be remarked that the model must by used with a solver able to cope with the non-linear and non-smooth characteristics of the model. In this specific case, a modified current injection-based power flow solver has been used. The solver is able to incorporate also non-reversible substations, which are the most common devices used currently for feeding DC systems. The effect of the on-board accumulation systems on the network efficiency will be analyzed using different real scenarios.


Introduction
The use of power flow algorithms for planning traction networks is a widely-accepted technique [1,2]. However, the use of accurate models of the network and the trains may result in very complicated simulations [3]. In [4], the authors proposed a methodology to perform a fast estimation of the aggregated railway power system and traffic performance. In [5], for instance, the authors tried to reduce the computational burden of the whole simulation, proposing a compression technique that reduces the number of necessary simulations for evaluating the performance of the traction network. In [6], the energetic macroscopic representation (EMR) was used in order to simplify the development of the train and network mathematical models. This approach has been use extensively. For instance, in [7], the authors compared different traction substation models. The models presented using the EMR were extremely accurate and precise [8], and they may describe every single part of the system and are appropriate for applications in which the transient behavior of the trains is critical, such as hardware in the loop applications [9]. A very detailed model of an electric train fed by supercapacitors connected to fuel cells was presented in [10]. Pseudo-static models have been traditionally less accurate, but much faster [11]. They are usually combined with power flow solvers that must be able to manage the observe that on-board and wayside energy storage are better solutions in terms of efficiency. On-board accumulation is slightly better in terms of efficiency because the losses in the DC system are slightly lower, but wayside energy storage can be more cost-effective since a single accumulation device can be shared by all the trains in the system. Similar conclusions were obtained in [26]; in this work, the use of wayside energy storage improved the overall efficiency by more than 20% and also improved the controllability and the voltage profile of the network. Similar conclusions were obtained in [24,27] supporting the previous conclusions. In [24], the use of energy storage at the substation level reduced the annual cost of energy by 13%. We would like to point out that the main contribution of this paper is the mathematical model; a detailed comparison between the infrastructure alternatives considering wayside and on-board energy storage and reversible substations is beyond the scope of the work, which is focused specifically on providing an accurate and efficient model to simulate on-board energy'storage.
What this work presents is a decoupled model that considers the train and the storage system as separate devices located at the same electrical point and interacting with each other and with the network. The authors did not consider the input filter, and the losses in this device can be estimated and added to the electromechanical and electrochemical device. This approach reduces the complexity of the mathematical model without reducing the accuracy when compared with previous proposed models [28]. The model has been embedded in a commercial software package that uses a modified current injection method to solve the power flow problem [12]; the solver allows combining regenerative trains with diode-based non-reversible substations. It must be remarked that even when there is the possibility of combining regenerative trains with Voltage Source Converter (VSC)-based reversible substations [29,30], in the vast majority of the cases, the existing infrastructures and also the new ones use non-controlled and non-reversible substations due to their robustness, reliability, and low cost [31]. The combination of the two technologies must be done carefully. It should be noticed that without a very precise coordination, the power regenerated by one train could be used to feed another train, but if there is no power demand when the train is regenerating, the voltage in the network will increase, and the substations will be blocked. In such cases, usually the voltage limit is reached, and most of the regenerated power must be burned in the rheostatic braking system. The use of on-board accumulation systems could alleviate this situation [32]. However, it will be demonstrated, the size and effectiveness of the storage system depend highly on the trains' power profile and the schedule.
The results obtained using this kind of simulations are very useful also to study how the accumulation systems are cycled in real applications. There are many on-going research projects that try to predict the life of the accumulation systems depending on the cycling. For instance, the common thought is that load leveling the charge and discharge of the battery packs can contribute to extending the life of the batteries. Recent studies revealed unexpected results. In [33], the authors cycled LiNiCoAl and LiFePO 4 modules with 18,650 cells 750 times in a period of six months, and they concluded that the degradation was higher with constant current loads than using dynamic pulse profiles. There is still much research to do in this direction, and every new storage technology requires these kinds of studies. The output of the pseudo-static simulation could be a valuable input for the researchers that try to reproduce real dynamic cycles in accumulation systems.
In the next section, the mathematical model of the train plus the energy storage device in different working modes is explained. Section 3 presents the results obtained using the proposed model in a real traction network. This section has four subsections: in the first three subsections, we describe the feeding infrastructure, the rolling stock, and the proposed scenarios. In the last subsection a deep analysis of the results is presented. Conclusions are stated in Section 4.

Mathematical Model of the Train and the Storage System
The proposed train model including the storage system has been simplified in order to lighten the computational burden. In [28], the storage system and the traction equipment of the train were considered as a single element; in this case, the storage system and the traction system separation will substantially simplify the solving procedure. This simplification in the model will not reduce the accuracy because both models are equivalent. The time savings in simulation were around 20%. The solving procedure in this case allows the decoupling between the storage system and the traction equipment, and this is the main reason for this simplification in the mathematical model. From now on, we will use the term train when referring to the traction and breaking device, and storage will be reserved for the accumulation device. Both devices will be modeled as separate, but coupled models.
The storage control will consider the train behavior, as well as the catenary voltage, so we will assume that the train and the storage will be traveling along the grid separated, but connected to the same catenary point.
The joint power and current can be easily calculated as the sum of the current of the train and the storage. For the sign criteria, we will consider that the train will have positive current and power in traction mode when it is absorbing power from the electric grid. In the other case, the storage will have positive power and current when it is in discharging mode. In Figure 1a,b, the schematic representation of the train and storage mathematical models is depicted in traction and braking modes, respectively. We only consider two situations that are the most common ones, the train in traction mode with the storage system discharging and the train in braking mode with the storage system in charging mode.

Train Plus Storage Working in Traction Mode
The behavior of the train and the storage in traction mode is summarized in Figure 1a, and the functions f t1 , g d1 , and g d2 are depicted in the curves presented in Figure 2a-c. In the next subsections, a detailed description of the train and the storage model will be presented.  The train model will read first the power reference provided by the external software; we will refer to this power as the mechanical reference power (P train mech,re f ). It is possible to use coupled models that generate the mechanical power reference during the electrical simulation [34]. However, in most of the cases, the train mechanical model and the electrical problem are solved separately. To obtain the train electrical reference power (P train re f ), the mechanical power will be divided by the efficiency of the electromechanical conversion that will include the motors, converters, and the rest of the equipment. This efficiency will be labeled as train efficiency in traction mode (ν train t ). This operation can be executed before launching the solver with all positive powers provided by the external software package, so the input file for the solver will already consider this efficiency according to the next expression: It must be pointed out that this electric power reference (P train re f ) is not the final electric power that the train will demand from the grid (P train ). The overcurrent protection can limit the power requested from the grid when the catenary voltage (V cat ) is very low for protection purposes.
Each train will have two configuration parameters to set the overcurrent protection (V train reg,t and V train min ); in Figure 2a, the function f t1 is depicted, and the mathematical expression of such function is expressed in (2). As can be observed, the actual power consumed by the train depends on the catenary voltage and the train electric power reference, as well as the parameters that define the overcurrent protection (P train = f t1 (P train re f , V train reg,t , V train min , V cat )). The maximum amount of power that can be extracted from the primary active power source of the storage system (ultracapacitors, batteries, flywheel, or whatever other technology used) is labeled as P STO max and depends on the level of charge of the storage system (E STO ) with the function g d1 (see Figure 2b) that is expressed in (3).
The maximum discharging power will be zero if the energy stored is zero, and it will increase in a linear way until a specific energy level is achieved (E STO reg,d ). The maximum discharging power will be constant for the interval between E STO reg,d and the maximum energy that the system can store (E STO max ). This maximum power (P STO max ) can be calculated at each instant once the storage energy level is updated, but before launching the power flow algorithm, because it is independent of the catenary voltage and the train behavior.
The maximum power available from the storage primary source (P STO max ) has to be multiplied by the electrochemical conversion process efficiency in discharging mode (ν STO d ) to obtain the maximum available electric power to be injected into the system (P STO max,av ). This operation can be also executed before launching the power flow solver.
Due to the storage system overcurrent protection, the maximum available power could be constrained if the voltage is too low by means of the function g d2 (see Figure 2c). The output of this function is the actual available power that the storage device can inject into the system considering already the catenary voltage constraint (P STO av ). To define the function g d2 , we need to define also the regulation parameters (V STO reg and V STO min ). The analytical expression of the function can be found in (5). Each train can have different regulation parameters. In the case of the storage system, these two parameters will be common for charging and discharging mode.

Coupling between the Train, Storage, and Network
Regarding the interaction between the train, the storage device, and the grid, many operational philosophies can be considered. In this particular case, one of the most popular working modes in real applications has been adopted, giving priority to the energy extraction from the storage system. When the train demands a given amount of power, if this power is available in the accumulation system, it is going to be extracted from it. If not, the train will extract the available power from the storage system, and the rest will be imported from the network. With this control philosophy, the actual power that the storage system is going to inject into the system is going to be calculated following the next expression: Due to the employed power flow solving procedure [11], we will use the previous iteration of the catenary voltage to calculate first the power of the train (P train ), then we will calculate the storage power (P STO ). The catenary net power, representing the whole set (train plus storage device), can be calculated as P cat = P train − P STO . Once the catenary power is calculated, a new iteration of the power flow will be launched to obtain a new value of the catenary voltage. The simulation will stop once the difference between all catenary voltages in all nodes in two successive iterations is lower than a specific tolerance.
Once the convergence is achieved, the energy level at the storage system must be updated using the following expression:

Train Plus Storage Working in Braking Mode
An analogous procedure can be used for describing the train and the storage device working in braking mode. Such behavior is summarized in Figure 1b and the functions f b1 , g c1 , and g c2 depicted in the curves represented in Figure 2d-f. In the next subsections, a detailed description of the train and the storage model will be presented.

Train Behavior (Positive Power in Traction Mode)
Again, the train model will read first the power reference provided by the external software (P train mech,re f ). To obtain the train electrical reference power (P train re f ), the mechanical power will be multiplied by the efficiency of the electromechanical conversion (see Equation (8)). This efficiency will be labeled as the train efficiency in braking mode (ν train b ).
The squeeze control can limit the power requested from the grid when the catenary voltage (V cat ) is very high for protection purposes. Each train will have two configuration parameters to set the squeeze control (V train reg,b and V train max ). In Figure 2d, the function f b1 is depicted. The mathematical expression of such a function is expressed in (9). As can be observed, the actual power that the train is going to regenerate depends on the catenary voltage and the train electric power reference, as well as the parameters that define the overcurrent protection

Storage Behavior (Negative Power in Braking Mode)
The maximum amount of power that can be injected into the primary active power source of the storage system (P STO max ) depends on the level of charge of the storage system (E STO ) with the function g c1 (see Figure 2e) that is expressed in (10). Again, this maximum power (P STO max ) can be calculated at each instant once the storage energy level is updated, but before launching the power flow algorithm, because it is independent of the catenary voltage and the train behavior.
The maximum power that could be injected into the storage referring to the catenary side of the storage system (P STO max,av ) has to be calculated using the electrochemical conversion process efficiency in charging mode (ν STO c ) according to (11). This operation can be also executed before launching the power flow solver.
The overcurrent protection of the storage system is considered also in braking mode using the function g c2 (see Figure 2f). The output of this function is the actual available power that the storage device can extract from the system considering already the catenary voltage constraint (P STO av ). To define the function g c2 , we need to define also the regulation parameters (V STO reg and V STO min ). The analytical expression of the function is similar to the one for the traction mode (see Equation (5)).

Coupling between the Train, Storage, and Network
To be consistent with the control philosophy previously described, the system will assign priority to the energy injection into the storage system. When the train regenerates a given amount of power, the storage system will try to use it (if possible) to increase its energy level. If not, the train will inject the maximum allowed power into the storage system, and the rest will be injected into the network. The actual power that the storage system is going to use for charging is going to be calculated following the next expression: The catenary net power, representing the whole set (train plus storage device) can be calculated as P cat = P train − P STO . At each instant, the energy level of the storage system can be calculated as follows:

Results' Analysis
In this section, we will analyze in depth the effect of the on-board accumulation system on the trains and network behavior in different scenarios. This section will have four subsections. In the first one, we will describe the basic feeding infrastructure, lines, and substations. In the second subsection, we will describe the rolling stock used in the simulations. The third subsection will be focused on the description of the different scenarios that will be analyzed in Section 4.

Feeding Infrastructure
The case study in this article will focus on the study of a real network consisting of two lines of 30.84 km and 36.93 km. The voltage level of the system was 3000 V. The simplified diagram of the network can be observed in Figure 3. Line blue is the longest; it had nine stops and four electrical nodes labeled as S1, S2, S3, and S4. The red line shared the first two electrical nodes with the blue line, and it had 17 stops and six electrical nodes labeled as S1, S2, S3, S4, S5, and S6. There were a total of eight electrical nodes and seven lines. Among the eight nodes, only three of them represented feeding substations; the rest were just nodes without any connection with the AC system. The three substations were placed in the nodes S3 and S5 of the red line and in the node S3 of the blue line. The three substations had the same characteristics. All of them were composed by a power transformer with rated power of 3 MW and a short circuit voltage of 5%; the no load output voltage of the rectifier was 3000 V, and the voltage at rated load (1000 A) was 2880 V. The equivalent impedance of each of the three substations in forward mode (AC to DC) was 270 mΩ. The equivalent impedance of the Energies 2019, 12, 2199 9 of 21 overhead conductor and the rails (return circuit) were respectively 28.605 mΩ/km and 7 mΩ/km. In Table 1, the lengths of the different segments of the red and blue lines can be observed.
Non-reversible substation

Rolling Stock
The train used in both lines was an electrical multiple unit (EMU). The whole unit was 2.940 m wide 4.265 m high and had a total length of 98.05 m, with an unladen weight of 157.3 t. The units were composed of five cars, the two ends having a driver's cabin and a normal floor. The middle car had a normal floor, while the other two cars had a low floor. The five cars were supported on two types of bogies, the trailer bogie and the tractor. The tractor bogie was always shared between two cars. The train was designed to use a standard Iberian track gauge (1668 mm) at a maximum speed of 120 km/h with almost 1000 passengers, although it could reach 160 km/h with minor modifications. The maximum total power of the train was 2.2 MW, and it had regenerative braking. In the base case, the trains were not going to be equipped with on-board energy storage systems. However, we considered the possibility of adding an on-board accumulator system based on a hybrid battery/ultracapacitor technology. The power profiles obtained during the simulation in the accumulation system would help the manufacturer to determine the percentage of energy that must be stored in the battery or ultracap parts. The electromechanical efficiency of the trains in traction and braking modes (ν train t , ν train b ), as well as the electrochemical efficiency of the storage system in charging and discharging modes (ν STO c , ν STO d ) were set to 0.95. The total accumulation system capacity (E STO max ) was 7 kWh, and the on-board energy storage device rated charging and discharging power (P STO rated,c , P STO rated,d ) was 1 MW. Regarding the protection curves of the trains and the storage elements, the minimum and the regulation voltage of the train in traction mode (V train min , V train reg,t ) were set to 1980 V and 2280 V. The same values have been selected for the minimum and regulation voltage of the energy storage system (V STO min , V STO reg ). In braking mode, the regulation voltage and the maximum voltage of the squeeze control (V train max , V train reg,b ) were 3300 V and 3600 V, respectively. In the cases in which the on-board accumulation system was activated, the system would be initialized with no charge.
In Table 2, the data summarizing the behavior of the trains in the different trips are collected. In the first column, we can observe the required mechanical power to complete the trip. It can be seen that the slope of the blue line was steeper because the difference between the power required for outward and return journeys was greater than on the red line. The average trip considering the two lines and both directions needed 229 kWh. The mechanical regeneration capacity in Column 2 is the available mechanical power that can be regenerated. Columns 3 and 4 contain the required electrical power and the electrical regeneration capacity considering already the efficiency of the electromechanical conversion. The electrical regeneration capacity was usually around 40% of the required electrical power, except in the S1 to S4 trip of the blue line. In this case, because the train ascended a steep slope, the regeneration capacity was much lower, around 22% of the required electrical power. In the fifth column, we can see the minimum electric consumption. This consumption was calculated as the required electrical power minus the electrical regeneration capacity. Off course, this is a theoretical consumption that considers that we took advantage of all electrical power available for regeneration. This is not true, mainly because of two reasons. First, part of the power that was available to be injected into the catenary was burned in the rheostatic braking system when the squeeze control was activated in order to maintain the catenary voltage below the maximum level. In addition, if the train was equipped with on-board accumulation, the efficiency of the electrochemical conversion during the charging and discharging process also reduced the percentage of available regenerated power that could be reused. For these reasons, we uses these minimum consumption figures as a theoretical ceiling to compare the different solutions, but we must be aware of the fact that we will not reach this theoretical ceiling.

Description of the Selected Scenarios
The authors developed eight different scenarios to study the influence of the accumulation system on the network, as can be observed in Table 3. There were four different paths for the trains, from S1 to S6 and from S6 to S1 for the red line and from S1 to S4 and from S4 to S1 for the blue line. Four different traffic densities were considered. Light traffic scenarios used a train headset of 50 min with 10 departures for each of the above-described routes. The medium traffic scenario considered a train headset of 35 min with 14 departures per route. The dense traffic scenario launched 24 trains per route with a headset of 20 min. Finally, the heavy traffic scenario launched 47 trains per route with a headset of 10 min. The simulation interval was very similar for all scenarios, and it went from eight hours and 18 min for the light traffic scenario to eight hours and 28 min for the heavy traffic scenario. Each of the four traffic densities were simulated without and with on-board accumulation systems without modifying the feeding infrastructure. The obtained results are presented in the next section and will clarify the effect of the on-board accumulation systems on the railway traction networks depending on the density of the traffic.

Analysis of the Results
The above-described scenarios were simulated, and the obtained results are analyzed in this subsection. With the software used for simulation, we can obtain time-varying series of each electrical variable in the train or in the feeding network. An example of this detailed analysis can be observed in Figure 4, in which Train Number 4 on the red line route starting from S1 is represented in Scenario L1. It must be noticed how on seven occasions, the train burned part of the regenerated energy using the rheostatic system due to the high catenary voltage, even when the train was equipped with an energy storage system and the storage capacity was not full. It can be observed also how the power extracted/injected into the catenary differed from the train mechanical reference due to the efficiency of the electromechanical conversion process, but also because part of the power was provided by the energy storage system. In Figure 5, the fourth train of the blue line starting from S1 is represented, in this case for the heavy traffic scenario (H1). It has to be remarked that the voltage level was much lower and the network receptivity higher. The burned power in this scenario was nearly zero since the overvoltage protection was only activated on four occasions with a very short duration.  The analysis of the time-varying curves was very interesting since it showed how electrical variables were correlated and it helped to understand the mathematical model proposed in the previous sections. However, in order to determine the effectiveness of the combination of the trains plus the infrastructure, another kind of analysis must be carried out. In Figure 6 is represented the power extracted from the AC network at substation S5 of the red line during the first four hours of simulation for all scenarios. Each subfigure represents two scenarios with the same train headset, but with and without the accumulation system. The correlation between the power extracted from the AC network was higher in the low traffic scenarios. Figure 7 depicts the energy extracted from the AC network at substation S5 in the eight possible scenarios. It must be noticed that the red solid line represents always the scenarios with energy storage. Obviously, scenarios with the same traffic were quite correlated. It should be pointed out that in the light, medium, and dense traffic scenarios, the energy consumption from the AC network in the cases with energy storage systems was a little bit lower. Paradoxically, that is not the case for the heavy traffic scenario in which the energy consumed from the AC network was higher when the trains were equipped with energy storage systems. As will be observed, in heavy traffic scenarios, it was more efficient to share the energy surplus with other trains than storing it in the on-board accumulation system.  In Figure 8, we represent the Marey diagrams of each scenario just for the red line. The horizontal axis represents time; in this case, we represented the first 80 min of simulation. In the vertical axis, we represent the position of the train. Solid black lines represent trains from S1 to S6, while dashed-dotted lines represent trains circulating from S6 to S1. Vertical red lines represent the instants at which all the substations were blocked at the same time. It must be noticed that in light, medium, and dense traffic scenarios, the percentage of instants in which all substations were blocked at the same time dropped drastically when we added on-board energy storage systems, improving this percentage by more than 10%. However, in the heavy traffic scenario, the percentage of blocking instants was already very low (5%) in the case without on-board accumulation. The installation of on-board accumulation in the heavy traffic scenario produced the blocking of all substations in only 0.8% of the cases. In the next paragraphs, we will analyze the aggregated results obtained for the eight scenarios from three points of view. First, we will present a full summary of the system. Second, a summary from the point of view of the trains will be analyzed, and finally, we will show a summary from the point of view of the substations. The use of aggregated data help us to measure the impact of the on-board storage solution in different traffic scenarios for the same feeding infrastructure. In Table 4, we summarize the most representative energies in the system; all the numbers are in MWh. In the first row, we have the electrical energy required by the trains. This energy is going to be the reference energy, since it only depends on the number of trains and does not depend on the equipment of the trains. This energy was the same with and without accumulation. For the light traffic scenario, the 40 trains required a total of 9.62 MWh. For the heavy scenario, the 188 trains required 45.2 MWh, the relation between the number of trains and the required electrical energy being linear. Something similar happened with the second and the third row, which represent the regeneration capacity and the minimum consumption, respectively. The regeneration capacity considers the electromechanical conversion efficiency, and the numbers obtained represent the electrical energy that could be used for injecting into the catenary or charging the storage system. The minimum consumption is a theoretical concept that cannot be reached since it is obtained by subtracting the regenerated capacity from the required electrical energy. Up to now, these numbers only depended on the traffic, and did not vary whether or not the trains had on-board energy storage equipment. The next row (fourth row) represents the energy demanded by the train at the catenary level. It differs from the electrical energy required by the trains for one reason in the case of trains without an accumulation system and two reasons in the case of trains with an on-board accumulation system. In the first case, the train demand was different from the electrical energy requirement because the over-current protection prevented the absorption of too much power in the case of low catenary voltage. In the second case, part of the electrical energy provided by the train can be provided by the on-board accumulation system. As can be observed, for the light traffic scenario without accumulation, the trains absorbed 0.72 MWh less than the required energy, so we would have delays in the system. For the heavy traffic scenario, the difference was 1.3 MW. However, in order to compare the scenarios, it is better to compare percentages with respect to the required electrical energy. These percentages can be found in Table 5. The authors are fully aware that this second table is redundant, but we think that it is important to analyze at the same time the data in MWh and in percentage with respect to the required electrical energy. The train demand in the scenario L0 was 92.6% of the required electrical energy, which means that 6.4% of the required energy cannot be provided. In order to distinguish the correlations between the different variables with the different traffic scenarios with and without energy storage and extract conclusions about the trends of the energy savings and consumption, the energetic data in percentage with respect to the total energy demanded by the trains are also represented in Figure 9. As can be observed, we had the same trends as the traffic increased for the cases with and without energy storage. However, even when the trend was the same, the values were very different, and it was necessary to add two different y-scales in order to compare the plots. The real energy demanded by the trains increased with the traffic, as well as the energy injected by the trains. That means that with denser traffic, the number of instants in which the overcurrent and overvoltage protections were activated was lower. We also reduced the energy provided by the substations and the grid losses when we increased the traffic.   Table 4. Black bars represent the cases without energy storage, and they are associated with the y-scale on the left. Red bars represent the cases with energy storage, and they are associated with the y-scale on the right.
The non-supplied energy in the case of the heavy traffic scenario without energy storage was only 0.6%. The non-supplied energy is represented in Row 9 of Tables 4 and 5. This non-supplied energy in percentage is a good index of the network congestion, as it can be observed that the more trains in the system, the less the network congestion. In general, we can state that the installation of on-board energy storage always reduced the amount of non-supplied energy. In the worst scenario (light traffic), the non-supplied energy was reduced from 6.4-1.7% when we added storage to the trains. Row 5 represents the energy regenerated by the trains that was actually injected into the feeding system. Obviously, this energy increased with the traffic, but it was always reduced when we added energy storage to the trains. In Row 6, the trains' net energy was obtained subtracting Row 5 (energy injected) from Row 4 (energy demanded). We can see that in the light, medium, and dense traffic scenario, the net energy was always lower when the trains were equipped with energy storage. However, the net energy was nearly the same in the scenarios H0 and H1. Row 7 represents the energy provided by the substations. It is very interesting to remark that even in the worst scenario (L0), the substations provided only 80% of the energy demanded by the trains, the rest being provided by other trains in braking mode. Again, for the light, medium, and dense traffic scenarios, the energy provided by the substations was reduced when we added the energy storage. In the heavy traffic scenario, this trend was inverted, and the energy provided by the substations was slightly higher when we added energy storage to the trains. This is because when we had many trains and they were very close to each other, it was more efficient to use the regenerated energy in other trains than the stored energy in the on-board accumulation system. This is coherent with the data obtained representing the losses in the feeding system (Row 10). In the heavy traffic scenario, the losses were higher when the trains were equipped with an energy storage system. Finally, Row 8 represents the energy burned in the rheostatic system of the train, in all cases, it suffered a significant reduction when we added the on-board accumulation system.
In Table 6, we can observe the trains' average behavior in the different trips and in all scenarios, as well as the average trip considering all possible routes in all scenarios. The table is split into four blocks. In the first one, we represent the total energy obtained from the catenary. It should be noticed that the slope of the blue line was steeper than the slope in the red line. That is the reason why the energy obtained from the catenary in the outward and return journeys was so different, as well as the energy injected into the catenary. It must be noted that for the average trip, the energy obtained from the catenary increased when we increased the traffic because the voltage was more stable (within the limits), and the non-supplied energy decreased. The minimum consumption (theoretical) for the average trip defined in Table 2 was 151 kWh, and this number was obtained subtracting the electrical regeneration capacity from the required electrical energy. The best scenario in terms of average net consumption was H0 (heavy traffic without accumulation), with a net consumption of 164 kWh, only 5.8% above the theoretical consumption. Adding accumulation to the heavy traffic scenario increased the average trip net consumption by 1 kWh. In the rest of the scenarios, adding accumulation always improved the average trip net consumption. For instance, in the light traffic scenario, the net consumption for the average trip was 183 kWh without accumulation and 167 kWh with accumulation. It is very important to remark that when we considered the net consumption in the average trip for comparing with the theoretical limit, we were not considering the non-supplied energy that was quite high in the light, medium, and dense traffic scenarios without accumulation. When we added accumulation to the system, the non-supplied energy in the average trip was always lower than 2% with respect to the required electrical energy for the average trip; and in the case of heavy traffic, lower than 1%. Table 7 contains the energy analysis of the substations for all scenarios; it is split into three blocks. The first block represents the energy injected into the DC traction system per each substation and the total. In the second block, we can observe the energy injected by each substation per train trip. Finally, this number is compared with the minimum theoretical consumption in Block 3. Again, we can observe that in the scenario H0, the total energy injected into all substations was only 9.3% above the minimum consumption.
Regarding the performance, both models (decoupled and integral) were equivalent from the point of view of the equations, and they provided the same result. However, the time savings on average when simulating the proposed decoupled model were around 20%. On average, the number of iterations of the decoupled model compared to the integral one was higher, from 6.5 iterations (average) to 7.9 iterations (average) per instant. The iterations using the decoupled model were faster. The average time to complete an iteration with the decoupled model was 0.39 ms, while the integral model spent 0.55 ms.

Conclusions
A simplified model of a storage system for power flow purposes was presented and tested. The proposed approach considered the train and the accumulation system as different devices placed at the same point of the railway traction network. Both devices were coupled and interacted with each other and with the electrical network. The proposed mathematical model reduced the computational burden of the simulation when it was compared with models that considered the train and the storage in an integral way. The model was embedded in a power flow solver, and extensive results were provided.
Analyzing the results obtained, it can be stated that the behavior of the system when adding this kind of technology cannot be easily predicted. It has been observed how reducing the train headset and increasing the traffic in the system, the non-supplied energy was reduced, as well as the number of instants at which all substations were blocked. In general, we could state that adding on-board accumulation systems to the trains reduced the burned energy in the rheostatic braking system due to the squeeze control activation and also the non-supplied power. However, this effect was highly correlated with the train headset. In order to perform comparisons with different traffic levels, all the energetic data must be normalized using as a rated value the total energy demanded by all trains in each specific scenario. In the cases studied in this work, we observed how the congestion in the network was inversely correlated with the traffic. That means that increasing the traffic, we alleviated the network, and the voltage profile was lower, so the activation of the overvoltage protection was less frequent. We observed also that with the heavy traffic scenario, the trains were able to absorb nearly 98% of the demanded energy when the on-board accumulation was installed and 77% without accumulation. The energy regenerated reached nearly 30% in the heavy traffic scenario with accumulation, when in the same situation without accumulation, it was only 9%, half of the injected regenerated power in the light traffic scenario with accumulation. As a general conclusion, it could be stated that up to a reasonable level, the increasing the traffic can alleviate the system, making it more efficient in relative terms. However, each specific infrastructure can have different behaviors, and it is important to make these kinds of studies during the infrastructure planning stage in order to make the correct decisions. Acknowledgments: The authors would like to thank CAFTurnkey & Engineering, especially Peru Bidaguren and Urtzi Armendariz for their support during the development of these models.

Conflicts of Interest:
The authors declare no conflict of interest. CAF TE provided the infrastructure data to perform the study. The Spanish Ministry of Science, Innovation and Universities had no role in the design of the study; in the collection, analyzes, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

P train
mech,re f Mechanical power provided by an external software package; it is positive in traction mode and negative in breaking mode. Storage system regulation energy in discharging and charging modes, respectively, used in the functions g d1 and g c1 . Storage system overcurrent protection parameters; see functions g c2 and g d2 .

Functions
f t1 , f b1 Functions for obtaining the actual power of the train in traction and braking mode, respectively; this is the real power absorbed from the grid. g d1 , g c1 Functions for obtaining the maximum amount of power that can be extracted from or injected into the storage system, respectively; this function relates the maximum power in the physical storage system depending on the energy level. g d2 , g c2 Functions for obtaining the available discharging and charging power in the storage system, respectively; this power considers the storage system's overall efficiency and the network conditions.

E STO
Energy in the storage system at a specific instant.

P train
Actual power that the train exchanges with the electrical traction network. P STO Actual power that the storage system exchanges with the electrical traction network.

P cat
Actual power that the set train plus storage system exchange with the electrical traction network; it is positive if it is absorbed from the grid. P train re f Train electrical reference power (considering the electromechanical efficiency); this power is obtained directly using the efficiency coefficients from the electrical machines and the input mechanical power; this power includes also the auxiliary power, and it refers to the catenary level. P STO rated,c , P STO rated,d Energy storage device rated power in charging and discharging modes, respectively. P STO max Maximum power at each instant that can be extracted from the storage system; this power is at the accumulation system side; it does not consider the electrochemical efficiency of the process. P STO max,av Maximum power at each instant that can be extracted from or injected into the storage system; it only depends on the energy stored and the electrochemical efficiency, not on the network conditions. P STO av Real available power that the storage system can inject into or extract from the system in discharging and charging mode. V cat Catenary voltage.