CN104218864B - A kind of frequency domain method of double-fed fan motor unit rotor side controller parameter identification - Google Patents

Abstract

The invention discloses the frequency domain method of a kind of double-fed fan motor unit rotor side controller parameter identification.First, rotor-side controller is determineddAxle andqDecoupler shaft model and parameter to be identified；Then in the reference signal of rotor-side controller, inject pseudo-random signal；And according to the input pseudo-random signal of synchronous acquisition and the rotor current of outputdAxle orqAxle component, the auto-power spectrum of calculating input signal and input signal and the crosspower spectrum of output signal, obtain the frequency sequence of the transmission function of rotor-side controller Decoupled Model accordingly；It is finally based on the Damped Nonlinear least square optimization identification each parameter of rotor-side controller.The present invention is based on the Decoupled Model identification each parameter of rotor-side controller, and experimental technique is simple；On input reference signal, stacking pseudorandom signal is as excitation, can carry out when not affecting double-fed fan motor unit and normally working, and method of testing is practical.

Description

Frequency domain method for identifying parameters of rotor side controller of doubly-fed wind turbine generator

Technical Field

The invention relates to a frequency domain method for identifying parameters of a rotor side controller of a doubly-fed wind turbine generator, and belongs to the field of power system modeling.

Background

Since the introduction of wind power generation into power systems, modeling problems, particularly dynamic modeling of wind turbines, have been the focus of research. The controller is an important component of the doubly-fed wind turbine generator, and the control mode and the model parameters of the controller have large influence on the dynamic characteristics of the doubly-fed wind turbine generator, so that the controller of the doubly-fed wind turbine generator needs to be accurately modeled.

In the existing literature, when identifying the parameters of the wind generating set controller, the identification is generally performed based on an original controller model using rotor voltage as an output quantity, however, there is coupling between d-axis and q-axis variables in the original model, if the parameters are directly identified based on the model, the number of experiments is large and the operation is complicated, and therefore, a d-axis and q-axis decoupling model of the controller using rotor current as an output quantity is proposed, and the parameters are identified accordingly.

In the existing parameter identification technology, a time domain identification method is mostly adopted. However, the time domain identification method usually requires a large disturbance experiment, and due to the characteristics of the power system, there is a certain risk that a large disturbance excitation is applied to the actual system, and the operation cannot be performed frequently.

Disclosure of Invention

The technical problem to be solved by the invention is to provide a frequency domain method for identifying parameters of a rotor side controller of a doubly-fed wind turbine generator, the current of a generator rotor is used as an output signal, a d-axis decoupling model and a q-axis decoupling model of the rotor side controller are realized, each parameter of the rotor side controller is identified, and the experimental method is simple; the pseudo-random binary signal is applied to the reference signal, the pseudo-random signal is superposed on the input reference signal to serve as excitation, the method can be carried out when the normal work of the double-fed wind turbine generator is not influenced, and the test method is practical.

The invention adopts the following technical scheme for solving the technical problems:

the invention provides a frequency domain method for identifying parameters of a rotor side controller of a doubly-fed wind turbine generator, which comprises the following steps of:

step 1, determining a d-axis and q-axis decoupling model and a parameter to be identified of a rotor side controller of a doubly-fed wind turbine generator, specifically:

the d-axis decoupling model of the rotor-side controller takes the voltage deviation amount as input and the d-axis component I of the rotor-side current_drAs an output, its transfer function H₁(s,θ₁) The expression of (a) is:

$H_1(s,{\mathrm{\θ}}_1)=\frac{I_{\mathrm{dr}}\left(s\right)}{V_{s\_\mathrm{ref}}-V_s}=\frac{\frac{K_{p3}{K}_{p4}}{L_r^{\′}/{\mathrm{\ω}}_B}{s}^2+\frac{K_{p4}{K}_{i3}+K_{p3}{K}_{i4}}{L_r^{\′}/{\mathrm{\ω}}_B}s+\frac{K_{i3}{K}_{i4}}{L_r^{\′}/{\mathrm{\ω}}_B}}{s(s^2+\frac{K_{p4}}{L_r^{\′}/{\mathrm{\ω}}_B}s+\frac{K_{i4}}{L_r^{\′}/{\mathrm{\ω}}_B})}---\left(1\right)$

in the formula, s is a Laplace operator; theta₁For the d-axis decoupling model to be identified, θ₁＝[K_p3,K_i3,K_p4,K_i4,L_r′]^T；I_dr(s) is the d-axis component I_drThe Ralsberg transform of (1); v_{s_ref}Is a stator voltage reference value; v_sIs the stator voltage; omega_BIs a reference angular frequency; k_p3And K_i3Proportional and integral gains, K, of the outer loop for d-axis voltage control, respectively_p4And K_i4Proportional and integral gains of the inner d-axis voltage control loop respectively; l is_r' is rotor transient inductance;

the q-axis decoupling model of the rotor side controller takes the active power deviation value as input and the q-axis component I of the rotor side current_qrAs an output, its transfer function H₂(s,θ₂) The expression of (a) is:

$H_2(s,{\mathrm{\θ}}_2)=\frac{I_{\mathrm{qr}}\left(s\right)}{P_{s\_\mathrm{ref}}-P_s}=\frac{\frac{K_{p2}{K}_{p1}}{L_r^{\′}/{\mathrm{\ω}}_B}{s}^2+\frac{K_{p2}{K}_{i1}+K_{p1}{K}_{i2}}{L_r^{\′}/{\mathrm{\ω}}_B}s+\frac{K_{i2}{K}_{i1}}{L_r^{\′}/{\mathrm{\ω}}_B}}{s(s^2+\frac{K_{p2}}{L_r^{\′}/{\mathrm{\ω}}_B}s+\frac{K_{i2}}{L_r^{\′}/{\mathrm{\ω}}_B})}---\left(2\right)$

in the formula, theta₂For the q-axis decoupling model to be identified, θ₂＝[K_p1,K_i1,K_p2,K_i2,L_r′]^T；I_qr(s) is the q-axis component I_qrThe Ralsberg transform of (1); p_{s_ref}Is an active power reference value; p_sActive power; k_p1And K_i1Respectively controlling the proportion and integral gain of an outer ring for the q-axis active power; k_p2And K_i2Proportional and integral gains of the q-axis active power control inner ring are respectively obtained;

step 2, in the reference value V of the stator voltage_{s_ref}Applied voltage pseudo-random signal V_PRBSAnd synchronously collecting V_PRBSAnd d-axis component I of rotor-side current_dr(ii) a At the active power reference value P_{s_ref}To which an active power pseudo-random signal P is applied_PRBSAnd synchronously collecting P_PRBSAnd q-axis component I of rotor-side current_qr；

Step 3, calculating a stator voltage reference value V_PRBSSelf-power spectrum G of_VV(ω) and V_PRBSAnd I_drCross power spectrum G of_VI(ω) to calculate a frequency series of transfer functions of the d-axis decoupling model of the rotor-side controllerWherein ω is the angular frequency;

step 4, calculating P_PRBSSelf-power spectrum G of_PP(ω) and P_PRBSAnd I_qrCross power spectrum G of_PI(ω) to calculate a frequency series of transfer functions of the rotor-side controller q-axis decoupling model

Step 5, replacing s in the step 1 with j omega to obtain a frequency domain transfer function H₁(ω,θ₁) And H₂(ω,θ₂) (ii) a By H₁(ω,θ₁) Fitting the frequency series of the transfer function in step 3To identify a parameter theta to be identified in a d-axis decoupling model₁(ii) a By H₂(ω,θ₂) Fitting the frequency series of the transfer function in step 4To identify a parameter theta to be identified in a q-axis decoupling model₂(ii) a Wherein j is an imaginary number;

step 6, identifying the hypothesis L before_r' known, output θ obtained in step 5₁In [ K ]_p3,K_i3,K_p4,K_i4]And theta₂In [ K ]_p1,K_i1,K_p2,K_i2]The result of the identification.

As a further optimization scheme of the invention, in step 3 or 4, the calculation process of the self-power spectrum is to calculate the autocorrelation function of the self-power spectrum, and then perform fourier transform on the calculated autocorrelation function, so as to obtain the self-power spectrum; the cross-power spectrum is calculated by firstly calculating the cross-correlation function of the two and then carrying out Fourier transform on the calculated cross-correlation function to obtain the cross-power spectrum.

As a further optimization scheme of the invention, the fitting in the step 5 adopts a damped nonlinear least square optimization method.

As a further optimization scheme of the invention, theta in the d-axis decoupling model is identified in step 5₁The objective function of the damping nonlinear least square optimization method is as follows:

$\min \underset{i=1}{\overset{n}{\mathrm{\Σ}}}{|{\hat{H}}_1\left(\mathrm{\ω}\right)-H_1(\mathrm{\ω},{\mathrm{\θ}}_1)|}^2$

wherein the parameter search range is theta₁∈[θ_1max,θ_1min]，θ_1maxAnd theta_1minRespectively is a parameter theta to be identified₁A maximum set value and a minimum set value; n is identification frequency window I_drThe total number of points of the transfer function sequence relative to the voltage deviation amount, i is an integer, i ∈ [1, n]。

As a further optimization scheme of the invention, theta in a q-axis decoupling model is identified in step 5₂The objective function of the damping nonlinear least square optimization method is as follows:

$\min \underset{k=1}{\overset{l}{\mathrm{\Σ}}}{|{\hat{H}}_2\left(\mathrm{\ω}\right)-H_2(\mathrm{\ω},{\mathrm{\θ}}_2)|}^2$

wherein the parameter search range is theta₂∈[θ_2max,θ_2min]，θ_2maxAnd theta_2minRespectively is a parameter theta to be identified₂A maximum set value and a minimum set value; l is the identification frequency window I_qrThe total number of points of the transfer function sequence relative to the voltage deviation value, k is an integer, k ∈ [1, l]。

As a further optimization scheme of the invention, the reference angular frequency omega in the step 1_B＝2πf_NWherein f is_NIs the nominal frequency of the system.

As a further optimization scheme of the invention, the rotor transient inductance in step 1Wherein L is_ssIs a stator inductance, L_ss＝L_sσ+L_m，L_rrIs the rotor inductance, L_rr＝L_rσ+L_m，L_sσAnd L_rσStator leakage inductance and rotor leakage inductance, L, respectively_mThe stator and the rotor are mutually inducted.

Compared with the prior art, the invention adopts the technical scheme that the generator rotor current is used as an output signal, so that a d-axis and q-axis decoupling model of the rotor side controller is realized, and the model structure is simple; the pseudo-random binary signal is applied to the reference signal, the test can be carried out when the normal work of the wind turbine generator is not influenced, the test method is practical, and the identification result is more in line with the actual running state of the fan.

Drawings

Fig. 1 is a structural diagram of a simulation test system for grid connection of a wind driven generator based on a DFIG.

Fig. 2 is a flow chart of the operation of the present invention.

Fig. 3 is an original model structure diagram of a doubly-fed wind turbine rotor-side controller.

Fig. 4 is a d-axis and q-axis decoupling model structure diagram of a doubly-fed wind turbine generator rotor side controller, wherein (a) is the d-axis decoupling model structure diagram; (b) the structure diagram of the q-axis decoupling model is shown.

Detailed Description

The technical scheme of the invention is further described in detail by combining the drawings and the specific embodiments:

in this embodiment, a doubly-fed induction machine (DFIG) is used as a basisThe Wind turbine generator system is directly connected to the grid, a simulation system is shown in FIG. 1, the simulation system is built in Matlab2012b software, all elements in the system are taken from Wind power plant calculation examples (Wind Farm (DFIG Phasor Model) Demo) carried by Matlab, element parameters are taken as default values, wherein β represents the pitch angle of a Wind turbine blade, and omega is_rRepresenting the DFIG rotor speed, I_rRepresenting the current in the rotor winding, V_dr、V_qrD-axis and q-axis components representing rotor excitation voltage, respectively; v_DCRepresenting the voltage of the DC side capacitor; v_gAnd I_gIs the output voltage and current of the grid-side converter; v_dg、V_qgRepresenting the d-axis and q-axis components of the grid-side converter output voltage, respectively.

As shown in fig. 2, a frequency domain method for identifying parameters of a rotor-side controller of a doubly-fed wind turbine generator specifically includes the following steps:

step 1, determining a d-axis and q-axis decoupling model and a parameter to be identified of a rotor side controller of the doubly-fed wind turbine generator.

The master model of the rotor-side controller is shown in FIG. 3, where K_p1And K_i1Respectively controlling the proportion and integral gain of an outer ring for the q-axis active power; k_p2And K_i2Proportional and integral gains of the q-axis active power control inner ring are respectively obtained; k_p3And K_i3Proportional and integral gains of the outer loop are controlled for the d-axis voltage respectively; k_p4And K_i4Proportional and integral gains of the inner d-axis voltage control loop respectively; i is_drAnd I_qrD-axis component and q-axis component of the rotor side current respectively; i is_dsAnd I_qsD-axis component and q-axis component of the stator side current respectively; i is_{dr_ref}And I_{qr_ref}Reference values for the d-axis component and the q-axis component of the rotor-side current, respectively; p_sAnd V_sActive power and stator voltage, respectively; p_{s_ref}And V_{s_ref}Reference values for active power and stator voltage, respectively; omega_sAnd ω_rThe angular speed of the power grid and the angular speed of the generator are respectively; l is_rrIs the rotor inductance, L_rr＝L_rσ+L_m，L_rσFor rotor leakage inductance, L_mIs the mutual inductance between stator and rotor windings.

Selecting d-axis component I of rotor-side current_drAnd q-axis component I_qrAs the output quantities, d-axis and q-axis decoupling models of the rotor-side controller are obtained as shown in (a) and (b) in fig. 4.

The transfer function of the d-axis decoupling model is:

$H_1(s,{\mathrm{\θ}}_1)=\frac{I_{\mathrm{dr}}\left(s\right)}{V_{s\_\mathrm{ref}}-V_s}=\frac{\frac{K_{p3}{K}_{p4}}{L_r^{\′}/{\mathrm{\ω}}_B}{s}^2+\frac{K_{p4}{K}_{i3}+K_{p3}{K}_{i4}}{L_r^{\′}/{\mathrm{\ω}}_B}s+\frac{K_{i3}{K}_{i4}}{L_r^{\′}/{\mathrm{\ω}}_B}}{s(s^2+\frac{K_{p4}}{L_r^{\′}/{\mathrm{\ω}}_B}s+\frac{K_{i4}}{L_r^{\′}/{\mathrm{\ω}}_B})}---\left(1\right)$

transfer function H of q-axis decoupling model₂(s,θ₂) Comprises the following steps:

$H_2(s,{\mathrm{\θ}}_2)=\frac{I_{\mathrm{qr}}\left(s\right)}{P_{s\_\mathrm{ref}}-P_s}=\frac{\frac{K_{p2}{K}_{p1}}{L_r^{\′}/{\mathrm{\ω}}_B}{s}^2+\frac{K_{p2}{K}_{i1}+K_{p1}{K}_{i2}}{L_r^{\′}/{\mathrm{\ω}}_B}s+\frac{K_{i2}{K}_{i1}}{L_r^{\′}/{\mathrm{\ω}}_B}}{s(s^2+\frac{K_{p2}}{L_r^{\′}/{\mathrm{\ω}}_B}s+\frac{K_{i2}}{L_r^{\′}/{\mathrm{\ω}}_B})}---\left(2\right)$

in the formula, s is a Laplace operator; theta₁For the d-axis decoupling model to be identified, θ₁＝[K_p3,K_i3,K_p4,K_i4,L_r′]^T；ω_BIs a reference angular frequency; l is_r' is rotor transient inductance; theta₂For the q-axis decoupling model to be identified, θ₂＝[K_p1,K_i1,K_p2,K_i2,L_r′]^T。

Wherein the reference angular frequency omega_B＝2πf_N，f_NIs the nominal frequency of the system; rotor transient inductanceL_ssIs a stator inductance, L_ss＝L_sσ+L_m，L_sσThe leakage inductance of the stator is obtained.

Step 2, in the reference value V of the stator voltage_{s_ref}Upper applied frequency bandwidth of [0.001,100]Pseudo-random binary voltage signal V with Hz and sampling frequency of 1000Hz_PRBSSynchronous measurement of V_PRBSAnd d-axis component I of rotor-side current_dr(ii) a At the active power reference value P_refUpper applied frequency bandwidth of [0.001,100]Pseudo-random binary active power signal P with Hz and sampling frequency of 1000Hz_PRBSSynchronous measurement of P_PRBSAnd q-axis component I of rotor-side current_qr。

Step 3, calculating V_PRBSSelf-power spectrum G of_VV(ω) and V_PRBSAnd I_drCross power spectrum G of_VI(ω) to calculate a frequency series of transfer functions of the d-axis decoupling model of the rotor-side controllerWhere ω is the frequency.

In this embodiment, V is calculated_PRBSWhen the spectrum is self-powered, V is calculated first_PRBSIs self-correlation function ofR_V(m) then to R_V(m) Fourier transform to obtain the self-power spectrum G_VV(ω), specifically as follows:

R_V(m)＝E[V_PBRS(t)V_PBRS(t+m)](3)

$G_{\mathrm{VV}}\left(\mathrm{\ω}\right)=\underset{m=-\∞}{\overset{\∞}{\mathrm{\Σ}}}{R}_V\left(m\right)e^{-\mathrm{j\ωm}{T}_s}---\left(4\right)$

in the formula, V_PBRS(t) is V_PBRSTime series of (1), V_PBRS(t + m) is V after delaying m time intervals_PBRSTime series of (E [ ·)]To take the mean, m is the delay, T_sIs the data sampling interval.

In this embodiment, V is calculated_PRBSAnd I_drWhen cross power spectrum is obtained, V is calculated first_PRBSAnd I_drCross correlation function R of_VI(m) then to R_VI(m) Fourier transform to obtain cross-power spectrum G_VI(ω), specifically as follows:

R_VI(m)＝E[V_PBRS(t)I_dr(t+m)](5)

$G_{\mathrm{VI}}\left(\mathrm{\ω}\right)=\underset{m=-\∞}{\overset{\∞}{\mathrm{\Σ}}}{R}_{\mathrm{VI}}\left(m\right)e^{-\mathrm{j\ωm}{T}_s}---\left(6\right)$

in the formula I_dr(t + m) is I delayed by m time intervals_dr(t) time series, I_dr(t) is I_drTime series.

And 4, step 4: calculating P_PRBSSelf-power spectrum G of_PP(ω) and P_PRBSAnd I_qrCross power spectrum G of_PI(omega) to derive a transfer function sequence of a rotor side controller q-axis decoupling model

In this embodiment, P is calculated_PRBSWhen the spectrum of the self-power is obtained, P is calculated first_PRBSIs the autocorrelation function R of_P(m) then to R_P(m) Fourier transform to obtain the self-power spectrum G_PP(ω), specifically as follows:

R_P(m)＝E[P_PBRS(t)P_PBRS(t+m)](7)

$G_{\mathrm{PP}}\left(\mathrm{\ω}\right)=\underset{m=-\∞}{\overset{\∞}{\mathrm{\Σ}}}{R}_P\left(m\right)e^{-\mathrm{j\ωm}{T}_s}---\left(8\right)$

in the formula, P_PBRS(t) is P_PBRSTime series of (1), P_PBRS(t + m) is P delayed by m time intervals_PBRSTime series, E [. C]To take the mean, m is the delay, T_sIs the data sampling interval.

In this embodiment, P is calculated_PRBSAnd I_qrWhen cross power spectrum of (2), first calculate P_PRBSAnd I_qrCross correlation function R of_PI(m) then to R_PI(m) Fourier transform to obtain cross-power spectrum G_PI(ω), specifically as follows:

R_PI(m)＝E[P_PBRS(t)I_qr(t+m)](9)

$G_{\mathrm{PI}}\left(\mathrm{\ω}\right)=\underset{m=-\∞}{\overset{\∞}{\mathrm{\Σ}}}{R}_{\mathrm{PI}}\left(m\right)e^{-\mathrm{j\ωm}{T}_s}---\left(10\right)$

in the formula I_qr(t + m) is I delayed by m time intervals_qr(t) time series, I_qr(t) is I_qrTime series.

And 5: replacing s in the step 1 by j omega to obtain a frequency domain transfer function H₁(ω,θ₁) And H₂(ω,θ₂) (ii) a By H₁(ω,θ₁) Fitting the sequence of transfer functions in step 3Accordingly, theta in the d-axis decoupling model is identified₁(ii) a By H₂(ω,θ₂) Fitting the sequence of transfer functions in step 4From this, identify theta in q-axis decoupling model₂。

The fitting method in this embodiment adopts a damped nonlinear least squares optimization method, and then:

theta in d-axis decoupling model₁The objective function of (a) is:

$\min \underset{i=1}{\overset{n}{\mathrm{\Σ}}}{|{\hat{H}}_1\left(\mathrm{\ω}\right)-H_1(\mathrm{\ω},{\mathrm{\θ}}_1)|}^2---\left(11\right)$

wherein the parameter search range is theta₁∈[θ_1max,θ_1min]，θ_1maxAnd theta_1minRespectively is a parameter theta to be identified₁A maximum set value and a minimum set value; n is identification frequency window I_drThe total number of points of the transfer function sequence relative to the voltage deviation amount, i is an integer, i ∈ [1, n]。

Theta in q-axis decoupling model₂The objective function of (a) is:

$\min \underset{k=1}{\overset{l}{\mathrm{\Σ}}}{|{\hat{H}}_2\left(\mathrm{\ω}\right)-H_2(\mathrm{\ω},{\mathrm{\θ}}_2)|}^2---\left(12\right)$

wherein the parameter search range is theta₂∈[θ_2max,θ_2min]，θ_2maxAnd theta_2minRespectively is a parameter theta to be identified₂A maximum set value and a minimum set value; l is the identification frequency window I_qrThe total number of points of the transfer function sequence relative to the voltage deviation value, k is an integer, k ∈ [1, l]。

Step 6, identifying the hypothesis L before_r' known, output θ obtained in step 5₁In [ K ]_p3,K_i3,K_p4,K_i4]And theta₂In [ K ]_p1,K_i1,K_p2,K_i2]The result of the identification.

In this embodiment, a calculation example of directly accessing a doubly-fed wind turbine generator to an infinite system as shown in fig. 1 is adopted, and system parameters are shown in table 1, where pu is a per unit value.

TABLE 1

Hypothesis before identification parameter L_r' As known, the search range for the remaining parameters is set to [ -80%, + 200% of the true value]. Fitting an actual measurement transfer function based on equation (1)Frequency sequence of (2), parameter theta of L-M optimization method₁＝[K_p3,K_i3,K_p4,K_i4]The identification results are shown in Table 2; fitting the measured transfer function based on equation (2)Frequency sequence of (2), parameter theta of L-M optimization method₂＝[K_p1,K_i1,K_p2,K_i2]The identification results are shown in Table 3. Wherein,

TABLE 2

Parameter(s)	Truth value	Search scope	The result of the recognition	Error/%)
					Kp3/pu	1.25	[0.25,2.5]	1.4561	16.4880
Ki3/1/s	300	[60,600]	298.8650	-0.3783
					Kp4/pu	0.3	[0.06,0.6]	0.2425	-19.1667
Ki4/1/s	8	[1.6,16]	5.3548	-33.0650

TABLE 3

Parameter(s)	Truth value	Search scope	The result of the recognition	Error/%)
					Kp1/pu	1	[0.2,2]	1.1429	14.2900
Ki1/1/s	100	[20,200]	109.2178	9.2178
					Kp2/pu	0.3	[0.06,0.6]	0.2562	-14.6000
Ki2/1/s	8	[1.6,16]	6.4615	-19.2312

From tables 2 and 3, it can be seen that L is assumed_r' in the known Presence, parameter K_i3、K_i1High identification accuracy of, K_p1～K_p4The parameter K of the recognition accuracy_i2And K_i4The identification error of (2) is the largest. However due to K_i2And K_i4The method belongs to the inner ring control parameters of a rotor side controller, and because the dynamic influence of the inner ring controller on the dynamic state of the double-fed wind turbine generator is small, the error of the inner ring controller is acceptable for analyzing the dynamic state of the double-fed wind turbine generator.

The above description is only an embodiment of the present invention, but the scope of the present invention is not limited thereto, and any person skilled in the art can understand that the modifications or substitutions within the technical scope of the present invention are included in the scope of the present invention, and therefore, the scope of the present invention should be subject to the protection scope of the claims.

Claims (7)

1. A frequency domain method for identifying parameters of a rotor side controller of a doubly-fed wind turbine generator is characterized by comprising the following steps of:

step 1, determining a d-axis and q-axis decoupling model and a parameter to be identified of a rotor side controller of a doubly-fed wind turbine generator, specifically:

the d-axis decoupling model of the rotor-side controller takes the voltage deviation amount as input and the d-axis component I of the rotor-side current_drAs an output, its transfer function H₁(s,θ₁) The expression of (a) is:

$H_1(s,{\mathrm{\θ}}_1)=\frac{I_{\mathrm{dr}}\left(s\right)}{V_{s\_\mathrm{ref}}-V_s}=\frac{\frac{K_{p3}{K}_{p4}}{L_r^{\′}/{\mathrm{\ω}}_B}{s}^2+\frac{K_{p4}{K}_{i3}+K_{p3}{K}_{i4}}{L_r^{\′}/{\mathrm{\ω}}_B}s+\frac{K_{i3}{K}_{i4}}{L_r^{\′}/{\mathrm{\ω}}_B}}{s(s^2+\frac{K_{p4}}{L_r^{\′}/{\mathrm{\ω}}_B}s+\frac{K_{i4}}{L_r^{\′}/{\mathrm{\ω}}_B})}---\left(1\right)$

in the formula, s is a Laplace operator; theta₁For the d-axis decoupling model to be identified, θ₁＝[K_p3,K_i3,K_p4,K_i4,L_r′]^T；I_dr(s) is the d-axis component I_drThe Ralsberg transform of (1); v_{s_ref}Is a stator voltage reference value; v_sIs the stator voltage; omega_BIs a reference angular frequency; k_p3And K_i3Proportional and integral gains, K, of the outer loop for d-axis voltage control, respectively_p4And K_i4Proportional and integral gains of the inner d-axis voltage control loop respectively; l is_r' is rotor transient inductance;

the q-axis decoupling model of the rotor side controller takes the active power deviation value as input and the q-axis component I of the rotor side current_qrAs an output, its transfer function H₂(s,θ₂) The expression of (a) is:

$H_2(s,{\mathrm{\θ}}_2)=\frac{I_{\mathrm{qr}}\left(s\right)}{P_{s\_\mathrm{ref}}-P_s}=\frac{\frac{K_{p2}{K}_{p1}}{L_r^{\′}/{\mathrm{\ω}}_B}{s}^2+\frac{K_{p2}{K}_{i1}+K_{p1}{K}_{i2}}{L_r^{\′}/{\mathrm{\ω}}_B}s+\frac{K_{i2}{K}_{i1}}{L_r^{\′}/{\mathrm{\ω}}_B}}{s(s^2+\frac{K_{p2}}{L_r^{\′}/{\mathrm{\ω}}_B}s+\frac{K_{i2}}{L_r^{\′}/{\mathrm{\ω}}_B})}---\left(2\right)$

in the formula, theta₂For the q-axis decoupling model to be identified, θ₂＝[K_p1,K_i1,K_p2,K_i2,L_r′]^T；I_qr(s) is the q-axis component I_qrThe Ralsberg transform of (1); p_{s_ref}Is a reference value of active power; p_sActive power; k_p1And K_i1Respectively controlling the proportion and integral gain of an outer ring for the q-axis active power; k_p2And K_i2Proportional and integral gains of the q-axis active power control inner ring are respectively obtained;

step 2, in the reference value V of the stator voltage_{s_ref}Applied voltage pseudo-random signal V_PRBSAnd synchronously collecting V_PRBSAnd d-axis component I of rotor-side current_dr(ii) a At reference value P of active power_{s_ref}To which an active power pseudo-random signal P is applied_PRBSAnd synchronously collecting P_PRBSAnd q-axis component I of rotor-side current_qr；

Step 3, calculating V_PRBSSelf-power spectrum G of_VV(ω) and V_PRBSAnd I_drCross power spectrum G of_VI(ω) to calculate a frequency series of transfer functions of the d-axis decoupling model of the rotor-side controllerWherein ω is the angular frequency;

step 4, calculating P_PRBSSelf-power spectrum G of_PP(ω) and P_PRBSAnd I_qrCross power spectrum G of_PI(ω) to calculate a frequency series of transfer functions of the rotor-side controller q-axis decoupling model

Step 5, replacing the Laplace operator s in the step 1 with j omega to obtain frequency domain transmissionTransfer function H₁(ω,θ₁) And H₂(ω,θ₂) (ii) a By H₁(ω,θ₁) Fitting the frequency series of the transfer function in step 3To identify a parameter theta to be identified in a d-axis decoupling model₁(ii) a By H₂(ω,θ₂) Fitting the frequency series of the transfer function in step 4To identify a parameter theta to be identified in a q-axis decoupling model₂(ii) a Wherein j is an imaginary number;

step 6, identifying the hypothesis L before_r' known, output θ obtained in step 5₁In [ K ]_p3,K_i3,K_p4,K_i4]And theta₂In [ K ]_p1,K_i1,K_p2,K_i2]The result of the identification.

2. The frequency domain method for identifying the parameters of the rotor side controller of the doubly-fed wind turbine generator set according to claim 1, wherein the self-power spectrum is calculated by firstly calculating an autocorrelation function of the doubly-fed wind turbine generator set, and then performing Fourier transform on the calculated autocorrelation function to obtain the self-power spectrum; the cross-power spectrum is calculated by firstly calculating the cross-correlation function of the two and then carrying out Fourier transform on the calculated cross-correlation function to obtain the cross-power spectrum.

3. The frequency domain method for identifying the parameters of the rotor side controller of the doubly-fed wind turbine generator set according to claim 1, wherein the fitting in the step 5 adopts a damped nonlinear least square optimization method.

4. The frequency domain method for identifying the parameters of the rotor-side controller of the doubly-fed wind turbine generator set according to claim 3, wherein the identification of θ in the d-axis decoupling model in the step 5 is performed₁The objective function of the damping nonlinear least square optimization method is as follows:

$\min \underset{i=1}{\overset{n}{\mathrm{\Σ}}}{|{\hat{H}}_1\left(\mathrm{\ω}\right)-H_1(\mathrm{\ω},{\mathrm{\θ}}_1)|}^2$

wherein the parameter search range is theta₁∈[θ_1max,θ_1min]，θ_1maxAnd theta_1minRespectively is a parameter theta to be identified₁A maximum set value and a minimum set value; n is identification frequency window I_drThe total number of points of the transfer function sequence relative to the voltage deviation amount, i is an integer, i ∈ [1, n]。

5. The frequency domain method for identifying the parameters of the rotor-side controller of the doubly-fed wind turbine generator set according to claim 3, wherein the identification of θ in the q-axis decoupling model in step 5 is performed₂And the damping nonlinear least square optimization method comprises the following objective functions:

$\min \underset{k=1}{\overset{l}{\mathrm{\Σ}}}{|{\hat{H}}_2\left(\mathrm{\ω}\right)-H_2(\mathrm{\ω},{\mathrm{\θ}}_2)|}^2$

wherein the parameter search range is theta₂∈[θ_2max,θ_2min]，θ_2maxAnd theta_2minRespectively is a parameter theta to be identified₂A maximum set value and a minimum set value; l is the identification frequency window I_qrThe total number of points of the transfer function sequence relative to the voltage deviation value, k is an integer, k ∈ [1, l]。

6. The frequency domain method for identifying the parameters of the rotor-side controller of the doubly-fed wind turbine generator set according to claim 1, wherein the reference angular frequency ω in step 1 is the frequency of the reference angular frequency ω_B＝2πf_NWherein f is_NIs the nominal frequency of the system.

7. The frequency domain method for identifying the parameters of the rotor-side controller of the doubly-fed wind turbine generator set according to claim 1, wherein the rotor transient inductance in the step 1 isWherein L is_ssIs a stator inductance, L_ss＝L_sσ+L_m，L_rrIs the rotor inductance, L_rr＝L_rσ+L_m，L_sσAnd L_rσStator leakage inductance and rotor leakage inductance, L, respectively_mThe stator and the rotor are mutually inducted.