CN104133989B - Meter and the wind power plant sequential export power calculation algorithms of icing loss - Google Patents

Abstract

The invention discloses under the situation that meter and icing lose, the computational methods of mountain area complexity landform leeward electric field power output consider the influence of Wind turbines and wind power plant inside cable and/or trolley line to power output while calculating；The result of gained is intended to the assessment of the power system Generation System Reliability accessed for wind power plant；The multimode probabilistic model that wind power plant is set up in cluster analysis is done to sequential export power, the reliability assessment of generating and transmitting system can be further used for.The method perfect influence factor of Power Output for Wind Power Field, particularly wind power plant icing congeal influence of the loss to Power Output for Wind Power Field, simulation ice-coating congeal loss result be actually consistent, the Evaluation accuracy of wind power plant reliability can be improved.

Description

Meter and the wind power plant sequential export power calculation algorithms of icing loss

Technical field

The invention belongs to field of power, more particularly to a kind of meter and the wind power plant sequential export power meter of icing loss Calculation method.

Background technology

As wind-powered electricity generation large-scale cluster formula is accessed, certain influence is also brought to power system.Therefore need to be to containing big rule The power system of mould wind power plant carries out reliability assessment.With being actually incorporated into the power networks for wind-powered electricity generation, many problems also begin to highlight Come.In February, 2011 to April causes Jiuquan, Guazhou County and body in Zhangjiakou Area, Hebei Province to recur three big rule due to cable fault The Wind turbines off-grid of mould, harsh weather can also affect greatly to the power output of wind power plant, such as Southwestern China portion area Annual December to Second Year March between, in the presence of the external elements such as temperature, humidity, circuit in wind power plant and Wind turbines can produce icing, cause that Wind turbines fault rate is high, power output weakens.

At present in the reliability assessment of electricity generation system, the model of wind power plant is defeated using multimode probabilistic model and sequential Go out power module.But do not counted in existing numerous models and above-mentioned factor influence.

By considering the factors such as icing loss, cable fault, analogue simulation obtains the power output mould of wind power plant to the present invention Type, temporal model can be used for the assessment of the reliability of electricity generation system, it is also possible to obtained after output power model is clustered The multimode probabilistic model of wind power plant, for generating and transmitting system reliability assessment.

For the use of theoretical icing model, major problem is that be input into, due to being influenceed by icing, the wind speed for measuring Data are not necessarily accurate.Additionally, the beginning and end in order to detect nature icing, in addition it is also necessary to which accurate air themperature and cloud base are high Angle value.And the acquisition of these parameters is very difficult, it may be said that icing theoretical modeling additionally depends on the development of this system with prediction. And icing model is often presented certain periodicity.A part be due to temperature cycles rise and fall round the clock, additionally, due to Disengaging caused by dynamics.So far, the detailed modeling of the process is also very limited, therefore the present invention is damaged for this icing Lose, after service data and icing loss data is obtained, the relation of power output and equivalent wind speed is built from external characteristics.Improve wind The influence factor of electric field power output.

Output power curve is obtained, the timing curve data can be used for the reliability assessment of temporal model, during one side pair The data of overture line are clustered, and can more be met actual and accurate multimode probabilistic model, are further used for hair transmission of electricity The assessment of system.

The content of the invention

In order to analyze influence of the icing loss to wind power plant sequential export power, the present invention proposes a kind of meter and icing is damaged The wind power plant sequential export power calculation algorithms of mistake, comprise the following steps：

Step 1, the data for gathering each Wind turbines of wind power plant, including：Wind turbines quantity N_WT, each height above sea levelWind speed { v in the T of wind power plant time period t=1,2 ...₁,v₂,…v_t,…,v_TAnd wind direction { d (v₁),d (v₂),…d(v_t),…,d(v_T)}；

Collection Wind turbines inherent parameters, including：Wind wheel radius r₀, sweeping area A_rotor, incision wind speed v_in, rated wind speed v_rated, cut-out wind speed v_out, rated power P_rated, hub height H；

Collection wind power plant inner member fault rate and repair rate, including：Wind turbines fault rate λ_WT, Wind turbines repair rate μ_WT, cable fault rate λ_CA, cable repair rate μ_CA；

Statistics during collection icing loss in each time period, including：Time period T during icing_iceInterior t_ice The actual wind speed of the i-th Fans of icing hourBlower fan power output under corresponding wind speedT_iceThe individual time period The loss electricity that congeals of interior i-th Fans

Step 2, by wind power plant time period t_ice=1 ... T_iceInterior wind speedStatistics, calculate Equivalent wind speed is calculated as follows during icing loses：

It is t_iceThe actual wind speed of the i-th Fans in the individual icing time period,It is the blower fan under corresponding wind speed Power output,It is t_iceThe loss electricity that congeals of the i-th Fans in the individual time period,It is t_iceThe individual time period The equivalent wind speed of interior i-th Fans,It is the blower fan power output under corresponding equivalent wind speed；Statistical equivalent wind speed probability Distribution, is simulated with piecewise function；

Step 3, carry out simulation calculation；The sequential Wind speed model of wind power plant is obtained with Monte Carlo sampling, it is considered to During icing loses, corresponding wind speed is modified, obtain sequential equivalent wind speed model；Obtain each with the sampling of Taka sieve is covered The time sequence status of Wind turbines and each cable and/or trolley line；Mountain area complexity relief model is set up, variant height wind-powered electricity generation is tried to achieve Wind speed at unit；The wake effect model of multi-overlapped is set up, the wind speed at each Wind turbines is corrected；With reference to Wind turbines The power characteristic of time sequence status and Wind turbines, tries to achieve the power output of each Wind turbines；With reference to cable and/or trolley line Time sequence status and wind power plant internal electric wiring, try to achieve the sequential export power of wind power plant.

Sequential Wind speed model uses Weibull fittings of distribution in the step 3, and Weibull probability distribution is as follows：

K is form parameter, and c is scale parameter, and v is wind speed；Form parameter and scale parameter are obtained by actual-structure measurement Go out；Using maximum likelihood estimate, it is assumed that wind farm wind velocity statistical sample is v={ v₁,v₂,v₃,…,v_n, n represents statistics sample This number, the then specific calculation expression that form parameter k and scale parameter c are solved using maximum likelihood estimate is as follows：

Form parameter k is solved using iterative method in formula, so as to try to achieve dimensional parameters c.

The mountain area complexity relief model is the place installed based on Lissaman wake effect models, different Wind turbines Difference, with different height above sea levels, and wind speed with altitude changes and changes, and causes the wind speed profile of wind power plant uneven, its change Change situation is represented with following formula：

In formula：v₀It is being highly h to be₀The wind speed for measuring, m/s；V (h) is the wind speed measured for h in height, m/s；α is wind Speed typically takes α=1/7 with height change factor；Thus obtain wind speed and increase as height increases；By lossless bernoulli side Cheng get, the wind speed decreased coefficient d of level terrain_FWith the wind speed decreased coefficient d of complicated landform_CRelation is as follows：

Then wind speed everywhere is：

Wherein, C_TIt is thrust coefficient, r is the radius influenceed by wake effect.

The wake effect model of the multi-overlapped is based on Lissaman wake effect models, it is assumed that be apart x_ij's WT_iAnd WT_jHeight above sea level be respectively h_i、h_j, then downstream Wind turbines WT_jBy upstream Wind turbines WT_iWind speed after influence Computing formula is as follows：

In formula：A_rotorIt is wind wheel sweeping area,A_shad.ijIt is the area of lap, r_iAnd r_jRespectively Upstream Wind turbines WT_iWith downstream Wind turbines WT_jWind wheel radius；Computing formula is tried to achieve by following formula：

Δ h=in formula | h_j-h_i|, it is that the height above sea level of two Wind turbines is poor；d_ijIt is upstream Wind turbines WT_iWith lower urticaria Group of motors WT_jBetween level interval；Wake effect is not paid attention in the situation of sowing：The Wind turbines failure of upstream, under Urticaria group of motors has no influence；Wind turbines power characteristic is as follows：

V, v_in, v_rated, v_outRespectively currently practical wind speed, incision wind speed, rated wind speed and cut-out wind speed, P_ratedIt is wind The rated power of machine.

Time sequence status { the S of Wind turbines is obtained in the step 3_WT.1,S_WT.2,…S_WT.t... .S_WT.T, wherein, some time The Wind turbines state at quarterS_WTi.t(i=1 ... N_WT, t=1,2, 3 ... T) states of the Wind turbines i in t is represented, its value is only 0 and 1, blower fan is represented respectively and is in stoppage in transit and operation shape State, N_WTIt is blower fan sum, T is simulation total time；Obtain the time sequence status { S of cable and/or trolley line_CL.1,S_CL.2,… S_CL.t,…S_CL.T, wherein, certain moment cable and/or trolley line state S_CLm.t(m=1,2,3 ... N_CL, t=1 ... T) and the state of cable and/or trolley line m in t is represented, its value is only 0 and 1, Cable is represented respectively and/or trolley line is in and stops transport and running status, N_CLIt is cable and/or trolley line number, when T is for simulation Between；

Bring the equivalent wind speed in step 2 into the P that exerts oneself that formula (9) tries to achieve t Wind turbines i_WTi.t, then try to achieve and combine wind It is P that the t Wind turbines i of group of motors state exerts oneself_WTi.tWith S_WTi.tProduct：P_WTi.t·S_WTi.t；

Try to achieve the t Wind turbines i with reference to Wind turbines state and cable and/or trolley line state exert oneself forSo this branch road power output sum isWind power plant Comprising a plurality of chain type branch road, after trying to achieve every branch road sum, addition obtains the output of wind power plant t.

The idiographic flow of the simulation calculation in the step 3, including：

1) the location of each Wind turbines of wind power plant coordinate is read in；

2) simulation time initialization, t is hourage, since t=1；

3) during judging whether to enter icing；According to the characteristics of ice-covering area, during 2 months December to next years were icing, meter Simulated time is calculated when 1 day 1 January, during 1- is icing in the 8760th hour the 1416th hour and 8017-, and with this It is criterion；

If 4) into during icing, using equivalent wind speed modeling wind speed；Otherwise, using Weibull distribution simulations Wind speed；

5) wind speed and direction of t hours is read in；

6) Monte Carlo sampling obtains the Wind turbines that normal operation and failure are stopped transport；

7) wind speed at each Wind turbines is calculated according to wake effect correlation formula (4)-(8)；

8) output of each Wind turbines actual power is calculated according to Wind turbines power characteristic formula (9)；

9) simulation time is set Y=8760 hours herein, i.e., 1 year, hourage t=t+1 judged whether simulation time is tied Beam, t>Y then terminates to calculate, t<Y is then transferred to calculation process 3).

It is perfect the beneficial effects of the present invention are there is provided a kind of method for setting up wind power plant Wind turbines icing model The influence factor of Power Output for Wind Power Field, particularly wind power plant icing are congealed influence of the loss to Power Output for Wind Power Field, simulation Icing congeal loss result be actually consistent, the Evaluation accuracy of wind power plant reliability can be improved.

Brief description of the drawings

Fig. 1 is the flow of simulation calculation；

Fig. 2 is wind farm wind velocity probability distribution statistical and fitted figure, and column curve is the distribution of actual count wind speed, is smoothed Curve is Weibull Probability Distribution Fittings；

The wind direction statistics of Fig. 3 wind power plants, rose circle diagram, N in figure, NNE, NE, ENE, E, ESE, SE, SSE, S, SSW, SW, WSW, W, WNW, NW, NNW, represent respectively north, northeast by north, northeast, northeast by east, east, southeast by east, the southeast, southeast by south, South, swbs, southwest, southwest by west, west, northwest by west, northwest, north-west by north；

Fig. 4 is the probability distribution of equivalent wind speed during icing as obtained by calculating loses, it can be seen that Weibull is distributed Cannot well be fitted, then be simulated using piecewise function；

Fig. 5 is the power curve of Wind turbines；

Fig. 6 is that wake effect model, v are blocked in complicated landform bottom point₀It is original wind speed, WT_iUpwindturbine, WT_jDownwindturbine is respectively upstream Wind turbines and the downstream Wind turbines (x influenceed by wake effect_i, y_i, z_i) (x_j, y_j, z_j) it is respectively its coordinate, d_ij, x_ijFor not Wei two Wind turbines horizontal axis away from, anterior-posterior horizontal distance, v_TBe by WT_iAfterwards with the Wind turbines at wind speed, v_j(x_ij) it is downstream Wind turbines WT_jBy upstream Wind turbines WT_iWind after influence Speed, A₀It is sweeping area, A_{X, j}It is Wind turbines WT_iIn Wind turbines WT_jGo out wake effect influence area, r_i(y_ij) it is the area Radius, A_shadowIt is preceding two-part lap area, angle, θ is the corresponding central angle of overlapping area, r₀It is wind wheel half Footpath；

Fig. 7 is corresponding multiple wake effect model, v₀It is original wind speed, 1#, 2# etc. are respectively Wind turbines numbering；

Fig. 8 is wind power plant Wind turbines layout, and transverse and longitudinal axle difference denotation coordination, black side's point represents that Wind turbines are installed Place, the height of Wind turbines is labeled in Wind turbines installation place；

Fig. 9 is cabling diagram；

Figure 10 is the sequential export power curve of wind power plant；

Figure 11 is one cabling diagram of branch road in wind power plant inside, and bus represents bus, and arrow represents power flow direction, WT Wind turbines are represented, WT1 ... WT9, WT10 represent Wind turbines numbering, and CL represents cable or trolley line CL1, CL10 representative Cable number.

Specific embodiment

Below in conjunction with the accompanying drawings, preferred embodiment is elaborated.It should be emphasized that the description below is merely exemplary , rather than in order to limit the scope of the present invention and its application.

The present invention proposes the wind power plant sequential export power calculation algorithms of a kind of meter and icing loss, including following step Suddenly：

Step 1, the data for gathering each Wind turbines of wind power plant, including：Wind turbines quantity N_WT, each height above sea levelWind speed { v in the T of wind power plant time period t=1,2 ...₁,v₂,…v_t,…,v_TAnd wind direction { d (v₁),d (v₂),…d(v_t),…,d(v_T)}；The wind direction statistics of wind power plant is illustrated in figure 3, wind direction is angularly divided into 16 directions, presents Go out rose circle diagram；

Collection Wind turbines inherent parameters, obtain under the parameter of Wind turbines：

Wind wheel radius r₀=45m

Sweeping area A_rotor=6232m²

Incision wind speed v_in=3m/s

Rated wind speed v_rated=12m/s

Cut-out wind speed v_out=25m/s

Rated power P_rated=2MW

Hub height H=60m

Wind power plant inner member fault rate and repair rate：Wind turbines fault rate λ_WT=0.012 times/year, Wind turbines are repaiied Multiple rate μ_WT=30 days；Cable fault rate λ_CA=0.008 times/year, cable repair rate μ_CA=12 days.

Statistics during collection icing loss in each time period, including：Time period T during icing_iceInterior t_ice The actual wind speed of the i-th Fans of icing hourBlower fan power output under corresponding wind speedT_iceThe individual time The loss electricity that congeals of the i-th Fans in section

Step 2, by wind power plant time period t_ice=1 ... T_iceInterior wind speedStatistics, calculate Equivalent wind speed is calculated as follows during icing loses：

It is t_iceThe actual wind speed of the i-th Fans in the individual icing time period,It is the blower fan under corresponding wind speed Power output,It is t_iceThe loss electricity that congeals of the i-th Fans in the individual time period,It is t_iceIn the individual time period The equivalent wind speed of the i-th Fans,It is the blower fan power output under corresponding equivalent wind speed；Statistical equivalent wind speed probability point Cloth, is simulated with piecewise function；

Step 3, the sequential Wind speed model that wind power plant is obtained with Monte Carlo sampling, it is considered to during icing loses, to phase The wind speed answered is modified, and obtains sequential equivalent wind speed model；Each Wind turbines and each cable are obtained with the sampling of Taka sieve is covered And/or the time sequence status of trolley line；Mountain area complexity relief model is set up, the wind speed at variant height Wind turbines is tried to achieve；Build The wake effect model of vertical multi-overlapped, corrects the wind speed at each Wind turbines；With reference to the time sequence status and wind-powered electricity generation of Wind turbines The power characteristic of unit, tries to achieve the power output of each Wind turbines；With reference to cable and/or the time sequence status and wind of trolley line Electric field internal electric wiring, tries to achieve the sequential export power of wind power plant.

Wherein, sequential Wind speed model uses Weibull fittings of distribution in step 3, is as shown in Figure 2 wind farm wind velocity probability Distribution statisticses and fitting, column curve are the distribution of actual count wind speed, and smoothed curve is Weibull Probability Distribution Fittings.

Weibull probability distribution is as follows：

K is form parameter, and c is scale parameter, and v is wind speed；Form parameter and scale parameter are obtained by actual-structure measurement Go out.Using maximum likelihood estimate, it is assumed that wind farm wind velocity statistical sample is v={ v₁,v₂,v₃,…,v_n, n represents statistics sample This number, the then specific calculation expression that form parameter k and scale parameter c are solved using maximum likelihood estimate is as follows：

Form parameter k can be solved using iterative method in formula, so as to try to achieve dimensional parameters c.

Equivalent wind speed is changed with piecewise function during icing, is calculated as follows：

It is t_iceThe actual wind speed of the i-th Fans in the individual icing time period,It is the blower fan under corresponding wind speed Power output,It is t_iceThe loss electricity that congeals of the i-th Fans in the individual time period,It is t_iceIn the individual time period The equivalent wind speed of the i-th Fans,It is the blower fan power output under corresponding equivalent wind speed.

Mountain area complexity relief model is that the place that different Wind turbines are installed is not based on Lissaman wake effect models Together, with different height above sea levels, and wind speed with altitude changes and changes, and causes the wind speed profile of wind power plant uneven.Specific mould Type is shown in Fig. 6, v₀It is original wind speed, WT_iUpwindturbine, WT_jDownwind turbine be respectively upstream Wind turbines and Downstream Wind turbines (the x influenceed by wake effect_i, y_i, z_i)(x_j, y_j, z_j) it is respectively its coordinate, d_ij, x_ijFor Wei not two wind-powered electricity generations The horizontal axis of unit away from, anterior-posterior horizontal distance, v_TIt is by WT_iAfterwards with the Wind turbines at wind speed, v_j(x_ij) it is lower urticaria Group of motors WT_jBy upstream Wind turbines WT_iWind speed after influence, A₀It is sweeping area, A_{X, j}It is Wind turbines WT_iIn wind turbine Group WT_jGo out wake effect influence area, r_i(y_ij) be the area radius, A_shadowIt is preceding two-part lap area, angle Degree θ is the corresponding central angle of overlapping area, r₀It is wind wheel radius.Its situation of change can have following formula to represent：

In formula：v₀It is being highly h to be₀The wind speed for measuring, m/s；V (h) is the wind speed measured for h in height, m/s；α is wind Speed typically takes α=1/7 with height change factor.Thus obtain wind speed and increase as height increases.By lossless bernoulli side Journey can be obtained, the wind speed decreased coefficient d of level terrain_FWith the wind speed decreased coefficient d of complicated landform_CRelation is as follows：

Then wind speed everywhere is：

Wherein, C_TIt is thrust coefficient, r is the radius influenceed by wake effect.

It, based on Lissaman wake effect models, is as shown in Figure 7 corresponding many that the wake effect model of multi-overlapped is Heavy-tailed flow effect model, v₀It is original wind speed, 1#, 2# etc. are respectively Wind turbines numbering.Assuming that being apart x_ijWT_iAnd WT_j's Height above sea level is respectively h_i、h_j, then downstream Wind turbines WT_jBy upstream Wind turbines WT_iWind speed computing formula after influence is such as Under：

In formula：A_rotorIt is wind wheel sweeping area,A_shad.ijIt is the area of lap, r_iAnd r_jRespectively Upstream Wind turbines WT_iWith downstream Wind turbines WT_jWind wheel radius；Computing formula can be tried to achieve by following formula：

Δ h=in formula | h_j-h_i|, it is that the height above sea level of two Wind turbines is poor, d_ijIt is upstream Wind turbines WT_iWith lower urticaria Group of motors WT_jBetween level interval.Wake effect is not paid attention in the situation of sowing：The Wind turbines failure of upstream, under Urticaria group of motors has no influence.Wind turbines power characteristic is as follows：

V, v_in, v_rated, v_outRespectively currently practical wind speed, incision wind speed, rated wind speed and cut-out wind speed, P_ratedIt is wind The rated power of machine.It is as shown in Figure 5 the power curve of Wind turbines, incision wind speed v_in=3m/s；Rated wind speed v_rated= 12m/s, cut-out wind speed v_out=25m/s, rated power P_rated=2MW.

Fig. 9 is cabling diagram, including power supply is Wind turbines, bus, transformer and cable connection, and connected mode is chain Formula is connected, and 10 typhoon group of motors is connected on every cable, totally 8 typhoon group of motors of cable 80, total installation of generating capacity 160MW.Wind-powered electricity generation Unit after boosting after step-up transformer by, through cable connection to bus, further accessing power network.With reference to cable and/or trolley line Time sequence status and wind power plant internal electric wiring, try to achieve the sequential export power of wind power plant, be as shown in Figure 10 wind power plant Sequential export power curve.

Time sequence status { the S of resulting Wind turbines_WT.1,S_WT.2,…S_WT.t... .S_WT.T, wherein, the wind turbine at certain moment Group stateS_WTi.t(i=1 ... N_WT, t=1,2,3 ... T) and represent wind In the state of t, its value is only 0 and 1 to group of motors i, blower fan is represented respectively is in and stop transport and running status, N_WTFor blower fan is total Number, T is simulation total time；Obtain the time sequence status { S of cable and/or trolley line_CL.1,S_CL.2,…S_CL.t,…S_CL.T, wherein, Certain moment cable and/or trolley line stateS_CLm.t(m=1,2,3 ... N_CL, t=1 ... T) represent the state of cable and/or trolley line m in t, its value is only 0 and 1, represent respectively cable and/or Trolley line is in stops transport and running status, N_CLIt is cable and/or trolley line number, T is simulated time.

With reference to the time sequence status and the power characteristic of Wind turbines of Wind turbines, the output work of each Wind turbines is tried to achieve Rate.Bring the equivalent wind speed in step 2 into the P that exerts oneself that formula (9) tries to achieve t Wind turbines i_WTi.t, with reference to the shape of Wind turbines State, it is P that final t Wind turbines i exerts oneself_WTi.tWith S_WTi.tProduct：P_WTi.t·S_WTi.t。

With reference to the time sequence status and wind power plant internal electric wiring of cable and/or trolley line, the sequential for trying to achieve wind power plant is defeated Go out power.Internal electric connection figure is shown in Fig. 9 and Figure 11

By taking Figure 11 as an example, work as t, Wind turbines i exert oneself forSo this branch road Power output sum isWind power plant includes a plurality of chain type branch road, tries to achieve every branch road sum Afterwards, it is added the output for obtaining wind power plant t.

The idiographic flow of the simulation calculation in step 3 as shown in figure 1, including：

1) the location of each Wind turbines of wind power plant coordinate is read in, as shown in figure 8, transverse and longitudinal axle represents seat respectively Mark, black side's point represents Wind turbines installation place, and the height of Wind turbines is labeled in Wind turbines installation place；

2) simulation time initialization, t is hourage, since t=1；

3) during judging whether to enter icing.According to the characteristics of ice-covering area, during 2 months December to next years were icing.Meter Simulated time is calculated when 1 day 1 January, during 1- is icing in the 8760th hour the 1416th hour and 8017-, and with this It is criterion；

If 4) into during icing, using equivalent wind speed modeling wind speed；Otherwise, using Weibull distribution simulations Wind speed.The asking for of equivalent wind speed sees formula (3), and the probability distribution graph of equivalent wind speed is shown in Fig. 4, it can be seen that Weibull distributions cannot Fitting well, is then simulated using piecewise function.

5) wind speed and direction of t hours is read in；

6) Monte Carlo sampling obtains the Wind turbines that normal operation and failure are stopped transport；

7) wind speed at each Wind turbines is calculated according to wake effect correlation formula (4)-(8)；Single wind energy conversion system wake flow Computation model-Lissaman models

8) output of each Wind turbines actual power is calculated according to Wind turbines power characteristic formula (9)；

9) simulation time is set Y=8760 hours, i.e., 1 year herein.Hourage t=t+1, judges whether simulation time is tied Beam, t>Y then terminates to calculate, t<Y is then transferred to calculation process 3).

The above, the only present invention preferably specific embodiment, but protection scope of the present invention is not limited thereto, Any one skilled in the art the invention discloses technical scope in, the change or replacement that can be readily occurred in, Should all be included within the scope of the present invention.Therefore, protection scope of the present invention should be with scope of the claims It is defined.

Claims (6)

1. the wind power plant sequential export power calculation algorithms that a kind of meter and icing lose, it is characterised in that comprise the following steps：

Step 1, the data for gathering each Wind turbines of wind power plant, including：Wind turbines quantity N_WT, each height above sea levelWind speed { v in the T of wind power plant time period t=1,2 ...₁,v₂,…v_t,…,v_TAnd wind direction { d (v₁),d (v₂),…d(v_t),…,d(v_T)}；

Collection Wind turbines inherent parameters, including：Wind wheel radius r₀, sweeping area A_rotor, incision wind speed v_in, rated wind speed v_rated, cut-out wind speed v_out, rated power P_rated, hub height H；

Collection wind power plant inner member fault rate and repair rate, including：Wind turbines fault rate λ_WT, Wind turbines repair rate μ_WT, Cable fault rate λ_CA, cable repair rate μ_CA；

Statistics during collection icing loss in each time period, including：Time period T during icing_iceInterior t_iceIcing The actual wind speed of the i-th Fans of hourBlower fan power output under corresponding wind speedT_iceIn the individual time period The loss electricity that congeals of i Fans

Step 2, by wind power plant time period t_ice=1 ... T_iceInterior wind speedStatistics, calculate icing Equivalent wind speed is calculated as follows during loss：

$\underset{t_{ice}=1}{\overset{T_{ice}}{\&Sigma;}}\underset{i=1}{\overset{N}{\&Sigma;}}P\left(v_{i.t_{ice}}\right)-\underset{t_{ice}=1}{\overset{T_{ice}}{\&Sigma;}}\underset{i=1}{\overset{N}{\&Sigma;}}{P}_{loss.i.t_{ice}}=\underset{t_{ice}=1}{\overset{T_{ice}}{\&Sigma;}}\underset{i=1}{\overset{N}{\&Sigma;}}P\left(V_{eq.i.t_{ice}}\right)---\left(3\right)$

It is t_iceThe actual wind speed of the i-th Fans in the individual icing time period,It is the blower fan output under corresponding wind speed Power,It is t_iceThe loss electricity that congeals of the i-th Fans in the individual time period,It is t_iceI-th in the individual time period The equivalent wind speed of Fans,It is the blower fan power output under corresponding equivalent wind speed；Statistical equivalent wind velocity distributing paremeter, It is simulated with piecewise function；

Step 3, carry out simulation calculation；The sequential Wind speed model of wind power plant is obtained with Monte Carlo sampling, it is considered to icing During loss, corresponding wind speed is modified, obtain sequential equivalent wind speed model；Each wind-powered electricity generation is obtained with Monte Carlo sampling The time sequence status of unit and each cable and/or trolley line；Mountain area complexity relief model is set up, variant height Wind turbines are tried to achieve The wind speed at place；The wake effect model of multi-overlapped is set up, the wind speed at each Wind turbines is corrected；With reference to the sequential of Wind turbines The power characteristic of state and Wind turbines, tries to achieve the power output of each Wind turbines；With reference to cable and/or trolley line when Sequence state and wind power plant internal electric wiring, try to achieve the sequential export power of wind power plant.

2. method according to claim 1, it is characterised in that sequential Wind speed model is using Weibull points in the step 3 Cloth is fitted, and Weibull probability distribution is as follows：

$f\left(v\right)=\left(\frac{k}{c}\right){\left(\frac{v}{c}\right)}^{k-1}\exp \&lsqb;-{\left(\frac{v}{c}\right)}^k\&rsqb;---\left(1\right)$

K is form parameter, and c is scale parameter, and v is wind speed；Form parameter and scale parameter draw by actual-structure measurement；Adopt With maximum likelihood estimate, it is assumed that wind farm wind velocity statistical sample is v={ v₁,v₂,v₃,…,v_n, n represents statistical sample Number, the then specific calculation expression that form parameter k and scale parameter c are solved using maximum likelihood estimate is as follows：

$\left\{\begin{array}{c}k={(\frac{\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{v}_i^k\ln v_i}{\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{v}_i^k}+\frac{1}{n}\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}\ln v_i)}^{-1}\\ c={\left(\frac{1}{n}\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{\left(v_i\right)}^k\right)}^{1/k}\end{array}\right.---\left(2\right)$

Form parameter k is solved using iterative method in formula, so as to try to achieve dimensional parameters c.

3. method according to claim 1, it is characterised in that the mountain area complexity relief model is based on Lissaman tails Stream effect model, the place that different Wind turbines are installed is different, and with different height above sea levels, and wind speed with altitude changes and becomes Change, cause the wind speed profile inequality of wind power plant, its situation of change to be represented with following formula：

$\frac{v\left(h\right)}{v_0}={\left(\frac{h}{h_0}\right)}^{\mathrm{\&alpha;}}---\left(4\right)$

In formula：v₀It is being highly h to be₀The wind speed for measuring, m/s；V (h) is the wind speed measured for h in height, m/s；α be wind speed with Height change factor, takes α=1/7；Thus obtain wind speed and increase as height increases；Obtained by lossless D.Bernolli equation, put down The wind speed decreased coefficient d of smooth landform_FWith the wind speed decreased coefficient d of complicated landform_CRelation is as follows：

$\frac{d_C}{d_F}={\left(\frac{v_0}{v}\right)}^2---\left(5\right)$

Then wind speed everywhere is：

$\begin{array}{c}v(x,h)=v\left(h\right)(1-d_C)\\ =v_0\&lsqb;1-(1-\sqrt{1-C_T}){\left(\frac{h_0}{h}\right)}^{2\mathrm{\&alpha;}}{\left(\frac{r_0}{r}\right)}^2\&rsqb;{\left(\frac{h}{h_0}\right)}^{\mathrm{\&alpha;}}\end{array}---\left(6\right)$

Wherein, C_TIt is thrust coefficient, r is the radius influenceed by wake effect.

4. method according to claim 1, it is characterised in that the wake effect model of the multi-overlapped be based on Lissaman wake effect models, it is assumed that be apart x_ijWT_iAnd WT_jHeight above sea level be respectively h_i、h_j, then downstream wind turbine Group WT_jBy upstream Wind turbines WT_iWind speed computing formula after influence is as follows：

$\begin{array}{c}{v}_j(x_{ij},h_j)=v_i\&lsqb;1-(1-\sqrt{1-C_T})\&CenterDot;{\left(\frac{h_i}{h_j}\right)}^{2\mathrm{\&alpha;}}\\ \&CenterDot;{\left(\frac{r_j}{r_i\left(x_{ij}\right)}\right)}^2\&CenterDot;\left(\frac{A_{shad.ij}}{A_{rotor}}\right)\&rsqb;{\left(\frac{h_j}{h_i}\right)}^{\mathrm{\&alpha;}}\end{array}---\left(7\right)$

In formula：A_rotorIt is wind wheel sweeping area, A_rotor=π r_j ²；A_shad.ijIt is the area of lap, r_iAnd r_jRespectively upstream Wind turbines WT_iWith downstream Wind turbines WT_jWind wheel radius；Computing formula is tried to achieve by following formula：

$\begin{array}{c}{A}_{shad.ij}=r_i^2\left(x_{ij}\right){\cos }^{-1}\left(\frac{r_i^2\left(x_{ij}\right)+d_{ij}^2+{\mathrm{\&Delta;h}}^2-r_j^2}{2r_i\left(x_{ij}\right)\&CenterDot;\sqrt{{\mathrm{\&Delta;h}}^2+d_{ij}^2}}\right)\\ +r_j^2{\cos }^{-1}\left(\frac{r_j^2+d_{ij}^2+{\mathrm{\&Delta;h}}^2-r_i^2\left(x_{ij}\right)}{2r_j\&CenterDot;\sqrt{{\mathrm{\&Delta;h}}^2+d_{ij}^2}}\right)\\ -r_i\left(x_{ij}\right)\&CenterDot;\sqrt{{\mathrm{\&Delta;h}}^2+d_{ij}^2}\\ \&CenterDot;\sin \&lsqb;{\cos }^{-1}\left(\frac{r_i^2\left(x_{ij}\right)+d_{ij}^2+{\mathrm{\&Delta;h}}^2-r_j^2}{2r_i\left(x_{ij}\right)\&CenterDot;\sqrt{{\mathrm{\&Delta;h}}^2+d_{ij}^2}}\right)\&rsqb;\end{array}---\left(8\right)$

Δ h=in formula | h_j-h_i|, it is that the height above sea level of two Wind turbines is poor；d_ijIt is upstream Wind turbines WT_iWith downstream wind turbine Group WT_jBetween level interval；Wake effect is not paid attention in the situation of sowing：The Wind turbines failure of upstream, to lower urticaria Group of motors has no influence；Wind turbines power characteristic is as follows：

$P\left(v\right)=\left\{\begin{array}{cc}0,& v>v_{out}orv<v_{in}\\ P_{rated}\frac{v^3-v_{in}^3}{v_{rated}^3-v_{in}^3},& v_{in}\&le;v\&le;v_{rated}\\ P_{rated},& v_{rated}\&le;v\&le;v_{out}\end{array}\right.---\left(9\right)$

V, v_in, v_rated, v_outRespectively currently practical wind speed, incision wind speed, rated wind speed and cut-out wind speed, P_ratedIt is blower fan Rated power.

5. method according to claim 1, it is characterised in that the time sequence status of Wind turbines are obtained in the step 3 {S_WT.1,S_WT.2,…S_WT.t... .S_WT.T, wherein, the Wind turbines state at certain momentS_WTi.t(i=1 ... N_WT, t=1,2,3 ... T) and represent Wind turbines i In the state of t, its value is only 0 and 1, blower fan is represented respectively and is in stoppage in transit and running status, N_WTIt is blower fan sum, T is Simulation total time；Obtain the time sequence status { S of cable and/or trolley line_CL.1,S_CL.2,…S_CL.t,…S_CL.T, wherein, certain moment Cable and/or trolley line stateS_CLm.t(m=1,2,3 ... N_CL, t= 1 ... T) state of cable and/or trolley line m in t is represented, its value is only 0 and 1, cable and/or trolley line is represented respectively In stoppage in transit and running status, N_CLIt is cable and/or trolley line number, T is simulated time；

Equivalent wind speed in step 2 is substituted into the P that exerts oneself that formula (9) tries to achieve t Wind turbines i_WTi.t,

$P\left(v\right)=\left\{\begin{array}{cc}0,& v>v_{out}orv<v_{in}\\ P_{rated}\frac{v^3-v_{in}^3}{v_{rated}^3-v_{in}^3},& v_{in}\&le;v\&le;v_{rated}\\ P_{rated},& v_{rated}\&le;v\&le;v_{out}\end{array}\right.---\left(9\right)$

V, v_in, v_rated, v_outRespectively currently practical wind speed, incision wind speed, rated wind speed and cut-out wind speed, P_ratedIt is blower fan Rated power；

The exerting oneself for t Wind turbines i tried to achieve again with reference to Wind turbines state is P_WTi.tWith S_WTi.tProduct：P_WTi.t· S_WTi.t；

Try to achieve the t Wind turbines i with reference to Wind turbines state and cable and/or trolley line state exert oneself forSo this branch road power output sum isWind power plant Comprising a plurality of chain type branch road, after trying to achieve every branch road sum, addition obtains the output of wind power plant t.

6. method according to claim 1, it is characterised in that the idiographic flow of the simulation calculation in the step 3, Including：

1) the location of each Wind turbines of wind power plant coordinate is read in；

2) simulation time initialization, t is hourage, since t=1；

3) during judging whether to enter icing；According to the characteristics of ice-covering area, during 2 months December to next years were icing, mould is calculated Pseudotime is when 1 day 1 January, during 1- the 8760th hour the 1416th hour and 8017- is icing, and as sentencing According to；

If 4) into during icing, using equivalent wind speed modeling wind speed；Otherwise, using Weibull distribution simulation wind speed；

5) wind speed and direction of t hours is read in；

6) Monte Carlo sampling obtains the Wind turbines that normal operation and failure are stopped transport；

7) wind speed at each Wind turbines is calculated according to wake effect correlation formula (4)-(8)；

$\frac{v\left(h\right)}{v_0}={\left(\frac{h}{h_0}\right)}^{\mathrm{\&alpha;}}---\left(4\right)$

In formula：v₀It is being highly h to be₀The wind speed for measuring, m/s；V (h) is the wind speed measured for h in height, m/s；α be wind speed with Height change factor, takes α=1/7；Thus obtain wind speed and increase as height increases；Obtained by lossless D.Bernolli equation, put down The wind speed decreased coefficient d of smooth landform_FWith the wind speed decreased coefficient d of complicated landform_CRelation is as follows：

$\frac{d_C}{d_F}={\left(\frac{v_0}{v}\right)}^2---\left(5\right)$

Then wind speed everywhere is：

$\begin{array}{c}v(x,h)=v\left(h\right)(1-d_C)\\ =v_0\&lsqb;1-(1-\sqrt{1-C_T}){\left(\frac{h_0}{h}\right)}^{2\mathrm{\&alpha;}}{\left(\frac{r_0}{r}\right)}^2\&rsqb;{\left(\frac{h}{h_0}\right)}^{\mathrm{\&alpha;}}\end{array}---\left(6\right)$

Wherein, C_TIt is thrust coefficient, r is the radius influenceed by wake effect；

Assuming that being apart x_ijWT_iAnd WT_jHeight above sea level be respectively h_i、h_j, then downstream Wind turbines WT_jBy upstream wind turbine Group WT_iWind speed computing formula after influence is as follows：

$\begin{array}{c}{v}_j(x_{ij},h_j)=v_i\&lsqb;1-(1-\sqrt{1-C_T})\&CenterDot;{\left(\frac{h_i}{h_j}\right)}^{2\mathrm{\&alpha;}}\\ \&CenterDot;{\left(\frac{r_j}{r_i\left(x_{ij}\right)}\right)}^2\&CenterDot;\left(\frac{A_{shad.ij}}{A_{rotor}}\right)\&rsqb;{\left(\frac{h_j}{h_i}\right)}^{\mathrm{\&alpha;}}\end{array}---\left(7\right)$

In formula：A_rotorIt is wind wheel sweeping area, A_rotor=π r_j ²；A_shad.ijIt is the area of lap, r_iAnd r_jRespectively upstream Wind turbines WT_iWith downstream Wind turbines WT_jWind wheel radius；Computing formula is tried to achieve by following formula：

$\begin{array}{c}{A}_{shad.ij}=r_i^2\left(x_{ij}\right){\cos }^{-1}\left(\frac{r_i^2\left(x_{ij}\right)+d_{ij}^2+{\mathrm{\&Delta;h}}^2-r_j^2}{2r_i\left(x_{ij}\right)\&CenterDot;\sqrt{{\mathrm{\&Delta;h}}^2+d_{ij}^2}}\right)\\ +r_j^2{\cos }^{-1}\left(\frac{r_j^2+d_{ij}^2+{\mathrm{\&Delta;h}}^2-r_i^2\left(x_{ij}\right)}{2r_j\&CenterDot;\sqrt{{\mathrm{\&Delta;h}}^2+d_{ij}^2}}\right)\\ -r_i\left(x_{ij}\right)\&CenterDot;\sqrt{{\mathrm{\&Delta;h}}^2+d_{ij}^2}\\ \&CenterDot;\sin \&lsqb;{\cos }^{-1}\left(\frac{r_i^2\left(x_{ij}\right)+d_{ij}^2+{\mathrm{\&Delta;h}}^2-r_j^2}{2r_i\left(x_{ij}\right)\&CenterDot;\sqrt{{\mathrm{\&Delta;h}}^2+d_{ij}^2}}\right)\&rsqb;\end{array}---\left(8\right)$

8) output of each Wind turbines actual power is calculated according to Wind turbines power characteristic formula (9)；

Wind turbines power characteristic is as follows：

$P\left(v\right)=\left\{\begin{array}{cc}0,& v>v_{out}orv<v_{in}\\ P_{rated}\frac{v^3-v_{in}^3}{v_{rated}^3-v_{in}^3},& v_{in}\&le;v\&le;v_{rated}\\ P_{rated},& v_{rated}\&le;v\&le;v_{out}\end{array}\right.---\left(9\right)$

V, v_in, v_rated, v_outRespectively currently practical wind speed, incision wind speed, rated wind speed and cut-out wind speed, P_ratedIt is blower fan Rated power

9) simulation time is set Y=8760 hours herein, i.e., 1 year, hourage t=t+1 judged whether simulation time terminates, t> Y then terminates to calculate, t<Y is then transferred to calculation process 3).