CN114626575A

CN114626575A - Reactive power planning method for receiving-end power grid containing high-permeability wind power and considering transient voltage stability

Info

Publication number: CN114626575A
Application number: CN202210116388.6A
Authority: CN
Inventors: 徐艳春; 蒋伟俊; 汪平; 孙思涵
Original assignee: China Three Gorges University CTGU
Current assignee: China Three Gorges University CTGU
Priority date: 2022-02-07
Filing date: 2022-02-07
Publication date: 2022-06-14

Abstract

The reactive power planning method for the receiving-end power grid containing high-permeability wind power and considering transient voltage stability comprises the following steps: based on the criterion of a multi-binary table, a transient voltage safety margin index is provided for evaluating the transient voltage stability of the system; based on the transient voltage safety margin index, a point distribution method of the dynamic reactive power compensation equipment is provided; establishing a differential dynamic reactive power compensation optimization model; and screening out the optimal scheme under each scene by an improved entropy weight optimal solution distance method, and then determining the final configuration scheme of the dynamic reactive power compensation equipment in the system. The method can comprehensively and comprehensively guide the configuration of the dynamic reactive power compensation device, and can save the economic cost to the maximum extent while ensuring the stable operation of the power grid.

Description

Reactive power planning method for receiving-end power grid containing high-permeability wind power and considering transient voltage stability

Technical Field

The invention relates to the technical field of reactive power planning of power grids, in particular to a receiving-end power grid reactive power planning method considering transient voltage stability and containing high-permeability wind power.

Background

Wind power is used as a main body of new energy power generation, the permeability of the wind power in a power grid is increasing year by year, and meanwhile, the problem of transient voltage is also increasingly highlighted. Reactive power planning is a mainstream approach for solving the voltage stability of a power grid, but most of the reactive power planning is oriented to the static voltage problem of the power grid at present, and aims to reduce the grid loss and the voltage deviation of a system while ensuring the optimal economic investment, and a planning method focusing on the transient voltage stability is rare. Therefore, it is necessary to establish a scientific and feasible reactive power planning method to improve the robustness of the power grid after the action of the large disturbance.

At present, the following three problems generally exist in a reactive power planning method considering the transient voltage stability of a power grid:

1) the planned dynamic reactive power source has single type, and the economic efficiency can not be considered while the system can obtain the optimal reactive power compensation effect;

2) the research background is mostly limited to the traditional power grid and does not conform to the current China power grid development pattern;

3) for a receiving end power grid containing DG (distributed generation), the influence of DG uncertainty on the planning is not considered by the existing reactive power planning method.

4) Due to the lack of a standard for measuring the transient voltage stability margin in detail, most of reactive power planning only focuses on the voltage instability problem caused by voltage sag, and the overvoltage problem caused by voltage sag is rarely considered.

In view of the fact that point placement and constant volume are indispensable core links in reactive power planning, many scholars have separately studied the reactive power planning, but still have many problems.

On the basis of a distribution method of dynamic reactive power compensation equipment, a dynamic reactive power optimization configuration method [ J ] based on an improved track sensitivity index, a power grid technology, 2012, 36 (2): 88-94.DOI:10.13335/j.1000-3673.pst.2012.02.020. ISSN: 1000 + 3673, CN: 11-2410/TM) defines an improved track sensitivity index, but the index does not truly reflect a track of reactive power increasing of a dynamic reactive power compensation device and accurately shows a support effect on bus voltage. The method comprises the following steps of taking account of reactive increment of a reactive compensation device on the basis of a dynamic reactive power optimization configuration method [ J ] power automation equipment 2016, 36 (9): 127-133.DOI:10.16081/j.issn.1006-6047.2016.09.019. ISSN: 1006-6047.CN: 32-1318/TM) considering grid characteristics, increasing the actual physical significance of the index, and enabling a calculation result to be more accurate, wherein the index has a strict requirement on sampling step length, and if the index is not properly arranged, the index is easy to be solved. In addition, a large amount of simulation data show that compared with a reactive voltage sensitivity index, the track sensitivity index based on the transient voltage stability margin has a remarkable advantage in point distribution. The document (treble. alternating current and direct current power grid transient voltage characteristic evaluation and control research [ D ]. Beijing: North China electric power university, 2020.DOI:10.27140/d.cnki. ghbbu.2020.000733.) is based on the track sensitivity index based on the transient voltage stability margin, but the distribution thought is lack of completeness. For the constant volume strategy of the dynamic reactive power compensation equipment, documents (Soxhlet, Liujiaqin, Jianweiong, and the like, large-scale new energy direct current delivery system phase modifier configuration research [ J ]. power automation equipment, 2019, 39 (9): 124-129.DOI:10.16081/J. epae.2019054. ISSN: 1006 and 6047, CN: 32-1318/TM.) explore the configuration scheme of the phase modifier in the large-scale new energy direct current delivery system, and the conclusion that the transient overvoltage problem of a delivery end system can be effectively solved by the low-voltage bus dispersedly arranged at each new energy collection station is obtained, and the documents do not specifically optimize the installation capacity and the access quantity of each phase modifier and have poor economy.

The document (Zhoushao, Tangfei, Liu Daojie, etc.. dynamic reactive compensation configuration method [ J ] considering reducing the risk of multi-feed direct current commutation failure [ high voltage technology, 2018, 44 (10): 3258-3265, DOI:10.13336/j.1003-6520, hve.20180925016, ISSN: 1003-.

Disclosure of Invention

Aiming at the defects of the existing power grid reactive power planning method, the invention provides the receiving end power grid reactive power planning method considering the transient voltage stability and containing the high-permeability wind power.

The technical scheme adopted by the invention is as follows:

the reactive power planning method for the receiving-end power grid containing high-permeability wind power and considering transient voltage stability comprises the following steps:

step 1: based on the criterion of a multi-binary table, a transient voltage safety margin index is provided for evaluating the transient voltage stability of the system;

step 2: based on the transient voltage safety margin index, a point distribution method of the dynamic reactive power compensation equipment is provided;

and step 3: establishing a differential dynamic reactive power compensation optimization model;

and 4, step 4: and screening out the optimal scheme under each scene by an improved entropy weight optimal solution distance method, and then determining the final configuration scheme of the dynamic reactive power compensation equipment in the system.

In the step 1, before the transient voltage safety margin index is provided, a wind power typical operation scene is constructed based on a wind power scene probability theory, as shown in formulas (1) to (6):

in formula (1), f (v) is a probability density function of wind speed; v is the wind speed; k is a shape parameter of wind speed distribution; c is a scale parameter, and the value of the scale parameter can be calculated by the wind speed mean value mu and the standard deviation sigma; e is a natural constant.

In the formula (2), P_wThe output power of the wind turbine generator is obtained; p_rIs windRated capacity of the motor set; v. of_ci、v_r、v_coRespectively the cut-in wind speed, the rated wind speed and the cut-out wind speed of the wind turbine generator;

in formulae (3) to (5), P₁、P₂、P₃The probabilities of the wind turbine generator operating in a zero output scene, an underrated output scene and a rated output scene are respectively shown.

The output power of the wind turbine generator is the output power of the wind turbine generator under the condition of underrated output.

The formula (6) gives a calculation method of the wind power output power by taking the underrated output scene as an example, and the wind power output power under other two typical scenes can be calculated sequentially by the same method.

Wind power output power in a zero output scene:

the output power of the wind turbine generator is the output power of the wind turbine generator in a zero-output scene.

Wind power output power under rated output scene:

the output power of the wind turbine generator is the output power of the wind turbine generator under a rated output scene.

In the step 1, the multi-binary-table criterion is an evaluation criterion for measuring whether the voltage of each bus of the system is unstable or not when the power system has large disturbance. By setting a voltage binary meter (V)_cr.1，T_cr.1)、(V_cr.2，T_cr.2)、…、(V_cr.n，T_cr.n) To request a certain bus voltage V_iBelow each predetermined threshold value V_cr.1、V_cr.2、…、V_cr.nMaximum duration T of_bRespectively not exceeding a prescribed time T_cr.1、T_cr.2、…、T_cr.nWhen the transient voltage of a certain bus meets the condition, the transient voltage of the bus is considered to be stable, otherwise, the transient voltage is unstable.

In step 1, the transient voltage safety margin index includes:

the transient voltage safety margin index of the bus is as follows:

in formula (7), σ is a step factor represented by formula (12); eta_iThe transient voltage safety margin indexes of the bus i under the constraint of the n low-voltage multi-binary tables and the m overvoltage binary tables; eta_i.d.nThe voltage sag safety margin index of the bus i under the constraint of the n low-voltage multi-binary tables is specifically shown in a formula (8); eta_i.r.mIs m pieces ofThe voltage sag index of the bus i under the constraint of the voltage multi-binary table is specifically shown in a formula (10):

in the formula (8), V_iNIs a rated voltage reference value; t is t_kAnd t'_kRespectively, the voltage of the bus i is lower than the voltage threshold value V in the low-voltage binary table in the dropping process_cr.d.kAnd above voltage threshold V during recovery_cr.d.kThe time of day; t is t_k+1And t'_k+1Respectively, the voltage of the bus i is lower than the voltage threshold value V in the low-voltage binary table in the dropping process_cr.d.k+1And above voltage threshold V during recovery_cr.d.k+1The time of day; v_i(t) is the transient voltage response curve of the bus i; k_nIs the weight coefficient of the nth threshold interval.

K_kThe weight coefficient of each threshold interval can be solved step by formula (9).

In the formula (9), T_cr.d.1Time threshold for the 1 st low voltage binary table; t is a unit of_cr.d.2Time threshold for 2 nd low voltage binary table; t is_cr.d.nA time threshold for the nth low voltage binary table; v_cr.d.1Is the voltage threshold of the 1 st low voltage binary table; v_cr.d.2Is the voltage threshold of the 2 nd low voltage binary table; v_cr.d.nIs the voltage threshold of the nth low voltage binary table; k₁The weight coefficient is the 1 st threshold interval; k₂The weight coefficient is the 2 nd threshold interval; k_nIs the weight coefficient of the nth threshold interval.

η_i.r.mThe voltage temporary rise index of the bus i under the constraint of m overvoltage multi-binary tables is specifically shown in a formula (10):

in the formula (10), V_cr.r.k+1The voltage threshold value of the (k + 1) th overvoltage binary table; t is t_i.r.kThe moment when the voltage of the bus i is higher than the voltage threshold value of the kth overvoltage binary table for the first time in the rising process after disturbance; t is t_i.r.k+1The moment when the voltage of the bus i is higher than the voltage threshold of the (k + 1) th overvoltage binary meter for the first time in the rising process after disturbance; t'_i.r.kThe time when the voltage of the bus i is equal to the voltage threshold value of the kth overvoltage binary meter for the first time in the descending process after the ascending is finished; m is a positive integer; t'_i.r.k+1The time when the voltage of the bus i is equal to the voltage threshold of the (k + 1) th overvoltage binary meter for the first time in the descending process after the ascending is finished; t is t_i.r.mThe moment when the voltage of the bus i is higher than the voltage threshold value of the mth overvoltage binary meter for the first time in the rising process after disturbance; t'_i.r.mThe moment when the voltage of the bus i is equal to the voltage threshold value of the mth overvoltage binary meter for the first time in the descending process after the ascending is finished; k_r.mIs {0 is less than or equal to V_i(t)≤V_cr.r.m+1}∩{t_i.r.m≤t≤t'_i.r.m+1The weight coefficient of the area between;

K_r.kis {0 is less than or equal to V_i(t)≤V_cr.r.k+1}∩{t_i.r.k≤t≤t'_i.r.k+1The magnitude of the weight coefficient of the region between (1) and (11) can be solved by the formula.

In the formula (11), k is a positive integer; t is_cr.r.1The time threshold value is the time threshold value in the 1 st overvoltage binary table; t is a unit of_cr.r.2Is the time threshold value in the 2 nd overvoltage binary table; t is_cr.r.kIs the time threshold value in the kth overvoltage binary table; t is_cr.r.k+1The time threshold value is the time threshold value in the (k + 1) th overvoltage binary table; k is_r.1Is {0 is less than or equal to V_i(t)≤V_cr.r.2}∩{t_i.r.1≤t≤t'_i.r.2The weight coefficient of the area between; k_r.kIs {0 is less than or equal to V_i(t)≤V_cr.r.k+1}∩{t_i.r.k≤t≤t'_i.r.k+1A weight coefficient of the area between; v_cr.r.1The voltage threshold value of the 1 st overvoltage binary table; v_cr.r.2Is the voltage threshold of the 2 nd overvoltage binary table; v_cr.r.kIs the voltage threshold of the kth overvoltage binary table; v_cr.r.k+1Is the voltage threshold of the (k + 1) th over-voltage binary table.

In the formula (12), e is a natural constant; t is time; omega is a kurtosis parameter and is used for characterizing the steepness degree of the function, and the omega is 1000 according to the invention.

Area voltage qualification rate index:

in formula (13), P_aThe voltage qualification rate indexes of the region a behind all N load buses of the whole system are considered, wherein the S operation modes of the system, T typical wind power scenes, M preselected fault sets and N load buses of the whole system are considered; p_lThe probability that the system operates in the operation mode l; n is a radical of_a.l.v.b.dIn order to increase the number of voltage-qualified buses in the area a when a d-type fault occurs at the load bus b in the operating mode and the wind power scene v, when eta of the bus i_i<When 1, the bus i is considered as a voltage qualified bus; n is a radical of_aThe number of all load buses in the area a; p_WT.vThe probability that wind power operates in a scene v; delta. for the preparation of a coating_l.v.dIn the running mode and the wind power scene v, the weight coefficient of the fault d is numerically equal to the probability of the fault d, and typical faults are independent of each other, so that the fault d has

For the probability of a fault d occurring on bus b, it is assumed here that a fault d occursThe probability of being equal at each bus of the system, i.e.

Regional voltage stability margin index:

in the formula (14), eta_aThe method comprises the following steps of considering S operation modes of a system, T typical wind power scenes, M preselected fault sets and voltage stability margin indexes of all N load bus rear areas a of the whole system; eta_i.l.v.d.bThe transient voltage safety margin index of a bus i in an area a is considered when a bus b has a fault d after the system runs in a mode I and a wind power scene v; p_lThe probability that the system operates in the operation mode l; p_WT.vThe probability that wind power operates in a scene v; delta_l.v.dThe weight coefficient of the fault d is under the condition of the operation mode and the wind power scene v;

the probability that the bus b has a fault d in the running mode and the wind power scene v is shown.

In the step 1, the reactive power planning is a means for improving the stability of the power grid and ensuring the safe and economic operation of the power grid by determining the optimal installation position and capacity of the reactive power compensation equipment on the basis of the power grid planning.

In the step 2, the point distribution method of the dynamic reactive power compensation equipment includes the following steps:

step 2.1: the key bus in the system is screened out according to formulas (15) to (18):

B_i.zs＝λ_iη'_i.risk (15)；

in the formula (15), B_i.zsIs the pivot value of the bus i; lambda [ alpha ]_iA weight coefficient, which is a bus i, whose value is defined by equation (16);

in the formula (16), D_iFor the degree of the bus i, reflecting the number of edges associated with the bus i;

is a weight coefficient, satisfies

In the present invention,

S_iactual injected power for bus i; s_baseThe reference value of the system power is 100 MVA;

the larger the value of the power transmitted or distributed by the bus i in the system is represented, the more power transmitted or distributed by the bus in the system is represented, and the more importance degree in the system is represented.

η'_i.riskA voltage instability risk factor for the bus i, whose value is defined by equation (17):

in the formula (17), N_l.v.d.iThe total number of the voltage instability buses in the dynamic subarea is determined when the bus i has a fault d in the operation mode and the wind power scene v; eta 'of'_l.v.g.dIn the l operation mode and the wind power scene v, when a bus i has a fault d, the transient voltage safety margin of a voltage instability bus g in a dynamic subarea is defined as eta'_l.v.g.d>1, the transient voltage of the bus g is unstable; p_lIs the probability that the system operates in the operating mode l; p_WT.vThe probability that wind power operates in a scene v; delta_l.v.dThe weight coefficient of the fault d under the operation mode l and the wind power scene v;

the probability that the bus i has a fault d under the operation mode l and the wind power scene v is shown; eta 'of'_i.riskAfter the bus i fails, the average size of the safety margin indexes of all the unstable buses in the dynamic partition is reflected, and the larger the value of the safety margin indexes is, the larger the influence degree of the bus on the fault is.

Determining a key busbar set of the system by calculating the pivot value of each busbar and arranging the pivot value in descending order according to a formula (18):

P_bus＝{i|sort{B_i.zs},i∈{1,2,...,N}} (18)；

in formula (18), sort { B_i.zsEach bus is according to B_i.zsThe values are sorted in descending order to form a set; b is_i.zsIs the pivot value of the bus i; i is the serial number of the bus; 1,2, N is the N load bus numbers within the system.

Step 2.2: and constructing a candidate bus bar set to be compensated by the formulas (19) to (23).

In the formula (19), SI1.i is a sensitivity index of the bus i based on the transient voltage safety margin of the bus; n is a radical of_l.v.d.iThe total number of the voltage instability buses in the dynamic subarea is determined when the bus i has a fault d in the operation mode and the wind power scene v; eta_g0.dBefore a dynamic reactive power compensation device is installed at a bus i, transient voltage safety margin indexes of a bus g when a system fails d after n low-voltage multi-binary tables and m overvoltage binary tables are considered; eta_g.dAfter a certain dynamic reactive power compensation device is installed at a bus i, when a system has a fault d, the transient voltage safety margin indexes of a bus g under the constraint of n low-voltage multi-binary tables and m overvoltage binary tables are obtained; delta Q_c.iIs the capacity of the dynamic reactive compensation equipment installed at the bus i.

If during the transient process a higher demand is placed on the fast voltage support capability of the dynamic reactive power compensation equipment, the system dynamic reactive power backup of the dynamic reactive power compensation equipment i defined by equation (20) can be used to further modify the si1. i.

Q_RTSi＝∫e^-t(Q_i-Q_i0)dt (20)；

In the formula (20), Q_RTSiThe method is characterized in that the method is a quantitative value for representing the dynamic reactive power reserve of the dynamic reactive power compensation equipment connected to a bus i; q_iThe reactive power is increased in the transient process of the dynamic reactive power compensation equipment connected to the bus i; q_i0The initial reactive power of the dynamic reactive power compensation equipment connected to the bus i is in steady-state operation; e.g. of a cylinder^-tThe introduced attenuation factor is used for quantifying the reactive power of the reactive compensation equipment at each moment in the transient process, and the faster the reactive power of the dynamic reactive compensation equipment is increased is more beneficial to the transient voltage stability of the system.

Based on the formula (19) and the formula (20), the sensitivity indexes of the transient voltage safety margin and the dynamic reactive response rate based on the bus shown in the formula (21) are defined.

In the formula (21), si2.i is a sensitivity index based on the transient voltage safety margin and the dynamic reactive response rate of the bus; max { Q }_RTSiI ∈ {1,2, …, N } } is for each Q_RTSiThe largest of the values; and SI1.i is a sensitivity index of the bus i based on the transient voltage safety margin of the bus.

A candidate busbar set based on the si1.i as shown in the formula (22) and a candidate busbar set based on the si2.i as shown in the formula (23) are further constructed.

C_SI1.bus＝{i|sort{SI1.i},i∈{1,2,...,N}} (22)；

In the formula (22), C_SI1.busIs a candidate bus set based on SI1. i; sort { SI1.i } is a set in which the buses are arranged in descending order according to the size of the SI1.i value.

C_SI2.bus＝{i|sort{SI2.i},i∈{1,2,...,N}} (23)；

In the formula (23), C_SI2.busIs a candidate bus set based on SI2. i; sort { SI2.i } is a set of buses arranged in descending order according to the size of the SI2.i value.

In the step 3, the objective function of the differentiated dynamic reactive power compensation optimization model is shown in formulas (24) to (26):

f₁(x)＝ω₁[(1-P_a)+η_a] (24)；

in the formula (24), P_aThe voltage qualification rate index of the area a; eta_aIs an index of the voltage stability margin of the region a.

min f＝{f₁(x),f₂(x)} (26)；

In formulae (24) to (26), f₁(x) And f₂(x) Respectively representing the dynamic reactive power compensation effect and the dynamic reactive power compensation economic cost for two sub-target functions to be optimized; omega₁、ω₂The optimized weight of the sub-objective function is satisfied with omega₁+ω₂＝1；1-P_aThe voltage instability rate index of the system is obtained; t is a unit of₂The operating life of the SVC; c_svc.uThe reactive compensation unit price of the SVC is achieved; q_svc.uIs the reactive compensation capacity of SVC; f_svc.uThe installation cost of SVC; t is₁The running age of the STATCOM; c_STATCOM.hThe reactive compensation unit price is the STATCOM reactive compensation unit price; q_STATCOM.hThe reactive compensation capacity of the STATCOM is obtained; f_STATCOM.hThe installation cost of the STATCOM; z_l.vThe number of compensation nodes for installing SVC under the operation mode l and the wind power scene v; h_l.vThe number of compensation nodes for installing the STATCOM in the running mode and the wind power scene v is increased; c_eIs the electricity price; ζ is the annual maximum load hours; p_l.v.lossThe network loss of the system is determined under the condition of the operation mode and the wind power scene v; minf denotes that the two sub-targeting functions are minimal.

In the step 3, the differentiated dynamic reactive power compensation is realized by: 1) dynamics of installation for different candidate nodesThe types of reactive compensation equipment are different; 2) a double optimization strategy is employed to enforce reactive planning. The double optimization strategy is an optimization idea which respectively considers two situations of optimal reactive compensation effect and optimal economic cost. When the optimal reactive compensation effect of the system during the transient process is considered, the optimization weight omega of the sub-target function is made₁:ω₂2:1, i.e. ω₁＝0.67，ω₂0.33; and when the economic cost is optimized, let omega₁:ω₂1:2, i.e.. omega₁＝0.33，ω₂＝0.67。

In the step 4, the improved entropy weight TOPSIS method is a decision evaluation method based on the entropy weight method, the coefficient of variation method and the TOPSIS method, and for s evaluation schemes, w evaluation indexes have the following evaluation steps:

step (1) of applying an extreme value method to enable an evaluation matrix A to be [ a ]_xy]_s×wNormalized to obtain a normalized evaluation matrix A '═ a'_xy]_s×wWherein: a'_xyThe values are as follows:

in formula (27), x is the current evaluation scheme; y is the current evaluation index; a is_xyIs the initial value of the index; a'_xyIs the new value of the index;

respectively corresponding to the evaluation index y and a_xyMaximum and minimum values of.

Step (2), calculating the coefficient of variation V according to the formula (28)_y；

In the formula (28), the reaction mixture is,

S_yare respectively a 'shown as a formula (29) and a formula (30)'_xyMean and standard deviation of (d).

Step (3) calculating the coefficient of variation V_yWeighted value W of each index_1y；

Step (4) of'_xyCalculating the information entropy E of the discrete distribution situation_y；

In the formula (32), E_yIs information entropy; s is the total number of evaluation protocols; x is the current evaluation scheme; a'_xyIs the new value of the index; ln is a logarithmic function based on natural constants.

Step (5), calculating the entropy E based on the information_yWeighted value W of each index_2y；

Step (6) of calculating the final weight W of each index_yConstructing a weighting matrix R of the method;

R＝(W_y×a'_xy)_s×w＝(r_xy)_s×w(35)；

in formula (35), R is a weighting matrix; w_yIs the final weight of the index; a'_xyIs the new value of the index; r is_xyThe index value in the weighting matrix; s is the total number of evaluation scenarios; w is the total number of evaluation indexes.

Step (7), determining the optimal scheme under each index by the weighting matrix R

And worst case scenario

In the formula (36), the reaction mixture is,

is an optimal scheme; r is_xyThe index value in the weighting matrix; x is the current evaluation scheme; y is the current evaluation index; w is the total number of evaluation indexes;

is an index r_xyMaximum value of (2).

In the formula (37), the reaction mixture is,

the best scheme is adopted; r is_xyThe index value in the weighting matrix; x is the current evaluation scheme; y is the current evaluation index; w is the total number of evaluation indexes;

is an index r_xyMinimum value of (1).

Step (ii) of(8) Sequentially calculating each scheme to be evaluated and the optimal scheme

And the worst scheme

European distance between

In the formula (38), the reaction mixture is,

for each solution to be evaluated and the optimal solution

The Euclidean distance between; y is the current evaluation index; w is the total number of evaluation indexes;

is an optimal scheme; r is_xyIs an index value in the weighting matrix.

In the formula (39), the compound represented by the formula (I),

for each solution to be evaluated and the worst solution

the best scheme is adopted; r is_xyIs an index value in the weighting matrix.

Step (9), calculating the closeness F of each scheme to be evaluated and the ideal scheme_x；

In step 4, the final configuration scheme is determined by formula (41):

in formula (41), PRO_finConfiguring a scheme for the final dynamic reactive power compensation equipment;

economic cost of the optimal scheme screened out under the operation mode and the wind power scene v; l is the operation mode of the system; s is the total number of the system operation modes; v is the running scene of wind power; and T is the total number of the typical wind power scenes.

The invention relates to a receiving-end power grid reactive power planning method considering transient voltage stability and containing high-permeability wind power, which has the following technical effects:

1) aiming at a receiving-end power grid containing high-permeability new energy, the influence of DG uncertainty on the receiving-end power grid needs to be considered during reactive power planning, the invention provides a differentiated configuration method of dynamic reactive power compensation equipment, constructs a set of transient voltage safety margin indexes capable of guiding reactive power planning of the power grid, and further provides a corresponding distribution method to reasonably determine a configuration scheme of the dynamic reactive power compensation equipment in a system.

2) The invention reflects the objective factor of wind power uncertainty in reactive power planning, and the planning result is more real and reasonable.

3) The transient voltage safety margin index based on the multiple binary table criterion can accurately evaluate the transient voltage stability of the system and provide guidance for reactive power planning.

4) The reactive power planning method not only can ensure the reactive power compensation effect, but also can save the economic cost to the maximum extent.

Drawings

Fig. 1 is an IEEE39 node system diagram.

FIG. 2 is a graph showing the results of comparison of sensitivity indexes.

Fig. 3 is a graph of compensated fault bus voltage waveforms.

Fig. 4 is a graph comparing the results of si2.i and si1. i.

Fig. 5 is a voltage waveform diagram of a key bus of a system part.

FIG. 6(a) is a voltage waveform diagram (uncompensated) of a full-network bus under each working condition;

fig. 6(b) is a voltage waveform diagram of the full-grid bus under each working condition (wind power zero output scene);

fig. 6(c) is a voltage waveform diagram of the full-grid bus under each working condition (wind power underrun output scene);

fig. 6(d) is a voltage waveform diagram (wind power rated output scene) of the full-grid bus under various working conditions;

FIG. 7 is a flow chart of the method of the present invention.

Detailed Description

As shown in fig. 7, the method for planning reactive power of the receiving-end power grid including high-permeability wind power considering transient voltage stabilization includes the following steps:

step 1: and on the basis of the multi-binary-table criterion, providing a transient voltage safety margin index for evaluating the transient voltage stability of the system.

Before a transient voltage safety margin index is provided, a wind power operation scene is reduced into three typical scenes of a zero output scene, an underrated output scene and a rated output scene based on a wind power scene probability theory, the operation probability and the output power of a wind turbine generator under each scene are calculated, and the specific process is shown in formulas (1) to (6):

in the formula (1), f (v) is a probability density function of wind speed; v is the wind speed; k is a shape parameter of wind speed distribution; c is a scale parameter, and the value of the scale parameter can be calculated by the wind speed mean value mu and the standard deviation sigma; e is a natural constant.

In the formula (2), P_wThe output power of the wind turbine generator is obtained; p_rThe rated capacity of the wind turbine generator is set; v. of_ci、v_r、v_coRespectively the cut-in wind speed, the rated wind speed and the cut-out wind speed of the wind turbine generator;

Wind power output power in a zero output scene:

Wind power output power under rated output scene:

the output power of the wind turbine generator is the output power of the wind turbine generator in a rated output scene.

And (3) constructing a bus transient voltage safety margin index, a regional voltage qualification rate index and a regional voltage stability margin index, wherein the specific processes are shown in formulas (7) to (14):

in the formula (7), wherein eta_iThe transient voltage safety margin indexes of the bus i under the constraint of the n low-voltage multi-binary tables and the m overvoltage binary tables; eta_i.d.nThe voltage sag safety margin index of the bus i under the constraint of the n low-voltage multi-binary tables is specifically shown in a formula (8); eta_i.r.mThe voltage temporary rise index of the bus i under the constraint of the m overvoltage multi-binary tables is specifically shown as a formula (10); σ is a step factor shown in equation (12).

In the formula (8), V_iNIs a foreheadA constant voltage reference value; t is t_kAnd t'_kRespectively, the voltage of the bus i is lower than the voltage threshold value V in the low-voltage binary table in the dropping process_cr.d.kAnd above voltage threshold V during recovery_cr.d.kThe time of day; k_kThe weight coefficient in each threshold interval can be solved step by the formula (9); t is t_k+1And t'_k+1Respectively, the voltage of the bus i is lower than the voltage threshold value V in the low-voltage binary table in the dropping process_cr.d.k+1And above voltage threshold V during recovery_cr.d.k+1The time of day; v_i(t) is the transient voltage response curve of the bus i; k is_nIs the weight coefficient of the nth threshold interval.

In the formula (9), T_cr.d.1Time threshold for the 1 st low voltage binary table; t is a unit of_cr.d.2Time threshold for the 2 nd low voltage binary table; t is a unit of_cr.d.nA time threshold for the nth low voltage binary table; v_cr.d.1Is the voltage threshold of the 1 st low voltage binary table; v_cr.d.2Is the voltage threshold of the 2 nd low voltage binary table; v_cr.d.nIs the voltage threshold of the nth low voltage binary table; k₁The weight coefficient is the 1 st threshold interval; k₂The weight coefficient is the 2 nd threshold interval; k_nIs the weight coefficient of the nth threshold interval.

In the formula (10), V_cr.r.k+1The voltage threshold value of the (k + 1) th overvoltage binary table; t is t_i.r.kThe moment when the voltage of the bus i is higher than the voltage threshold of the kth overvoltage binary meter for the first time in the rising process after disturbance; t is t_i.r.k+1The moment when the voltage of the bus i is higher than the voltage threshold of the (k + 1) th overvoltage binary meter for the first time in the rising process after disturbance; t'_i.r.kThe bus i voltage is equal to the k over-voltage for the first time in the descending process after the ascending is finishedThe moment of pressing the binary meter voltage threshold; m is a positive integer; t'_i.r.k+1The time when the voltage of the bus i is equal to the voltage threshold of the (k + 1) th overvoltage binary meter for the first time in the descending process after the ascending is finished; t is t_i.r.mThe moment when the voltage of the bus i is higher than the voltage threshold value of the mth overvoltage binary meter for the first time in the rising process after disturbance; t'_i.r.mThe moment when the voltage of the bus i is equal to the voltage threshold value of the mth overvoltage binary meter for the first time in the descending process after the ascending is finished; k_r.mIs {0 is less than or equal to V_i(t)≤V_cr.r.m+1}∩{t_i.r.m≤t≤t'_i.r.m+1The weight coefficient of the area between; k_r.kIs {0 is less than or equal to V_i(t)≤V_cr.r.k+1}∩{t_i.r.k≤t≤t'_i.r.k+1The magnitude of the weight coefficient of the region between (1) and (11) can be solved by the formula.

In the formula (11), k is a positive integer; t is_cr.r.1The time threshold value is the time threshold value in the 1 st overvoltage binary table; t is a unit of_cr.r.2Is the time threshold value in the 2 nd overvoltage binary table; t is a unit of_cr.r.kIs the time threshold value in the kth overvoltage binary table; t is_cr.r.k+1The time threshold value is the time threshold value in the (k + 1) th overvoltage binary table; k_r.1Is {0 is less than or equal to V_i(t)≤V_cr.r.2}∩{t_i.r.1≤t≤t'_i.r.2The weight coefficient of the area between; k_r.kIs {0 is less than or equal to V_i(t)≤V_cr.r.k+1}∩{t_i.r.k≤t≤t'_i.r.k+1A weight coefficient of the area between; v_cr.r.1The voltage threshold value of the 1 st overvoltage binary table; v_cr.r.2Is the voltage threshold of the 2 nd overvoltage binary table; v_cr.r.kA voltage threshold value of a kth overvoltage binary table; v_cr.r.k+1Is the voltage threshold of the (k + 1) th overvoltage binary table.

In formula (13), P_aThe voltage qualification rate indexes of the region a behind all N load buses of the whole system are considered, wherein the S operation modes of the system, T typical wind power scenes, M preselected fault sets and N load buses of the whole system are considered; p is_lThe probability that the system operates in the operation mode l; n is a radical of_a.l.v.b.dIn order to increase the number of voltage-qualified buses in the area a when a d-type fault occurs at the load bus b in the operating mode and the wind power scene v, when eta of the bus i_i<1, the bus i is considered to be a voltage qualified bus; n is a radical of_aThe number of all load buses in the area a; p_WT.vThe probability that wind power operates in a scene v; delta_l.v.dIn order to ensure that the weight coefficient of the fault d is numerically equal to the probability of the fault under the condition of the operation mode and the wind power scene v, the typical faults are independent of each other, so that the fault d has

In order to ensure the probability of the fault d of the bus b under the condition of the operation mode and the wind power scene v, the probability of the fault d occurring at each bus of the system is assumed to be equal, namely

In the formula (14), eta_aThe method comprises the following steps of considering S operation modes of the system, T typical wind power scenes, M preselected fault sets and voltage stability margin indexes of all N load bus rear areas a of the whole system; eta_i.l.v.d.bThe transient voltage safety margin index of a bus i in an area a is considered when a bus b has a fault d after the system runs in a mode I and a wind power scene v; p_lThe probability that the system operates in the operation mode l; p_WT.vThe probability that wind power operates in a scene v; delta. for the preparation of a coating_l.v.dThe weight coefficient of the fault d under the operation mode l and the wind power scene v;

the probability that the bus b has a fault d under the operation mode l and the wind power scene v is shown.

Step 2: and providing a distribution method of the dynamic reactive power compensation equipment based on the transient voltage safety margin index.

Before the distribution, the key buses in the system need to be screened out according to the formulas (15) to (18).

B_i.zs＝λ_iη'_i.risk (15)；

In the formula (15), B_i.zsIs the pivot value of the bus; lambda [ alpha ]_iA weight coefficient, which is a bus i, whose value is defined by equation (16); eta 'of'_i.riskThe value of the voltage instability risk factor for the bus i is defined by equation (17).

is a weight coefficient, satisfies

In the invention

S_iActual injected power for bus i; s. the_baseThe reference value of the system power is taken as 100 MVA;

In the formula (17), N_l.v.d.iThe total number of the voltage instability buses in the dynamic subarea is determined when the bus i has a fault d in the operation mode and the wind power scene v; eta 'of'_l.v.g.dIn the l operation mode and the wind power scene v, when a bus i has a fault d, the transient voltage safety margin of a voltage instability bus g in a dynamic subarea is defined as eta'_l.v.g.d>1, the transient voltage of the bus g is unstable; p_lThe probability that the system operates in the operation mode l; p_WT.vThe probability that wind power operates in a scene v; delta_l.v.dThe weight coefficient of the fault d under the operation mode l and the wind power scene v;

P_bus＝{i|sort{B_i.zs},i∈{1,2,...,N}} (18)；

in formula (18), sort { B_i.zsIs that each bus is according to B_i.zsThe values are sorted in descending order to form a set; b is_i.zsIs the pivot value of the bus i; i is the serial number of the bus; 1,2, N is the number of N load buses in the system.

And then, constructing a candidate bus bar set to be compensated by the formulas (19) to (23).

Q_RTSi＝∫e^-t(Q_i-Q_i0)dt (20)；

In the formula (20), Q_RTSiThe method is characterized in that the method is a quantitative value for representing the dynamic reactive power reserve of the dynamic reactive power compensation equipment connected to a bus i; q_iThe reactive power is increased in the transient process of the dynamic reactive power compensation equipment connected to the bus i; q_i0The initial reactive power of the dynamic reactive power compensation equipment connected to the bus i is in steady-state operation; e.g. of the type^-tThe introduced attenuation factor is used for quantifying the reactive power of the reactive compensation equipment at each moment in the transient process, and the faster the reactive power of the dynamic reactive compensation equipment is increased is more beneficial to the transient voltage stability of the system.

Based on the formula (19) and the formula (20), the sensitivity index of the bus-based transient voltage safety margin and the dynamic reactive response rate is defined as shown in the formula (21).

In the formula (21), si2.i is a sensitivity index based on the transient voltage safety margin and the dynamic reactive response rate of the bus; max { Q }_RTSiI ∈ {1,2, …, N } } is each Q_RTSiThe largest of the values; and SI1.i is a sensitivity index of the bus i based on the transient voltage safety margin of the bus.

An SI1. i-based candidate busbar set as shown in formula (22) and an SI2. i-based candidate busbar set as shown in formula (23) are further constructed.

C_SI1.bus＝{i|sort{SI1.i},i∈{1,2,...,N}} (22)；

In the formula (22), C_SI1.busIs a candidate bus set based on SI1. i; sort { SI1.i } is a set of buses arranged in descending order according to the size of the SI1.i value.

C_SI2.bus＝{i|sort{SI2.i},i∈{1,2,...,N}} (23)；

And step 3: establishing a differential dynamic reactive power compensation optimization model, and solving a Pareto optimal solution set under each scene through a multi-objective grey wolf algorithm (MOGWO).

Establishing an objective function of a differential dynamic reactive power compensation optimization model, as shown in formulas (24) to (26):

f₁(x)＝ω₁[(1-P_a)+η_a] (24)；

min f＝{f₁(x),f₂(x)} (26)；

In formulae (24) to (26), f₁(x) And f₂(x) Respectively representing the dynamic reactive power compensation effect and the dynamic reactive power compensation economic cost for two sub-target functions to be optimized; omega₁、ω₂The optimization weight of the sub-objective function is satisfied₁+ω₂＝1；1-P_aThe voltage instability rate index of the system is obtained; t is₂The operating life of the SVC; c_svc.uThe reactive compensation unit price of the SVC is achieved; q_svc.uIs the reactive compensation capacity of SVC; f_svc.uThe installation cost of SVC; t is₁The running age of the STATCOM; c_STATCOM.hThe reactive compensation unit price is STATCOM; q_STATCOM.hThe reactive compensation capacity of the STATCOM is obtained; f_STATCOM.hThe installation cost of the STATCOM; z_l.vThe number of compensation nodes for installing SVC under the operation mode l and the wind power scene v; h_l.vThe number of compensation nodes for installing the STATCOM in the running mode and the wind power scene v is increased; c_eIs the electricity price; ζ is the maximum annual load hours; p_l.v.lossThe network loss of the system is determined under the condition of the operation mode and the wind power scene v; minf denotes that the two sub-targeting functions are minimal.

The method comprises the following main steps of solving a Pareto optimal solution set under each scene through a multi-objective wolf algorithm (MOGWO):

1) setting a population number N of wolfs corresponding to a capacity of a dynamic reactive power compensation device_pMaximum iteration number I, dimension D, and upper and lower limits U of capacity of dynamic reactive power compensation equipment_bAnd L_bNumber of outside population Archive N_ArAnd related control parameters;

2) randomly generating the wolf individuals meeting the conditions, calculating the objective function value of each individual, determining the current optimal individual, and updating Archive;

3) selecting three wolfs of alpha, beta and delta in Archive based on a roulette method, sequentially updating the positions of the rest wolfs according to the positions of the wolfs, and calculating corresponding objective function values of the wolfs;

4) finding out a new non-dominant individual according to the objective function value of each wolf, and updating Archive;

5) repeating the step 3) and the step 4) until the iteration is carried out for I times, if the number of the external population reaches the upper limit in the process, randomly eliminating part of individuals from the crowded group according to the grouping result of the objective function of each individual in the Archive until the number of the external population Archive is equal to N_Ar；

6) And outputting the gray wolf position in the Archive, wherein the position represents the solved Pareto optimal solution set.

The improved entropy weight TOPSIS method is used for evaluating s evaluation schemes and w evaluation indexes, and the steps are shown in formulas (27) to (40):

step (1) of applying an extreme method to set the evaluation matrix A to [ a ═ a-_xy]_s×wNormalized to obtain a normalized evaluation matrix A '═ a'_xy]_s×wWherein a'_xyThe values are as follows:

in formula (27), x is the current evaluation scheme; y is the current evaluation index; a is a_xyIs the initial value of the index; a'_xyIs the new value of the index;

respectively corresponding to the evaluation index y and a_xyA maximum value and a minimum value of (c).

In the formula (28), the reaction mixture is,

In the formula (32), E_yIs the information entropy; s is the total number of evaluation scenarios; x is the current evaluation scheme; a'_xyIs a new value of the index; ln is a logarithmic function based on natural constants.

Step (5) calculating information-based entropy E_yWeighted value W of each index_2y；

R＝(W_y×a'_xy)_s×w＝(r_xy)_s×w (35)；

And worst case scenario

In the formula (36), the reaction mixture is,

is an index r_xyMaximum value of (2).

In the formula (37), the reaction mixture is,

is an index r_xyMinimum value of (1).

And the worst scheme

European distance between

In the formula (38), the reaction mixture is,

for each solution to be evaluated and the optimal solution

is an optimal scheme; r is_xyIs an index value in the weighting matrix.

In the formula (39), the compound represented by the formula (I),

for each solution to be evaluated and the worst solution

The euclidean distance between them; y is the current evaluation index; w is the total number of evaluation indexes;

is the worst scheme; r is_xyIs an index value in the weighting matrix.

Step (9) calculating the closeness F of each scheme to be evaluated and the ideal scheme_x；

Determining the final configuration scheme of the dynamic reactive compensation equipment in the system according to the formula (41):

Example (b):

the parameters involved in the invention are as follows:

1) the invention sets the following low-voltage multi-binary meter and over-voltage binary meter:

[0.80p.u.,10s ], [0.75p.u.,1s ], [0.7p.u.,0.2s ], [0.65p.u.,0.1s ], [1.1p.u.,1.2s ], [1.15p.u.,0.1s ]. At this time, the number n of the low-voltage multi-binary meters and the number m of the overvoltage multi-binary meters take values of 4 and 2 respectively;

2) considering the entire system as the area to be studied, i.e. N_a＝N＝39；

3) In reactive power planning, only the normal operation mode of a power system is taken as an example, and the operation working conditions of wind power under three typical scenes are considered, namely S is 1; t is 3; p_WT.1＝0.1561；P_WT.2＝0.6837；P_WT.3＝0.1602；

4) Three-phase short-circuit faults and single-phase earth faults form a preselected fault set, and the probability that each fault occurs at each part of the system is assumed to be equal, namely M is 2; delta_l.v.1＝0.07；δ_l.v.2＝0.93；

5) The bus weighting coefficients calculated from the network topology and the load flow information shown in fig. 1 by the formula (16) are shown in table 1:

TABLE 1 weighting factor of each busbar

6) From the results of the system N-1 fault scan, the bus voltage instability risk factors calculated by equation (17) are shown in Table 2:

TABLE 2 Risk factors for Voltage instability of the buses

7) The economic parameters and the investment costs of the dynamic reactive power compensation system are given in Table 3, and C_e0.617 yuan/(KW · h); ζ 3600 h;

TABLE 3 economic parameters and investment costs of dynamic reactive power compensation equipment

8)U_i.min、U_i.maxRespectively taking 0.95p.u. and 1.05 p.u.; q_i.min＝0；Q_i.maxDepending on the type of dynamic reactive compensation equipment. According to a typical model (Liu Zhen ya, Zhang Qin Ping, Wang Yating, etc.) recommended by the great Motor Special Committee of the Chinese Motor engineering society, research on reactive compensation measures for improving the safety and stability level of a New Ganqing 750kV sending-end power grid in northwest [ J]China electrical engineeringAcademic, 2015, 35 (5): 1015-1022.DOI 10.13334/j.0258-8013.pcsee.2015.05.001. ISSN: 0258-8013, CN: 11-2107/TM) and the relevant enterprise standards of the national grid company for STATCOM (national grid company. Q/GDW 241.1-2008 section 1 of the chain static synchronous compensator: function standard guide rule [ S ]]Beijing: the reliability of STATCOM equipment with the capacity of 300Mvar and above needs to be verified by the national grid company 2008, DL/T1215.1-2013), so the invention aims at STATCOM, Q_i.maxSet to 300 Mvar; and for SVC, Q_i.maxTaking the mass as 200 Mvar;

9) in MOGWO algorithm, N_p100; i ═ 50, D ═ 3, and 4; for STATCOM, U_b300 for SVC, U_b＝200；L_b＝0；N_Ar＝20。

And (3) carrying out simulation analysis on the receiving-end power grid containing the high-permeability wind power shown in the figure 1 based on a PSD-BPA and MATLAB combined simulation platform, and verifying the rationality and the effectiveness of the provided reactive power planning method. The system N-1 fault scanning result shows that the influence on the key bus and the whole-network bus is the most serious after the head end of the line 16-17 generates the three-phase short-circuit fault, so the invention sets that the head end of the 5 th cyclic line 16-17 generates the three-phase short-circuit fault, the breaker at the head end of the following 4.5 cyclic lines acts, and the breaker at the tail end of the 5 th cyclic line acts, and the optimal access position of the dynamic reactive power compensation is explored on the basis of the three-phase short-circuit fault.

The distribution point of the invention selects 100Mvar SVC, and examines the voltage supporting condition of the SVC to the key bus and the whole network bus by connecting the SVC to each bus in sequence, and because the buses 30-39 are generator buses and the voltage amplitude is constant, the buses can not be considered when the distribution point is carried out. To show the superiority of the point distribution method, the sensitivity index of equation (15) based on the bus transient voltage safety margin is compared with the conventional method 1 and the conventional method 2, and the comparison result is shown in fig. 2.

As can be seen from fig. 2, the bus sensitivity ranking results obtained by the three methods are different, and the fault bus voltage waveform and the system transient voltage safety index when the candidate nodes screened by the method of the present invention are the bus 22, the bus 21, and the bus 23 … …, and the best candidate nodes screened by the

prior methods

1 and 2 are the bus 15, the bus 14, the bus 23 … …, the bus 8, the bus 26, and the bus 28 … … 100Mvar, respectively, and the SVC is sequentially connected to the bus 22, the bus 15, and the bus 8 are shown in fig. 3 and table 4, respectively.

TABLE 4 System transient Voltage safety index

As can be seen from fig. 3 and table 4, when the dynamic reactive power compensation equipment is distributed based on the method of the present invention, during the transient process, the faulty bus and all buses in the system can obtain better compensation effect, which highlights the advantages of the method of the present invention. In addition, candidate bus bars capable of providing fast voltage support can be further selected from the candidate bus bar set constructed by the formula (17), and fig. 4 compares the sensitivity indexes of the bus bars based on the formula (15) and the formula (17).

From the size of SI2.i of each bus bar in FIG. 4, equation (19) constructs candidate bus bar set C consisting of bus bar 19, bus bar 23, and bus bar 22 … …_SI2.bus. FIG. 5 compares SVC connected to C respectively_SI1.busAnd C_SI2.busAfter the top-ranked bus, the voltage waveform of the critical bus of the system part during the transient process. As can be seen from FIG. 5, compare C_SI1.busSelecting C_SI2.busIn the bus distribution point, the key bus of the system can show better transient voltage characteristics in the early stage of the transient process, and the reasonability and the correctness of the sensitivity index based on the bus transient voltage safety margin and the dynamic reactive response speed are verified.

Aiming at the bus with the front centralized sequencing candidate buses, the STATCOM is installed on the bus, the voltage of the whole network bus is supported by virtue of the quick response capability of the STATCOM during the transient process, and SVCs are sequentially installed on the rest buses according to the sequencing sequence, so that the economic cost of reactive power compensation is emphasized. Tables 5 and 6 are based on C, respectively_SI1.busAnd C_SI2.busAnd the system voltage stability margin index under the preset compensation scheme is shown.

TABLE 5 base on C_SI1.busPreset compensation scheme of

TABLE 6 bases on C_SI2.busPreset compensation scheme of

Tables 5 and 6 show that when the reactive compensation capacity is constant, the STATCOM is installed on the bus which is arranged in the candidate bus in a centralized and front-ranked mode, so that the transient voltage stability of the system is better, the scientificity of the reactive planning strategy is proved, and a foundation is laid for subsequent constant volume.

According to the screened C_SI1.busAnd C_SI2.busBased on the double optimization strategy, the optimal capacity of the dynamic reactive power compensation equipment installed on each bus under different scenes is solved by using an MOGWO algorithm, a Pareto optimal solution set of the wind power under a rated output scene is given in a table 7, and an evaluation result of each optimal scheme after being evaluated by an improved entropy weight TOPSIS method is shown in a table 8.

TABLE 7 optimal configuration scheme of dynamic reactive power compensation equipment in rated output scene

TABLE 8 evaluation results of each optimal solution under rated output scenario

Table 8 shows that the better the reactive compensation, the higher the corresponding economic cost. And comprehensively considering, selecting a scheme 3 with the maximum closeness as an optimal configuration scheme of the dynamic reactive power compensation equipment in the wind power rated output scene.

The optimal configuration scheme of the dynamic reactive power compensation equipment in the wind power zero output scene and the underrated output scene is shown in table 9.

Configuration scheme of table 9 dynamic reactive power compensation equipment in different wind power scenes

To determine the final configuration scheme of the dynamic reactive power compensation equipment in the system, table 10 summarizes the optimal scheme evaluation results of each wind power scene.

TABLE 10 evaluation results of optimal scenarios under each wind power scenario

And (4) determining the optimal scheme of the wind power under the zero output scene as the final configuration scheme of the dynamic reactive power compensation equipment of the system by using the data of the table 10 according to the formula (37). Fig. 6(a) -6 (d) show the transient voltage waveforms of the full-grid bus in various operating conditions under the final configuration scheme. Fig. 6(a) -6 (d) illustrate that all buses of the system can maintain the transient voltage stability in all working conditions under the effect of the final configuration scheme, which highlights the applicability of the constant volume method of the invention.

Claims

1. The reactive power planning method of the receiving-end power grid containing high-permeability wind power and considering transient voltage stability is characterized by comprising the following steps of:

step 2: based on the transient voltage safety margin index, a distribution method of the dynamic reactive power compensation equipment is provided;

2. The reactive power planning method for the receiving-end power grid containing the high-permeability wind power considering the transient voltage stability according to claim 1, wherein the method comprises the following steps: in the step 1, before the transient voltage safety margin index is provided, a wind power typical operation scene is constructed based on a wind power scene probability theory, as shown in formulas (1) to (6):

in the formula (1), f (v) is a probability density function of wind speed; v is the wind speed; k is a shape parameter of wind speed distribution; c is a scale parameter, and the value of the scale parameter can be calculated by the wind speed mean value mu and the standard deviation sigma; e is a natural constant;

in formulae (3) to (5), P₁、P₂、P₃The probabilities of the wind turbine generator operating in a zero output scene, an underrated output scene and a rated output scene are respectively set;

the output power of the wind turbine generator is the output power of the wind turbine generator under the condition of underrated output;

wind power output power in a zero output scene:

the output power of the wind turbine generator is the output power of the wind turbine generator in a zero-output scene;

wind power output power under rated output scene:

3. The method of claim 1 for wind power with high permeability considering transient voltage stabilizationThe reactive power planning method for the receiving-end power grid is characterized by comprising the following steps: in step 1, the multiple binary table criterion refers to: by setting a voltage binary meter (V)_cr.1，T_cr.1)、(V_cr.2，T_cr.2)、…、(V_cr.n，T_cr.n) To request a certain bus voltage V_iBelow each predetermined threshold value V_cr.1、V_cr.2、…、V_cr.nMaximum duration T of_bRespectively not exceeding a prescribed time T_cr.1、T_cr.2、…、T_cr.nWhen the transient voltage of a certain bus meets the condition, the transient voltage of the bus is considered to be stable, otherwise, the transient voltage is unstable.

4. The reactive power planning method for the receiving-end power grid containing the high-permeability wind power considering the transient voltage stability according to claim 1, wherein the method comprises the following steps: in step 1, the transient voltage safety margin indicator includes:

the transient voltage safety margin index of the bus is as follows:

in formula (7), σ is a step factor represented by formula (12); eta_iThe transient voltage safety margin indexes of the bus i under the constraint of the n low-voltage multi-binary tables and the m overvoltage binary tables; eta_i.d.nThe voltage sag safety margin index of the bus i under the constraint of the n low-voltage multi-binary tables is specifically shown in a formula (8); eta_i.r.mThe voltage transient index of the bus i under the constraint of the m overvoltage multi-binary tables is specifically shown in a formula (10):

in the formula (8), V_iNIs a nominal voltage reference value; t is t_kAnd t'_kRespectively, the voltage of the bus i is lower than the voltage threshold value V in the low-voltage binary table in the dropping process_cr.d.kAnd recoveredAbove voltage threshold V in stroke_cr.d.kThe time of day; t is t_k+1And t'_k+1Respectively, the voltage of the bus i is lower than the voltage threshold value V in the low-voltage binary table in the dropping process_cr.d.k+1And above voltage threshold V during recovery_cr.d.k+1The time of day; v_i(t) is the transient voltage response curve of the bus i; k_nA weight coefficient for the nth threshold interval;

K_kthe weight coefficient of each threshold interval can be solved step by the formula (9);

in the formula (9), T_cr.d.1Time threshold for the 1 st low voltage binary table; t is_cr.d.2Time threshold for the 2 nd low voltage binary table; t is_cr.d.nA time threshold for the nth low voltage binary table; v_cr.d.1Is the voltage threshold of the 1 st low voltage binary table; v_cr.d.2Is the voltage threshold of the 2 nd low voltage binary table; v_cr.d.nIs the voltage threshold of the nth low voltage binary table; k₁A weight coefficient of a 1 st threshold interval; k is₂A weight coefficient for the 2 nd threshold interval; k_nA weight coefficient for the nth threshold interval;

η_i.r.mthe voltage transient index of the bus i under the constraint of the m overvoltage multi-binary tables is specifically shown in a formula (10):

in the formula (10), V_cr.r.k+1The voltage threshold value of the (k + 1) th overvoltage binary table; t is t_i.r.kThe moment when the voltage of the bus i is higher than the voltage threshold of the kth overvoltage binary meter for the first time in the rising process after disturbance; t is t_i.r.k+1The moment when the voltage of the bus i is higher than the voltage threshold of the (k + 1) th overvoltage binary meter for the first time in the rising process after disturbance; t'_i.r.kFor bus i voltage to end up risingThe moment of the first time being equal to the voltage threshold of the kth overvoltage binary table in the later descending process; m is a positive integer; t'_i.r.k+1The time when the voltage of the bus i is equal to the voltage threshold of the (k + 1) th overvoltage binary meter for the first time in the descending process after the ascending is finished; t is t_i.r.mThe moment when the voltage of the bus i is higher than the voltage threshold value of the mth overvoltage binary meter for the first time in the rising process after disturbance; t'_i.r.mThe moment when the voltage of the bus i is equal to the voltage threshold value of the mth overvoltage binary meter for the first time in the descending process after the ascending is finished;

K_r.mis {0 is less than or equal to V_i(t)≤V_cr.r.m+1}∩{t_i.r.m≤t≤t'_i.r.m+1A weight coefficient of the area between;

K_r.kis {0 is less than or equal to V_i(t)≤V_cr.r.k+1}∩{t_i.r.k≤t≤t'_i.r.k+1The weight coefficient of the region between (11) can be solved by the formula;

in the formula (11), k is a positive integer; t is_cr.r.1The time threshold value is the time threshold value in the 1 st overvoltage binary table; t is_cr.r.2Is the time threshold value in the 2 nd overvoltage binary table; t is_cr.r.kIs the time threshold value in the kth overvoltage binary table; t is a unit of_cr.r.k+1The time threshold value is the time threshold value in the (k + 1) th overvoltage binary table; k is_r.1Is {0 ≦ V_i(t)≤V_cr.r.2}∩{t_i.r.1≤t≤t'_i.r.2The weight coefficient of the area between; k_r.kIs {0 is less than or equal to V_i(t)≤V_cr.r.k+1}∩{t_i.r.k≤t≤t'_i.r.k+1The weight coefficient of the area between; v_cr.r.1The voltage threshold value of the 1 st overvoltage binary table; v_cr.r.2Is the voltage threshold of the 2 nd overvoltage binary table; v_cr.r.kA voltage threshold value of a kth overvoltage binary table; v_cr.r.k+1The voltage threshold value of the (k + 1) th overvoltage binary table;

in the formula (12), e is a natural constant; t is time; omega is a kurtosis parameter and is used for representing the steepness degree of the function;

area voltage qualification rate index:

in formula (13), P_aThe voltage qualification rate indexes of the region a behind all N load buses of the whole system are considered, wherein the S operation modes of the system, T typical wind power scenes, M preselected fault sets and N load buses of the whole system are considered; p_lIs the probability that the system operates in the operating mode l; n is a radical of_a.l.v.b.dIn order to increase the number of voltage-qualified buses in the area a when a d-type fault occurs at the load bus b in the operating mode and the wind power scene v, when eta of the bus i_i<When 1, the bus i is considered as a voltage qualified bus; n is a radical of_aThe number of all load buses in the area a; p_WT.vThe probability that wind power operates in a scene v; delta_l.v.dIn the running mode and the wind power scene v, the weight coefficient of the fault d is numerically equal to the probability of the fault d, and typical faults are independent of each other, so that the fault d has

For the probability of the bus b having the fault d, it is assumed that the probability of the fault d occurring at each bus of the system is equal, i.e. the probability of the fault d occurring at each bus of the system is equal

Regional voltage stability margin index:

in the formula (14), eta_aThe method comprises the following steps of considering S operation modes of the system, T typical wind power scenes, M preselected fault sets and voltage stability margin indexes of all N load bus rear areas a of the whole system; eta_i.l.v.d.bThe transient voltage safety margin index of a bus i in an area a is considered when a bus b has a fault d after the system runs in a mode I and a wind power scene v; p_lThe probability that the system operates in the operation mode l; p is_WT.vThe probability that wind power operates in a scene v; delta_l.v.dThe weight coefficient of the fault d is under the condition of the operation mode and the wind power scene v;

5. The reactive power planning method for the receiving-end power grid containing the high-permeability wind power considering the transient voltage stability according to claim 1, wherein the method comprises the following steps: in the step 2, the point distribution method of the dynamic reactive power compensation equipment includes the following steps:

B_i.zs＝λ_iη'_i.risk (15)；

is a weight coefficient, satisfies

S_iActual injected power for bus i; s_baseIs a reference value of the system power;

the bus i is characterized by the transmission or distribution power of the bus i in the system, and the larger the value of the bus i is, the more the bus i transmits or distributes power in the system, the greater the importance degree in the system is;

η'_i.riska voltage instability risk factor for bus i, whose value is defined by equation (17):

in the formula (17), N_l.v.d.iThe total number of the voltage instability buses in the dynamic subarea is determined when the bus i has a fault d in the operation mode and the wind power scene v; eta 'of'_l.v.g.dIn the l operation mode and the wind power scene v, when a bus i has a fault d, the transient voltage safety margin of a voltage instability bus g in a dynamic subarea is defined as eta'_l.v.g.d>1, the transient voltage of the bus g is unstable; p_lIs the probability that the system operates in the operating mode l; p_WT.vThe probability that wind power operates in a scene v; delta_l.v.dThe weight coefficient of the fault d is under the condition of the operation mode and the wind power scene v;

the probability that the bus i has a fault d in the operation mode and the wind power scene v is shown; eta 'of'_i.riskAfter the bus i breaks down, the average size of safety margin indexes of all the instability buses in the dynamic partition is reflected, and the larger the value of the safety margin indexes is, the larger the influence degree of the bus on the fault is;

P_bus＝{i|sort{B_i.zs},i∈{1,2,...,N}} (18)；

in formula (18), sort { B_i.zsIs that each bus is according to B_i.zsThe values are sorted in descending order to form a set; b is_i.zsIs the pivot value of the bus i; i is the serial number of the bus; 1,2, wherein N is the serial number of N load buses in the system;

step 2.2: constructing a candidate bus set to be compensated by formulas (19) to (23);

in the formula (19), SI1.i is a sensitivity index of the bus i based on the transient voltage safety margin of the bus; n is a radical of_l.v.d.iThe total number of the voltage instability buses in the dynamic subarea is determined when the bus i has a fault d in the operation mode and the wind power scene v; eta_g0.dBefore a dynamic reactive power compensation device is installed at a bus i, transient voltage safety margin indexes of a bus g when a system fails d after n low-voltage multi-binary tables and m overvoltage binary tables are considered; eta_g.dAfter a certain dynamic reactive power compensation device is installed at a bus i, when a system has a fault d, the transient voltage safety margin indexes of a bus g under the constraint of n low-voltage multi-binary tables and m overvoltage binary tables are obtained; delta Q_c.iIs the capacity of the dynamic reactive power compensation equipment installed at the bus i;

if during the transient process, a higher requirement is made on the fast voltage support capability of the dynamic reactive power compensation equipment, the system dynamic reactive power backup of the dynamic reactive power compensation equipment i defined by the formula (20) can be used for further correcting the SI1. i;

Q_RTSi＝∫e^-t(Q_i-Q_i0)dt (20)；

in the formula (20), Q_RTSiThe method is characterized in that the method is a quantitative value for representing the dynamic reactive power reserve of the dynamic reactive power compensation equipment connected to a bus i; q_iThe reactive power is increased in the transient process of the dynamic reactive power compensation equipment connected to the bus i; q_i0For steady-state operation, the dynamic reactive compensation connected to the bus iInitial reactive power of the device; e.g. of the type^-tThe dynamic reactive compensation equipment is used for quantizing the reactive power increased and generated by the reactive compensation equipment at each moment in the transient process, and the faster the reactive power increased and generated by the dynamic reactive compensation equipment is, the more beneficial the transient voltage stability of the system is;

based on the formula (19) and the formula (20), defining a sensitivity index of the transient voltage safety margin and the dynamic reactive response rate based on the bus shown in the formula (21);

in the formula (21), si2.i is a sensitivity index based on the transient voltage safety margin and the dynamic reactive response rate of the bus; max { Q }_RTSiI ∈ {1,2, …, N } } is for each Q_RTSiThe largest of the values; SI1.i is a sensitivity index of the bus i based on the transient voltage safety margin of the bus;

further constructing an SI1. i-based candidate busbar set as shown in formula (22) and an SI2. i-based candidate busbar set as shown in formula (23);

C_SI1.bus＝{i|sort{SI1.i},i∈{1,2,...,N}} (22)；

in the formula (22), C_SI1.busIs a candidate bus set based on SI1. i; sort { SI1.i } is a set formed by arranging buses in descending order according to the SI1.i value;

C_SI2.bus＝{i|sort{SI2.i},i∈{1,2,...,N}} (23)；

6. The reactive power planning method for the receiving-end power grid containing high-permeability wind power considering transient voltage stabilization according to claim 1 is characterized in that: in the step 3, the objective function of the differentiated dynamic reactive power compensation optimization model is shown in formulas (24) to (26):

f₁(x)＝ω₁[(1-P_a)+η_a] (24)；

in the formula (24), P_aThe voltage qualification rate index of the area a; eta_aThe voltage stability margin index of the area a is obtained;

min f＝{f₁(x),f₂(x)} (26)；

in formulae (24) to (26), f₁(x) And f₂(x) Respectively representing the dynamic reactive power compensation effect and the dynamic reactive power compensation economic cost for two sub-target functions to be optimized; omega₁、ω₂The optimization weight of the sub-objective function is satisfied₁+ω₂＝1；1-P_aThe voltage instability rate index of the system is obtained; t is a unit of₂The operating life of the SVC; c_svc.uThe reactive compensation unit price of the SVC is achieved; q_svc.uIs the reactive compensation capacity of SVC; f_svc.uThe installation cost of SVC; t is₁The running age of the STATCOM; c_STATCOM.hThe reactive compensation unit price is STATCOM; q_STATCOM.hThe reactive compensation capacity of the STATCOM is obtained; f_STATCOM.hThe installation cost of the STATCOM; z_l.vThe number of compensation nodes of the SVC is installed in the l operation mode and the wind power scene v; h_l.vThe number of compensation nodes for installing the STATCOM in the running mode and the wind power scene v is increased; c_eIs the electricity price; ζ is the annual maximum load hours; p_l.v.lossThe network loss of the system is determined under the condition of the operation mode and the wind power scene v; minf denotes that the two sub-objective functions are minimal.

7. The reactive power planning method for the receiving-end power grid containing the high-permeability wind power considering the transient voltage stability according to claim 1, wherein the method comprises the following steps: in the step 4, the improved entropy weight TOPSIS method is specifically as follows:

for s evaluation schemes, w evaluation indices, whose evaluation steps are shown below:

step (1) of applying an extreme method to set the evaluation matrix A to [ a ═ a-_xy]_s×wStandard of meritThe normalized evaluation matrix A ' ═ a ' was obtained '_xy]_s×wWherein: a'_xyThe values are as follows:

respectively corresponding to the evaluation index y and a_xyMaximum and minimum values of;

In the formula (28), the reaction mixture is,

S_yare respectively a ' shown as a ' and 30 '_xyMean and standard deviation of;

In formula (32), E_yIs the information entropy; s is the total number of evaluation scenarios; x is the current evaluation scheme; a'_xyIs the new value of the index; ln is a logarithmic function with a natural constant as a base number;

R＝(W_y×a'_xy)_s×w＝(r_xy)_s×w (35)；

in formula (35), R is a weighting matrix; w is a group of_yIs the final weight of the index; a'_xyIs the new value of the index; r is_xyThe index value in the weighting matrix; s is the total number of evaluation scenarios; w is the total number of evaluation indexes;

And worst case scenario

In the formula (36), the reaction mixture is,

is an index r_xyMaximum value of (1);

in the formula (37), the reaction mixture is,

is the worst scheme; r is_xyThe index value in the weighting matrix; x is the current evaluation scheme; y is the current evaluation index; w is the total number of evaluation indexes;

is an index r_xyMinimum value of (d);

step (8), calculating each scheme to be evaluated and the optimal scheme in sequence

And the worst scheme

European distance between

In the formula (38), the reaction mixture is,

for each solution to be evaluated and the optimal solution

is an optimal scheme; r is_xyThe index value in the weighting matrix;

in the formula (39), the compound represented by the formula (I),

for each solution to be evaluated and the worst solution

is the worst scheme; r is_xyThe index value in the weighting matrix;

8. The reactive power planning method for the receiving-end power grid containing the high-permeability wind power considering the transient voltage stability according to claim 1, wherein the method comprises the following steps: in step 4, the final configuration scheme is determined by formula (41):

economic cost of the optimal scheme screened out under the operation mode and the wind power scene v; l is the operation mode of the system; s is the total number of the system operation modes; v is the running scene of wind power; t is the total number of typical wind power scenes.

9. Differentiation dynamic reactive power compensation optimization model, its characterized in that: the objective function of this model is shown in equations (24) to (26):

f₁(x)＝ω₁[(1-P_a)+η_a] (24)；

in the formula (24), P_aThe voltage qualified rate index of the area a is obtained; eta_aThe voltage stability margin index of the area a is obtained;

min f＝{f₁(x),f₂(x)} (26)；

in formulae (24) to (26), f₁(x) And f₂(x) Respectively representing the dynamic reactive power compensation effect and the dynamic reactive power compensation economic cost for two sub-objective functions to be optimized; omega₁、ω₂The optimization weight of the sub-objective function is satisfied₁+ω₂＝1；1-P_aThe voltage instability rate index of the system is obtained; t is₂The operating life of the SVC; c_svc.uThe reactive compensation unit price of the SVC is achieved; q_svc.uIs the reactive compensation capacity of SVC; f_svc.uThe installation cost of SVC; t is₁The running age of the STATCOM; c_STATCOM.hThe reactive compensation unit price is STATCOM; q_STATCOM.hThe reactive compensation capacity of the STATCOM is obtained; f_STATCOM.hThe installation cost of the STATCOM; z_l.vThe number of compensation nodes for installing SVC under the operation mode l and the wind power scene v; h_l.vThe number of compensation nodes for installing the STATCOM in the running mode and the wind power scene v is increased; c_eIs the electricity price; ζ is the maximum annual load hours; p_l.v.lossThe network loss of the system is determined under the condition of the operation mode and the wind power scene v; minf denotes that the two sub-targeting functions are minimal.

10. Improved entropy weighted TOPSIS process characterized by: the method evaluates s evaluation schemes and w evaluation indexes, and comprises the following steps:

step 1), applying an extreme value method to enable an evaluation matrix A to be [ a ]_xy]_s×wNormalized to obtain a normalized evaluation matrix A '═ a'_xy]_s×wWherein a'_xyThe values are as follows:

step 2), calculating the coefficient of variation V according to a formula (28)_y；

In the formula (28), the reaction mixture is,

S_yare respectively a ' shown as a ' and 30 '_xyMean and standard deviation of;

step 3), calculating the coefficient of variation V_yWeighted value W of each index_1y；

Step 4) of'_xyCalculating the information entropy E of the discrete distribution situation_y；

In the formula (32), E_yIs the information entropy; s is the total number of evaluation scenarios; x is the current evaluation scheme; a'_xyIs the new value of the index; ln is a logarithmic function with a natural constant as a base number;

step 5), calculating the entropy E based on the information_yWeighted value W of each index_2y；

Step 6), calculating the final weight of each indexHeavy W_yConstructing a weighting matrix R of the method;

R＝(W_y×a'_xy)_s×w＝(r_xy)_s×w (35)；

step 7), determining the optimal scheme under each index by the weighting matrix R

And worst case scenario

In the formula (36), the reaction mixture is,

is an optimal scheme; r is a radical of hydrogen_xyThe index value in the weighting matrix; x is the current evaluation scheme; y is the current evaluation index; w is the total number of evaluation indexes;

is an index r_xyMaximum value of (1);

in the formula (37), the reaction mixture is,

the best scheme is adopted; r is a radical of hydrogen_xyThe index value in the weighting matrix; x is the current evaluation scheme; y is the current evaluation index; w is the total number of evaluation indexes;