Power system time sequence production simulation method considering equipment shutdown

Abstract

The invention relates to the technical field of power system production simulation, in particular to a power system time sequence production simulation method considering equipment shutdown. Comprising the steps of obtaining power grid information of a power system, and establishing a power grid structure model taking minimum power generation cost as a target function; performing mathematical modeling on operation constraints of various power supplies in a power grid according to operation principles of various power supplies; based on a power grid, a power supply, a load sequence and a new energy power generation output time sequence, establishing a power system operation mode optimization model suitable for new energy time sequence production simulation; and based on an operational research optimization algorithm, solving the optimization model to obtain power supply generated output data under the per-time section of the power system, and calculating the maximum consumption amount of the new energy. According to the method, time sequence production simulation of a power system is realized under the condition of considering two outage states of equipment failure and maintenance, and the problem of operation index evaluation of a power supply base and a direct-current transmission channel is solved.

Description

Technical Field

[0001] The present invention relates to the technical field of power system production simulation, and in particular to a chronological production simulation method for power systems considering equipment outages. Background Art

[0002] Production simulation is an important means for power system simulation analysis, especially for balance analysis. It can guide production, planning, and even market transactions by calculating indicators such as power supply gaps and unit utilization hours over a period of time. In traditional power systems, the installed capacity of new energy is relatively low, and the uncertainty of the power system mainly comes from power loads, with relatively little uncertainty. Usually, production simulation is based on continuous load curves. Such methods have low data volume requirements and exhibit good analysis accuracy and speed performance when the grid digital operation level and new energy penetration rate are relatively low. However, traditional chronological production simulation may have such a problem: due to the lack of consideration of transmission network power flow constraints, the new energy consumption capacity obtained from chronological production simulation may be overly optimistic. For example, at a certain moment, all new energy in the system is consumed, but at this time, the node voltage and branch power of some power flow sections may exceed the reasonable range, which makes the evaluation of the new energy consumption capacity of the system exceed the maximum new energy capacity that the actual power system can bear.

[0003] With the continuous increase in new energy penetration, the strong randomness and volatility of new energy output will lead to random changes and high-frequency fluctuations in power balance requirements. The grid operation mode changes frequently and rapidly, and the operation states of conventional units change frequently (such as rapid power ramping and frequent start-stop). However, the stochastic production simulation method based on continuous load curves is difficult to consider constraints such as ramp rates and unit start-stop. Therefore, it is urgent to utilize the massive data input and stronger computing power brought by the digital transformation of the power grid, consider more detailed operation constraints (especially chronological coupling constraints) and operation scenarios, and further innovate production simulation technology to meet the needs of new power systems, providing a panoramic and accurate decision-making reference for grid operation and planning.

[0004] Chronological production simulation can consider the chronological characteristics of new energy power generation output and the relevant constraints of unit operation by conducting hour-by-hour power balance of the power system at a long time scale and a fine time scale. It can model the flexibility of the power system, simulate the power and electricity balance process, and become an important simulation analysis tool in the power industry, providing reference opinions for grid companies, power generation enterprises, power users, and even energy and power policy-making institutions for operation and planning decisions. However, due to the large number of integer variables and complex system operation constraints (some are cross-period coupling constraints) in chronological production simulation, and in practical applications, different scenarios often need to be analyzed, the direct solution calculation is extremely time-consuming, resulting in limited application. In recent years, a large number of solution techniques for improving the calculation efficiency of chronological production simulation have been proposed and played an important role in practical applications. However, different solution techniques will simplify the model differently, resulting in corresponding application limitations when using different solution techniques. Therefore, considering the spatio-temporal distribution characteristics of renewable energy in deserts, gobi, and wastelands, proposing a chronological production simulation method for power systems considering equipment outages is of great significance for the further application of chronological production simulation. Summary of the Invention

[0005] Aiming at the deficiencies of the existing technology, the present invention intends to propose a chronological production simulation and power source optimization configuration method for power systems that simultaneously considers two outage situations: equipment maintenance and faults. The method includes:

[0006] Step S1: Obtain the grid information of the power system and establish a grid structure model with minimizing the generation cost as the objective function;

[0007] Step S2: Based on the operating principles of various power sources, conduct mathematical modeling of the operating constraints of various power sources in the grid;

[0008] Step S3: Based on the grid, power sources, load sequences, and new energy power generation output time series, establish an optimization model for the operation mode of the power system adapted to new energy chronological production simulation;

[0009] Step S4: Based on the optimization algorithm of operations research, solve the optimization model to obtain the power generation output data of power sources under each time section of the power system, and calculate the maximum accommodation amount of new energy.

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

[0011] (1) Realize the chronological production simulation of the power system considering two outage states: equipment faults and maintenance, and solve the problem of evaluating the operation indicators of power source bases and DC transmission channels.

[0012] (2) Further solve the problem of optimizing the configuration of the capacity of power source bases through the chronological production simulation method of the power system. Description of the Drawings

[0013] Figure 1 It is a schematic flowchart of the chronological production simulation method for a power system considering equipment outages in the embodiments of the present invention; Specific embodiments

[0014] The chronological production simulation method for a power system considering equipment outages in this embodiment is basically as Figure 1 shown, and the specific steps are as follows:

[0015] Step S1: Obtain the grid information of the power system and establish a grid structure model.

[0016] Step S2: According to the operating principles of various power sources, mathematically model the operating constraints of thermal power, hydropower, energy storage, and nuclear power sources in the grid;

[0017] Step S3: Based on the grid, power sources, load sequences, and new energy power generation output time series, establish an optimized model for the operating mode of the power system suitable for new energy chronological production simulation;

[0018] Step S4: Based on the optimization algorithm of operations research, solve the optimization model to obtain the power generation output data of the power sources at each time section of the power system, and calculate the maximum accommodation capacity of new energy.

[0019] The grid structure model in Step S1 specifically includes the following contents:

[0020] S11: Objective function of chronological production simulation

[0021] The main objective of the model is to optimize the operation of the power system by minimizing the generation cost. The objective function f total mainly includes the investment cost f inv , the operating cost f ope , the penalty term f ne,lost for renewable energy curtailment, and the penalty term f l,lost for load power outage. Considering that the cash value changes over time, calculate the equivalent annual cost of the planning scheme, and the calculation is as follows:

[0022] minf total = f inv + f ope + f ne,lost + f l,lost (1)

[0023] S12: Discount calculation

[0024] Convert the funds at different times into the equivalent amount at the current time, and this amount is called the present value. This conversion is called discount calculation, and the present value is also called the discounted value. The present value of the funds occurs at the beginning of the first year. Convert the funds into the equivalent amount at a future time, and this amount is called the future value. The future value of the funds is sometimes also called the final value. The final value of the funds occurs at the end of the last year. Convert the funds into the amount of equal payment per period, usually one year per period, and this amount is called the equivalent annual value. The equivalent annual value of the funds occurs at the end of each year.

[0025] Calculating the future value F from the present value P is called the calculation of principal and interest. Assuming the interest rate is i, the relationship between the future value F at the end of the nth year and the present value P is:

[0026] F n =P(1 + i) n (2)

[0027] Among them, (1 + i) n is the coefficient of the sum of principal and interest for a lump-sum payment.

[0028] When calculating using the above formula, it should be noted that the P value occurs at the beginning of the first year, and the future value occurs at the end of the nth year.

[0029] Calculating the present value P from the future value F is called discount calculation. From the relationship formula between the future value F and the present value P, we can get:

[0030]

[0031] Among them, is called the single-payment discount factor, which is the reciprocal of the single-payment compound amount factor; F n is the future value at the end of the nth year relative to the present value P.

[0032] Adopt the equivalent annual value method to calculate the standby-related costs, that is, convert the funds into the amount of equal payment per period, usually one year per period, and this amount is the equivalent annual value. The equivalent annual value of the funds occurs at the end of each year.

[0033] Calculating the future value F from the annual value A is called the calculation of the sum of principal and interest of the equivalent annual value. When the cash flow of the annual value A occurs at the end of each year from t = 1 year to t = n years, the future value F at the end of the nth year is equal to the sum of the future values of each A value in these n cash flows, that is:

[0034] F n =A + A(1 + i) + A(1 + i) 2 +… + A(1 + i) n-1 (4)

[0035] This is a problem of summing a geometric series, whose common ratio is 1 + i. Therefore, from the summation formula of the geometric series, we can get:

[0036]

[0037] Among them, is the equal annual value of principal and interest factor. This factor reflects the relationship between the equal annual value A for n years and the future value F at the end of the nth year.

[0038] The calculation of obtaining the present value P from the equal annual value A is called the present value calculation of the equal annual value.

[0039] From this, it can be obtained that:

[0040]

[0041] Among them, becomes the present value factor of the equal annual value.

[0042] Using the sinking fund factor, the amount of equal annual savings from now until the end of the nth year to pay for the expenses in the nth year can be calculated.

[0043] Using the equal annual value method to allocate the investment cost of the spare capacity r to the m (m = 1, 2,..., M) sampling years:

[0044]

[0045] Among them, f Σ (r) is the total investment cost of the spare capacity r; χ is the discount rate; J is the equipment life cycle. Thus, during the sampling process, f m (r) is a constant, m = 1, 2,..., M.

[0046] S13: Investment cost calculation

[0047] The specific calculation formula of the investment cost f in the objective function inv is as follows:

[0048]

[0049] In the formula, f inv,w , f inv,pv , f inv,g , f inv,ess are the investment costs of wind power, photovoltaic, thermal power units and energy storage; are the capacity investment cost coefficients of wind power, photovoltaic, thermal power units and energy storage; E w , E pv , E g , E ess are the installed capacities of wind power, photovoltaic, thermal power units and energy storage. Among them, the energy storage is calculated according to the power capacity, and the energy storage power capacity and energy capacity are uniformly taken as 1:2; T w , T pv , Tg , T ess is the service life of wind power, photovoltaic, thermal power units and energy storage, and χ is the discount rate.

[0050] S14: Operating cost calculation

[0051] The operating cost f in the objective function ope is calculated as follows:

[0052]

[0053] In the formula, f ope,w , f ope,pv , f ope,g , f ope,ess are the operating costs of the receiving-end wind power, photovoltaic and thermal power units, and f ope,L is the operating cost of the sending-end wind-solar-storage; is the operation and maintenance cost coefficient of the receiving-end wind power, photovoltaic and thermal power units, is the fuel cost coefficient of the receiving-end thermal power unit, c rp is the receiving-end landing comprehensive electricity price; P w , P pv , P g are the powers of the receiving-end wind power, photovoltaic and thermal power units; T is the operation period, taking 8760 h, and Δt is the simulation step size, taking 1 h; P flow23 is the active power transmitted between nodes 2 and 3.

[0054] S15: Renewable energy curtailment penalty term

[0055] According to the operation simulation results, the curtailment amount of renewable energy is counted, and the renewable energy curtailment penalty term is calculated as follows:

[0056]

[0057] In the formula, is the renewable energy curtailment penalty coefficient.

[0058] S16: Load power outage penalty term

[0059] According to the operation simulation results, the power outage amounts of each level of load are counted, and combined with the power outage penalty coefficient of the secondary load, the overall power outage penalty term of the system is obtained. The size of the penalty term not only reflects the economic losses brought by the power outage to users, but also can be used as an important reference index for optimizing the power system planning and operation strategy. The load power outage penalty term is calculated as follows:

[0060]

[0061] In the formula, is the power outage penalty coefficient of the secondary load.

[0062] In step S2 of this embodiment, the established operation constraint mathematical model is specifically as follows:

[0063] S21: Power supply reliability constraint

[0064]

[0065] In the formula: P loss is the load shedding amount in the system, representing the total load that needs to be reduced to ensure the safe and reliable power supply. The unit is megawatt (MW). f LPSP is the maximum allowable load shedding percentage, which is taken as 5% in this model and is an index to measure the power supply reliability of the system. is the simulated annual load data of the system, with the unit of megawatt (MW).

[0066] S22: Energy storage system operation constraint

[0067] The operation constraints of the energy storage system are divided into energy storage capacity change limit, energy storage state constraint, and energy storage periodic constraint. The specific expressions are as follows:

[0068] Expression of energy storage capacity change limit:

[0069]

[0070] In the formula: E ess,sk is the maximum energy storage capacity of the energy storage system, with the unit of megawatt-hour (MWh), which controls the charging and discharging rates of the energy storage system. ΔE sk is the energy change amount of the energy storage system in two consecutive time steps, reflecting the charging and discharging behavior of the energy storage system. t s,hour is the energy storage charging and discharging time, which is set to 2h in this embodiment.

[0071] Expression of energy storage state constraint:

[0072] SOC min ×E ess,sk ≤E sk ≤SOC max ×E ess,sk (14)

[0073] In the formula: SOC min and SOC max respectively represent the minimum and maximum state of charge of the energy storage system, which are the boundaries for the safe operation of the energy storage system. In this model, they are taken as 5% and 95% respectively. E ess,sk is the energy capacity of the energy storage system, with the unit of megawatt-hour (MWh), representing the maximum amount of electricity that the energy storage system can store. E sk is the energy state of the energy storage system at each time step, reflecting the charging and discharging process of the energy storage system.

[0074] Energy storage periodic constraint expression:

[0075] E sk (1) = E sk (N) (15)

[0076] S23: Spinning reserve constraint

[0077]

[0078] Where: Pgmax(i) is the upper limit of the active power output of unit i. In this embodiment, P re is the load reserve coefficient, which is taken as 5% in this model. X(i, t) is a binary variable representing the operating state of unit i at time t. When X(i, t) equals 1, it means unit i is in the operating state at time t. When X(i, t) equals 0, it means unit i is in the shutdown state at time t.

[0079] S24: Conventional unit operation constraint

[0080] The conventional unit operation constraints include: unit output power limit, unit ramp rate constraint, minimum on / off time constraint, and unit status logic constraint.

[0081] Unit output power limit expression:

[0082] P gmin(i) X(i, t) ≤ P g (i, t) ≤ P gmax(i) X(i, t) (17)

[0083] Where: P gmin (i) is the minimum technical output of unit i, and P gmax (i) is the maximum technical output of unit i.

[0084] Unit ramp rate constraint expression:

[0085]

[0086] Where: P up (i) and P down (i) are the maximum allowable uphill rate and maximum allowable downhill rate of unit i, respectively.

[0087] Minimum on / off time constraint expression:

[0088]

[0089] Where: Y(i,t) and Z(i,t) respectively represent the binary variables of the start-up and shutdown states of unit i at time t. When Y(i,t) = 1, it means the unit is starting at time t; when Y(i,t) = 0, it means the unit is not in the start-up state. When Z(i,t) = 1, it means the unit is stopping its operation at time t; when Z(i,t) = 0, it means the unit is not performing a shutdown operation; T on and T off are respectively the minimum continuous operation time and the minimum continuous shutdown time of the unit, and are taken as 2h in this model.

[0090] Unit status logical constraint expression:

[0091]

[0092] S25: Renewable energy output constraint

[0093]

[0094] Where: and respectively represent the theoretical maximum output of wind power and photovoltaic power at the nth port.

[0095] S26: Transmission line capacity constraint

[0096]

[0097] Where: In this embodiment, it is considered that the sending ends 1 and 2 are at the same site, so there is no DC transmission line constraint. P flow23 is the active power between the sending end node 2 and the receiving end node 3. P flow34 is the active power between the receiving end nodes 3 and 4. P max is the maximum active power of the DC transmission channel, and is taken as 10000MW in the model established in this embodiment.

[0098] S27: Node power balance constraint

[0099]

[0100] The above formula is the power balance formula for 3 nodes at the sending end and the receiving end, that is: the power balance between the sending end node 2 and the receiving end node 3, the injection power of the receiving end node 3 plus the power received from the sending end nodes 1 and 2 is equal to zero, and the injection power of the receiving end node 4 plus the power received from the receiving end node 3 is equal to zero; P flow23 is the active power transmitted from the sending end nodes 1 and 2 to the receiving end node 3; P wp (n) and P pv (n) respectively represent the wind power and photovoltaic power output at the nth port; E sk(n) represents the energy storage of the nth port; P load,n represents the load of the nth port.

[0101] S28: Converter station capacity constraint

[0102]

[0103] In the formula, the meanings of the variables are the same as those in S27, is half of the maximum active power of the DC transmission channel, that is, 5000 MW.

[0104] The specific implementation process of step S3 in this embodiment is as follows:

[0105] S31: In the capacity optimization configuration model considering economy, the objective function used is the same as the objective function in step S2, and the specific function content can refer to formulas (1) to (11). The newly added target variables to be optimized are as follows: E w,1 is the installed capacity of the wind turbine generator at port 1, E pv,1 is the installed capacity of the photovoltaic generator at port 1, E ess,1 is the installed capacity of the energy storage at port 1, E esp,1 is the total installed capacity at port 1, E w,2 is the installed capacity of the wind turbine generator at port 2, E pv,2 is the installed capacity of the photovoltaic generator at port 2, E ess,2 is the installed capacity of the energy storage at port 2, E esp,2 is the total installed capacity at port 2, E esp is the total installed capacity of port 1 and port 2. Then there is:

[0106] E esp = E esp,1 + E esp,2 (25)

[0107] S32: Add the constraint conditions of the capacity optimization configuration system

[0108] (1) Consider the output constraint of conventional units for faults and maintenance outages

[0109] On the basis of S24, considering the fault state of thermal power units, the original output constraint conditions of conventional units are updated:

[0110] P gmin(i) X(i,t) ≤ P g (i,t) ≤ P gmax(i) X(i,t).*X g,i (26)

[0111] In the formula, X g,i is the fault state sequence of the thermal power unit at the ith node.

[0112] (2) Consider the renewable energy output constraints considering faults and maintenance outages

[0113] Based on S25, update the original renewable energy output constraint conditions. The updated expression is as follows:

[0114]

[0115] In the formula, X wp,n represents the wind power output status sequence of the nth node, and X pv,n represents the photovoltaic output status sequence of the nth node.

[0116] (3) Consider the transmission line capacity constraints considering faults and maintenance outages

[0117] Based on S26, update the transmission line capacity constraints. The updated expression is as follows:

[0118]

[0119] In the formula, X line1 , X line2 are the line status sequences of the transmission lines between ports 2, 3 and ports 3, 4 respectively.

[0120] (4) Consider the energy storage configuration capacity constraints

[0121] Since the energy storage configuration capacity should meet the regulations of the local power grid energy bureau, among which the requirements for the guaranteed projects of the Inner Mongolia Autonomous Region Energy Bureau are: the energy storage configuration capacity should not be less than 15% of the total installed capacity of wind and light, and the energy storage duration should be more than 2 hours. Therefore, in this embodiment, the energy storage capacity configuration constraint is set to be greater than or equal to 15% of the total wind and light capacity.

[0122]

[0123] In the formula, α represents the set minimum energy storage ratio. In this embodiment, α is 15%, E esp,w is the total installed capacity of wind power, and E esp,pv is the total installed capacity of photovoltaic power.

[0124] (5) Transmission channel annual utilization hours ≥ 4500 hours constraint (additional constraint)

[0125] This embodiment clearly requires a power source optimization configuration scheme with the constraint that the annual utilization hours of the transmission channel ≥ 4500 hours. Therefore, when proposing a recommended scheme for power source optimization configuration, in addition to the above conventional constraint conditions, it is also necessary to additionally consider the constraint condition that the annual utilization hours of the DC transmission channel ≥ 4500 hours, as shown in the following formula.

[0126]

[0127] Where P flow23 It is the power flowing through the DC transmission channel connecting the sending-end converter station and the receiving-end converter station.

[0128] In step S4 of this embodiment, the new energy time-series production simulation model can be mathematically reduced to solving a mixed integer linear programming problem, and its mathematical model is abbreviated as follows:

[0129]

[0130] Where x represents the set of variables to be optimized; f(x) represents the optimization objective function; g i (x)≥0 represents the inequality constraint set; h j (x) = 0 represents an equality constraint set; f(x), g i (x) and h j (x) are all linear expressions of decision variables.

[0131] The core algorithm for solving mixed integer programming is the branch and bound method (B&B), whose basic idea is to search the space of all feasible solutions (limited number) of optimization problems with constraints. For a mixed integer programming problem A with maximization, the integer variables are relaxed into continuous variables to obtain the corresponding programming problem B. Starting from solving problem B, if its optimal solution does not meet the integer condition of A, then the optimal objective function of B must be the upper bound of the optimal objective function of A, and the objective function value of any integer feasible solution of A is the lower bound. According to the optimization result of problem B, the variable value range is divided by integers to obtain branch problems B1 and B2 of B. Solve B1 and B2, and continue to delineate branches according to the results. Iterate in sequence until a solution that meets the constraints of A is generated in the branch problem, which is the optimal solution. The branch and bound method is to divide the feasible domain of B into sub-regions, continuously branch, prune and delimit. After finding a better feasible integer solution, update the lower bound, gradually increase the lower bound and reduce the upper bound, and finally obtain the optimal solution of problem A.

[0132] For nonlinear mixed integer programming problems such as unit power output optimization, while the branch and bound method is reasonably applied to split the original mixed integer programming problem, the interior point method or exterior point method is also used for nonlinear optimization calculation. In order to ensure the speed of solution and the global optimality of the results, the optimization calculation code needs to be written scientifically and reasonably, and various algorithms should be organically combined.

[0133] Common integer programming models include the knapsack problem, set covering, packing and partitioning problems, at least K of L constraints must be satisfied, a function that takes L values, the fixed cost problem, If-Then constraints, and piecewise linear functions. The following mainly introduces If-Then constraints.

[0134] In many applications, if the constraint f(x1, x2,..., x n ) > 0 is satisfied, then the constraint g(x1, x2,..., x n ) ≥ 0 must also be satisfied. To ensure this, a 0-1 variable y can be introduced. When f(x1, x2,..., x n ) > 0, y = 0, and then it is required that when y = 0, g(x1, x2,..., x n ) ≥ 0. Thus, this requirement can be expressed as:

[0135]

[0136] where M represents a sufficiently large constant that should ensure that all (x1, x2,..., x n ) that satisfy the other constraints in the problem also satisfy f(x1, x2,..., x n ) ≤ M and -g(x1, x2,..., x n ) ≥ M.

[0137] It can be seen that if f > 0, then necessarily y = 0. Then, from the constraint conditions, it can be known that -g < 0 or g > 0, which is the required result.

[0138] Let the matrix Α ∈ R m×n , the column vector b ∈ R m , b ∈ R n , x ∈ R n , and the set of variable subscripts Then the general mixed integer programming problem (Mixed Integer Programming, MIP) can be expressed as:

[0139]

[0140] where l j , u j ∈ R ∪ {±∞}. Denote the set of subscripts of 0-1 variables as B = {j ∈ I | l j = 0 and u j = 1}, and the set of subscripts of continuous variables as C = N / I.

[0141] If Then the optimization problem in Equation (33) is called a linear programming problem (LP); if I = N, it is called an integer programming problem (IP); if I = B, it is called a mixed binary programming problem (MBP); if I = B = N, it is called a binary programming problem (BP); if the integer constraints are removed, the linear programming problem can be obtained:

[0142]

[0143] To obtain the optimal solution of the integer programming, the simplest algorithm is to use the enumeration method (which means a huge amount of computation). To minimize the amount of computation as much as possible, the branch and bound method is proposed, which is actually a kind of implicit enumeration method.

[0144] Consider the general linear mixed-integer programming problem. Let its feasible region be S and the optimal objective value be Z * , and the feasible region of its linear programming relaxation problem is P.

[0145] Let S = S1 ∪ S2 ∪ … ∪ S k be a decomposition of the feasible region S. Let z i = min{c T x: x ∈ S i}, i = 1, …, k, z * = min{c T x: x ∈ S}, then

[0146] Let S = S1 ∪ S2 ∪ … ∪ S k be a decomposition of the feasible region S. Let z i = min{c T x: x ∈ S i}, i = 1, …, k, z * = min{c T x: x ∈ S}, An upper bound of the expression z i (that is, ), z i is a lower bound of z i (that is, z i ≤ z i ). Then is an upper bound of z * , is a lower bound of z * .

[0147] The branch-and-bound method based on linear programming problems consists of a series of iterations. Let the feasible domain of the mixed integer programming problem be S, and the feasible domain of its relaxation problem be J. Initially, let П1 = {J}, and assume that a candidate solution has been obtained. The corresponding target value is (If there is no candidate solution, then ). In the subsequent iterations, Δ is decomposed into a series of pairwise disjoint polyhedral sets J1,…,J k . Assume that at the kth iteration, a decomposition Π of J has been obtained k ={J1,…,J k},satisfy:

[0148] (1){J1,…,J k}isI k There are two disjoint polyhedral sets in the set (each of which can be represented by a set of inequalities) and J = J1∩J2∩…∩J k .

[0149] (2) All feasible solutions to the MIP problem are contained in the set J1,…,J k If S i =S∩J i , i=1,...,k, then S i is the feasible domain of the ith subproblem of the mixed integer programming problem, and {S i :i=1,...,k} constitutes a decomposition of the feasible domain S of the mixed integer programming problem.

[0150] Then from Π k Choose a subproblem to solve, let’s start with Π k Select subproblem J i (That is, the feasible domain is J i The i-th subproblem of this problem has a feasible domain of S i The linear programming relaxation problem of the ith MIP subproblem of . The subproblems correspond to the feasible domain one by one, and the two will no longer be distinguished in the future), and the linear programming method is used to solve it. The following situations may occur:

[0151] 1) Polyhedral set J i is an empty set. Then we can solve the subproblem J i From k The removal of branches from the tree is called pruning due to infeasibility.

[0152] 2) The optimal value of this problem At this time, in the feasible region J i There will not be a solution δ + (k) A better feasible solution, subproblem J i From kDeleted from it, which is called pruning according to the boundary.

[0153] 3) The optimal solution x of this problem i is a feasible solution to the mixed-integer programming problem, that is This indicates that the i-th sub-problem z i = min{c T x|x ∈ S i} has been solved. At this time, there must be Therefore, it is necessary to update the candidate solution δ - (k) and That is, let Then delete the sub-problem J i from Π k This is called pruning according to optimality.

[0154] In the above three pruning cases, we all get Then enter the (k + 1)-th iteration.

[0155] If it is not the above case, then the optimal solution x of the sub-problem Δ j does not satisfy the integer constraint requirements of the mixed-integer programming problem. Without loss of generality, let x i = {ξ1,…,ξ n}, ξ j (j ∈ I) is a fraction, then define:

[0156]

[0157]

[0158] And let Then enter the (k + 1)-th iteration.

[0159] The above is the basic iterative process of the branch and bound method. The iteration will continue until before it terminates. Therefore, the branch and bound method is essentially an enumeration method and may generate a huge amount of computation. However, if a suitable search strategy is adopted and combined with a good bound, the amount of computation can also be greatly saved. In the search, there are the following two basic problems to be solved:

[0160] (1) In the case where multiple ξ j (j ∈ I) are non-integer, which variable should be selected as the branching variable. This problem is called the branching problem.

[0161] (2) Which problem should be solved first in the branching problem. This choice is called the node selection problem.

[0162] For the first problem, that is, if there are two or more integer variables in the optimal solution of a sub-problem that are fractional, if it is known in advance which variable is more important, then the most important variable should be selected for branching. Otherwise, an estimate should be made for the integer variables taking fractional values to see which one is more important. This method of selecting branching variables is called the branching strategy.

[0163] For the second problem, that is, which problem to solve first among the leaf nodes of the branch, the commonly used rules include last in, first out (LIFO) or depth-first search (DFS), first in, first out (FIFO) or breadth-first search (BFS) rules, or a compromise between the two. These methods of selecting nodes are called the node selection strategy.

[0164] (1) Branching strategy. Now denote the feasible region of the MIP problem as X. The basic algorithm for selecting branching variables is as follows:

[0165] Input: The current sub-problem Q and the optimal solution of its linear programming relaxation

[0166] Output: The subscript j of an integer variable taking a fractional value, that is, j ∈ I, but

[0167] Basic algorithm for selecting branching variables:

[0168] Let F denote the set of subscripts that can be used as branching variables, that is, let

[0169] For all j ∈ F, calculate an evaluation value s j ∈ R, let p = argmax k∈F {s k} (if there are multiple, then select one of them in a certain way), then select x p as the branching variable and return the subscript value J.

[0170] Commonly used methods for selecting branching variables include the most "fractional" variable branching method (most infeasible branching), the least infeasible fractional value variable branching method (least infeasible branching), the pseudo-cost branching method (pseudo-cost branching), the strong branching method (strong branching), and the hybrid branching method.

[0171] (2) Node selection strategy. In the branch-and-bound method, after solving a node (sub-problem), the next step is to select a leaf node in the current search tree as the sub-problem to continue the branch-and-bound solution process. In the branch-and-bound search tree of MIP, the selection of nodes usually requires achieving the following two opposite goals: finding a good MIP feasible solution to improve the original bound (primal bound, which is the upper bound for minimization problems and the lower bound for maximization problems), which can be used to prune the search tree relatively quickly; improving the global dual bound (dual bound, which is the lower bound for minimization problems and the upper bound for maximization problems).

[0172] In addition to finding an MIP feasible solution through the original heuristic method, when solving the linear programming relaxation problem of a node, if the corresponding relaxed optimal solution is also MIP feasible, then an MIP feasible solution is found. In actual calculations, such a situation usually occurs deep in the search tree. Therefore, in order to find an MIP feasible solution more quickly, it is natural to adopt the depth-first search strategy (DFS). However, a method like depth-first completely ignores the second goal because the nodes with the smallest lower bounds are usually near the root node of the search tree. In order to enhance the global dual bound as much as possible, the best-first search strategy should be used. For minimization problems, this method always selects the leaf node with the smallest objective function value to solve. In order to achieve the above two goals simultaneously, a strategy that combines these two methods is generated, called the best-first search strategy with plunging.

[0173] A variant of the best-first search strategy is the best-estimate search strategy. Instead of selecting the node with the best dual bound to solve, this method estimates the optimal objective function value of the corresponding node in the MIP feasible region and then selects the node with the best estimate to solve. The optimal objective function value of a node in the MIP feasible region is estimated by the dual bound of the node, the degree to which the integer variables take fractional values, and the pseudo-cost values of the corresponding variables. The goal of the best-estimate search strategy is to find a good or even optimal feasible solution as early as possible. If this strategy is combined with the depth-first strategy, the best-estimate search strategy with plunging is obtained. In addition, the above three methods can also be used to generate different search strategies.

Claims

1. A chronological production simulation method for a power system considering equipment outages, characterized in that Including: Step S1: Obtain the grid information of the power system and establish a grid structure model with the goal of minimizing the generation cost. Step S2: Based on the operating principles of various power sources, mathematically model the operating constraints of thermal power, hydropower, energy storage, and nuclear power sources in the grid. Step S3: Based on the grid, power sources, load sequence, and new energy generation output time series, establish an optimization model for the operation mode of the power system adapted to the new energy time-series production simulation. Step S4: Based on the optimization algorithm of operations research, solve the optimization model to obtain the power generation output data of the power sources at each time section of the power system, and calculate the maximum accommodation capacity of the new energy.

2. The method according to claim 1, characterized in that, The grid structure model in Step S1 specifically includes the following: S11: The objective function f of the sequential production simulation total , there is minf total = f inv + f ope + f ne,lost + f l,lost (1) f inv is the investment cost, f ope is the operating cost, f ne,lost is the penalty term for curtailment of renewable energy, f l,lost is the penalty term for load power outage; S12: Discount calculation It is called the calculation of the sum of principal and interest to calculate the future value F from the present value P. Assuming the interest rate is i, the relationship between the future value F at the end of the nth year and the present value P is: F n = P(1 + i) n (2) Among them, (1 + i) n is the coefficient of the sum of principal and interest for a lump-sum payment; When calculating using the above formula, it should be noted that the P value occurs at the beginning of the first year, and the future value occurs at the end of the nth year. Calculating the present value P from the future value F is called discount calculation. From the relationship formula between the future value F and the present value P, we can get: Among them, is called the single-payment discount factor and is the reciprocal of the single-payment compound-amount factor; F n is the future value at the end of the nth year relative to the present value P; The equal annual value method is used to calculate the standby-related costs, and the equal annual value of the funds occurs at the end of each year. Calculating the future value F from the annual value A is called the calculation of the sum of principal and interest of the equal annual value. When the cash flow of the annual value A occurs at the end of each year from t = 1 year to t = n years, the future value F at the end of the nth year is equal to the sum of the future values of each A value in these n cash flows, that is: Among them, is the equal annual value of principal and interest coefficient; Calculating the present value P from the equal annual value A is called the present value calculation of the equal annual value. From this, we can get: Among them, becomes the present value factor of equal annual value; Use the sinking fund coefficient to calculate the amount of equal annual savings from now until the end of the nth year in order to pay the expenses in the nth year. Use the equal annual value method to allocate the investment cost of the standby capacity r to the m (m = 1, 2,..., M) sampling years: where f Σ (r) is the total investment cost of the spare capacity r; χ is the discount rate; J is the equipment life cycle, such that f m (r) is a constant, m = 1, 2,..., M; S13: Investment cost calculation Investment cost f inv The specific calculation formula is as follows: Where f inv,w , f inv,pv , f inv,g , f inv,ess are the investment costs of wind power, photovoltaic, thermal power generation units and energy storage; are the capacity investment cost coefficients of wind power, photovoltaic, thermal power generation units and energy storage; E w , E pv , E g , E ess are the installed capacities of wind power, photovoltaic, thermal power generation units and energy storage. Among them, the energy storage is calculated according to the power capacity, and the energy storage power capacity and energy capacity are uniformly taken as 1:2; T w , T pv , T g , T ess are the service lives of wind power, photovoltaic, thermal power generation units and energy storage, and χ is the discount rate; S14: Operating cost calculation Operating cost f ope The specific calculation formula is as follows: where, f ope,w , f ope,pv , f ope,g , f ope,ess are the operating costs of the receiving-end wind power, photovoltaic and thermal power units, and f ope,L is the operating cost of the sending-end wind-solar-storage; is the operation and maintenance cost coefficient of the receiving-end wind power, photovoltaic and thermal power units, is the fuel cost coefficient of the receiving-end thermal power unit, c rp is the receiving-end landing comprehensive electricity price; P w , P pv , P g are the powers of the receiving-end wind power, photovoltaic and thermal power units; T is the operating period, and Δt is the simulation step size; P flow23 is the active power transmitted between nodes 2 and 3; S15: Renewable energy curtailment penalty term According to the operation simulation results, count the curtailment amount of renewable energy, and calculate the renewable energy curtailment penalty term as follows: In the formula, is the penalty coefficient for the curtailment of renewable energy; S16: Load power outage penalty term According to the operation simulation results, count the power outage amounts of each level of load, and combine the power outage penalty coefficient of the secondary load to obtain the overall power outage penalty term of the system as follows: In the formula, is the power outage penalty coefficient for secondary loads.

3. According to the method described in claim 1, the operating constraint mathematical model established in step S2 is specifically as follows: S21: Power supply reliability constraint Where: P loss is the load shedding amount in the system, representing the total load that needs to be reduced to ensure the safe and reliable power supply, with the unit of megawatt, f LPSP is the maximum allowable load shedding percentage, is the annual load data of the simulated system, with the unit of megawatt; S22: Energy storage system operating constraint The operating constraints of the energy storage system are divided into energy storage capacity change limit, energy storage state constraint, and energy storage periodic constraint; Expression of energy storage capacity change limit: Where: E ess,sk is the maximum energy storage capacity of the energy storage system, in megawatt-hours; ΔE sk is the energy change of the energy storage system in two consecutive time steps, and t s,hour is the charge and discharge time of the energy storage; Expression of energy storage state constraint: SOC min ×E ess,sk ≤E sk ≤SOC max ×E ess,sk (13) Where: SOC min and SOC max represent the minimum and maximum state of charge of the energy storage system respectively, E ess,sk is the energy capacity of the energy storage system, in megawatt-hours, and E sk is the energy state of the energy storage system at each time step; Expression of energy storage periodic constraint: E sk (1) = E sk (N) (14) S23: Spinning reserve constraint where: Pgmax(i) is the upper limit of the active power output of unit i, P re is the load reserve coefficient, and X(i,t) represents the binary variable of the operating state of unit i at time t. When X(i,t) equals 1, it means that unit i is in the operating state at time t, and when X(i,t) equals 0, it means that unit i is in the shutdown state; S24: Conventional unit operating constraint The operating constraints of conventional units are: unit output power limit, unit ramp rate constraint, minimum on / off time constraint, unit status logic constraint; Expression of unit output power limit: P gmin(i) X(i,t) ≤ P g (i,t) ≤ P gmax(i) X(i,t) (16) Where: P gmin (i) is the minimum technical output of unit i, P gmax (i) is the maximum technical output of unit i; Expression of unit ramp rate constraint: Where: P up (i) and P down (i) are respectively the maximum uphill rate and the maximum downhill rate allowed for unit i; Expression of minimum on / off time constraint: Where: Y(i,t) and Z(i,t) respectively represent binary variables of the start-up and shutdown states of unit i at time t. When Y(i,t) = 1, it means the unit is starting up at time t; when Y(i,t) = 0, it means the unit is not in the start-up state at time t. When Z(i,t) = 1, it means the unit is stopping operation at time t; when Z(i,t) = 0, it means the unit is not performing a shutdown operation at time t; T on and T off are respectively the minimum continuous operation time and the minimum continuous shutdown time of the unit; Expression of unit status logic constraint: S25: Renewable energy output constraint Where: and respectively represent the theoretical maximum output of wind power and photovoltaic power at the nth port; S26: Transmission line capacity constraint Where: P flow23 is the active power between the sending end node 2 and the receiving end node 3, P flow34 is the active power between the receiving end nodes 3 and 4, P max is the maximum active power of the DC transmission channel; S27: Node power balance constraint The above formula is the power balance formula for three nodes at the sending end and the receiving end, that is: the power balance between node 2 at the sending end and node 3 at the receiving end, the injected power of receiving end node 3 plus the power received from sending end nodes 1 and 2 equals zero, and the injected power of receiving end node 4 plus the power received from receiving end node 3 equals zero; P flow23 is the active power transmitted from sending end nodes 1 and 2 to both ends of receiving end node 3; P wp (n) and P pv (n) respectively represent the wind power and photovoltaic output of the nth port; E sk (n) represents the energy storage of the nth port; P load,n represents the load of the nth port; S28: Converter station capacity constraint 4. The method according to claim 3, wherein The specific implementation process of step S3 is as follows: S31: Adopt the same objective function as in step S2, and additionally add the following target variables to be optimized: E w,1 is the installed capacity of the wind turbine generator for Port 1, E pv,1 is the installed capacity of the photovoltaic unit for Port 1, E ess,1 is the installed capacity of the energy storage for Port 1, E esp,1 is the total installed capacity for Port 1, E w,2 is the installed capacity of the wind turbine generator for Port 2, E pv,2 is the installed capacity of the photovoltaic unit for Port 2, E ess,2 is the installed capacity of the energy storage for Port 2, E esp,2 is the total installed capacity for Port 2, E esp is the total installed capacity of Port 1 and Port 2, then there is: E esp = E esp,1 + E esp,2 (24) S32: Add the constraint conditions of the capacity optimization configuration system (1) Consider the output constraints of conventional units during faults and maintenance outages On the basis of S24, considering the fault states of thermal power units, the original output constraint conditions of conventional units are updated as follows: P gmin(i) X(i,t) ≤ P g (i,t) ≤ P gmax(i) X(i,t).*X g,i (25) where X g,i is the fault status sequence of the thermal power unit at the i-th node; (2) Consider the output constraints of renewable energy during faults and maintenance outages On the basis of S25, the original output constraint conditions of renewable energy are updated as follows: where X wp,n represents the wind power output status sequence of the nth node, and X pv,n represents the photovoltaic power output status sequence of the nth node; (3) Consider the capacity constraints of transmission lines during faults and maintenance outages On the basis of S26, the capacity constraints of transmission lines are updated as follows: wherein, X line1 , X line2 are respectively the line state sequences of the transmission lines between port 2 and port 3 and between port 3 and port 4; (4) Consider the capacity constraints of energy storage configuration where α represents the set minimum energy storage ratio, E esp,w is the total installed capacity of wind power, and E esp,pv is the total installed capacity of photovoltaic power; (5) Annual utilization hours constraint of transmission channels Where P flow23 is the power flowing through the HVDC transmission channel connecting the sending converter station and the receiving converter station.

5. The method according to claim 4, wherein In step S4, during the process of solving the optimization model based on the optimization algorithm of operations research, a best-first search strategy that takes into account depth is adopted.

Images