Parameters Identification of Equivalent Model of Permanent Magnet Synchronous Generator (PMSG) Wind Farm Based on Analysis of Trajectory Sensitivity

As wind farms have great influences on power system stability, it is essential to develop an adaptive as well as robust equivalent model of it. In this paper, a detailed equivalent model of PMSG wind farm and initialization method is developed. The trajectory sensitivity of parameters is analyzed. Then, the key parameters are estimated using improved Genetic Learning Particle Swarm Optimization (GLPSO) hybrid algorithm with phasor measurement unit (PMU). The description and generalization capability, stability for parameter identification of the equivalent model under wake effects, and when some wind turbines are off-line or wind speed is unknown after an event are analyzed. The maximum differences between the values of estimated parameters and their real ones are less than 10% for the proportional magnification constant of DC voltage controller Kp2 and grid side current controller Kp3. The convergence rate and global optimization performance of the improved GLPSO hybrid algorithm is 0.5 times higher than the classical particle swarm optimization algorithm (PSO) and genetic algorithm (GA).


Introduction
Connecting a large scale PMSG wind farm to grid brings about great influence on angle [1,2], voltage [3,4] and frequency [5,6] stability of power system. If detailed model of PMSG wind farm is built to simulate all wind turbines, the complexity and simulation time will greatly increase [7]. Therefore, the establishment of a simplified, accurate and effective wind farm equivalent model is essential for analysis and control of power system [8].
Wind farm can be equivalent to a single machine [9,10] or multi-machines [11,12]. The equivalent methods mainly include aggregation method using detailed physical parameters [13] or parameters estimation method using PMU data in the point of common interconnection. Yingchen Zhang et al. report that the output from the standard model cannot match the measurements of he DFIG wind farm [14]. J. Brochu et al. carry out short circuit test on field and show that for stability and electromagnetic transient studies, single machine representation of wind farm after disturbance provide good and adequate results [15]. Jin Lin et al. analyze the wide-area trajectory sensitivity for parameters of constant speed induction generator wind farm. GA is used to estimate the parameters [16]. Xueyang Cheng et al. estimate the key parameters of the equivalent model of DFIG wind farm using PSO gradient search optimization with PMU data [17]. Yinfeng Wang et al. use improved GA to strengthen the capability of global optimization [18]. L. P. Kunjumuhammed et al. pointed out that aggregation Energies 2020, 13, 4607 2 of 17 methods must be improved to enhance precision [19]. Dong-Eok Kim et al. formulate a linear model for DFIG wind farm. The recursive least square method is used to estimate the parameters [20]. Yuhao Zhou et al. propose an aggregation method of collecting system. The robustness of the equivalent model for DFIG wind farm when some wind turbines are offline is analyzed [21]. Jian Zhang et al. consider the control system limiting units of grid side (GSC) and machine side converters (MSC) of DFIG wind turbine in the parameters estimation process. As an outcome, the precision are greatly improved [22]. Xueping Pan et al. propose a simplified model of PMSG wind turbine, in which the mechanical subsystem is neglected. The combination of PSO method and Levenberg-Marquardt method is used to estimate the simplified model parameters. Furthermore, the simplified PMSG wind turbine model and the classical load model are paralleled to illustrate the generalized synthesis load in the distribution network [23]. However, the permanent magnet flux linkage of PMSG is not identified. Moreover, the parameters of PMSG wind farm are identified with load parameters. Peipeng Zhou et al. formulate the frequency-domain impedance equivalent model of PMSG wind farms to cope with sub-synchronization and sup-synchronization oscillation problems [24]. Parameter sensitivity of wind generator input admittance to wind speed is analyzed. The input admittance and the participation factor for the dominant oscillation mode are used as the indexes of dividing groups. Based on the circuit equations, node-by-node aggregation is proposed. However, this approach still belongs to the aggregation method. After long-time running, the dynamical characteristics may be changed due to the change of parameters of each PMSG wind turbine. As a result, this method cannot be applicable any longer.
At present, a large number of researches focus on DFIG wind farm. The literature on equivalent model of PMSG wind farm is very rare. In this paper, in view of the fact that the traditional aggregation method cannot solve the problem of parameter variation during long-time running of wind farms, and manufacturers regard control system parameters as commercial secrets, we propose that the time-varying and control parameters with high sensitivity are identified using improved GLPSO hybrid algorithm while the values of other parameters are fixed to aggregated or typical values.
The organization of this paper is as follows. In Section 2, the initialization method of the equivalent model of a PMSG wind farm is proposed. In Section 3, the trajectory sensitivity analysis of time-varying and control system parameters is performed. In Section 4, parameters identification method is formulated. In Section 5, simulation cases are conducted. Conclusions are summarized in Section 6.

Initialization Method for the Equivalent Model of PMSG Wind Turbine
Suppose the steady state voltage magnitude, injected active and reactive power of common interconnection point of wind farm is denoted by V 0 , P 0 and Q 0 . Then current magnitude squared injecting into the grid from wind farm is given by: The active power loss on the coupling and filter inductor is given by: where R g is the resistance of coupling and filter inductor. The stator loss of PMSG is: where R s is the PMSG stator resistance. i qs is the q-axis current of stator of PMSG. The electromagnetic power of PMSG is: where E is the electromotive force of PMSG. ω e is the electrical angular velocity of rotor. ψ PM is permanent magnet flux linkage of PMSG. P EM can be rewritten as: Equations (3) and (4) are substituted into (5). Then Equation (6) can be obtained: Generally, the optimal torque control is adopted. The reference of electro-magnetic torque is generally set as: where k opt is the optimal torque-speed coefficient, which can be obtained according to the power-speed curve of PMSG wind turbine. The electromagnetic torque of PMSG is given by: In steady state, the electromagnetic torque is equal to that of the control command. That is: Substituting (7) and (8) into (9), results in: Equation (10) is substituted into (6). Then Equation (11) can be obtained: When solving (11), resistance R s can be neglected first. Thus, (12) can be obtained. Then, according to (13), the estimation of steady state rotational electric angular speed of rotor ω e0 can be obtained.
Taking ω e0 as the initial value, the Newton-Raphson approach is adopted to compute the rotational electric angular speed of rotor by solving Equation (11).

The Necessity of Trajectory Sensitivity Analysis
In the equivalent model, there are many parameters. If all parameters are identified, it will inevitably consume a lot of computational time, and the identification results will be very dispersive. As a consequence, the equivalent model is very difficult to apply. Therefore, only the key parameters Energies 2020, 13, 4607 4 of 17 that have strong impact on power system dynamics need to be identified. As direct current (DC) capacitors, filter inductors and equivalent distribution network impedance are static components, they can be fixed to nameplate or technical manuals values, and are not identified. The specific aggregation methods can be found in [18]. Moreover, PMSG is a rotating element. After long-term operation, its stator resistance, reactance and flux linkage may differ enormously from nameplate data and should be taken as candidate parameters to be identified. Although the inertial time constant is less time-varying, its value is usually not listed in the nameplate and technical manual. Therefore, it is taken as a candidate parameter to be identified as well. Controller parameters are generally regarded as commercial secrets by manufacturers. Their values cannot be found in nameplate and technical manuals. Therefore, controller parameters should also be taken as candidate parameters to be identified. However, low sensitivity parameters in candidate parameters to be identified are generally difficult to identify accurately and can be fixed to typical values. Therefore, sensitivity analysis of candidate identified parameters is carried out in this paper to discern the high sensitivity parameters.
The relationship between sensitivity analysis and parameter identification is that the former guides the latter. That is, sensitivity analysis tells which parameters should be identified. If all the candidate identification parameters are identified without sensitivity analysis, the convergence rate of the program will be very slow, and the identification results will be very dispersive.

Algorithm of Trajectory Sensitivity Analysis
The computing method of trajectory sensitivity and its average value in theory is shown in Appendix A. However, the trajectory sensitivity and its average value are very hard to compute analytically according to Equation (A1)-(A3) in Appendix A. However, the numerical solution can be obtained using simulation software. The algorithm is as follows.
(1) Firstly, a simulation example is built to set the candidate identification parameters to the nameplate or technical manuals values. If they cannot be found on the nameplate and technical manual, they are set as the typical values. (2) Then, the power flow is computed. After that, the transient analysis and computation are conducted. As a result, P 0 (t) and Q 0 (t) are recorded. (3) Increase the value of k th candidate identification parameter by 5% based on the value set in Step (1). Repeat step (2) to get P k (t) and Q k (t). (4) The trajectory sensitivity of the k th candidate identification parameter are computed as (5) The average trajectory sensitivities of the k th candidate identification parameter for active as well as reactive power are computed as where N is the number of sample data. (6) Repeat Steps (2)-(5) to compute the trajectory sensitivity and average trajectory sensitivity of all candidate identification parameters.

A Case of Trajectory Sensitivity Analysis
The trajectory sensitivity of the candidate identification parameters of the equivalent model is investigated using the simulation system shown in Figure 1. A round-rotor PMSG wind turbine with rated active power of 1.5 MW is connected to the 575 V bus. The parameters of physical and control system of the PMSG wind turbine are shown in Tables 1 and 2, respectively. Doctor Jiayang Ruan, who is the author of many published papers on wind power generation, developed the detailed PMSG wind turbine model according to Ref. [25] which can be downloaded from Ref. [26]. We modify it to the phasor model using MATLAB/SIMULINK. Energies 2020, 13, x FOR PEER REVIEW 5 of 18 rated active power of 1.5 MW is connected to the 575 V bus. The parameters of physical and control system of the PMSG wind turbine are shown in Tables 1 and 2, respectively. Doctor Jiayang Ruan, who is the author of many published papers on wind power generation, developed the detailed PMSG wind turbine model according to Ref. [25] which can be downloaded from Ref. [26]. We modify it to the phasor model using MATLAB/SIMULINK.  A three-phase short circuit fault with 8 and 1 ohms transition resistance, 4 cycles duration are respectively set to the midpoint of line L1. The simulation time lasts for 0.2 s. The voltage magnitude of bus B0 dips about 30% when the transition resistance is set to 8 ohm and 80% when the transition resistance is set to 1 ohm. The candidate identification parameters are increased by 5% in turn according to the values in Tables 1 and 2.
The average trajectory sensitivity of each candidate identification parameter for 8-ohm transition resistance is shown in Table 3   A three-phase short circuit fault with 8 and 1 ohms transition resistance, 4 cycles duration are respectively set to the midpoint of line L1. The simulation time lasts for 0.2 s. The voltage magnitude of bus B0 dips about 30% when the transition resistance is set to 8 ohm and 80% when the transition resistance is set to 1 ohm. The candidate identification parameters are increased by 5% in turn according to the values in Tables 1 and 2.
The average trajectory sensitivity of each candidate identification parameter for 8-ohm transition resistance is shown in Table 3. Clearly, the average trajectory sensitivity of stator reactance L d , L q , inertia time constant T j , proportional and integral magnifications of MSC control loop K p1 and K i1 are all 0. Therefore, they are difficult to identify. As a result, T j K p1 and K i1 can be fixed as typical values while L d and L q can be fixed as the values on nameplate or technical manuals. The trajectory sensitivities of active and reactive power for other candidate identification parameters are shown in Figures 2 and 3, respectively. It can be seen from the first figure in Figures 2a and 3a that, when the trajectory sensitivity of active (reactive) power for stator resistance R s is negative, that for permanent magnet flux linkage ψ PM is positive. Moreover, when the trajectory sensitivity of active (reactive) power for R s is at the trough, that for the permanent magnet flux linkage ψ PM is at the peak. That is, the phase of trajectory sensitivity of active and reactive power for the stator resistance R s is opposite to the flux linkage ψ PM . This is because the stator resistance is an energy consuming element, while the permanent magnet is the key element of energy conversion. Their characteristics are different. Since the phase of trajectory sensitivity of both active and reactive power for resistance and flux linkage are opposite, increasing the value of one parameter versus decreasing the value of the other parameter can lead to the same output active and reactive power changes. Therefore, it is difficult to accurately identify the two parameters at the same time. Further research shows that when the resistance and flux are involved in parameter identification process at the same time, their identification results are very dispersive and far away from their true values. That is, there exists the problem of parameter identifiability. Therefore, the feasible scheme is to fix the stator resistance as the value on the nameplate and make it not participate in identification.     The average trajectory sensitivity of each candidate identification parameter for 1-ohm transition resistance is shown in Table 4. Clearly, the average trajectory sensitivity of stator reactance, inertia time constant, proportional and integral amplification constants of MSC control loop are still 0. Moreover, the average trajectory sensitivity of R s , ψ PM , K p2 , K i2 , K p3 and K i3 decrease dramatically compared to the 8-ohm transition resistance scenario. This is because the current limit unit of GSC control loop is activated when voltage magnitude drops too much.
I g , I gmin and I gmax are the currents of the limiting unit in the controller of GSC, its minimum and maximum respectively. I m , I mmin and I mmax are the currents of the limiting unit in the controller of MSC, its minimum and maximum respectively. In this paper, only ψ PM , K p2 , K i2 , K p3 , K i3 are identified using PMU data.

Parameters Identification Algorithm
In [27], a novel GLPSO hybrid algorithm with cascade structure using a single exemplary is proposed. However, in GLPSO hybrid algorithm, the inertia weight coefficient ω, acceleration Energies 2020, 13, 4607 8 of 17 coefficient c and crossover and mutation probabilities P c and P m are all constant. In view of this deficiency, the improved GA (IGA) and improved PSO (IPSO), which are derived from standard GA and PSO, are combined with GLPSO to further accelerate the convergence rate in this paper. The procedures of the improved GLPSO hybrid algorithm are shown in Figure 4. It includes the following 5 parts.

Configurations of Case
As shown in Figure 5, the detailed model of WECC benchmark system is used to test the capability of the proposed method. The physical and control system parameters of each PMSG wind turbine and the step up transformer are shown in Table 1 and 2 respectively. The system base capacity is chosen to be B e q = 30 S S = MVA when all PMSG wind turbines are online. A three-phase shortcircuit fault with 8 Ω transition resistance lasting for 4 cycles is set at the middle L1. The simulation time last for 0.2 s. Population size and maximum iteration of GLPSO is set to 50 and 15 respectively. The parameters identification interval is shown in Table 5. (1) The exemplar, velocity and location for each particle as well as population size, maximimal iterations and other parameters are initialized. The fitness value of each exemplar and particle is computed and the best particle is saved.
(2) The optimization parameters of IGA are updated according to (19), (20) [18]: where P i c is crossover probability at the i th iteration. P c0 is the initial crossover probability. N is maximal number of iterations.
The recommended value range of mutation probability P m is 0.001-0.1. The P m is adjusted according to (20) [14].
where τ is a binary denoting whether reseting P i m to the initial value P 0 m . (3) Each exemplar is updated using IGA operation with the genetic materials gbest and pbest. Then, the objective function is computedfor each exemplar. After that, the optimal value and corresponding exemplar are saved.
(4) After IGA operation, the inertia weight coefficient ω, acceleration coefficient c are updated as follows.
where ω max and ω min are the upper and lower bound of ω. c max and c min are the upper and lower bound of c.
The IPSO method is adopted to update all the particles using the new exemplar. The gbest and pbest as well as objective function of every particle are computed. Then the optimal fitness value and corresponding particle is saved. If the end conditions are satisfied, the program is terminated. Otherwise, go to Step 2. The new iteration cycle starts.

Configurations of Case
As shown in Figure 5, the detailed model of WECC benchmark system is used to test the capability of the proposed method. The physical and control system parameters of each PMSG wind turbine and the step up transformer are shown in Tables 1 and 2 respectively. The system base capacity is chosen to be S B = S eq = 30 MVA when all PMSG wind turbines are online. A three-phase short-circuit fault with 8 Ω transition resistance lasting for 4 cycles is set at the middle L1. The simulation time last for 0.2 s. Population size and maximum iteration of GLPSO is set to 50 and 15 respectively. The parameters identification interval is shown in Table 5.

Convergence Rate and Searching Capability of Improved GLPSO
Set wind speed according to Table 5 in [22]. The value of objective function for 3 methods are shown in Figure 6. Clearly, the global searching capability of improved GLPSO is much higher than canonical PSO and GA. The reason is that the capability and convergence rate of GLPSO is much higher than PSO and GA, as well as the optimization coefficients are updated adaptively.

Convergence Rate and Searching Capability of Improved GLPSO
Set wind speed according to Table 5 in [22]. The value of objective function for 3 methods are shown in Figure 6. Clearly, the global searching capability of improved GLPSO is much higher than canonical PSO and GA. The reason is that the capability and convergence rate of GLPSO is much higher than PSO and GA, as well as the optimization coefficients are updated adaptively.

Convergence Rate and Searching Capability of Improved GLPSO
Set wind speed according to Table 5 in [22]. The value of objective function for 3 methods are shown in Figure 6. Clearly, the global searching capability of improved GLPSO is much higher than canonical PSO and GA. The reason is that the capability and convergence rate of GLPSO is much higher than PSO and GA, as well as the optimization coefficients are updated adaptively.

Dispersion of Identified Parameters
The results of identified parameters under various wind speeds and wake effect scenarios are shown in Table 6. Clearly, the dispersion of p2 K , p3 K is low. Other parameters have a certain degree of dispersion. Assume that the wind turbines 5,9,13, and 17 in Figure 5 are offline. The system base

Dispersion of Identified Parameters
The results of identified parameters under various wind speeds and wake effect scenarios are shown in Table 6. Clearly, the dispersion of K p2 , K p3 is low. Other parameters have a certain degree of dispersion. Assume that the wind turbines 5, 9, 13, and 17 in Figure 5 are offline. The system base capacity is chosen to be S B = S eq = 30 × 0.8 MVA. As shown in Table 7, the dispersion of K p2 and K p3 is still low. As shown in Table 5 in [22], the wake effect is considered by setting wind speeds of each PMSG wind turbine different from each other. As shown in the last row in Tables 6 and 7, the precision of parameter identification results for wake effect is slightly lower than other scenarios due to the large difference of operating states of each PMSG wind turbine. However, it is still satisfactory.

Verification of Descriptive Capability
The fitting curves when the wind speed of each turbine is set to 8 m/s, 9 m/s, 10 m/s, and 11 m/s are shown in Figure 7. Clearly, the output active and reactive power of the equivalent model overlaps those of wind farm respectively. Therefore, the equivalent model can describe dynamic characteristics of the wind farm after disturbance with very high precision.

Verification of Generalization Capability
Equivalent model identified with 9 m/s wind speed is utilized to describe wind farm with 7 m/s, 12 m/s wind speed, under wake effect and 20% wind turbines are off line respectively. As shown in Figure 8, the deviation of output active and reactive power curves of the equivalent model from those of wind farm is very small, although the wind speeds differ much. As shown in Table 8, the goodness of fit [18] are very close to 100%. It concludes that the generalization capability of the proposed equivalent model and parameters identification strategy is high enough. is still low. As shown in Table 5 in [22], the wake effect is considered by setting wind speeds of each PMSG wind turbine different from each other. As shown in the last row in Tables 6 and 7, the precision of parameter identification results for wake effect is slightly lower than other scenarios due to the large difference of operating states of each PMSG wind turbine. However, it is still satisfactory.

Verification of Descriptive Capability
The fitting curves when the wind speed of each turbine is set to 8 m/s, 9 m/s, 10 m/s, and 11 m/s are shown in Figure 7. Clearly, the output active and reactive power of the equivalent model overlaps those of wind farm respectively. Therefore, the equivalent model can describe dynamic characteristics of the wind farm after disturbance with very high precision.

Verification of Generalization Capability
Equivalent model identified with 9 m/s wind speed is utilized to describe wind farm with 7 m/s, 12 m/s wind speed, under wake effect and 20% wind turbines are off line respectively. As shown in Figure 8, the deviation of output active and reactive power curves of the equivalent model from those of wind farm is very small, although the wind speeds differ much. As shown in Table 8, the goodness of fit [18] are very close to 100%. It concludes that the generalization capability of the proposed equivalent model and parameters identification strategy is high enough.

Verification of Generalization Capability
Equivalent model identified with 9 m/s wind speed is utilized to describe wind farm with 7 m/s, 12 m/s wind speed, under wake effect and 20% wind turbines are off line respectively. As shown in Figure 8, the deviation of output active and reactive power curves of the equivalent model from those of wind farm is very small, although the wind speeds differ much. As shown in Table 8, the goodness of fit [18] are very close to 100%. It concludes that the generalization capability of the proposed equivalent model and parameters identification strategy is high enough.
(a)   The results of parameters identification when wind speed is unknown are shown in Table 9. Clearly, they are consistent with those in Tables 6 and 7. It is worth pointing out that the identification results of ψ PM in Tables 6, 7 and 9 deviate a lot from the true value 1.1884. When the three-phase short-circuit transition resistance is set to 1 ohm, we have carried out the program of parameters identification many times. It is found that although other parameters cannot be identified precisely, K p3 and ψ PM can be identified precisely. For all the programs, the mismatch between the identification results for ψ PM and the true value 1.1884 are less than 5%. Further, the transient process time of power system is very short, generally within 10 s. The dynamic characteristics of PMSG wind farm depend on the regulating speed of GSC converter, which is generally in millisecond. As shown in Figures 7 and 8, the transient process of PMSG wind farm is very short, far less than 0.2 s. In such a short period of time, the wind speed can be regarded as a constant value. The proposed method has good robustness and adaptability for the randomness of wind speed.

Conclusions
We show that the sensitivity of inertia time constant T j , stator reactance L d , L q and proportional and integral amplification constant of current controller of MSC K p1 , K i1 are very small after trajectory sensitivity analysis. Moreover, the phase of trajectory sensitivity for stator resistance R s is opposite to permanent magnet flux linkage ψ PM . Simulation results using the WECC benchmark system indicate that the capability of the proposed improved GLPSO hybrid algorithm are much higher canonical PSO and GA. Further, the approach proposed is preferable to those in [20,24]. Because in [20] the parameters of PMSG wind farm are identified with load parameters. Since the sensitivity for parameters of load are much higher than parameters of a PMSG wind turbine, it is very difficult to identify the parameters of a PMSG wind turbine precisely if they are identified at the same time. Furthermore, the permanent magnet flux linkage of PMSG is not identified in [20], which has a significant impact on power system dynamics. Moreover, since the approach in [24] belongs to the aggregation method, detailed information of the whole wind farm is necessary. Furthermore, the method in [24] cannot be applicable after long-term operation of wind farm since the parameters of each PMSG wind turbine may deviate greatly from the values on the nameplate. Therefore, the proposed method is preferable to the on in [24] as well. and Development Plan "Important Scientific Instruments and Equipment Development"] grant number [2016YFF0102200].

Conflicts of Interest:
The authors declare no conflict of interest.

Trajectory Sensitivity Analysis
In this paper, the analysis method of trajectory sensitivity and average trajectory sensitivity in reference [13] is adopted. The dynamics of power system is given by: dx dt = f(x, y, u, θ) y = g(x, u, θ) where state variables, stator and rotor flux linkage and rotor speed for example are denoted by x. The input variables, wind speed, voltage of common interconnection point for example are denoted by u. Vector y represents output variables, such as output active power, reactive power, etc. Vector θ represents parameters, such as equivalent distribution network impedance, stator and rotor impedance, excitation reactance, etc. The trajectory sensitivity of parameters reflects the change degree of output variables when the parameters change. In a certain operating state, the trajectory sensitivity of the output y j (t) relative to the parameter θ i is defined as S y j /θ i (t) = ∂y j (t)/y j (t) y j (θ r ,θ i0 +∆θ i ,t)−y j0 (θ r ,θ i0 ,t) y j0 (θ r ,θ i0 ,t) where y j0 is the output when ith parameter value is set to θ i0 . ∆θ i is the increment on θ i . θ r are the parameters in the model except θ i . The average trajectory sensitivity is defined as where t 1 and t 2 are the time before and after the disturbance, which can be reasonably selected according to the output curve samples data.