Reactive Power Planning Method for Receiving-End Power Grid with High-Penetration Wind Power Considering Transient Voltage Stability

Abstract

A reactive power planning method for a receiving-end power grid with high-penetration wind power considering transient voltage stability, comprising the following steps: Based on multiple binary table criteria, a transient voltage safety margin index is proposed to evaluate the transient voltage stability of the system; Based on the transient voltage safety margin index, a method for locating dynamic reactive power compensation devices is proposed; A differential dynamic reactive power compensation optimization model is established; Through the improved entropy weight TOPSIS method, the best solutions in each scenario are screened out, and then the final configuration plan of the dynamic reactive power compensation devices in the system is determined. The method of the present invention can guide the configuration of dynamic reactive power compensation devices more comprehensively and comprehensively, while ensuring the stable operation of the power grid, maximizing the saving of economic costs.

Description

Technical Field

[0001] The present invention relates to the technical field of reactive power planning for power grids, and particularly relates to a method for reactive power planning of a receiving-end power grid with high-penetration wind power considering transient voltage stability. Background Art

[0002] As the main body of new energy power generation, the penetration rate of wind power in the power grid is increasing year by year. At the same time, the resulting transient voltage problems are becoming increasingly prominent. Reactive power planning is the mainstream way to solve the voltage stability of the power grid. However, at present, most of it is oriented to the static voltage problems of the power grid, aiming to ensure the optimal economic investment while reducing the system power loss and voltage deviation, and there are few planning methods focusing on transient voltage stability. Therefore, it is necessary to establish a set of scientific and feasible reactive power planning methods to improve the robustness of the power grid after being affected by large disturbances.

[0003] At present, the reactive power planning methods considering the transient voltage stability of the power grid generally have the following three types of problems:

[0004] 1) The types of dynamic reactive power sources in the planning are single, and it is impossible to take into account the economy while ensuring that the system can obtain the best reactive power compensation effect;

[0005] 2) The research backgrounds are mostly limited to traditional power grids, which do not conform to the current power grid development pattern in China;

[0006] 3) For the receiving-end power grid with DG (distributed generation), the existing reactive power planning methods do not consider the influence of DG uncertainty on the planning.

[0007] 4) Due to the lack of a standard that can measure the transient voltage stability margin in detail, most reactive power planning only focuses on the voltage instability problem caused by voltage sag, while the overvoltage problem caused by voltage swell is rarely considered.

[0008] In view of the fact that site selection and capacity determination are indispensable core links in reactive power planning, many scholars have respectively carried out research on them, but there are still many problems.

[0009] In the method of locating dynamic reactive power compensation equipment, the literature (Huang Hongyang, Yang Fenyan, Xu Zheng, et al. An improved trajectory sensitivity index-based dynamic reactive power optimization configuration method [J]. Power System Technology, 2012, 36(2): 88-94. DOI: 10.13335 / j.1000-3673.pst.2012.02.020. ISSN: 1000-3673, CN: 11-2410 / TM) defined an index based on an improved trajectory sensitivity, but this index does not truly reflect the trajectory of the additional reactive power of the dynamic reactive power compensation device and does not accurately show the support effect on the bus voltage. The literature (Huang Xiaoqing, Ruan Chicheng, Zou Jiaxin, et al. A dynamic reactive power optimization configuration method considering grid characteristics [J]. Electric 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) considered the reactive power increment of the reactive power compensation device on this basis, adding the actual physical meaning of this index and making its calculation results more accurate. However, this index has relatively strict requirements for the sampling step. If set improperly, it is easy to cause the index to have no solution. In addition, a large amount of simulation data shows that compared with the reactive power-voltage sensitivity index, the trajectory sensitivity index based on the transient voltage stability margin has significant advantages in locating. The literature (Cui Jiehao. Research on transient voltage characteristics evaluation and control of AC-DC power grids [D]. Beijing: North China Electric Power University, 2020. DOI: 10.27140 / d.cnki.ghbbu.2020.000733.) used the trajectory sensitivity index based on the transient voltage stability margin as the basis for locating, but its locating ideas lack completeness. For the rating strategy of dynamic reactive power compensation equipment, the literature (Suo Zhiwen, Liu Jianqin, Jiang Weiyong, et al. Research on the configuration of synchronous condensers in large-scale new energy DC transmission systems [J]. Electric Power Automation Equipment, 2019, 39(9): 124-129. DOI: 10.16081 / j.epae.201909054. ISSN: 1006-6047, CN: 32-1318 / TM) explored the configuration scheme of synchronous condensers in large-scale new energy DC transmission systems and concluded that installing synchronous condensers separately on the low-voltage buses of each new energy collection station can effectively solve the transient overvoltage problem of the sending-end system. However, this literature did not specifically optimize the installed capacity and access quantity of each synchronous condenser, and its economy is poor.

[0010] Literature (Zhou Shihao, Tang Fei, Liu Dichen, etc. Method for configuring dynamic reactive power compensation considering reducing commutation failure risk of multi-infeed HVDC systems [J]. High Voltage Engineering, 2018, 44(10): 3258-3265. DOI: 10.13336 / j.1003-6520.hve.20180925016. ISSN: 1003-6520, CN: 42-1239 / TM) quantitatively optimized the configuration capacity of dynamic reactive power compensation devices based on a decomposition-based multi-objective optimization algorithm, but it was for reducing the commutation failure risk of multi-infeed HVDC power grids. SUMMARY OF THE INVENTION

[0011] In view of the deficiencies of existing power grid reactive power planning methods, the present invention provides a reactive power planning method for a receiving-end power grid with high-penetration wind power considering transient voltage stability, which can more comprehensively and thoroughly guide the configuration of dynamic reactive power compensation devices and maximize cost savings while ensuring the stable operation of the power grid.

[0012] The technical solution adopted by the present invention is as follows:

[0013] A reactive power planning method for a receiving-end power grid with high-penetration wind power considering transient voltage stability, comprising the following steps:

[0014] Step 1: Based on a multi-binary table criterion, a transient voltage security margin index is proposed to evaluate the transient voltage stability of the system;

[0015] Step 2: Based on the transient voltage security margin index, a method for locating dynamic reactive power compensation devices is proposed;

[0016] Step 3: Establish a differential dynamic reactive power compensation optimization model;

[0017] Step 4: Through an improved entropy weight TOPSIS method, the best solutions in each scenario are screened out, and then the final configuration scheme of dynamic reactive power compensation devices in the system is determined.

[0018] In the said Step 1, before proposing the transient voltage security margin index, first based on the wind power scenario probability theory, a typical wind power operation scenario is constructed, as shown in formulas (1) to (6):

[0019]

[0020] In formula (1), f(v) is the probability density function of wind speed; v is the wind speed magnitude; k is the shape parameter of the wind speed distribution; c is the scale parameter, and its value can be calculated from the wind speed mean μ and standard deviation σ; e is the natural constant.

[0021]

[0022] In formula (2), P wis the output power of the wind turbine; P r is the rated capacity of the wind turbine; v ci and v r and v co are the cut-in wind speed, rated wind speed, and cut-out wind speed of the wind turbine respectively;

[0023]

[0024]

[0025]

[0026] In equations (3) to (5), P 1 and P 2 and P 3 are the probabilities of the wind turbine operating in the zero-output scenario, under-output scenario, and rated-output scenario respectively.

[0027]

[0028] is the magnitude of the output power of the wind turbine in the under-output scenario.

[0029] Taking the under-output scenario as an example, equation (6) gives the calculation method of the wind power output. Similarly, the wind power outputs in the other two typical scenarios can be calculated in turn.

[0030] Wind power output in the zero-output scenario:

[0031]

[0032] is the magnitude of the output power of the wind turbine in the zero-output scenario.

[0033] Wind power output in the rated-output scenario:

[0034]

[0035] is the magnitude of the output power of the wind turbine in the rated-output scenario.

[0036] In step 1, the multi-binary table criterion is an evaluation criterion used to measure whether the bus voltages in the power system are unstable when a large disturbance occurs in the power system. By setting the voltage binary tables (V cr.1 , T cr.1 ), (V cr.2 , T cr.2 ), …, (V cr.n , T cr.n ) to require the bus voltage Vi Lower than each preset threshold value V cr.1 、V cr.2 、…、V cr.n The longest duration T b Does not exceed the specified time T cr.1 、T cr.2 、…、T cr.n respectively. When the transient voltage of a certain bus meets this condition, it is considered that the transient voltage of the bus is stable; otherwise, the transient voltage is unstable.

[0037] In the said step 1, the transient voltage safety margin index includes:

[0038] The transient voltage safety margin index of the bus:

[0039]

[0040] In formula (7), σ is the step factor shown in formula (12); η i Is the transient voltage safety margin index of bus i under the constraints of n low-voltage multi-binary tables and m over-voltage binary tables; η i.d.n Is the voltage sag safety margin index of bus i under the constraints of n low-voltage multi-binary tables, as shown in formula (8) specifically; η i.r.m Is the voltage swell index of bus i under the constraints of m over-voltage multi-binary tables, as shown in formula (10) specifically:

[0041]

[0042] In formula (8), V iN Is the rated voltage reference value; t k And t' k Are respectively the moments when the voltage of bus i is lower than the voltage threshold V in the low-voltage binary table during the voltage drop process cr.d.k And higher than the voltage threshold V during the recovery process cr.d.k ; t k+1 And t' k+1 Are respectively the moments when the voltage of bus i is lower than the voltage threshold V in the low-voltage binary table during the voltage drop process cr.d.k+1 And higher than the voltage threshold V during the recovery process cr.d.k+1 ; V i (t) is the transient voltage response curve of bus i; K n Is the weight coefficient of the nth threshold interval.

[0043] K k Is the weight coefficient of each threshold interval, which can be solved step by step by formula (9).

[0044]

[0045] In formula (9), T cr.d.1 is the time threshold of the first low-voltage binary table; T cr.d.2 is the time threshold of the second low-voltage binary table; T cr.d.n is the time threshold of the nth low-voltage binary table; V cr.d.1 is the voltage threshold of the first low-voltage binary table; V cr.d.2 is the voltage threshold of the second low-voltage binary table; V cr.d.n is the voltage threshold of the nth low-voltage binary table; K 1 is the weight coefficient of the first threshold interval; K 2 is the weight coefficient of the second threshold interval; K n is the weight coefficient of the nth threshold interval.

[0046] η i.r.m is the voltage sag index of bus i under the constraint of m over-voltage multi-binary tables, as shown in formula (10) specifically:

[0047]

[0048] In formula (10), V cr.r.k+1 is the voltage threshold of the (k + 1)th over-voltage binary table; t i.r.k is the moment when the voltage of bus i first exceeds the voltage threshold of the kth over-voltage binary table during the rising process after the disturbance; t i.r.k+1 is the moment when the voltage of bus i first exceeds the voltage threshold of the (k + 1)th over-voltage binary table during the rising process after the disturbance; t' i.r.k is the moment when the voltage of bus i first equals the voltage threshold of the kth over-voltage binary table during the falling process after the rising ends; m is a positive integer; t' i.r.k+1 is the moment when the voltage of bus i first equals the voltage threshold of the (k + 1)th over-voltage binary table during the falling process after the rising ends; t i.r.m is the moment when the voltage of bus i first exceeds the voltage threshold of the mth over-voltage binary table during the rising process after the disturbance; t' i.r.m is the moment when the voltage of bus i first equals the voltage threshold of the mth over-voltage binary table during the falling process after the rising ends; K r.m is the weight coefficient of the region between {0 ≤ V i (t) ≤ V cr.r.m+1} ∩ {t i.r.m ≤ t ≤ t' i.r.m+1};

[0049] K r.k is the weight coefficient of the region between {0 ≤ V i (t) ≤ V cr.r.k+1} ∩ {t i.r.k ≤ t ≤ t' i.r.k+1The weight coefficient of the region between} can be solved by formula (11).

[0050]

[0051] In formula (11), k is a positive integer; T cr.r.1 is the time threshold in the first overvoltage binary table; T cr.r.2 is the time threshold in the second overvoltage binary table; T cr.r.k is the time threshold in the k-th overvoltage binary table; T cr.r.k+1 is the time threshold in the (k + 1)-th overvoltage binary table; K r.1 is the weight coefficient of the region between {0 ≤ V i (t) ≤ V cr.r.2} ∩ {t i.r.1 ≤ t ≤ t' i.r.2}; K r.k is the weight coefficient of the region between {0 ≤ V i (t) ≤ V cr.r.k+1} ∩ {t i.r.k ≤ t ≤ t' i.r.k+1}; V cr.r.1 is the voltage threshold of the first overvoltage binary table; V cr.r.2 is the voltage threshold of the second overvoltage binary table; V cr.r.k is the voltage threshold of the k-th overvoltage binary table; V cr.r.k+1 is the voltage threshold of the (k + 1)-th overvoltage binary table.

[0052]

[0053] In formula (12), e is the natural constant; t is the time; ω is the kurtosis parameter, used to characterize the steepness of the function. In the present invention, ω takes 1000.

[0054] Regional voltage qualification rate index:

[0055]

[0056] In formula (13), P a is the voltage qualification rate index of region a considering S operating modes of the system, T typical wind power scenarios, M preselected fault sets, and all N load buses of the entire system; P l is the probability that the system operates in operating mode l; N a.l.v.b.d is the number of voltage qualified buses in region a when a fault of type d occurs at load bus b under operating mode l and wind power scenario v. When η i < 1 for bus i, bus i can be considered a voltage qualified bus; N a is the number of all load buses in region a; PWT.v The probability of wind power operating in scenario v; δ l.v.d is the weight coefficient of fault d under the l operating mode and the wind power scenario v, and numerically equals the probability of its occurrence. Since each typical fault is independent, we have is the probability of fault d occurring at bus b. Here, it is assumed that the probability of fault d occurring at each bus in the system is equal, that is

[0057] Regional voltage stability margin index:

[0058]

[0059] In formula (14), η a is the voltage stability margin index of region a considering S operating modes of the system, T typical wind power scenarios, M preselected fault sets, and all N load buses in the entire system; η i.l.v.d.b is the transient voltage safety margin index of bus i in region a when fault d occurs at bus b considering the system operating in the l operating mode and the wind power scenario v; P l is the probability of the system operating in operating mode l; P WT.v is the probability of wind power operating in scenario v; δ l.v.d is the weight coefficient of fault d under the l operating mode and the wind power scenario v; is the probability of fault d occurring at bus b under the l operating mode and the wind power scenario v.

[0060] In step 1, reactive power planning is a means to improve the stability of the power grid, ensure its safe and economic operation by determining the optimal installation location and capacity of reactive power compensation equipment on the basis of power grid planning.

[0061] In step 2, the method for distributing dynamic reactive power compensation equipment includes the following steps:

[0062] Step 2.1: Screen the key buses in the system according to formulas (15) to (18):

[0063] B i.zs = λ i η' i.risk (15);

[0064] In formula (15), B i.zs is the central value of bus i; λ i is the weight coefficient of bus i, and its value is defined by formula (16);

[0065]

[0066] In formula (16), D iis the degree of bus i, reflecting the number of edges associated with bus i; is the weight coefficient, satisfying In the present invention, S i is the actual injected power of bus i; S base is the reference value of the system power, which is taken as 100 MVA in the present invention; characterizes the magnitude of the power transmitted or distributed by bus i in the system. The larger its value, the more power the bus transmits and distributes in the system, and the greater its importance in the system.

[0067] η' i.risk is the voltage instability risk factor of bus i, and its value is defined by formula (17):

[0068]

[0069] In formula (17), N l.v.d.i is the total number of voltage instability buses in the dynamic partition of bus i when a fault d occurs in bus i under the l operating mode and the wind power scenario v; η' l.v.g.d is the transient voltage safety margin of the voltage instability bus g in the dynamic partition of bus i when a fault d occurs in bus i under the l operating mode and the wind power scenario v. When η' l.v.g.d > 1, it can be considered that the transient voltage of bus g is unstable; P l is the probability of the system operating in the l operating mode; P WT.v is the probability of the wind power operating in the v scenario; δ l.v.d is the weight coefficient of the fault d under the l operating mode and the wind power scenario v; is the probability of a fault d occurring in bus i under the l operating mode and the wind power scenario v; η' i.risk reflects the average magnitude of the safety margin index of all unstable buses in the dynamic partition of bus i after a fault occurs in bus i. The larger its value, the greater the degree of influence of the fault on the bus.

[0070] By calculating the central values of each bus and sorting them in descending order according to formula (18), the key bus set of the system is determined:

[0071] P bus ={i|sort{B i.zs},i∈{1,2,...,N}} (18);

[0072] In formula (18), sort{B i.zs} is the set composed of each bus sorted in descending order according to the value of B i.zs ; B i.zsis the central value of bus i; i is the bus number; 1, 2,..., N are the numbers of N load buses in the system.

[0073] Step 2.2: Construct the candidate bus set to be compensated according to Formulas (19) to (23).

[0074]

[0075] In Formula (19), SI1.i is the sensitivity index of bus i based on the transient voltage security margin of the bus; N l.v.d.i is the total number of voltage instability buses in the dynamic partition when a fault d occurs at bus i under the l operation mode and wind power scenario v; η g0.d is the transient voltage security margin index of bus g when a fault d occurs in the system after considering the constraints of n low-voltage multi-binary tables and m over-voltage binary tables before installing the dynamic reactive power compensation device at bus i; η g.d is the transient voltage security margin index of bus g under the constraints of n low-voltage multi-binary tables and m over-voltage binary tables when a fault d occurs in the system after installing a certain dynamic reactive power compensation device at bus i; △Q c.i is the capacity of the dynamic reactive power compensation device installed at bus i.

[0076] If higher requirements are put forward for the fast voltage support ability of the dynamic reactive power compensation device during the transient process, the system dynamic reactive power reserve of the dynamic reactive power compensation device i defined by Formula (20) can be used to further correct SI1.i.

[0077] Q RTSi = ∫e -t (Q i -Q i0 )dt (20);

[0078] In Formula (20), Q RTSi is the quantization value representing the dynamic reactive power reserve of the dynamic reactive power compensation device connected to bus i; Q i is the reactive power incremented by the dynamic reactive power compensation device connected to bus i during the transient process; Q i0 is the initial reactive power of the dynamic reactive power compensation device connected to bus i during steady-state operation; e -t is the introduced attenuation factor used to quantify the reactive power incremented by the reactive power compensation device at each moment during the transient process. The faster the reactive power compensation device increments reactive power, the more beneficial it is to the transient voltage stability of the system.

[0079] Based on Formulas (19) and (20), define the sensitivity index based on the transient voltage security margin of the bus and the dynamic reactive power response rate as shown in Formula (21).

[0080]

[0081] In Equation (21), SI2.i is the sensitivity index based on the transient voltage security margin and dynamic reactive power response rate of the bus; max{Q RTSi , i∈{1, 2, …, N}} is the maximum value among each Q RTSi value; SI1.i is the sensitivity index of the transient voltage security margin of bus i based on the bus.

[0082] Furthermore, a candidate bus set based on SI1.i as shown in Equation (22) and a candidate bus set based on SI2.i as shown in Equation (23) are constructed.

[0083] C SI1.bus = {i|sort{SI1.i}, i∈{1, 2,..., N}} (22);

[0084] In Equation (22), C SI1.bus is the candidate bus set based on SI1.i; sort{SI1.i} is the set formed by arranging each bus in descending order according to the value of SI1.i.

[0085] C SI2.bus = {i|sort{SI2.i}, i∈{1, 2,..., N}} (23);

[0086] In Equation (23), C SI2.bus is the candidate bus set based on SI2.i; sort{SI2.i} is the set formed by arranging each bus in descending order according to the value of SI2.i.

[0087] In Step 3, the objective function of the differential dynamic reactive power compensation optimization model is as shown in Equations (24) to (26):

[0088] f 1 (x) = ω 1 [(1 - P a ) + η a (24);

[0089] In Equation (24), P a is the voltage qualification rate index of area a; η a is the voltage stability margin index of area a.

[0090]

[0091] min f = {f 1 (x), f 2 (x)} (26);

[0092] In formulas (24) to (26), f 1 (x) and f 2 (x) are two sub-objective functions to be optimized, respectively representing the dynamic reactive power compensation effect and the economic cost of dynamic reactive power compensation; ω 1 、ω 2 are the optimization weights of the sub-objective functions, satisfying ω 1 +ω 2 =1; 1 - P a is the system voltage instability rate index; T 2 is the operation life of the SVC; C svc.u is the unit price of reactive power compensation of the SVC; Q svc.u is the reactive power compensation capacity of the SVC; F svc.u is the installation cost of the SVC; T 1 is the operation life of the STATCOM; C STATCOM.h is the unit price of reactive power compensation of the STATCOM; Q STATCOM.h is the reactive power compensation capacity of the STATCOM; F STATCOM.h is the installation cost of the STATCOM; Z l.v is the number of compensation nodes for installing the SVC in the l operation mode and the wind power scenario v; H l.v is the number of compensation nodes for installing the STATCOM in the l operation mode and the wind power scenario v; C e is the electricity price; ζ is the annual maximum load hours; P l.v.loss is the network loss of the system in the l operation mode and the wind power scenario v; minf represents the minimum of the two sub-objective functions.

[0093] In step 3, the differential dynamic reactive power compensation is reflected in: 1) For different candidate nodes, the types of dynamic reactive power compensation devices installed are not the same; 2) A dual optimization strategy is adopted to implement reactive power planning. The dual optimization strategy is an optimization idea that respectively considers two situations of the best reactive power compensation effect and the optimal economic cost. When considering that the system can obtain the best reactive power compensation effect during the transient process, let the optimization weight ω 1 : ω 2 =2:1, that is, ω 1 =0.67, ω 2 =0.33; while when focusing on the optimal economic cost, then let ω 1 : ω 2 =1:2, that is, ω 1 =0.33, ω 2 =0.67.

[0094] In 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. For s evaluation schemes and w evaluation indicators, its evaluation steps are as follows:

[0095] Step (1): Use the extreme value method to standardize the evaluation matrix A = [a xy s×w to obtain the standardized evaluation matrix A' = [a' xy s×w , where: a' xy takes the following values:

[0096]

[0097] In formula (27), x is the current evaluation scheme; y is the current evaluation indicator; a xy is the initial value of the indicator; a' xy is the new value of the indicator; correspond to the maximum and minimum values of a xy under the evaluation indicator y, respectively.

[0098] Step (2): Calculate the coefficient of variation V y according to formula (28);

[0099]

[0100] In formula (28), S y are the mean and standard deviation of a' xy shown in formulas (29) and (30), respectively.

[0101]

[0102]

[0103] Step (3): Calculate the weight value W y of each indicator based on the coefficient of variation V 1y ;

[0104]

[0105] Step (4): Calculate its information entropy E xy from the discrete distribution of a' y ;

[0106]

[0107] In formula (32), E y is the information entropy; s is the total number of evaluation schemes; x is the current evaluation scheme; a' xy ​​is the new value of the index; ln is the natural logarithm function with the natural constant as the base.

[0108] Step (5), calculate the weight value W of each index based on the information entropy E y under 2y ;

[0109]

[0110] Step (6), calculate the final weight W of each index y , and construct its weighted matrix R;

[0111]

[0112] R = (W y × a') xy ) s×w = (r xy ) s×w (35);

[0113] In formula (35), R is the weighted matrix; W y is the final weight of the index; a' xy is the new value of the index; r xy is the index value in the weighted matrix; s is the total number of evaluation schemes; w is the total number of evaluation indicators.

[0114] Step (7), determine the optimal scheme and the worst scheme

[0115]

[0116] In formula (36), is the optimal scheme; r xy is the index value in the weighted matrix; x is the current evaluation scheme; y is the current evaluation indicator; w is the total number of evaluation indicators; is the maximum value in the index r xy .

[0117]

[0118] In formula (37), is the worst scheme; r xy is the index value in the weighted matrix; x is the current evaluation scheme; y is the current evaluation indicator; w is the total number of evaluation indicators; is the minimum value in the index r xy .

[0119] Step (8), calculate the difference between each scheme to be evaluated and the optimal scheme and the worst scheme Euclidean distance between

[0120]

[0121] In formula (38), is the Euclidean distance between each scheme to be evaluated and the optimal scheme ; y is the current evaluation index; w is the total number of evaluation indexes; is the optimal scheme; r xy is the index value in the weighted matrix.

[0122]

[0123] In formula (39), is the Euclidean distance between each scheme to be evaluated and the worst scheme ; y is the current evaluation index; w is the total number of evaluation indexes; is the worst scheme; r xy is the index value in the weighted matrix.

[0124] Step (9), calculate the closeness degree F between each scheme to be evaluated and the ideal scheme x ;

[0125]

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

[0127]

[0128] In formula (41), PRO fin is the final configuration scheme of the dynamic reactive power compensation equipment; is the economic cost of the best scheme screened under the l operation mode and the wind power scenario v; l is the operation mode of the system; S is the total number of system operation modes; v is the operation scenario of the wind power; T is the total number of typical wind power scenarios.

[0129] For a reactive power planning method for a receiving-end power grid with high-penetration wind power considering transient voltage stability according to the present invention, the technical effects are as follows:

[0130] 1) For a receiving-end power grid with high-penetration new energy, when performing reactive power planning, it is necessary to consider the influence of DG uncertainty on it. The present invention proposes a configuration method for differential dynamic reactive power compensation equipment, constructs a set of transient voltage safety margin indexes that can guide the reactive power planning of the power grid, and then proposes a corresponding point-placement method to reasonably determine the configuration scheme of the dynamic reactive power compensation equipment in the system.

[0131] 2) The present invention incorporates the objective factor of the uncertainty of wind power into the reactive power planning, making the planning results more realistic and reasonable.

[0132] 3) The transient voltage security margin index based on multiple binary table criteria proposed by the present invention can not only accurately evaluate the transient voltage stability of the system, but also provide guidance for reactive power planning.

[0133] 4) The reactive power planning method of the present invention can not only ensure the reactive power compensation effect, but also save economic costs to the greatest extent. Description of the Drawings

[0134] Figure 1 It is a diagram of the IEEE 39-bus system.

[0135] Figure 2 It is a diagram of the comparison results of sensitivity indexes.

[0136] Figure 3 It is a diagram of the voltage waveform of the fault bus after compensation.

[0137] Figure 4 It is a diagram of the comparison results of SI2.i and SI1.i.

[0138] Figure 5 It is a diagram of the voltage waveforms of some key buses of the system.

[0139] Figure 6(a) is a diagram of the voltage waveforms of all buses in the network under various working conditions (uncompensated);

[0140] Figure 6(b) is a diagram of the voltage waveforms of all buses in the network under various working conditions (wind power zero output scenario);

[0141] Figure 6(c) is a diagram of the voltage waveforms of all buses in the network under various working conditions (wind power under-output scenario);

[0142] Figure 6(d) is a diagram of the voltage waveforms of all buses in the network under various working conditions (wind power rated output scenario);

[0143] Figure 7 It is a flowchart of the method of the present invention. Detailed Embodiments

[0144] As Figure 7 shown, the method for reactive power planning of a receiving-end power grid with high-penetration wind power considering transient voltage stability includes the following steps:

[0145] Step 1: Based on multiple binary table criteria, a transient voltage security margin index is proposed to evaluate the transient voltage stability of the system.

[0146] Before presenting the transient voltage security margin index, based on the probability theory of wind power scenarios, the wind power operation scenarios are reduced to three typical scenarios: zero-output scenario, under-output scenario, and rated-output scenario. Then, calculate the operation probability and output power of the wind turbine in each scenario. The specific process can be seen in formulas (1) to (6):

[0147]

[0148] In formula (1), f(v) is the probability density function of wind speed; v is the wind speed magnitude; k is the shape parameter of the wind speed distribution; c is the scale parameter, and its value can be calculated from the mean wind speed μ and standard deviation σ; e is the natural constant.

[0149]

[0150] In formula (2), P w is the output power of the wind turbine; P r is the rated capacity of the wind turbine; v ci , v r , v co are the cut-in wind speed, rated wind speed, and cut-out wind speed of the wind turbine respectively.

[0151]

[0152]

[0153]

[0154] In formulas (3) to (5), P 1 , P 2 , P 3 are the probabilities of the wind turbine operating in the zero-output scenario, under-output scenario, and rated-output scenario respectively.

[0155]

[0156] Taking the under-output scenario as an example, formula (6) gives the calculation method of wind power output. Similarly, the wind power output in the other two typical scenarios can be calculated in turn.

[0157] Wind power output in the zero-output scenario:

[0158]

[0159] is the output power magnitude of the wind turbine in the zero-output scenario.

[0160] Wind power output in the rated-output scenario:

[0161]

[0162] is the output power of the wind turbine under the rated output scenario.

[0163] Construct the bus transient voltage safety margin index, regional voltage qualification rate index, and regional voltage stability margin index. The specific process is shown in formulas (7) to (14) as follows:

[0164]

[0165] In formula (7), where η i is the transient voltage safety margin index of bus i under the constraints of n low-voltage multi-binary tables and m over-voltage binary tables; η i.d.n is the voltage sag safety margin index of bus i under the constraints of n low-voltage multi-binary tables, as shown in formula (8) specifically; η i.r.m is the voltage swell index of bus i under the constraints of m over-voltage multi-binary tables, as shown in formula (10) specifically; σ is the step factor shown in formula (12).

[0166]

[0167] In formula (8), V iN is the rated voltage reference value; t k and t' k are the moments when the voltage of bus i is lower than the voltage threshold V cr.d.k in the low-voltage binary table during the voltage drop process and higher than the voltage threshold V cr.d.k during the recovery process respectively; K k is the weight coefficient within each threshold interval, which can be solved step by step by formula (9); t k+1 and t' k+1 are the moments when the voltage of bus i is lower than the voltage threshold V cr.d.k+1 in the low-voltage binary table during the voltage drop process and higher than the voltage threshold V cr.d.k+1 during the recovery process respectively; V i (t) is the transient voltage response curve of bus i; K n is the weight coefficient of the nth threshold interval.

[0168]

[0169] In formula (9), T cr.d.1 is the time threshold of the first low-voltage binary table; T cr.d.2 is the time threshold of the second low-voltage binary table; T cr.d.n is the time threshold of the nth low-voltage binary table; V cr.d.1 is the voltage threshold of the first low-voltage binary table; V cr.d.2is the voltage threshold of the second low-voltage binary table; V cr.d.n is the voltage threshold of the nth low-voltage binary table; K 1 is the weight coefficient of the first threshold interval; K 2 is the weight coefficient of the second threshold interval; K n is the weight coefficient of the nth threshold interval.

[0170]

[0171] In Equation (10), V cr.r.k+1 is the voltage threshold of the (k + 1)th overvoltage binary table; t i.r.k is the moment when the voltage of bus i first exceeds the voltage threshold of the kth overvoltage binary table during the rising process after the disturbance; t i.r.k+1 is the moment when the voltage of bus i first exceeds the voltage threshold of the (k + 1)th overvoltage binary table during the rising process after the disturbance; t' i.r.k is the moment when the voltage of bus i first equals the voltage threshold of the kth overvoltage binary table during the falling process after the rising ends; m is a positive integer; t' i.r.k+1 is the moment when the voltage of bus i first equals the voltage threshold of the (k + 1)th overvoltage binary table during the falling process after the rising ends; t i.r.m is the moment when the voltage of bus i first exceeds the voltage threshold of the mth overvoltage binary table during the rising process after the disturbance; t' i.r.m is the moment when the voltage of bus i first equals the voltage threshold of the mth overvoltage binary table during the falling process after the rising ends; K r.m is for {0 ≤ V i (t) ≤ V cr.r.m+1} ∩ {t i.r.m ≤ t ≤ t' i.r.m+1} the weight coefficient of the area between; K r.k is for {0 ≤ V i (t) ≤ V cr.r.k+1} ∩ {t i.r.k ≤ t ≤ t' i.r.k+1} the weight coefficient of the area between, and its numerical value can be solved by Equation (11).

[0172]

[0173] In Equation (11), k is a positive integer; T cr.r.1 is the time threshold in the first overvoltage binary table; T cr.r.2 is the time threshold in the second overvoltage binary table; T cr.r.k is the time threshold in the kth overvoltage binary table; T cr.r.k+1 is the time threshold in the (k + 1)th overvoltage binary table; K r.1 is for {0 ≤ V i(t) ≤ V cr.r.2} ∩ {t i.r.1 ≤ t ≤ t' i.r.2}; the weight coefficient of the region between; K r.k is {0 ≤ V i (t) ≤ V cr.r.k+1} ∩ {t i.r.k ≤ t ≤ t' i.r.k+1}; the weight coefficient of the region between; V cr.r.1 is the voltage threshold of the first overvoltage binary table; V cr.r.2 is the voltage threshold of the second overvoltage binary table; V cr.r.k is the voltage threshold of the k-th overvoltage binary table; V cr.r.k+1 is the voltage threshold of the (k + 1)-th overvoltage binary table.

[0174]

[0175] In Equation (12), e is the natural constant; t is time; ω is the kurtosis parameter, used to characterize the steepness of the function. In the present invention, ω takes 1000.

[0176]

[0177] In Equation (13), P a is the voltage qualification rate index of region a considering the S operating modes of the system, T typical wind power scenarios, M preselected fault sets, and all N load buses of the entire system; P l is the probability that the system operates in operating mode l; N a.l.v.b.d is the number of voltage qualified buses in region a when a fault of type d occurs at load bus b under operating mode l and wind power scenario v. When η i < 1 for bus i, it can be considered that bus i is a voltage qualified bus; N a is the total number of load buses in region a; P WT.v is the probability that wind power operates in scenario v; δ l.v.d is the weight coefficient of fault d under operating mode l and wind power scenario v, and numerically equals its occurrence probability. Since the typical faults are independent of each other, so is the probability that a fault of d occurs at bus b under operating mode l and wind power scenario v. Here, it is assumed that the probability of fault d occurring at each bus of the system is equal, that is

[0178]

[0179] In Equation (14), η ais the voltage stability margin index of area a considering S operating modes of the system, T typical wind power scenarios, M preselected fault sets, and all N load buses in the entire system; η i.l.v.d.b is the transient voltage security margin index of bus i in area a when a fault d occurs at bus b considering the system in operating mode l and wind power scenario v; P l is the probability of the system operating in operating mode l; P WT.v is the probability of the wind power operating in scenario v; δ l.v.d is the weight coefficient of fault d under operating mode l and wind power scenario v; is the probability of a fault d occurring at bus b under operating mode l and wind power scenario v.

[0180] Step 2: Based on the transient voltage security margin index, propose a method for placing dynamic reactive power compensation devices.

[0181] Before placing the devices, key buses in the system need to be screened according to formulas (15) to (18).

[0182] B i.zs = λ i η' i.risk (15);

[0183] In formula (15), B i.zs is the central value of the bus; λ i is the weight coefficient of bus i, and its value is defined by formula (16); η' i.risk is the voltage instability risk factor of bus i, and its value is defined by formula (17).

[0184]

[0185] In formula (16), D i is the degree of bus i, reflecting the number of edges associated with bus i; is the weight coefficient, satisfying In the present invention S i is the actual injection power of bus i; S base is the reference value of the system power, which is taken as 100 MVA in the present invention; characterizes the magnitude of the power transmitted or distributed by bus i in the system. The larger its value, the more power the bus transmits and distributes in the system, and the greater its importance in the system.

[0186]

[0187] In formula (17), N l.v.d.i is the total number of voltage instability buses in the dynamic partition of bus i when a fault d occurs at bus i under operating mode l and wind power scenario v; η'l.v.g.d When a fault d occurs on bus i under the l operation mode and the wind power scenario v, the transient voltage safety margin of the voltage-unstable bus g in its dynamic partition. When η' l.v.g.d > 1, it can be considered that the transient voltage of bus g is unstable; P l is the probability that the system operates in the l operation mode; P WT.v is the probability that the wind power operates in the v scenario; δ l.v.d is the weight coefficient of the fault d under the l operation mode and the wind power scenario v; is the probability that a fault d occurs on bus i under the l operation mode and the wind power scenario v; η' i.risk reflects the average size of the safety margin indexes of all unstable buses in the dynamic partition of bus i after a fault occurs on it. The larger its value, the greater the degree of influence of the fault on this bus.

[0188] By calculating the central values of each bus and sorting them in descending order according to formula (18), the key bus set of the system is determined:

[0189] P bus = {i|sort{B i.zs}, i ∈ {1, 2,..., N}} (18);

[0190] In formula (18), sort{B i.zs} is the set composed of each bus sorted in descending order according to the value of B i.zs ; B i.zs is the central value of bus i; i is the bus number; 1, 2,..., N are the numbers of N load buses in the system.

[0191] After that, the candidate bus set to be compensated is constructed by formulas (19) to (23).

[0192]

[0193] In formula (19), SI1.i is the sensitivity index of bus i based on the transient voltage safety margin of the bus; N l.v.d.i is the total number of voltage-unstable buses in the dynamic partition of bus i when a fault d occurs under the l operation mode and the wind power scenario v; η g0.d is the transient voltage safety margin index of bus g when a fault d occurs in the system after considering the constraints of n low-voltage multi-binary tables and m over-voltage binary tables before installing the dynamic reactive power compensation device at bus i; η g.d is the transient voltage safety margin index of bus g under the constraints of n low-voltage multi-binary tables and m over-voltage binary tables when a fault d occurs in the system after installing a certain dynamic reactive power compensation device at bus i; △Q c.i is the capacity of the dynamic reactive power compensation equipment installed at bus i.

[0194] When higher requirements are put forward for the fast voltage support ability of dynamic reactive power compensation equipment during the transient process, the system dynamic reactive power reserve of dynamic reactive power compensation equipment i defined by formula (20) can be used to further correct SI1.i.

[0195] Q RTS i = ∫e -t (Q i -Q i0 )dt (20);

[0196] In formula (20), Q RTSi is the quantization value representing the dynamic reactive power reserve of the dynamic reactive power compensation equipment connected to bus i; Q i is the reactive power incremented by the dynamic reactive power compensation equipment connected to bus i during the transient process; Q i0 is the initial reactive power of the dynamic reactive power compensation equipment connected to bus i during steady-state operation; e -t is the introduced attenuation factor, which is used to quantify the reactive power incremented by the reactive power compensation equipment at each moment during the transient process. The faster the dynamic reactive power compensation equipment increments reactive power, the more beneficial it is to the transient voltage stability of the system.

[0197] Based on formula (19) and formula (20), a sensitivity index of the transient voltage security margin and dynamic reactive power response rate based on the bus is defined as shown in formula (21).

[0198]

[0199] In formula (21), SI2.i is the sensitivity index of the transient voltage security margin and dynamic reactive power response rate based on the bus; max{Q RTSi , i ∈ {1, 2,..., N}} is the maximum value among each Q RTSi value; SI1.i is the sensitivity index of the transient voltage security margin based on bus i.

[0200] Furthermore, a candidate bus set based on SI1.i as shown in formula (22) and a candidate bus set based on SI2.i as shown in formula (23) are constructed.

[0201] C SI1.bus = {i|sort{SI1.i}, i ∈ {1, 2,..., N}} (22);

[0202] In formula (22), C SI1.bus is the candidate bus set based on SI1.i; sort{SI1.i} is the set composed of each bus sorted in descending order according to the value of SI1.i.

[0203] CSI2.bus ={i|sort{SI2.i}, i ∈ {1, 2, ..., N}} (23);

[0204] In formula (23), C SI2.bus is the candidate bus set based on SI2.i; sort{SI2.i} is the set composed of each bus arranged in descending order according to the value of SI2.i.

[0205] Step 3: Establish a differential dynamic reactive power compensation optimization model and solve the Pareto optimal solution set under each scenario through the multi-objective grey wolf optimizer (MOGWO).

[0206] Establish the objective function of the differential dynamic reactive power compensation optimization model as shown in formulas (24) - (26):

[0207] f 1 (x) = ω 1 [(1 - P a ) + η a (24);

[0208] In formula (24), P a is the voltage qualification rate index of area a; η a is the voltage stability margin index of area a.

[0209]

[0210] min f = {f 1 (x), f 2 (x)} (26);

[0211] In formulas (24) - (26), f 1 (x) and f 2 (x) are two sub-objective functions to be optimized, respectively representing the dynamic reactive power compensation effect and the economic cost of dynamic reactive power compensation; ω 1 , ω 2 are the optimization weights of the sub-objective functions, satisfying ω 1 + ω 2 = 1; 1 - P a is the system voltage instability rate index; T 2 is the operation years of SVC; C svc.u is the reactive power compensation unit price of SVC; Q svc.u is the reactive power compensation capacity of SVC; F svc.u is the installation cost of SVC; T 1 is the operation years of STATCOM; C STATCOM.h is the reactive power compensation unit price of STATCOM; Q STATCOM.h is the reactive power compensation capacity of STATCOM; FSTATCOM.h is the installation cost of STATCOM; Z l.v is the number of compensation nodes for installing SVC under the l operation mode and wind power scenario v; H l.v is the number of compensation nodes for installing STATCOM under the l operation mode and wind power scenario v; C e is the electricity price; ζ is the annual maximum load hours; P l.v.loss is the network loss of the system under the l operation mode and wind power scenario v; minf represents the minimum of two sub-objective functions.

[0212] The main steps to solve the Pareto optimal solution set in each scenario by the multi-objective grey wolf algorithm (MOGWO) are as follows:

[0213] 1) Set the grey wolf population size N corresponding to the capacity of the dynamic reactive power compensation device p , the maximum number of iterations I, the dimension D, the upper and lower limits U b and L b of the capacity of the dynamic reactive power compensation device, the number N Ar of the external population Archive, and the relevant control parameters;

[0214] 2) Randomly generate grey wolf individuals that meet the conditions, calculate the objective function values of each individual, determine the current optimal individual, and update Archive;

[0215] 3) Select three leading wolves α, β, and δ from Archive based on the roulette method, update the positions of the remaining grey wolves in turn by the positions of the leading wolves, and calculate their corresponding objective function values;

[0216] 4) According to the objective function values of each grey wolf, find the new non-dominated individuals and update Archive;

[0217] 5) Repeat steps 3) and 4) until the iteration reaches I times. If the number of the external population reaches the upper limit during this process, randomly eliminate some individuals from the crowded group according to the grouping results of the objective functions of each individual in Archive until the number of the external population Archive is equal to N Ar ;

[0218] 6) Output the positions of the grey wolves in Archive, and this position represents the solved Pareto optimal solution set.

[0219] Step 4: Screen out the best solutions in each scenario through the improved entropy weight method for evaluating the distance from the ideal solution, and then determine the final configuration plan of the dynamic reactive power compensation device in the system.

[0220] Evaluate s evaluation schemes and w evaluation indicators by using the improved entropy weight TOPSIS method. The steps are shown in formulas (27) to (40):

[0221] Step (1): Use the extreme value method to standardize the evaluation matrix A = [a xy s×w to obtain the standardized evaluation matrix A' = [a' xy s×w , where the value of a' xy is as follows:

[0222]

[0223] In Equation (27), x is the current evaluation scheme; y is the current evaluation index; a xy is the initial value of the index; a' xy is the new value of the index; respectively correspond to the maximum and minimum values of a xy under the evaluation index y.

[0224] Step (2): Calculate the coefficient of variation V y according to Equation (28);

[0225]

[0226] In Equation (28), S y are respectively the mean and standard deviation of a' xy as shown in Equation (29) and Equation (30).

[0227]

[0228]

[0229] Step (3): Calculate the weight value W y of each index based on the coefficient of variation V 1y ;

[0230]

[0231] Step (4): Calculate its information entropy E xy from the discrete distribution of a' y ;

[0232]

[0233] In Equation (32), E y is the information entropy; s is the total number of evaluation schemes; x is the current evaluation scheme; a' xy is the new value of the index; ln is the logarithmic function with the natural constant as the base.

[0234] Step (5): Calculate based on the information entropy E y ​​The weight value W of each of the following indicators 2y ;

[0235]

[0236] Step (6), calculate the final weight W of each indicator y , and construct its weighted matrix R;

[0237]

[0238] R = (W y × a' xy ) s×w = (r xy ) s×w (35);

[0239] In formula (35), R is the weighted matrix; W y is the final weight of the indicator; a' xy is the new value of the indicator; r xy is the indicator value in the weighted matrix; s is the total number of evaluation schemes; w is the total number of evaluation indicators.

[0240] Step (7), determine the optimal scheme and the worst scheme

[0241]

[0242] In formula (36), is the optimal scheme; r xy is the indicator value in the weighted matrix; x is the current evaluation scheme; y is the current evaluation indicator; w is the total number of evaluation indicators; is the maximum value in the indicator r xy .

[0243]

[0244] In formula (37), is the worst scheme; r xy is the indicator value in the weighted matrix; x is the current evaluation scheme; y is the current evaluation indicator; w is the total number of evaluation indicators; is the minimum value in the indicator r xy .

[0245] Step (8), calculate the Euclidean distance between each scheme to be evaluated and the optimal scheme and the worst scheme

[0246]

[0247] In formula (38), is the Euclidean distance between each scheme to be evaluated and the optimal scheme ; y is the current evaluation index; w is the total number of evaluation indexes; is the optimal scheme; r xy is the index value in the weighted matrix.

[0248]

[0249] In formula (39), is the Euclidean distance between each scheme to be evaluated and the worst scheme ; y is the current evaluation index; w is the total number of evaluation indexes; is the worst scheme; r xy is the index value in the weighted matrix.

[0250] Step (9), calculate the closeness degree F of each scheme to be evaluated and the ideal scheme x ;

[0251]

[0252] Determine the final configuration scheme of the dynamic reactive power compensation device in the system by formula (41):

[0253]

[0254] In formula (41), PRO fin is the final configuration scheme of the dynamic reactive power compensation device; is the economic cost of the best scheme selected under the l operation mode and the wind power scenario v; l is the operation mode of the system; S is the total number of system operation modes; v is the operation scenario of wind power; T is the total number of typical wind power scenarios.

[0255] Example:

[0256] The parameters related to the present invention are as follows:

[0257] 1) The present invention sets the following low-voltage multi-binary table and over-voltage binary table:

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

[0259] 2) Regarding the entire system as the area to be studied, that is, N a = N = 39;

[0260] 3) When performing reactive power planning, only the normal operation mode of the power system is taken as an example, and the operating conditions of wind power under three typical scenarios are considered, that is, S = 1; T = 3; P WT.1 = 0.1561; P WT.2 = 0.6837; P WT.3 = 0.1602;

[0261] 4) A preselected fault set is composed of three-phase short-circuit faults and single-phase ground faults, and it is assumed that the probabilities of each fault occurring at various locations in the system are equal, that is, M = 2; δ l.v.1 = 0.07; δ l.v.2 = 0.93;

[0262] 5) From the network topology structure and power flow information shown by Figure 1 , the weight coefficients of each bus calculated by formula (16) are shown in Table 1:

[0263] Table 1 Weight coefficients of each bus

[0264]

[0265] 6) According to the results of the system N-1 fault scan, the voltage instability risk factors of each bus calculated by formula (17) are shown in Table 2:

[0266] Table 2 Voltage instability risk factors of each bus

[0267]

[0268] 7) The economic parameters and investment costs of the dynamic reactive power compensation equipment are shown in Table 3. In addition, C e = 0.617 yuan / (kW·h); ζ = 3600h;

[0269] Table 3 Economic parameters and investment costs of the dynamic reactive power compensation equipment

[0270]

[0271] 8) U i.min , U i.max are taken as 0.95 p.u. and 1.05 p.u. respectively; Q i.min = 0; Q i.maxIt is determined by the type of dynamic reactive power compensation equipment. According to the typical model recommended by the Large Electric Machine Subcommittee of the Chinese Society for Electrical Engineering (Liu Zhenya, Zhang Qiping, Wang Yating, et al. Research on reactive power compensation measures to improve the safety and stability level of the 750 kV sending-end power grid in the new Gansu-Qinghai region of Northwest China [J]. Proceedings of the CSEE, 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 State Grid Corporation of China regarding STATCOM (State Grid Corporation of China. Q / GDW 241.1—2008 Chain-type Static Synchronous Compensator Part 1: Functional Specification Guide [S]. Beijing: State Grid Corporation of China, 2008. DL / T 1215.1-2013), the reliability of STATCOM equipment with a capacity of 300 Mvar and above remains to be verified. Therefore, this invention focuses on STATCOM, Q i.max is set to 300 Mvar; while for SVC, Q i.max is taken as 200 Mvar;

[0272] 9) In the MOGWO algorithm, N p = 100; I = 50, D = 3 and 4; for STATCOM, U b = 300, for SVC, U b = 200; L b = 0; N Ar = 20.

[0273] Based on the combined simulation platform of PSD-BPA and MATLAB, the Figure 1 receiving-end power grid with high-penetration wind power shown is simulated and analyzed to verify the rationality and effectiveness of the proposed reactive power planning method. The system N-1 fault scanning results show that the three-phase short-circuit fault at the head of line 16-17 has the most serious impact on the key buses and all the buses in the power grid. Therefore, this invention sets a three-phase short-circuit fault at the head of line 16-17 in the 5th cycle, the circuit breakers at the heads of the lines in the subsequent 4.5 cycles act, and the circuit breaker at the end of the line in the 5th cycle acts, based on which the best access position of dynamic reactive power compensation is explored.

[0274] In this invention, 100 Mvar SVC is selected for the placement points. By connecting it to each bus in turn, its voltage support for the key buses and all the buses in the power grid is investigated. Since buses 30 to 39 are generator buses and their voltage amplitudes are constant, these buses can be not considered during the placement points. To demonstrate the superiority of the proposed placement point method, the sensitivity index based on the transient voltage safety margin of the bus in Equation (15) is now compared with the existing method 1 and the existing method 2, and the comparison results are as Figure 2 shown.

[0275] Analysis Figure 2 It can be seen that the sensitivity ranking results of each bus obtained by the three methods are not exactly the same. The candidate nodes screened by the method of the present invention are bus 22, bus 21, bus 23... The best candidate nodes selected by the existing method 1 and the existing method 2 are bus 15, bus 14, bus 23... and bus 8, bus 26, bus 28... respectively. When the 100 Mvar SVC is sequentially connected to bus 22, bus 15, and bus 8, the fault bus voltage waveforms and the system transient voltage security indexes are respectively as Figure 3 shown in Table 4.

[0276] Table 4 System Transient Voltage Security Indexes

[0277]

[0278]

[0279] From Figure 3 Table 4, it can be seen that when the dynamic reactive power compensation equipment is arranged based on the method of the present invention, during the transient process, better compensation effects can be obtained for both the fault bus and all buses in the system, highlighting the advantages of the method of the present invention. In addition, in the constructed candidate bus set, the candidate buses that can provide fast voltage support can be further selected by Equation (17), Figure 4 and the sensitivity indexes of each bus based on Equation (15) and Equation (17) are compared.

[0280] According to Figure 4 the magnitudes of SI2.i of each bus in SI2.bus , the candidate bus set C composed of bus 19, bus 23, bus 22... is constructed by Equation (19) SI2.bus . Figure 5 The voltage waveforms of some key buses in the system during the transient process are compared when the SVC is connected to the bus with the highest ranking in C SI1.bus and C SI2.bus respectively. From Figure 5 it can be seen that compared with C SI1.bus , when the buses in C SI2.bus are selected for bus placement, the key buses of the system can show better transient voltage characteristics in the early stage of the transient process, which verifies the rationality and correctness of the sensitivity indexes of the present invention based on the transient voltage security margin of the bus and the dynamic reactive power response rate.

[0281] For the buses with the highest ranking in the candidate bus set, the present invention installs a STATCOM on this bus to support the voltage of all network buses during the transient process with its fast response ability, and SVCs are installed on the remaining buses in the order of ranking, aiming to focus on the economic cost of reactive power compensation. Tables 5 and 6 are respectively based on C SI1.bus and C SI2.bus, showing the system voltage stability margin index under the preset compensation scheme.

[0282] Table 5 is based on C SI1.bus of the preset compensation scheme

[0283]

[0284] Table 6 is based on C SI2.bus of the preset compensation scheme

[0285]

[0286] Tables 5 and 6 show that when the reactive power compensation capacity is certain, installing the STATCOM on the buses with higher rankings in the candidate bus set is more beneficial to the transient voltage stability of the system, which verifies the scientificity of the reactive power planning strategy of the present invention and lays a foundation for subsequent capacity determination.

[0287] According to the selected C SI1.bus and C SI2.bus , based on the aforementioned dual optimization strategy, the MOGWO algorithm is used to solve the optimal capacity of the dynamic reactive power compensation devices installed on each bus under different scenarios. Table 7 gives the Pareto optimal solution set of the wind power under the rated output scenario, and the evaluation results of each optimal scheme after being evaluated by the improved entropy weight TOPSIS method are shown in Table 8.

[0288] Table 7 Optimal configuration scheme of dynamic reactive power compensation devices under the rated output scenario

[0289]

[0290] Table 8 Evaluation results of each optimal scheme under the rated output scenario

[0291]

[0292] Table 8 shows that the better the reactive power compensation effect, the higher the corresponding economic cost. Considering comprehensively, Scheme 3 with the largest closeness is selected as the best configuration scheme of the dynamic reactive power compensation devices under the wind power rated output scenario.

[0293] The best configuration schemes of the dynamic reactive power compensation devices under the wind power zero output scenario and the underrated output scenario are shown in Table 9.

[0294] Table 9 Configuration schemes of dynamic reactive power compensation devices under different wind power scenarios

[0295]

[0296] To determine the final configuration scheme of the dynamic reactive power compensation devices in the system, Table 10 summarizes the evaluation results of the best schemes for each wind power scenario.

[0297] Evaluation results of the optimal solutions under each wind power scenario in Table 10

[0298]

[0299] Using the data in Table 10, the optimal solution for the wind power in the zero-output scenario is determined by Equation (37) as the final configuration plan of the system dynamic reactive power compensation equipment. Figures 6(a) to 6(d) The transient voltage waveforms of all the network buses under each working condition in the final configuration plan are shown. Figures 6(a) to 6(d) It shows that under the action of the final configuration plan, all the buses of the system can maintain transient voltage stability under each working condition, which highlights the applicability of the fixed-capacity method of the present invention.

Claims

1. A reactive power planning method for a receiving-end power grid with high-penetration wind power considering transient voltage stability, characterized in that it includes the following steps: Step 1: Based on the multi-binary table criterion, a transient voltage security margin index is proposed to evaluate the transient voltage stability of the system; Step 2: Based on the transient voltage security margin index, a method for locating dynamic reactive power compensation devices is proposed; Step 3: Establish a differential dynamic reactive power compensation optimization model; Step 4: Through the improved entropy weight TOPSIS method, the best solutions in each scenario are screened out, and then the final configuration plan of the dynamic reactive power compensation devices in the system is determined; In the said Step 1, before proposing the transient voltage security margin index, first based on the wind power scenario probability theory, typical wind power operation scenarios are constructed as shown in formulas (1) to (6): (1)； In Equation (1), is the probability density function of wind speed; v is the magnitude of wind speed; k is the shape parameter of wind speed distribution; c is the scale parameter, and its value can be calculated from the mean wind speed μ and the standard deviation σ ; e is the natural constant; (2)； In formula (2), P w is the output power of the wind turbine; is the rated capacity of the wind turbine; , , are the cut-in wind speed, rated wind speed, and cut-out wind speed of the wind turbine, respectively; ; ; (3)； (4)； (5)； In formulas (3) to (5), P 1 , P 2 , P 3 are the probabilities of the wind turbine operating in the zero-output scenario, under-output scenario, and rated-output scenario, respectively; (6)； is the output power of the wind turbine under the under-output scenario; Wind power output power in the zero-output scenario: is the output power of the wind turbine under the zero-output scenario; Wind power output power in the rated-output scenario: is the output power of the wind turbine under the rated output scenario; In the said Step 1, the transient voltage security margin index includes: Transient voltage security margin index of the bus: (7)； In formula (7), is the step factor shown in formula (12); is n low-voltage multi-binary tables and m the transient voltage security margin index of the bus under the constraint of i over-voltage binary tables; is n the voltage sag security margin index of the bus under the constraint of i low-voltage multi-binary tables, as specifically shown in formula (8); is m the voltage swell index of the bus under the constraint of i over-voltage multi-binary tables, as specifically shown in formula (10): (8)； In formula (8), is the rated voltage reference value; and are respectively the moments when the busbar i voltage is lower than the voltage threshold in the low-voltage binary table during the voltage dip V cr.d.k and higher than the voltage threshold V cr.d.k during the recovery process; and are respectively the moments when the busbar i voltage is lower than the voltage threshold in the low-voltage binary table during the voltage dip V cr.d.k+1 and higher than the voltage threshold V cr.d.k+1 during the recovery process; is the transient voltage response curve of the busbar i ; is the weight coefficient of the n th threshold interval; is the weight coefficient for each threshold interval and can be solved step by step using formula (9); (9)； In Equation (9), is the time threshold of the first low-voltage binary table; is the time threshold of the second low-voltage binary table; is the n time threshold of the is the voltage threshold of the first low-voltage binary table; is the voltage threshold of the second low-voltage binary table; is the n voltage threshold of the is the weight coefficient of the first threshold interval; is the weight coefficient of the second threshold interval; is the n weight coefficient of the For m the voltage sag index of the bus under the constraint of the overvoltage multi-binary table i is shown in the specific formula (10) as follows: (10)； In formula (10), is the k voltage threshold of the ( t i.r.k +1)-th overvoltage binary table; i is the moment when the k bus voltage is first higher than the voltage threshold of the bus i during the rising process after the disturbance; k is the moment when the t' i.r.k bus voltage is first higher than the voltage threshold of the i +1)-th overvoltage binary table during the rising process after the disturbance; k is the moment when the m bus voltage is first equal to the voltage threshold of the t' i.r.k+1 bus i during the falling process after the rising ends; k is the moment when the t i.r.m bus voltage is first equal to the voltage threshold of the i +1)-th overvoltage binary table during the falling process after the rising ends; m is the moment when the t' i.r.m bus voltage is first higher than the voltage threshold of the i during the rising process after the disturbance; m is the moment when the For the weight coefficient in the region between {0 ≤ V i ( t ) ≤ V cr.r.m+1}} ∩ { t i.r.m ≤ t ≤ t' i.r.m+1} K r.k For {0 ≤ V i ( t ) ≤ V cr.r.k+1}} ∩ { t i.r.k ≤ t ≤ t' i.r.k+1} the weight coefficient of the area in between, the numerical value of which can be solved by formula (11); (11)； In formula (11), k is a positive integer; T cr.r.1 is the time threshold in the first overvoltage binary table; T cr.r.2 is the time threshold in the second overvoltage binary table; T cr.r.k is the k th time threshold in the overvoltage binary table; T cr.r.k+1 is the k +1 th time threshold in the overvoltage binary table; is the weight coefficient for the region between {0 ≤ V i ( t ) ≤ V cr.r.2} ∩ { t i.r.1 ≤ t ≤ t' i.r.2}; K r.k is the weight coefficient for the region between {0 ≤ V i ( t ) ≤ V cr.r.k+1} ∩ { t i.r.k ≤ t ≤ t' i.r.k+1}; is the voltage threshold of the first overvoltage binary table; is the voltage threshold of the second overvoltage binary table; is the voltage threshold of the k th overvoltage binary table; is the voltage threshold of the k +1 th overvoltage binary table; (12)； In Equation (12), e is the natural constant; t is time; ω is the kurtosis parameter, which is used to characterize the steepness of the function; Regional voltage qualification rate index: (13)； In formula (13), is the voltage qualification rate index of area a after considering S system operation modes, T a typical wind power scenario, M a preselected fault set, and all N load buses in the entire system; is the probability that the system operates in operation mode l ; is l under the operation mode and wind power scenario v , the number of voltage qualified buses in area a when a fault of type b occurs at load bus d . When the i of bus is i , the bus N a is considered a voltage qualified bus; is the probability that wind power operates in scenario v ; is l under the operation mode and wind power scenario v , the weight coefficient of fault d , which is numerically equal to its occurrence probability. Since each typical fault is independent, ; is the probability that a fault occurs at bus b . Here, it is assumed that the probability of fault d occurring at each bus in the system is equal, that is d ; ; Regional voltage stability margin index: (14)； In formula (14), is the voltage stability margin index of area a considering S system operation modes, T typical wind power scenarios, M a preselected fault set, and all N load buses in the whole system; is the transient voltage security margin index of the bus in area a when a fault occurs in l the system operation mode, wind power scenario v and the bus b after the d fault occurs; i is the probability that the system operates in operation mode ; l is the probability that the wind power operates in scenario ; v is the weight coefficient of fault under l the system operation mode and wind power scenario v ; d is the probability that a fault occurs in bus under l the system operation mode and wind power scenario v ; b when the fault d occurs; In the said Step 2, the method for locating dynamic reactive power compensation devices includes the following steps: Step 2.1: Screen out the key buses in the system according to formulas (15) to (18): (15)； In formula (15), is the central value of the busbar i ; is the weight coefficient of the busbar i , and its value is defined by formula (16); (16)； In formula (16), is the degree of the busbar i , reflecting the number of edges associated with the busbar i ; , are weight coefficients, satisfying ; S i is the actual injection power of the busbar i ; S base is the reference value of the system power; characterizes the magnitude of the power transmitted or distributed by the busbar i in the system. The larger its value, the more power the busbar transmits and distributes in the system, and the greater its importance in the system; is the voltage instability risk factor of the busbar i and its value is defined by formula (17): (17)； In formula (17), is the total number of voltage - unstable buses in its dynamic partition when l in the operation mode and wind - power scenario v the bus i has a fault d ; is the transient voltage safety margin of the voltage - unstable bus in its dynamic partition when l in the operation mode and wind - power scenario v the bus i has a fault d . When g , it can be considered that the bus has transient voltage instability; g ; is the probability that the system operates in the operation mode l ; is the probability that the wind power operates in the scenario v ; is the weight coefficient of the fault l in the operation mode and wind - power scenario v ; d is the probability that the bus has a fault l in the operation mode and wind - power scenario v ; i has a fault d ; reflects the average size of the safety - margin index of all unstable buses in its dynamic partition after the bus i has a fault. The larger its value, the greater the degree of influence of the fault on the bus. By calculating the central values of each bus and arranging them in descending order according to formula (18), the key bus set of the system is determined: (18)； In formula (18), sort{ B i.zs} is a set formed by arranging each busbar in descending order according to the B i.zs value; is the central value of busbar i ; is the number of the busbar; is the N load busbar numbers in the system; Step 2.2: Construct a set of candidate buses to be compensated from formulas (19) to (23); (19)； In formula (19), is the busbar i sensitivity index of transient voltage security margin based on the busbar; is at l operation mode and wind power scenario v Under the condition, when a fault occurs in busbar i the total number of voltage - unstable busbars in its dynamic partition; d is before installing the dynamic reactive power compensation device at busbar Considering i low - voltage multi - binary tables and n and m constraints of over - voltage binary tables, when a fault occurs in the system d the transient voltage security margin index of busbar g ; is after installing a certain dynamic reactive power compensation device at busbar i When a fault occurs in the system d at this time, n low - voltage multi - binary tables and m constraints of over - voltage binary tables, the transient voltage security margin index of busbar g ; is the capacity of the dynamic reactive power compensation device installed at busbar i ; When higher requirements are put forward for the fast voltage support ability of dynamic reactive power compensation equipment during the transient process, the system dynamic reactive power reserve of the dynamic reactive power compensation equipment defined by formula (20) can be used to i further correct; ​ (20)； In formula (20), Q RTSi is the quantization value of the dynamic reactive power reserve of the dynamic reactive power compensation device connected to the bus i location; Q i is the reactive power increased by the dynamic reactive power compensation device connected to the bus i location during the transient process; Q i0 is the initial reactive power of the dynamic reactive power compensation device connected to the bus i location during steady-state operation; e -t is the introduced attenuation factor, which is used to quantify the reactive power increased by the reactive power compensation device at each moment during the transient process. The faster the dynamic reactive power compensation device increases reactive power, the more beneficial it is to the transient voltage stability of the system. Based on formulas (19) and (20), a sensitivity index based on the transient voltage security margin and dynamic reactive power response rate of the bus as shown in formula (21) is defined; (21)； In formula (21), is the sensitivity index based on the transient voltage security margin of the busbar and the dynamic reactive power response rate; is the largest among the Q RTSi values; is the sensitivity index of the transient voltage security margin of the busbar i based on the busbar; Further construct the candidate bus set shown in formula (22) based on and the candidate bus set shown in formula (23) based on ; (22)； In formula (22), is the candidate bus set based on ; is the set formed by arranging each bus in descending order according to the value. (23)； In formula (23), is the candidate bus set based on ; is the set formed by arranging each bus in descending order according to the value.

2. The reactive power planning method for a receiving-end power grid with high-penetration wind power considering transient voltage stability according to claim 1, characterized in that: In the said step 1, the multi-binary table criterion refers to: by setting voltage binary tables ( V cr.1 , T cr.1 ), ( V cr.2 , T cr.2 ), …, ( V cr.n , T cr.n ) to require that the voltage V i of a certain bus is lower than each preset threshold value V cr.1 , V cr.2 , …, V cr.n and the longest duration T b does not exceed the specified time T cr.1 , T cr.2 , …, T cr.n respectively. When the transient voltage of a certain bus meets this condition, it is considered that the transient voltage of the bus is stable; otherwise, the transient voltage is unstable.

3. The reactive power planning method for a receiving-end power grid with high-penetration wind power considering transient voltage stability according to claim 1, characterized in that: In the said Step 3, the objective function of the differential dynamic reactive power compensation optimization model is as shown in formulas (24) to (26): (24)； In formula (24), is the voltage qualification rate index of area a; is the voltage stability margin index of area a; (25)； (26)； In formulas (24) to (26), and are two sub-objective functions to be optimized, respectively representing the dynamic reactive power compensation effect and the economic cost of dynamic reactive power compensation; ω 1 , ω 2 are the optimization weights of the sub-objective functions, satisfying ω 1 + ω 2 = 1; is the system voltage instability rate index; T 2 is the operating life of the SVC; is the unit price of reactive power compensation of the SVC; is the reactive power compensation capacity of the SVC; is the installation cost of the SVC; T 1 is the operating life of the STATCOM; is the unit price of reactive power compensation of the STATCOM; is the reactive power compensation capacity of the STATCOM; is the installation cost of the STATCOM; is at l operating mode and wind power scenario v the number of compensation nodes for installing the SVC; is at l operating mode and wind power scenario v the number of compensation nodes for installing the STATCOM; C e is the electricity price; is the annual maximum load hours; is at l operating mode and wind power scenario v the network loss of the system; indicates that the two sub-objective functions are minimized.

4. The reactive power planning method for a receiving-end power grid with high-penetration wind power considering transient voltage stability according to claim 1, characterized in that: In the said Step 4, the improved entropy weight TOPSIS method is specifically as follows: For s evaluation scenarios w evaluation metrics, the evaluation steps are as follows: Step (1), use the extreme value method to normalize the evaluation matrix A = s×w to obtain the normalized evaluation matrix A'= s×w , where: takes the following values:​​​ (27)； In formula (27), x is the current evaluation plan; y is the current evaluation index; is the initial value of the index; is the new value of the index; and correspond to the maximum and minimum values of the evaluation index y respectively under ; Step (2), calculate the coefficient of variation according to formula (28) V y ; (28)； In formula (28), and are the average value and the standard deviation of shown in formula (29) and formula (30) respectively; (29)； (30)； Step (3), calculate the weight values of each index based on the coefficient of variation V y ; ; (31)； Step (4), calculate its information entropy from the discrete distribution of ; (32)； In formula (32), is the information entropy; s is the total number of evaluation schemes; x is the current evaluation scheme; is the new value of the index; ln is the logarithmic function with the natural constant as the base; Step (5), calculate the weight values of each index based on information entropy ;​ (33)； Step (6), calculate the final weights of each index , and construct its weighted matrix R ; (34)； (35)； In formula (35), R is the weighting matrix; is the final weight of the index; is the new value of the index; is the index value in the weighting matrix; s is the total number of evaluation schemes; w is the total number of evaluation indicators; Step (7), determine the optimal solution and the worst solution for each index from the weighted matrix R ; and the worst solution ; (36)； In Equation (36), is the optimal solution; is the index value in the weighting matrix; x is the current evaluation solution; y is the current evaluation index; w is the total number of evaluation indices; is the index the maximum value in; (37)； In Equation (37), is the worst solution; is the index value in the weighted matrix; x is the current evaluation solution; y is the current evaluation index; w is the total number of evaluation indices; is the index the minimum value in; Step (8), calculate the Euclidean distances between each solution to be evaluated and the optimal solution and the worst solution respectively , ; (38)； In formula (38), is the Euclidean distance between each scheme to be evaluated and the optimal scheme ; y is the current evaluation index; w is the total number of evaluation indexes; is the optimal scheme; is the index value in the weighting matrix; (39)； In formula (39), is the Euclidean distance between each solution to be evaluated and the worst solution ; y is the current evaluation index; w is the total number of evaluation indices; is the worst solution; is the index value in the weighted matrix; Step (9), calculate the closeness degree of each alternative to be evaluated and the ideal alternative ; (40)。 5. The reactive power planning method for a receiving-end power grid with high-penetration wind power considering transient voltage stability according to claim 1, characterized in that: In the said Step 4, the final configuration plan is determined by formula (41): (41)； In formula (41), is the final configuration scheme of the dynamic reactive power compensation device; is the economic cost of the best scheme screened under l the operation mode and the wind power scenario v ; l is the operation mode of the system; S is the total number of the operation modes of the system; v is the operation scenario of the wind power; T is the total number of the typical wind power scenarios.

Images