Research on the Influence of Traction Load on Transient Stability of Power Grid Based on Parameter Identification

Abstract

The traction load of electrified railways is a special user of the power system, which will have a significant impact on and pose a challenge to the power system after grid connection. Considering the impact and fluctuation of traction load connected to the power grid, this paper proposed a method to study the influence of traction load on the transient stability of the power grid based on parameter identification. Firstly, according to the operation characteristics of traction load, a comprehensive traction load model based on the measurement-based method of induction motor parallel power function model is selected, and the objective function of parameter identification is determined. Then, the specific steps of using the improved gray wolf optimizer (IGWO) algorithm to achieve the parameter identification of the traction comprehensive load model are introduced, including chaotic map, nonlinear convergence factor, and individual position update. Next, with the IEEE39 bus system as the network background, the load parameters identified by the measured data of a line of the Beijing–Shanghai high-speed railway are imported into the Power System Analysis Software Package (PSASP) 7.0 for transient simulation to analyze the impact of traction load on the transient stability of the power grid. Finally, three typical load models are compared with the traction load model to draw relevant conclusions, and it is verified that the traction load characteristics are different from the general load model, which provides reference for the connection of electrified railway traction load to the power grid and the economic and technical construction of the power system.

1. Introduction

With the continuous progress of industrialization, the loads in the power system have become more complex and diverse, and the number of nonlinear loads has sharply increased. Their operation will cause serious harm to the power quality of the grid, affecting the operation and management of the power system.

In reference [1], the Brazilian distribution system is used as an example to demonstrate that the penetration of distributed generation directly affects the power quality and will cause some harmonics to increase. In reference [2], the operating characteristics of two typical nonlinear load simulation models are established to conclude that the presence of a large number of nonlinear harmonic sources in the distribution network is the main factor affecting power quality. In reference [3], the grey clustering analysis algorithm is used to classify common nonlinear impact loads, and a power quality management plan is proposed for high-power variable frequency induction motors connected to the power grid. In reference [4], the Monte Carlo simulation is applied to calculate probabilistic power flow, and the matching degree between new energy output power and load power are found to have a significant impact on the voltage imbalance, power factor, and voltage deviation of the power system.

Electrified railway traction load, as one of the typical nonlinear impact loads [5], mainly consists of electric locomotive traction drive system and its auxiliary control power supply load. The traction load does not have energy supply capacity itself, and needs to obtain power from the power grid, it is characterized by the large range of load fluctuation and randomness; negative sequence current and harmonics [6]; low power factor and large energy consumption, etc. The development of technology has made possible massive transportation systems driven by renewable sources. The use of rail transportation is expected to increase [7], and the problems of power quality and stability fluctuation caused by it to the power grid cannot be ignored.

At present, the analysis of the influence of traction load on the safety and stability of the power system mainly revolves around the impact on voltage and the transient stability of the regional power grid caused by traction after grid connection.

In reference [8], the effect of the traction load model of DC locomotives on the transient stability of power grid after grid connection is analyzed from multiple perspectives, such as via load model comparison and the proportion of traction load in the power grid by using the time–domain simulation method. In reference [9], a modeling method of the impact load based on the power system simulation software BPA China version is proposed and applied to build simulation examples to analyze the effects of different types of impact loads on the transient stability of the power grid. In reference [10], the influence of concentrated integration with electrified railway traction load and photovoltaic power station on the stability are analyzed and evaluated. In reference [11], the influence of the commissioning of Lhasa–Nyingchi Railway on the security and stability of the weak power grid in the Tibet region are analyzed from three aspects: power flow analysis, stability analysis, and short-circuit capacity analysis.

By synthesizing the current research status at home and abroad, research on the safety and stability of traction load on the power grid mainly focused on the system power quality problems caused by traction load after grid connection [12], such as negative sequence current and harmonics, etc. However, there was little consideration given to the influence of the comprehensive load characteristics on the system operation status after grid connection. As an important load, traction load needs to be systematically analyzed for its dynamic characteristics on the transient stability of the power system [13], so as to provide support for the future planning, construction, and development of the relevant power and railway industries.

In view of the above problems, considering the dynamic and static characteristics of traction load and its impact and fluctuation after grid connection, this paper proposed an analysis method of the influence of traction load on transient stability of power grid based on parameter identification. By using the improved gray wolf optimizer, the measured data collected during the operation of high-speed railway traction load are combined with the theories of transient stability analysis to establish an accurate parameter model, which can better describe the relevant characteristics of traction load.

In PSASP, the IEEE39 node system is used as a simulation example to integrate the traction load model into power system, and the traction load model is compared with three typical load models for experiments by setting the corresponding fault disturbances; the effects of the traction load connection on the transient stability of the power system are analyzed in three dimensions: the relative power angle change of generators, change in bus voltage, and fault limit clearing time under different disturbances, which provided a reference value for the power grid and electrified railway operation.

2. Traction Comprehensive Load Model of Electrified Railway

2.1. Model Structure

The traction load of the electrified railway mainly refers to the traction load of electric locomotive, including traction drive system and auxiliary power supply control system. According to its characteristics, we selected the traction comprehensive load model of the induction motor parallel power function model proposed in reference [14] by using the measurement-based method. This model took into account the components of the traction power supply system comprehensively and established a “dynamic + static” load model based on three-phase power to accurately describe the dynamic and static characteristics of the traction load.

In the power system calculation, the induction motor model is the most commonly used dynamic load model. According to different parameter quantities, such as dynamic response accuracy and voltage stability index accuracy, the induction motor transient model is usually divided into electromagnetic transient, electromechanical transient, mechanical transient, and voltage transient models. The accuracy characteristics of the four types of model are shown in

Table 1. Comparison of four induction motor models.

Induction Motor Model	Electromagnetic Transient Model	Electromechanical Transient Model	Mechanical Transient Model	Voltage Transient Model
Active power response accuracy	Better	Better	Good	Bad
Reactive power response accuracy	Better	Better	Bad	Good
Stability accuracy	Better	Better	Better	Good
Parameter complexity	High	High	Low	Low

.

For static load characteristics, two forms of polynomial model Equation (1) and power function model Equation (2) are often used for description.

$$\left\{\begin{array}{l}P=P_0{\left(\frac{U}{U_0}\right)}^{p_v}{\left(\frac{f}{f_0}\right)}^{p_f}\\ Q=Q_0{\left(\frac{U}{U_0}\right)}^{q_v}{\left(\frac{f}{f_0}\right)}^{q_f}\end{array}\right.$$

(1)

where $P$ is the active power, $Q$ is the reactive power, $U$ is the voltage value, $U_0$ is the voltage base value, $f$ is the frequency value, $f_0$ is the frequency base value, $p_v$ and $q_v$ are, respectively, load active and reactive voltage coefficients, $p_f$ and $q_f$ are, respectively, load active and reactive frequency characteristic coefficients.

$$\left\{\begin{array}{l}P=P_0[a_p{\left(\frac{U}{U_0}\right)}^2+b_p\left(\frac{U}{U_0}\right)+c_p]\\ Q=Q_0[a_q{\left(\frac{Q}{Q_0}\right)}^2+b_q\left(\frac{Q}{Q_0}\right)+c_q]\end{array}\right.$$

(2)

where $a_p$, $b_p$, and $c_p$ are the percentage of active power of constant impedance load, constant current load, and constant power load as a percentage of total active power, respectively.

According to reference [15], the dynamic behavior trends of the electromagnetic transient model and the electromechanical transient model are basically consistent. In order to increase the calculation accuracy and reduce the number of parameters to be identified, we have chosen a third-order electromechanical transient motor model for the dynamic part, the state equation of which is shown in Equation (3). Meanwhile, the power function model with fewer parameters and suitable accuracy is chosen for the static part.

$$\left\{\begin{array}{l}{T}^{\prime }\frac{d{E}^{\prime }}{dt}=-E^{\prime }+C\cdot U\cos \delta \\ \frac{d\delta }{dt}=-C\cdot \frac{U\sin \delta }{T^{\prime }{E}^{\prime }}+\omega -{\omega }_s\\ M\frac{d\omega }{dt}=-\frac{U\cdot E^{\prime }\sin \delta }{X^{\prime }}-T_m\end{array}\right.$$

(3)

In Equation (3), $E^{\prime }=\sqrt{E_q^{\prime 2}+E_d^{\prime 2}}$, $\delta =-\mathrm{a}\mathrm{r}\mathrm{c}\mathrm{t}\mathrm{a}\mathrm{n}\left(E_d^{\prime }/E_q^{\prime }\right)$, $C=\left(X-X^{\prime }/X\right)$, $T^{\prime }=X^{\prime }{T}_{d0}^{\prime }/X$, where $E^{\prime }$ is the voltage behind transient electromotive force, $E_d^{\prime }$ and $E_q^{\prime }$ are the direct and quadrature axis rotor voltage, respectively, $\delta$ is the angle behind transient reactance, $\omega$ is the shaft speed, ${\omega }_s$ is the synchronous speed, $X$ is the open circuit reactance, $X^{\prime }$ is the transient reactance, $T_{d0}^{\prime }=\left(X_m+X_r\right)/R_r$ is the rotor open-circuit time constant, $X_m$ is the excitation reactance, $X_r$ is the rotor reactance, $R_r$ is the rotor resistance, $M=2H/{\omega }_s$ is the motor inertia time constant, $H$ is the inertia constant, $T_m$ is the applied load torque.

The structure of traction comprehensive load model for dynamic induction motor load parallel power function static load is shown in Figure 1.

2.2. Mathematical Description

The dynamic part of the load mainly included some major power components, such as high-power traction motors, power electronic converter devices, and other devices, which are represented by a third-order model of induction motor in polar coordinate form [16], and the equation of power is as follows:

$$\left\{\begin{array}{l}{P}_d=-\frac{U{E}^{\prime }}{X^{\prime }}\sin \delta \\ Q_d=\frac{U\left(U-E^{\prime }\cos \delta \right)}{X^{\prime }}\end{array}\right.$$

(4)

where $P_d$ and $Q_d$ are the active and reactive power of the dynamic part of the load, respectively.

The static part of the load mainly included the auxiliary and control circuits in the electric locomotive, such as lighting, heating, ventilation, water tank, and other equipment, which are expressed in the form of power function. Since frequency changes are generally not considered during normal operation of the power system, and frequency changes can be negligible in certain transient processes. Therefore, Equation (2) can be written as follows:

$$\left\{\begin{array}{l}{P}_s=P_{s0}{\left(\frac{U}{U_0}\right)}^{p_v}\\ Q_s=Q_{s0}{\left(\frac{U}{U_0}\right)}^{q_v}\end{array}\right.$$

(5)

where $P_s$ and $Q_s$ are, respectively, the active and reactive power of the static part of the load. $P_{s0}$ and $Q_{s0}$ are the active and reactive power under steady-state power flow before disturbance, respectively. $U_0$ is the steady-state voltage. $P_U,P_f,Q_U,Q_f$ represent the static characteristics of the load, and are the objects of parameter identification for load model.

The total power of the comprehensive load can be expressed as:

$$\left\{\begin{array}{l}P=P_s+P_d=P_{s0}{\left(\frac{U}{U_0}\right)}^{p_v}-\frac{U{E}^{\prime }}{X^{\prime }}\sin \delta \\ Q=Q_s+Q_d=Q_{s0}{\left(\frac{U}{U_0}\right)}^{q_v}+\frac{U\left(U-E^{\prime }\cos \delta \right)}{X^{\prime }}\end{array}\right.$$

(6)

As can be seen above, the parameters to be determined are composed of 5 dynamic load parameters ${\theta }_d=\left[T^{\prime },X^{\prime },C,M,T_m\right]$ and 2 static load parameters ${\theta }_s=\left[p_v,q_v\right]$, constituting a total of 7 parameters to be identified, and the set is represented as $\theta =\left[T^{\prime },X^{\prime },C,M,T_m,p_v,q_v\right]$.

The objective function of parameter identification is a positive monotonically increasing function of the output error, so the objective function in this model is the minimum value of the error sum of squares of the measured sample data and model calculated data:

$$E^{*}=\min \epsilon \left(\theta \right)=\frac{1}{N}\min {\displaystyle \sum _{k=1}^N\left[{\left(P_{Mk}-P_{Ck}\right)}^2+{\left(Q_{Mk}-Q_{Ck}\right)}^2\right]}$$

(7)

where $P_{Mk}$ and $Q_{Mk}$ are the true values of the measured active and reactive power; $P_{Ck}$ and $Q_{Ck}$ are the active and reactive power calculated by the requested model; and $N$ is the number of sample data groups.

2.3. Analysis Process

The analysis of transient stability problem of power system mainly studied whether the system can enter into a new stable operation state or has the ability to recover to the original stable operation state after a large disturbance under the originally stable initial operation conditions. In this paper, the analysis process of the impact of traction load on transient stability of the power grid after grid connection is shown in Figure 2.

We chose a typical IEEE39 bus system as an example. The system power flow under steady state is first calculated as the initial value of transient calculation, and the effects of generator governor, voltage regulator, and PSS are been considered. The parameters of electrical components such as generators, transformers, and power lines are used as typical parameters. The improved gray wolf optimizer is used to import and identify the traction load model from the actual measured data of high-speed railway operation. By comparing the transient stability of traction load with three other typical load models, the transient stability effects of traction load access to the grid are indirectly analyzed and relevant conclusions are drawn.

3. Parameter Identification Based on the Improved Gray Wolf Optimizer

3.1. Principle of Gray Wolf Optimizer

The gray wolf optimizer (GWO) is an algorithm proposed by Mirjalili [17], an Australian scholar, to imitate the process of siege and predation by gray wolves. The gray wolves can be divided into four levels according to their population status, which are arranged according to the pyramid pattern, and their number decreases layer by layer, as shown in Figure 3.

At the top of the pyramid are the leaders and managers of the gray wolf population, called α wolves, who are responsible for formulating the action guidelines and directions of the population; at the second level are the deputy leaders of the wolf population, called β wolves, who have partial decision-making and leadership power; at the third level are the coordinators of the group, called δ wolves, who are responsible for assisting α and β wolves; at the lowest level are the executors of the group, called ω wolves, who obey the orders of the other three levels of gray wolves and balance the internal relationship within the group.

When the gray wolf pack rounds up prey, α wolves are generally responsible for leading the pack, β wolves are in charge of command, δ wolves complete coordination, and ω wolves surround at the end. The position of the prey corresponds to the optimal solution of the objective function [18], and the wolf pack moves according to the location of the prey. The α, β, and δ wolves that surround the prey first correspond to the three solutions with the highest fitness, while the other wolves automatically adjust the surround range according to the positions of α, β, and δ wolves. After each prey update, the wolves will recalculate the fitness and automatically upgrade the three wolves closest to the prey to α, β, δ wolves, so as to approach the optimal solution [19].

The mechanism of individual wolves tracking prey orientation [20] is shown in Figure 4.

3.2. Improved Gray Wolf Optimizer Algorithm

The overall structure of the GWO algorithm is simple, requires fewer parameters to be adjusted, and is easy to implement [21]. In addition, it can automatically adjust the convergence factor and the information feedback mechanism, achieving a balance between global search and local development. However, like other swarm intelligence optimization algorithms, the GWO algorithm inevitably suffers from the problem of falling into the local optimum easily, and will lose the diversity of the population in the later stage of the search, which reduces the optimization accuracy.

In this paper, we took advantage of the stochastic and ergodic characteristics of chaotic systems, introduced a nonlinear adaptive convergence factor, and combined the idea of the typical population optimization algorithm, Particle Swarm Optimization (PSO), to update the positions of individual wolves, so as to balance the local search ability of the algorithm.

3.2.1. PWLCM Chaotic Mapping

In solving function optimization problems, the GWO algorithm usually uses randomly generated data as the initial population information, which will make it difficult to retain the diversity of the population and result in unsatisfactory optimization search results. However, chaotic motion is characterized by randomness, traversal, and regularity, which can enable the algorithm to easily escape from the local optimal solution, thus maintaining the population diversity and improving the global search capability [22].

When selecting chaotic mapping, two important characteristics of chaotic mapping must be considered, namely “simplicity” and “ergodicity”. Compared with other one-dimensional chaotic systems, the piecewise linear chaotic mapping is relatively uniform in phase distribution and simple in terms of equations, which satisfies these above two characteristics [23]. Therefore, in this paper, we adopted the chaotic property of PWLCM (also known as piecewise) chaotic mapping to replace the random initialization of the gray wolf population, so as to achieve a more uniform distribution of the gray wolf population in the search space, its iterative formula is described as follows:

$${\mathrm{x}}_{t+1}=\left\{\begin{array}{ll}\frac{x_1}{p}& x_t\in \left(0,p\right)\\ \frac{1-x_t}{1-p}& x_t\in \left(p,1\right)\end{array}\right.$$

(8)

where $p\in \left(0,1\right),x_t\in \left[0,1\right],\mathrm{a}\mathrm{n}\mathrm{d}n=1,2,\dots$

Different dimensional parameters were set in Matlab to simulate PWLCM chaotic mapping between (0, 1), and the results are shown in Figure 5. By comparison, at p = 0.4, the parameter variables have relatively good randomness, ergodicity, and regularity. Therefore, in this paper, the PWLCM chaotic mapping at p = 0.4 was taken to initialize the gray wolf population uniformly.

3.2.2. Nonlinear Convergence Factor

The GWO algorithm suffers from the problem of balancing the global and local searching ability, where the global searching ability plays an important role in the richness of the population and the local searching ability determines the convergence speed of the algorithm. Therefore, the variation in the convergence factor is crucial [24]. In the traditional GWO algorithm, the influence coefficient C, the convergence factor a, and convergence coefficient A is shown in the following equation:

$$C=2\cdot r_1$$

(9)

$$a=2-2\times \frac{x_k}{T_{\max }}$$

(10)

$$A=2a\cdot r_2-a$$

(11)

where $x_k$ is the current iteration ordinal number, $T_{max}$ is the maximum iteration ordinal number, $r_1$ and $r_2$ are random numbers between [0, 1]. Obviously, this convergence factor shows a linear decreasing trend and cannot achieve the purpose of global fast optimization. For this reason, an adaptive nonlinear convergence factor was introduced in this paper, whose mathematical model is shown in the following equation.

$$a=\frac{2}{\ln 2}\times \ln \left(2-{\left(\frac{t}{T_{\max }}\right)}^e\right)$$

(12)

To verify its effectiveness, we listed the convergence factor update strategies in references [25,26] in Equations (13) and (14). The dynamic changes of the four convergence factors are shown in Figure 6.

$$a=2\times \cos \left(\frac{t}{T_{\max }}\right)$$

(13)

$$a=2-2\left(\frac{1}{e-1}\times \left(e^{\frac{t}{T_{\max }}}-1\right)\right)$$

(14)

From the figure, it can be seen that compared with other update strategies, in the early stage, a should slowly decrease with the increase in iteration number to maximize the global search, which is conducive to avoid falling into the local optimum. In the middle and later iteration stage, a should converge quickly to search finely in a small range to ensure the accuracy of the optimization search.

3.2.3. Change Individual Location Update Method

In the traditional GWO algorithm, the position of gray wolf individuals is updated as:

$$X\left(t+1\right)=\frac{X_{\alpha }\left(t\right)+X_{\beta }\left(t\right)+X_{\delta }\left(t\right)}{3}$$

(15)

where $X_{\alpha }$, $X_{\beta }$, and $X_{\delta }$, respectively, represent the best position of α, β, and δ wolves.

In the process of individual position iteration, the equation only considers two factors, the current individual position and the historical optimal position, which leads to the lack of diversity in the gray wolf population.

In PSO, the optimal positions of all particles and the global optimum of the population can be preserved, which has good effects on faster convergence and avoiding prematurely falling into the local optimal solution. Therefore, in order to further improve the search capability and convergence speed of GWO algorithm, we combined with the idea of PSO algorithm and introduced the historical optimal solutions experienced by individual gray wolves into the position update formula [27] so that they could retain their own optimal position information. The improved position update formula is as follows:

$$X\left(t+1\right)=\frac{1}{n}\left(\begin{array}{l}{X}_{\alpha }+X_{\beta }+X_{\delta }+b_1\cdot r_1\left(X_{\alpha }-X_t\right)+\\ b_2\cdot r_2\left(X_{\beta }-X_t\right)+b_3\cdot r_3\left(X_{\delta }-X_t\right)\end{array}\right)$$

(16)

where $b_1,b_2,b_3,r_1,r_2,\mathrm{and}{r}_3$ are random numbers between [0, 1], $X_t$ is a random individual in the population, and $n$ is the number of leader layers in the gray wolf population, which is taken to be 3 in this paper. The individuals with better diversity generated in the iterative process can prevent the phenomenon of early maturation.

3.2.4. Algorithm Flow

The process of using the IGWO algorithm to identify the parameters of the traction load model is shown in Figure 7.

The specific steps of the algorithm are shown in

Table 2. Steps of the improved gray wolf population optimization algorithm.

Step	Description
1	Set the parameters of the algorithm, including the number of groups (population size) of parameters to be identified N, the number of parameters to be identified in each group (dimension) Dim, convergence factor a, control parameter A, swing factor C, the maximum number of iterations Tmax, and the objective function, set the current number of iterations t = 0
2	Initialize the individual positions of the gray wolf population using PWLCM chaotic mapping, generate N gray wolf individuals with dimension Dim
3	Calculate the fitness value fi corresponding to the current position of all gray wolf individuals (in this paper, we took the objective function value as the fitness value) and sort them to derive the first three fitness values in this calculation process, correspond them to the α, β, and δ wolves, and the corresponding positional information is recorded as Xα, Xβ, and Xδ
4	Calculate and update the value of the nonlinear convergence factor a according to (12), and then calculate the values of C and A according to Equations (9) and (11)
5	Update the gray wolf population individual locations based on (16), calculate new population fitness values, and reorder to update α, β, and δ values
6	Determine whether the algorithm has reached the maximum number of iterations Tmax. If the iteration reaches the maximum value, then output the population global optimal fitness value fbest and its corresponding population position coordinates Xbest, where fbest is the target solution of the optimization objective function and Xbest is the optimal set of identification parameters. Otherwise, go to step 3 to continue the iterative operation

.

3.3. Identification Results of Traction Load Model Parameters

The actual measured data of a line operation of the Beijing–Shanghai high-speed railway was taken as an example; the parameters of the traction comprehensive load model were identified by using the IGWO algorithm proposed in this paper. Through calculation, we obtained the model active power and reactive power fitting response curves, and compared them with the original data, as shown in Figure 8.

From Figure 8, it can be seen that the fitting curves of the calculated values of active power and reactive power of the traction load model are basically consistent with the data curves of the actual measured values, indicating that the identified traction load model can track the measured power response curve well.

According to the research in reference [28], the relative error coefficient can be used as an evaluation indicator for load model:

$$\epsilon ={\frac{\left(\frac{1}{N}{\displaystyle \sum _{k=1}^N{\left[y(k)-\stackrel{\wedge }{y(k)}\right]}^2}\right)}{{\left(\frac{1}{N}{\displaystyle \sum _{k=1}^N{\left[y(k)\right]}^2}\right)}^{\frac{1}{2}}}}^{\frac{1}{2}}$$

(17)

where $N$ is the total number of samples used for estimation and $y\left(k\right)$ and $\widehat{y}\left(k\right)$ are the measured and simulated (active and reactive) power for the kth sample, respectively. If the value is less than 0.05, the dynamic load model is said to be acceptable.

The parameter identification and result errors are shown in

Table 3. Results and errors of parameter identification.

Identification Parameters	Identification Results
T′	1.0645
X′	0.8516
C	1.1332
M	0.6967
Tm	0.5745
pv	1.1437
qv	1.0969
${\epsilon }_p$	1.626 × 10−2
${\epsilon }_q$	1.074 × 10−2
Obj.E*(θ)	3.841 × 10−4

.

From Table 3, it can be concluded that the identified parameter results are within a reasonable range, the relative error coefficient of active power ${\epsilon }_p$ is 1.63%, the relative error coefficient of reactive power ${\epsilon }_q$ is 1.07%, and both active power and reactive power errors are less than 0.05, which proves the validity of the parameters of the traction load model identified using the IGWO algorithm. Figure 9 shows the variation in relative error coefficients with population size for four algorithms of GWO, PSO, Sparrow Search Algorithm (SSA), and IGWO.

From Figure 9, it can be seen that with the increase in population size, the IGWO algorithm has a lower error value and smaller changes. The results indicate that the IGWO algorithm has better error stability and parameter identification precision, which can improve the accuracy of load modeling.

The variety rule of the global optimal solution during the iteration process is reflected by the global optimal fitness value, which represents the convergence characteristics of the algorithm. Figure 10 plots the iterative variation curves of the optimal fitness values of the four algorithms.

From Figure 10, it can be concluded that with the continuous iteration of the algorithm, IGWO converges faster than GWO, the fitness value decreases, and both SSA and IGWO demonstrate better search accuracy and convergence speed than GWO and PSO. However, it is easy for SSA to fall into local optimal solution problems, and the IGWO algorithm proposed in this paper has better computational accuracy and convergence.

4. Example Analysis

4.1. Simulation Network and Parameter Settings

In this study, we improved the IEEE39 bus system and used it as an analysis model. After repeated testing, this system can meet the operation of high-power loads. Connecting traction loads to this testing network will not cause large-scale power blackouts, so its ability to restore stability can be further analyzed.

There were 19 loads in the system, of which the load bus BUS21 was the main research object in the traction load transient calculation, and the loads corresponding to the rest of the buses were set as the constant impedance loads commonly used in the stability calculation by setting up a three-phase short circuit grounding fault at the exit of bus BUS21 to simulate the grid transient process, the short circuit moment of 1.00 s, fault removal moment of 1.23 s. The system unit system calculations chose the per unit value and the generator was the sixth-order synchronous generator model, and took into account the role of the automatic excitation regulators and speed regulators. The system structure is shown in Figure 11, with the fault location marked by the red arrow.

The load of BUS21 was set to be different load models, including the following: the identified traction load model; the constant impedance load model; the traditional comprehensive load model (40% constant power load in parallel with 60% constant impedance load model [29]); and the large industrial load model (induction motor in parallel with 40% constant impedance load model). Among them, the large industrial load model selected the typical parameters recommended by the IEEE Task Force Load Modeling Working Group [30], which are shown in

Table 4. Typical parameters of induction motors recommended by IEEE Task Force Load Modeling Working Group.

Parameters	Numerical Value	Parameters	Numerical Value
Rs	0.013	A	1
xs	0.067	B	0
xm	3.800	D	0
Rr	0.009	Inertia Constant H	1.5
xr	0.170	Load Ratio k	0.8

, and the units of the parameters in the table are all per unit value. The rest of the system loads were modeled as constant impedance loads.

After a three-phase short circuit grounding fault, we compared the changes in the relative power angle curves between the generators and the changes in the BUS21 bus voltage when the BUS21 bus load was the above four different load models.

4.2. Transient Stability Analysis

The stability of the power system mainly consists of three aspects: power angle stability, voltage stability, and frequency stability. When all three meet the stability conditions, the system is stable; otherwise, the system is unstable. In this paper, the specific criteria for transient stability are as follows, in which the voltage stability criterion adopts the engineering experience criterion [31]:

• Power angle stability: after a power system failure, the relative angle between any main generating units within the system does not exceed 180° with reduced amplitude oscillation.
• Voltage stability: after the disappearance of the power system fault, the bus voltage of each pilot in the system is restored to more than 0.8 per unit value, and the bus voltage is restored to below 0.75 per unit value of the time does not exceed 1 s.
• Frequency stability: after the power system failure, through the adoption of protective measures, the system frequency does not change substantially, and can be restored to the normal range in a short period of time without affecting the safe operation of the power system.

In order to compare the transient calculation parameters of the traction load with the other three load models more intuitively, we analyzed the generator relative power angle variation curves and the load bus voltage variation curves of the four load models under the fault calculation example described above.

4.2.1. Generator Relative Power Angle Variation Curves

Figure 12 shows the power angle variation curves of BUS35 generator relative to BUS36 generator when connected to four different load models.

From Figure 12, it can be seen that during the fault period, except for the industrial load model, the relative power angle amplitude of the traction load model is the largest, and the curve oscillation trend is close to that of the constant impedance load model; in the 3rd second after the fault is cleared, the traction load model oscillates normally with reduced amplitude like the other load models, and then it slowly decays to be within the range of the normal relative power angle.

4.2.2. Load Bus Voltage Variation Curves

Figure 13 shows the voltage variation curves of the BUS21 load bus when four different load models are connected.

From Figure 13, it can be seen that before and after the fault disappears, the voltage curves of the traction load model are significantly different from the other three load models, and, at the same time, the voltage drop increased during the fault period are not consistent, and the sensitivity of the models to the voltage disturbances was also different. The specific analysis results are as follows:

• After the 1.0 s fault occurs, the traction load voltage drop amplitude during the fault period is the smallest, and the traction load voltage drop is the largest for the traditional comprehensive load, which indicates that compared with other load models, the power quality changes resulting from the voltage drop characteristics caused by the traction load model during the fault period have less impact on the power equipment.
• After the 1.0 s fault occurs, in terms of bus voltage recovery speed, the constant impedance load and traditional comprehensive load model are faster, followed by the traction load model, and the large industrial load model is the slowest. To summarize, it can be seen that traction load is more sensitive to voltage disturbances, the frequency of bus voltage changes is lower, and there will be a short recovery process. After the disappearance of the fault in 1.23 s, compared with the traction load model, the bus voltage recovery of the other three load models is more stable, while the traction load model needs to go through several cycles to recover to the per unit value of the original bus voltage, which has the problem of relatively unfavorable voltage power quality.

The transient stability analysis calculation statements under the four load models are summarized as shown in

Table 5. Four load models’ transient stability analysis calculation statement.

Load Model	Stable Situation	Maximum Power Angle Difference	Time for Bus Voltage to Return to 0.8 p.u.	Maximum Frequency Difference	Amplitude of BUS21 Voltage Drop during Fault
Constant impedance loads	steadiness	91.377	1.75s	0.01255	87.43%
Traditional comprehensive load	steadiness	109.962	1.84s	0.01761	92.50%
Large industrial load	steadiness	183.226	2.22s	0.01954	78.93%
Traction load	steadiness	127.837	2.19s	0.01682	74.97%

, which records the important analysis parameters in the transient stability calculation operations.

From Table 5, it can be seen that the traction load has the smallest voltage drop in amplitude during fault period, but other integrated parameters are between the large industrial load and the conventional comprehensive load. Therefore, the traction load model is more prone to losing stability compared to the constant impedance model and conventional load model, and more likely to remain stable compared to the large industrial load model.

4.2.3. Fault Limit Removal Time Calculation

A summary of the fault limit clearing time for different load models under different faults obtained by PSASP simulation is shown in

Table 6. Different fault limit clearing times under different load models.

Load Model	Constant Impedance Load	Traditional Comprehensive Load	Large Industrial Load	Traction Load
Single-phase short circuit grounding	0.8460 s	0.4700 s	0.4374 s	0.6792 s
Two-phase short circuit	0.3474 s	0.2810 s	0.2920 s	0.2924 s
Two-phase short circuit grounding	0.3180 s	0.2730 s	0.2744 s	0.2974 s
Three-phase short circuit	0.3229 s	0.2396 s	0.2472 s	0.2820 s
Three-phase short circuit grounding	0.3414 s	0.2325 s	0.2471 s	0.2825 s
Three-phase disconnection	-	0.9634 s	1.2280 s	2.7106 s

.

From Table 6, it can be seen that the limit clearing time of the traction load model under any fault is between that of the constant impedance model and the traditional load model. Therefore, the transient characteristics of the traction load are different from the general load model. If the traditional load model is used to describe traction model, it will be too conservative, which is not conducive to the economy of the relay protection equipment; if the constant impedance model is used, it will be too optimistic, which is not conducive to the safe operation of the relay protection equipment.

5. Discussion

From the perspective of the dynamic characteristics of the load model, we studied the transient stability calculation of the power system connected to the traction load model. We used the IGWO algorithm to identify the model parameters based on measured data, and analyzed the impact of different load models in the IEEE39 bus system on the transient stability of the power system after grid connection. It was found that the traction load model exhibits different performance from typical load models.

The identification results indicated that the traction comprehensive load model differs greatly from the parameters recommended by the IEEE load modeling working group, and both the physical meaning and mathematical expression have changed, which is consistent with the description in reference [32], that is, the model parameters have a certain degree of “Inexplicability of physical meaning”. However, this paper concluded that the relative power angle amplitude variation in traction load model is better than that of the industrial load model, which is different from the conclusion in reference [8], that the power angle variation characteristics of traction load model are better than the constant impedance load model. This phenomenon may be due to the complexity of traction load itself, and there are usually certain differences in the identification parameters of the traction load in different time and regions. Therefore, we should establish a comprehensive load model that can accurately reflect the dynamic characteristics of traction load based on the measured data at the research site.

At present, research on traction comprehensive load models mainly focuses on the mechanistic load model of “induction motor parallel traction motor and constant impedance” [8,32,33]. We proposed the load model of third-order electromechanical transient induction motor parallel power function, which can achieve the improvement of calculation speed through the simplification of the model, and the improved IGWO algorithm is used for parameter identification, which can reflect the dynamic characteristics of traction load well. This provides a certain reference value for the model structure and parameter identification optimization algorithm selection of load modeling, lays a foundation for further research on the impact of traction load on the safe and stable operation of power system, and is of great significance for the economic and technical aspects of power system construction.

It is worth mentioning that the impact power on the power grid generated when the traction load drives into or out of the traction substation or passing neutral section at a certain instant cannot be ignored. The limitation of this paper is that the modeling did not consider the effects of other factors such as the negative sequence current caused by load impact characteristics.

6. Conclusions

In this paper, we identified the parameters of traction comprehensive load model based on the proposed IGWO algorithm and analyzed the effects of different load models in the IEEE39 bus system on the transient stability of the power system. The conclusions are as follows:

• During the fault process, except for the industrial load model, the generator relative power angle amplitude of the traction load model is the largest, and the curve oscillation trend is close to that of the constant impedance load model; when the fault disappears, the traction load oscillates at a reduced amplitude like other load models, and then slowly decays to within the range of the normal relative power angle;
• The speed of voltage recovery stability of the traction load model is between that of the constant impedance load and large industrial load model, and its voltage drop amplitude is the lowest, and although there is no obvious oscillatory process, it requires multiple cycles to recover to the original stable state;
• Under all fault cases, the fault limit clearing time of the traction load is between that of the constant impedance load model and the traditional load model, which is numerically closer to the constant impedance load, but in a real situation, if the constant impedance approach is used to describe the traction load, it is not conducive to the smooth operation of the relay protection equipment of the power system;
• Traction load is different from other typical load models, and its dynamic load characteristics should be fully considered during the modeling process in order to better simulate the actual situation of the traction load after the connection of the grid to the power grid, which is beneficial to the economy and technology of the power grid’s construction and provides a certain reference for subsequent research;
• Further research can be divided into three aspects.
  • Studying the coupling between load characteristics and impact characteristics; the dynamic part can adopt an electromagnetic transient model and consider components such as power electronic converters to establish a more accurate traction comprehensive load model;
  • Connecting a traction load model to the actual topology of a regional power grid for analysis;
  • Developing a traction load modeling and parameter identification system to make the transient stability analysis of power systems more convenient and intuitive.