Power system uncertainty optimization scheduling method and system considering security risk

Abstract

The invention relates to a power system uncertainty optimization scheduling method and system considering safety risks, and the method comprises the steps: constructing a risk measurement model considering faults and system states at the same time, and enabling the risk measurement model to be used for quantifying and evaluating the safety of a system; constructing a scheduling optimization model, taking minimization of ground state power generation cost, expected power generation cost and expected wind curtailment load shedding cost as a target function, and taking an alternating current power flow constraint, a unit output upper and lower limit constraint, a unit climbing constraint and a system comprehensive risk constraint under an original scene set and a fault as constraint conditions; reducing the original scene set by adopting a scene reduction method based on a Wasserstein distance; solving a scheduling optimization model based on the subtracted optimal typical scene set and the corresponding probability distribution; and performing system scheduling according to a scheduling scheme obtained according to a solving result. The method has a more obvious indication effect on the safety state of the system, meanwhile, the calculation speed of the scheduling optimization model is increased, and the engineering application requirement is greatly met.

Description

Technical Field

[0001] The present invention relates to the technical field of power system optimization and dispatching, and in particular to a power system uncertainty optimization and dispatching method and system taking into account safety risks. Background Art

[0002] Under the background of new power system construction, affected by the uncertainty of renewable energy output and the uncertainty of system component failure, the system state and cross-section power flow change dramatically, which in turn affects or destroys the stable operation of the system. In order to make the economic operation of the power system more stable, it is usually necessary to adopt a corresponding risk model that takes into account the uncertainty of renewable energy output or the uncertainty of system component failure to quantify and evaluate the static safety risk of the system. If the risk value of the system exceeds the limit, it is necessary to carry out corresponding optimization scheduling that takes into account the safety risk so that the risk value of the system is reduced to meet the requirements and prevent problems before they happen. At present, risk scheduling methods can be divided into two categories: preventive control risk scheduling and corrective control risk scheduling. The main idea of ​​preventive control risk scheduling is to adjust the current output plan to make the risk index meet the constraints under various expected scenarios (including renewable energy output and failure scenarios). Its advantage is that the model is simple, but because it does not consider the readjustment of the output plan under each scenario, it is difficult to obtain the optimal solution that meets all constraints. The main idea of ​​corrective control risk scheduling is to introduce measures such as output plan readjustment, load shedding, wind and solar abandonment in a single expected scenario to improve the flexibility of the risk scheduling plan and its reliability under the expected scenario. Compared with preventive and control risk scheduling, the conservatism of the scheduling plan is greatly reduced, and the effectiveness and feasibility of risk scheduling are improved.

[0003] The existing power system uncertainty optimization scheduling taking into account safety risks mainly has the following problems: the measurement of safety risks in the existing risk scheduling process is unreasonable. The risk measurement based on CVaR conditional risk has strict theoretical support, but it involves complex calculations due to integral operations; when the existing risk scheduling method uses random optimization methods to solve, the number of scenarios is large and the calculation speed is slow, which is difficult to meet the requirements of engineering applications. Summary of the invention

[0004] In view of the shortcomings of the prior art, the present invention provides a method and system for optimizing the scheduling of power system uncertainty taking into account safety risks, with the aim of improving the rationality of safety risk measurement while improving the scheduling optimization model to increase the solution speed, reduce the solution difficulty, and meet the requirements of engineering applications.

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

[0006] The present invention provides a power system uncertainty optimization dispatching method taking into account safety risks, comprising:

[0007] Constructing a risk measurement model that takes into account both faults and system status to quantify and evaluate system safety; the system status is determined by the voltage of each node in the system and the transmission power value of each line;

[0008] Construct a dispatch optimization model with the objective function of minimizing the base-state power generation cost, the expected power generation cost, and the expected wind curtailment and load shedding cost, and with the constraints of the original scenario set and the AC power flow constraints under faults, the upper and lower limits of the unit output, the unit ramp constraints, and the system comprehensive risk constraints; the original scenario set includes the base-state scenario before dispatch and the new energy uncertainty scenario set obtained by sampling the new energy output distribution;

[0009] The original scene set is reduced by using a scene reduction method based on the wasserstein distance to obtain an optimal typical scene set and a corresponding probability distribution;

[0010] The scheduling optimization model is solved based on the optimal typical scenario set and the corresponding probability distribution to obtain a scheduling plan, and the system is scheduled according to the scheduling plan.

[0011] Further technical solutions are:

[0012] The construction of the risk measurement model includes:

[0013] Establish the system overall voltage over-limit risk index And the system overall branch power flow over-limit risk index

[0014] in:

[0015] m is the number of faults;

[0016] p s,j is the probability of the jth device N-1 being out of service, i.e., failure j, in the system;

[0017] α V,j is the system voltage over-limit risk index under fault j, N is the number of nodes in the system, α V,i,j is the voltage over-limit risk of node i under fault j:

[0018]

[0019] In the above formula, is the voltage V at node i under fault j i,j The probability density function of V min,i , V max,i are the lower and upper limits of the voltage at node i respectively;

[0020] α P,j is the system branch power flow over-limit risk index under fault j, N l is the total number of system lines, α P,i,j is the over-limit risk of line i under fault j:

[0021]

[0022] In the above formula, represents the transmission power P of line i under fault j i,j The probability density function of max,i represents the maximum transmission power of line i;

[0023] Based on the α V , α P , α V,j , α P,j Obtain the system comprehensive safety risk index α under fault j s,j , and the system comprehensive security risk index α taking into account all system states s :

[0024]

[0025] Among them, the weight parameter σ1, σ2, are greater than or equal to 0,

[0026] The expression of the risk measurement model is discretized, so that when calculating α s,j , α s In the process, the probability density function The integral operation is converted into a discretization operation.

[0027] The objective function of the scheduling optimization model is:

[0028]

[0029] Where T represents the total number of time periods; N g Indicates the number of thermal power generators; N s Indicates the number of scenes in the original scene set; N k Indicates the number of N-1 outages; N d Indicates the number of loads; N w Indicates the number of wind turbines; Δd i,t,k , Δw i,t,k They represent the actual output of the i-th thermal power generator in the k-th N-1 outage base state scenario during period t; C LS , C WS They represent load shedding cost and wind abandonment cost respectively; f i(·) represents the power generation cost function of the i-th thermal power generator with respect to output; Δd i,t,s,k , Δw i,t,s,k They represent the actual output of the i-th thermal power generator, the load shedding amount of the i-th load, and the wind abandonment amount of the i-th wind turbine generator in the t-th period in the s-th scenario under the k-th N-1 outage;

[0030] Among them, the probability of the kth N-1 outage occurring is k=1,...N k ; When k = 0, represents the probability of the ground state scenario; N l is the total number of lines; k 、pro l are the probability of failure of line k and line l respectively.

[0031] The AC power flow constraints of the scheduling optimization model are:

[0032]

[0033] Among them, p ij,t,s,k ,q ij,t,s,k are the active power and reactive power flowing from node i to node j in the sth scenario during period t under the kth N-1 outage; v i,t,s,k is the voltage of node i in time period t in the sth scenario under the kth N-1 outage; θ ij,t,s,k is the voltage phase angle difference between node i and node j in period t in the sth scenario under the kth N-1 outage; and is the reactance and susceptance of the line between node i and node j under the kth N-1 outage; They represent the active and reactive outputs of the i-th thermal power generator in the s-th scenario during the k-th N-1 outage, respectively; They represent the active and reactive outputs of the i-th wind turbine in the s-th scenario during the t-th period under the k-th N-1 outage, respectively; They represent the i-th active load and reactive load in the t-th period in the s-th scenario under the k-th N-1 outage; Δp g,i,t,s,k , Δq g,i,t,s,k They respectively represent the active output adjustment and reactive output adjustment of the i-th thermal power generator in the k-th N-1 outage scenario during period t.

[0034] The system comprehensive risk constraint of the scheduling optimization model is:

[0035]

[0036] The above formula indicates that the total risk of the system in T time periods is less than the set risk threshold TOL, and the weight parameter prob s Represents the probability of scene s appearing in the original scene set;

[0037] Among them, g v (v i,t,s,k ) is the voltage v of node i under the kth N-1 outage in scenario s during period t i,t,s,k The risk value of crossing the limit is:

[0038]

[0039] In the above formula, v min 、v max They are v i,t,s,k The lower and upper limits of

[0040] g l (p l,t,s,k ) is the power flow p of line l under the kth N-1 outage in scenario s during period t l,t,s,k The risk value of crossing the limit is:

[0041]

[0042] In the above formula, -p l,max 、p l,max They are v i,t,s,k The lower and upper limits of .

[0043] The scene reduction method based on the wasserstein distance is used to reduce the original scene set to obtain the optimal typical scene set and the corresponding probability distribution, including:

[0044] Build an optimization problem model:

[0045]

[0046] Where N represents the number of original scene sets before reduction; p i Indicates the corresponding scene The probability of k i A 0-1 variable indicating whether the i-th scene is a classic scene; x i,j A 0-1 variable indicating whether the original scene i is assigned to the classic scene j;

[0047] Determine the number of classic scenarios, including:

[0048] (1) Set the number of typical scenarios K = 2 and give a convergence threshold tol;

[0049] (2) Solve the optimization problem model and find K classic scenarios based on the values ​​of the 0-1 variables And calculate the probability of the corresponding classic scene [prob1,...,prob K ]:

[0050]

[0051] Among them, x j,i A 0-1 variable indicating whether the original scene j is assigned to the classic scene i;

[0052] (3) If (f K-1 -f K ) / f K-1 ≤tol, the calculation ends, and the number of classic scenes K at this time is the selected optimal value; if (f K-1 -f K ) / f K-1 >tol, then K=K+1, jump to step (2), until the convergence condition is met, and output the calculation result.

[0053] Solving the scheduling optimization model based on the optimal typical scenario set and the corresponding probability distribution includes:

[0054] According to the calculation results of solving the optimization problem model, the objective function of the scheduling optimization model and the system comprehensive risk constraint are transformed to obtain a model after scenario reduction, and its objective function is as follows:

[0055]

[0056] Among them, prob s Represents the probability of scene s appearing in the classic scene set after elimination.

[0057] 9. The method according to claim 8, characterized in that solving the model after scene reduction comprises:

[0058] Take λ,λ≥0 to relax the system comprehensive risk constraint to the objective function and obtain the relaxed model:

[0059]

[0060] AC power flow constraints, output upper and lower limit constraints, and ramp constraints under all faults and scenarios

[0061] Among them, wv 1,i,t,s,k 、wv 2,i,t,s,k 、wv 3,i,t,s,k 、wv 4,i,t,s,k , and wp 1,i,t,s,k 、wp 2,i,t,s,k 、wp 3,i,t,s,k 、wp 4,i,t,s,k All of them are continuous variables introduced by piecewise function continuation of the system comprehensive risk constraint;

[0062] For a given λ,λ>0, the relaxed model is decomposed into (N k +1)(N s +1) mutually independent sub-problems are solved separately. The details of each sub-problem are as follows:

[0063]

[0064] st fault k, AC power flow constraints under scenario s, output upper and lower limit constraints, climbing constraints

[0065] The present invention also provides a power system uncertainty optimization dispatching system taking into account safety risks, comprising:

[0066] A risk measurement module, which is used to construct a risk measurement model that takes into account both faults and system status to quantify and evaluate system safety; the system status is determined by the voltage of each node in the system and the transmission power value of each line;

[0067] A scheduling optimization model construction module is used to construct a scheduling optimization model, which takes the minimization of the base state power generation cost, the expected power generation cost, and the expected wind curtailment and load shedding cost as the objective function, and takes the AC power flow constraints, the upper and lower limit constraints of the unit output, the unit climbing constraints, and the system comprehensive risk constraints for the original scenario set and faults as the constraint conditions; the original scenario set includes the base state scenario before scheduling and the new energy uncertainty scenario set obtained by sampling the new energy output distribution;

[0068] A scene reduction module, which uses a scene reduction method based on wasserstein distance to reduce the original scene set to obtain an optimal typical scene set and a corresponding probability distribution;

[0069] A solution module solves the scheduling optimization model based on the optimal typical scenario set and the corresponding probability distribution, obtains a scheduling plan, and schedules the system with the scheduling plan.

[0070] The beneficial effects of the present invention are as follows:

[0071] The present invention provides a new idea for uncertainty optimization dispatch calculation under new power systems, and has the following advantages:

[0072] According to the SRI risk expression, the present invention obtains the node voltage over-limit risk index under the fault and the system voltage over-limit risk index under the fault, and then quantifies the node and system overall voltage over-limit and line flow and system flow overall over-limit risks caused by the uncertainty of new energy output under the fault, calculates the corresponding risks under all faults, and takes the probability weighted sum of the risk values ​​of all faults and the base state as the safety risk that takes into account the uncertainty of the system state and the occurrence of faults, thereby establishing a risk measurement model, which has a more obvious indicative effect on the system safety state and effectively reduces the economic cost of the system.

[0073] The present invention converts the integral operation in the risk measurement model calculation process into a discretization operation, thereby avoiding complex integral operations.

[0074] The constraints of the scheduling optimization model of the present invention include AC power flow constraints under various scenarios and faults, and system comprehensive risk constraints. At the same time, second-order cone relaxation and piecewise function continuation technology are introduced to convert them into convex constraints, which greatly improves the solvability and computational efficiency of the model. In addition, the scenario reduction method based on the wasserstein distance reduces the number of scenarios on the basis of reasonably characterizing the uncertainty of new energy output, solving the problem of unreasonable definition of the scene set spacing in the existing scenario reduction. By defining the distance between scene sets through the supremum of the difference between the objective function before and after the scenario reduction, the change in calculation accuracy before and after the scenario reduction is more intuitive.

[0075] In view of the particularity of the scheduling optimization model after scene reduction, the present invention adopts Lagrangian decomposition and parallel computing technology to quickly solve the model, thereby improving the calculation speed and realizing the rapid solution of the model.

[0076] Other features and advantages of the present invention will be set forth in the following description or may be learned by practicing the present invention. BRIEF DESCRIPTION OF THE DRAWINGS

[0077] Figure 1 Schematic diagram of the process of the embodiment of the present invention.

[0078] Figure 2 These are 200 original wind power scenarios and predicted scenario diagrams provided in the verification experiment of the embodiment of the present invention.

[0079] Figure 3 It is the change of the Wasserstein distance between the original scene and the typical scene when K=2-20 in the verification experiment of the embodiment of the present invention.

[0080] Figure 4 These are typical scenarios and corresponding probabilities after K=15 is eliminated in the verification experiment of the embodiment of the present invention.

[0081] Figure 5These are typical scenarios and corresponding probabilities after K=10 is eliminated in the verification experiment of the embodiment of the present invention;

[0082] Figure 6 These are typical scenarios and corresponding probabilities after K=5 is eliminated in the verification experiment of the embodiment of the present invention.

[0083] Figure 7 It is a topology diagram of the IEEE39 node system in the verification experiment of the embodiment of the present invention. DETAILED DESCRIPTION

[0084] The specific implementation of the present invention is described below with reference to the accompanying drawings.

[0085] See also Figure 1 , a power system uncertainty optimization scheduling method taking into account security risks in this embodiment includes:

[0086] S1. Construct a risk measurement model that takes into account both faults and system status to quantify and evaluate system safety; the system status is determined by the voltage of each node in the system and the transmission power value of each line.

[0087] In this embodiment, it is considered that there are m devices in the system (only lines and units are considered) and N-1 outages occur, and the probability of the jth device N-1 outage, that is, the probability of fault j occurring is p j , then according to the Synthesis Risk Index (SRI) expression, the node voltage over-limit (branch power flow over-limit) risk index under fault j and the system voltage over-limit (branch power flow over-limit) risk index under the fault can be obtained, and then the node and system overall voltage over-limit (line power flow and system power flow over-limit) risk caused by the uncertainty of new energy output under fault j can be quantified. Calculate the corresponding risks under all m faults, that is, take the probability weighted sum of the risk values ​​of all faults and the base state as the safety risk that takes into account the uncertainty of system state and fault occurrence. The branch includes lines and transformers; the node includes busbars and units and loads on the busbars.

[0088] Specifically, the construction of the risk measurement model includes:

[0089] For a power grid consisting of N nodes (i.e., busbars and the units and loads on the busbars), there are a total of m devices (i.e., lines and units) in the system that have N-1 outages. j (j=1,...,m) represents the failure probability of the jth device, then the system base state probability p s,0 : The probability p of the jth device N-1 being out of service, i.e., failure j, in the system s,j for:

[0090] According to the SRI risk expression, the voltage over-limit risk α of node i under fault j is V,i,j for:

[0091]

[0092] Among them, V i,j is the voltage of node i under fault j, is its probability density function; V min,i , V max,i Represent the lower and upper limits of the voltage at node i respectively. The comprehensive voltage over-limit risk of node i It is the weighted sum of the probability of voltage over-limit risk of node i in each state. This indicator comprehensively measures the over-limit risk of node i. The larger its value means the greater the over-limit risk of the node: Accordingly, the system voltage over-limit risk index under fault j is the sum of the voltage over-limit risks of all nodes in the system, that is: N is the number of nodes in the system. Considering m faults, the overall voltage over-limit risk index of the system is α V Then it is the weighted sum of the base state probability and fault probability of the system voltage over-limit risk under the base state and each fault:

[0093] The branch power flow over-limit risk index and system comprehensive safety risk index under fault conditions include:

[0094] The calculation of branch power flow over-limit risk index is similar to that of node voltage over-limit risk. The over-limit risk index of line i under fault j is:

[0095]

[0096] Among them, P max,i represents the maximum transmission power of line i; represents the probability density function of the transmission power of line i under fault j. If the total number of system lines is N l , the system branch power flow over-limit risk index under fault j is: Considering m N-1 outage faults, the overall branch power flow over-limit risk index α P , the corresponding base state probability and fault probability weighted sum of the branch power flow over-limit risk of the system under the base state and various faults is:

[0097] Based on the above indicator α V , α P , α V,j , α P,j , obtain the system comprehensive safety risk index α under fault j s,j , and the system comprehensive security risk index α taking into account all system statess :

[0098]

[0099] Among them, the weight parameter σ1, σ2, are greater than or equal to 0,

[0100] For the above risk measurement model, the calculation of each risk index involves the integral operation of the probability density function, and this nonlinear operation will undoubtedly affect the solvability and computational complexity of the model. Therefore, this embodiment performs linearization processing on it. The specific process is as follows:

[0101] General Rewritten as: in, V i,j The value range of is divided into Z parts, and the bandwidth of each part is h, then: Among them, V i,j,z Represents the midpoint value (or left and right endpoint values) of the zth value interval.

[0102] Suppose the sample set of voltage at node i under fault j is {V i,j} 1×Ns , 1×Ns represents the sample set size, formula f(V i,j ) can actually be obtained through statistical calculation of the sample set. Similarly, V i,j The value range of is divided into Z parts, and the bandwidth of each part is h, then: Among them, V i,j,z Indicates V i,j The midpoint value of the zth interval, n z represents the number of sample points falling in the zth value interval. Then the formula It can be written as: For n that falls in the zth value interval z Samples Respectively represent V i,j The first and nth values ​​distributed in the zth interval z samples, it can be approximately considered that: According to the homogeneous additivity of linear functions, in the function g(V i,j ) has: The half width near the segmentation point is ( Approaching 0 + ) has: Then, the formula and Substitution We can get:

[0103]

[0104] The indicators of the discretized risk model are as follows:

[0105] Node voltage over-limit risk indicator:

[0106]

[0107] Among them, N Σ represents the number of samples, V i,j,n Indicates the nth V i,j sample.

[0108] Branch power flow over-limit risk index:

[0109]

[0110] Among them, N Σ represents the number of samples, P i,j,n Indicates the nth P i,j sample;

[0111]

[0112] System comprehensive security risk indicators:

[0113]

[0114] Therefore, this embodiment transforms the integral operation in the risk indicator into a discretization operation, avoids complex integral operations, and greatly improves the solvability and computational efficiency of the model.

[0115] S2. Construct a dispatch optimization model with the objective function of minimizing the base state power generation cost, expected power generation cost, and expected wind curtailment and load shedding cost, and the constraints of the original scenario set and the AC power flow constraints under faults, the upper and lower limits of the unit output, the unit ramp constraints, and the system comprehensive risk constraints; the original scenario set includes the base state scenario before dispatch and the new energy uncertainty scenario set obtained by sampling the new energy output distribution. The base state scenario is the scenario under the new energy output at the current moment before the system is re-dispatched.

[0116] The objective function of the scheduling optimization model is:

[0117]

[0118] Where T represents the total number of time periods; N g Indicates the number of thermal power generators; N s Indicates the number of scenes in the original scene set; N kIndicates the number of N-1 outages; N d Indicates the number of loads; N w Indicates the number of wind turbines; Δd i,t,k , Δw i,t,k They represent the actual output of the i-th thermal power generator in the k-th N-1 outage base state scenario during period t; C LS , C WS They represent load shedding cost and wind abandonment cost respectively; f i (·) represents the power generation cost function of the i-th thermal power generator with respect to output; Δd i,t,s,k , Δw i,t,s,k They represent the actual output of the i-th thermal power generator, the load shedding amount of the i-th load, and the wind abandonment amount of the i-th wind turbine generator in the t-th period in the s-th scenario under the k-th N-1 outage;

[0119] Among them, the probability of the kth N-1 outage occurring is k=1,...N k ; When k = 0, represents the probability of the ground state scenario; N l is the total number of lines; k 、pro l are the probability of failure of line k and line l respectively.

[0120] The AC power flow constraint is:

[0121]

[0122] Among them, p ij,t,s,k ,q ij,t,s,k are the active power and reactive power flowing from node i to node j in the sth scenario during period t under the kth N-1 outage; v i,t,s,k is the voltage of node i in time period t in the sth scenario under the kth N-1 outage; θ ij,t,s,k is the voltage phase angle difference between node i and node j in period t in the sth scenario under the kth N-1 outage; and is the reactance and susceptance of the line between node i and node j under the kth N-1 outage; They represent the active and reactive outputs of the i-th thermal power generator in the s-th scenario during the k-th N-1 outage, respectively; They represent the active and reactive outputs of the i-th wind turbine in the s-th scenario during the t-th period under the k-th N-1 outage, respectively; They represent the i-th active load and reactive load in the t-th period in the s-th scenario under the k-th N-1 outage; Δp g,i,t,s,k , Δq g,i,t,s,kThey respectively represent the active output adjustment and reactive output adjustment of the i-th thermal power generator in the k-th N-1 outage scenario during period t.

[0123] The comprehensive risk constraint of the system is:

[0124]

[0125] The above formula indicates that the total risk of the system in T time periods is less than the set risk threshold TOL, and the weight parameter prob s Represents the probability of scene s appearing in the original scene set;

[0126] Among them, g v (v i,t,s,k ) is the voltage v of node i under the kth N-1 outage in scenario s during period t i,t,s,k The risk value of crossing the limit is:

[0127]

[0128] In the above formula, v min 、v max They are v i,t,s,k The lower and upper limits of

[0129] g l (p l,t,s,k ) is the power flow p of line l under the kth N-1 outage in scenario s during period t l,t,s,k The risk value of crossing the limit is:

[0130]

[0131] In the above formula, -p l,max 、p l,max They are v i,t,s,k The lower and upper limits of .

[0132] The upper and lower limits of the unit output are:

[0133]

[0134] Among them, P g,max , P g,min , Q g,max , Q g,min They respectively represent the upper and lower limits of active output and reactive output of thermal power generators.

[0135] The unit climbing constraints:

[0136]

[0137] Among them, R gIndicates the unit climbing rate; ΔT indicates the length of time per unit period.

[0138] In order to facilitate the solution, the AC power flow constraint is transformed into a convex constraint by introducing a second-order cone relaxation, including:

[0139] Since the AC power flow constraints contain quadratic terms and nonlinear terms such as cosine and sinine, they are usually difficult to solve. Therefore, it is necessary to introduce the second-order cone relaxation technique and set the intermediate variables: Then it must satisfy:

[0140] (r ij,t,s,k ) 2 +(t ij,t,s,k ) 2 =u i,t,s,k u j,t,s,k , relaxed to: (r ij,t,s,k ) 2 +(t ij,t,s,k ) 2 ≤u i,t,s, k u j,t,s,k ;

[0141] General Substituting the original AC power flow constraint into the base state AC power flow constraint after the second-order cone relaxation is introduced:

[0142]

[0143] (r ij,t,s,k ) 2 +(t ij,t,s,k ) 2 ≤u i,t,s,k u j,t,s,k

[0144] For the comprehensive risk constraint of the system, the piecewise function continuation technology is introduced to transform it into a convex constraint, including:

[0145] Since the constraints formed by piecewise functions cannot be solved directly by the solver, the above piecewise linear functions can be converted into linear constraints by introducing 0-1 variables, as shown in the formula As an example, introduce the continuous variable wv 1,i,t,s,k 、wv 2,i,t,s,k 、wv 3,i,t,s,k 、wv 4,i,t,s,k and 0-1 variable zv 1,i,t,s,k 、zv 2,i,t,s,k 、zv 3,i,t,s,k ,Mode Transformed into the following linear constraints:

[0146] g v (v i,t,s,k )=wv 1,i,t,s,kg v (0)+wv 2,i,t,s,k g v (v min )+wv 3,i,t,s,k g v (v max )+wv 4,i,t,s,k g v (2)

[0147] v i,t,s,k =wv 1,i,t,s,k 0+wv 2,i,t,s,k ·v min +wv 3,i,t,s,k ·v max +wv 4,i,t,s,k 2

[0148] v i,t,s,k =wv 1,i,t,s,k 0+wv 2,i,t,s,k ·v min +wv 3,i,t,s,k ·v max +wv 4,i,t,s,k 2

[0149] wv 1,i,t,s,k +wv 2,i,t,s,k +wv 3,i,t,s,k +wv 4,i,t,s,k =1

[0150] z 1,i,t,s,k +zv 2,i,t,s,k +zv 3,i,t,s,k =1

[0151] 0≤wv 1,i,t,s,k ≤zv 1,i,t,s,k

[0152] 0≤wv 2,i,t,s,k ≤zv 1,i,t,s,k +zv 2,i,t,s,k

[0153] 0≤wv 3,i,t,s,k ≤zv 2,i,t,s,k +zv 3,i,t,s,k

[0154] 0≤wv 4,i,t,s,k ≤zv 3,i,t,s,k

[0155] Similarly, It can be transformed into the following linear constraint:

[0156]

[0157] v i,t,s,k =wv 1,i,t,s,k·0+wv 2,i,t,s,k ·in min +wv 3,i,t,s,k ·in max +wv 4,i,t,s,k ·2

[0158] wv 1,i,t,s,k +wv 2,i,t,s,k +wv 3,i,t,s,k +wv 4,i,t,s,k =1

[0159] from 1,i,t,s,k +v 2,i,t,s,k +v 3,i,t,s,k =1

[0160] 0≤wv 1,i,t,s,k ≤zv 1,i,t,s,k

[0161] 0≤wv 2,i,t,s,k ≤zv 1,i,t,s,k +v 2,i,t,s,k

[0162] 0≤wv 3,i,t,s,k ≤zv 2,i,t,s,k +v 3,i,t,s,k

[0163] 0≤wv 4,i,t,s,k ≤zv 3,i,t,s,k

[0164] p l,t,s,k =wp 1,l,t,s,k ·(-1000)+wp 2,l,t,s,k ·(-p l,max )+wp 3,l,t,s,k ·p l,max +wp 4,l,t,s,k ·1000

[0165] wp 1,l,t,s,k +wp 2,l,t,s,k +wp 3,l,t,s,k +wp 4,l,t,s,k =1

[0166] zp 1,l,t,s,k +zp 2,l,t,s,k +zp 3,l,t,s,k =1

[0167] 0≤wp 1,l,t,s,k ≤zp 1,l,t,s,k

[0168] 0≤wp 2,l,t,s,k ≤zp 1,l,t,s,k +zp 2,l,t,s,k

[0169] 0≤wp3,l,t,s,k ≤zp 2,l,t,s,k +zp 3,l,t,s,k

[0170] 0≤wp 4,l,t,s,k ≤zp 3,l,t,s,k

[0171] Among them, wp 1,i,t,s,k 、wp 2,i,t,s,k 、wp 3,i,t,s,k 、wp 4,i,t,s,k is a continuous variable; zp 1,i,t,s,k 、zp 2,i,t,s,k 、zp 3,i,t,s,k It is a 0-1 variable.

[0172] S3. Use a scene reduction method based on wasserstein distance to reduce the original scene set to obtain an optimal typical scene set and a corresponding probability distribution.

[0173] The original scene set includes a large-scale new energy uncertainty scene set obtained by sampling the new energy output distribution. Due to its large number, the random optimization calculation amount is usually too large to be calculated. Therefore, this embodiment uses a scene reduction method to simplify the original scene set. The probability characteristics of the original scene set are approximated to the greatest extent with the smallest possible scene set. The reduced scene set is called a typical scene.

[0174] The following steps are involved:

[0175] For a given set of original scenes and the number of typical scenes, an optimization problem model is established:

[0176]

[0177] Where N represents the number of original scene sets before reduction; p i Indicates the corresponding scene The probability of k i A 0-1 variable indicating whether the i-th scene is a classic scene; x i,j A 0-1 variable indicating whether the original scene i is assigned to the classic scene j.

[0178] For the above 0-1 planning problem with a given number of classic scenarios, solvers such as Cplex and Gurobi can be directly called to solve it. However, there is still a lack of effective methods for determining the optimal number of classic scenarios. Therefore, this embodiment further proposes the following method for determining the optimal number of typical scenarios, including:

[0179] (1) Set the number of typical scenarios K = 2, give the convergence threshold tol, and set the function value of the first iteration f1 = +∞;

[0180] (2) Solve the optimization problem model and find K classic scenarios based on the values ​​of the 0-1 variables And calculate the probability of the corresponding classic scene [prob1,...,prob K ]:

[0181]

[0182] Among them, x j,i A 0-1 variable indicating whether the original scene j is assigned to the classic scene i;

[0183] (3) If (f K-1 -f K ) / f K-1 ≤tol, the calculation ends, and the number of classic scenes K at this time is the selected optimal value; if (f K-1 -f K ) / f K-1 >tol, then K=K+1, jump to step (2), until the convergence condition is met, and output the calculation result.

[0184] S4, solving the scheduling optimization model based on the optimal typical scenario set and the corresponding probability distribution, obtaining a scheduling plan, and scheduling the system with the scheduling plan. Specifically, it includes the following steps:

[0185] According to the calculation results of solving the optimization problem model, the objective function of the scheduling optimization model and the system comprehensive risk constraint are transformed to obtain a model after scenario reduction, and its objective function is as follows:

[0186]

[0187] Among them, prob s Represents the probability of scene s appearing in the classic scene set after elimination.

[0188] At the same time, the system comprehensive risk constraint is rewritten, and the other constraints remain unchanged. The model after the scenario reduction is solved, including:

[0189] Take λ,λ≥0 to relax the system comprehensive risk constraint to the objective function and obtain the relaxed model:

[0190]

[0191] AC power flow constraints, output upper and lower limit constraints, and ramp constraints under all faults and scenarios

[0192] For a given λ,λ>0, the relaxed model is decomposed into (N k +1)(N s +1) mutually independent sub-problems are solved separately. The details of each sub-problem are as follows:

[0193]

[0194] st fault k, AC power flow constraints under scenario s, output upper and lower limit constraints, climbing constraints

[0195] The following proves the duality of the relaxed problem and the original problem:

[0196] Assume that the feasible domain of the relaxed model is S1, and the feasible domain of the original model is S2. Then when λ>0, we must have:

[0197]

[0198] set up:

[0199]

[0200] When the solution is in S2 (i.e. x∈S2), we have:

[0201]

[0202] When λ>0, we must have:

[0203]

[0204] If and only if the following holds:

[0205]

[0206] The two questions are equivalent:

[0207]

[0208] Therefore, the original problem and maxL(λ) are dual problems. The solution to the original problem can be obtained by solving maxL(λ). Due to the particularity of the L(λ) structure, this embodiment specifically uses a subgradient algorithm combined with parallel computing technology to quickly solve the model. The specific algorithm is as follows:

[0209] 1) Give a larger initial value of λ (λ>0) so that the solution of the first iteration is within the feasible region S2; give a convergence threshold tol.

[0210] 2) Use Gurobi and other solvers combined with parallel computing technology to solve the subproblem L(λ) s,k (s=0,...,K;k=0,...,N k ).

[0211] 3) Solve for the subgradient g:

[0212]

[0213] 4) If g≤0 and gλ≤tol, the calculation ends; if g>0, λ=λ+Δλ, and jump to step 2); if g≤0 and gλ>tol, λ=λ-Δλ, and jump to step 2).

[0214] In order to further verify the feasibility and effectiveness of the method of this embodiment, an experimental example is used to further illustrate it.

[0215] See also Figure 2 , which are 200 wind power original scenarios and predicted scenarios provided in this experimental example, use the open source wind farm data of Eastern Wind Dataset of the American Renewable Energy Laboratory to construct random optimization scenarios, among which the actual wind power output and predicted output data of wind farm No. 3980 are taken as an example to verify the effectiveness of the scenario generation and reduction technology proposed in this embodiment. Starting from the predicted value of wind farm No. 3980 at 0:00 on June 6, 2004, 200 initial scenarios with a length of 12 hours are generated as follows Figure 1 As shown in the figure, the black bold line represents the forecast scenario from 0:00 on June 6, 2004 to 12:00 on June 6, 2004. As can be seen from the figure, the 200 wind power scenarios are distributed around the forecast scenario, and in each scenario, the output difference between adjacent time periods is not too large, which is relatively reasonable. This is because the Markov chain technology retains the time series correlation of the historical wind power data, and because the forecast error is considered to be a conditional probability distribution based on the forecast value (actual value), it avoids the irrationality of output fluctuations in adjacent time periods in a single scenario.

[0216] After generating 200 original scenes, scene reduction was performed using a typical number of scenes K = 2 to 20, see Figure 3 The change of the objective function value, i.e., the Wasserstein distance between the original scene and the typical scene, is plotted when K = 2 to 20. As can be seen from the figure, as K increases, the Wasserstein distance between the original scene and the typical scene decreases, but the decrease gradually decreases.

[0217] See also Figure 4 , Figure 5 , Figure 6 The classic scenes and corresponding scene probabilities after elimination are given when K=5, K=10, and K=15 respectively. As can be seen from the figure, when K=15 and K=10, the distribution of the typical scenes is very similar to the original scenes, and basically most of the features of the original scenes are retained. When K=5, the typical scenes retain the distribution of the original scenes to a great extent. Given a convergence threshold of 0.03, according to the optimal typical scene number determination algorithm of this embodiment, the optimal number of classic scenes K=5 is selected to simplify the calculation amount for random optimization calculation.

[0218] See also Figure 7 This is the topology diagram of the IEEE 39-node system used in this embodiment. The system includes 39 buses, 10 thermal power generators, and in addition, in order to consider the uncertainty of the output of new energy, a wind turbine is added to node 1.

[0219] A comparative example is constructed to illustrate the superiority of this embodiment and the following indicators are used to evaluate the prediction accuracy:

[0220] The calculation example is calculated using the traditional deterministic method that does not consider new energy and failure uncertainty, the stochastic optimization model based on CVaR conditional probability risk, and the method proposed in this embodiment, mainly considering the three evaluation indicators of total system cost, risk value, and calculation time. The calculation method and evaluation indicators are shown in Table 1 below:

[0221] Table 1 Calculation method and evaluation index

[0222]

[0223]

[0224] As can be seen from the above table, the deterministic method does not consider the uncertainty of wind power output and fault uncertainty, so its model scale is small, relatively simple and the solution speed is extremely fast. In addition, the voltage and branch transmission power are strictly limited during the optimization process, so the risk cannot be calculated and the risk value is 0. For the random optimization method based on the CVaR conditional risk index, the CVaR conditional risk index characterizes the expected value of the part of the quantile corresponding to the confidence level under a given confidence level, that is, the expected value of the worst case. Therefore, when the confidence level is larger, the risk value becomes larger. In order to reduce the risk, it is necessary to continuously adjust the thermal power output or even cut the load, which leads to a rapid increase in the total cost of the system. The risk index of the method proposed in the present invention characterizes the expected value of the system at this time. Compared with the method based on the CVaR conditional risk index, the meaning of its characterization is clearer, and it can better indicate the current system safety risk state, and has a lower economic cost. At the same time, combined with the fast parallel calculation method based on Lagrangian decomposition proposed in the present invention, it has obvious advantages in computational efficiency, which fully verifies the effectiveness and superiority of the method of the present invention.

[0225] This embodiment also provides a power system uncertainty optimization dispatching system taking into account security risks, which is characterized by comprising:

[0226] A risk measurement module, which is used to construct a risk measurement model that takes into account both faults and system status to quantify and evaluate system safety; the system status is determined by the voltage of each node in the system and the transmission power value of each line;

[0227] A scheduling optimization model construction module is used to construct a scheduling optimization model, which takes the minimization of the base state power generation cost, the expected power generation cost, and the expected wind curtailment and load shedding cost as the objective function, and takes the AC power flow constraints, the upper and lower limit constraints of the unit output, the unit climbing constraints, and the system comprehensive risk constraints for the original scenario set and faults as the constraint conditions; the original scenario set includes the base state scenario before scheduling and the new energy uncertainty scenario set obtained by sampling the new energy output distribution;

[0228] A scene reduction module, which uses a scene reduction method based on wasserstein distance to reduce the original scene set to obtain an optimal typical scene set and a corresponding probability distribution;

[0229] A solution module solves the scheduling optimization model based on the optimal typical scenario set and the corresponding probability distribution, obtains a scheduling plan, and schedules the system with the scheduling plan.

[0230] Those skilled in the art can understand that the above description is only a preferred embodiment of the present invention and is not intended to limit the present invention. Although the present invention is described in detail with reference to the aforementioned embodiments, those skilled in the art can still modify the technical solutions recorded in the aforementioned embodiments or replace some of the technical features therein by equivalents. Any modification, equivalent replacement, improvement, etc. made within the spirit and principle of the present invention shall be included in the protection scope of the present invention.

Claims

1. A method for optimizing power system uncertainty dispatching taking into account safety risks, characterized in that: include: Constructing a risk measurement model that takes into account both faults and system status to quantify and evaluate system safety; the system status is determined by the voltage of each node in the system and the transmission power value of each line; Construct a dispatch optimization model with the objective function of minimizing the base-state power generation cost, the expected power generation cost, and the expected wind curtailment and load shedding cost, and with the constraints of the original scenario set and the AC power flow constraints under faults, the upper and lower limits of the unit output, the unit ramp constraints, and the system comprehensive risk constraints; the original scenario set includes the base-state scenario before dispatch and the new energy uncertainty scenario set obtained by sampling the new energy output distribution; The original scene set is reduced by using a scene reduction method based on the wasserstein distance to obtain an optimal typical scene set and a corresponding probability distribution; The scheduling optimization model is solved based on the optimal typical scenario set and the corresponding probability distribution to obtain a scheduling plan, and the system is scheduled according to the scheduling plan.

2. The power system uncertainty optimization dispatching method taking into account security risks according to claim 1 is characterized in that: The construction of the risk measurement model includes: Establish the system overall voltage over-limit risk index And the system overall branch power flow over-limit risk index in: m is the number of faults; p s,j is the probability of the jth device N-1 being out of service, i.e., failure j, in the system; α V,j is the system voltage over-limit risk index under fault j, N is the number of nodes in the system, α V,i,j is the voltage over-limit risk of node i under fault j: In the above formula, is the voltage V at node i under fault j i,j The probability density function of V min,i , V max,i are the lower and upper limits of the voltage at node i respectively; α P,j is the system branch power flow over-limit risk index under fault j, N l is the total number of system lines, α P,i,j is the over-limit risk of line i under fault j: In the above formula, represents the transmission power P of line i under fault j i,j The probability density function of max,i represents the maximum transmission power of line i; Based on the α V , α P , α V,j , α P,j Obtain the system comprehensive safety risk index α under fault j s,j , and the system comprehensive security risk index α taking into account all system states s : Among them, the weight parameter σ1, σ2, are both greater than or equal to 0, 3. The method according to claim 2, characterized in that The expression of the risk measurement model is discretized, so that when calculating α s,j , α s In the process, the probability density function The integral operation is converted into a discretization operation.

4. The method according to claim 1, characterized in that The objective function of the scheduling optimization model is: Where T represents the total number of time periods; N g Indicates the number of thermal power generators; N s Indicates the number of scenes in the original scene set; N k Indicates the number of N-1 outages; N d Indicates the number of loads; N w Indicates the number of wind turbines; Δd i,t,k , Δw i,t,k They represent the actual output of the i-th thermal power generator in the k-th N-1 outage base state scenario during period t; C LS , C WS They represent load shedding cost and wind abandonment cost respectively; f i (·) represents the power generation cost function of the i-th thermal power generator with respect to output; Δd i,t,s,k , Δw i,t,s,k They represent the actual output of the i-th thermal power generator, the load shedding amount of the i-th load, and the wind abandonment amount of the i-th wind turbine generator in the t-th period in the s-th scenario under the k-th N-1 outage; Among them, the probability of the kth N-1 outage occurring is k=1,N...;k=0, represents the probability of the ground state scenario; N l is the total number of lines; k 、pro l are the probability of failure of line k and line l respectively.

5. The method according to claim 4, characterized in that The AC power flow constraints of the scheduling optimization model are: Among them, p ij,t,s,k ,q ij,t,s,k are the active power and reactive power flowing from node i to node j in the sth scenario during period t under the kth N-1 outage; v i,t,s,k is the voltage of node i in time period t in the sth scenario under the kth N-1 outage; θ ij,t,s,k is the voltage phase angle difference between node i and node j in period t in the sth scenario under the kth N-1 outage; and is the reactance and susceptance of the line between node i and node j under the kth N-1 outage; They represent the active and reactive outputs of the i-th thermal power generator in the s-th scenario during the k-th N-1 outage, respectively; They represent the active and reactive outputs of the i-th wind turbine in the s-th scenario during the t-th period under the k-th N-1 outage, respectively; They represent the i-th active load and reactive load in the t-th period in the s-th scenario under the k-th N-1 outage; Δp g,i,t,s,k , Δq g,i,t,s,k They respectively represent the active output adjustment and reactive output adjustment of the i-th thermal power generator in the k-th N-1 outage scenario during period t.

6. The method according to claim 5, characterized in that The system comprehensive risk constraint of the scheduling optimization model is: The above formula indicates that the total risk of the system in T time periods is less than the set risk threshold TOL, and the weight parameter prob s Represents the probability of scene s appearing in the original scene set; Among them, g v (v i,t,s,k ) is the voltage v of node i under the kth N-1 outage in scenario s during period t i,t,s,k The risk value of crossing the limit is: In the above formula, v min 、v max They are v i,t,s,k The lower and upper limits of g l (p l,t,s,k ) is the power flow p of line l under the kth N-1 outage in scenario s during period t l,t,s,k The risk value of crossing the limit is: In the above formula, -p l,max 、p l,max They are v i,t,s,k The lower and upper limits of .

7. The method according to claim 4, characterized in that The scene reduction method based on the wasserstein distance is used to reduce the original scene set to obtain the optimal typical scene set and the corresponding probability distribution, including: Build an optimization problem model: Where N represents the number of original scene sets before reduction; p i Indicates the corresponding scene The probability of k i A 0-1 variable indicating whether the i-th scene is a classic scene; x i,j A 0-1 variable indicating whether the original scene i is assigned to the classic scene j; Determine the number of classic scenarios, including: (1) Set the number of typical scenarios K = 2 and give a convergence threshold tol; (2) Solve the optimization problem model and find K classic scenarios based on the values ​​of the 0-1 variables And calculate the probability of the corresponding classic scene [prob1,...,prob K ]: Among them, x j,i A 0-1 variable indicating whether the original scene j is assigned to the classic scene i; (3) If (f K-1 -f K ) / f K-1 ≤tol, the calculation ends, and the number of classic scenes K at this time is the selected optimal value; if (f K-1 -f K ) / f K-1 >tol, then K=K+1, jump to step (2), until the convergence condition is met, and output the calculation result.

8. The method according to claim 7, characterized in that Solving the scheduling optimization model based on the optimal typical scenario set and the corresponding probability distribution includes: According to the calculation results of solving the optimization problem model, the objective function of the scheduling optimization model and the system comprehensive risk constraint are transformed to obtain a model after scenario reduction, and its objective function is as follows: Among them, prob s Represents the probability of scene s appearing in the classic scene set after elimination.

9. The method according to claim 8, characterized in that Solving the scene reduction model includes: Take λ,λ≥0 to relax the system comprehensive risk constraint to the objective function and obtain the relaxed model: AC power flow constraints, output upper and lower limit constraints, and ramp constraints under all faults and scenarios Among them, wv 1,i,t,s,k 、wv 2,i,t,s,k 、wv 3,i,t,s,k 、wv 4,i,t,s,k , and wp 1,i,t,s,k 、wp 2,i,t,s,k 、wp 3,i,t,s,k 、wp 4,i,t,s,k All of them are continuous variables introduced by piecewise function continuation of the system comprehensive risk constraint; For a given λ,λ>0, the relaxed model is decomposed into (N k +1)(N s +1) mutually independent sub-problems are solved separately. The details of each sub-problem are as follows: st fault k, AC power flow constraints, output upper and lower limit constraints, and climbing constraints under scenario s.

10. An uncertainty optimization dispatching system for power systems taking into account safety risks, characterized in that: include: A risk measurement module, which is used to construct a risk measurement model that takes into account both faults and system status to quantify and evaluate system safety; the system status is determined by the voltage of each node in the system and the transmission power value of each line; A scheduling optimization model construction module is used to construct a scheduling optimization model, which takes the minimization of the base state power generation cost, the expected power generation cost, and the expected wind curtailment and load shedding cost as the objective function, and takes the AC power flow constraints, the upper and lower limit constraints of the unit output, the unit climbing constraints, and the system comprehensive risk constraints for the original scenario set and faults as the constraint conditions; the original scenario set includes the base state scenario before scheduling and the new energy uncertainty scenario set obtained by sampling the new energy output distribution; A scene reduction module, which uses a scene reduction method based on wasserstein distance to reduce the original scene set to obtain an optimal typical scene set and a corresponding probability distribution; A solution module solves the scheduling optimization model based on the optimal typical scenario set and the corresponding probability distribution, obtains a scheduling plan, and schedules the system with the scheduling plan.

Images