CN105809297B - A kind of thermal power plant's environmental economy dispatching method based on multiple target differential evolution algorithm - Google Patents

Abstract

The present invention discloses a kind of thermal power plant's environmental economy dispatching method based on multiple target differential evolution algorithm, comprising the following steps: establishes using generate electricity network minimal and the minimum regulation goal of pollutant discharge amount, using generator capacity and power-balance as constraint condition thermal power plant environmental economy scheduling model；Multiple target differential evolution algorithm is recycled to optimize model, obtain optimal Pareto disaggregation, wherein multiple target differential evolution algorithm is scanned for using differential variation operator, access times and accumulation performance based on operator each when making a variation several times recently when each mutation operation select mutation operator, and ensure the convergence of disaggregation and the uniformity of distribution using non-dominated ranking, domination the methods of number and hypervolume contribution amount；Decision finally is carried out using Fuzzy Set Theory, selects compromise solution as final scheduling scheme from Pareto solution concentration.The characteristic that there is the method for the present invention precision height, Pareto forward position disaggregation to be evenly distributed with fast convergence rate, and it is easy to Project Realization.

Description

Thermal power plant environment economic dispatching method based on multi-target differential evolution algorithm

Technical Field

The invention relates to the technical field of power system scheduling, in particular to a thermal power plant environment economic scheduling method based on a multi-objective differential evolution algorithm.

Background

The economic dispatching of the power system is a dispatching scheme for solving the lowest power generation cost on the basis of meeting the operation constraint condition of the power system. In the case of a thermal power generator, a large amount of pollution gases such as sulfur oxides, nitrogen oxides, and carbon dioxide, or greenhouse gases are emitted during power generation. If the pollutant discharge amount in the power generation process is considered, the original single-target economic dispatching is changed into multi-target environmental economic dispatching. The difficulty in formulating a scheduling scheme is increased due to the conflicting cost and pollutant emission objectives.

The methods for solving the environmental economic dispatch can be generally divided into two categories, namely a mathematical programming method and an intelligent optimization method. The mathematical programming method mainly adopts means such as a constraint method, a weight coefficient sum method and the like to convert a multi-objective scheduling problem into a single-objective scheduling problem, and obtains a pareto solution set through multiple solving. The pareto solution set can be obtained by single solution of the intelligent optimization method, the pareto solution set has the characteristics of simplicity in implementation, high optimization efficiency, suitability for solving nonlinear problems and the like, and is widely applied to solving of environmental economic scheduling problems, such as genetic algorithm, particle swarm algorithm, bacterial foraging algorithm, random search and the like. The differential evolution algorithm is an intelligent optimization technology with simplicity, parallelism and strong robustness, carries out optimization in a search space through special variation operation, and is an effective means for solving the economic dispatching of the environment. In addition, no matter a mathematical programming method or an intelligent optimization method, after a pareto solution set is obtained, how to make a reasonable decision to select a certain compromise solution from the solution set as a final implementation scheme is also one of the keys for realizing the environmental-economic scheduling of the thermal power plant.

At present, some researches are carried out on the environmental economic dispatching of the thermal power plant by the academic community, and partial results are obtained, such as ' power system environmental economic dispatching based on multi-objective evolutionary algorithm ' by zhuyingsheng et al, and ' power environmental economic dispatching adopting multi-objective evolutionary algorithm based on decomposition ' (power environmental economic dispatching adopting multi-objective improved differential evolutionary algorithm) ', and ' environmental economic power generation dispatching adopting multi-objective improved differential evolutionary algorithm ' by hu and et al.

Disclosure of Invention

Aiming at the problems, the invention aims to provide a thermal power plant environment economic dispatching method based on a multi-objective differential evolution algorithm, overcomes the defect that the conventional method is difficult to obtain a pareto solution set with optimal convergence and distribution uniformity when solving the problem, and carries out reasonable decision to realize selection of a final implementation scheme.

In order to solve the technical problem, the invention adopts the technical scheme that:

a thermal power plant environment economic dispatching method based on a multi-target differential evolution algorithm comprises the following steps:

step 1: establishing mathematical model for environmental economic dispatch of thermal power plant

The mathematical model of the economic dispatch of the thermal power plant comprises an objective function and a constraint condition, wherein the objective function is described by the following formula:

wherein K is the number of generators in the thermal generator set, P_kFor the kth generator to output active power, F_k(P_k) The active power output for the kth generator is P_kThe required electricity generation cost of time, E_k(P_k) The output active power of the ith generator is P_kThe discharge amount of pollutants, K ═ 1, 2., K, the power generation cost F_k(P_k) The calculation expression of (a) is:

F_k(P_k)＝a_k+b_kP_k+c_kP_k ²，

wherein, a_k、b_kAnd c_kThe cost coefficient of the kth generator, K is 1,2, and K, describes the cost of fuel, labor, equipment maintenance and the like required by the generator to generate electricity, and is a known parameter for a specific generator;

pollutant emission E_k(P_k) The calculation expression of (a) is:

E_k(P_k)＝α_k+β_kP_k+γ_kP_k ²+ξ_kexp(λ_kP_k)，

wherein, α_k、β_k、γ_k、ξ_kAnd λ_kThe characteristic parameters of pollutant emission of the kth generator are respectively, and K is 1,2, a.

The constraint conditions of the environmental economic dispatching of the thermal power plant comprise generator capacity constraint and power balance constraint; the generator capacity constraint is described by:

wherein,andrespectively outputting a minimum value and a maximum value of active power of a kth generator, wherein K is 1, 2. The power balance constraint is that the total output power of the generator set is equal to the sum of network loss and load power, and the expression is as follows:

wherein,the total output active power P of K generators in the thermal generator set_DFor system load power, P_LPower is lost to the network. P_DIs a known parameter, P_LCan be calculated by the B coefficient method as follows:

wherein, B_ij、B_0iAnd B₀₀For the network loss factor, i 1,2, a, K, j 1,2, K, a parameter known for the particular power system;

step 2: obtaining various parameters required in the model in the step 1, including: total load power P of system in power grid_DCalculating the parameter B of the network power loss_ij、B_0iAnd B₀₀Minimum and maximum output active power of each generator in thermal generator setAndcost parameter a of generator set_k、b_kAnd c_iAnd pollutant discharge amount parameter α of generator set_k、β_k、γ_k、ξ_kAnd λ_k，i＝1,2,...,K，j＝1,2,...,K，k＝1,2,...,K；

And step 3: the mathematical model in the step 1 is optimized and solved by adopting a multi-objective differential evolution algorithm to obtain a pareto solution setN is the number of solution concentration scheduling schemes, and solution concentration P^g,iRepresents the ith scheduling scheme, i 1, 2.., N;

and 4, step 4: on the basis of the pareto solution set obtained in the step 3, the satisfaction degree of each scheduling scheme is calculated by using a fuzzy set theory method, and the calculation expression is as follows:

wherein, mu_nRepresenting the overall satisfaction of the nth scheduling scheme in the pareto solution set,representing the satisfaction of the cost objective in the nth scheduling scheme,represents the satisfaction of the pollutant discharge amount in the nth scheduling scheme, wherein N is 1, 2.Andcan be calculated from the following equations:

wherein, F_nAnd E_nRespectively representing the power generation cost and the pollutant discharge amount of the nth scheduling scheme, F_minAnd F_maxRespectively representing the minimum and maximum values of the electricity generation cost in all the N scheduling schemes, E_minAnd E_maxRespectively representing the minimum value and the maximum value of pollutant discharge amount in all N scheduling schemes;

and 5: the satisfaction degree mu of the scheduling scheme obtained in the step 4_nThe numerical values are compared and sequenced, the scheduling scheme with the maximum satisfaction numerical value is selected as the final scheduling scheme, and the numerical value of the selected scheduling scheme is corresponding to the pareto frontier distribution map drawn according to the method.

Further, the various model parameters obtained in step 2 depend on specific thermal power plant generator sets, and for some thermal power plant generator sets, the index item ξ in the scheduling objective function value is not considered when the scheduling objective function value is calculated_kexp(λ_kP_k) Then ξ can be ordered_kAnd λ_k0, j 1,2, K, which uses the parameters obtained to model a specific thermal power plant.

Further, the step 3 comprises the following steps:

step 3.1, initializing parameters NP, W, C, G, α, setting an iteration counter initial value G to 1 and setting an iteration counter initial value G to q_l＝0，n_l0,1, 2, 3; initializing a solution set of NP initial solutions, also called candidate solutions, denoted P^g＝{P^g,1,P^g,2,···,P^g,NPAt each initial solution }It represents one of the scheduling schemes that is,representing the output active power of the kth generator in the ith scheduling scheme, initiallyIn thatIs randomly selected, i 1,2,., NP, K1, 2., K;

step 3.2, the candidate solution may not satisfy the power balance constraint, and the solution which does not satisfy the constraint condition is corrected by adopting a heuristic method, wherein the correction step is as follows: for candidate solutionsi 1, 2.. NP, the network loss is first calculated using the B-coefficient methodAnd calculating the candidate solution constraint violation according to:

if it is notIf not, randomly selecting one generator K from K generators to output power of the generator K, wherein the generator K belongs to {1,2Increase in current valueThen according to the generator capacity constraint pairBoundary constraint processing is carried out, namely the boundary constraint processing is set as a value rangeA boundary value or a random value of (1); recalculating network losses on this basisIf the corrected candidate solution still violates the power balance constraint, the correction process is repeated until a certain number of repairs is reached orEnding when the value is less than the minimum value of the constraint violation;

step 3.3, calculate candidate P first^g,iCost of each generatorAnd discharge of pollutantsK1, 2,. K; then calculates candidate P^g,iTotal cost ofAnd total amount of pollutant emissionsi＝1,2,...,NP；

Step 3.4, for solution set P^gThe ith candidate solution P in^g,iSelecting an operator r from the 3 mutation operators of the differential evolution algorithm, wherein the operator r belongs to {1,2,3} pair P^g,iPerforming mutation operation to generate variable V^g,iUpdating the using times n of the r-th operator in the latest W mutation operations at the same time_rAnd 3 differential operators are respectively described as:

V^g,i＝P^g,r1+α·(P^g,r1-P^g,r1)，

V^g,i＝P^g,r1+α·(P^g,r2-P^g,r3)+α·(P^g,r4-P^g,r5)，

V^g,i＝P^g,i+rand(0,1)·(P^g,r1-P^g,i)+α·(P^g,r2-P^g,r3)，

where α is the parameter for performing the corresponding operation, P^g,r1、P^g,r2、P^g,r3、P^g,r4And P^g,r5Are candidate solutions participating in mutation operations in differential evolution operators, and are selected from a solution set P^gIs randomly selected but different from each other, and rand (0,1) is a random number between 0 and 1; the selection method of the operator r comprises the following steps: if there are unused operators in the 3 operators, then randomly selecting one of the unused operators, and if all 3 operators are used, then selecting according to the following formula:

wherein n is_lThe number of times of using the l-th operator in the last W mutation operations, q_lC is a constant for the cumulative performance of operator l in the last W mutation operations. If P^g,iAnd (3) performing operation by using the mutation operator of the 1 st type or the mutation operator of the 2 nd type, and further performing the following operation: first generating a random integer k_randE {1, 2.., K }, and then for P^g,iEach variable of (1)Performing the crossover operation described by the following equation to generate a new candidate solution U^g,i：

Wherein,andeach represents a variable V^g,iAnd U^g,iβ is the parameter required for crossover operation if P^g,iOperating by using the 3 rd operator, the new candidate solution U is obtained^g,iValue of V^g,i；

Step 3.5, if candidate solution U^g,iIf the generator capacity constraint and the power balance constraint are not satisfied, the method described in the step 3.2 is adopted for correction, and a corrected candidate solution U is calculated^g,iCorresponding cost value and pollutant discharge amount; then put U^g,iAdding to solution set P^gIn this case P^gThe number of inner candidate solutions is NP + 1;

step 3.6, according to the new solution set P^gDeleting 1 worst solution to maintain the population number at NP; the deleting step is as follows: performing non-dominant ranking according to the cost values and pollutant discharge amount of all candidate solutions, and performing P^gIs divided into₁,F₂,...,F_LWaiting for the front edge; if the number of leading edge layers L is greater than 1 and F_LIf only 1 solution exists, the solution is directly deleted; if the number of leading edge layers F_LGreater than 1 and F_LIf there are multiple solutions, F is deleted_LThe solution with the most number of times dominated by other candidate solutions; if the number L of the leading edge layers is equal to 1, deleting the solution with the minimum value of the over-volume contribution of each candidate solution, wherein the over-volume contribution of each candidate solution is defined as the area of the area surrounded by all the candidate solutions and the reference point minus the area of the area surrounded by the remaining candidate solutions and the reference point after the candidate solutions are removed, and each target value of the reference point is the maximum value of the target value in all the candidate solutions;

step 3.7, calculating the operator effect gamma after the r mutation operator is applied in step 3.4_rAnd updating the accumulated performance q of each operator in the last W mutation operations_l，l＝1,2,3，γ_rCalculated using the formula:

γ_r＝1/m₁·1/(m₂+1)+m₃，

wherein m is₁For the solution candidate U generated in step 3.4 using the r-th mutation operator^g,iLeading edge sequence number, m, of non-bad sort in step 3.6₂For the solution candidate U generated in step 3.4 using the r-th mutation operator^g,iThe number of times dominated by other candidate solutions; m is₃Is a candidate solution U^g,iThe amount of the over-volume contribution of (c); cumulative performance q_lThe calculation method comprises the following steps: for all the recent W mutation operations, the W effect values are sorted from small to large, and the W-th effect is sorted as R_wThen, the delay ranking value DR of the w-th effect is calculated according to the following formula_wNamely:

DR_w＝(W-R_w+1)·D^W-w，

DR calculated after application of W-times operator_wAnd (4) further calculating the area under the AUC curve of each operator, and calculating the accumulated performance q of each operator according to the formula_l，l＝1,2,3：

Wherein, AUC_lThe area of the region enclosed by the ROC curve of the ith operator;

step 3.8, repeating the step 3.4 to the step 3.7 until P^gPerforming mutation operation on each candidate solution;

step 3.9, adding 1 to the G value of the iteration number counter, and repeating the step 3.4 to the step 3.8 until the G value of the iteration number counter reaches the maximum iteration number G; then to the solution set P^gPerforming a non-inferior ranking, P^gAll non-dominant solutions in (a) constitute a pareto solution setN is a solution setThe number of middle solutions.

Compared with the prior art, the invention has the following beneficial effects:

the method is simple to implement, high in operation speed and strong in robustness, realizes power balance constraint in the searching process by adopting an easy-to-implement heuristic mode, considers network loss in power balance, and is calculated and solved by a B coefficient method, and realizes generator capacity constraint in the searching process by a simple boundary constraint method;

in the multi-objective differential evolution algorithm, the convergence of a pareto solution set is ensured by adopting the measures of cooperative search of a plurality of differential evolution operators, inferior sequencing, domination frequency statistics and the like, wherein the differential evolution operator selection is realized according to the use frequency and the accumulated performance of each operator in a plurality of last mutation operations; ensuring the uniformity of the distribution of the pareto solution set by adopting an over-volume contribution strategy, and providing an optimal candidate scheme for subsequent decision;

and thirdly, on the basis of obtaining the pareto solution set with the optimal convergence and distribution uniformity, providing a scheduling method with the highest satisfaction degree for a decision maker by using a fuzzy set theory, and realizing the comprehensive optimal effect of the power generation cost and pollutant discharge amount of the thermal power plant.

In a word, the capacity constraint and the power balance constraint of the thermal generator set are realized by adopting a simple constraint processing method; the convergence and the distribution uniformity of the pareto solution sets in the searching process are ensured by adopting multi-differential operator collaborative searching and strategies based on super-volume contribution and the like, and the defect that the distribution uniformity of the solution sets is difficult to accurately keep by using a crowding degree method in a general multi-objective algorithm is overcome; the method adopts fuzzy set theory to select the scheme with the most satisfaction degree from the pareto solution set as the implementation scheme, and has strong practicability.

Drawings

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

FIG. 2 is a graph of a pareto frontier profile and selected trade-off obtained by an embodiment of the present invention;

FIG. 3 is a compromise between pareto front profile and selection obtained by the GDE3 method;

FIG. 4 is a graph of the convergence of the optimal cost target value of the present invention with the GDE3 method;

fig. 5 is a convergence curve of the optimal pollutant emission target value of the GDE3 method according to the present invention.

Detailed Description

The invention will be further described in detail with reference to examples of embodiments shown in the drawings to which, however, the invention is not restricted.

A thermal power plant environment economic dispatching method based on a multi-objective differential evolution algorithm is characterized by firstly establishing a thermal power plant environment economic dispatching mathematical model and acquiring model parameters of a specific thermal power plant generator set; then solving the model by utilizing a multi-target differential evolution algorithm to obtain a pareto frontier solution set; and calculating the satisfaction degree of each candidate scheduling scheme in the pareto frontier solution set on the basis, and selecting the scheduling scheme with the maximum satisfaction degree as the final scheduling scheme. The flow chart of the invention is shown in fig. 1, and specifically comprises the following steps:

step 1, establishing a mathematical model for environmental economic dispatching of a thermal power plant

The mathematical model of the economic dispatch of the thermal power plant comprises an objective function and a constraint condition, wherein the objective function is described by the following formula:

wherein K is the number of generators in the thermal generator set, P_kFor the kth generator to output active power, F_k(P_k) Is a k-th stationThe output active power of the generator is P_kThe required electricity generation cost of time, E_k(P_k) The active power output for the kth generator is P_kThe discharge amount of pollutants, K ═ 1, 2., K, the power generation cost F_k(P_k) The calculation expression of (a) is:

wherein, a_k、b_kAnd c_kThe cost coefficient of the kth generator, K is 1,2, and K, describes the cost of fuel, labor, equipment maintenance and the like required by the generator to generate electricity, and is a known parameter for a specific generator; pollutant emission E_k(P_k) The calculation expression of (a) is:

wherein, α_k、β_k、γ_k、ξ_kAnd λ_kThe characteristic parameters of pollutant emission of the kth generator are respectively, and K is 1, 2.

The constraint conditions of the environmental economic dispatching of the thermal power plant comprise generator capacity constraint and power balance constraint; the generator capacity constraint is described by:

wherein,andrespectively generate power for the kth stationMinimum and maximum machine output active power, K1, 2, K, which is a known parameter for a particular generator; the power balance constraint is that the total output power of the generator set is equal to the sum of network loss and load power, and the expression is as follows:

wherein,the total output active power P of K generators in the thermal generator set_DFor system load power, P_LPower is lost for the network; p_DIs a known parameter, P_LCan be calculated by the B coefficient method as follows:

wherein, B_ij、B_0iAnd B₀₀For the network loss factor, i 1,2, a, K, j 1,2, K, a parameter known for the particular power system;

for some thermal power plant environment economic scheduling problems, if the network loss constraint P is not considered_LThen the power balance constraint need only be satisfiedWithout calculating the loss by using B coefficient method, P can be made at this time_LIs 0. The mathematical model is the basis for multi-target scheduling in the subsequent steps;

step 2, obtaining various parameters in the model, including the total load power P of the system in the power grid_DCalculating the parameter B of the network power loss_ij、B_0iAnd B₀₀Minimum and maximum output active power of each generator in thermal generator setAndcost parameter a of generator set_k、b_kAnd c_iAnd pollutant discharge amount parameter α of generator set_k、β_k、γ_k、ξ_kAnd λ_k，i＝1,2,...,K，j＝1,2,...,K，k＝1,2,...,K。

The various model parameters obtained in the step depend on specific thermal power plant generator sets, and for some thermal power plant generator sets, if the index items ξ in the scheduling objective function values are not considered in calculation of the scheduling objective function values_kexp(λ_kP_k) Then ξ can be ordered_kAnd λ_kIs 0, j 1,2, a., K1, 2, a., K; the step of using the obtained parameters to instantiate a model of the specific thermal power plant;

and step 3: optimizing and solving the mathematical model by adopting a multi-target differential evolution algorithm to obtain a pareto solution setN is the number of solution concentration scheduling schemes, and solution concentration P^g,iRepresents the ith scheduling scheme, i 1, 2.., N;

this step results in a pareto solution set for the economic environment scheduling of a particular thermal power plant, which solution set represents several superior solutions, but since there is no direct goodness relationship between the solutions, i.e. for both solutions in the solution set, say a and B, solution a is superior to solution B on one objective function, but inferior to solution B on the other objective function; since the number of solutions in the pareto solution set is limited, and in order to obtain more excellent solutions to ensure the performance of the final selection of the compromise solution, the obtained pareto solution set should have good convergence and distribution uniformity, that is, the pareto solution set should approach the real pareto frontier of the problem as much as possible, and all solutions should be distributed uniformly as much as possible at the pareto frontier;

and 4, step 4: on the basis of the pareto solution set obtained in the step 3, the satisfaction degree of each scheduling scheme is calculated by using a fuzzy set theory method, and the calculation expression is as follows:

wherein, mu_nRepresenting the overall satisfaction of the nth scheduling scheme in the pareto solution set,representing the satisfaction of the cost objective in the nth scheduling scheme,represents the satisfaction of the pollutant discharge amount in the nth scheduling scheme, wherein N is 1, 2.Andcan be calculated from the following equations:

wherein, F_nAnd E_nRespectively representing the power generation cost and the pollutant discharge amount of the nth scheduling scheme, F_minAnd F_maxRespectively representing the minimum and maximum values of the electricity generation cost in all the N scheduling schemes, E_minAnd E_maxRespectively representing the minimum value and the maximum value of pollutant discharge amount in all N scheduling schemes;

in the pareto solution set obtained in the step 3, the solutions cannot be directly compared, and the overall quality degree of each solution is quantitatively described by adopting satisfaction in the step, so that the method is a basis for obtaining a compromise solution in the next step;

and 5: on the basis of obtaining the satisfaction degrees of all the scheduling schemes in the step 4, selecting the scheduling scheme with the maximum satisfaction degree as a final scheduling scheme;

this step selects the solution with the greatest satisfaction from the pareto solution set as the final scheduling scheme, which represents an optimal compromise between cost objectives and pollutant emissions, based on the satisfaction values of the solutions in the pareto solution set.

Further, the specific process of performing optimization solution on the mathematical model by using the multi-objective differential evolution algorithm in the step 3 is as follows:

step 3.1, initializing parameters NP, W, C, G, α, setting an iteration counter initial value G to 1 and setting an iteration counter initial value G to q_l＝0，n_l0,1, 2, 3; initializing a solution set of NP initial solutions, also called candidate solutions, denoted P^g＝{P^g,1,P^g,2,···,P^g,NPAt each initial solution }It represents one of the scheduling schemes that is,representing the output active power of the kth generator in the ith scheduling scheme, initiallyIn thatIs randomly selected, i 1,2,., NP, K1, 2., K;

step 3.2, the candidate solution may not satisfy the power balance constraint, and for the condition that the constraint strip is not satisfiedThe solution of the element is corrected by adopting a heuristic method, and the correction steps are as follows: for candidate solutionsi 1, 2.. NP, the network loss is first calculated using the B-coefficient methodAnd calculating the candidate solution constraint violation according to:

if it is notIf not, randomly selecting one generator K from K generators to output power of the generator K, wherein the generator K belongs to {1,2Increase in current valueThen according to the generator capacity constraint pairBoundary constraint processing is carried out, namely the boundary constraint processing is set as a value rangeA boundary value or a random value of (1); recalculating network losses on this basisIf the corrected candidate solution still violates the power balance constraint, the correction process is repeated until a certain number of repairs is reached orValue less than aboutEnding at a minimum value of the amount of the violation;

step 3.3, calculate candidate P first^g,iCost of each generatorAnd discharge of pollutantsK1, 2,. K; then calculates candidate P^g,iTotal cost ofAnd total amount of pollutant emissionsi＝1,2,...,NP；

Step 3.4, for solution set P^gThe ith candidate solution P in^g,iSelecting an operator r from the 3 mutation operators of the differential evolution algorithm, wherein the operator r belongs to {1,2,3} pair P^g,iPerforming mutation operation to generate variable V^g,iUpdating the using times n of the r-th operator in the latest W mutation operations at the same time_r(ii) a The 3 differential operators are respectively described as:

V^g,i＝P^g,r1+α·(P^g,r1-P^g,r1)，

V^g,i＝P^g,r1+α·(P^g,r2-P^g,r3)+α·(P^g,r4-P^g,r5)，

V^g,i＝P^g,i+rand(0,1)·(P^g,r1-P^g,i)+α·(P^g,r2-P^g,r3)，

where α is the parameter for performing the corresponding operation, P^g,r1、P^g,r2、P^g,r3、P^g,r4And P^g,r5Are candidate solutions participating in mutation operations in differential evolution operators, and are selected from a solution set P^gAre randomly selected but different from each other, and rand (0,1) is random between 0 and 1Counting; the selection method of the operator r comprises the following steps: if there are unused operators in the 3 operators, then randomly selecting one of the unused operators, and if all 3 operators are used, then selecting according to the following formula:

wherein n is_lThe number of times of using the l-th operator in the last W mutation operations, q_lC is a constant for the cumulative performance of operator l in the last W mutation operations. If P^g,iAnd (3) performing operation by using the mutation operator of the 1 st type or the mutation operator of the 2 nd type, and further performing the following operation: first generating a random integer k_randE {1, 2.., K }, and then for P^g,iEach variable of (1)Performing the crossover operation described by the following equation to generate a new candidate solution U^g,i：

Wherein,andeach represents a variable V^g,iAnd U^g,iβ is the parameter required for crossover operation if P^g,iOperating by using the 3 rd operator, the new candidate solution U is obtained^g,iValue of V^g,i；

Step 3.5, if candidate solution U^g,iIf the generator capacity constraint and the power balance constraint are not satisfied, the method described in the step 3.2 is adopted for correction, and a corrected candidate solution U is calculated^g,iCorrespond toCost value and pollutant emissions; then put U^g,iAdding to solution set P^gIn this case P^gThe number of inner candidate solutions is NP + 1;

step 3.6, according to the new solution set P^gThe 1 worst solution is deleted, so that the population number is maintained at NP. The deleting step is as follows: performing non-dominant ranking according to the cost values and pollutant discharge amount of all candidate solutions, and performing P^gIs divided into₁,F₂,...,F_LWaiting for the front edge; if the number of leading edge layers L is greater than 1 and F_LIf only 1 solution exists, the solution is directly deleted; if the number of leading edge layers F_LGreater than 1 and F_LIf there are multiple solutions, F is deleted_LThe solution with the most number of times dominated by other candidate solutions; if the number L of the leading edge layers is equal to 1, deleting the solution with the minimum value of the over-volume contribution of each candidate solution, wherein the over-volume contribution of each candidate solution is defined as the area of the area surrounded by all the candidate solutions and the reference point minus the area of the area surrounded by the remaining candidate solutions and the reference point after the candidate solutions are removed, and each target value of the reference point is the maximum value of the target value in all the candidate solutions;

step 3.7, calculating the operator effect gamma after the r mutation operator is applied in step 3.4_rAnd updating the accumulated performance q of each operator in the last W mutation operations_l，l＝1,2,3。γ_rCalculated using the formula:

γ_r＝1/m₁·1/(m₂+1)+m₃，

wherein m is₁For the solution candidate U generated in step 3.4 using the r-th mutation operator^g,iLeading edge sequence number, m, of non-bad sort in step 3.6₂For the solution candidate U generated in step 3.4 using the r-th mutation operator^g,iThe number of times dominated by other candidate solutions; m is₃Is a candidate solution U^g,iThe amount of the over-volume contribution of (c); cumulative performance q_lThe calculation method comprises the following steps: for all the recent W mutation operations, the W effect values are sorted from small to large, and the W-th effect is sorted as R_wThen, the w-th effect is calculated according to the following formulaCorresponding delayed sorting value DR_wNamely:

DR_w＝(W-R_w+1)·D^W-w，

DR calculated after application of W-times operator_wAnd (4) value and drawing an AUC curve of each operator by using an area under curve method, and calculating the accumulated performance q of each operator according to the following formula on the basis_l，l＝1,2,3：

Wherein, AUC_lThe area of the region enclosed by the AUC curve of the ith operator;

step 3.8, repeating step 3.4-step 3.7 until P^gPerforming mutation operation on each candidate solution;

step 3.9, adding 1 to the G value of the iteration number counter, and repeating the steps 3.4-3.8 until the G value of the iteration number counter reaches the maximum iteration number G; then to the solution set P^gPerforming a non-inferior ranking, P^gAll non-dominant solutions in (a) constitute a pareto solution setN is a solution setThe number of middle solutions.

Examples

As shown in fig. 1, by applying the method provided by the present invention, an IEEE30 node standard power system is selected as an experimental object to perform an environmental economic dispatch experiment, which includes the following specific steps:

1) establishing a mathematical model of the environmental economic dispatching of the thermal power plant, wherein aiming at an IEEE30 node system, the mathematical model can be integrally described as follows:

F_k(P_k)＝a_k+b_kP_k+c_kP_k ²，

E_k(P_k)＝α_k+β_kP_k+γ_kP_k ²+ξ_kexp(λ_kP_k)，

the generator capacity constraint and the power balance constraint should be satisfied simultaneously, namely:

wherein:

in the model, K is the number of the thermal generators in the generator set, P_kFor the kth thermal generator to output active power, F_k(P_k) The active power output for the kth thermal power generator is P_kThe required cost of electricity generation a_k、b_kAnd c_iCost factor of the ith thermal generator, E_k(P_k) The output active power of the ith thermal generator is P_kDischarge amount of pollutants of hour, α_k、β_k、γ_k、ξ_kAnd λ_kRespectively are pollutant discharge characteristic parameters of a kth thermal power generator,andrespectively outputting the minimum value and the maximum value, P, of the active power of the kth generator_DFor system load power, P_LFor network power loss, B_ij、B_0iAnd B₀₀A network loss factor, i 1,2, a., K, j 1,2, a., K1, 2, a., K;

2) obtaining various parameters in the model, wherein the number K of generators in an IEEE30 node standard system is 6, and the total load P of the system_D283.4MW, the parameters in the model and the network loss coefficients are shown in Table 1 and Table 2, respectively;

TABLE 1 IEEE30 node Standard System Generator set data

TABLE 2 IEEE30 node Standard System network loss coefficients

3) Initializing parameters NP-50, W-20, C-3, G-200, α -0.5, β -0.9, setting G-1, q_l＝0，n_l0,1, 2, 3; initialization solution set P^g＝{P^g,1,P^g,2,···,P^g,50Each candidate solutionIt represents one of the scheduling schemes that is,in thatIn the step (2), the random selection is carried out,i＝1,2,...,50,k＝1,2,...,6；

4) the candidate solution may not satisfy the power balance constraint formula, and the solution which does not satisfy the constraint condition is corrected by adopting a heuristic method, wherein the correction step is as follows: for candidate solutions1,2, 50, first the network loss coefficient B is used_ij、B_0iAnd B₀₀Calculating network lossAnd calculating the candidate solution constraint violation according to:

if it is notIf the power is not zero, one generator k is randomly selected from 6 generators, wherein the generator k belongs to {1, 2.. multidot.6 }, and the output power of the generator k is outputIncrease in current valueThen according to the generator capacity constraint type pairBoundary constraint processing is carried out, namely the boundary constraint processing is set as a value rangeA boundary value or a random value of (1); recalculating network losses on this basisIf the revised candidate solution still violatesPower balance constraint, repeating the above correction process until reaching a certain repair number orEnding when the value is less than the minimum value of the constraint violation;

5) for each solution P in the candidate solution set^g,i1, 2.., 50, calculating the cost of each generatorAnd discharge of pollutants1,2, 6; then calculate the plan P^g,iTotal cost ofAnd total amount of pollutant emissions

6) For solution set P^gThe ith candidate solution P in^g,iSelecting an operator r e {1,2,3} pair P from the following 3 differential algorithm mutation operators^g,iPerforming an operation to generate a variable V^g,iAnd updating the using times n of the r-th operator in the last 20 mutation operations_r。

V^g,i＝P^g,r1+α·(P^g,r1-P^g,r1)，

V^g,i＝P^g,r1+α·(P^g,r2-P^g,r3)+α·(P^g,r4-P^g,r5)，

V^g,i＝P^g,i+rand(0,1)·(P^g,r1-P^g,i)+α·(P^g,r2-P^g,r3)，

Wherein α is 0.5, P^g,r1、P^g,r2、P^g,r3、P^g,r4And P^g,r5Are candidate solutions for the operation of the differential evolution operator, which are represented by a solution set P^gAre randomly selected but different from each other, and rand (0,1) is a random number between 0 and 1. The operator r selection method comprises the following steps: if there are unused operators in the 3 operators, randomly selecting one of the unused operators, and if all of the 3 operators are used, selecting according to the following formula:

wherein, C is 3, n_lThe number of uses of the l-th operator in the last 20 mutation operations, q_lThe accumulated performance of operator l in all mutation operations of the last 20 times. If P^g,iMutation by using the 1 st or 2 nd operator generates a random integer k_randE {1, 2.., 6}, and then for P^g,iEach variable of (1)Performing a crossover operation to generate a new candidate solution U^g,iNamely:

wherein, β is 0.9,andeach represents a variable V^g,iAnd U^g,iThe kth variable of (1); if P^g,iCarrying out mutation by adopting a 3 rd operator, and then carrying out new candidate solution U^g,iValue of V^g,i；

7) If candidate solution U^g,iIf the constraint formula is not satisfied, correcting according to the method described in the step 4), and calculating the corrected candidate solution U^g,iCorresponding cost value and pollutant discharge amount; then put U^g,iAdding to solution set P^gIn (i) P^gThe number of inner candidate solutions is 51;

8) according to the new solution set P^gThe 1 worst solution is deleted, and the population number is maintained at 50. The deleting step is as follows: performing non-inferior ranking according to the cost values and pollutant discharge amount of all candidate solutions, and performing P^gIs divided into₁,F₂,...,F_LWaiting for the front edge; if the number of leading edge layers L is greater than 1 and F_LIf only 1 solution exists, the solution is directly deleted; if the number of leading edge layers F_LGreater than 1 and F_LIf there are multiple solutions, F is deleted_LThe solution with the most number of times dominated by other candidate solutions; if the number L of the leading edge layers is equal to 1, deleting the solution with the minimum value of the over-volume contribution of each candidate solution, wherein the over-volume contribution of each candidate solution is defined as the area of the area surrounded by all the candidate solutions and the reference point minus the area of the area surrounded by the remaining candidate solutions and the reference point after the candidate solutions are removed, and each target value of the reference point is the maximum value of the target value in all the candidate solutions;

9) calculating the operator effect gamma after the r mutation operator is applied in the step 6)_rAnd updating the cumulative performance of each operator in the last 20 mutation operations, namely q_lAnd l is 1,2, 3. Effect value gamma_rCalculated using the formula:

γ_r＝1/m₁·1/(m₂+1)+m₃，

wherein m is₁For step 6), using the r-th mutation operator to generate a candidate solution U^g,iLeading edge sequence number m of non-inferior sorting in step 8)₂The solution candidate U generated by the r mutation operator in the step 6)^g,iThe number of times dominated by other candidate solutions; m is₃Is a candidate solution U^g,iThe over-volume contribution of (c); cumulative performance q_lThe calculation method comprises the following steps: for all mutation operations of the last 20 times, 20 effect values are sorted from small to large, and the w-th effect is sorted as R_wThen the w-th effect is calculated according to the following formulaDelayed sorting value DR_wNamely:

DR_w＝(20-R_w+1)·D^20-w，

DR calculated after 20 times of operator application_wThe value and area under the curve method draws the AUC curve of each operator, and on the basis, the accumulated performance q of each operator is calculated according to the following formula_l，l＝1,2,3：

Wherein, AUC_lThe area of the region enclosed by the AUC curve of the ith operator;

10) repeating the steps 6) to 9) until P^gPerforming mutation operation on each candidate solution;

11) adding 1 to the value of the iteration counter g, and repeating the steps 6) to 10) until g is 200; then to the solution set P^gPerforming a non-inferior ranking, P^gAll non-dominant solutions in (a) constitute a pareto solution setN is a solution setThe number of middle solutions.

12) Computing by fuzzy set theory methodSatisfaction degree mu of each candidate scheduling scheme_nFor the nth scheme, μ_nThe calculation method comprises the following steps:

wherein, F_nAnd E_nRespectively representing the power generation cost and the pollutant discharge amount of the nth scheduling scheme, F_minAnd F_maxRespectively representing the minimum and maximum values of the electricity generation costs in all the scheduling schemes, E_minAnd E_maxRespectively representing the minimum and maximum values of pollutant emissions in all scheduling schemes.

13) μ from all scheduling schemes calculated in step 12)_nAnd selecting the scheduling scheme with the maximum satisfaction as the final scheduling scheme.

Computer simulations were performed on this example and compared to the well-known multi-objective differential evolution algorithm GDE 3. Table 3, table 4 and table 5 show the optimal cost scheduling result, optimal pollutant emission scheduling result and compromise scheduling result of the method of the present invention and the GDE3 method, respectively. Fig. 2 and 3 show the pareto frontier distribution and selection trade-off obtained by the method of the present invention and the GDE3 method, respectively, and fig. 4 and 5 show the optimal cost scheduling and optimal emission scheduling target value convergence curves of the method of the present invention and the GDE3 method, respectively. From the pareto frontier profile obtained in the method of the present invention, as shown in fig. 2, the corresponding compromise point location is found.

Aiming at the optimal cost scheduling, the results in the table 3 show that the cost of the scheduling scheme is 606.0130$/h, which is superior to 606.0323$/h of the scheduling scheme of the GDE3 method; aiming at the optimal emission scheduling, the results in Table 4 show that the pollutant emission amount of the scheduling scheme is 0.1942t/h, which is slightly superior to 0.1943t/h of the scheduling scheme of the GDE3 method; for the compromise solution, the results in Table 3 show that the scheduling scheme of the invention has the cost of 614.5114$/h and the pollutant discharge amount of 0.2016t/h, and is superior to the scheduling scheme of the GDE3 method. Meanwhile, fig. 2 and 3 show that the pareto fronts obtained by the method of the present invention are more uniformly distributed than the pareto fronts obtained by the GDE3 method; the convergence curves of fig. 4 and 5 also show that the method of the present invention also has better convergence than the GDE3 method. The results fully show that the method has the characteristics of high precision, uniform pareto frontier solution set distribution and high convergence rate in the problem of environmental economic dispatching of the thermal power plant.

TABLE 3 optimal cost scheduling results for the method of the present invention and GDE3 method

TABLE 4 optimal emissions scheduling results for the method of the present invention and GDE3 method

TABLE 5 compromise of scheduling results for the method of the present invention and the GDE3 method

The foregoing is only a preferred embodiment of the present invention, and it should be noted that the present invention is not limited to the above-mentioned embodiment, and further modifications can be made without departing from the principle of the present invention, and these modifications should also be construed as the protection scope of the present invention.

Claims (1)

1. A thermal power plant environment economic dispatching method based on a multi-objective differential evolution algorithm is characterized by comprising the following steps:

step 1: establishing mathematical model for environmental economic dispatch of thermal power plant

The mathematical model of the economic dispatch of the thermal power plant comprises an objective function and a constraint condition, wherein the objective function is described by the following formula:

wherein K is the number of generators in the thermal generator set, P_kFor the kth generator to output active power, F_k(P_k) The active power output for the kth generator is P_kThe required electricity generation cost of time, E_k(P_k) The output active power of the ith generator is P_kThe discharge amount of pollutants, K ═ 1, 2., K, the power generation cost F_k(P_k) The calculation expression of (a) is:

F_k(P_k)＝a_k+b_kP_k+c_kP_k ²，

wherein, a_k、b_kAnd c_kA cost factor for the kth generator, K being 1,2, a, K, describing the cost of fuel, labor and equipment maintenance required to generate the generator, is a known parameter for the particular generator;

pollutant emission E_k(P_k) The calculation expression of (a) is:

E_k(P_k)＝α_k+β_kP_k+γ_kP_k ²+ξ_k exp(λ_kP_k)，

wherein, α_k、β_k、γ_k、ξ_kAnd λ_kThe characteristic parameters of pollutant emission of the kth generator are respectively, and K is 1,2, a.

The constraint conditions of the environmental economic dispatching of the thermal power plant comprise generator capacity constraint and power balance constraint; the generator capacity constraint is described by:

wherein,andrespectively outputting a minimum value and a maximum value of active power of a kth generator, wherein K is 1, 2. The power balance constraint is that the total output power of the generator set is equal to the sum of network loss and load power, and the expression is as follows:

wherein,the total output active power P of K generators in the thermal generator set_DFor system load power, P_LPower is lost for the network; p_DIs a known parameter, P_LCalculated by the B coefficient method as follows:

wherein, B_ij、B_0iAnd B₀₀For the network loss factor, i 1,2, a, K, j 1,2, K, a parameter known for the particular power system;

step 2: obtaining various parameters required in the model in the step 1, including: total load power P of system in power grid_DCalculating the parameter B of the network power loss_ij、B_0iAnd B₀₀Minimum and maximum output active power of each generator in thermal generator setAndcost parameter a of generator set_k、b_kAnd c_iAnd pollutant discharge amount parameter α of generator set_k、β_k、γ_k、ξ_kAnd λ_k，i＝1,2,...,K，j＝1,2,...,K，k＝1,2,...,K；

And step 3: the mathematical model in the step 1 is optimized and solved by adopting a multi-objective differential evolution algorithm to obtain a pareto solution setN is the number of solution concentration scheduling schemes, and solution concentration P^g,iRepresents the ith scheduling scheme, i 1, 2.., N;

and 4, step 4: on the basis of the pareto solution set obtained in the step 3, the satisfaction degree of each scheduling scheme is calculated by using a fuzzy set theory method, and the calculation expression is as follows:

wherein, mu_nRepresenting the overall satisfaction of the nth scheduling scheme in the pareto solution set,representing the satisfaction of the cost objective in the nth scheduling scheme,represents the satisfaction degree of the pollutant discharge amount in the nth scheduling scheme, wherein N is 1, 2.Andare calculated by the following formulas, respectively:

wherein, F_nAnd E_nRespectively representing the power generation cost and the pollutant discharge amount of the nth scheduling scheme, F_minAnd F_maxRespectively representing the minimum and maximum values of the electricity generation cost in all the N scheduling schemes, E_minAnd E_maxRespectively representing the minimum value and the maximum value of pollutant discharge amount in all N scheduling schemes;

and 5: the satisfaction degree mu of the scheduling scheme obtained in the step 4_nComparing and sequencing the numerical values, selecting the scheduling scheme with the maximum satisfaction degree value as a final scheduling scheme, and corresponding the numerical value of the selected scheduling scheme to the pareto frontier distribution map drawn according to the method;

the various model parameters obtained in the step 2 depend on specific thermal power plant generator sets, and for some thermal power plant generator sets, if the index items ξ in the scheduling objective function values are not considered when calculating the scheduling objective function values_k exp(λ_kP_k) Then let ξ_kAnd λ_k0, j 1,2, a., K1, 2, a., K, which instantiates a model for a specific thermal power plant using the obtained parameters;

the step 3 comprises the following steps:

step 3.1, initializing parameters NP, W, C, G, α, setting an iteration counter initial value G to 1 and setting an iteration counter initial value G to q_l＝0，n_l0,1, 2, 3; initializing a solution set of NP initial solutions, also called candidate solutions, denoted P^g＝{P^g,1,P^g,2,···,P^g,NPAt each initial solution }It represents one of the scheduling schemes that is,representing the output active power of the kth generator in the ith scheduling scheme, initiallyIn thatIs randomly selected, i 1,2,., NP, K1, 2., K;

step 3.2, the candidate solution may not satisfy the power balance constraint, and the solution which does not satisfy the constraint condition is corrected by adopting a heuristic method, wherein the correction step is as follows: for candidate solutionsFirstly, the network loss is calculated by using a B coefficient methodAnd calculating the candidate solution constraint violation according to:

in thatWhen the output power is not zero, randomly selecting one generator K from K generators to form a generator K belonging to {1, 2.. multidot.K }, and outputting the output powerIncrease in current valueThen according to the generator capacity constraint pairBoundary constraint processing is carried out, namely the boundary constraint processing is set as a value rangeA boundary value or a random value of (1); recalculating network losses on this basis

If the corrected candidate solution still violates the power balance constraint, the correction process is repeated until a certain number of repairs is reached orEnding when the value is less than the minimum value of the constraint violation;

step 3.3, calculate candidate P first^g,iCost of each generatorAnd discharge of pollutants Then calculates candidate P^g,iTotal cost ofAnd total amount of pollutant emissions

Step 3.4, for solution set P^gThe ith candidate solution P in^g,iSelecting an operator r from the 3 mutation operators of the differential evolution algorithm, wherein the operator r belongs to {1,2,3} pair P^g,iPerforming mutation operation to generate variable V^g,iUpdating the using times n of the r-th operator in the latest W mutation operations at the same time_rAnd 3 differential operators are respectively described as:

V^g,i＝P^g,r1+α·(P^g,r1-P^g,r1)，

V^g,i＝P^g,r1+α·(P^g,r2-P^g,r3)+α·(P^g,r4-P^g,r5)，

V^g,i＝P^g,i+rand(0,1)·(P^g,r1-P^g,i)+α·(P^g,r2-P^g,r3)，

where α is the parameter for performing the corresponding operation, P^g,r1、P^g,r2、P^g,r3、P^g,r4And P^g,r5Are candidate solutions participating in mutation operations in differential evolution operators, and are selected from a solution set P^gIs randomly selected, but is different from each other, and rand (0,1) is a random number between 0 and 1; the selection method of the operator r comprises the following steps: if there are unused operators in the 3 operators, then randomly selecting one of the unused operators, and if all 3 operators are used, then selecting according to the following formula:

wherein n is_lThe number of times of using the l-th operator in the last W mutation operations, q_lThe accumulated performance of the operator l in the last W mutation operations is shown, and C is a constant; if P^g,iAnd (3) performing operation by using the mutation operator of the 1 st type or the mutation operator of the 2 nd type, and further performing the following operation: first generating a random integer k_randE {1, 2.., K }, and then for P^g,iEach variable of (1)Performing the crossover operation described by the following equation to generate a new candidate solution U^g,i：

Wherein,andeach represents a variable V^g,iAnd U^g,iβ is the parameter required for crossover operation if P^g,iOperating by using the 3 rd operator, the new candidate solution U is obtained^g,iValue of V^g,i；

Step 3.5, if candidate solution U^g,iIf the generator capacity constraint and the power balance constraint are not satisfied, the method described in the step 3.2 is adopted for correction, and a corrected candidate solution U is calculated^g,iCorresponding cost value and pollutant discharge amount; then put U^g,iAdding to solution set P^gIn this case P^gThe number of inner candidate solutions is NP + 1;

step 3.6, according to the new solution set P^gDeleting 1 worst solution to maintain the population number at NP; the deleting step is as follows: performing non-dominant ranking according to the cost values and pollutant discharge amount of all candidate solutions, and performing P^gIs divided into₁,F₂,...,F_LA leading edge; if the number of leading edge layers L is greater than 1 and F_LIf only 1 solution exists, the solution is directly deleted; if the number of leading edge layers F_LGreater than 1 and F_LIf there are multiple solutions, F is deleted_LThe solution with the most number of times dominated by other candidate solutions; if the number L of the leading edge layers is equal to 1, deleting the solution with the minimum value of the over-volume contribution of each candidate solution, wherein the over-volume contribution of each candidate solution is defined as the area of the area surrounded by all the candidate solutions and the reference point minus the area of the area surrounded by the remaining candidate solutions and the reference point after the candidate solutions are removed, and each target value of the reference point is the maximum value of the target value in all the candidate solutions;

step 3.7, calculating the operator effect gamma after the r mutation operator is applied in step 3.4_rAnd updating the accumulated performance q of each operator in the last W mutation operations_l，l＝1,2,3，γ_rCalculated using the formula:

γ_r＝1/m₁·1/(m₂+1)+m₃，

wherein m is₁For the solution candidate U generated in step 3.4 using the r-th mutation operator^g,iLeading sequence of non-bad sequence in step 3.6Number m₂For the solution candidate U generated in step 3.4 using the r-th mutation operator^g,iThe number of times dominated by other candidate solutions; m is₃Is a candidate solution U^g,iThe amount of the over-volume contribution of (c); cumulative performance q_lThe calculation method comprises the following steps: for all the recent W mutation operations, the W effect values are sorted from small to large, and the W-th effect is sorted as R_wThen, the delay ranking value DR of the w-th effect is calculated according to the following formula_wNamely:

DR_w＝(W-R_w+1)·D^W-w，

DR calculated after application of W-times operator_wAnd (4) further calculating the area under the AUC curve of each operator, and calculating the accumulated performance q of each operator according to the formula_l，l＝1,2,3：

Wherein, AUC_lThe area of the region enclosed by the AUC curve of the ith operator;

step 3.8, repeating the step 3.4 to the step 3.7 until P^gPerforming mutation operation on each candidate solution;

step 3.9, adding 1 to the G value of the iteration number counter, and repeating the step 3.4 to the step 3.8 until the G value of the iteration number counter reaches the maximum iteration number G; then to the solution set P^gPerforming a non-inferior ranking, P^gAll non-dominant solutions in (a) constitute a pareto solution setN is a solution setThe number of middle solutions.