Optimization to POD parameters of DFIGs based on the 2nd order eigenvalue sensitivity of power systems

Here, the analytical model of the 2nd order eigenvalue sensitivity is improved, and applied to optimize the parameters of the power oscillation dampers of the doubly-fed induction generators. To solve the 2nd order eigenvalue sensitivity, the eigenvector sensitivity is required whose difﬁculty lies in insufﬁcient constraints, for which two constraints about the magnitude and the angle of eigenvector elements are newly introduced, and the normalization conditions are derived to solve the eigenvector sensitivity. The 2nd order eigenvalue sensitivity is applied to analyse the change of the low-frequency oscillation modes with the varying of parameters in the power systems with the doubly-fed induction generators. To improve the low-frequency oscillation modes, the optimization model of the power oscillation damper parameters is newly formulated based on the 1st and the 2nd order eigenvalue sensitivities and solved by the interior point method. The numerical results are provided to verify the effectiveness of the constraints proposed; the advantages of the 2nd order eigenvalue sensitivity in the eigen-analysis and parameters optimization, and the robustness against wind ﬂuctuation.


INTRODUCTION
With the integration of the doubly-fed induction generators (DFIGs), the small signal stability of the power systems may be affected because of the volatility of the wind power and the reduction of the system inertia. The low frequency oscillation (LFO) modes may be affected by the interactions between the DFIGs and synchronous generators (SGs) [1][2], and the oscillation modes related to the transfer shaft of the DFIGs may be introduced [3]. In [4][5][6][7], the impacts of the wind power penetration and the locations of the DFIGs on the small signal stability of power systems are studied based on the eigen-analysis and time domain simulation. The LFO modes of power systems may be damped by the power system stabilizer (PSS) installed at the SGs, or recently, the power oscillation damper (POD) at the DFIGs [8,9]. With the signals from the transmission lines and the SGs inputted to the POD [10,11], the LFO modes may be damped by the POD whose control effects are affected by the parameters setting. The heuristic algorithms, e.g. the whale optimisation algorithm, imperialist competitive algorithms and particle swarm optimization (PSO) algorithm [12][13][14]  while the performance and the convergence characteristics may be limited in the large power systems. In [15], based on the transfer functions of the open loops in the power systems, the phase compensation required for the POD design is obtained, while the integration of the DFIGs may increase the difficulty in deriving the transfer functions. Compared with the methods above, the eigenvalue sensitivity models have the advantages of the high computational efficiency and easy implementation. With the gradient information with respect to parameters provided by the eigenvalue sensitivity, the PSS parameters may be optimized in [16,17], which is also suitable for the POD design. In [18], the effectiveness of the POD installed at the DFIGs under different wind speeds is studied based on the eigenvalue sensitivity with respect to the wind penetration.
The 1 st order eigenvalue sensitivity only shows the linear function between the oscillation modes and the parameters, whose accuracy is limited. The truncation errors of the eigenvalue sensitivity models may be reduced by applying the successive linear programming [19], whose computational efficiency is affected by the multi-step predictions, or introducing the 2nd order eigenvalue sensitivity, while the eigenvector sensitivity required is difficult to solve. To solve IET Gener. Transm. Distrib. 2021;15:1123-1135.
wileyonlinelibrary.com/iet-gtd the problem, the method based on the 2nd order perturbation theory and the method with the normalization conditions are proposed in existing researches. Though the derivatives of the high-dimensional state matrix are avoided by the former method, all the eigenvalues and eigenvectors are required to derive the analytical expression of the 2nd order eigenvalue sensitivity in [21][22][23], which is complicated. In [24][25][26], with the normalization conditions, the eigenvector sensitivity is solved, based on which, the 2nd order eigenvalue sensitivity may be obtained. Since only the eigenvectors corresponding to the eigenvalues of interest are involved, the method with the normalization conditions is easy to apply. However, there are few literatures involving the causes of the insufficient constraints in the eigenvalue sensitivity model, which may the basis for proposing the normalization conditions. And the normalization conditions may affect the accuracy of the eigenvector sensitivity.
Based on the normalization conditions about the magnitude and the angle of eigenvector elements, the analytical expression of the 2nd order eigenvalue sensitivity is derived in this paper. To damp the LFO modes, the POD is installed at the DFIGs and tuned based on the 1st and the 2nd order eigenvalue sensitivities. The main contributions of this paper are summarised as follow, 1. Two constraints about the magnitude and the angle of eigenvector elements are newly proposed, based on which, the normalization conditions are derived to solve the eigenvector sensitivity and the 2nd order eigenvalue sensitivity. The case study shows the effectiveness of the proposed normalization conditions. 2. The change of LFO modes in the wind power systems is analysed based on the eigenvalue sensitivity. The case study shows that, by introducing the 2nd order eigenvalue sensitivity, the accuracy of the sensitivity model is improved. 3. To damp the LFO modes, the optimization problem of the POD parameters is modelled based on the 1st and the 2nd order eigenvalue sensitivities and solved with the interior point method. The case study shows that the POD may improve the damping of the power system effectively. It is shown in this work that the optimization results of the POD parameters may be more reliable with the introduction of the 2nd order eigenvalue sensitivity.
This paper is organised as follow: Section 2 proposes the new normalization conditions about the magnitude and the angle of eigenvector elements to solve the eigenvector sensitivity, based on which, the 2nd order eigenvalue sensitivity may be obtained. Section 3 explains the modelling of the DFIGs and POD, and formulates the state equations of the wind power systems. Section 4 formulates the optimization model of the POD based on the 1st and 2nd order eigenvalue sensitivities. In Section 5, the case study based on the New England 39-bus test system with 2 DFIGs shows the effectiveness of the proposed method, the advantages of the 2nd order eigenvalue sensitivity and the robustness of the POD. The conclusions of this paper are given in Section 6.

The 1st order eigenvalue sensitivity
With the linearization to differential/algebraic equations, the power systems are given as, where A sys is the state matrix; x is the state variable; s is the Laplace operator and Δ is the increment. And the eigenvalue λ is obtained by, The 1st order eigenvalue sensitivity with respect to the configurational or control parameter α is expressed by, where v, w are the right and left eigenvectors corresponding to λ 1 , respectively.

2nd order eigenvalue sensitivity
The 1st order eigenvalue sensitivity shows the linear function between the oscillation modes and parameters. However, the truncation error with the 1st order eigenvalue sensitivity may yield error in analysing and damping the oscillation modes. Hence, the 2nd order eigenvalue sensitivity is applied here.
With the parameter perturbation, we may find the 1st and the 2nd order eigenvalue sensitivity, but the calculation effort is huge. An analytical expression of the eigenvalue sensitivities is more convenient for actual application. The 2nd order eigenvalue sensitivity with respect to the parameters α and μ is defined as, where E is the identity matrix. It is seen that to find the 2nd order eigenvalue sensitivity, the eigenvector sensitivity is needed, which may be obtained with the derivatives of Equation (2), But Equation (5) can not be solved because (A sys -λE) is singular, which means some constraints are ignored. To solve the eigenvector sensitivity, a normalization condition based on the magnitude of eigenvectors is proposed in [24], where the superscript * is the conjugate transpose. However, the Equations (5) and (6) still can not be solved because the real/imaginary parts of the eigenvector sensitivity are required, which means more normalization conditions are needed.
In [25], since the effectiveness of Equation (6) is independent of the value of Im[v*(∂v /∂α)], the normalization condition is extended by assuming Im[v*(∂v/∂α)] equal to 0, However, the effect of the value of Im[v*(∂v /∂α)] on the constraints of Equation (5) is not considered, which may affect the results of the eigenvector sensitivity.
In this paper, the causes of the insufficient constraints are analysed, which may be the basis for the normalization conditions. From Equation (2), it is found that there are multiple eigenvectors corresponding to one eigenvalue in the eigenanalysis models, which means some constraints about the eigenvectors are not involved and may cause the problem in solving the eigenvector sensitivity. Thus, the additional constraints to determine the eigenvector are required, based on which, the normalization conditions may be improved.
First, with the rearranges of the real/imaginary parts of Equation (5), the number of normalization conditions required may be analysed based on rank of the coefficient matrix, Assuming A sys is a n-dimensional matrix, the rank of the 2ndimensional matrix A sum is 2n-2, which means 2 normalization conditions are needed to solve the eigenvector sensitivity.
Then, to determine the eigenvector, the constraints of eigenvectors corresponding to λ are studied based on the eigenvector elements. Assuming v i is the ith element of v, Equation (2) is expressed as, The following equation is obtained by multiplying a complex number k∠θ by both sides of (9).
where k is not equal to 0.

FIGURE 1 Impacts k and θ on v i
It is found that the eigenvectors change with the varying of k and θ, while all the eigenvectors meet the constraints of Equation (9). Thus, to determine the eigenvector, the effects of k and θ on the eigenvectors need to be analyzed, which may be studied based on the elements. Figure 1 shows v i in the complex plane corresponding to different values of k and θ, in which the magnitude and angle of v i are controlled by k and θ, respectively.
It should be noted that the effects of k and θ on each element of v are same, e.g. the multiple of the magnitude changing and the angle of rotation. Thus, by taking the constraints of v i into Equation (2), v may be determined when the unique value of v i is obtained. And the eigenvector sensitivity may be solved with the normalization conditions based on the constraints of v i . Considering the constraints of the magnitude and angle, the normalization conditions are improved as follow,

the magnitude of v i
By controlling k, the magnitude of v i may be kept unchanged, where c is a constant.
With the derivatives of Equation (11) with respect to α, the normalization condition is derived as, Since the ratio between the magnitudes of v and v i remains unchanged with the varying of k, the constraint of Equation (11) may be replaced by maintaining the magnitude of v to be equal to 1. Thus, the constraints of Equation (12) is same as Equation (6), which may explain the normalization condition based on the eigenvector elements.

the angle of v i
By multiplying a complex number with the angle of θ, the angle of v i is maintained to be equal to 0, which may be expressed by the constraints of the imaginary part of v i , With the derivatives of Equation (13) with respect to α, the 2nd normalization condition is derived as, By substituting Equations (12) and (14) into Equation (8), the rank of the coefficient matrix is 2n, which means the eigenvector sensitivity and the 2nd order eigenvalue sensitivity may be solved. Since the normalization conditions are proposed based on the constraints of Equations (11) and (13), the right/left eigenvectors should meet the constraints provided by them.

Modelling of DFIG
In the study of the LFO modes of wind power systems, the DFIGs and the power systems may be considered as two open loop subsystems, respectively. The active/reactive power outputs of the DFIGs (P DFIG , Q DFIG ) are provided from the DFIGs for the power systems while the amplitude/angle of voltage of the DFIGs (U DFIG , θ DFIG ) are controlled by power systems. Figure 2(a) shows the configuration of the DFIGs [27], which is consisted with the wind turbine (WT), the pitch angle controller (PAC, the 2nd order), the transfer shaft (the 3rd order), the induction generator (IG, the 4th order), the rotorside converter (RSC, the 4th order) and its controller, the GSC (the 4th order) and its controller, the transformer, the dc capacitor (the 1st order) and the POD (the 3rd order). Figure 2(b) shows the control strategies of the DFIGs. The energy of the wind is converted into the mechanical energy of the DFIGs by the WT, whose captured power (P t ) may be expressed as, where P wloss is the lost power of the wind, C P is the coefficient of the power utilization which is affected by the rotor speeds of the WT (ω t) and the pitch angle (β). To maintain the optimal rotor speeds of the WT, β is controlled by the PAC, which also regulates P t and the torque of the WT (T t ).
The WT is connected to the rotor of the IG though the transfer shaft which increases ω t to a higher rotor speeds of the IG (ω r). The transfer shaft is described by the two-lumped masses model, whose variables may be obtained by, To achieve the decoupling control of active/reactive powers of the DFIGs, the vector control technique is applied. The daxis is oriented based on the stator voltage of the IG and q-axis is assumed to be 90 • ahead of the d-axis. The RSC controls the active/reactive powers of the stator of the IG (P s and Q s ), which are regulated by the direct/quadrature currents of the rotor of the IG (I rd and I rq ), respectively. The GSC maintains the voltage of the dc capacitor (U dc ) and controls the reactive power from the convectors to the grid (Q g ), which are regulated by the direct/quadrature currents of the GSC (I gd and I gq ), respectively. U dc is determined by the power imbalance of the converters. Under the stator-voltage-oriented frame, the quadrature stator voltage of the IG is equal to 0, and the other voltages of the stator/rotor of the IG are given by the electromagnetic transient of the IG, where I and U are the current and voltage, respectively; L ss , L rr are the inductances of the stator and rotor of the IG, respectively; L m is the magnetizing reactance of the IG; R s and R r are the resistances of the stator and rotor of the IG respectively and the subscripts sd, sq, rd, rq denote the direct/quadrature variables of the stator and rotor of the IG, respectively.

Modelling of power oscillation damper
To damp the LFO modes, the POD may be installed at the DFIGs, which shares the similar structure with the PSS at the SGs. In this paper, the POS is the POD is modelled based on the PSS1A which is easily to tune. In existing researches, the POD is usually installed at the control loops of convectors in the DFIGs, for example, the active/reactive power control loops of the RSC and the dc voltage/reactive power control loops of the GSC. Though the convectors of the DFIGs do not yield the LFO modes for the fast response, the LFO modes are still affected by the convectors since the power outputs of the DFIGs are controlled by those. Considering the active power of the DFIGs is controlled by the active power control loop of the RSC directly, the POD installed at this locations may be more effective in damping the LFO modes than those at the other locations in the DFIG, for the LFO modes may be more sensitivity to the active power than the reactive power and the dc voltages.
As shown in Figure 3, with the rotor speeds of the SGs strongly associated with the LFO modes adopted as the input signal (ω in), the output signal of the POD (P POD ) is provided for the RSC controller. K POD sets the gain of the POD, which determines the damping effect on the LFO modes; T h is the time constant of the washout filter which ensures the output of

State equations of power systems with DFIGs
Combining the models of the DFIGs, the POD, the SGs (the 4th order) and the grid, the state equations of wind power systems are given by (18), where A, B, C and D are the coefficient matrices, y denote the algebraic variables.
There are three frames in the wind power systems, e.g. the x/y frame in the grid, the d/q frame in the DFIGs and the frame oriented based on the power angle in the SGs. Since the angle of the d-axis in the DFIGs leading the x-axis in the grid is equal to the angle of the stator voltage of the IG (θ IG,s ), the interface equations of the current injected to the grid from the DFIGs (I DFIG ) are given by Equation (19), where the subscribes x, y, d and q denote the x/y and d/q frames. The interface equations of the SGs may be obtained by replacing θ IG,s by the power angle of the SGs. [ To simplify the linear models of power systems, the constant power loads in the power systems, which are usually considered as the P∖Q buses in the power flow program, are replaced by the constant impedance loads. It helps to reduce the number of algebraic variables in the state equations while only the DFIGs and SGs may provide the currents injected to the grid.

Robustness of POD
Based on the modelling of power systems with DFIGs, the close loop with unknown disturbance is shown in Figure 4. On the basis of the small gain theorem and robust stabilisation, the controller robustness may be quantified by the infinity norm of M, whose marginal value is 2 [13].

LFO modes of wind power systems
The LFO modes are strongly associated with the SGs, whose electromechanical loop participation ratio (ρ SG ) is much greater than 1 and the frequency (f) is within 0.1 to 2.5 Hz, where p is the participation factor and the subscript SG,MT denotes the mechanical transient of the SGs. In this paper, the LFO modes with the weakest damping ratio may be defined as the critical mode of wind power systems, λ c = σ c +jω c .

Optimization model of POD parameters
With the rotor speeds of the SGs strongly associated with critical mode applied as the input signal, the LFO modes may be damped by the POD, and the control parameters of the POD, e.g., K POD , T d1 and T d2 , may be optimized [9]. To find the optimal solution within the ranges of the POD parameters, the damping effects of the POD are formulated as the objective function given by Equation (23), where ξ c is the damping ratio of the critical mode.
Since the Hopf bifurcation may happen when the oscillation modes are close to the imaginary axis, the limitation of σ c is considered as the inequality constraint given, where σ max is the marginal value equal to −0.15.
The limitations of the POD parameters are given as, (25)

Parameters optimization with interior point method
The optimization model may be solved by the interior point method. The inequality constraints may be transformed into the equality constraints with the introduction of the slack variables (u, l), The object function may be transformed into the Lagrange function, where z T and o T are the Lagrange multipliers.
According to the Kuhn-Tucker theorem, the optimal solution is obtained when the derivatives of the Lagrange function with respect to the variables and the Lagrange multipliers are equal to 0.
The equations may be solved by the Newton-Raphson method, whose difficulty lies in the Jacobi matrix and the Hessian matrix. Based on the 1st and the 2nd order eigenvalue sensitivities, λ c may be formulated by the POD parameters, where the subscript 0 denotes the initial value. Thus, the 1st and the 2nd order derivatives of ξ c and σ c with respect to the POD parameters may be obtained which are required by the Jacobi matrix and the Hessian matrix.
(28) Figure 5 shows the flowchart of the optimization of the POD parameters based on the 1st and the 2nd order eigenvalue sensitivities.

5
RESULTS AND DISCUSSIONS Figure 6 shows the New England 39-bus test system with 2 DFIGs, which is widely applied to study the LFO modes. The SGs are represented by the fourth-order model, which are equipped with the automatic voltage regulator. The modelling of the DFIGs is given in Section 3, whose parameters may be found in [3]. There are three areas in the test system which are connected by four interconnectors, and two wind farms are aggregated to DFIG 1 and DFIG 2 respectively. The area 1 and area 2 are connected by the interconnector between bus 1 and bus 39 and the interconnector between bus 3 and bus 4. The area 3 is connected to area 1 by the interconnector between bus 16 and bus 17, which is connected to area 2 by the interconnector between bus 14 and bus 15. A wind farm with 90 DFIGs are aggregated to the DFIG 1 located at the bus 9, and the other wind farm with 60 DFIGs are aggregated to the DFIG 2 located at the bus 14. υ w1 is equal to 11 m/s and υ w2 is equal to 9 m/s.

Impact of DFIGs integrated on LFO modes
The LFO modes in the power systems without and with the DFIGs are given in Tables 1 and 2 respectively, which are defined in Section 4.1. With the DFIGs integrated, the margin of the damping ratio is reduced from 0.0188 to 0.0180. Among 9 LFO modes existing in power systems, the LFO mode λ 35,36 with the weakest damping ratio may be considered as the critical mode (in bold). More detailed analysis based on the participation factors shows that λ 35,36 is strongly associated with the SG 1 and SG 10 . Table 3 shows the 1st order eigenvalue sensitivities of λ 35,36 with respect to the parameters of the PODs installed at the DFIG 1 and DFIG 2 respectively. It is seen that λ 35,36 is more

FIGURE 7
Effects of υ w1 on ζ 35,36 with and without POD sensitive to the POD installed at the DFIG 1 than that installed at the DFIG 2 , which may be affected by the locations, capacities and parameters of the DFIGs. To avoid the interactions between the PODs at different wind farms, the POD is only installed at the DFIG 1 in this paper. The POD parameters in the initial point are given as follow: K POD = 10, T d1 = 0.5s, T d2 = 0.13s.

Effects of wind speeds on critical mode
With υ w1 varying from 8 to 14 m/s, Figure 7 shows the curves of the damping ratio ζ 35,36 with and without the POD installed at the DFIG 1 . It is found that, under various wind speeds, ζ 35,36 is improved with the POD installed, which means the POD is effective in damping the critical mode. Since the active power of the DFIGs is controlled to damp the critical mode by the POD, ζ 35,36 become more sensitive to υ w1 with the application of the POD, which affects the captured power of the DFIGs.

5.2
Validation of the 2nd order eigenvalue sensitivity

Accuracy of eigenvector sensitivity
As shown in Figure 8, the eigenvector sensitivity with respect to T d2 is obtained based on the proposed method and the parameter perturbation, respectively. Due to the space limitations, only the sensitivities of 20 elements of v with respect to T d2 are given. The eigenvector sensitivities obtained based on the both methods are equal, which proves the effectiveness of the proposed method. The initial values of the POD parameters are given in Section 5.1. Figure 9 shows the sensitivities of 20 elements of v 0 when Im[(v 0 )*(∂v/∂T d2 )] is equal to 0 and 0.5, respectively. It is seen that the result of ∂v/∂T d2 changes with the varying of Im[(v)*(∂v/∂T d2 )], which means the assumption about Im[(v)*(∂v/∂T d2 )] may affect the eigenvector sensitivity.

Accuracy of the 2nd eigenvalue sensitivity
The comparison between curves of λ 35 and ζ 35,36 computed based on the 1st order eigenvalue sensitivity only (blue line) and those obtained with the introduction of the 2nd order eigenvalue sensitivity (red line) is shown in Figure 10, with T d2 varying from 0.10 to 0.16 s. The actual value (black line) is obtained based on the continuously changing parameter method by increasing T d2 from 0.10 to 0.16 s in the steps of 0.00006 s. The initial values of the POD parameters are given in Section 5.1.
Near the initial point, the curves of λ 35 and ξ 35,36 are tangent to the actual value, which means the eigenvalue sensitivities are effective. It may be seen that, though the errors of λ 35 and ξ 35,36 are based on the eigenvalue sensitivities are inevitable, the errors are reduced obviously by applying the 2nd order eigenvalue sensitivity, since the nonlinear structure of the power systems is considered.
Apart from it, there is a turning of the actual value of λ 35 , which corresponds to the change of the trend of ξ 35,36 . It may be seen that the result of ξ 35,36 based on the 1st order eigenvalue sensitivity shows the wrong trend with T d2 increasing from 0.1 to about 0.12 s, while the result of ξ 35,36 may follow the change of the actual value with the introduction of the 2nd order eigenvalue sensitivity. This mistake is due to the linear model of the 1st order eigenvalue sensitivity, which may be unignorable in the study of oscillation modes and affect the parameters optimization.
To quantitatively study the accuracy and the computational performance, Table 4 shows the maximum errors and computation times corresponding to the methods above. The result obtained by the continuously changing parameter method may be considered as the actual value, which is more accurate than those based on the eigenvalue sensitivities. However, the continuously changing parameter method requires 19.056 s, whose computational efficiency is low. By introducing the 2nd order eigenvalue sensitivity, the maximum errors of λ 35 and ξ 35,36 are reduced by 26.3% and 66.1% respectively, and the computation time only increases from 0.799 to 0.863 s. Thus, the eigenvalue sensitivity model is more computationally economic than the continuously changing parameters method, whose accuracy is improved obviously by introducing the 2nd order eigenvalue sensitivity.

Validation to control effect of POD with time-domain analysis
Based on the method explained in Section2, the 1st and the 2nd order eigenvalue sensitivities of λ 35 with respect to the POD parameters are given in Table 5.
The parameters of the POD a are optimized with the linear programming based on the 1st order eigenvalue sensitivity only, while the POD b is tuned with the introduction of the 2nd order eigenvalue sensitivity. For comparison, the POD c is tuned based on the PSO algorithm. The parameters and computation times corresponding to the PODs are given in Table 6. The ranges of the POD parameters are given as follow: K POD [5,15], T d1 [0.4s, 0.6s]. T d2 [0.05s, 0.20s].
Compared with the POD a tuned based on the 1st order eigenvalue sensitivity, the computation time corresponding to the POD b increases by 0.058 s with the introduction of the 2 nd order eigenvalue sensitivity. However, the POD c tuned based on the PSO algorithm requires more than 30× of the computation time compared with the PODs tuned based on the eigenvalue sensitivities. Thus, the optimization with eigenvalue sensitivity has an obvious advantage in computation efficiency than the heuristic methods. Table 7 shows the critical modes of power systems without the POD, equipped with the POD a , POD b and POD c respectively. The POD a , POD b and POD c are effective in damping the critical modes. The POD c damps the critical modes most effectively since the corresponding ζ 35,36 is the largest and λ 35,36 is the left-most, which means that the errors of the eigenvalue sensitivities may affect the parameters optimization. However, with the introduction of the 2nd order eigenvalue sensitivity, the damping effect of the POD b is close to the POD c , which is improved obviously compared with the POD a tuned based on the 1st order eigenvalue sensitivity only.
Since the critical mode λ 35,36 is strongly associated with the SG 1 and SG 10 , the damping ratio may be extracted from the rotor dynamics of the SGs. As shown in Figure 11, the oscillations of power angle difference (δ 1,10 ) and rotor speeds of the SG 1 (ω 1 ) damped by POD a (blue line), POD b (red line) and POD c (black line) are compared to study the damping effect of the PODs. A three-phase short-circuit fault occurs at bus 3 at 0.3 s, which is cleared at 0.4 s.
It is seen that the oscillations of δ 1,10 and ω 1 corresponding to the POD a is the most obvious among the PODs, whose amplitude is the largest. And the oscillations damped by the POD b and POD c are close in about 15 s, which correspond to the critical modes in Table 7. Thus, with the introduction of the 2 nd order eigenvalue sensitivity, the optimization based on the eigenvalue sensitivity model becomes more reliable, which has a higher computation efficiency than the PSO algorithm.

Robustness of POD against wind fluctuation
Based on the models in Section 3.4, the robustness indexes of the POD a , POD b and POD c are shown in Table 8. Compared with the marginal value equal to 2, the robustness of the PODs may be acceptable. More analysis based on time domain simulation is given to study the impacts of the wind speeds on the POD. As shown in Figure 12, the damping effects of the POD are studied based on the oscillations of δ 1,10 when υ w1 is equal to 8 and 14 m/s, respectively. The POD parameters are same as the POD b which is tuned based on the 1st order and the 2nd order eigenvalue sensitivities when υ w1 is equal to 11 m/s. It is seen that the POD damps the oscillations of δ 1,10 effectively under different wind speeds and the damping effects under the higher wind speed may be more obvious.

CONCLUSIONS
The constraints about the magnitude and the angle of eigenvector elements are proposed to derive the analytical expression of the 2nd order eigenvalue sensitivity. To study the LFO modes, the state equations of power systems with the DFIGs are formulated, and the critical modes are found. To improve the small signal stability of wind power systems, the DFIGs are equipped with the POD whose parameters are optimized based on the 1st and the 2nd order eigenvalue sensitivities. The numerical results show that, 1. The eigenvector sensitivity obtained based on the proposed method is close to the actual value, which proves the effectiveness of the normalization conditions. 2. By applying the 2nd order eigenvalue sensitivity, the errors of the predictions based on the eigenvalue sensitivity are reduced, which means the accuracy of the eigenvalue sensitivity is improved. 3. With the introduction of the 2nd order eigenvalue sensitivity, the optimization models become more reliable, which are more computationally economic than that based on the PSO method. 4. The POD may damp the LFO modes under different wind speeds.
In the following studies, the interaction between the PODs installed in the different DFIGs may be studied.