Urban Rail System Modeling and Simulation Based on Dynamic Train Density

Abstract

To further improve the simulation calculation ability of urban rail traction systems during the peak operation period and provide an accurate and reliable simulation tool for the subsequent train schedule and energy storage system design, a multi-train circuit model with a bilateral power supply was established in this paper, and a power calculation algorithm based on dynamic train density was designed. The circuit topology in the model can be dynamically adjusted according to the number of trains to improve the operation rate. Based on the spatial and electrical data of a real section of the subway, the urban rail circuit model was built on the MATLAB platform, and the actual operation data of the subway was imported for verification. The experimental results show that the multi-train model can accurately reflect the influence of voltage fluctuations on the traction system under different train running conditions, and the results fit the actual operation conditions. By comparing the influence of different train intervals on the RBE (regenerative braking energy) utilization, the results show that the optimal RBE utilization rate can be achieved by adjusting the train interval in the peak period.

1. Introduction

Urban rail transit is an efficient, safe, and environmentally friendly urban transportation mode, but it also faces problems such as high energy consumption, high operating costs, and large environmental impacts [1]. In order to improve the energy utilization efficiency and reduce the energy consumption of urban rail transit, the research on energy-saving technology of rail transit has become an important topic [2]. To reduce the comprehensive cost of design and verification of energy-saving technologies, such as timetable optimization [3] and energy storage system strategy [4], and to improve the efficiency and safety of the development process, it is necessary to establish an accurate urban rail equivalent simulation model. Regenerative braking technology [5] is the most basic energy-saving technology in urban rail systems. The way urban rail is currently modeled depends on how the regenerative braking energy is used. Chen and Hou et al. [6,7], based on the kinematic model of the train, used an evolutionary algorithm or a dynamic planning method to improve the acceleration and deceleration overlap time of the opposing train, so as to achieve the maximum utilization of RBE. The model can quickly complete the strategy verification under a long sampling step. However, the overlap time study ignores the coupling relationship between multiple trains and traction substations, so it is impossible to analyze the electrical quantities of each node of the traction system under this mutual influence. The train simulation model based on the MATLAB/Simulink motor module [8,9] can realize the comprehensive modeling of the motor, power electronics, and traction system on the basis of kinematics. However, a single-motor model cannot simulate realistic magnitude urban rail power [10], while increasing the number of motors increases the simulation time significantly. In order to improve the power simulation speed, current mainstream studies, such as the study by Saleh and Zhu et al. [11,12], usually equate the train as a current source with the traction system components to form an equivalent circuit model and use the current injection method to perform the current calculation in order to obtain the RBE and traction system power flow distribution. Compared with the literature [9], Khodaparastan et al. [13] used the equivalent circuit model to effectively improve the simulation speed and time span, which can realize the simulation analysis for 24 h. Zhu [14] investigated the problem of parameter tuning of the traction system and the energy storage system under different train operation scenarios. Sun [15] further combined the circuit model with the schedule optimization, equating the overlap time method to the overlap current method, and obtained a more accurate direct utilization of RBE through iterative current calculation, but the iterative calculation requires a lot of time. Chen et al. [16] proposed a method of calculating the current distribution without iteration in the system with fewer nodes and obtained the higher-order equations on the required train voltages by collating the node equations, and then input the circuit parameters to obtain the accurate train voltage, which further improved the simulation speed; however, during the peak period of operation, the interval time between trains decreased, resulting in an increase in the frequency of changes in the number of trains in the interval. In order to study the railway traffic quality under different train densities, Janusz Woch constructed the SOUT (System for Evaluation of Track Array) system [17]. In the research based on the SOUT system, Jeremi Rychlewski [18] analyzed the influence of train delay on traffic quality in the peak period, which further proves the importance of intelligent management. The delay of the train may lead to longer timetable confusion. For the urban rail traction system, the temporary change of the operation strategy during the operation of the train may bring an extremely unstable influence to the traction system.

In previous studies, the research on the urban rail model coupled with kinematics and electrics has basically been based on the fixed circuit topology. The purpose is to improve the utilization rate of RBE under a specific train condition, but this ignores the dynamic process of the train entering and leaving the interval during the peak period of operation, in this dynamic process; the power flow calculation based on the current injection method should also switch the circuit topology and the current injection position with the actual position of the train to avoid affecting the power flow calculation and RBE utilization analysis of the subsequent short-interval train entering the interval. The main contribution is to establish a fixed number of train models, and on this basis, an urban rail traction power calculation algorithm based on dynamic train density is designed. During the non-peak period, the algorithm can control the model to switch the circuit topology to a single train model. During the peak period of operation, the circuit topology is switched by calculating the number of trains in the interval in real time, and the current injection position is dynamically synchronized with the actual position of the train, which improves the dynamic characteristics of the power system and the accuracy of power flow calculation. The dynamic train density model is constructed on the MATLAB/Simulink platform, which proves the applicability of the model in full-time simulation. Under the same basic timetable, the simulation of different train interval time is set up. The results show that the control peak train interval time is 149 s and 179 s, which can achieve 25–32% RBE utilization.

The rest of the section is structured as follows: Section 2 presents the basic train model, traction substation, and contact network model. Section 3 describes the effect of dynamic train movements in the peak section on the model structure and the corresponding calculation method. Section 4 shows the model validation and simulation results. Section 5 summarizes the whole paper.

2. Urban Rail Traction System Model

2.1. Train Model

Subway vehicles consist of power cars and trailers, each of which is equipped with four traction motors. The electromagnetic torque generated by the traction motors is transmitted to the wheels through a gearbox to produce forward motion of the wheels, and the total traction force of the train is equal to the sum of the forces generated by all the traction motors. As shown in Figure 1, the train kinematic model is modeled in reverse [13], where the train number i and the corresponding displacement data x_i are used to obtain the train velocity v and acceleration a through the derivation operation, and Newton’s second law (Equations (1) and (2)) is used to compute the traction force F_T and braking force F_B of the train.

$$v=\frac{dx}{dt}$$

(1)

$$F_T-F_B-F_W=m\frac{dv}{dt}$$

(2)

where F_W is the sum of the basic resistance per unit F_w0 and the additional resistance per unit F_wj, where the additional resistance consists of the ramp resistance and the curve resistance, which are calculated as Equations (3) and (4) [19], respectively.

$$F_{w0}=2.7551+0.014v+0.00075v^2$$

(3)

$$F_{wj}=g\sin \theta +g\frac{700}{R}$$

(4)

In order to simplify the complexity of the train system, reduce the amount and time of calculation, and improve the operational efficiency and stability of the model, the traction power is calculated directly using Equation (5) to calculate the traction force, and the braking power is calculated by Equation (6), where η_m is the motor efficiency and η_i is the inverter efficiency.

$$P_t=\frac{F_{T}v}{{\eta }_m{\eta }_i}$$

(5)

$$P_b=F_{B}v{\eta }_m{\eta }_i$$

(6)

In the current calculation, the train is modeled as a controlled current source, and its controlled signal is the current result of dividing the traction power calculated by the kinematic model with the train voltage. As shown in the parallel part of the controlled current source model in Figure 1, the on-board braking resistor is modeled as a series connection of a diode, an IGBT, and a braking resistor, and when the pantograph voltage VTrain exceeds the braking resistor activation threshold, the comparator generates a PWM control signal g to control the on-off of the IGBT to maintain the stability of the pantograph voltage.

2.2. Varying Catenary Resistance and Rail Resistance Models

As an example, if there are two trains in the section, according to the distance between the train and the traction substation, the variable catenary resistance can be expressed as Equations (7)–(9):

$$R_{uf}=(D_{ts2}-x_{t1})\cdot \Delta r$$

(7)

$$R_{us}=(x_{t1}-x_{t2})\cdot \Delta r$$

(8)

$$R_{ut}=(x_{t2}-D_{ts1})\cdot \Delta r$$

(9)

If there is only one train in the section, the variable catenary resistance is expressed as Equations (10) and (11):

$$R_{uf}=(D_{ts2}-x_{t1})\cdot \Delta r$$

(10)

$$R_{us}=(x_{t1}-D_{ts1})\cdot \Delta r$$

(11)

where Δr is the unit catenary resistance; D_ts1 and D_ts2 are the positions of the first and second traction stations; x_t1 and x_t2 are the real-time positions of train 1 and train 2, respectively; R_uf is the forward catenary resistance of the upline train; R_us is the catenary resistance between the upline double train or the backward catenary resistance of the upline train; R_ut is the backward catenary resistance of the upline train that leans backward; and in the same way, R_df, R_ds, and R_dt correspond to the forward, middle and backward catenary resistance of the downline train. The rail resistance is replaced by the unit resistance Δrail. Taking R_uf, R_ruf as an example, in order to intuitively reflect the real-time voltage fluctuation results of the catenary and rail, the catenary resistance and rail resistance are modeled as controlled voltage sources in Simulink, as shown in Figure 2; the control signal has a high degree of freedom, and the required algorithm can be designed on the basis of the basic control signal to expand the primary control function. Since the rail resistance is coupled with the trend of the catenary resistance, the two will be uniformly referred to as the catenary resistance in the subsequent discussion in order to better describe their change process.

2.3. Bilateral Power Supply Section Traction System Modeling

The bilateral power supply traction system consists of two traction substations and the contact network, as shown in Figure 3. The traction substations step down the 35 kV AC MV network through two 12-pulse transformers, each of which is Δ-connected on the network side, as well as Δ-connected and Y-connected on the valve side to reduce the harmonic component, and then rectify the full-wave rectifier to nominal 1500 V DC, which will be supplied to the contact network through the feeder circuit breakers. According to the number ratio of subway traction substation to passenger station in a city, the passenger station is set in the bilateral power supply system on average.

A section of the subway power supply in a certain city was selected as the research target. This section has three stations. The first and third stations are traction substations, which are set as Traction Substation 1 and Traction Substation 2, respectively. The second station is a step-down substation. Since it does not have energy interaction with the DC traction network under study, it will not be discussed in the model. Considering that the research is based on the peak period, to facilitate subsequent dynamic train density research, a circuit model topology with the highest train density was built, as shown in Figure 4. Due to the limitation of train safety distance and the ratio between the number of traction substations and passenger stations, the maximum number of trains in the same direction in the section is two; therefore, train 1 and train 2 were set up as double trains in the downline in the section, and train 3 and train 4 were set up as double trains in the upline in the section. U_ts1, U_ts2 and I_ts1, I_ts2 are, respectively, the output voltage and output current of the Traction Substation 1 and 2.

3. Circuit Topology Analysis and Dynamic Train Density Power Calculation Algorithm

3.1. Peak Period Traction System Circuit Topology Analysis

The specific change process of the number of trains in the peak area is shown in Figure 5. At t = t1, both train_j and train_j + 1 exist in area i, the catenary resistance is calculated by Equations (7)–(9), and P_Train_j and P_Train_j + 1 are injected into the corresponding trains into the controlled current source to carry out the calculation of the current. And when t = t2 moment, train_j leaves region i, and only train_j + 1 exists in the section; at this time, the peak state in the section changes to the low peak state in a shorter time range, and the catenary resistance should also be switched to the calculation of Equations (10) and (11). At the same time, the current injection of train_j + 1 should be switched to the position of P_Train_j at the moment of t1, which will provide the calculation of the subsequent import of train_j + 2 at the moment of t3 with a spare current source. According to the above analysis, if there is only one train on the one-way line, from the perspective of circuit calculation, the circuit topology can be constructed only by considering the resistance of the two lines before and after the train. However, in the circuit model, we need to think about how to deal with the third catenary resistance R_ut in this case. There are two solutions: the first solution is to retain R_ut and allocate the resistance of R_us to R_ut equally; the second solution is to set the equivalent current source model of R_ut and the second train to zero. The number of nodes on the line will be reduced and the circuit topology will be changed, and then the resistance of the line behind the train will be imported into R_us. In the circuit analysis at moments t1 and t2, as illustrated in Figure 5, it is evident that regardless of whether it is during the peak density of t1 or the lower density of t2, the value of R_uf is consistently determined by the position of the leading train within the interval. The value of R_us is influenced by both trains at t1 and is solely determined by the leading train at t2. On the other hand, R_ut is consistently determined by the position of the trailing train within the interval. Therefore, considering the coupling relationship between different catenary resistances and each train, the second scheme was used to design the subsequent algorithm, which can avoid the non-synchronous condition that the number of trains changes; the number of catenary resistances in the circuit topology remained unchanged; and the motion phenomenon and circuit calculation of the train were further unified.

The number of trains and catenary resistance are summarized as shown in

Table 1. Catenary resistance variations.

Direction	Variable Catenary Resistance	State
upline	$R_{ut}+R_{us}=R_{us},R_{uf}$	S1
	$R_{ut},R_{us},R_{uf}$	S2
	$R_{ut}+R_{us}=R_{us},R_{uf}$	S3
	$R_{Li}$	S4
downline	$R_{dt}+R_{ds}=R_{ds},R_{df}$	S5
	$R_{dt},R_{ds},R_{df}$	S6
	$R_{dt}+R_{ds}=R_{ds},R_{df}$	S7
	$R_{Li}$	S8

. Among them, S1 and S5 represent the state of only one train in the one-way interval, S2 and S6 represent the state of two trains, S3 and S7 represent the first train in the interval leaving the interval, S4 and S8 represent the state of no train in the interval, and R_Li is the maximum resistance of the line in the interval. Taking the unidirectional urban rail system as a background example, when in the off-peak period, there is only one train coming in one direction, so only the S1 and S5 states may occur, the R_ut and the rear train current source model is continuously set to zero, and R_us and R_uf are used for calculation. When the train leaves the zone in states S4 and S8, set R_us and all train current source models to zero, the value of R_uf is set to the maximum resistance of the line R_Li, and off-peak period urban rail system modeling can be realized.

When in the peak period, the single-direction train state is in S3, S7 and S2, S6 cycle. At this time, the switching of the catenary resistance, as shown in the above analysis of Figure 5, is consistent with the number of trains and the number of states to give the R_us, and the corresponding train current source model of the zero signal to switch the topology of the circuit can be realized in the peak period of the urban rail system modeling.

3.2. Dynamic Train Density Power Calculation Algorithm

Considering the applicability of the model in the whole period, a power calculation algorithm of traction system based on train density in the interval is proposed on the circuit model of fixed topology. The algorithm consists of two modules: the upper train number control module and the lower power flow calculation module. The control flow is shown in Figure 6. The spatial and electrical data behavior is very important to the development of the algorithm. The relationship between the train position x_tj and the selected substation position D_{ts_i} determines whether the train can participate in the power calculation of the lower module. The upper module receives the train diagram data and counts by judging whether the train position is in the selected interval [D_{ts_i}, D_{ts_i+2}]. According to the number of trains, the judgment of the catenary resistance state in Table 1 is completed, and the zero signal and the resistance signal are generated and transmitted to the lower module. The lower layer module switches the circuit topology according to the zeroing signal. If the current is in the non-peak period (states S1, S5, S4, and S8), shown in Figure 6, the circuit topology is switched to the non-peak period topology. The power of the train is divided by the train voltage to obtain the basic signal of the train current source model, which is imported into the model for power calculation. If the current is the peak period (states S2 and S6), shown in Figure 6, the power of the train x_tj and x_tj+1 is processed as above, and the current source models before and after are imported, respectively, and the power calculation is carried out. After the power calculation of the current time step is completed, the program will return to the train number judgment step. If the number of trains is reduced during the peak period (states S3 and S7), the upper module can provide the circuit topology switching signal to the lower module in advance to ensure that the train power data can be imported into the correct train current source model, thus ensuring the accuracy of power calculation during non-peak and peak periods. After the circuit topology is switched, the power of the train x_tj+1 is processed with the above and imported into the front current source for power calculation, which realizes the application of the model in the whole period. The main steps are as follows:

• Import train diagram data, train parameters, and traction system parameters.
• Judge whether the train distance x_tj is in the selected interval [D_{ts_i}, D_{ts_i+2}] and count.
• According to the total number of trains, and the resistance state of the line shown in Table 1, the zero-setting signal and impedance signal are generated.
• If the state is S1, S5, S4, S8, the zero-signal switching circuit is used as a non-peak topology.
  • The basic signal of the train current source model is obtained by dividing the power of the train and the voltage of the train and import the model for power calculation.
  • After the calculation of the current time step is completed, return to step 2.
• If the state is S2 and S6, the zero-signal switching circuit is used as the peak high-train-density topology.
  • After the power of the train x_tj and x_tj+1 is processed as above, the power calculation is carried out by introducing the current source model before and after.
  • If the states are S3 and S7, the zero-signal switching circuit is used as the peak low-train-density topology.
  • After the power of the train x_tj+1 is processed as above, it is imported into the front current source for power calculation.
  • After the calculation of the current time step is completed, return to step 2.
• Cycle from steps 2 to 5 until the end of the simulation.

4. Simulation Verification and Analysis

4.1. Dynamic Train Density-Based Traction System Model Verification

In order to verify the proposed dynamic train density metro rail circuit model in detail, the spatial location distribution of a power supply section in a metro line project was taken as the research object, and the train density traction system model was built in MATLAB/Simulink 2022b. The relevant parameters of the traction substation and the rectifier transformer parameters are shown in

Table 2. Basic parameters of the traction substation.

Parameter	Value
Rated capacity of net side	2 × 25 MVA
Rated capacity of value side	2 × 3000 kVA
Net-side rated voltage	35 kV
Value side rated voltage	1210 V
No-load output voltage	1680 V
Equivalent internal resistance	0.03 Ω/km
Equivalent inductance	0.0025 F

, and the basic data of the subway B-type car [20] are shown in

Table 3. Basic data of subway B-type vehicles.

Parameter	Value
Train formation	4M2T
Average load	52.4 t
Average passenger	870
Total train weight	241.4 t
Max acceleration	1.0 m·s−2
Unit resistance	2.7551 + 0.014v + 0.00075 v2
Braking resistor limiting voltage	1800 V
Braking resistor starting voltage	1730 V

. The line parameters of traction system are shown in

Table 4. Catenary and rail parameters.

Parameter	Value
Distance between FDG and WJS	730 m
Distance between WJS and MWD	760 m
Catenary unit resistance	0.013 Ω/km
Rail line unit resistance	0.019 Ω/km

.

To verify the accuracy of the model, actual measurements with 2000 samples per second were used as a benchmark for comparison. The items for comparison included the traction substation output voltage and train current. Part of the operating map data was imported into the train kinematics model, and the output voltage and train current obtained after model calculation are shown in Figure 7. During the periods of 8–22 s and 26–40 s, the demand current rose with the acceleration process of the train. During the period 22–26 s, the current dropped to a smaller value due to the reduced traction of the train. After 40 s, the train entered the braking phase, and the regenerative braking generated a large amount of current, resulting in a rapid increase in the output voltage of the traction substation in Figure 7c. At the same time, the braking resistor was activated and stabilized the voltage at its starting voltage threshold at 50 s. The simulation and experimental results show that the traction system-related data obtained from the simulation are very close to the actual current and voltage data. Figure 7d–f show the comparison between the simulation and actual situation of another section of the train running for a single time, and the simulation and actual operation knot further shows that the traction calculation model has a high degree of reduction.

The operating diagram data [20] for the time period from 17:15 to 17:22 was selected as shown in Figure 8a, and the time window was set to 400 s for dynamic train density model validation. Setting Figure 8b as the operating map curve under the synchronous condition of up and down trains was used to simulate the distribution of RBE utilization under the unoptimized operating map. Figure 8c is the operating diagram with the down train interval delayed by 30 s on the basis of the synchronized condition of the up and down trains to simulate the RBE utilization in the case that the trains cannot run exactly according to the schedule during the peak operation period.

According to the running diagram shown in Figure 8a, the dynamic train density model was verified. The change trend of the resistance of the upward catenary is shown in Figure 9a–c, and the change trend of the resistance of the downward catenary is shown in Figure 9d–f. Figure 9g–l represent the train entry signs in the upline and downline sections, respectively. From 0 to 180 s, the upline train in the interval was only one train, and the front- and rear- catenary resistance R_uf and R_us of the train can normally realize the change of train position. Between 180 and 190 s, the upline train 2 entered the interval, and the three catenary resistors were enabled at the same time. After 190 s, train 1 left the interval, while R_uf did not change to 0 but rose to 0.18 Ω, indicating that train 2 successfully realized position switching in circuit calculation. After train 2 left the interval in 340 s, train 3 repeated the calculation process of train 1. In summary, the proposed model effectively achieves the design goals of circuit topology changes and train cycle calculations during peak hours. The trend of the resistance of the downline was the same as that of the upline.

In the traction system current calculation based on the current injection method, the power of the train is usually injected and then divided by the real-time voltage values at the two ends of the controlled current sources as current injection. Figure 10a,b show the real-time injected power of the train current source in the down and uplines, respectively, and Figure 10c,d show the signs of the presence of the respective train numbers in the sections in the down and uplines, respectively. By comparing the train presence signs, it can be observed that in the upline, when train 1 was alone in the section (5–160 s), the system automatically allocated power to the first current source for current calculation. When train 2 was in the section (160–180 s), the traction power of train 2 was arranged to the second current source. When train 1 left the section at 190 s, the power injection of train 2 was automatically switched to the first current source; thus, the injected power of the second current source was abruptly reduced to 0, which provided a current injection location for train 3 to enter at the following 310 s. The above simulation cases show that the proposed dynamic train density model is able to accomplish the synchronization of the injected current position with the train kinematics in the case of short train time intervals.

The output voltage U_ts and current I_ts of the traction substation comparing the dynamic train density model with the fixed topology model are shown in Figure 11. According to the train operation diagram shown in Figure 8a, it can be seen that in 0–40 s, the upline train 1 was close to U_Ts1 and started before the downline train 1, so the network voltage impact of U_Ts1 in Figure 11a was larger than that of U_Ts2 in Figure 11c. In 166–172 s, the dynamic train density model can simulate that the second train of the upline enters into the zone; the first train RBE was absorbed by the second train, which effectively suppressed the fluctuation of U_Ts1 and U_Ts2 in this period, while the fixed topology model was unable to simulate the traction system under the condition of changing number of trains in the peak period. After 250 s, as the number of trains in the interval changes frequently, the dynamic train density model can provide the traction system simulation results under complex train conditions. Figure 11b,d show the output currents of the traction substation under the dynamic train density model and fixed topology model models, respectively. The dynamic train density model can simulate the high load frequency of the traction substation in the case of frequent starting and stopping of trains during the peak period, and it demonstrates the feasibility of a good multi-train dynamic traction simulation.

4.2. Analysis of Operation Diagram and RBE Utilization Rate

In order to study the calculation effect of the proposed model on the RBE utilization rate at different intervals, multiple rounds of simulation were carried out based on the three operation diagrams shown in Figure 8. The RBE utilization rate is calculated by Equation (12):

$$\eta =\frac{{\displaystyle \sum _{j=1}^2{\displaystyle {\int }_{j\_s}^{j\_e}{U}_{dc\_j}{i}_{dc\_j}dt}-{\displaystyle \sum _{i=1}^4{\displaystyle {\int }_{t\_ts}^{t\_te}{P}_{t}dt}}+{\displaystyle \sum _{k=1}^6{\displaystyle {\int }_{t\_bs}^{t\_be}{i}_{re}^2\Delta R_{L}dt}}}}{{\displaystyle \sum _{i=1}^4{\displaystyle {\int }_{t\_bs}^{t\_be}{P}_{b}dt}}}$$

(12)

where the total RBE energy is calculated as the integral sum of all train braking powers P_b, the total catenary resistance loss is calculated as the sum of the losses of the six sections, and the direct utilization of RBE by the train reduces the total energy of the traction substation, so the RBE utilization rate η is calculated by the difference between the integral of the total train traction power P_t and the integral of the power of all the traction stations P_{dc_j}, divided by the total RBE energy.

Figure 12 shows the time-domain distribution of the influence of different train intervals on the utilization rate of RBE on the basis of the preset train diagram. With 159 s as the interval center term and 10 s as the tolerance, a total of eight different interval time differences were extended in two directions, thus forming nine train interval sequences. In Figure 12, the interval of 159 s is expressed as sequence 5, and the interval of 149 s and 179 s are expressed as sequence 4 and sequence 7. Case 1, case 2, and case 3 represent the RBE utilization distribution under the operation diagram shown in Figure 8a–c, respectively. By comparing case 2 and case 1, it can be found that under ideal conditions, the unified optimization of train diagram can effectively improve the overall RBE utilization rate by 10%. In the case of a fixed first train diagram, by modifying the interval time of subsequent trains to 139, 149, or 179 s, the braking and traction train conditions can be overlapped, as shown in the additional 149 and 179 s of the train diagram curve shown in Figure 8b, thereby continuing to increase the RBE utilization rate to a maximum of 32%.

Case 3 is the distribution of RBE utilization rate obtained by simulating the operation of train delays during the peak operation period. Compared with the ideal result of case 2, the interval time of the second train in the actual train diagram was increased to 210 s, resulting in the subsequent 300–400 s. The RBE utilization rate was less than that of case 2. In addition, at about 200 s, sequence 4 and sequence 7 represent the influence of overlapping conditions of multiple trains in the same direction and overlapping conditions of multiple trains in the opposite direction on the utilization rate of RBE, respectively. The same-direction utilization was affected by the change of resistance loss of catenary, and the RBE utilization rate was reduced by 5% compared with the opposite overlap condition. If the train interval is less than 149 s and more than 179 s, the overlapping conditions cannot be utilized, so the RBE utilization rate is greatly reduced. At this time, the energy storage system is required to indirectly utilize and realize the indirect utilization of RBE energy and maintain the stability of the traction network voltage.

5. Conclusions

Aiming at the problem of train density variation in urban rail rush hours, as the existing research does not analyze the dynamic circuit topology in detail, a multi-train system circuit model based on the dynamic train density calculation algorithm was proposed in this paper. Taking a certain subway project as the object, the time domain simulation of multi-train dynamic operation in a fixed area was carried out, and the time domain distribution of RBE utilization under different interval times was analyzed. The results show that the RBE utilization can be increased by up to 32% by optimizing the interval time of 149 s and 179 s in the ideal case of rush hours. By comparing the simulation results of the existing models with fixed number of trains in the whole period, the proposed dynamic train density model can expand the applicability of the existing model in the rush hours, and an effective simulation tool for the subsequent research on the combination of energy storage system strategy design and timetable optimization can be provided. However, the study focuses on realizing the dynamic train density function, and factors such as rail leakage current were not taken into account, so the calculation results have a slight error compared with the actual ones, and only the modeling in the bilateral power supply zones was considered. In the future, we will focus on improving the calculation accuracy of the model and expanding the model into a multi-zone simulation environment to better simulate the actual train operation. We will also conduct the study of cooperative voltage stabilization of multiple energy storage systems on the multi-zone model.