CN112186766B

CN112186766B - Multistage decomposition method based on Benders decomposition

Info

Publication number: CN112186766B
Application number: CN202011112195.0A
Authority: CN
Inventors: 杨楠; 王璇; 李宏圣; 黎索亚; 叶迪; 黄禹; 董邦天
Original assignee: China Three Gorges University CTGU
Current assignee: China Three Gorges University CTGU
Priority date: 2017-11-30
Filing date: 2017-11-30
Publication date: 2022-07-01
Anticipated expiration: 2037-11-30
Also published as: CN112186766A; CN107947164B; CN112186765A; CN112186765B; CN107947164A

Abstract

A multistage decomposition method based on Benders decomposition is used for decomposing problems of a day-ahead scheduling model of unit combination decision, and comprises a UC decision main problem in a basic scene, a network safety checking sub problem in the basic scene and a robust checking sub problem in a worst scene after decomposition. The method comprises the following steps when the problem decomposition is carried out: step 1, modeling and solving a UC main problem; step 2, modeling and solving a safety subproblem in a basic scene; and 3, modeling and solving the safety sub-problem under the uncertain scene.

Description

Multistage decomposition method based on Benders decomposition

Technical Field

The invention discloses a multistage decomposition method based on Benders decomposition, and relates to the field of power system scheduling.

Background

Wind power and photovoltaic are pollution-free and green renewable energy sources, are widely distributed, have high energy density and are suitable for large-scale development, so that wind power and photovoltaic power generation technologies are highly valued by countries in the world. However, because the output of the power grid has the characteristics of randomness and volatility, the large-scale access to the power grid brings great challenges to the traditional scheduling method, and therefore, the research on the day-ahead scheduling method of the power system under the large-scale access of various new energy resources has important theoretical value and practical significance.

At present, a plurality of experts and scholars research the day-ahead scheduling problem of an electric power system under new energy access from different angles, but generally only a single uncertainty variable is considered, however, an actual electric power system comprises multiple uncertainty factors such as wind power output, photovoltaic output and load prediction error, and the existing day-ahead scheduling only considering single uncertainty is obviously difficult to guarantee the decision effectiveness and influence the economical efficiency of system operation. Therefore, considering the influence of multiple uncertainty factors in the scheduling problem has become a hot point for the research of experts in recent years.

Disclosure of Invention

Aiming at the defects of the existing method, the invention provides an improved robust scheduling method comprehensively considering uncertainty and correlation of load, wind power and photovoltaic.

The technical scheme adopted by the invention is as follows:

a power system day-ahead robust scheduling method considering multiple uncertainties and correlations comprises the following steps:

1) the method comprises the following steps of (1) taking multiple random factors into consideration for day-ahead robust scheduling modeling;

2) solving a worst scenario;

3) solving the model by a Benders decomposition method;

in the step 1), firstly, the day-ahead robust scheduling modeling of multiple random factors is considered; then carrying out day-ahead scheduling modeling in a basic scene; then carrying out day-ahead scheduling modeling in an uncertain scene;

in the step 2), firstly, probability density function modeling of random factors is carried out, then samples are sampled, and then orthogonal transformation matrix derivation is carried out; then, solving a worst scene;

in the step 3), firstly modeling a UC main problem, and then modeling a safety sub-problem in a basic scene; and then modeling a safety sub-problem considering uncertainty.

In the step 1), multiple random factors comprise wind power, photovoltaic and load prediction error factors.

In the step 1), the unit combination decision is divided into a basic scene and a worst scene for modeling respectively.

The basic scenario is based on the uncertainty factor power predicted value, the minimum total operation cost of the system is taken as a target, and various conventional constraint conditions of the system in a deterministic environment are considered.

The worst scenario is based on the maximum fluctuation output of the uncertain power output, and the uncertain constraint condition is considered.

In step 2), the worst scene is sought, and then robustness check is performed by directly using the worst scene.

Firstly, probability density functions of wind power, photovoltaic and load are respectively constructed by utilizing nonparametric kernel density estimation, then, a Latin hypercube sampling is utilized to generate samples, finally, a Cholesky decomposition method is adopted to convert the random samples with correlation into mutually independent random samples, and the worst scene is determined on the basis of the mutually independent random samples.

If the number of load samples based on the historical data is n, constructing a probability density model of the load based on a non-parametric kernel density estimation method as follows:

in the formula, phi (P)_d) Is a probability density function of the load; k (P)_dL) is a kernel function; p_dmIs the m-th sample value in the load sample, and l is the bandwidth.

Selecting a Gaussian function as a kernel function of the load probability density model, and solving the bandwidth l to obtain a probability density function phi (P) of the system load_d)。

Calculating the probability density function phi (P) of the wind power output_w) And a photovoltaic output probability density function phi (P)_v). As shown in formula (2) and formula (3).

In the formula, phi (P)_w) Is a probability density function of the load; k (P)_wL) is a kernel function; p_wmAnd l is the bandwidth, wherein the m is the sample value in the wind power output sample.

In the formula, phi (P)_v) Is a probability density function of the load; k (P)_vL) is a kernel function; p_vmAnd l is the bandwidth, wherein the value is the mth sample value in the photovoltaic output sample.

Latin hypercube sampling is adopted for hierarchical sampling.

Setting the sampling scale as N, Y_m＝F_m(X_m) Denotes the m-th random variable X_mThe specific sampling process of the probability density function of (1) is as follows: will be interval [0,1]Equally dividing the sampling value into N equal parts, selecting the middle value of each subinterval, and obtaining the sampling value through the inverse function of the middle value

And obtaining a sample matrix of all the random variables after sampling of all the random variables is completed.

The method comprises the steps of describing the correlation among loads, wind power and photovoltaic power generation by utilizing a correlation coefficient matrix, and setting a sample matrix obtained by Latin hypercube sampling as W ═ W₁,w₂,…w_l]^TThe matrix of correlation coefficients isC_w：

The elements of the matrix can be obtained by equation (5):

in the formula:

and

are respectively an input variable w_iAnd w_jThe standard deviation of (a); c_ov(w_i,w_j) As an input variable w_iAnd w_jThe covariance of (a). Matrix of correlation coefficients C_wIf it is a positive definite matrix, Cholesky decomposition may be performed on the coefficient matrix:

C_w＝GG^T (6)

in the formula: g is a lower triangular matrix, in which the elements can be found by equation (7):

assuming an orthogonal matrix B, the input random variable W with correlation can be converted into an uncorrelated random variable Y:

Y＝BW (8)

matrix C of correlation coefficients due to uncorrelated random variables Y_YIs an identity matrix I, and thus:

C_Y＝ρ(Y,Y^T)＝ρ(BW,W^TB^T)＝Bρ(W,W^T)B^T＝BC_WB^T＝I (9)

further, it can be obtained from the formula (6):

C_Y＝BC_WB^T＝BGG^TB^T＝(BG)(BG)^T＝I (10)

derived from the above equation:

B＝G^-1 (11)

on the premise that the input uncertainty W having correlation is known, it can be changed to an uncorrelated random variable Y by orthogonal transformation matrix.

And converting the load, wind power and photovoltaic output matrixes with correlation into independent matrixes through the obtained orthogonal transformation matrix, eliminating the correlation among the independent matrixes, and further solving the worst scene through a method of linear superposition of the worst scene.

Solving the model by adopting a Benders decomposition method; decomposing an original problem into a main problem and two sub-problems, wherein the main problem is a UC decision main problem under a basic scene, and the two sub-problems are respectively a network safety check sub-problem under the basic scene and a robust check sub-problem under a worst scene;

the method comprises the following steps: firstly, modeling a UC main problem; then, modeling a safety sub-problem in a basic scene; then, modeling a safety subproblem considering uncertainty;

the method comprises the following specific steps:

1): UC Main problem modeling

The UC major problem includes the target function (1) and the constraint equations (3) through (9), as well as all the Benders cuts generated.

2): safety sub-problem modeling under basic scene

The network security syndrome problem under the basic scene is shown as formula (27), and the network security of the scheduling scheme is ensured by checking the load flow out-of-limit condition of the UC main problem.

In the formula: lambda [ alpha ]_1,l,t,λ_2,l,tIs a dual variable of the network security constraint; v. of_l,tIs the relaxation variable.

Introducing a relaxation variable v into a subproblem_l,tThe role of the method is to temporarily relieve the network security constraint by using a relaxation variable when the constraint condition cannot be met so as to ensure that the subproblem is always solved. V if final optimization is obtained_l,tIf the maximum value is greater than the given safety threshold value, the optimal unit combination scheme obtained by the main problem cannot meet the network safety constraint, and therefore the Benders cut as follows needs to be returned:

3): uncertainty-aware safety sub-problem modeling

The safety check model under the uncertain scene is shown as a formula (29), and the safety check model is checked in the worst scene

And

whether the safety constraint can be met.

S.t.

0≤v_1,lt,v_2t,v_3t(29) In the formula:

and

the actual values of the unit output, the wind power output, the photovoltaic output and the load under the worst scene are obtained. Lambda [ alpha ]_1,it,λ_2,it,μ_1,it,μ_2,it,η_1,it,η_2,itIs a dual variable of a rotation standby constraint, a unit capacity constraint and a climbing constraint.

And if the safety constraint cannot be met, returning to the Benders cut formula (30) to the main problem, wherein the main problem serves as a constraint condition, so that the unit combination and output scheme can be adaptively adjusted on the worst scene to meet the network safety constraint and the power balance.

Through the steps, the power system day-ahead robust scheduling considering multiple uncertainties and correlations is completed.

The invention has the following technical effects:

1) compared with the traditional method, on one hand, the method can effectively account for multiple uncertain factors and the correlation thereof, thereby improving the economical efficiency of system operation while ensuring the robustness of a scheduling decision scheme; on the other hand, the improved robust scheduling model effectively avoids solving the unit combination problem in a massive scene, and has higher solving efficiency.

2) The robust optimization method is improved by introducing a Cholesky decomposition method, so that the worst scene is directly determined without multi-scene load flow calculation, and the applicability and the solving efficiency of the robust optimization model are effectively improved.

Drawings

FIG. 1 is a block diagram of the overall modeling of the present invention;

FIG. 2 is a flow chart of a worst case scenario extraction of the present invention;

FIG. 3 is a block diagram of the general concept of the algorithm of the present invention;

FIG. 4 illustrates wind and photovoltaic active power output curves according to an embodiment of the present invention;

FIG. 5 is a table of worst scenario contrast at time 1;

FIG. 6 is a table of comparison results of scheduling costs from day-ahead;

fig. 7 is a table of comparison of calculated efficiencies.

Detailed Description

step 1: day-ahead scheduling modeling in basic scenarios

The method is based on a robust optimization thought, and a day-ahead scheduling model considering wind power, photovoltaic and load prediction error uncertainty is constructed. The method divides the unit combination decision into a basic scene and a worst scene for respectively modeling, wherein the basic scene takes uncertainty factor power prediction value as a basis, takes the minimum total operation cost of the system as a target, and simultaneously considers various conventional constraint conditions of the system in a deterministic environment, thereby ensuring the economy of the scheduling decision; in the worst scene, the maximum fluctuation output of the uncertain power output is taken as the basis, and the uncertain constraint condition is considered, so that the robustness of the decision scheme in the basic scene in the uncertain environment is ensured. The model block diagram is shown in fig. 1.

Step 1.1: day-ahead scheduling modeling in a basic scenario: including objective functions and constraints

1) Determining an objective function containing fuel cost and unit start-stop cost of a conventional unit:

wherein:

in the formula:

the active output of the unit i at the moment t is represented;

representing the state of the unit i at the time t;

for the cost of starting up and stopping the machine, a_i、b_i、c_iIs a parameter of the generating cost function of the unit.

2) Constraints contained in the model are established, such as:

(1) system power balance constraint, under the condition of not counting network loss:

in the formula:

and the predicted values of the wind power output and the load power at the moment t and the photovoltaic power generation output are obtained.

(2) And (3) output constraint of a conventional unit:

in the formula: p is_i ^minAnd P_i ^maxAnd respectively representing the maximum and minimum active output of the thermal power generating unit i.

(3) Unit ramp rate constraint

Ascending and climbing rate constraint

Descending and climbing rate constraint

In the formula: UR_i,DR_iAnd the climbing power of the conventional unit is limited.

(4) Conventional unit start-stop time constraints

In the formula:

for the on-off time, T, of the unit i at time T_on,iT_off,iAre minimum boot and downtime constraints.

(5) Network security constraints

In the formula:

for maximum current constraints of the line, SF_l,mIs the node power transfer factor. U (m), D (m), W (m) and V (m) are respectively the bus bar sets where the conventional unit, the load, the wind power and the photovoltaic are located.

Step 1.2: day-ahead scheduling modeling in uncertain scenes

The day-ahead unit combination and the unit output plan obtained in the basic scene should ensure the robustness of the system in an uncertain environment, so the day-ahead plan needs to be checked by using system constraint conditions in the uncertain environment. It should be noted that in the actual calculation, the worst scenario may be first sought, and then the worst scenario may be directly utilized for robustness check.

Establishing constraints in the model:

in the formula:

and

the actual values of wind power output, photovoltaic output and load under uncertain conditions.

(2) And (3) output constraint of a conventional unit:

in the formula:

is uncertainThe conventional unit under the circumstances actually goes out.

(3) Unit ramp rate constraint

Uphill climb rate constraint

Descending and climbing rate constraint

(4) Rotational reserve capacity constraint

In the formula:

the device is used for the positive and negative rotation of the conventional unit for standby.

(5) Network security constraints

In the formula:

for maximum current constraint of the line, SF_l,mIs the node power transfer factor. U (m), D (m), W (m) and V (m) are respectively a bus bar set where a conventional unit, a load, wind power and photovoltaic are located.

Step 2: worst scenario solution

As can be seen from step 1, in the actual calculation, the worst scenario may be first sought, and then the worst scenario may be directly used for robustness check. Therefore, the invention provides a method for rapidly calculating the worst scene suitable for multiple correlation randomness factors based on the Cholesky decomposition theory. Firstly, probability density functions of wind power, photovoltaic and load are respectively constructed by utilizing nonparametric kernel density estimation, then, a Latin hypercube sampling is utilized to generate samples, finally, a Cholesky decomposition method is adopted to convert the random samples with correlation into mutually independent random samples, and the worst scene is determined on the basis of the mutually independent random samples.

Step 2.1: probability density function modeling of random factors

If the number of load samples based on historical data is n, constructing a probability density model of the load based on a nonparametric kernel density estimation method as follows:

The invention selects a Gaussian function as a kernel function of the load probability density model and solves the bandwidth l to obtain the probability density function phi (P) of the system load_d)。

The probability density function phi (P) of the wind power output can be obtained by the same method_w) And a photovoltaic output probability density function phi (P)_v). As shown in formula (17) and formula (18).

In the formula, phi (P)_v) Is a probability density function of the load; k (P)_vL) is a kernel function; p_vmIs a photovoltaicThe mth sample value in the output samples, l, is the bandwidth.

Step 2.2: sample sampling

Latin hypercube sampling is a hierarchical sampling method which has a sample memory function and can avoid the sampling of already appeared samples. Let the sampling scale be N, Y_m＝F_m(X_m) Denotes the m-th random variable X_mIs determined. The specific sampling process is as follows: will be interval [0,1]Equally dividing the sampling value into N equal parts, selecting the middle value of each subinterval, and obtaining the sampling value through the inverse function of the middle value

Step 2.3: orthogonal transformation matrix derivation

The method comprises the steps of describing the correlation among loads, wind power and photovoltaic power generation by utilizing a correlation coefficient matrix, and setting a sample matrix obtained by Latin hypercube sampling as W ═ W₁,w₂,…w_l]^TThe matrix of correlation coefficients is C_w：

The elements of the matrix can be obtained by equation (20):

in the formula:

and

are respectively an input variable w_iAnd w_jStandard deviation of (d); c_ov(w_i,w_j) As an input variable w_iAnd w_jThe covariance of (a).

As can be seen by definition, the correlation coefficient matrix C of the present invention_wIf it is a positive definite matrix, Cholesky decomposition may be performed on the coefficient matrix:

C_w＝GG^T (6)

in the formula: g is a lower triangular matrix in which the elements can be found by equation (22).

Y＝BW (8)

C_Y＝ρ(Y,Y^T)＝ρ(BW,W^TB^T)＝Bρ(W,W^T)B^T＝BC_WB^T＝I (9)

further, it can be obtained from the formula (21):

C_Y＝BC_WB^T＝BGG^TB^T＝(BG)(BG)^T＝I (10)

derived from the above equation:

B＝G^-1 (11)

Step 2.4: solving a worst scene;

and (4) converting the load, wind power and photovoltaic output matrixes with correlation into independent matrixes through the orthogonal transformation matrix obtained in the step (2.3), eliminating the correlation among the independent matrixes, and further solving the worst scene through a method of linearly overlapping the worst scene. A worst case scenario retrieval flow diagram is shown in fig. 2.

And step 3: model solution

The invention provides a multistage decomposition algorithm based on Benders decomposition to solve a model.

The traditional Benders decomposition method generally decomposes an original problem into two primary and secondary problems, and decomposes the original problem into a main problem and two sub-problems in consideration of the particularity of the model of the invention, wherein the main problem is a main UC decision problem under a basic scene, and the two sub-problems are respectively as follows: the network safety checking sub-problem under the basic scene and the robust checking sub-problem under the worst scene. The whole frame is shown in figure 1.

Step 3.1: UC Main problem modeling

The main UC problem includes the objective function (1) and constraint equations (3-9) and all the generated Benders cuts.

Step 3.2: safety sub-problem modeling under basic scene

Introducing a relaxation variable v in the subproblem_l,tThe role of the method is to temporarily relieve the network security constraint by using a relaxation variable when the constraint condition cannot be met so as to ensure that the subproblem is always solved. V obtained if final optimization_l,tIf the maximum value is greater than the given safety threshold value, the optimal unit combination scheme obtained by the main problem cannot meet the network safety constraint, and therefore the Benders cut as follows needs to be returned:

step 3.3: uncertainty-aware safety sub-problem modeling

And

whether the safety constraint can be met.

S.t.

0≤v_1,lt,v_2t,v_3t (29)

In the formula:

and

And if the safety constraint cannot be met, returning to the Benders cut formula (30) to the main problem, wherein the Benders cut formula serves as a constraint condition, so that the unit combination and the output scheme can be adaptively adjusted on the worst scene to meet the network safety constraint and the power balance.

And 4, step 4: establishing an operation mode;

in order to comparatively analyze the effectiveness and the correctness of the scheduling model established by the invention, the following two operation modes are established:

mode 1: robust day-ahead scheduling that accounts for multiple uncertainties but does not take correlations into account.

Mode 2: a robust day-ahead scheduling model that takes into account multiple uncertainties and their correlations.

In particular, the invention verifies the correctness of the proposed model with a modified IEEE-118 node system. The system comprises 54 conventional thermal power generating units, 3 wind power plants and 1 photovoltaic power station. The rated power of the wind power plant is respectively 100MW, 200MW and 250MW, and the wind power plant is positioned at nodes No. 5, 9 and 48; the capacity of the photovoltaic power station is 300MW and is positioned at the No. 20 node. The active power output curve of wind power and photovoltaic is shown in fig. 5. The normal rotation standby requirement of a conventional unit in the system is the maximum load of the system8% of the total system load, the negative spin reserve requirement is 2% of the minimum system load. The threshold values of the problem check values of the safety check sub-problems are all 10^-3MWh. And adopting Latin hypercube sampling to sample the load, the wind power and the photovoltaic output 100 times per hour, and totaling 2400 groups of samples. The related calculations are all completed on an Intel core i5-4460 processor 3.20GHz, 8G memory computer, and Matlab and Cplex 12.5 are adopted to program and solve the example.

According to the worst scenario solution method provided by the present invention, the worst scenario comparison result of the uncertainty factors of the mode 1 and the mode 2 is shown in fig. 5.

As can be seen from fig. 5, if the probabilistic correlation between uncertainty factors is not considered, the worst scenario obtained by simple linear superposition is conservative, which will affect the economy of the future scheduling decision, and it can be seen from the precise calculation of the method herein that the worst scenario is actually impossible to occur due to the existence of the correlation between uncertainty factors.

The day-ahead scheduling schemes in mode 1 and mode 2 are calculated separately, and the results are shown in fig. 6.

As can be seen from fig. 6, compared to the conventional day-ahead scheduling method without considering the probability correlation of multiple uncertainty factors, the method proposed in the present invention reduces the day-ahead scheduling cost of the present example by about 38650 $. The method has the advantages that the method fully considers the probability correlation among wind power, photovoltaic and loads, and effectively avoids extreme scenes which are not possible to occur in the robustness checking process, so that the operation cost of the system is reduced while the robustness of the system is ensured by the day-ahead scheduling decision.

In order to compare the advantages of the method of the invention compared with the traditional method, the standard Benders decomposition method and the method of the invention are respectively adopted to simulate the same example, and the calculation efficiency is compared. The comparative results are shown in FIG. 7.

As can be seen from fig. 7: compared with the calculation efficiency, the method provided by the invention is remarkably improved, and the worst scene solving method based on Cholesky decomposition effectively avoids solving problems of unit combination and load flow calculation under a large number of scenes, so that the calculation efficiency is improved by 411.5%.

The present invention has been described in terms of the preferred embodiments, but the above embodiments are not intended to limit the present invention in any way, and all technical solutions obtained by substituting equivalents or equivalent variations fall within the scope of the technical solutions of the present invention.

Claims

1. A multistage decomposition method based on Benders decomposition is characterized in that the method is used for decomposing problems of a day-ahead scheduling model of unit combination decision, and the decomposed problems comprise a UC decision main problem in a basic scene, a network safety checking sub problem in the basic scene and a robust checking sub problem in a worst scene;

the method comprises the following steps when the problem decomposition is carried out:

step 1, modeling and solving a UC main problem;

step 2, modeling and solving a safety subproblem in a basic scene;

step 3, modeling and solving a safety subproblem under the uncertain scene;

in step 2, the sub-problem of security in the basic scene is as formula (27), and the network security of the scheduling scheme is ensured by checking the out-of-limit situation of the power flow of the main problem of UC;

in the formula: lambda [ alpha ]_1,l,t,λ_2,l,tIs a dual variable of the network security constraint, v_l,tIn order to be a function of the relaxation variable,

representing the maximum current constraint, SF, of the line_l,mRepresents the power transfer factor of the node, U (m), D (m), W (m) and V (m) are respectivelyAnd the conventional generator set, the load, the wind power and the photovoltaic are located on the bus collection.

2. A multi-stage decomposition method based on Benders' decomposition according to claim 1,

when the constraint condition can not be met, temporarily relieving the network security constraint by using a relaxation variable so as to ensure that the subproblems are always solved; v obtained if final optimization_l,tIf the maximum value is greater than the given safety threshold value, the optimal unit combination scheme obtained by the main problem cannot meet the network safety constraint, and therefore the Benders cut as follows needs to be returned:

in the formula (I), the compound is shown in the specification,

representing the active power output, SF, of the unit i at time t_l,mFor the node power transfer factor, U (m), W (m) and V (m) are respectively the bus collection where the conventional generator set, the wind power and the photovoltaic are located.

3. The Benders decomposition-based multi-stage decomposition method according to claim 1, wherein in step 3, the safety subproblem model under uncertain scenarios is shown as formula (29), which verifies that in worst scenario

And

whether the safety constraint can be met or not can be judged,

S.t.

0≤v_1,lt,v_2t,v_3t (29)

in the formula:

and

actual values of the output, wind power output, photovoltaic output and load of the unit under the worst scene; lambda_1,it,λ_2,it,μ_1,it,μ_2,it,η_1,it,η_2,itIs a rotational standby constraintThe dual variables of the unit capacity constraint and the climbing constraint, U (m), D (m), W (m) and V (m) are bus sets where a conventional unit, a load, wind power and photovoltaic are respectively located; SF_l,mIn order for the power transfer factor of the node,

for maximum current constraint of the line, P_i ^minAnd P_i ^maxRespectively representing the minimum and maximum active power output of the thermal power generating unit i,

and

for reserve in positive and negative rotation of conventional units, UR_iAnd DR_iAnd the climbing power of the conventional unit is limited.

4. A multistage decomposition method based on Benders decomposition according to claim 3, characterized in that if the safety constraint cannot be satisfied, Benders' cut (30) is returned to the main problem as a constraint condition so that the unit combination and contribution scheme adaptively adjusts to satisfy the network safety constraint and power balance on the worst scenario;

and

for reserve for positive and negative rotation of conventional units, UR_iAnd DR_iFor conventional units climbing power limitation, P_i ^minAnd P_i ^maxRespectively representing the minimum and maximum active power output of the thermal power generating unit i,

representing the state of the unit i at the time t;

and the active output of the unit i at the moment t is shown.