Multi-train modeling and simulation integrated with traction power supply solver using simplified Newton–Raphson method

Multi-train modeling and simulation plays a vital role in railway electrification during operation and planning phase. Study of peak power demand and energy consumed by each traction substation needs to be determined to verify that electrical energy flowing in its railway power feeding system is appropriate or not. Gauss–Seidel, conventional Newton–Raphson, and current injection methods are well-known and widely accepted as a tool for electrical power network solver in DC railway power supply study. In this paper, a simplified Newton–Raphson method has been proposed. The proposed method employs a set of current-balance equations at each electrical node instead of the conventional power-balance equation used in the conventional Newton–Raphson method. This concept can remarkably reduce execution time and computing complexity for multi-train simulation. To evaluate its use, Sukhumvit line of Bangkok transit system (BTS) of Thailand with 21.6-km line length and 22 passenger stopping stations is set as a test system. The multi-train simulation integrated with the proposed power network solver is developed to simulate 1-h operation service of selected 5-min headway. From the obtained results, the proposed method is more efficient with approximately 18 % faster than the conventional Newton–Raphson method and just over 6 % faster than the current injection method.


Introduction
In the recent decades, demand growth in public transport systems has increased rapidly. Several cities across the world have planned to develop their own urban mass transit systems or to extend their existing routes to cover every street corner. Most urban metro systems require DC traction power supply to energize their rail vehicles [1][2][3]. The third rail conductor in DC power feeding systems is typically used for urban metros with the standard DC supply voltage of 750 V. At higher voltage level, 1500 VDC or 3000 VDC, the overhead catenary feeding configuration is more appropriate. It is necessary to characterize electrical performance and power loading at traction substations for planning, designing, and operation of mass rapid transit. Multi-train system simulation [4][5][6][7] integrated with a power network solver is a potential tool to exhibit power supply performances.
DC railway power flow calculation has been continually developed. Some may consider that DC railway power flow is a reduced version of AC power flow. As AC power flow, Gauss-Seidel and Newton-Raphson methods [8][9][10] are both well-known and widely accepted. In DC railway power systems, these two methods have been commonly employed in case of non-linear traction power load. The nature of DC railway power system is as simple as DC linear circuits unless traction power load model is taken into account. From the literature [11][12][13][14] and also proof by simulation experiences, the current injection method or alternatively current-vector iterative method (CIM) is more efficient than any others. However, for a large-scale DC power network, several hundred nodes or up to a thousand nodes, incomplete Cholesky conjugate gradient (ICCG) method [15] to handle large sparse matrices is preferred.
In this paper, the well-known Newton-Raphson power flow method has been revised and therefore simplified. A set of current-balanced equations governing each node is selected instead of a set of conventional power-balanced equations. According to [16], the simplified Newton-Raphson method for AC power flow calculation has been proposed. With good performance and its simplicity, it is expected to have advantages over the above-mentioned methods.
This paper is organized into seven sections. Section 2 gives a brief in train movement and performance calculations. Section 3 is a review of DC railway power systems and its power flow solution methods. Section 4 describes Newton-Raphson methods for DC railway power systems in both power and current expressions. Section 5 illustrates multi-train system simulation integrated with power network solvers. Simulation results and conclusion remarks are in Sects. 6 and 7, respectively.

Train movement and performance calculation
The key dynamic variables of train movement are position, velocity, and acceleration rate. During single train motion, the relationships among these variables are only subject to the straightforward kinematic equation according to Newton's second law of motion [17][18][19]. Considering Fig. 1, a train vehicle is moving up on an inclined rail surface. This motion can be expressed mathematically by using the free body diagram describing all the forces acting on the train as shown in Eq. 1.
where F T denotes the tractive effort of the train, F R denotes the resistance force of the train, F RR denotes the rolling resistance force of the train, F dyna denotes the aerodynamic drag force of the train, F grad denotes the gradient force of the train, M eff denotes the total effective mass of the train, and a denotes the train acceleration.
There are resistance forces, F R , opposing the train motion. The resistance forces can be categorized into three types: frictional forces, commonly called ''rolling resistance'', air resistance or dynamic drag forces, and gravitational or gradient forces as described in Eq. 2.

Rolling resistance
The rolling resistance is the resistance to motion of rotating parts. It is mainly caused by frictional torques (bearing torques, gear teeth friction, brake pads). Equation 3 is the mathematical representation of the rolling resistance.
where W is the axle load, f R is the rolling resistance coefficient, f 0 and f 1 are two arbitrary constants, and v is the train velocity.

Aerodynamic drag resistance
The motion of a train takes place in the air and the force exerted by air on the train will influence the motion. The aerodynamic resistance force results from three basic effects: (i) the pressure difference in front and behind the train due to the separation of the air flow and the vortex creation behind the vehicle, (ii) skin friction representing the surface roughness of the vehicle body, and (iii) internal flow of air entering the internal parts of the vehicle. It is common to express the aerodynamic drag force in the basic form as Eq. 4: where q air is air density (kg/m 3 ), C d is an aerodynamic drag coefficient, A F is the projected frontal area of the train, and v air is the speed of air relative to the train body.

Gradient force
The gradient force on a slope will act in the opposite direction of the train for uphill. The positive and negative signs are for the downhill and uphill motions, respectively. The gradient force is a constant force as long as the slope is constant. Equation 5 is the mathematical representation of the gradient force.
whereg is the gravitational constant (9.81 m/s 2 );h is the slope angle. The train's aerodynamic and rolling resistance forces are the properties of each train vehicle. They can only vary with the train speed and can be separated from the gradient force for convenient use. The combination of these two resistance forces is called the drag force. Different  [20] is commonly used in Eq. 6.
where a, b, and c are drag-force coefficients. It is common to use different values for open or tunnel situations. Trains and their movements can be modeled mathematically. When the simulation starts, the first train is released on the track. It is in the powering mode to accelerate the train to reach the maximum train travel speed. When the maximum speed is reached, the train is switched to operate in many different ways due to the travel purpose e.g., energy saving operation. However, most of the train operating mode before the braking is always the coasting mode. The train's power can be computed by multiplying the tractive effort and the train speed. The tractive effort can be calculated by using the Newton's second law of motion which considers impacts of dynamic drag, gradient force, train's mass, and acceleration or deceleration as described previously.

Railway power supply study
The railway traction power supply system has several configuration features different to a normal power system, which are summarized as follows.

DC railway power supply system
The commonly used voltages for DC railways are 600 V, 750 V, and 1.5 kV for urban, interurban metros, and regional system, respectively [1][2][3]. Overhead catenary is typically used for light rail system at 600-800 V and for conventional interurban or regional systems at 1.5 or 3 kV. Because of the large currents involved compared to highvoltage AC system, the DC copper contact wire is made from heavier gage material. DC railway power supply system has several configuration features different from a normal DC power system. However, there are some simplifications of the power network modeling. DC feeding system feature includes a three-phase bridged silicon rectifier for conversion from alternating to direct current. Figure 2 is an example showing the structure of a DC feeding circuit connected to the nearest rectifier substation.
The DC railway power supply system is not a simple DC linear circuit due to two causes of non-linearity. The first is the rectifier substation that does not allow the current flowing in the negative direction. The second is the traction power of the train. The equivalent circuit of DC railway power supply systems is summarized in Fig. 3. The need for efficient DC railway power flow calculation is that although the total number of trains during service hours and the number of rectifier substations in DC mass transit systems can be up to a hundred nodes, it is computational burden when multi-train system simulation is considered. For a typical metro train service between 6.00 a.m. and 10.00 p.m. (16 h), the simulator with a time step of 0.5 s needs to recall the DC power flow solver for 2 9 16 9 3,600 = 115,200 times.
As shown in Fig. 3, the rectifier substation is modeled by Norton's equivalent source [13,14] in which I ss and R ss represent the Norton's short-circuit current and the Norton's resistance, respectively; R U1 and R U2 are the conductor resistances of the up-track sections; R D1 and R D2 are the conductor resistances of the down-track sections; P U1 and P D1 are the power consumptions of the running trains on the up-track and the down-track. The diode placed at the substation terminal is used to prevent any negative current flowing into the substation. The train model is represented by a controlled current model, I T = P T /V T . Hence, the DC power flow equation at bus k can be described as follows: where V k is the voltage at bus k, I SS,k is the short-circuit capacity of the substation at bus k, P k is the power load of the train connected at bus k, G k,i is an element k,i of the bus conductance matrix.
It should be noted that the bus conductance matrix G can be formulated as it is for the bus admittance matrix Y in complex power system analysis [8][9][10].

Gauss-Seidel DC railway power flow solution
The power flow equations as shown in Eq. 7 can be rearranged in order to obtain the updated voltage equation at a given bus k. This updated equation of iteration h ? 1 as expressed in Eq. 8 needs to be re-calculated bus-by-bus at each iteration. To accelerate the solution convergence, the accelerating factor, a f , can be applied. This is similar to that of AC power flow [8][9][10] as shown in Eq. 9.

Conventional Newton-Raphson DC railway power flow solution
Similar to AC power flow calculation [8][9][10], the updated voltage is calculated using Taylor series expansion of the power mismatches, as shown in Eq. 10. With first-order derivatives of Eq. 10, the Jacobian matrix, can be formulated. For DC railway power system, Eqs. 11 and 12 (for diagonal and off-diagonal elements, respectively) are used to compute elements of the Jacobian matrix. Therefore, the updated voltage at bus k of the h ? 1 iteration can be found in Eq. 13.

Current injection method for DC railway power flow solution
This method is based on the current-balance equation at each bus rather than the power-balanced equation [8][9][10]. From Eq. 7, it can be re-written into the current form as shown in Eq. 14. Although the power of trains and the rectifier property are non-linear, Eq. 14 can be transformed into a linear equation at each iteration when the initial or previously known voltage of all buses is assumed. This implies that P T,k /V k is iteratively arbitrary.

Current injection method for DC railway power flow solution
The simplified Newton-Raphson AC power flow method has been proposed in [16]. With its remarkable features, this technique can also be applied to solve DC railway power flow problems. It is based on the current-balance equations. Starting with Eq. 14 as the CIM, the updated voltage is considered by using Taylor series expansion of the current mismatches as shown in Eq. 16.
In comparison with the conventional Newton-Raphson method, Eq. 18 is simpler than Eq. 11. There is no loop (summation form) to compute the diagonal elements. In addition, Eq. 19 to compute the off-diagonal elements is constant throughout the calculation (voltage independent). It is computed only once before the iterative process is started. The Jacobian matrix of the simplified Newton-Raphson method needs to re-compute only the diagonal elements while Eq. 12 is voltage-dependent. All elements of the Jacobian matrix of the conventional Newton-Raphson method must be re-calculated at each iteration which causes slow execution times.

Multi-train system simulation
The multi-train system simulator (MTS) with discrete time update is used as the main simulation core. It can be programmed in various programming software, e.g., FOR-TRAN, C/C??, JAVA, MATLAB, etc. The metro train service can be implemented by using MTS. The interaction between MTS and power network solvers can be demonstrated as follows.
In Fig. 4, there are four main blocks: (i) System data, (ii) Multi-train simulator (main program), (iii) Network capture, and (iv) Power network solver. At each discrete time

Multi-train system simulator (main program)
Train  update, the multi-train simulator (main program) is used to simulate position and power consumption of all trains. The change in train positions and powers causes voltage variation in the power feeding system. The network capture will be called every time update to prepare bus data and line data for the power network solver, signal (a). After receiving signal (b) the multi-train simulator (main program) calls the power network solver, and signal (c) for power flow calculation. The power network solver sends signal (d) to the network capture in order to receive bus data and line data. After receiving the necessary data from the network capture, signal (e), the power network solver can perform the power flow calculation using the selected method. The voltage solution obtained by the power network solver is sent back to the multi-train simulator (main program), signal (f) and can be used to evaluate the train performances for the next time update (if required). This repetitive process will be performed until the stop time is reached.
To clarify the idea of MTS incorporating the network capture and power network solver, a system snapshot taken at a particular time step is depicted and illustrated in Fig. 5 to exhibit how the processes of MTS works.

Test system
A dense metro train service with a 5-minute headway is modeled for the simulation tests as shown in Fig. 6. It is Sukhumvit line (light green line) of Bangkok Transit System, so-called BTS Sky Train [21]. It consists of 22 passenger platforms and 10 rectifier substations. The trains receive electrical power from the 3rd conductor at 750 V. The rectifier substations are operated by BTS operator at no-load rectifier substation voltage of 900 V. The technical data [22] of BTS's EMUs (electric multiple units) are described in Fig. 7. The trains were all assumed to be identical. Figure 8 shows the distance-time curves from these tests using the multi-train system simulator described in the previous section. It exhibits the BTS service for onehour operation from the beginning, 6.00 a.m.. For more details, the speed-time trajectory for the first train, on the up-track is selected and shown in Fig. 9.
The test system is the BTS-Sukhumvit line of 21.6 km long. The test-case scenario is an hour operation starting from 6.00 a.m. having uniform 5-minute headway. The train's acceleration rate is set as 1.0 m/s 2 . The travel time of a one-way running train is 30 min and 36.5 s. The speed limit is assumed at 80 km/h for the entire route.

Simulation results
The test is concerned with DC metro train service. The system is examined by the multi-train system simulation coded in the MATLAB programming environment developed by the School of Electrical Engineering, Suranaree University of Technology, Thailand, to study the BTS-Sukhumvit line's train service, with uniform 5-minute headway. The effectiveness of SNR (simplified Newton-Raphson power flow method) compared with CNR (conventional Newton-Raphson power flow method), GSM (Gauss-Seidel power flow method), GSA (accelerated Gauss-Seidel power flow method), and CIM (current injection method) has been examined.
This test was performed on a Mac-book pro (Intel Core i5-2.8 GHz, DDR3 1600 MHz-4 GB) with MATLAB 7. With 1 9 10 -4 p.u. equally applied to the relative termination criterion (maximum power mismatch for the CNR, maximum current mismatch for the SNR, and maximum voltage error for the GSM, the GSA, and the CIM), their power flow solutions are compared. It reveals that the results obtained by the four power flow methods are exactly the same, but the total number of iterations required and execution times are different depending on their individual performances. Figure 10 shows the voltage profiles of the first train and the first rectifier substation. The power drawn from the first rectifier substation can be depicted in Fig. 11. It presents almost 3-MW of the peak power drawn from the first rectifier substation.
The convergence of the power flow methods for each test system is an essential indication to examine how the solution sequence moves toward the true solution and to show that the generated sequence is bounded. This roughly describes the rate of error reduction only and cannot be used to judge the computational speed of the calculation.
Thus, the execution times need to be observed carefully and also to be compared. In addition, the execution times for the four power flow methods applied to the test system are recorded and presented in Table 1.
The test was performed repeatedly for 30 trials per method. This can evaluate the effectiveness of each method. On the assessment of the overall execution time, it is perceived that the SNR is the fastest method while the CIM is the second and the GSM comes last. The average execution times are 8.19, 7.17, 6.72, 11.46, and 8.75 s for the CNR, CIM, SNR, GSM, and GSA, respectively. The optimal accelerating factor of the GSA is 1.37. This factor was obtained by varying the value between 1.0 and 1.6. Figure 12 summarizes the optimal tuning of the accelerating factor. However, with the optimal accelerating factor, the GSA leads to 8.75 s for the average execution time, it is 2.03 s slower than the best (SNR).
The complication of this study is due to the change in a total number of buses in the DC railway power supply  system. The BTS-Sukhumvit line has 10 rectifier substations. This implies that the minimum number of buses in the entire system is 10. After starting the service, the total bus number increases as the number of running trains increases. Figure 13 describes the total number of buses in the entire system at each particular time step.
The test result shows that the total number of buses can be varied in between 10 and 24. The maximum number of total buses is 24. This situation occurs when there is a total of 14 running trains in service. As can be seen in Fig. 13, in the steady train service, the total bus number of the entire system is in between 21 and 24.

Discussion
To exhibit the convergence rate of each power flow method, comparisons among convergence curves at a selected time, 2,500.0 s, of the four methods are shown in Figs. 14 and 15. The convergence curves describe their own convergence property. As can be seen, the quadratic convergence can be found in both Newton-Raphson power  The execution time of the proposed algorithm is faster than that of the conventional Newton-Raphson method. This execution time depends mainly on the amount of floating-point operations (FLOP) if these two methods are assumed to be performed in the same software on the same computing machine. It assumes that the other steps of these two Newton-Raphson methods are exactly the same, therefore, the Jacobian updating step dominates the overall execution time. The Jacobian updating formulae for the  conventional Newton-Raphson methods are described in Eqs. 11 and 12. It has a loop in Eq. 11 and its FLOP per iteration is N (total number of buses). In contrast, Eqs. 18 and 19 give the updating equations for the proposed method. Remarkably, there is no loop in both equations. Equation 18 may spend a few FLOP per iteration only and is independent from the total bus number. This is the main support of the proposed method to gain faster execution times compared to the conventional method.  Power network solver for DC railway power supply is one of the most essential parts in the multi-train simulation in order to analyze, simulate, design, and control the steadystate train service. Although there exist several powerful power flow solvers based on the conventional Newton-Raphson method, their problem formulation gives complication due to the need to calculate derivatives in the Jacobian matrix. The simplified Newton-Raphson method uses non-linear current mismatch equations instead of the commonly used power mismatches to simplify overall equation complexity. In addition, a total number of multiplication operations required by the conventional method is linearly proportional to the size of the Jacobian matrix, while that of the proposed method is nearly constant. This means that the calculation time of the conventional method increases more rapidly as the total bus number increases than that of the proposed method does. This can lead to improvement of multi-train simulation software development in fast computational speed and less memory usage.
To investigate the effectiveness of the proposed method, BTS-Sukhumvit line in Bangkok, Thailand is selected as the test system. It consists of 22 passenger stopping stations, 10 rectifier substation, and the total line length is 21.6-km. The rated voltage of the power feeding system is 750 VDC. The test-case scenario is assumed at the speed limit of 80 km/h with the 5-minute headway. The simulation is conducted with one-hour service operation. The train service with the given headway spends just over 30 min to reach the steady train service operation. The result shows that the total bus number of the entire system can vary in between 21 and 24 at the steady train service. As a result, the SNR shows the best performance among the four methods in the execution time. It is about 18 % and 6 % faster than the CNR and the CIM, respectively.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http:// creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.