Data-driven wasserstein fuzzy set based active distribution network distribution robust day-ahead scheduling method

Abstract

The application discloses a kind of based on data-driven wasserstein fuzzy set active distribution network distribution robust day-ahead scheduling method, realize the construction compact DRO fuzzy set, and effectively reduce the conservativeness of day-ahead scheduling result;The method constructs the wasserstein condition generated adversarial network model (CWGAN-GP) based on gradient penalty norm, for wind, light output day-ahead scene generation, and proposes the abnormal sample identification method of improved DBSCAN clustering combined with CNN-BiGRU automatic encoder, to improve the credibility of generated scene set;Adopt the compact boundary of data support set determined based on non-parametric kernel density estimation (NKDE) confidence interval, and combined with wasserstein metric to construct DRO fuzzy set.The application compared with the existing DRO day-ahead scheduling method, realizes the deep combination of data-driven method and DRO model, effectively reduces the conservativeness of fuzzy set, and improves the economy of day-ahead scheduling scheme and the adaptability of coping with new energy output uncertainty under the premise of guaranteeing decision robustness.

Description

Technical Field

[0001] This invention relates to the field of active distribution network day-ahead scheduling, specifically a data-driven Wasserstein fuzzy set-based active distribution network day-ahead scheduling method that considers the uncertainty of renewable energy output and the coordinated regulation of "source, grid, load, and storage". Background Technology

[0002] With the large-scale integration of flexible loads and distributed generation, the "high-voltage and high-efficiency" characteristics of the power system are becoming increasingly prominent, and the operation and control methods are becoming more complex. The inherent intermittency and randomness of distributed generation output, represented by wind and solar power, make it difficult to predict accurately. In cases of drastic fluctuations in actual wind and solar output, the planned day-ahead dispatch may fail to guarantee the economic efficiency and stability of system operation, and in severe cases, may even trigger operational accidents such as power flow exceeding limits and load shedding.

[0003] Therefore, considering the uncertainty of renewable energy output, formulating a reasonable and reliable day-ahead dispatch plan for the active distribution network through coordinated regulation of "source, grid, load, and storage" is of significant practical importance. Commonly used uncertainty handling methods primarily include stochastic optimization and robust optimization. Stochastic optimization requires prior assumptions about the probability distribution of uncertain variables and characterizing their uncertainty through simulation of numerous scenarios. On the one hand, the true distribution of uncertain variables is often difficult to obtain effectively in practice; on the other hand, a large number of computational scenarios can severely reduce the efficiency of problem solving. Robust optimization does not require assumptions about the probability distribution of uncertain variables; instead, it controls the fluctuation range of uncertain variables by setting boundary parameters of a closed set.

[0004] However, robust optimization only considers the worst-case scenario, resulting in overly conservative solutions. Partially robust optimization, on the other hand, combines the characteristics of stochastic and robust optimization to find feasible solutions that optimize the expected value of the objective function under the worst-case probability distribution of uncertain variables, thus effectively balancing the economy and robustness of the decision-making results. The construction methods of partially robust fuzzy sets generally include four types: the moment information method, the KL divergence method, the discrete scenario method, and the Wasserstein distance method. Considering that the moment information method cannot effectively guarantee the convergence of the fuzzy set to the true distribution, and that the KL divergence method and the discrete scenario method are only applicable to the discrete distribution of uncertain variables, the Wasserstein distance method not only fully utilizes sample data information but also exhibits good out-of-sample performance for both discrete and continuous distribution uncertain variables.

[0005] Therefore, this invention proposes an active distribution network day-ahead scheduling method based on data-driven Wasserstein fuzzy sets. By utilizing data generation and anomaly identification techniques, a highly reliable sample set of wind and solar day-ahead output is constructed. Based on this, a DRO fuzzy set is constructed by combining NKDE and Wasserstein metrics, thus proposing an active distribution network day-ahead scheduling method based on data-driven Wasserstein fuzzy sets. Summary of the Invention

[0006] This invention provides an active distribution network day-ahead scheduling method based on data-driven Wasserstein fuzzy sets, achieving the construction of a compact DRO fuzzy set and effectively reducing the conservatism of day-ahead scheduling results. The method constructs a Wasserstein Conditional Generative Adversarial Network (CWGAN-GP) model based on gradient penalty norm for generating day-ahead scenarios for wind and solar power output. It also proposes an anomaly sample identification method combining a CNN-BiGRU autoencoder and improved DBSCAN clustering to enhance the credibility of the generated scenario set. Based on this, a confidence interval based on nonparametric kernel density estimation (NKDE) is used to determine the compact boundary of the data support set, and a DRO fuzzy set is constructed using Wasserstein metrics. Furthermore, an active distribution network day-ahead scheduling method based on data-driven Wasserstein fuzzy sets is proposed.

[0007] This invention is achieved using the following technical solution:

[0008] Step 1: Initialize the parameters of the simulation system, including the reference power, reference voltage, branch and node parameters, gas turbine, energy storage device, distributed photovoltaic and wind turbine, demand response parameters and node set.

[0009] Step 2: Train the CWGAN-GP model based on historical predicted and measured data of distributed photovoltaic and wind turbine output.

[0010] Step 3: Based on the day-ahead forecast values ​​of distributed photovoltaic and wind turbine output, generate N sets of day-ahead scenarios for wind and solar power output using the CWGAN-GP model trained in Step 2.

[0011] Step 4: An autoencoder based on a CNN-BiGRU architecture is used to perform feature dimensionality reduction on the N sets of day-ahead wind and solar power output scenes generated in Step 3. Based on the dimensionality reduction results, an improved DBSCAN algorithm is used to identify abnormal samples. Specifically, the improved DBSCAN algorithm utilizes an improved Sparrow Search Algorithm (ISSA) to adaptively optimize the DBSCAN parameters.

[0012] Step 5: Based on the wind and solar power output day-ahead generated scene set after Step 4 and the known day-ahead predicted values, calculate the wind and solar power output prediction error scene set.

[0013] Step 6: The probability distribution of the prediction error scenario set obtained in Step 5 is fitted using NKDE, and this is used as the central empirical distribution of the Wasserstein fuzzy set. Confidence interval optimization is then used to select the boundary of the support set (i.e., the prediction error scenario set obtained in Step 5) to construct the DRO fuzzy set. Based on this, a two-stage active distribution network day-ahead scheduling model is proposed, realizing the organic combination of data-driven methods and DRO.

[0014] Step 7: Based on strong duality theory, the two-stage DRO model constructed in Step 6 is transformed into a mixed integer linear programming problem (MILP) to achieve an efficient solution.

[0015] Beneficial effects

[0016] Compared with existing DRO day-ahead scheduling methods, this invention achieves a deep integration of data-driven methods and DRO models, effectively reducing the conservatism of fuzzy sets, and improving the economy of day-ahead scheduling schemes and their adaptability to uncertainties in renewable energy output while ensuring decision robustness. Attached Figure Description

[0017] Figure 1 This is a flowchart of an active distribution network day-ahead scheduling method based on data-driven Wasserstein fuzzy sets, which is involved in this invention.

[0018] Figure 2 This is a structural diagram of the CGAN model involved in this invention;

[0019] Figure 3 This is a schematic diagram of the CNN-BiGRU autoencoder structure involved in this invention;

[0020] Figure 4 This is a flowchart of the method for identifying anomalies in wind and solar power output samples based on ISSA-DBSCAN, which is involved in this invention. Detailed Implementation

[0021] 1. The construction and training of the CWGAN-GP model involved in step 2 are as follows:

[0022] The CGAN model mainly consists of a generator and a discriminator. By introducing conditional information, it combines supervised learning with GAN to generate samples that satisfy specified features. Its basic structure is as follows: Figure 2As shown. Considering the gradient vanishing and mode collapse problems in the traditional CGAN model during training, this invention introduces Wasserstein distance and gradient penalty norm into the discriminator loss function to improve the model training stability. The resulting loss functions of the generator and discriminator in the CWGAN-GP model are shown in equations (69) and (70):

[0023] Loss G =-E x'～P(x') [D(x'|c)] (69)

[0024]

[0025] In the formula, Loss G and Loss D Let G and D represent the loss functions of the generator and discriminator, respectively; c represents the input conditional information; E x～P(x) (·) and E x′～P(x′) (·) represent the expected values ​​of the true data distribution P(x) and the generated data distribution P(x′), respectively; This represents a random sample value on the line connecting the real data x and the generated data x', satisfying... This represents the introduced gradient penalty norm.

[0026] The training process of the CWGAN-GP model is as follows:

[0027] Step 1: x the predicted wind and solar power output values pre The condition variable c is concatenated with the random noise vector z and input into the generator G to generate the daytime wind and solar power output samples G(z|c).

[0028] Step 2: After concatenating the actual wind and solar power output values ​​x and the generated sample G(z|c) with the condition variable c, input them into the discriminator model D and output the corresponding discrimination results.

[0029] Step 3: Feed the discrimination results from Step 2 back to the generator G and discriminator D models to update the network parameters.

[0030] Step 4: After reaching the maximum number of iterations, the model stops training and extracts the current generator G as the scene generation model. Input the day-ahead predicted value to generate the corresponding day-ahead samples of wind and solar power output.

[0031] 2. The construction of the anomaly identification model based on CNN-BiGRU-AE feature dimensionality reduction and ISSA-DBSCAN clustering in step 4. The specific details are as follows:

[0032] (1) Feature dimensionality reduction method based on CNN-BiGRU autoencoder

[0033] An automatic encoder mainly consists of two parts: an encoder and a decoder. Its basic working principle is shown in equations (71)-(72).

[0034]

[0035]

[0036] In the formula, φ and W represent the encoder and decoder models, respectively; en and W de These represent the encoder and decoder network parameters, respectively; X represents the input data. Z represents the reconstructed output data; Z represents the dimensionality reduction feature. Encoder The high-dimensional input data X is mapped to a low-dimensional vector space to obtain the dimensionality-reduced features Z; the decoder φ reconstructs the input Z into the output based on the dimensionality-reduced features. The loss function of an autoencoder is designed to minimize the difference between the encoder's original input X and the decoder's reconstructed output. The mean square error between them is shown in equation (73).

[0037]

[0038] In the formula, N is the number of training samples; X i and Let be the i-th training input sample and its corresponding reconstructed output.

[0039] Considering the significant variable and temporal coupling characteristics between the generated wind and solar power output samples, this invention employs an autoencoder based on an improved CNN-BiGRU hybrid network design to fully extract the deep-level features between the wind and solar power output sequences. Its structure is as follows: Figure 3 As shown.

[0040] The daytime wind and solar power output samples X = [x1,…x2] generated by CWGAN are used to calculate the daytime wind and solar power output samples. i ,…x N ] T As input, where x i ∈R 2×24The encoder consists of a CNN layer, a pooling layer, and a BiGRU layer. First, the CNN layer extracts the variable coupling features between wind and solar power outputs, reducing the input data dimension from 2×24 to 1×24, thus achieving dimensionality reduction. Second, the pooled data is then input into the BiGRU layer, where a bidirectional GRU structure is used to extract the temporal coupling features between time series data, achieving dimensionality reduction in the time dimension. Finally, a Softmax activation function layer normalizes the output of the BiGRU layer, resulting in the dimensionality-reduced feature result f. i =[f1,f2,…f L ], i∈[1,N],f i ∈R 1×L L represents the length of the dimensionality reduction feature, which is determined by the number of neurons in that layer.

[0041] The decoder employs a symmetrical structure similar to the encoder, first using the dimensionality-reduced features f output by the encoder. i As input, the data is passed through a BiGRU layer and an upsampling layer to achieve initial data reconstruction, and then through a deconvolution layer to obtain the final data reconstruction result. i∈[1,N].

[0042] (2) Anomaly identification method based on ISSA-DBSCAN

[0043] 1) DBSCAN clustering algorithm

[0044] The DBSCAN clustering algorithm is based on the spatial density distribution of sample points. It clusters points in high-density areas into a single cluster and identifies outliers in low-density areas as noise points. The implementation steps are as follows:

[0045] Step 1: Traverse all data samples, calculate the Euclidean distance between sample i and other samples, and compare the Euclidean distance with the neighborhood radius parameter l. Eps The size of the sample is used to obtain the number of sample points N(i) in the neighborhood of sample i, and then the sample points are divided into core sample points X. Core and boundary sample point X Border The division rules are shown in equation (74).

[0046]

[0047] In the formula, δ N This is the minimum threshold for the number of sample points in the neighborhood of the core point.

[0048] Step 2: For any core sample point If sample x j Located at a distance neighborhood radius l Eps If inside, then x is considered j arrive It is based on the direct attainability of density; if in x j ,…,x p If there is a direct density-accessible relationship between them, then x is called x p arrive It is based on density reachability; furthermore, if x q and x p All about If a density reachability relationship exists, then x is called x p and x q It is density-connectable. Sample points that satisfy any of the above relationships can be connected to... They are grouped into the same cluster. This is achieved by traversing all core points X. Core This allows the formation of multiple densely connected clusters with arbitrary shapes.

[0049] Step 3: Sample points that do not belong to any cluster are classified as outlier noise points X. Noise And as an identified anomalous sample.

[0050] 2) Improve SSA

[0051] The Sparrow Optimization Algorithm (SSA) is an intelligent algorithm based on two different behaviors of sparrows: foraging and anti-predation. This algorithm exhibits excellent performance in terms of running speed, convergence accuracy, and global optimization ability. In SSA, sparrows are divided into three roles: discoverers, followers, and alerters. In each iteration, the fitness values ​​of all individuals are calculated and optimally ranked. The sparrows with the highest ranking are selected as discoverers, responsible for providing the population with foraging directions and areas. Followers will follow and monitor the discoverers, and once they detect that the discoverer has found a better foraging location, they will immediately leave their current location to compete for food. Alerters, when they detect approaching danger, will alert the other individuals, and the sparrow population will immediately take anti-predation actions.

[0052] The location update formula for the discoverer is shown in equation (75).

[0053]

[0054] In the formula, Let Iter be the position of the i-th sparrow in the n-th dimension of the search space during the t-th iteration; maxα represents the maximum number of iterations; α and Q are both random numbers, where α∈(0,1] and Q follows a standard normal distribution; L is a 1×D vector of all 1s, where D is the total dimension of the search space; ST is the safety value, ST∈[0.5,1]; R is the warning value, R∈[0,1]. When R<ST, it means that there is no danger for the population to forage, and the discoverer will continue searching along that area; when R≥ST, it means that a sparrow has discovered a predator and sounded the alarm, and the population will immediately take anti-predation actions, with all individuals moving to other areas to continue foraging.

[0055] The position update formula for followers is shown in equation (76).

[0056]

[0057] In the formula, and Let A and B represent the best and worst foraging positions of the population in the (t+1)th and tth iterations, respectively; A is a 1×D matrix, where each element is randomly assigned a value of 1 or -1, and satisfies A + =A T (AA T ) -1 S represents the total population of sparrows. When i > S / 2, it indicates that the i-th follower, due to its low fitness value, cannot obtain food and needs to forage in other areas. When i ≤ S / 2, it indicates that the i-th follower can continue foraging near the optimal location.

[0058] The location update formula for the early warning system is shown in equation (77).

[0059]

[0060] In the formula, β is the step size control parameter, which is a random number following a standard normal distribution. K is a random number taking values ​​in the range [-1, 1], and f i f is the fitness value of the current i-th sparrow individual; g and f w These are the fitness values ​​of the best and worst individuals in the current population, respectively; e is an infinitesimally small real number, used to avoid the denominator being zero.

[0061] Considering that the purpose of combining SSA with DBSCAN is to reasonably determine the hyperparameter values ​​and improve the identification accuracy of wind and solar power output anomalies as much as possible, this section proposes two evaluation metrics, norm_dis and noise_dis, to improve the fitness function of the SSA algorithm, as shown in equations (78) and (79).

[0062]

[0063]

[0064] In the formula, norm_dis represents the intra-class average distance of all normal samples, that is, the average of the sum of the distances from each normal sample point to the other normal sample points; noise_dis represents the inter-class average distance of all abnormal samples to normal samples, that is, the average of the sum of the distances from each abnormal sample point to all normal sample points; N norm N represents the number of normal samples; noise Represents the number of outlier samples; ||x i -x j ||2 represents the sample point x i With x j The L2 norm distance between them.

[0065] The fitness function in the improved SSA is shown in equation (80).

[0066] Fitness=λ·norm_dis-(1-λ)·noise_dis (80)

[0067] In the formula, λ is the weighting factor, λ∈[0,1]. The fitness function is the weighted result of the indices norm_dis and noise_dis; the smaller the value, the better the fitness function is obtained in the current SSA optimization. Eps and δ N The better the DBSCAN algorithm is at identifying abnormal samples, the better it is at identifying abnormal samples.

[0068] 3) Anomaly identification method based on ISSA-DBSCAN

[0069] Based on the dimensionality reduction results of the wind and solar power generation scene set by the CNN-BiGRU autoencoder, the DBSCAN algorithm based on ISSA parameter adaptive optimization is used for clustering. Outlier noise points in the clustering results are identified as anomalous samples, thus achieving anomaly identification for the wind and solar power generation scene set. The implementation process of the ISSA-DBSCAN-based anomalous sample identification method for wind and solar power generation scene set is as follows: Figure 4 As shown. The specific steps are as follows:

[0070] Step 1: Input the dimensionality reduction results of the wind and solar power generation scene set obtained by the CNN-BiGRU autoencoder, and initialize the ISSA parameters, including: DBSCAN parameters l Eps and δ N The search range, the sparrow population size S, and the maximum number of iterations Iter max And randomly generate the initial position of the population in the search space.

[0071] Step 2: For each individual sparrow, optimize the parameters based on its current position (i.e., the optimal value). and The DBSCAN clustering algorithm is executed to obtain the results of the identified normal and abnormal samples.

[0072] Step 3: Based on the DBSCAN identification results, calculate the fitness function values ​​of all individuals according to equation (80), sort them in ascending order, and determine the optimal fitness value of the population in the current iteration. and worst fitness value and the corresponding optimal individual and the worst individual

[0073] Step 4: Update the positions of the discoverers, followers and early warning birds in the sparrow population according to equations (75)-(77), and set the iteration number t = t + 1.

[0074] Step 5: If t≤Iter max Then repeat steps 2 and 3 until the ISSA algorithm termination condition is met. Output the final optimal individual position x. best That is, the optimal combination of DBSCAN parameters and

[0075] Step 6: Based on the optimal parameter combination determined by ISSA, execute the DBSCAN clustering algorithm to obtain the abnormal sample identification results of the wind and solar power generation scene set.

[0076] 3. Step 6 involves the construction of fuzzy sets based on NKDE and Wasserstein metrics, and the construction of a two-stage active distribution network day-ahead scheduling model. The specific details are as follows:

[0077] (1) Construction of fuzzy sets based on NKDE and Wasserstein metric

[0078] The NKDE method does not require prior assumptions about the probability distribution model followed by random variables; instead, it directly estimates the empirical distribution function of random variables using known data samples. Therefore, compared with parameter estimation methods, NKDE can make fuller use of data sample information, and the fitted distribution more closely matches the actual characteristics. Assumptions If there are N samples of wind and solar power output prediction deviations, then the general form of NKDE is shown in equation (81).

[0079] )

[0080] In the formula, These are uncertain parameters, namely, the prediction deviations for wind and solar power output; is the probability density function of uncertain parameters fitted by NKDE; N is the sample size; h is the window width; K h(·) is the kernel function for a window width of h, and the Gaussian kernel function K is usually used. h (x), whose formula is shown in equation (82).

[0081]

[0082] Therefore, the uncertain parameter The probability density function can be expressed as equation (83).

[0083]

[0084] In the formula, It is an uncertain parameter The probability density function.

[0085] The uncertain parameters can be obtained by integrating the probability density function of equation (83). The cumulative probability distribution function is shown in equation (84).

[0086]

[0087] In the formula, F NKDE (x) is an uncertain parameter The cumulative probability distribution function.

[0088] Based on the prediction bias probability distribution obtained from NKDE, this invention proposes to use a confidence level of 1-α. Bilateral quantiles are used to construct the boundary of the uncertainty set, thereby minimizing the conservatism of the uncertainty set while ensuring full utilization of data sample information. The uncertainty set Ξ constructed in this way is shown in equations (85)-(87).

[0089]

[0090]

[0091]

[0092] In the formula, ψ and It is the prediction error bilateral quantiles, ψ It is an uncertain parameter The lower bound value, It is an uncertain parameter The upper bound of Ξ; The uncertain set, i.e. the support set.

[0093] Based on the support set Ξ constructed above, this invention uses the Wasserstein metric to quantify the empirical distribution of wind and solar power output prediction biases. With the true distribution The degree of deviation between them is constructed as an empirical distribution. Centered on, ε N A Wasserstein sphere with radius (β) is used as a sub-Bruker bar fuzzy set Γ. By adjusting the radius ε N The magnitude of (β) can directly affect the performance of the fuzzy set, thereby controlling the conservatism of the optimization results. The expression of the fuzzy set Γ is shown in equation (88).

[0094]

[0095]

[0096] In the formula, Γ is the constructed Wasserstein fuzzy set; M(Ξ) is the uncertain parameter. All possible distributions that exist within its support set Ξ; ε N (β) is the radius of the Wasserstein fuzzy sphere, which is related to the sample size N and the confidence level β, and satisfies Typically, p = 1 is chosen, i.e., L1 norm distance is used; uncertain parameters Follows the true distribution Prediction bias sample Follows empirical distribution for and They follow a joint probability distribution.

[0097] In addition, radius ε N (β) can be obtained by solving the optimization problems shown in equations (90) and (91).

[0098]

[0099]

[0100] In the formula, It is the average value of the data sample; C is a constant, which can be calculated by equation (23); ρ is a non-negative parameter; This is the nth prediction bias sample.

[0101] (2) Construction of a two-stage active distribution network day-ahead scheduling model

[0102] 1) Objective function

[0103] The objective function of the two-stage active distribution network day-ahead scheduling model consists of two stages: pre-scheduling and rescheduling. In the pre-scheduling stage, based on the day-ahead forecasts of wind and solar power output, a day-ahead scheduling plan is formulated for each resource in the system, including switch reconfiguration plans, gas turbine start-up and shutdown and output plans, energy storage charging and discharging plans, demand response plans, power purchase plans, and the amount of wind and solar power curtailment. In the rescheduling stage, considering the disturbances to system operation caused by the uncertainty of actual wind and solar power output, an intraday adjustment plan for flexible resources in the system is formulated, including adjustments to gas turbine output, energy storage charging and discharging, demand response, and the actual amount of wind and solar power curtailment within the scheduling period. The objective function of the two-stage active distribution network day-ahead scheduling model is the minimum sum of the pre-scheduling cost and the expected value of the rescheduling cost under the worst-case distribution p, as shown in equation (92).

[0104]

[0105] In the formula, T is the scheduling period; C Grid C SW C Ess C IS C PG C IDR and C Cur These are, respectively, the main grid power purchase cost during the pre-dispatch phase, the switching operation cost, the energy storage operation and maintenance cost, the gas turbine unit start-up and shutdown cost, the gas turbine unit operating cost, the demand response cost, and the total penalty cost for curtailment of solar and wind power; c price c is the electricity purchase cost coefficient. sw c is the cost coefficient for switching operations. ess c represents the energy storage operation and maintenance cost coefficient. su and c sd These are the start-up and shutdown cost coefficients for the gas turbine; a, b, and c are the operating cost coefficients for each item of the gas turbine; c shift and c rup These are the load transfer and cost reduction coefficients, respectively; c cur,wt and c cur,pv These are the penalty cost coefficients for wind curtailment and solar curtailment, respectively. and As a 0-1 identifier variable representing the change in the state of a branch switch, if This indicates that the branch switch ij changes from an open state to a closed state during time period t; otherwise, it is 0. Similarly, if This indicates that the branch switch ij changes from the closed state to the open state during time period t; otherwise, it is 0. To represent the operating status of a gas turbine unit using 0-1 identifier variables, if This indicates that unit j is in operation; if This indicates that unit j is in a stopped state; and As a 0-1 identifier variable representing changes in the operating state of a gas turbine unit, if This indicates that unit j transitions from a stopped state to a started state during time period t; otherwise, it is 0. Similarly, if This indicates that unit j changed from the start-up state to the shutdown state during time period t; otherwise, it is 0. The active power output of gas turbine j during time period t; The switching power between the root node and the main network during time period t; and These represent the charging and discharging power of energy storage device j during time period t; Let be the active load transfer amount at node j at time t; Let be the active power load reduction at node j at time t; and The abandoned power of distributed photovoltaic power j and wind turbine j during time period t are respectively; L is the set of branches; N pv N wt N Ess N Shift N Rup and N PG These are the sets of nodes belonging to distributed photovoltaic, distributed wind turbines, energy storage devices, transferable loads, load reduction, and gas turbines, respectively.

[0106] Rescheduling phase objective function This includes: energy storage charging and discharging rescheduling costs, gas turbine output rescheduling costs, demand response rescheduling costs, and wind and solar curtailment rescheduling penalty costs, as shown in Equation (93).

[0107]

[0108] In the formula, It is the cost of energy storage charging, discharging, and redistribution. It is the cost of rescheduling gas turbine output. It is the cost of demand response rescheduling. It is the penalty cost of redistributing wind and solar power that is being curtailed; It is the cost coefficient for energy storage operation and maintenance adjustments; It is the cost coefficient for adjusting the output of the gas turbine; and This refers to the adjusted energy storage charging and discharging power. This refers to the adjusted output of the gas turbine. and It refers to the adjusted active power load transfer and reduction amounts; and This refers to the adjusted power of wind and solar power curtailment.

[0109] 2) Constraints related to the pre-scheduling phase

[0110] ① Current constraints

[0111] The original nonlinear power flow constraints based on Distflow are transformed into linear constraints through second-order cone relaxation, as shown in equations (94)-(98).

[0112]

[0113]

[0114]

[0115]

[0116] ||[2P ij,t 2Q ij,t u i,t -i ij,t ] T ||2≤u i,t +i ij,t (98)

[0117] In the formula, δ(j) represents the set of terminal nodes of the branch with j as the first terminal node; λ(j) represents the set of terminal nodes of the branch with j as the first terminal node; P ij,t and Q ij,t Represent the active and reactive power on branch ij during time period t, respectively; r ij and x ij Let I represent the resistance and reactance of branch ij, respectively. ij,t and i ij,t U represents the amplitude and square of the current in branch ij during time period t, respectively; i,t and u i,t P represents the magnitude and square of the voltage at node i during time period t; j,t and Q j,t These represent the active and reactive power injected into node j during time period t, respectively. and Let represent the active load and reactive load of node j in time period t, respectively; This represents the reactive power output of gas turbine j during time period t.

[0118] For branches containing openable switches, this paper uses the Big-M method to relax the constraint equation (96) of the relevant branches and replaces it with constraint equation (99) and equation (100).

[0119]

[0120]

[0121] In the formula, z ij,tLet z be a 0-1 identifier variable representing the state of branch switch ij during time period t. ij,t When z = 1, the branch switch ij is in the closed state; when z = 1, the branch switch ij is in the closed state. ij,t When = 0, the branch switch ij is in the open state; M is any large positive real number.

[0122] ② Network topology constraints

[0123]

[0124]

[0125]

[0126] -M·z ij,t ≤F ij,t ≤M·z ij,t (104)

[0127] In the formula, F ij,t W represents the virtual power of branch ij during time period t. j The virtual power supply output power is unlimited and can be any large real number; other nodes are set with virtual loads of unit 1; n Node and n Sub These represent the total number of system nodes and the number of substation nodes, respectively; N Sub Represents the set of substation nodes, N\N Sub This represents the set of all nodes except for the substation nodes.

[0128] ③Switch action count constraint

[0129]

[0130]

[0131]

[0132]

[0133] In the formula, This indicates the maximum number of times a single switch can be activated within a scheduling period. This indicates the maximum number of all switching actions within the scheduling period.

[0134] ④ Energy storage operation constraints

[0135]

[0136]

[0137]

[0138]

[0139]

[0140] In the formula, Let be a 0-1 flag variable representing the charging and discharging state of Ess at node j in time period t. When, it indicates that Ess is in a charging state. When this occurs, it indicates that Ess is in a discharge state; and These are the maximum charging power and maximum discharging power of Ess, respectively. S represents the remaining battery power of node j at time period t; Ess,min and S Ess,max These are the upper and lower limits of the Ess state of charge, respectively; η Ess,ch and η Ess,dis These are the Ess charge / discharge efficiencies, respectively. Let σ be the Ess capacity at node j; Ess Let be the self-discharge rate of Ess.

[0141] ⑤ Gas turbine constraints

[0142]

[0143]

[0144]

[0145] In the formula, and These are the maximum and minimum limits for the output of the gas turbine, respectively. and These represent the maximum ramp rate and ramp rate of the gas turbine.

[0146] ⑥ Demand response constraints

[0147] • Transferable load

[0148]

[0149]

[0150]

[0151] In the formula, and Let be the active and reactive load transfer amounts at node j at time t, respectively. The power transferred into the node is recorded as a positive value, and the power transferred out of the node is recorded as a negative value. Let be the load power factor of node j; is the maximum allowable transfer coefficient of active power load of the node; constraint (119) indicates that the total load of each node remains unchanged during the scheduling cycle before and after load transfer.

[0152] Interruptible load

[0153]

[0154]

[0155]

[0156] In the formula, and These represent the active and reactive load reductions at node j at time t, respectively. This represents the maximum permissible reduction factor for the active power load at a node. This represents the maximum allowable reduction in total load for node j within the scheduling period.

[0157] ⑦ Power exchange constraints with the upstream power grid

[0158]

[0159]

[0160]

[0161] In the formula, The reactive power exchanged between the root node and the main network during time period t; This represents the maximum limit of active power exchanged between the root node and the main network during time period t. This represents the maximum capacity of the root node substation. Furthermore, since this paper focuses on the operation optimization and scheduling problem within the regional distribution network, it is stipulated that the root node is not allowed to send power back to the upper-level power grid.

[0162] ⑧ Distributed DG output constraints

[0163]

[0164]

[0165] In the formula, and These are the day-ahead forecasts for the output of distributed photovoltaic and wind turbines, respectively.

[0166] ⑨ Operational safety constraints

[0167]

[0168] In the formula, U max and U min These are the upper and lower limits of the node voltage amplitude, respectively; Imax S is the upper limit of the current amplitude in branch ij; ij,max This represents the maximum current carrying capacity of branch ij.

[0169] 3) Constraints related to the rescheduling phase

[0170] To mitigate the impact of wind and solar power output forecasting deviations on system stability during the rescheduling phase, it is necessary to reschedule flexible resources within the system, such as gas turbines, loads that can be reduced or transferred, and energy storage. Assume the random vectors of the wind and solar power output forecasting deviations are as follows: and Then the uncertainty vector of the total prediction bias for wind and solar power can be obtained. As shown in equation (129).

[0171]

[0172] This invention employs an affine strategy as an adjustment strategy for the wind and solar power output prediction deviations during the rescheduling phase. After affine transformation, the expressions for the variables in the rescheduling phase are shown in equations (130)-(136).

[0173]

[0174]

[0175]

[0176]

[0177]

[0178]

[0179]

[0180] In the formula, This is the set of affine adjustment factors for rescheduled variables, used to characterize the response of the system's flexible resources to wind and solar output forecast deviations. To reasonably control the response range, the following is specified... The value range of any element in the set is [-1, 1]. This is a flexible resource rescheduling scheme based on an affine policy. Furthermore, all 0-1 state variables during the rescheduling phase... All of these will remain consistent with the pre-scheduling phase and will not be adjusted further.

[0181] 4. The dual transformation solution method involved in step 7 is as follows:

[0182] After affine transformation, the objective function of the two-stage sub-Bruker model can be converted into a compact form, as shown in equation (137):

[0183]

[0184]

[0185] In the formula, c T and d T These are the coefficient vectors of the objective functions for the pre-scheduling and rescheduling phases, respectively; where, and x, The functional relationship between them is determined by the affine strategy. Equation (138) represents the relevant constraints of the pre-scheduling and rescheduling phases, where A and E are coefficient matrices, b and This is the corresponding parameter vector.

[0186] Since the probability distribution p of the uncertain variables in the support set Ξ is unknown, the model cannot be solved directly. However, according to the strong duality theory, the upper bound problem of the expected objective function under the worst distribution in the rescheduling stage can be transformed into a lower bound problem, as shown in Equation (139).

[0187]

[0188] In the formula, γ is the dual variable; inf{·} represents the infimum function.

[0189] The objective function of the rescheduling stage after dual transformation is merged with the objective function of the outer prescheduling stage to obtain a robust optimization model equivalent to the original DRO model under the support set Ξ, as shown in Equation (140).

[0190]

[0191] However, the objective function (140) also belongs to the min-sup problem at this time, which is difficult to solve directly. Therefore, the auxiliary variable σ is introduced. n (n=1,2,…N), the above min-sup problem is further transformed into equations (141) and (142).

[0192]

[0193]

[0194] Since the objective function (141) after the transformation by introducing auxiliary variables is also linear, and in the constraint (142), In the interval and The above is about variables. convex function, and It is about Since the expression is a linear function, this model belongs to a convex programming problem, and its optimal solution must lie in the uncertain variables. The lower bound ψ, the upper bound and sample set Obtained from [location].

[0195] After the above processing, the two-stage distributed day-ahead scheduling model of the data-driven Wasserstein active distribution network constructed in this invention can be finally transformed into a mixed integer linear programming problem that is easy to solve, as shown in equation (143).

[0196]

[0197] It should be noted that, in this document, relational terms such as "first" and "second" are used only to distinguish one entity or operation from another, and do not necessarily require or imply any such actual relationship or order between these entities or operations. Furthermore, the terms "comprising," "including," or any other variations thereof are intended to cover non-exclusive inclusion, such that a process, method, article, or apparatus that comprises a list of elements includes not only those elements but also other elements not expressly listed, or elements inherent to such process, method, article, or apparatus.

[0198] Although embodiments of the invention have been shown and described, it will be understood by those skilled in the art that various changes, modifications, substitutions and alterations can be made to these embodiments without departing from the principles and spirit of the invention, the scope of which is defined by the appended claims and their equivalents.

Claims

1. A data-driven Wasserstein fuzzy set based active power distribution network distribution robust day-ahead scheduling method, characterized in that: The CWGAN-GP model is used to generate wind and light output day-ahead scenarios, and the abnormal sample identification method based on CNN-BiGRU automatic encoder and improved DBSCAN clustering is used to improve the credibility of the generated scenario set; on this basis, the DRO fuzzy set based on NKDE and Wasserstein metric is constructed; the specific steps are as follows: Step 1: initialize the system parameters of the example, including the reference power, the reference voltage, the branch and node parameters, the gas turbine, the energy storage device, the distributed photovoltaic and wind turbine, the demand response parameters and the node set; Step 2: train the CWGAN-GP model according to the historical prediction data and the historical measured data of the distributed photovoltaic and wind turbine output; Step 3: generate N sets of wind and light output day-ahead scenarios based on the day-ahead prediction value of the distributed photovoltaic and wind turbine output using the CWGAN-GP model trained in step 2; Step 4: use the automatic encoder based on the CNN-BiGRU structure to reduce the dimensionality of the N sets of wind and light output day-ahead scenarios generated in step 3; based on the dimensionality reduction result, use the improved DBSCAN algorithm to realize abnormal sample identification; in the improved DBSCAN algorithm, use the improved sparrow search algorithm (ISSA) to realize adaptive optimization of the DBSCAN parameters; Step 5: calculate the wind and light output prediction error scenario set according to the wind and light output day-ahead scenario set processed in step 4 and the known day-ahead prediction value; Step 6: use NKDE to fit the probability distribution of the prediction error scenario set obtained in step 5 as the center experience distribution of the Wasserstein fuzzy set, and use the confidence interval optimization to select the boundary of the prediction error scenario set obtained in step 5 to construct the DRO fuzzy set; based on this, a two-stage active distribution network distribution robust day-ahead scheduling model is proposed to realize the organic combination of data-driven method and DRO; Step 7: based on the strong duality theory, the two-stage active distribution network distribution robust day-ahead scheduling model constructed in step 6 is converted into a mixed integer linear programming problem (MILP) to realize effective solution.

2. The data-driven Wasserstein fuzzy set based active power distribution network distribution robust day-ahead scheduling method according to claim 1, characterized in that: The construction and training of the CWGAN-GP model in step 2 are as follows: The CGAN model mainly consists of a generator and a discriminator. By introducing conditional information, it combines supervised learning with GAN to generate samples that meet specified characteristics. To address the gradient vanishing and mode collapse problems in traditional CGAN models during training, the Wasserstein distance and gradient penalty norm are introduced into the discriminator loss function to improve model training stability. The loss functions of the generator and discriminator in the CWGAN-GP model are shown in equations (1) and (2): Loss G = -E x'～P(x') [D(x'|c)] (1) where Loss G and Loss D denote the generator G and discriminator D loss functions, respectively; c denotes the input conditional information; E x～P(x) (·) and E x′～P(x′) (·) denote the expected values of the real data distribution P(x) and the generated data distribution P(x'), respectively; denotes a randomly sampled value on the line connecting the real data x and the generated data x', satisfying denotes the introduced gradient penalty norm; The training process of the CWGAN-GP model is as follows: Step 1: wind, light power forecast values x pre As a condition variable c, and with random noise vector z spliced, input into the generator G, to generate wind, light day-ahead power sample G(z|c); Step 2: concatenate the actual wind and light output values x and the generated samples G(z|c) with the conditional variables c respectively, and input them into the discriminator model D, and output the corresponding discrimination results; Step 3: feedback the discrimination results in step 2 to the generator G and discriminator D model to update the network parameters; Step 4: After reaching the maximum number of iterations, the model stops training and extracts the current generator G as the scene generation model; input the day-ahead prediction value to generate the corresponding wind and light output day-ahead sample.

3. The data-driven Wasserstein fuzzy set based active power distribution network distribution robust day-ahead scheduling method of claim 1, wherein: The construction of the abnormal sample identification model based on CNN-BiGRU automatic encoder feature dimension reduction and ISSA-DBSCAN clustering in step 4 is as follows: (1) Feature dimension reduction method based on CNN-BiGRU automatic encoder The automatic encoder is mainly composed of an encoder (Encoder) and a decoder (Decoder), and its basic working principle is shown in equations (3)-(4): wherein, and φ represent the encoder and decoder model, respectively; W en and W de represent the encoder and decoder network parameters, respectively; X represents the input data; represents the reconstructed output data; Z represents the reduced dimension feature; the encoder maps the high-dimensional input data X to a low-dimensional vector space and obtains the reduced dimension feature Z; the decoder φ reconstructs the input Z into the output The loss function of the autoencoder is designed to minimize the mean square error between the original input X of the encoder and the reconstructed output of the decoder, as shown in equation (5): where N is the number of training samples; X i and is the i-th training input sample and its corresponding reconstructed output; Considering the obvious variable coupling and time coupling characteristics between the generated wind and light output samples, an automatic encoder based on a CNN-BiGRU hybrid network is used to fully extract deep features between wind and light output sequences. The wind and light day-ahead output sample X = [x1, … x i ,…x N ] T generated by the CWGAN is taken as input, where x i ∈R 2×24 ; the encoder is sequentially composed of a CNN layer, a pooling layer, and a BiGRU layer; first, the CNN layer is used to extract variable coupling features between wind and light outputs, so that the input data dimension is reduced from 2x24 to 1x24, realizing feature dimension reduction of variables; second, the data processed by the pooling layer is input into the BiGRU layer, and the bidirectional GRU structure is used to extract forward and backward time sequence coupling features of the time sequence, realizing feature dimension reduction of the time dimension; finally, the output of the BiGRU layer is normalized by a Softmax activation function layer, so as to obtain the dimension-reduced feature result f i = [f1, f2, … f L ], i ∈ [1, N], f i ∈R 1×L , and L represents the dimension-reduced feature length, which is determined by the number of neurons in the layer. The decoder adopts a similar symmetrical structure as the encoder, and first reconstructs the dimension-reduced features f i As input, sequentially pass through the BiGRU layer and the upsampling layer to realize preliminary data reconstruction, and then pass through a deconvolution layer to obtain the final data reconstruction result (2) Abnormal sample identification method based on ISSA-DBSCAN 1) DBSCAN clustering algorithm The DBSCAN clustering algorithm is based on the density distribution of sample points in space, and the points in high-density areas are clustered into a cluster, and the outlier points in low-density areas are identified as noise points. The implementation steps are as follows: Step 1: Traverse all data samples, calculate the Euclidean distance between sample i and other samples, compare the Euclidean distance with the size of the neighborhood radius parameter l Eps , obtain the number of sample points N(i) in the neighborhood of sample i, and divide it into core sample points X Core and boundary sample points X Border , the division rule is shown in equation (6): In the formula, δ N is the minimum threshold value of the number of sample points in the core point neighborhood; Step 2: For any core sample point If sample x j Located at a distance neighborhood radius l Eps If inside, then x is considered j arrive It is based on the direct attainability of density; if in x j ,…,x p If there is a direct density-accessible relationship between them, then x is called x p arrive It is based on density reachability; furthermore, if x q and x p All about If a density reachability relationship exists, then x is called x p and x q It is density-connectable; sample points that satisfy any of the above relationships can be connected to... Classify them into the same cluster; by traversing all core points X Core This allows the formation of multiple densely connected clusters with arbitrary shapes; Step 3: The sample points not belonging to any cluster are classified as outlier noise points X Noise and as recognized abnormal samples; 2) Improved SSA The sparrow optimization algorithm is an intelligent algorithm based on the foraging and anti-predation behaviors of sparrow groups. This algorithm has excellent performance in terms of running speed, convergence accuracy, and global optimization ability. In SSA, sparrows are divided into three types: discoverers, followers, and warners. In each iteration, the fitness values of all individuals are calculated and sorted, and the top sparrows are selected as discoverers, which provide the foraging direction and area for the population. Followers follow and monitor the discoverers, and if they detect a better foraging location, they will immediately leave their current location to compete for food. Warners will alert the rest of the individuals when they detect danger, and the sparrow population will immediately take anti-predation behavior. The position update formula of the discoverer is shown in equation (7): wherein, is the position of the ith sparrow in the nth dimension of the search space in the tth iteration; Iter max is the maximum number of iterations; both a and Q are random numbers, where a e (0, 1], and Q follows a standard normal distribution; L is a 1 x D all-one vector, where D is the total dimension of the search space; ST is a safety value, ST e [0.5, 1]; R is a warning value, R e [0, 1]; when R < ST, it means that the population foraging is not in danger, and the discoverer will continue searching in this area; when R > ST, it means that a sparrow has discovered a predator and sent out an alarm, and the population will immediately make an anti-predation behavior, and all individuals will go to other areas to continue foraging; The position update formula of the follower is shown in equation (8): wherein, and respectively represent the best and worst foraging positions of the population in the (t+1)th and tth iterations; A is a 1 x D matrix, wherein each element is randomly assigned as 1 or -1, and satisfies A + = A T (AA T ) -1 ; S is the total number of sparrow population, when i > S / 2, it indicates that the ith follower cannot obtain food due to the low fitness value and needs to go to other areas to forage; when i ≤ S / 2, it indicates that the ith follower can continue to forage near the best position. The position update formula of the warner is shown in equation (9): where β is a step size control parameter, is a random number obeying the standard normal distribution; K is a random number taking values in [-1, 1], f i is the fitness value of the current i-th Sparrow individual; f g and f w are the fitness values of the current best individual and the worst individual in the population, respectively; e is an infinitesimal real number, which acts to avoid a zero denominator; Considering that the purpose of combining SSA with DBSCAN is to reasonably determine the parameter values and improve the identification accuracy of wind and light output abnormal samples as much as possible, two evaluation indicators norm_dis and noise_dis are proposed to improve the fitness function of the SSA algorithm, as shown in equations (10) and (11): In the formula, norm_dis represents the intra-class average distance of all normal samples, that is, the average of the sum of the distances from each normal sample point to the other normal sample points; noise_dis represents the inter-class average distance of all abnormal samples to normal samples, that is, the average of the sum of the distances from each abnormal sample point to all normal sample points; N norm N represents the number of normal samples; noise Represents the number of outlier samples; ||x i -x j ||2 represents the sample point x i With x j L2 norm distance between them; The fitness function of the improved SSA is shown in equation (12): Fitness = λ · norm_dis - (1 - λ) · noise_dis (12) In the formula, λ is a weight factor, λ ∈ [0, 1]; the fitness function is a weighted result of the indexes norm_dis and noise_dis, and the smaller the value is, the better the parameters l Eps and δ N The better the effect of abnormal sample identification of the DBSCAN algorithm is. 3) Abnormal sample identification method based on ISSA-DBSCAN Based on the dimension reduction results of the wind and light output generation scene set by the CNN-BiGRU automatic encoder, the DBSCAN algorithm based on ISSA parameter adaptive optimization is used for clustering, and the outlier noise points in the clustering results are identified as abnormal samples, thereby realizing the abnormal identification of the wind and light output generation scene set. The specific steps are as follows: Step 1: input the dimensionality reduction result of wind, light output generated scene set obtained by the CNN-BiGRU automatic encoder, initialize the ISSA parameters, including: the DBSCAN parameter l Eps and the optimization range of δ N , the sparrow population size S, the maximum iteration number Iter max , and randomly generate the initial position of the population in the search space; Step 2: For each sparrow individual, according to its current position, the DBSCAN clustering algorithm is executed to obtain the recognized normal sample and abnormal sample results; Step 3: Based on the DBSCAN identification result, the fitness function value of all individuals is calculated according to formula (12), and is sorted in ascending order to determine the optimal fitness value of the population in the current iteration and the worst fitness value and the corresponding optimal individual and the worst individual Step 4: According to formula (7) - formula (9), the positions of the discoverer, follower and early warning person in the sparrow population are updated, and the iteration number t is set to t+1; Step 5: If t≤ Iter max then repeat the operation of Step 2 and Step 3 until the ISSA algorithm termination condition is met; output the final optimal individual position x best i.e. the DBSCAN optimal parameter combination and Step 6: Based on the optimal parameter combination determined by ISSA, the DBSCAN clustering algorithm is executed to obtain the abnormal sample recognition result of the wind and light output generation scenario set.

4. The data-driven Wasserstein fuzzy set based active power distribution network distribution robust day-ahead scheduling method of claim 1, wherein: The construction of fuzzy set based on NKDE and Wasserstein metric and the construction of two-stage active distribution network distribution robust day-ahead scheduling model in step 6; The specific content is as follows: (1) Fuzzy set construction based on NKDE and Wasserstein metric The NKDE method does not need to assume the probability distribution model of the random variable in advance, but directly estimates the empirical distribution function of the random variable using the known data sample; Therefore, compared with the parameter estimation method, NKDE can more fully utilize the data sample information, and the fitted distribution is more in line with the actual characteristics; Assume For N samples of wind, light output prediction bias, the general form of NKDE is shown in equation (13): wherein is the uncertainty parameter, i.e. the wind, light output prediction error; is the probability density function of the uncertainty parameter fitted by the NKDE; N is the sample size; h is the window width; K h is the kernel function with window width h, usually a Gaussian kernel K h (x) as shown in equation (14): Thus, the probability density function of the uncertain parameter can be expressed as equation (15): wherein is the probability density function of the uncertain parameter . The cumulative probability distribution function of the uncertain parameter can be obtained by integrating the probability density function of equation (15), as shown in equation (16): where F NKDE (x) is the cumulative probability distribution function of the uncertain parameter x Based on the predicted bias probability distribution obtained by the NKDE, the confidence level of 1-α is used The bilateral quantile constructs the boundary of the uncertainty set, and thus reduces the conservativeness of the uncertainty set as much as possible while ensuring the full use of data sample information. The uncertainty set Ξ constructed in this way is shown in equations (17)-(19). In the formula, Ψ and It is the prediction error bilateral quantiles, Ψ It is an uncertain parameter The lower bound value, It is an uncertain parameter The upper bound of Ξ; The uncertain set, i.e., the support set; Based on the constructed support set Ξ, the Wasserstein metric is used to quantify the deviation between the empirical distribution of wind and power output prediction bias and the true distribution , and a Wasserstein ball with the empirical distribution as the center and ε N (β) as the radius is constructed as the distribution robust fuzzy set Γ; by adjusting the size of the radius Γ N (β), the performance of the fuzzy set can be directly affected, thereby controlling the conservatism of the optimization result, and the expression of the fuzzy set Γ is shown in equations (20) and (21): where Γ is the constructed Wasserstein fuzzy set; M(Ξ) is the uncertainty parameter all distributions possible within its support set Ξ; ε N (β) is the Wasserstein fuzzy ball radius, related to the sample size N and the confidence level β, and satisfies typically take p = 1, i.e., the L1 norm distance; uncertainty parameter subject to the true distribution predicted bias sample subject to the empirical distribution is and subject to the joint probability distribution; Further, the radius ε N (β) can be obtained by solving the optimization problem shown in equations (22) and (23): wherein is the mean value of the data samples; C is a constant, which can be obtained by equation (23); p is a non-negative parameter; is the n-th predicted bias sample; (2) Construction of two-stage active distribution network distribution robust day-ahead scheduling model 1) Objective function The objective function of the two-stage active distribution network distribution robust day-ahead scheduling model is composed of pre-scheduling and rescheduling two stages; The pre-scheduling stage formulates the day-ahead scheduling scheme of each resource of the system according to the day-ahead prediction value of wind and light output, including switch reconstruction plan, gas turbine start-stop and output plan, energy storage charging and discharging plan, demand response plan, power purchase plan and wind and light abandonment amount in the scheduling period; The rescheduling stage formulates the day-to-day adjustment scheme of the flexible resources of the system according to the disturbance of the actual output uncertainty of wind and light to the system operation, including gas turbine output adjustment amount, energy storage charging and discharging adjustment amount, demand response adjustment amount and actual wind and light abandonment amount in the scheduling period; The minimum sum of pre-scheduling cost and expected value of rescheduling cost under the worst distribution p is taken as the objective function of the two-stage distribution robust day-ahead scheduling model, as shown in formula (24): In the formula, T is a scheduling period; C Grid , C SW , C Ess , C IS , C PG , C IDR and C Cur are respectively a pre-scheduling stage main grid electricity purchasing cost, a switch action cost, an energy storage operation and maintenance cost, a gas turbine unit start and stop cost, a gas turbine unit operation cost, a demand response cost and a total penalty cost of abandoned light and wind; c price is a purchasing cost coefficient; c sw is a switch action cost coefficient; c ess is an energy storage operation and maintenance cost coefficient; c su and c sd are respectively a gas turbine start and stop cost coefficient; a, b and c are respectively a gas turbine corresponding operation cost coefficient of each term. c shift and c rup are the load shifting and curtailment cost coefficients, respectively;c cur,wt and c cur,pv are the wind and solar curtailment cost coefficients, respectively; and are 0-1 indicator variables representing the branch switch state change, i.e., if then it means that the branch switch ij changes from open to close at time period t, otherwise it is 0; similarly, if then it means that the branch switch ij changes from close to open at time period t, otherwise it is 0; are 0-1 indicator variables representing the gas turbine unit operation state, i.e., if then it means that the unit j is in operation, otherwise it is 0; if then it means that the unit j is in shutdown, otherwise it is 0; and are 0-1 indicator variables representing the gas turbine unit operation state change, i.e., if then it means that the unit j changes from shutdown to startup at time period t, otherwise it is 0; similarly, if then it means that the unit j changes from startup to shutdown at time period t, otherwise it is 0; is the active power output of the gas turbine j at time period t; P t Grid is the exchange power between the root node and the main grid at time period t; and are the charging and discharging power of the energy storage device j at time period t, respectively; is the active load shifting amount of node j at time t; is the active load curtailment amount of node j at time t; and are the curtailment power of the distributed photovoltaic and wind turbine j at time period t, respectively; L is the branch set; N pv , N wt , N Ess , N Shift , N Rup and N PG are the node sets of distributed photovoltaic, distributed wind turbine, energy storage device, transferable load, curtailable load and gas turbine, respectively; Rescheduling stage objective function The rescheduling stage objective function includes the energy storage charging and discharging rescheduling cost, the gas turbine output rescheduling cost, the demand response rescheduling cost and the wind and solar curtailment rescheduling penalty cost, as shown in equation (25): In the formula, is the energy storage charging and discharging rescheduling cost, is the gas turbine output rescheduling cost, is the demand response rescheduling cost, is the wind and light curtailment rescheduling penalty cost; is the energy storage operation and maintenance adjustment cost coefficient; is the gas turbine output adjustment cost coefficient; and are the adjusted energy storage charging and discharging powers; is the adjusted gas turbine output; and are the adjusted active load transfer and curtailment amounts; and are the adjusted wind and light curtailment powers; 2) Pre-scheduling stage related constraints ①Power flow constraints The original nonlinear power flow constraints based on Distflow are processed by second-order cone relaxation to convert them into linear constraints, as shown in formula (26) - formula (30): ||[2P ij,t 2Q ij,t u i,t -i ij,t ] T ||2≤u i,t +i ij,t (30) where δ(j) denotes the set of branch end nodes with j as the head end node; λ(j) denotes the set of branch head end nodes with j as the end node; P ij,t and Q ij,t denote the real and reactive power on branch ij at time t, respectively; r ij and x ij denote the resistance and reactance of branch ij, respectively; I ij,t and i ij,t denote the magnitude and square of the current on branch ij at time t, respectively; U i,t and u i,t denote the magnitude and square of the voltage at node i at time t; P j,t and Q j,t denote the real and reactive power injection at node j at time t, respectively; and denote the real and reactive load at node j at time t, respectively; denotes the reactive power output of gas turbine j at time t. For branches containing breakable switches, this paper uses Big-M method to relax the constraint formula (28) of the related branches, and replaces it with constraint formula (31) and formula (32): where z ij,t is a 0-1 indicator variable representing the state of branch switch ij at time t, z ij,t = 1 if branch switch ij is closed, and z ij,t = 0 if branch switch ij is open; and M is any large positive real number. ②Network topology constraints - M · z ij,t ≤ F ij,t ≤ M · z ij,t (36) In the formula, F ij,t is the virtual power of the branch ij for the time period t; W j is the unlimited virtual power output, which can be any large real number, and the load values of other nodes are set to be virtual loads of unit 1; n Node and n Sub are the total number of nodes and the number of substation nodes respectively; N Sub represents the set of substation nodes, N\N Sub represents the set of all nodes except the substation nodes; ③Switch action frequency constraints In the formula, represents the upper limit of the number of single switch actions within a scheduling period, represents the upper limit of the total number of all switch actions within a scheduling period; ④Energy storage operation constraints wherein, is a 0-1 flag variable representing the state of charge and discharge of Ess at node j at time t, which represents that Ess is in the state of charging when , and represents that Ess is in the state of discharging when ; and are the maximum charging power and the maximum discharging power of Ess, respectively; is the residual capacity of Ess at node j at time t; S Ess,min and S Ess,max are the upper and lower limit values of the state of charge of Ess, respectively; η Ess,ch and η Ess,dis are the charging and discharging efficiencies of Ess, respectively; is the capacity of Ess at node j; σ Ess is the self-discharge rate of Ess; ⑤Gas turbine constraints wherein and are the maximum and minimum limits of the gas turbine output; and are the maximum and minimum ramp rates of the gas turbine. ⑥Demand response constraints ●Transferable load wherein, and are the active and reactive load transfer of node j at time t, respectively, with positive values for power transfer into the node and negative values for power transfer out of the node; is the power factor of the load of node j; is the maximum allowed active load transfer coefficient of node j; constraint (51) indicates that the total load of each node remains unchanged before and after the load transfer during the dispatching period. ●Interruptible load wherein, and are the active and reactive load curtailment of node j at time t, respectively; is the maximum allowed curtailment factor of the active load of the node; is the maximum allowed value of the total amount of curtailment of the load of node j in the dispatching period. ⑦Exchange power with the upper-level power grid constraint wherein Qmax(t) is the reactive power exchanged by the root node with the main grid for the time period t; Pmax(t) is the maximum limit of active power exchanged by the root node with the main grid for the time period t; Cmax is the maximum capacity of the root node substation; it is provided that the root node is not allowed to send power back to the superior grid. ⑧Distributed DG output constraint wherein, and are the day-ahead forecast values of distributed photovoltaic and wind farm power output, respectively; ⑨Operation safety constraint where U max and U min are the upper and lower limits of the node voltage magnitude; I max is the upper limit of the branch ij current magnitude; S ij,max is the upper limit of the branch ij current magnitude; 3) Rescheduling stage related constraints To mitigate the impact of wind and solar power output prediction error on the stable operation of the system in the rescheduling stage, flexible resources such as gas turbines, reducible and transferable loads, and energy storage need to be rescheduled. Assuming that the wind and solar power output prediction error random vectors are and The total wind and solar power prediction error uncertainty vector can be obtained as as shown in equation (61). Affine strategy is used as the adjustment strategy of the rescheduling stage in response to wind and light output prediction deviation; After affine transformation, the expression of the rescheduling stage variable is as shown in formula (62) - formula (68): In the formula, is a set of rescheduling variable affine adjustment factors, used to represent the degree of response of the system flexible resources to the wind and light output prediction deviation; in order to reasonably control the response range, it is stipulated that Any element in the range [-1, 1] Flexible resource rescheduling scheme according to affine strategy; in addition, all 0-1 state variables in rescheduling phase will not be adjusted, which are consistent with pre-scheduling phase.

Images