CN105425611B - A kind of turbine-generator units Excitation System Parameter Identification of Synchronous method - Google Patents

Abstract

The invention discloses a kind of turbine-generator units Excitation System Parameter Identification of Synchronous method, for accurately obtaining turbine-generator units parameters of excitation system.Phantom is set up according to turbine-generator units excitation system, then set up according to this analogue system and use the real system output weighted error quadratic sum with identification system output as the object function of parameter identification, use the method for optimizing of present invention design to solve object function and obtain optimal control parameter.The turbine-generator units Excitation System Parameter Identification of Synchronous method of present invention design, uses a kind of novel heuristic value optimization object function, may search for less target function value, can obtain more accurate identified parameters.More accurate identified parameters makes identification system output coincide preferably with real system output.

Description

A kind of turbine-generator units Excitation System Parameter Identification of Synchronous method

Technical field

The invention belongs to parameters of electric power system optimisation technique field, send out more particularly, to a kind of water wheels Group of motors Excitation System Parameter Identification of Synchronous method.

Background technology

Under or accident condition properly functioning in power system, excitation system of hydrogenerator plays important Effect.It has control voltage, controls the distribution of reactive power, improves hydrogenerator parallel running Stability, improve the ability of the stability of power system.But power department is frequently used allusion quotation in the past Shape parameter, it is difficult to the real dynamic characteristic reflecting power system.And hydrogenerator in power system The model of group excitation system and the order of accuarcy of parameter are directly connected to the safe and stable operation water of system Flat.Therefore, the excitation system running scene carries out identification experiment, enters according to the data of collection in worksite Row Excitation System Parameter Identification of Synchronous is a very important job.

Excitation System Parameter Identification of Synchronous method conventional in China's power system mostly is frequency domain or time domain identification side Method, including fast Fourier transform (fast Fourier transform, FFT), method of least square (Least Squares method, LSE), piecewise linearity polynomial function method (Piecewise liner polynomial function,PLPE).Above-mentioned discrimination method principle is clear, simple and easy to do, the most at China's electric power System is on the actual application.But, these methods are the parameter identifications of linear system, but It is that excitation system often exists some nonlinear elements, such as amplitude limit link etc., therefore in identification non-thread Deficiency is yet suffered from during property excitation system.If will identification excitation system in aforementioned manners, need to make line Propertyization processes, and the result of institute's identification cannot accurately reflect the dynamic characteristic of nonlinear system.People is had for this Propose to be used for intelligent optimization algorithm the parameter identification of excitation system.Particle swarm optimization algorithm (particle Swarm optimization, PSO), gravitation search algorithm (gravitational search algorithm, GSA) parameter identification of excitation system it is applied to, the shortcoming overcoming frequency domain method and time domain method, Can effectively pick out systematic parameter.But PSO and GSA search optimal solution during still Shortcomings, is easily trapped into local optimum, causes finally cannot searching global optimum.

Summary of the invention

For the deficiency of traditional method, the present invention proposes a kind of based on novel heuristic value Turbine-generator units Excitation System Parameter Identification of Synchronous method, can effectively pick out systematic parameter, have very Good practical value.

To achieve these goals, the invention provides a kind of turbine-generator units parameters of excitation system to distinguish Knowledge method, comprises the steps:

Step (1): set up turbine-generator units excitation system phantom, determine parameter to be identified. Turbine-generator units excitation system structure is as it is shown in figure 1, described system includes PID controller, puts Big unit, exciter, hydrogenerator, measuring unit.Concrete, this excitation system is an allusion quotation The feedback control system of type, the hydrogenerator set end voltage arrived by measuring unit measurement is with given Reference voltage compares, and obtains the side-play amount of system output, and this side-play amount produces control through PID controller Signal processed, then act on exciter after amplifying unit amplifies, it is achieved the regulation to excitation voltage, Reach to regulate further the purpose of hydrogenerator set end voltage.Conventional controller includes, PID controller And the modified model PID controller such as Fractional Order PID Controller, fuzzy controller.For convenience of saying Bright, the present invention is using PID controller as Excitation Controller.V in Fig. 1_refFor reference voltage, V_cFor PID controller exports, V_RExport for amplifying unit, V_FExport for exciter, V_tFor hydrogenerator machine Terminal voltage, V_SExport for sensor.Wherein k_A,τ_AIt is respectively amplifying unit gain and time constant, k_E,τ_E It is respectively exciter gain and time constant, k_G,τ_GIt is respectively hydrogenerator gain and time constant, k_S,τ_SIt is respectively sensor gain and time constant.The parameter vector needing identification is $\widehat{\mathrm{\θ}}=\[k_A,{\mathrm{\τ}}_A,k_E,{\mathrm{\τ}}_E,k_G,{\mathrm{\τ}}_G,k_S,{\mathrm{\τ}}_S\];$

Step (2): gather actual excitation system dynamic process data.Actual excitation system is carried out electricity Pressure step disturbance test, gathers actual excitation system dynamic process data, and dynamic process data include respectively Link exports, i.e. amplifying unit output V_R, exciter output V_F, hydrogenerator set end voltage V_tWith Sensor output V_S；

Step (3): set up Excitation System Parameter Identification of Synchronous object function.Use real system output and distinguish The weighted error quadratic sum of knowledge system output is as the object function of parameter identification.Object function defines such as Under:

${\mathrm{minf}}_{WMSE}\left(\widehat{\mathrm{\θ}}\right)=\underset{k=1}{\overset{N_s}{\Σ}}\underset{j=1}{\overset{s}{\Σ}}{w}_j{\left(y_j\right(k)-{\hat{y}}_j(k\left)\right)}^2$

Wherein N_sExporting sampling number for system, s is that system exports number, w=[w₁,w₂,w₃,w₄] power Weight,It it is parameter to be identified.Under identical systems inputs, for actual system System is output as y_j(k)∈{V_R(k),V_F(k),V_t(k),V_S(k)}；Identification system is output as ${\hat{y}}_j\left(k\right)\∈\{{\hat{V}}_R\left(k\right),{\hat{V}}_F\left(k\right),{\hat{V}}_t\left(k\right),{\hat{V}}_S\left(k\right)\}.$ Wherein,It is the function of system parameter to be identified, when Parameter vectorWhen changing, the analogue system utilizing step (1) to set up obtains four systems output, I.e. four groups curves of output, with discrete series representation be:Meter Calculate corresponding target function value, by minimization object function, system of asking for parameter to be identified；

Step (4): use heuristic value to solve the object function of Excitation System Parameter Identification of Synchronous, Obtain unidentified system parameter.Identification step is as follows:

Step 1: algorithm initialization: arrange algorithm parameter, including population size N, maximum iteration time T, individual random search quantity N_l, eliminate range coefficient σ, skip threshold p；Determine that search is to be identified The span of parameter is [B_L,B_U], concrete k_A∈[k_A,min,k_A,max], τ_A∈[τ_A,min,τ_A,max], k_E∈[k_E,min,k_E,max], τ_E∈[τ_E,min,τ_E,max]k_G∈[k_G,min,k_G,max], τ_G∈[τ_G,min,τ_G,max], k_S∈[k_S,min,k_S,max], τ_S∈[τ_S,min,τ_S,max], i.e. B_L=[k_A,min,τ_A,min,k_E,min,τ_E,min,k_G,min,τ_G,min,k_S,min,τ_S,min] represent excitation system parameter to be identified Little value, B_U=[k_A,max,τ_A,max,k_E,max,τ_E,max,k_G,max,τ_G,max,k_S,max,τ_S,max] represent excitation system ginseng to be identified The maximum of number.At solution space [B_L,B_UThe position vector of all individualities in random initializtion colony in], One group position vector is expressed as X_i=[k_A,τ_A,k_E,τ_E,k_G,τ_G,k_S,τ_S].Arranging maximum iteration time is T, Make current iteration number of times t=0；

Step 2: calculate individual target function valueAnd find colony Object function minima, the individuality with minimum target functional value is defined as current optimum individual X_B(t)；

Step 3: to all individual X_i(t) (i=1 ..., N) carry out individual random search；

The individual searching times l=0 of Step 3.1: order；

Step 3.2: look around a positionCalculate $X_i^{play}\left(t\right),i=1,\mathrm{...},N:$

$X_i^{play}\left(t\right)=X_i\left(t\right)+rand\·{\mathrm{\ϵ}}_{play}$

Rand is random number between (0,1), ε_playFor looking around step-length, ε_play=0.1 | | B_U-B_L||；

Step 3.3: calculate next current location

$\left\{\begin{array}{ccc}{X}_i^{self}\left(t\right)=X_i\left(t\right)+rand\&CenterDot;\frac{X_i^{play}\left(t\right)-X_i\left(t\right)}{\left|\right|X_i^{play}\left(t\right)-X_i\left(t\right)\left|\right|}\&CenterDot;{\mathrm{\&epsiv;}}_{step}& if& f\left(X_i^{play}\right(t\left)\right)<f\left(X_i\right(t\left)\right)\\ X_i^{self}\left(t\right)=X_i\left(t\right)& if& f\left(X_i^{paly}\right(t\left)\right)\&GreaterEqual;f\left(X_i\right(t\left)\right)\end{array}\right.$

Rand is random number between (0,1), ε_stepFor inertia step-length, ε_step=0.2 | | B_U-B_L||；

Step 3.4:l=l+1, if l is < N_l, go to Step 3.2；Otherwise, Step 4 is gone to；

Step 4: update individual position vector X according to individual location updating formula_i(t), i=1 ..., N:

$\left\{\begin{array}{c}{\mathrm{\&delta;}}_i=|c_1\&CenterDot;X_B\left(t\right)-X_i\left(t\right)|\\ X_i(t+1)=X_B\left(t\right)+c_2\&CenterDot;{\mathrm{\&delta;}}_i\end{array}\right.$

Wherein, δ_iDistance vector with current optimum individual individual for middle i-th, random number c₁=2 rand, c₂=(2 rand-1) exp (-10 t/T)；It can thus be appreciated that c₁For the random number between (0,2), represent The charisma of excellent individuality, works as c₁> 1 time, represent current optimum individual power of influence strengthen, otherwise weaken； c₂For dynamic random number；

Step 5: judge individual the need of being eliminated and reinitializing:

Step 5.1: if i-th individuality meets formula, this individuality is eliminated and reinitializes:

$F_i^t>F_{ave}^t+\mathrm{\&omega;}\&CenterDot;(F_{ave}^t-F_{\min }^t),i=1,...,N$

Wherein,It is the t meansigma methods for population all individual goals functional value,It it is minimum mesh Offer of tender numerical value, ω be one with iterations the parameter of linear increment,Value model Enclose for [-σ, σ]；

Step 5.2: the individual initialization being eliminated:

X_i=rand (1, D) × (B_U-B_L)+B_L

Wherein, D is position vector dimension, D=8；

Step 6: judge whether that continuous p is not moved for current optimum individual position, if it is, Think population extinction, according to the population that formula following formula inverting reconstruct is new:

$X_i=X_B+rand\&times;\frac{R^2}{{\mathrm{\&delta;}}_i},i=1,2,...,N$

Wherein R is radius of inversion, R=0.1 | | B_U-B_L||；Rand is random number between (0,1), p For skip threshold；

Step 7:t=t+1, if t > T, algorithm terminates, and exports the most current optimum individual of optimal location vector Position is as whole solution；Otherwise, Step 2 is proceeded to.Described optimal location vector is system ginseng to be identified Number vector.

Compared with prior art, when utilizing the method for the invention identification parameters of excitation system, can search Rope, to less target function value, can obtain more accurate identified parameters.More accurate identified parameters makes Obtain identification system output and coincide preferable with real system output.

Accompanying drawing explanation

Fig. 1 is that turbine-generator units excitation system of the present invention transmits function block diagram；

Fig. 2 is that the identification system output of the method for the invention is compared with real system output；

Fig. 3 is that the identification system output of tradition GSA method is compared with real system output.

Detailed description of the invention

In order to make the purpose of the present invention, technical scheme and advantage clearer, below in conjunction with accompanying drawing And embodiment, the present invention is further elaborated.Should be appreciated that described herein specifically Embodiment only in order to explain the present invention, is not intended to limit the present invention.Additionally, it is disclosed below Just may be used as long as technical characteristic involved in each embodiment of the present invention does not constitutes conflict each other To be mutually combined.

For effect of the present invention being described, below using a certain turbine-generator units excitation system as the present invention's The inventive method is further detailed by objective for implementation:

Step (1): set up turbine-generator units excitation system phantom, determine parameter to be identified. Excitation system structure chart is as it is shown in figure 1, include PID control unit, amplifying unit, exciter, water Turbine generator, measuring unit.V in Fig. 1_refFor reference voltage, V_cExport for PID controller, V_RFor Amplifying unit exports, V_FExport for exciter, V_tFor hydrogenerator set end voltage, V_SFor sensor Output.Wherein k_A,τ_AIt is respectively amplifying unit gain and time constant, k_E,τ_EIt is respectively exciter gain And time constant, k_G,τ_GIt is respectively hydrogenerator gain and time constant, k_S,τ_SIt is respectively sensor Gain and time constant.The parameter vector needing identification isTo water wheels Generating unit excitation analogue system carries out voltage step test, and gathering simulation system exports, including amplifying Unit output V_R, exciter output V_F, hydrogenerator set end voltage V_tV is exported with sensor_S.Electricity Pressure disturbance V_ref=1, simulation time is 3 seconds, and simulation step length is 0.01 second；

Step (2): at voltage disturbance V_refWhen=1, actual excitation system is carried out voltage step disturbance examination Testing, gather real system output, the sampling time is 3 seconds, and the sampling interval is 0.01 second.Dynamic process Data include that each link exports, it may be assumed that amplifying unit output V_R, exciter output V_F, hydrogenerator Set end voltage V_tV is exported with sensor_S；

Step (3): set up Excitation System Parameter Identification of Synchronous object function.Use real system output and distinguish The weighted error quadratic sum of knowledge system output is as the object function of parameter identification.Object function defines such as Under:

${\mathrm{minf}}_{WMSE}\left(\widehat{\mathrm{\&theta;}}\right)=\underset{k=1}{\overset{N_s}{\&Sigma;}}\underset{j=1}{\overset{s}{\&Sigma;}}{w}_j{\left(y_j\right(k\left){\hat{y}}_j\right(k\left)\right)}^2$

Wherein w=[w₁,w₂,w₃,w₄] it is each link weight.It it is ginseng to be identified Number.According to system output, the sensitivity of parameter being calculated weight, the weight obtained is: w=[0.3259, 0.4923,0.1346,0.0472].Concrete process list of references: the Li Chaoshun of calculating. Hydropower Unit controls system System identification and fault diagnosis research [D]. the Central China University of Science and Technology, 2010.DOI:10.7666/d.d152664.

Under identical systems inputs, real system is output as y_j(k)∈{V_R(k),V_F(k),V_t(k),V_S(k)}；Identification System is output as ${\hat{y}}_j\left(k\right)\&Element;\{{\hat{V}}_R\left(k\right),{\hat{V}}_F\left(k\right),{\hat{V}}_t\left(k\right),{\hat{V}}_S\left(k\right)\}.$ Wherein,It it is system parameter to be identified Function, works as parameter vectorWhen changing, the analogue system utilizing step (1) to set up obtains four System exports, and i.e. four groups curves of output with discrete series representation are: Calculate corresponding target function value, by minimization object function, system of asking for parameter to be identified.

Step (4): use heuristic value to solve the object function of Excitation System Parameter Identification of Synchronous, Obtain unidentified system parameter.

Step 1: algorithm initialization: arrange algorithm parameter, including population size N, maximum iteration time T, individual random search quantity N_l, eliminate range coefficient σ, skip threshold p；Determine that search is to be identified The span of parameter is [B_L,B_U], concrete k_A∈[k_A,min,k_A,max], τ_A∈[τ_A,min,τ_A,max], k_E∈[k_E,min,k_E,max], τ_E∈[τ_E,min,τ_E,max]k_G∈[k_G,min,k_G,max], τ_G∈[τ_G,min,τ_G,max], k_S∈[k_S,min,k_S,max], τ_S∈[τ_S,min,τ_S,max], i.e. B_L=[k_A,min,τ_A,min,k_E,min,τ_E,min,k_G,min,τ_G,min,k_S,min,τ_S,min] represent excitation system parameter to be identified Little value, B_U=[k_A,max,τ_A,max,k_E,max,τ_E,max,k_G,max,τ_G,max,k_S,max,τ_S,max] represent excitation system ginseng to be identified The maximum of number.At solution space [B_L,B_UThe position vector of all individualities in random initializtion colony in], One group position vector is expressed as X_i=[k_A,τ_A,k_E,τ_E,k_G,τ_G,k_S,τ_S].Arranging maximum iteration time is T, Make current iteration number of times t=0；

Algorithm parameter is arranged: population scale N=30, maximum iteration time T=400, eliminates range coefficient σ=0.01, individual random search number N_l=2, jump threshold values is p=100, other value default setting.Its Middle B_L=[5,0.05,0.5,0.2,0.5,0.5,0.5,0.05], B_U=[20,0.2,2,0.8,2,2,2,0.2]；

Step 2: calculate individual target function valueCalculate colony's mesh Scalar functions minima, the individuality with minimum target functional value is defined as current optimum individual X_B(t)；

Step 3: to all individual X_i(t) (i=1 ..., N) carry out individual random search:

The individual searching times l=0 of Step 3.1: order；

Step 3.2: look around a positionCalculate $X_i^{play}\left(t\right),i=1,...,N:$

$X_i^{play}\left(t\right)=X_i\left(t\right)+rand\&CenterDot;{\mathrm{\&epsiv;}}_{play}$

Rand is random number between (0,1), ε_playFor looking around step-length, ε_play=0.1 | | B_U-B_L||；

Step 3.3: calculate next current location

$\left\{\begin{array}{ccc}{X}_i^{self}\left(t\right)=X_i\left(t\right)+rand\&CenterDot;\frac{X_i^{play}\left(t\right)-X_i\left(t\right)}{\left|\right|X_i^{play}\left(t\right)-X_i\left(t\right)\left|\right|}\&CenterDot;{\mathrm{\&epsiv;}}_{step}& if& f\left(X_i^{play}\right(t\left)\right)<f\left(X_i\right(t\left)\right)\\ X_i^{self}\left(t\right)=X_i\left(t\right)& if& f\left(X_i^{paly}\right(t\left)\right)\&GreaterEqual;f\left(X_i\right(t\left)\right)\end{array}\right.$

Rand is random number between (0,1), ε_stepFor inertia step-length, ε_step=0.2 | | B_U-B_L||；

Step 3.4:l=l+1, if l is < N_l, go to Step 3.2；Otherwise, Step 4 is gone to；

Step 4: update individual position vector X according to individual location updating formula_i(t), i=1 ..., N；

$\left\{\begin{array}{c}{\mathrm{\&delta;}}_i=|c_1\&CenterDot;X_B\left(t\right)-X_i\left(t\right)|\\ X_i(t+1)=X_B\left(t\right)+c_2\&CenterDot;{\mathrm{\&delta;}}_i\end{array}\right.$

δ_iDistance vector with current optimum individual individual for middle i-th, random number c₁=2 rand, c₂=(2 rand-1) exp (-10 t/T)；It can thus be appreciated that c₁For the random number between (0,2), represent The charisma of excellent individuality, works as c₁> 1 time, represent current optimum individual power of influence strengthen, otherwise weaken； c₂For dynamic random number；

Step 5: judge individual the need of being eliminated and reinitializing:

Step 5.1: if i-th individuality meets formula, this individuality is eliminated and reinitializes:

$F_i^t>F_{ave}^t+\mathrm{\&omega;}\&CenterDot;(F_{ave}^t-F_{\min }^t),i=1,...,N$

Wherein,It is the t meansigma methods for population all individual goals functional value,It it is minimum mesh Offer of tender numerical value, ω be one with iterations the parameter of linear increment,Value model Enclose for [-σ, σ]；

Step 5.2: the individual initialization being eliminated:

X_i=rand (1, D) × (B_U-B_L)+B_L

Wherein, D is position vector dimension, D=8；

Step 6: judge whether that continuous p is not moved for current optimum individual position, if it is, Think population extinction, according to the population that formula following formula inverting reconstruct is new:

$X_i=X_B+rand\&times;\frac{R^2}{{\mathrm{\&delta;}}_i},i=1,2,...,N$

Wherein R is radius of inversion, R=0.1 | | B_U-B_L||；Rand is random number between (0,1), p For skip threshold；

Step 7:t=t+1, if t > T, algorithm terminates, and exports the most current optimum individual of optimal location vector Position is as whole solution；Otherwise, Step 2 is proceeded to.Described optimal location vector is system ginseng to be identified Number vector.

For illustrate effect of the present invention, by based on tradition GSA algorithm Excitation System Parameter Identification of Synchronous method with The inventive method contrasts.Wherein, the parameter of GSA algorithm is set to: gravitational constant G₀=30, Attenuation quotient β=9, population scale N=30, maximum iteration time T=400.

In order to test stability and the effectiveness of the method for the invention, parameter identification experiment is repeated 20 times, Result of the test statistical analysis is as shown in table 1 table 2.Table 1 is the excitation that different discrimination method identification obtains The meansigma methods of systematic parameter and standard deviation.Table 2 is that different discrimination method searches for the target function value obtained Statistical result, including minima, maximum, meansigma methods and standard deviation.

As can be seen from Table 1 and Table 2, for turbine-generator units Excitation System Parameter Identification of Synchronous problem, The method of the invention obtains less target function value, and caulocarpic standard deviation is less, calculates Method has higher stability.

The optimum identified parameters statistical result of table 1

Table 2 optimal objective function statistical result

Fig. 2 with Fig. 3 is respectively the phantom output that the parameter using distinct methods identification to obtain is corresponding With real system output contrast, correlation curve includes that amplifying unit exports V_R, exciter output V_F, water Turbine generator set end voltage V_tV is exported with sensor_S.Figure it is seen that pass through the inventive method The excitation system curve of output that identification obtains is identical with the curve of output of real system, illustrate based on The identification system of the method for the invention the most highly levels off to real system.Can from Fig. 3 curve of output To find out, the system obtained by GSA identification is the most variant with real system.Experimental result explanation The inventive method has more preferable identification effect than GSA, and algorithm performance is more excellent.

As it will be easily appreciated by one skilled in the art that and the foregoing is only presently preferred embodiments of the present invention, Not in order to limit the present invention, all made within the spirit and principles in the present invention any amendment, etc. With replacement and improvement etc., should be included within the scope of the present invention.

Claims (5)

1. a turbine-generator units Excitation System Parameter Identification of Synchronous method, it is characterised in that described method Comprise the steps:

Step (1): set up turbine-generator units excitation system phantom, determines parameter to be identified: Described turbine-generator units excitation system includes PID controller, amplifying unit, exciter, water wheels Electromotor, measuring unit；In described turbine-generator units excitation system, measured by measuring unit To hydrogenerator set end voltage compare with given reference voltage, obtain the side-play amount of system output, This side-play amount produces control signal through PID controller, then acts on after amplifying unit amplifies and encourage Magnetomechanical, it is achieved the regulation to excitation voltage, thus regulation hydrogenerator set end voltage further；Need The parameter vector wanting identification isWherein k_A,τ_AIt is respectively amplifying unit Gain and time constant, k_E,τ_EIt is respectively exciter gain and time constant, k_G,τ_GIt is respectively water wheels to send out Motor gain and time constant, k_S,τ_SIt is respectively sensor gain and time constant；

Step (2): gather actual excitation system dynamic process data: actual excitation system is carried out electricity Pressure step disturbance test, gathers actual excitation system dynamic process data, and dynamic process data include respectively Link exports, i.e. amplifying unit output V_R, exciter output V_F, hydrogenerator set end voltage V_tWith Sensor output V_S；

Step (3): set up Excitation System Parameter Identification of Synchronous object function, uses real system output and distinguishes The weighted error quadratic sum of knowledge system output defines such as the object function of parameter identification, object function Under:

${\mathrm{minf}}_{WMSE}\left(\widehat{\mathrm{\&theta;}}\right)=\underset{k=1}{\overset{N_s}{\&Sigma;}}\underset{j=1}{\overset{s}{\&Sigma;}}{w}_j{\left(y_j\right(k)-{\hat{y}}_j(k\left)\right)}^2$

Wherein N_sExporting sampling number for system, s is that system exports number, w_jFor weight,It is parameter to be identified, y_j(k)∈{V_R(k),V_F(k),V_t(k),V_S(k) } it is actual System exports；Export for identification system；

Step (4): use object function in heuristic value solution procedure (3), it is thus achieved that treat Identification system parameter；Wherein said step (3) specifically includes following sub-step:

Step 1: algorithm initialization: arrange algorithm parameter, including population size N, maximum iteration time T, individual random search quantity N_l, eliminate range coefficient σ, skip threshold p；Determine that search is to be identified The span of parameter is [B_L,B_U], concrete k_A∈[k_A,min,k_A,max], τ_A∈[τ_A,min,τ_A,max], k_E∈[k_E,min,k_E,max], τ_E∈[τ_E,min,τ_E,max]k_G∈[k_G,min,k_G,max], τ_G∈[τ_G,min,τ_G,max], k_S∈[k_S,min,k_S,max], τ_S∈[τ_S,min,τ_S,max], i.e. B_L=[k_A,min,τ_A,min,k_E,min,τ_E,min,k_G,min,τ_G,min,k_S,min,τ_S,min] represent excitation system parameter to be identified Little value, B_U=[k_A,max,τ_A,max,k_E,max,τ_E,max,k_G,max,τ_G,max,k_S,max,τ_S,max] represent excitation system ginseng to be identified The maximum of number；At solution space [B_L,B_UThe position vector of all individualities in random initializtion colony in], One group position vector is expressed as X_i=[k_A,τ_A,k_E,τ_E,k_G,τ_G,k_S,τ_S]；Make current iteration number of times t=0；

Step 2: calculate the target function value F of each individuality_i ^t=f_WMSE(X_i(t)), i=1 ..., N, and find Target population function minimum, the individuality with minimum target functional value is defined as current optimum individual X_B(t)；

Step 3: to all individual X_i(t), i=1 ..., N carries out individual random search, updates X_i(t)；

Step 4: update individual position vector X according to individual location updating formula_i(t), i=1 ..., N:

$\left\{\begin{array}{c}{\mathrm{\&delta;}}_i=|c_1\&CenterDot;X_B(t)-X_i(t\left)\right|\\ X_i(t+1)=X_B\left(t\right)+c_2\&CenterDot;{\mathrm{\&delta;}}_i\end{array}\right.$

δ_iDistance vector with current optimum individual individual for middle i-th, random number c₁=2 rand, c₂=(2 rand-1) exp (-10 t/T), rand are random number between (0,1)；

Step 5: judge individual the need of being eliminated and reinitializing:

Step 5.1: if i-th individuality meets formula, this individuality is eliminated and reinitializes:

$F_i^t>F_{ave}^t+\mathrm{\&omega;}\&CenterDot;(F_{ave}^t-F_{min}^t),i=1,\mathrm{...},N$

Wherein,It is the t meansigma methods for population all individual goals functional value,It it is minimum mesh Offer of tender numerical value, ω be one with iterations the parameter of linear increment,Value model Enclose for [-σ, σ]；

Step 5.2: the individual initialization being eliminated:

X_i=rand (1, D) × (B_U-B_L)+B_L

Wherein, D is position vector dimension；

Step 6: judge whether that continuous p is not moved for current optimum individual position, if it is, Think population extinction, the population that inverting reconstruct is new according to the following formula:

$X_i=X_B+rand\&times;\frac{R^2}{{\mathrm{\&delta;}}_i},i=1,2,\mathrm{...},N$

Wherein R is radius of inversion, R=0.1 | | B_U-B_L||；Rand is random number between (0,1), p For skip threshold；

Step 7:t=t+1, if t > T, algorithm terminates, and exports current optimum individual position as whole solution, Current optimum individual position is the parameters of excitation system that identification draws；Otherwise, Step 2 is proceeded to.

2. the method for claim 1, it is characterised in that described step Step 3 specifically includes Following sub-step:

The individual searching times l=0 of Step 3.1: order；

Step 3.2: look around a positionCalculateI=1 ..., N；

$X_i^{play}\left(t\right)=X_i\left(t\right)+rand\&CenterDot;{\mathrm{\&epsiv;}}_{play}$

Rand is random number between (0,1), ε_playFor looking around step-length；

Step 3.3: calculate next current location

$\left\{\begin{array}{ccc}{X}_i^{self}\left(t\right)=X_i\left(t\right)+rand\&CenterDot;\frac{X_i^{play}\left(t\right)-X_i\left(t\right)}{\left|\right|X_i^{play}\left(t\right)-X_i\left(t\right)\left|\right|}\&CenterDot;{\mathrm{\&epsiv;}}_{step}& if& f\left(X_i^{play}\right(t\left)\right)<f\left(X_i\right(t\left)\right)\\ X_i^{self}\left(t\right)=X_i\left(t\right)& if& f\left(X_i^{play}\right(t\left)\right)\&GreaterEqual;f\left(X_i\right(t\left)\right)\end{array}\right.$

Rand is random number between (0,1), ε_stepFor inertia step-length；

Step 3.4:l=l+1, if l is ＜ N_l, go to Step 3.2；Otherwise, Step 4 is gone to.

3. method as claimed in claim 2, it is characterised in that in Step 3.2 ε_play=0.1 | | B_U-B_L||。

4. method as claimed in claim 2, it is characterised in that in Step 3.3 ε_step=0.2 | | B_U-B_L||。

5. method as claimed in claim 2, it is characterised in that D=8 in Step 5.2.