Novel method for power system synthesis load identification

Funding information National Natural Science Foundation of China, Grant/Award Number: 51577017; Fundamental Research Funds for the Central Universities, Grant/Award Number: 2019CDXYDQ0010 Abstract Power load parameters are one of the most important parameters for power system model. The load has the characteristics of time-varying, diversity and non-linearity. Therefore, obtaining the load parameters accurately is crucial to the safe and stable operation of the power system. The research in this paper is based on the synthesis load model which takes into account the impedance of distribution network. The multi-fault curves identification method and chaos-whale algorithm are adopted to improve the accuracy and practicability of load parameter identification. First, the algorithm proposed in this paper applies chaos algorithm to whale algorithm, and uses the randomness and ergodicity of chaos algorithm to make up for the shortcomings of the whale algorithm that the solution accuracy is not high and it is easy to fall into local optimum. Secondly, the identification method of multifault curves is to select several different faults for parameter identification, and the average value of each group of results obtained is regarded as the final result. Then compared with single curve identification, the multi-fault curves method is more superior. Finally, the accuracy of the results is verified by other faults. Simulation results show that the method is reliable to a certain extent.


INTRODUCTION
The mathematical model of each element in the power system is the basis of system analysis and calculation, and the accuracy of the simulation calculation result depends on the quality of the model. The detailed mathematical models and modelling techniques of generators, speed control systems, transformers and other elements in power system analysis have been well developed [1]. But for a long time, as one of the important components of the power system, the load has the characteristics of complexity, time-varying and non-linearity, which determines the difficulty of establishing its mathematical model [2,3].Therefore, it is urgent to improve the accuracy of load model parameters.
From the perspective of identification algorithms, many scholars have made great contributions to the parameter identification of load models. The early parameter identification algorithms were mainly linear methods, such as the least-squares method [4]. These methods were mainly applicable to the analysis and research of linear models. With the improvement of the load model structure, the synthesis load model considering the impedance of the distribution network has been widely used [5,6]. Therefore, the linear methods used to identify the synthe-  [7]. A simplified model identification strategy based on parameter sensitivity analysis was proposed in the literature [8], and genetic algorithm was used for parameter identification. In literature [9], the crossover operation in a genetic algorithm was introduced into particle swarm optimization, which improved the convergence of traditional particle swarm optimization. Then the improved particle swarm optimization algorithm was applied to load identification. Through the sensitivity analysis of the load model parameters, literature [10] used a traditional genetic algorithm and the effectiveness and accuracy of the method were proved by simulation experiments. In recent years, with the development of power systems, measurement techniques and load structures have also continued to change. For the composite load model with distributed generation (CMPLDWG) model, scholars have adopted Levenberg-Marquardt algorithm [11], deep learning method [12], double-deep Q-learning method [13], the support vector machine (SVM) approach [14], and so on. Nowadays, many novel swarm intelligence algorithms continue to appear. Among them, the whale optimization algorithm is applied in all directions of the power system, such as transient stability constrained optimal power flow [15], parameter identification of PV model [16] and parameter identification of isolated Wind-Diesel power systems [17]. Therefore, this paper adopts the whale optimization algorithm to identify load parameters.
From the perspective of identification methods, the methods used in the above-mentioned literature are all single-curve identification methods, so the correctness of the results of the above-mentioned literature depends on the choice of intelligent algorithms. It was observed that multiple solutions were obtained when parameter identification was carried out many times for a certain disturbance data, which indicated that the solution of parameter identification was not unique. It was also found that if there were N disturbance data, a single curve identification of these data had N load models, and the parameter values of these models were not the same. As the recorded data and the number of identifications increased, many models were obtained. However, in the actual simulation calculation process, we usually use a more widely used load model to perform simulation calculations under multiple faults. The multicurves method in literature [18] was to identify the fault curves of different lines and average the obtained parameter values. But the power fitting effect of this literature is not very good.
The load is a part of the distribution network, and the correct selection of load model parameters can improve the accuracy of the simulation calculation of the distribution network. The main motivation and research objectives of our paper are described as follows: (1) To reduce the time of parameter identification and improve the accuracy of parameter identification, chaos strategy is applied to whale optimization algorithm from the aspect of identification algorithm. (2) To solve the problem of choosing a set of more typical and practical parameters from multiple sets of model parameters, this paper proposes a new identification method-multi-fault curves identification method.
The article is divided as follows: Section 2 introduces the structure of the synthesis load model. In Section 3, the chaoswhale optimization algorithm is proposed. Section 4 uses the chaos-whale optimization algorithm to identify the parameters of the load model from a single-curve and multi-fault curves respectively. Section 5 is the further simulation of the larger system. Section 6 is the conclusion.

THE SYNTHESIS LOAD MODEL
The distribution network is mainly composed of distribution transformers, distribution lines, reactive power compensation and loads. The existing load model does not consider the influence of the distribution system, the synthesis load model proposed in the literature [4] considering the distribution network is closer to the actual situation of the power grid, which regards the transformer and line between the load and the bus as the impedance of the distribution network. And the induction motor, static load and reactive power compensation of the distribution network are connected to the virtual bus [5]. The structure of the synthesis load model is shown in Figure 1.

The static load model
The static load model represents the functional relationship between the load power and the bus voltage at a certain time [19]. The polynomial expression of the static load is shown as follows: where P 0 , Q 0 are respectively the initial active and reactive power of the distribution network branches;Z P is proportional coefficients of the constant impedance in the static active load; I P is proportional coefficients of the constant current in the static active load; P P is proportional coefficients of the constant power in the static active load; Z Q is proportional coefficients of the constant impedance in the static reactive load; I Q is proportional coefficients of the constant current in the static reactive load; Q Q is proportional coefficients of the constant power in the static reactive load.

The dynamic load model
Dynamic load models usually use induction motor models. The equivalent circuit of the induction motor is shown in Figure 2: In Figure 2, R S represents the stator resistance, X S is the stator reactance, R r is the rotor resistance, X r is the rotor reactance, and X m is the excitation reactance.
Induction motor models are generally divided into first-order, third-order, and fifth-order models. The first-order induction motor model considers the mechanical transient process. The where, E d ′ andE q ′ are respectively the d-axis and q-axis transient potential; I d and I q are respectively the d-axis and q-axis stator current; is the speed of the rotor; T E is the electromagnetic torque; T M is the mechanical torque.
The whole synthesis load model involves too many parameters. The parameters to be identified are summarized in Table 1.
In Table 1, P MP is the proportion of dynamic load, S 0 represents the initial slip of the induction motor.

CHAOS-WHALE ALGORITHM
Whale optimization algorithm (WOA) is a swarm intelligence optimization algorithm that imitates whale predation proposed by Griffith University scholars Seyedali Mirjalili and Andrew Lewis [20]. The WOA mainly has three behaviours of randomly searching, surrounding and attacking prey. The basic whale optimization algorithm has low accuracy and is easy to fall into local optimal. Chaos optimization strategy has randomness, ergodicity and sensitivity to initial conditions to achieve global optimization, so chaos optimization strategy is introduced into the whale optimization algorithm to make up for the shortcomings of whale optimization algorithm to improve the ability of global optimization.

Initial population
To improve the dependence of the WOA algorithm on the initial population, Sine chaotic mapping is used to replace the random generation method, which makes the distribution of the initial population more reasonable and increases the diversity of the initial population. Sine chaotic map population initialization steps are as follows: 2. The chaotic variables and population variables of the next iteration are as follows: Among them, i is the size of the population, X max and X min are the upper and lower limits of the population variables respectively; a′ is the control coefficient, which can be taken as [0,4].

Randomly looking for prey
Based on the value of | A |, the whale will choose one from the behaviour of randomly looking for prey or shrinking to encircle the prey as the next action. When | A | > 1, the random search for prey is conducted. This is a global search process, and the search equation is as follows: where, t represents the current iteration; X rand is a randomly selected position from the current population; r is a random number in [0,1].
In the basic whale algorithm, a is a linearly decreasing value from 2 to 0, which ensures the exploration ability at the early stage of the iteration and the development ability at the later stage. However, in the process of the load model parameter identification, there are a lot of parameters to identify, which requires the algorithm to have a higher exploration ability, and ensures that the algorithm performs as many global searches as possible in the early stage. The linear convergence factor is changed to a nonlinear convergence factor. The equation is as follows: where, Max_iter is the maximum number of iterations.

Shrinking surrounding prey
The behaviour of shrinking and encircling prey and attacking prey are both local optimization. The inertia weights in the particle swarm optimization algorithm are used to the whale algorithm. This article introduces the weights of the adaptive adjustment to balance the global searchability and the local optimization ability to prevent falling into local optimum. Therefore, this paper proposes an adaptive adjustment weight for nonlinear changes. The weighting equation is as follows: Each whale position in the search space represents a solution. It is now assumed that the prey position is the optimal position vector. The individual whales in the group shrink to the optimal position to update their position. The mathematical equation for shrinking and surrounding the prey is as follows: Where, X * (t ) is the position vector of the best solution obtained.

Attacking prey
Whales attack in the form of a spiral, which is based on the distance between the head whale and its prey.
Among them, D* represents the distance between the whale and its prey, b is the constant of the logarithmic spiral shape, and l is the random number in [−1,1].
It is worth noting that the head whale approaches its prey in a contracted enclosure while swimming along a spiral path. Although this behaviour occurs at the same time, let p be the selection probability, that is, when p is less than p 0 , the first expression in Equation (20) is executed and when p is greater than p 0 , the second expression in Equation (20) is executed, thereby updating the position of the head whale.
Where p is a random number between [0,1].

Chaos search optimization for the best individual
To further improve the solution accuracy of the whale optimization algorithm, the best optimization results of the above algorithm are retained, and then a chaotic optimization search is carried out for them [16]. The steps are as follows: 1. Set the maximum number of chaotic searches and reserve the best optimal individuals in the population; 2. The chaotic optimization search is performed on the best optimal individuals, and the formula is as follows: Among them, X is the best optimal individuals in the population, and X 2 is the result of chaos optimization. 3. Calculate the objective function values f(X) and f(X 2 ) of X and X 2 . If f(X) > f(X 2 ), then the optimal individual is X 2 , otherwise, the optimal individual is X.

The flowchart of the algorithm
The flowchart of the improved whale optimization algorithm is shown in Figure 3.

The selection of parameters to be identified
Experimental simulations were carried out in PSD-BPA and MATLAB. The small power grid of the IEEE30 structure shown in Figure 4 was built.
The synthesis load model mainly includes static load model, dynamic load model and distribution network impedance. The static load model in PSD-BPA uses a polynomial form, and the dynamic load model uses a third-order induction motor model.  Since there are many parameters to be identified, this will increase the identification time and reduce the identification accuracy. According to the sensitivity of the parameter [21,22], for the motor model, we choose the parameter with high sensitivity as the final parameter to be identified.
We used the perturbation method for sensitivity analysis. The perturbation method was based on the simulation system, and the dynamic response curve of the load model was observed after a certain parameter was adjusted by the control variable method. The step length of the perturbation method was set as 5%, and the fault was set as a permanent three-phase shortcircuit fault. To improve accuracy, the median method was used to calculate relative value of the trajectory sensitivity [22]. The calculation formula is as follows: The sensitivity curves of each parameter of the motor are shown in Figure 5.
As can be seen from Figure 5, the parameter R r , P MP and S 0 have a high sensitivity, so we take them as the final parameters  to be identified The parameters with less sensitivity are selected as the non-key parameters, and the typical value is adopted, as shown in Table 2.
After the remaining nine parameters are selected as the final identification parameters, the range of the identification parameters needs to be determined. The choice of identification range will also have a greater influence on the identification results. Considering the rationality of the parameters, the parameter identification range is shown in Table 3.

Multi-fault curves identification method
The load is time-varying, so it is now assumed that the load composition ratio is approximately the same, that is to say, the load characteristics are similar. On the premise of this assumption,

Selection of multi-fault curves
Three types of faults are set between bus 10 and bus 21, which include: Fault 1: A permanent three-phase short-circuit fault Fault 2: A permanent single-phase short-circuit fault Fault 3: A phase-to-phase short-circuit fault.
The voltage curves of the three faults are shown in Figure 6 respectively.
In order to prove the applicability of the identification results, the selected test faults are not only the fault between bus 10 and bus 21, but also the faults of other lines. The three tested faults are considered for the simulation, which include: Test 1: A phase-to-phase short-circuit fault on the line between bus 25 and bus 21. Test 2: A permanent single-phase short-circuit fault on the line between bus 2 and bus 6. Test 3: An instantaneous single-phase short-circuit fault on the line between bus 10 and bus 21.
The voltage curves of these faults are shown in Figure 7:

Determination of model parameters and objective function
The above three kinds of fault curves fitting will produce three sets of parameters, and the final model parameters can use the parameter averaging method. The specific equation is as Among them, A is the final parameter to be identified, and A1, A2, and A3 are the parameters identified under the above three fault conditions. The objective function for the synthesis load parameter identification is as follows: is the error function; P si is the measured active power; P mi is the fitted active power; Q si is the measured reactive power; Q mi is the fitted reactive power.

Analysis of simulation results
To visually show the final fitting effect of single-curve identification and multi-fault curves identification, the fitting situation of the identification result and the measured result is given by the following chart. The parameter results of single-curve identification and multi-fault curves identification are listed in the following tables.
To compare the results of single-curve identification and multi-fault curves identification more intuitively and clearly, the power fitting curves of the three faults are shown in Figure 8 respectively.
From Tables 4-6, it can be concluded that the error of the multi-fault curves identification results is smaller than the singlecurve identification error in the three fault situations. It can be seen from Figure 8 that the multi-fault curves identification effect is significantly better than the single-curve identification    The power fitting curves. (a) Active power fitting curves for test1, (b) reactive power fitting curves for test 1, (c) active power fitting curves for test 2, (d) reactive power fitting curves for test 2, (e) active power fitting curves for test 3, (f) reactive power fitting curves for test 3 in the case of the permanent three-phase short-circuit fault. The multi-fault curves identification effect is better at the peak or trough of the curve in the case of the permanent single-phase short-circuit fault. The multi-fault curves recognition effect of the first pendulum is better than that of single-curve in the case of the phase-to-phase short-circuit fault.
To further prove that the method has certain practicability, the three test faults are tested using the identification results of multi-fault curves. The effect of active and reactive power fitting is shown below: In order to see the fitting effect of Figure 9 more intuitively, the parameter results of multi-fault curves identification and the errors of the three test faults are listed in the Table 7.
The errors of the above three test faults are respectively 0.02599, 0.014357 and 0.009013. It can be seen from the error results and the power fitting curves that the fitting effect of the multi-fault curves is good, which proves the practicability and accuracy of this method.

FURTHER SIMULATION OF A LARGER NETWORK
To further prove that the method has certain practicability, simulations of the larger network and the distribution network are implemented.

5.1
The simulation of IEEE118 system In section 4, only the IEEE30 system was simulated. In order to prove that the identification algorithm and method proposed in this paper are applicable to a larger network, simulation will be carried out in IEEE118 system. The IEEE118 system is shown in Figure 10. The basic steps are the same as those in Section 4. Three types of faults are set between bus 49 and bus 54, which include: The errors of single-curve identification and multi-fault curves identification are shown in Table 8.
From Figure 11, it can be seen that regardless of single-curve identification or multi-fault curves identification, the curve trend is roughly the same. However, the effect of single-curve identification is not very good at peaks and valleys compared to multi-fault curves identification. Therefore, the method proposed in this paper is also applicable to the parameter identification of the synthesis load of a bigger system. It can be seen from Table 8 that the error of the multi-fault curves identification result is much smaller than that of the single-curve identification result.

5.2
The simulation of IEEE300 system The simulation results of the IEEE 30 system and the IEEE 118 system show that the multi-fault curves identification has certain applicability. In order to prove that the method can be applied to larger systems, the IEEE300 bus system is simulated. Three faults are set between bus 35 and bus 72. which include: Fault 7: A permanent three-phase short-circuit fault Fault 8: A permanent single-phase short-circuit fault Fault 9: A phase-to-phase short-circuit fault.
The power fitting curves of the three faults are shown respectively in Figure 12.
It can be concluded from the Figure 12 that although the curves of the two identification methods are similar, the multifault curves identification method is more similar to the actual curve. No matter what kind of fault, during the occurrence of the fault, the peak value of single-curve identification is much larger than the actual result, and the valley value is much smaller than the actual result. The results of multi-fault curves identification in peak and valley values are much better than that of single-curve identification in any fault case.

The simulation of the unbalanced 3-phase distribution feeder
In order to further prove the applicable scope of the method proposed in this paper, this part will be simulated in the power distribution system. The unbalanced 3-phase distribution feeder is shown in Figure 13.
Three faults are set between bus 632 and bus 671. which include: Fault 10: A permanent three-phase short-circuit fault Fault 11: A permanent single-phase short-circuit fault Fault 12: A phase-to-phase short-circuit fault.
Since this identification is a simulation of the three-phase unbalanced system and each phase has different loads, the final identification results are the fitting results of each phase respectively. For each phase, the power fitting curves of the three faults are shown in Figures 14-16, respectively. identification is more consistent with the actual situation than the single-curve identification.

5.4
The simulation of the 118 bus radial distribution feeder In Section 5.3, the unbalanced 3-phase distribution feeder was simulated. To verify that the method proposed in this paper can be applied to a larger distribution system, the selected distribution system is the 118 bus radial distribution feeder [23,24], and is shown in Figure 17.
Three types of faults are set between bus 49 and bus 50, which include: To compare the results of single-curve identification and multi-fault curves identification more intuitively and clearly, the power fitting curves are shown in Figure 18 respectively. It can be clearly observed from Figure 18 that in the case of the permanent three-phase short-circuit fault, the effect of multi-fault curves identification is better during the process of peak drop. In the case of the permanent single-phase shortcircuit fault, the multi-fault curves identification in the first few swing cycles is closer to the actual curve. In the case of the phase-to-phase short-circuit fault, the multi-fault curves identification has a better fitting effect at peak and valley values.

CONCLUSION
In order to improve the accuracy and applicability of parameter identification, this paper proposes a multi-fault curves identification method, and the chaos-whale optimization algorithm is applied to the identification of multi-fault curves to minimize the mean square error between the calculated power value and the measured power value. Through the simulation calculation of IEEE 30 system, IEEE 118 system, IEEE 300 system and distribution network, including the unbalanced 3-phase distribution feeder and the 118 bus radial distribution feeder, selecting a series of short-circuit faults for simulation calculation and comparing with single-curve identification. It is found that this method is suitable for different types of net-works, and the result of simulation in IEEE standard system is better. Compared with single-curve identification, the method proposed in this paper is more effective in the fault stage. This method improves the accuracy of the simulation calculation of the distribution network, and ensures the safe and reliable operation of the power grid.