Power grid day-ahead scheduling method considering new energy frequency modulation-exit dynamics

Abstract

The invention discloses a power grid day-ahead scheduling method and system considering a new energy frequency modulation-exit dynamic state, and the method is characterized in that the power grid day-ahead scheduling method based on multiple times of frequency drop and the new energy frequency modulation-exit dynamic state comprises the following steps: S1, considering the uncertainty of wind power output, establishing a high-proportion new energy frequency modulation-exit dynamic power grid day-ahead scheduling model; s2, realizing differentiation of the power grid frequency response model, and establishing a dynamic frequency security constraint model considering frequency secondary drop; and S3, dividing the problem of the optimization model into a main problem and a sub-problem by using a Benders decomposition algorithm, and connecting the main problem and the sub-problem by using Benders decomposition to carry out iterative solution. The system dynamic frequency model including the frequency multi-drop process provided by the invention has obvious advantages in guaranteeing the new energy power grid frequency safety problem, and in addition, the new energy frequency modulation-exit dynamic power grid day-ahead scheduling method is considered, so that the dynamic frequency safety problem in the new energy uncertain environment can be effectively solved, and the new energy frequency modulation-exit dynamic power grid day-ahead scheduling method can be effectively optimized. The method is suitable for a high-proportion new energy grid-connected power grid dispatching scene.

Description

Technical Field

[0001] The present invention relates to a power grid day-ahead dispatching method, and in particular to a power grid day-ahead dispatching method considering the frequency modulation-exit dynamics of new energy sources Background Art

[0002] With the advancement of the "dual carbon" goals and the construction of new power systems, the proportion of new energy power generation connected to the grid, represented by wind power, has increased year by year. However, most new energy power generation is connected to the grid through power electronic interfaces and is decoupled from the grid frequency, resulting in weak inertia support capabilities; in addition, in order to ensure its economy, most new energy operates in MPPT mode, resulting in weak frequency support capabilities, posing a great threat to the frequency security of new energy power grids. The realization of reasonable inertia and backup configuration through day-ahead scheduling can effectively alleviate frequency security issues. Therefore, it is urgent to study the day-ahead scheduling method for new energy power grids that takes into account the dynamic security of system frequency.

[0003] However, there are many challenges in the research of day-ahead dispatch of renewable energy power grid considering the dynamic safety of system frequency: first, the frequency dynamics of renewable energy power grid includes multiple frequency drop links related to frequency regulation and withdrawal, and is closely related to the system frequency dynamics and the frequency regulation dynamics of renewable energy, and the model characterization is extremely complex; second, the frequency dynamics of renewable energy power grid is described by a series of differential equations, which is difficult to be compatible with the algebraic optimization model of day-ahead dispatch; finally, the inherent random fluctuation characteristics of renewable energy not only threaten the steady-state power balance of the day-ahead dispatch process, but also cause the uncertainty of transient power disturbances, threatening dynamic frequency safety. Relevant scholars have conducted some research on the day-ahead dispatch of renewable energy power grid considering the dynamic safety of frequency. Most of them use simplified wind turbine models to model the system frequency dynamics, ignoring the secondary drop problem of system frequency caused by the frequency regulation and withdrawal characteristics of the wind turbine itself.

[0004] In terms of uncertainty response, most attention is paid to the impact of random volatility of renewable energy on steady-state power balance, while there is little research on methods to deal with random power fluctuations under frequency safety constraints.

[0005] In summary, existing day-ahead scheduling research is still insufficient in addressing dynamic frequency security issues in uncertain environments.

[0006] This paper proposes a day-ahead dispatching method for power grids that considers the frequency modulation-exit dynamics of renewable energy. This paper fully considers the system frequency dynamics under the whole process of renewable energy frequency modulation-exit, and responds to the impact of renewable energy uncertainty on steady-state power balance and transient power disturbances based on multi-scenario stochastic optimization, which can effectively ensure the dynamic frequency security of renewable energy power grids from the dispatching level. This paper proposes a day-ahead dispatching model for power grids that considers the frequency modulation-exit dynamics of renewable energy, which effectively responds to the dynamic frequency security issues under the uncertain environment of renewable energy; considering the frequency modulation-exit characteristics of renewable energy, a system frequency dynamic model including multiple frequency drop processes is established, and the effective compatibility of the frequency dynamic model and the dispatching model is achieved based on model difference and Benders decomposition algorithm; based on multi-scenario stochastic optimization, the impact of renewable energy uncertainty on steady-state power balance and transient power disturbances is taken into account, effectively responding to the random fluctuation characteristics of renewable energy. Summary of the invention

[0007] Symbols and definitions:

[0008] A. Indices, Sets, and Matrices:

[0009] t,i,ij Time, node, transmission line index k Algorithm Iteration Index

[0010] B. Variables

[0011]

[0012]

[0013] C. Constants

[0014]

[0015] The purpose of the present invention is to provide a power grid day-ahead dispatching method considering the frequency modulation-exit dynamics of new energy. The method proposes a new power grid day-ahead dispatching method, based on the uncertainty of wind power output, constructs a frequency response model under the participation of the whole process of wind power frequency modulation-exit dynamics, and uses the difference method to convert the differential form of the power grid frequency response model into a discretized algebraic form; by considering the secondary drop of frequency, a dynamic frequency safety constraint is established; by using the Benders algorithm, the optimization scheduling problem is decomposed into a main problem and a sub-problem, wherein the main problem is a unit combination solution problem containing conventional constraints, which mainly deals with variables related to the conventional constraints of the unit combination, and the sub-problem is a problem containing dynamic frequency constraints and frequency safety verification, which mainly deals with variables related to the dynamic frequency constraints, and then uses the Benders cut to connect the main problem and the sub-problem for iterative solution until the optimal solution is obtained. The method can effectively deal with the dynamic frequency safety problem in the uncertain environment of new energy, and is suitable for the power grid dispatching scenario with a high proportion of new energy grid connection.

[0016] To achieve the above object, the present invention provides a power grid day-ahead dispatching method considering the frequency modulation-exit dynamics of renewable energy, characterized in that: based on a dispatching method considering the frequency modulation-exit dynamics of renewable energy and multiple frequency drops, the steps are as follows:

[0017] S1. Considering the uncertainty of wind power output, a day-ahead dispatching model of the power grid with high proportion of new energy frequency regulation and withdrawal dynamics is established;

[0018] S2. Differentiate the power grid frequency response model and establish a dynamic frequency safety constraint model that takes into account the secondary frequency drop;

[0019] S3. Use the Benders decomposition algorithm to divide the optimization model into a main problem and sub-problems, and use Benders cut to connect the main problem and sub-problems for iterative solution.

[0020] Optimized, high-proportion renewable energy frequency regulation-exit dynamic grid day-ahead dispatch model

[0021] (1) Wind power output uncertainty modeling

[0022] Currently in the field of wind power research, wind speed is generally regarded as obeying a two-parameter Weibull distribution, and its probability density function is shown in formula (1).

[0023]

[0024] Where: x represents the wind speed; k and c represent the shape parameters corresponding to the Weibull distribution.

[0025] (2) Wind speed distribution description method considering time series fluctuations

[0026] When forecasting wind power output, since the k value and c value can only simply describe the distribution probability of wind speed, it is difficult to obtain a wind speed that conforms to the time series and the wind turbine output at a specific moment. Therefore, it is necessary to consider the time series fluctuation of wind speed.

[0027] This section starts from the law of wind speed change and studies the wind speed fluctuation characteristics. It is found that its fluctuation characteristics are reflected in the wind speed change. The difference between the wind speed values ​​of two adjacent points in the wind speed sequence is the change, as shown in formula (2). Based on the expression of the measured wind speed and Δv, we can get Figure 1 Distribution of wind speed variation.

[0028] Δv={v2-v1,v3-v2,…,v n -v n-1} (2)

[0029] The probability distribution g(x) of the wind speed change is represented by the Δv probability distribution P in formula (3).

[0030] g(x)=P{ζ∈Δv} (3)

[0031] As can be seen from the figure above, the wind speed change is a single-peak distribution across the zero point, but for sequences less than 0, the two-parameter Weibull distribution is difficult to describe, so the overall data is simply processed to make all values ​​in the entire sequence greater than or equal to zero. The wind speed change value is:

[0032] Δv∈[Δv min ,Δv max ] (4)

[0033] Perform Δv processing on the sequence, add the absolute value of the minimum wind speed change to the entire sequence, and get the sequence Δv t As shown in formula (5), at this time, all values ​​in the sequence are greater than or equal to zero.

[0034] Δv t =Δv+|Δv min | (5)

[0035] Δv t ∈[0,|Δv min |+|Δv max |] (6)

[0036] At this time, the sequence Δv t It can satisfy the two-parameter Weibull distribution, given by the sequence Δv t The resulting distribution is g t (x):

[0037]

[0038] Where x represents the change in wind speed

[0039] According to Δv t The numerical relationship with Δv modifies the distribution of wind speed variation expressed by Weibull distribution:

[0040] g(x)=g t (x+|Δv min |) (8)

[0041] The corrected Weibull distribution of wind speed variation g t (x) is:

[0042]

[0043] The wind speed value is regarded as the combined effect of the speed and change at the previous point. Therefore, the wind speed at the previous moment needs to be considered in the calculation process, which is called the basic wind speed. According to the wind speed distribution, the influence of the basic wind speed and the change on the amplitude of the wind speed change is defined as the time series volatility of the wind speed. Let the basic wind speed be n, the Weibull distribution of the wind speed be f(x), and the Weibull distribution of the wind speed change be g(x). The simple joint expression of the basic wind speed and the wind speed change combined with the wind speed distribution is:

[0044] f(x)·g(xn)(10)

[0045] Integrating this expression gives the cumulative probability density function:

[0046]

[0047] Through the description of the probability density properties:

[0048]

[0049] From the above formula, it is not difficult to find that formula (11) does not meet the conditions of the probability density function. Therefore, by introducing an amplification factor m, in order to make it meet the conditions of the probability density function, according to formula (12), the amplification factor m needs to meet formula (13):

[0050]

[0051] When f(x), g(x) and real number n are determined, the amplification factor m can be obtained from formula (14):

[0052]

[0053] In summary, the wind speed distribution function considering time series fluctuations is shown in formula (15), where n representing the basic wind speed is a changing real number, and the amplification factor m also changes with the change of n.

[0054]

[0055] (3) Wind speed scenario reduction

[0056] Based on the change in wind speed and the historical wind speed data, the wind speed time series is considered. The 24-hour wind speed data under 100 scenarios are obtained through Latin hypercube sampling, and the probability of each scenario is the same.

[0057] The 100 scenes sampled are reduced by the Kantorovich distance scene reduction technique, and finally the wind speed data of the required 5 scenes are obtained. The specific process of the Kantorovich distance scene reduction technique is as follows.

[0058] For the generated simulated wind speed data, first calculate the scene wind speed difference L(ξ i ,ξ j );For any sceneξ i , find the scene with the smallest wind speed difference, denoted as ξ s , and calculate the wind speed difference of this scene and ξ i The product of the scene probabilities m(ξ i ,ξ s )p(ξ i ); then find the scene pair with the smallest product (ξ i ,ξ s ):min{m(ξ i ,ξ s )p(ξ i )}, update the scene probability value, p(ξ s )=p(ξ s )+p(ξ i ), and i Cut from the scene; finally update the number of scenes, n1 = n1-n i , n i The number of pairs of scenes with the smallest product is found, and the process is repeated until the set target number of scenes is reached.

[0059] (4) Objective function

[0060] The objective function is the operating cost F of the conventional unit g , start-stop cost F sw and the wind turbine curtailment cost F cut The sum is the smallest, that is

[0061]

[0062] Where: N G is the number of conventional units; C1 is the operating cost coefficient of conventional units; ρ s represents the probability of the sth scenario; represents the operating power of the i-th conventional unit at time t in the s-th scenario; represents the start / stop status of the i-th conventional unit at time t in the s-th scenario; C on，i represents the startup cost of conventional units; N W is the number of wind turbines; C cut is the wind turbine abandonment penalty coefficient; represents the maximum operating power of the jth wind farm under the sth scenario; It represents the operating power of the j-th wind farm at time t in the s-th scenario.

[0063] (5) Constraints

[0064] a) Power balance constraints

[0065] The total power generation and load balance of the operating units:

[0066]

[0067] Where: represents the load power at time t in the sth scenario.

[0068] b) Power constraints of thermal power units

[0069] The output of thermal power units is within a certain range:

[0070] Where: and They represent the upper and lower limits of active and reactive power of the ith thermal power unit in the sth scenario respectively; It represents the output of the i-th thermal power unit at time t in the s-th scenario.

[0071] c) Climbing constraints of thermal power units

[0072] The output of thermal power units is within a certain range

[0073]

[0074] Where: and It represents the maximum output increase and decrease per unit time of the i-th thermal power unit in the s-th scenario; and It represents the increased output when the i-th thermal power unit is started and the reduced output when it is shut down in the s-th scenario.

[0075] d) Wind turbine power constraints

[0076]

[0077] Where: is a Boolean variable, indicating the start and stop status of the j-th wind farm.

[0078] e) Start and stop time constraints for thermal power units

[0079]

[0080] Where: and Indicates the continuous start-up and shutdown hours of the i-th thermal power unit; and Represents the minimum start-up and shutdown hours of the i-th thermal power unit.

[0081] f) Power flow constraints

[0082] Since the power system restoration based on dynamic frequency security mainly considers the active power of the system, the relationship between reactive power and voltage is ignored in the line power flow calculation. In the DC power flow calculation, the line active power flow is expressed as formula (22).

[0083]

[0084] Where: θ s,mk,t is the voltage difference at both ends of line mk at time t; G s,mk and B s,mk are the conductance and susceptance of branch mk respectively.

[0085] In order to simplify the calculation, three simplifying assumptions are made:

[0086] (1) In a normally operating power system, the node voltage is close to the rated voltage, which can be considered as U m =U m =1;

[0087] (2) The voltage difference between the two ends of the line is θ m -θ k is small, so we can assume that sinθ mk =θ mk , cosθ mk =1;

[0088] (3) In a high-voltage transmission network, the line resistance is much smaller than its reactance and can be ignored, i.e., G mk =0.

[0089] In summary, equation (22) can be simplified to equation (23), and the network DC power flow constraint of the system can be expressed by equation (24).

[0090]

[0091] Where: P s,mk,t represents the line transmission power between node m and node k at time t; θ s,m,t represents the phase angle of node m at time t; θ s,k,t represents the phase angle of node k at time t; x s,mk Represents the reactance of branch mk.

[0092]

[0093] In the formula: the direction from node m to node k is defined as the positive direction; Indicates the maximum line capacity of branch mk.

[0094] Preferably, the grid frequency response model differentiation and power system dynamic frequency security constraints:

[0095] (1) Dynamic frequency response model of wind turbine power grid

[0096] The differential form of the rotor motion equation is the physical basis of the dynamic response of the power grid frequency, as shown in equation (25).

[0097]

[0098] Where: Δf is the frequency change; H i , H j and D are the equivalent system inertia constant of conventional units, the equivalent system inertia constant of wind turbines and the damping coefficient respectively; ΔP L is the external power disturbance; ΔP g is the power variation of the conventional unit with primary frequency regulation; ΔP w is the power change of the wind turbine with primary frequency regulation.

[0099] (2) Differential method to construct dynamic frequency safety constraints

[0100] In order to further study the grid dispatch model with dynamic frequency security constraints, the unit combination is considered as a multi-period optimization problem, and the frequency differential equations in each period are converted into discretized algebraic forms through difference equations. The model is simplified to consist of only simple algebraic equations, and the time step involved in the difference equation is set to 0.1s.

[0101] First, the rotor motion equation and the relationship between conventional units, wind turbines and frequency are expressed using differential equations, as shown in formula (26).

[0102]

[0103] After the differential equation is established, the differential equation is processed by the difference method, using finite differences to replace differentials and finite difference quotients to replace derivatives, thereby changing the differential equation into a difference equation (algebraic equation) and converting the model containing differential equations and algebraic equations into a problem containing only algebraic equations, which is convenient for solving. The specific conversion steps are as follows:

[0104] First, the first-order differentials of the frequency change, wind power change and thermal power change are transformed to obtain the differential formula, as shown in formula (27).

[0105]

[0106] Then, the dynamic frequency constraint is derived based on the differential formula of the three variables, and the relationship between the frequency changes at each moment is obtained, as shown in formula (28).

[0107]

[0108] Where: Δf t,n It represents the frequency change in the nth discrete period at time t; It represents the total change of thermal power in the nth discrete period at time t; represents the total change of wind power in the nth discrete period at time t; ΔP t L represents the disturbance power at time t.

[0109] The change of wind power is shown in formula (29):

[0110]

[0111] The change of thermal power is shown in formula (30):

[0112]

[0113] Based on the above formula, set the frequency deviation Δf t,n and the quasi-steady-state frequency deviation Δf ss The maximum value of Strictly ensure the safety of the power grid frequency during the dispatch process, and set the maximum As the upper limit of the dynamic frequency and considering the uncertainty scenario of wind speed output, a complete dynamic frequency constraint is formed as shown in formula (31).

[0114]

[0115] In summary, the dynamic frequency security constraint (31) is incorporated into the day-ahead dispatch model of the renewable energy power system in Chapter 1, and a day-ahead dispatch model of the renewable energy power grid with a high proportion under the dynamic frequency security constraint is obtained.

[0116] (3) Mechanical power model

[0117] The mechanical power model in the per-unit value system is as follows:

[0118]

[0119] Where: P m Capture mechanical power for the wind turbine; C p is the wind energy capture coefficient. To facilitate the linearization of the subsequent model, C p Use the simplified calculation formula above; C p,n is the rated wind energy capture coefficient; v is the actual wind speed; v n is the rated wind speed; β is the pitch angle; λ is the tip speed ratio; λ n is the rated tip speed ratio; ω r is the fan speed; r,n is the rated speed of the fan.

[0120] Formula (32) is Taylor expanded at the initial operating point of the wind turbine, and solved by incremental partial derivative. The mechanical power increment can be obtained by integration as shown in formula (33).

[0121] ΔP m =k ω Δω r +k v Δv (33)

[0122] Where k ω , k v After integration, Δω r , Δv corresponding coefficient.

[0123] (4) Electromagnetic power model

[0124] When analyzing the electromagnetic power output by the fan, the generator part of the fan can be equivalent to the first-order inertia link, and the electromagnetic power output by the generator can be expressed as:

[0125]

[0126] Where: P e P is the electromagnetic power output of the fan; ref is the power setting value of the fan electrical part; α is the first-order inertia constant; P ωr is the MPPT active reference value calculated based on the speed; ΔP is the additional active increment of the fan; K c It is the proportional coefficient of the maximum power tracking curve of the wind turbine.

[0127] The fan electromagnetic power P can be obtained e The increment expression is:

[0128]

[0129] In summary, when considering the complete curve of the secondary frequency drop, the change in wind power is as shown in formula (36).

[0130]

[0131] Replacing the wind power variation in equation (29) with equation (36), we obtain the dynamic frequency constraint considering the secondary frequency drop.

[0132] Preferably, Benders decomposition algorithm

[0133] (1) Model decomposition

[0134] Since wind power is introduced into the power system dispatching problem to participate in the primary frequency regulation and the secondary frequency drop process, and the differential equation of the frequency is differentiated, the dynamic frequency constraint is too large to be solved. Therefore, the optimization model is considered to be decomposed into Figure 5 As shown in the figure, there is a main problem and a sub-problem. The main problem is the unit combination solution problem with conventional constraints, and the sub-problem is the frequency safety verification problem with dynamic frequency constraints. The main problem and the sub-problem are connected through Benders cut.

[0135] The overall problem to be solved is shown in formula (37).

[0136]

[0137] Wherein: x represents the variable related to the conventional constraint of the unit combination, y represents the variable related to the dynamic frequency constraint, the constraint Ax ≥ b represents the conventional constraint of the unit combination, which is composed of equations (17) to (24), and the constraint Ex + Fy ≥ h represents the frequency constraint, that is, equation (31).

[0138] (2) Benders decomposition algorithm

[0139] a) Model optimization process

[0140] The model optimization process mentioned in the present invention consists of the following five steps:

[0141] 1). After initialization, find the initial master problem containing only conventional constraints:

[0142]

[0143] 2). Solve the master problem with only conventional constraints to obtain a solution x* that satisfies only conventional constraints. Then introduce a penalty variable

[0144] Quantify ε into the subproblem and substitute the obtained solution x* into the subproblem:

[0145]

[0146] 3) Dualize the subproblems to obtain:

[0147]

[0148] 4). Look at the value of the penalty variable in the dual subproblem to determine whether the solution obtained by the conventional constraint satisfies the frequency constraint.

[0149] Ⅰ. If satisfied, the optimal solution is obtained and output

[0150] Ⅱ. If not satisfied, proceed to step 5

[0151] 5). By adding constraints to the main problem, the conventional constraint solutions that just do not meet the frequency constraint are eliminated, and the process is repeated until the optimal solution is obtained.

[0152] BRIEF DESCRIPTION OF THE DRAWINGS

[0153] Figure 1 This is the Benders decomposition scheduling model flow chart

[0154] Figure 2 It is the IEEE 24-node system adopted in the embodiment of the present invention;

[0155] Figure 3 is the load curve adopted by the embodiment of the present invention and the wind farm output curve considering dynamic frequency constraints;

[0156] Figure 4 is the number of thermal power units started in 24 periods in consideration of dynamic frequency constraints adopted in the embodiment of the present invention;

[0157] Figure 5 are the maximum frequency deviation and steady-state frequency deviation in each time period obtained by each method in the embodiments of the present invention;

[0158] Figure 6 is the frequency change rate in each time period obtained by each method in the embodiments of the present invention;

[0159] Figure 7 is the maximum deviation of the frequency variation in each time period under the two conditions obtained by the methods of the embodiments of the present invention;

[0160] Figure 8 is the steady-state frequency deviation in each time period under two conditions obtained by the methods of the embodiments of the present invention;

[0161] Fig. 9 is the frequency change rate of each time period in two cases obtained by each method of the embodiments of the present invention;

[0162] Fig.10 is the time domain trajectory of the frequency variation used in the embodiment of the present invention;

[0163] Fig.11 It is the wind power output without considering the frequency constraint and the dynamic frequency constraint adopted in the embodiment of the present invention;

[0164] Fig.12 is the maximum frequency deviation in each time period under the three scenarios obtained by the methods of the embodiments of the present invention;

[0165] Fig.13 is the wind speed scenario taking into account the time series obtained after reduction obtained by each method of the embodiments of the present invention;

[0166] Fig.14 : is a time domain trajectory diagram of frequency variation under three scenarios obtained by various methods of the embodiments of the present invention;

[0167] Fig.15is the number of thermal power units started under the six scenarios obtained by the methods of the embodiments of the present invention;

[0168] Fig.16 is the number of Benders decomposition iterations and convergence error obtained by the method of the embodiment of the present invention; DETAILED DESCRIPTION

[0169] The present invention is further described below in conjunction with the accompanying drawings and specific embodiments.

[0170] Unless otherwise defined, technical or scientific terms used in the present invention shall have the common meanings understood by one having ordinary skills in the field to which the present invention belongs.

[0171] The purpose of the present invention is to provide a power grid day-ahead dispatching method that considers the frequency modulation-exit dynamics of renewable energy. The method proposes a new power grid day-ahead dispatching method. Based on the uncertainty of wind power output, a frequency response model is constructed under the participation of the whole process of wind power frequency modulation-exit dynamics, and the differential method is used to convert the differential form of the power grid frequency response model into a discretized algebraic form; by considering the secondary drop of frequency, a dynamic frequency safety constraint is established; by using the Benders algorithm, the optimization scheduling problem is decomposed into a main problem and a sub-problem, wherein the main problem is a unit combination solution problem containing conventional constraints, which mainly deals with variables related to conventional constraints of unit combination, and the sub-problem is a problem containing dynamic frequency constraints and frequency safety verification, which mainly deals with variables related to dynamic frequency constraints, and then the main problem and the sub-problem are connected by Benders cut for iterative solution until the optimal solution is obtained. The method is decomposed into problems S1, S2, and S3.

[0172] In the embodiment study, the present invention appropriately transforms the 24-node system and constructs a high-proportion wind power grid dispatching model under dynamic frequency constraints, thereby verifying the effectiveness of the day-ahead dispatching model under dynamic frequency security constraints proposed by the present invention.

[0173] A. Test system and parameter settings

[0174] This section makes appropriate modifications to the 24-node system, such as Figure 2 As shown in the figure, it contains 24 buses and 34 branches. It contains 10 thermal power units, located at bus i = {1, 2, 7, 13, 15, 16, 18, 21, 22, 23}, and 6 wind farms, located at bus w = {3, 5, 13, 17, 21, 23}, with capacities of 150MW, 200MW, 200MW, 100MW, 100MW, 50MW, respectively.

[0175] B. Scheduling scheme based on dynamic frequency constraints

[0176] Considering the day-ahead dispatch time scale, the planned dispatch cycle is 24 hours. The base frequency of the example system is 50Hz, the maximum frequency change limit is -0.5Hz, and the power disturbance in each period is set to the wind turbine with the maximum output in the current period offline.

[0177] Depend on Figure 3 and Figure 4 It can be seen that the number of thermal power units in operation is closely related to the load power and wind power output in each period. As the power load continues to increase, it is difficult for wind power output to maintain power balance. At this time, thermal power units begin to start up to maintain the power balance of the power grid.

[0178] According to the obtained optimal unit conditions, the maximum frequency deviation, steady-state frequency deviation and frequency change rate in each period can be obtained, such as Figure 5 and Figure 6 shown.

[0179] Depend on Figure 5 and Figure 6 It can be seen that the maximum frequency deviation, steady-state frequency deviation and frequency change rate in each time period all meet the limit values. Therefore, considering the dynamic frequency security constraint in the day-ahead dispatch model of the power grid can consider economic dispatch on the basis of ensuring the safe operation of the power system.

[0180] C. Effectiveness of scheduling scheme with dynamic frequency safety constraints

[0181] In order to verify the effectiveness of adding dynamic frequency constraints to the power grid dispatch problem proposed in this paper, this section analyzes and compares the dispatch models in two cases:

[0182] Model 1 represents the power grid dispatching model without considering the system frequency security constraints; Model 2 represents the power grid dispatching model considering the system dynamic frequency security constraints.

[0183] It can be seen from the figure that the number of thermal power units started in each period of model 1 without considering frequency constraints is less than that of model 2 with considering dynamic frequency constraints. The maximum frequency change in each period under the two different models is as follows: Figure 7 shown.

[0184] Table 1 Maximum frequency deviation in 24 hours

[0185]

[0186] It can be seen from the figure that the maximum frequency change in each period of model 1 is greater than that of model 2. Taking the 1h period and 11h period as examples, it can be seen from Table 1 that the maximum frequency change is reduced by 0.42999Hz and 0.30802Hz respectively, and the lowest frequency point of model 2 meets the limit value, indicating that considering the dynamic frequency safety constraint can ensure the safe operation of the system. The steady-state frequency deviation of each period under the two sets of models is as follows Figure 8 shown.

[0187] Table 2 Steady-state frequency deviation in 24 periods

[0188]

[0189] It can be seen from the figure that the steady-state frequency deviation of each period of Model 1 without considering frequency constraints is greater than that of Model 2 considering dynamic frequency constraints. Taking the 1h period with the smallest load power and the 11h period with the largest load power as examples, the steady-state frequency deviations under the two different models are compared through Table 2. It is found that compared with Model 1, the steady-state frequency deviation of Model 2 in the 1h period and 11h period are reduced by 0.19591Hz and 0.16303Hz respectively, and the frequency stability of the system is improved. The frequency change rate of each period under the two groups of models is shown in Fig. 9 shown.

[0190] Table 3 Frequency change rate in 24 periods

[0191]

[0192] As can be seen from the figure, the frequency change rate of each period of Model 1 without considering frequency constraints is greater than that of Model 2 considering dynamic frequency constraints. Taking the 1h period with the smallest load power and the 11h period with the largest load power as examples, the frequency change rate under the two groups of different models is compared through Table 3. It is found that compared with Model 1, the steady-state frequency deviation of Model 2 in the 1h period and 11h period is reduced by 0.21615Hz / s and 0.2579Hz / s respectively, and the frequency stability of the system is improved.

[0193] D. Effectiveness of the scheduling scheme considering the secondary frequency drop

[0194] In a transmission network with a high proportion of wind turbine penetration, the proportion of synchronous units is greatly reduced, and the system inertia level and frequency regulation capability are also greatly weakened. In order to verify the necessity and effectiveness of adding frequency secondary drop to dynamic frequency constraints, this section analyzes frequency constraints in three different cases:

[0195] Scenario 1: Only conventional constraints on the day-ahead dispatch of the power system; Scenario 2: Dynamic frequency constraints are added on the basis of Scenario 1; Scenario 3: Secondary drop frequency constraints are added on the basis of Scenario 2.

[0196] Taking the 1-hour period as an example, Table 4 shows the dynamic frequency indicators under three different scenarios. The specific frequency changes within 30 seconds are as follows: Fig.10 shown.

[0197] Table 4 Dynamic frequency indexes for 1h period under three scenarios

[0198]

[0199] It can be seen from Table 4 that the lowest point of the primary frequency modulation in scenario 1 is -0.807Hz, exceeding the limit of 0.307Hz; the lowest point of the secondary frequency modulation is -0.994Hz, exceeding the limit of 0.494Hz; the steady-state frequency is -0.317Hz, exceeding the limit of 0.117Hz; the lowest point of the primary frequency modulation in scenario 2 is -0.438Hz, within the limit; the lowest point of the secondary frequency modulation is -0.591Hz, exceeding the limit of 0.091Hz; the steady-state frequency is -0.167Hz, within the limit; the lowest point of the primary frequency modulation in scenario 3 is -0.412Hz, within the limit; the lowest point of the secondary frequency modulation is -0.287Hz, within the limit; the steady-state frequency is -0.155Hz, within the limit.

[0200] In scenario 1, the dynamic frequency constraint was not considered, and the frequency change rate was -0.5303Hz / s, which was greater than the set frequency change rate limit. After considering the dynamic frequency, the frequency change rates of scenario 2 and scenario 3 were -0.258Hz / s and -0.2426Hz / s, respectively, both meeting the set frequency change rate limit. Compared with not considering the dynamic frequency constraint, the frequency change rate was reduced by about 50%, greatly improving the frequency stability of the system.

[0201] Depend on Fig.11 It can be seen that the wind power output under dynamic frequency constraints is significantly lower than that without considering frequency constraints. Fig.11 ,It is found that the maximum frequency deviation is significantly smaller when the wind power output is low, ,which indicates that when the wind power generation is high, the system inertia is reduced, and at the same time the ,secondary drop of the frequency is affected, resulting in a larger frequency deviation.

[0202] contrast Fig.12 From the maximum frequency deviation under different scenarios, it can be seen that in scenario 1, all periods exceed the maximum frequency limit, while in scenario 2, some cases are within the maximum frequency limit, and scenario 3 fully meets the maximum frequency limit. This shows that in a power grid with a high proportion of wind power penetration, if only dynamic frequency constraints are considered, the secondary frequency drop problem may also be caused by the withdrawal of wind power from frequency regulation, resulting in frequency exceeding the limit; therefore, considering the secondary frequency drop in the dynamic frequency constraint can effectively ensure the frequency security of the power system. The rationality and feasibility of the dispatching model are verified, and the effectiveness of considering dynamic frequency security constraints in power grid dispatching is fully demonstrated.

[0203] Table 5 Comparison of power generation costs under three scenarios

[0204]

[0205] As shown in Table 5, in order to fully meet the frequency constraints of the entire frequency regulation process, the total power generation cost of scenario 3 increased by 3.5% compared with scenario 1. After considering the frequency constraints, the power generation cost increased slightly. However, compared with not considering the frequency constraints, the dynamic frequency indicators of scenario 3 are all within the safe range, ensuring the safe operation of the power system.

[0206] E. Effectiveness of dispatching scheme considering wind power uncertainty

[0207] According to the generation and reduction of wind power uncertainty scenarios in Section 1, five typical wind power output scenario sets are obtained as follows: Fig.13 As shown, the probability of each scenario is listed in Table 6.

[0208] Table 6 Probability of wind power scenarios

[0209]

[0210] For the above wind power output scenarios, under random optimization conditions, the frequency changes of five scenarios are analyzed to explore the impact of the randomness of wind power output on frequency. The discussion is divided into the following six scenarios:

[0211] Scenario 1 to Scenario 5 represent Fig.13 The wind power output of S1-S5 is shown in Figure 6, and scenario 6 represents the wind power output of the five scenarios through probability superposition. In order to consider the impact of wind turbine output on grid frequency, the power disturbance in each period is set to the wind turbine with the maximum output in the current period being offline.

[0212] Table 7 Dynamic frequency constraint indicators for 1h period under six scenarios

[0213]

[0214] Taking the 1h period as an example, it can be seen from Table 7 that the dynamic frequency constraint indicators of scenarios 1 to 6 are all within the limit range. At the same time, it is not difficult to find that the frequency changes are uneven under different wind power generation scenarios, especially in scenarios 2 and 4, although they do not exceed the maximum frequency change limit, they are close to the limit; when considering the comprehensive scenario, the dynamic frequency constraint indicator is not as good as the frequency indicator in some scenarios, but the frequency change is relatively average, and there will be no situation of approaching or exceeding the limit, so that the frequency safety and stability are better reflected.

[0215] The wind power output of three scenarios is selected, and the specific frequency changes within 30 seconds under three different scenarios are given, such as Fig.14 shown.

[0216] Depend on Fig.14It can be seen that all three scenarios can meet the safety requirements of dynamic frequency constraints. Compared with scenario 4, the overall frequency change in scenario 1 is smaller and the frequency safety is higher.

[0217] Depend on Fig.15 It can be seen that due to the large number of thermal power units in operation in scenarios 1, 3, and 5, the inertia is relatively sufficient, the system dynamic frequency response is improved, and the frequency fluctuation amplitude is also suppressed. Compared with scenarios 1 and 4, the frequency change in scenario 6 is more even, which better reflects the frequency safety and stability.

[0218] In addition, since the power disturbance is set to the maximum output wind turbine in the current period, Fig.15 It can be seen that the output of wind turbines in scenario 1 is relatively small, while the output of wind turbines in scenario 4 is large, resulting in large frequency changes. It can be concluded that the uncertainty of wind turbine output poses a certain threat to the frequency security of the power grid. By using stochastic optimization to consider the impact of wind power uncertainty on the power grid frequency based on dynamic frequency constraints, the security and stability of the power grid frequency can be better guaranteed when a high proportion of wind power is connected.

[0219] F. Effectiveness of Benders Algorithm

[0220] In order to verify the effectiveness of the Benders algorithm, the objective function ε of the subproblem in Section 3 is taken as the convergence error. When it approaches 0, it means that the unit commitment optimization result at this time satisfies the frequency constraint, which proves the effectiveness of the Benders algorithm. Fig.15 As shown, when the threshold ε = 0.01, the Benders decomposition converges after iterations, which takes 133 seconds.

[0221] Fig.16 The convergence error diagram of different iteration generations of the model solved by the Benders algorithm shows that in the initial iteration, the convergence error is as high as 0.4. As the iteration continues, the convergence error gradually decreases, and finally the preset convergence accuracy is reached after 14 iterations.

[0222] There is a significant difference in the computational performance between the conventional solution method and the solution method using Benders decomposition for the day-ahead scheduling model calculation considering dynamic frequency security constraints. The specific situation is shown in Table 8.

[0223] Table 8 Comparison of solution rates of the two methods

[0224]

[0225] Table 8 compares the solution speed when the step size is 0.01s. It can be seen from the table that for the day-ahead dispatch model with dynamic frequency security constraints, the solution time using the Benders algorithm is 65.9% shorter than the direct solution, and the solution rate is greatly improved. In summary, it can be shown that the Benders algorithm can effectively solve the day-ahead dispatch model of a high proportion of new energy power grid considering dynamic frequency security constraints proposed in this chapter.

Claims

1. A method and system for day-ahead dispatching of a power grid considering the frequency modulation and withdrawal dynamics of new energy sources, characterized in that: The day-ahead dispatching method of the power grid based on multiple frequency drops and the frequency regulation-exit dynamics of renewable energy is as follows: S1. Considering the uncertainty of wind power output, a day-ahead dispatching model of the power grid with high proportion of new energy frequency regulation and withdrawal dynamics is established; S2. Differentiate the power grid frequency response model and establish a dynamic frequency safety constraint model that takes into account the secondary frequency drop; S3. Use the Benders decomposition algorithm to divide the optimization model into a main problem and sub-problems, and use Benders cut to connect the main problem and sub-problems for iterative solution.

2. The method and system for day-ahead dispatching of a power grid considering the frequency modulation and withdrawal dynamics of new energy sources according to claim 1, characterized in that: Considering the uncertainty of wind power output, a method for day-ahead dispatching of power grid with dynamic frequency regulation and withdrawal of renewable energy is established: (1) Wind power output uncertainty modeling Currently in the field of wind power research, wind speed is generally regarded as obeying a two-parameter Weibull distribution, and its probability density function is shown in formula (1). Where: x represents the wind speed; k and c represent the shape parameters corresponding to the Weibull distribution. (2) Wind speed distribution description method considering time series fluctuations When forecasting wind power output, since the k value and c value can only simply describe the distribution probability of wind speed, it is difficult to obtain a wind speed that conforms to the time series and the wind turbine output at a specific moment. Therefore, it is necessary to consider the time series fluctuation of wind speed. This section studies the wind speed fluctuation characteristics based on the wind speed variation law. It is found that its fluctuation characteristics are reflected in the wind speed variation. The difference between the wind speed values ​​of two adjacent points in the wind speed sequence is the variation, as shown in formula (2). Based on the measured wind speed and the expression of Δv, the distribution of wind speed variation in Figure 1 can be obtained. Δv={v2-v1,v3-v2,…,v n -v n-1 } (2) The probability distribution g(x) of the wind speed change is represented by the Δv probability distribution P in formula (3). g(x)=P{ζ∈Δv} (3) As can be seen from the figure above, the wind speed change is a single-peak distribution across the zero point, but for sequences less than 0, the two-parameter Weibull distribution is difficult to describe, so the overall data is simply processed to make all values ​​in the entire sequence greater than or equal to zero. The wind speed change value is: Δv∈[Δv min ,Δv max ](4) Perform Δv processing on the sequence, add the absolute value of the minimum wind speed change to the entire sequence, and get the sequence Δv t As shown in formula (5), at this time, all values ​​in the sequence are greater than or equal to zero. Δv t =Δv+|Δv min | (5) Δv t ∈[0,|Δv min |+|Δv max |](6) At this time, the sequence Δv t It can satisfy the two-parameter Weibull distribution, given by the sequence Δv t The resulting distribution is g t (x): Where x represents the change in wind speed According to Δv t The numerical relationship with Δv modifies the distribution of wind speed variation expressed by Weibull distribution: g(x)=g t (x+|Δv min |) (8) The corrected Weibull distribution of wind speed variation g t (x) is: The wind speed value is regarded as the combined effect of the speed and change at the previous point. Therefore, the wind speed at the previous moment needs to be considered in the calculation process, which is called the basic wind speed. According to the wind speed distribution, the influence of the basic wind speed and the change on the amplitude of the wind speed change is defined as the time series volatility of the wind speed. Let the basic wind speed be n, the Weibull distribution of the wind speed be f(x), and the Weibull distribution of the wind speed change be g(x). The simple joint expression of the basic wind speed and the wind speed change combined with the wind speed distribution is: f(x)·g(xn) (10) Integrating this expression gives the cumulative probability density function: Through the description of the probability density properties: From the above formula, it is not difficult to find that formula (11) does not meet the conditions of the probability density function. Therefore, by introducing an amplification factor m, in order to make it meet the conditions of the probability density function, according to formula (12), the amplification factor m needs to meet formula (13): When f(x), g(x) and real number n are determined, the amplification factor m can be obtained from formula (14): In summary, the wind speed distribution function considering time series fluctuations is shown in formula (15), where n representing the basic wind speed is a changing real number, and the amplification factor m also changes with the change of n. (3) Wind speed scenario reduction Based on the change in wind speed and the historical wind speed data, the wind speed time series is considered. The 24-hour wind speed data under 100 scenarios are obtained through Latin hypercube sampling, and the probability of each scenario is the same. The 100 scenes sampled are reduced by the Kantorovich distance scene reduction technique, and finally the wind speed data of the required 5 scenes are obtained. The specific process of the Kantorovich distance scene reduction technique is as follows. For the generated simulated wind speed data, first calculate the scene wind speed difference L(ξ i ,ξ j );For any sceneξ i , find the scene with the smallest wind speed difference, denoted as ξ s , and calculate the wind speed difference of this scene and ξ i The product of the scene probabilities m(ξ i ,ξ s )p(ξ i ); then find the scene pair with the smallest product (ξ i ,ξ s ):min{m(ξ i ,ξ s )p(ξ i )}, update the scene probability value, p(ξ s )=p(ξ s )+p(ξ i ), and i Cut from the scene; finally update the number of scenes, n1 = n1-n i , n i The number of pairs of scenes with the smallest product is found, and the process is repeated until the set target number of scenes is reached. (4) Objective function The objective function is the operating cost F of the conventional unit g , start-stop cost F sw and the wind turbine abandonment cost F cut The sum is the smallest, that is Where: N G is the number of conventional units; C1 is the operating cost coefficient of conventional units; ρ s represents the probability of the sth scenario; represents the operating power of the i-th conventional unit at time t in the s-th scenario; represents the start / stop status of the i-th conventional unit at time t in the s-th scenario; C on，i represents the startup cost of conventional units; N W is the number of wind turbines; C cut is the wind turbine abandonment penalty coefficient; represents the maximum operating power of the jth wind farm under the sth scenario; It represents the operating power of the j-th wind farm at time t in the s-th scenario. (5) Constraints a) Power balance constraints The total power generation and load balance of the operating units: Where: represents the load power at time t in the sth scenario. b) Power constraints of thermal power units The output of thermal power units is within a certain range: Where: and They represent the upper and lower limits of active and reactive power of the ith thermal power unit in the sth scenario respectively; It represents the output of the i-th thermal power unit at time t in the s-th scenario. c) Climbing constraints of thermal power units The output of thermal power units is within a certain range Where: and It represents the maximum output increase and decrease per unit time of the i-th thermal power unit in the s-th scenario; and It represents the increased output when the i-th thermal power unit is started and the reduced output when it is shut down in the s-th scenario. d) Wind turbine power constraints Where: is a Boolean variable, indicating the start and stop status of the j-th wind farm. e) Start and stop time constraints for thermal power units Where: and Indicates the continuous start-up and shutdown hours of the i-th thermal power unit; and Represents the minimum start-up and shutdown hours of the i-th thermal power unit. f) Power flow constraints Since the power system restoration based on dynamic frequency security mainly considers the active power of the system, the relationship between reactive power and voltage is ignored in the line power flow calculation. In the DC power flow calculation, the line active power flow is expressed as formula (22). Where: θ s,mk,t is the voltage difference at both ends of line mk at time t; G s,mk and B s,mk are the conductance and susceptance of branch mk respectively. In order to simplify the calculation, three simplifying assumptions are made: (1) In a normally operating power system, the node voltage is close to the rated voltage, which can be considered as U m =U m =1; (2) The voltage difference between the two ends of the line is θ m -θ k is small, so we can assume that sinθ mk =θ mk , cosθ mk =1; (3) In a high-voltage transmission network, the line resistance is much smaller than its reactance and can be ignored, i.e., G mk =0. In summary, equation (22) can be simplified to equation (23), and the network DC power flow constraint of the system can be expressed by equation (24). Where: P s,mk,t represents the line transmission power between node m and node k at time t; θ s,m,t represents the phase angle of node m at time t; θ s,k,t represents the phase angle of node k at time t; x s,mk Represents the reactance of branch mk. In the formula: the direction from node m to node k is defined as the positive direction; Indicates the maximum line capacity of branch mk.

3. The method and system for day-ahead dispatching of a power grid considering the frequency modulation and withdrawal dynamics of new energy sources according to claim 1, characterized in that: Differentiation of the grid frequency response model: (1) Dynamic frequency response model of wind turbine power grid The differential form of the rotor motion equation is the physical basis of the dynamic response of the power grid frequency, as shown in equation (25). Where: Δf is the frequency change; H i , H j and D are the equivalent system inertia constant of conventional units, the equivalent system inertia constant of wind turbines and the damping coefficient respectively; ΔP L is the external power disturbance; ΔP g is the power variation of the conventional unit with primary frequency regulation; ΔP w is the power change of the wind turbine with primary frequency regulation. (2) Differential method to construct dynamic frequency safety constraints In order to further study the grid dispatch model with dynamic frequency security constraints, the unit combination is considered as a multi-period optimization problem, and the frequency differential equations in each period are converted into discretized algebraic forms through difference equations. The model is simplified to consist of only simple algebraic equations, and the time step involved in the difference equation is set to 0.1s. First, the rotor motion equation and the relationship between conventional units, wind turbines and frequency are expressed using differential equations, as shown in formula (26). After the differential equation is established, the differential equation is processed by the difference method, using finite differences to replace differentials and finite difference quotients to replace derivatives, thereby changing the differential equation into a difference equation (algebraic equation) and converting the model containing differential equations and algebraic equations into a problem containing only algebraic equations, which is convenient for solving. The specific conversion steps are as follows: First, the first-order differentials of the frequency change, wind power change and thermal power change are transformed to obtain the differential formula, as shown in formula (27). Then, the dynamic frequency constraint is derived based on the differential formula of the three variables, and the relationship between the frequency changes at each moment is obtained, as shown in formula (28). Where: Δf t,n It represents the frequency change in the nth discrete period at time t; It represents the total change of thermal power in the nth discrete period at time t; represents the total change of wind power in the nth discrete period at time t; ΔP t L represents the disturbance power at time t. The change of wind power is shown in formula (29): The change of thermal power is shown in formula (30): Based on the above formula, set the frequency deviation Δf t,n and the quasi-steady-state frequency deviation Δf ss The maximum value of Strictly ensure the safety of the power grid frequency during the dispatch process, and set the maximum As the upper limit of the dynamic frequency and considering the uncertainty scenario of wind speed output, a complete dynamic frequency constraint is formed as shown in formula (31). In summary, the dynamic frequency security constraint (31) is incorporated into the day-ahead dispatch model of the renewable energy power system in Chapter 1, and a day-ahead dispatch model of the renewable energy power grid with a high proportion under the dynamic frequency security constraint is obtained.

4. The method and system for day-ahead dispatching of a power grid considering the frequency modulation and withdrawal dynamics of new energy sources according to claim 1, characterized in that: Dynamic frequency safety constraint model considering secondary drop (1) Mechanical power model The mechanical power model in the per-unit value system is as follows: Where: P m Capture mechanical power for the wind turbine; C p is the wind energy capture coefficient. To facilitate the linearization of the subsequent model, C p Use the simplified calculation formula above; C p,n is the rated wind energy capture coefficient; v is the actual wind speed; v n is the rated wind speed; β is the pitch angle; λ is the tip speed ratio; λ n is the rated tip speed ratio; ω r is the fan speed; r,n is the rated speed of the fan. Formula (32) is Taylor expanded at the initial operating point of the wind turbine, and solved by incremental partial derivative. The mechanical power increment can be obtained by integration as shown in formula (33). ΔP m =k ω Give r +k v Δv (33) Where k ω , k v After integration, Δω r , Δv corresponding coefficient. (2) Electromagnetic power model When analyzing the electromagnetic power output by the fan, the generator part of the fan can be equivalent to the first-order inertia link, and the electromagnetic power output by the generator can be expressed as: Where: P e P is the electromagnetic power output of the fan; ref is the power setting value of the fan electrical part; α is the first-order inertia constant; P ωr It is the MPPT active reference value calculated based on the speed; ΔP is the additional active power increment of the fan; K c It is the proportional coefficient of the maximum power tracking curve of the wind turbine. The fan electromagnetic power P can be obtained e The increment expression is: In summary, when considering the complete curve of the secondary frequency drop, the change in wind power is as shown in formula (36). Replacing the wind power variation in equation (29) with equation (36), we obtain the dynamic frequency constraint considering the secondary frequency drop.

5. The method and system for day-ahead dispatching of a power grid considering the frequency modulation and withdrawal dynamics of new energy sources according to claim 1, characterized in that: Benders decomposition algorithm (1) Model decomposition Since wind power is introduced into the power system dispatching problem to participate in the primary frequency regulation and the secondary frequency drop process, and the differential equation of the frequency is differentiated, the dynamic frequency constraint is too large and difficult to solve. Therefore, the optimization model is considered to be decomposed into a main problem and a sub-problem as shown in Figure 5. The main problem is the unit combination solution problem with conventional constraints, and the sub-problem is the frequency safety verification problem with dynamic frequency constraints. The main problem and the sub-problem are then connected through Benders cut. The overall problem to be solved is shown in formula (37). Wherein: x represents the variable related to the conventional constraint of the unit combination, y represents the variable related to the dynamic frequency constraint, the constraint Ax ≥ b represents the conventional constraint of the unit combination, which is composed of equations (17) to (24), and the constraint Ex + Fy ≥ h represents the frequency constraint, that is, equation (31). (2) Benders decomposition algorithm a) Model optimization process The model optimization process mentioned in the present invention consists of the following five steps: 1). After initialization, find the initial master problem containing only conventional constraints: 2). Solve the main problem with only conventional constraints to obtain a solution x* that only satisfies conventional constraints. Then introduce a penalty variable ε into the subproblem and substitute the obtained solution x* into the subproblem: 3) Dualize the subproblems to obtain: 4). Look at the value of the penalty variable in the dual subproblem to determine whether the solution obtained by the conventional constraint satisfies the frequency constraint. Ⅰ. If satisfied, the optimal solution is obtained and output Ⅱ. If not satisfied, proceed to step 5 5). By adding constraints to the main problem, the conventional constraint solutions that just do not meet the frequency constraint are eliminated, and the process is repeated until the optimal solution is obtained.

Images