Abstract
Dynamic economic and environmental load dispatch (DEED) aims to determine the amount of electricity generated from power plants during the planning period to meet load demand while minimizing energy consumption costs and environmental pollution emission indicators subject to the operation requirements. This planning problem is usually expressed using a nonsmooth cost function, taking into account various equality and inequality constraints such as valve-point effects, operational limits, power balance, and ramp rate limits. This paper presents DEED models developed for a system consisting of thermal units, wind power generators, photovoltaic (PV) generators, and energy storage (ES). A selection hyper-heuristic algorithm is proposed to solve the problems. Three heuristic mutation operators formed a low-level operator pool to direct search the solution space of DEED. The high level of SHHA evaluates the performances of the low-level operators and dynamically adjusts the chosen probability of each operator. Simulation experiments were carried out on four systems of different types or sizes. The results verified the effectiveness of the proposed method.
1. Introduction
With the rapid growth of modern industry and global economy, people’s demand for electricity in production and life is constantly increasing. Therefore, how to dispatch power production according to the forecasted load demand and improve the generation efficiency of generating units has always been an important research issue in power system operational planning. Traditional economic load dispatch aims to determine the operation status of each unit in a power plant within the planned time and minimize the energy consumption cost while meeting the demand for power load and the operational requirements of generating units. With the gradual attention of human beings to the ecological environment, the harmful substances such as sulfur oxides, nitrogen oxides, and carbon dioxide emitted from fossil fuels to the atmosphere in the combustion process have gradually attracted attention. While reducing energy consumption costs, further reducing environmental pollution emission indicators has become a new requirement for optimizing load dispatch.
For dispatching tasks, energy consumption cost and environmental pollution emissions are two key indicators and minor improvements can bring significant economic and environmental benefits [1]. The optimization of these two indicators is essentially contradictory. Economic and emission load dispatch (EED) is a multiobjective optimization problem that considers these two conflicting indicators at the same time [2]. In addition, because the mathematical expressions of the two indexes are nonconvex and nonsmooth objective functions, it is difficult to solve the optimization problem. The energy crisis caused by the long-term consumption of fossil energy and the global demand for clean energy promotes the development of renewable energy. Wind power, solar power, and plug-in electric vehicles have been widely developed in the last decade and might become the main alternatives in the power system in the near future [3, 4]. But at the same time, the introduction of these renewable energy sources will bring more complex planning and uncertainty to the power grid production [5].
Heuristic optimization algorithms are widely used to solve nonsmooth optimization problems with high nonlinearity and discontinuity of input and output, such as the economic and environmental dispatch problem, due to their strong global search ability and low requirements for some mathematical properties of the objective function. An improved bacterial foraging algorithm with a parameter automation strategy and crossover operation is proposed to improve computational efficiency and solve the static/dynamic EED [6]. In [7], an improved particle swarm optimization algorithm is applied to EED including solar photovoltaic generation. In [8], a bilevel dynamic economic emission dispatch model is proposed considering wind power integration. Also, robust optimization and a nested solution method are proposed for EED considering wind power uncertainty. A dispatch model combining environmental awareness and risk perception is proposed in [9], and nondominated sorting genetic algorithm-II embedded constraints handling techniques are devised for the problem. In [10], an enhanced multiobjective bee colony optimization algorithm is proposed for solving short-term hydro-thermal-wind complementary scheduling considering the uncertainty of wind power. In [2], a dynamic economic emission dispatch (DEED) incorporating renewable energy and energy storage is formulated and an adaptive robust method is proposed to solve the problem.
There are many differences in EED models that have been solved in the existing research. Some scholars focus on the dynamic economic emission dispatch (DEED) problems which consider only thermal generation, single-objective DEED with renewable energy generation (REG), and the multiobjective static DEED with REG. This paper formulates a DEED model including thermal generation, wind energy generation, solar photovoltaic (PV) generation, and energy storage (ES). A selection hyper-heuristic algorithm is proposed to solve the Pareto optimal solution of the problem. Three operators formed a low-level operator pool to directly search the solution space. The high-level structure is responsible for managing low-level operators. The operators are evaluated and ranked by analyzing the results of the nondominated solutions. Then, the chosen probability of each low-level operator is determined by the selection mechanism, which is provided for the next call. The effectiveness of the optimization strategy is validated through a series of simulations. The results are compared with those of existing state-of-the-art optimization techniques.
The rest of this paper is organized as follows. Section 2 presents the DEED mathematical model. The application of the proposed hyper-heuristic method is described in Section 3. To illustrate its effectiveness, the proposed method was tested on a series of systems in Section 4. Conclusions and future research are discussed in Section 5.
2. Problem Formulation
In this section, we will discuss the mathematical model of the DEED problem, which is used to distribute the output power of all generators economically and environmentally during the dispatching period and satisfies a series of equality and inequality constraints.
2.1. Objective Function
Compared with the static EED, which is scheduled at a certain hour, DEED adopts dynamic scheduling for a load cycle of 24 hours. Consider a power system with thermal power generators, wind power generators, and PV generators over time intervals. The basic objective function of the DEED problem can be formulated as the following optimal function:where and are the number of generator and time interval, respectively. denotes the total fuel costs over the dispatch period. , , and are cost functions of the thermal power generator, wind power generator, and solar power generator, respectively. is the time horizon. , , and are the power output of thermal power generator, wind power generator, and PV generator, respectively. , , , , and are environmental cost coefficients.
When the valve-point effect (VPE) is considered, the total fuel costs of thermal power generators can be written as follows:where is the minimum power outputs of generator . , , and are the fuel cost coefficients of generator . and are the coefficients of generator reflecting the valve-point effects. Valve-point effect is the phenomenon of rapid increase consumption caused by valve drawing effects during start-up of thermal power generators.
The total economic costs of wind power generators can be written as follows:where and are the schedule and actual power of wind generator at time interval . denotes the cost function of generator (payment to wind farm operators for actual use of wind power), which can be calculated by the cost coefficients and the wind probability density function. is the direct cost coefficient of generator . denotes the penalty cost function for not using all available power of generator , which has a linear correlation with the difference between available power and actual power. is the underestimating penalty cost coefficient. is the probability distribution function of wind power. denotes the reserve costs associated with wind power uncertainty (penalties for estimating available wind power). is the rated power of generator , and is the reserve cost coefficient.
The total economic costs of solar power generators can be written as follows:where and are the schedule and actual power of solar generator at time interval . denotes the weighted cost of generator based on solar irradiance, which can be calculated by the cost coefficients and the solar probability density function. and are the penalty cost function for not using all available power of generator and the reserve costs associated with power uncertainty, respectively. is the operating cost coefficient of solar power operator. is the rated power of generator , and is the upper limit of output of the corresponding operator. and are the underestimating penalty cost coefficient and the reserve cost coefficient, respectively. is the probability distribution function of solar power.
2.2. Constraints
(1)Real Power Balance Constrains In order to ensure the safe and stable operation of generators, each generator needs to run in a specified power range: where and are the minimum and maximum power outputs of thermal generator , respectively. and denote the charge and discharge power of the ES unit. is the power capacity of ES unit.(2)Power Operating Constrains In order to ensure the safe and stable operation of generators, each generator needs to run in a specified power range: where and are the minimum and maximum power outputs of thermal generator , respectively. and denote the charge and discharge power of the ES unit. is the power capacity of the ES unit.(3)Generator Ramp Rate Constrains The power output of the generators should not fluctuate too much during the interval time interval. Ramp rate constraints are used to restrict the adjusting range of the generator power: where and are the ramp-up and ramp-down rate constraints of generator , respectively.(4)Spinning Reserve Requirements Spinning reserve requirements (SRRs) refer to the readily available capacity, which can compensate for the generation gap caused by generator faults or power fluctuation in a short time. SRR constraints are frequently applied in the unit commitment problem and economic dispatch problem [11–13]. It can be employed to provide the operating reserves to balance the influence brought by the renewables uncertainty and forecast error [14]. where and are the spinning reserves for the one-hour and 10-minute compensation time.(5)Energy Storage SOC Constraints The state of charge (SOC) denotes an ES’s real-time power state, which is the variable that actually operates in the scheduling process. It is calculated according to the state at the previous time interval:where is the SOC of the ES at time interval which lies between 0 and 1. and are the minimum and maximum limits of SOC, respectively.
3. Multiobjective Selection Hyper-Heuristic Algorithm
In this paper, a selection hyper-heuristic algorithm is used to solve the Pareto optimal solution of the DEED problem. The method consists of two levels, and the low-level structure contains an algorithm pool, which is composed of three heuristic operators to directly search the solution space. In the management of the high-level structure, the low-level operators are evaluated and ranked by calculating the results of nondominated solutions. Then, the chosen probability of each low-level operator is calculated through the selection mechanism, which is provided for the next call.
3.1. The Framework of Selection Hyper-Heuristic Algorithm
Hyper-heuristic algorithms are widely used in the field of automatic algorithm design [15]. Hyper-heuristics can be defined as “automated methods for selecting or generating heuristics to solve hard computing search problems.” [16] This is different from most implementations of metaheuristic methods that operate directly on solution spaces. It was originally used in 1997 to describe a protocol that combines several AI methods in the context of automatic theorem proving. The proposed multiobjective selection hyper-heuristic algorithm is shown in Figure 1. The low-level operator pool consists of three heuristic operators, and the high level of the proposed algorithm consists of an evaluation strategy and a selection mechanism.
The low-level operator pool plays an important role in searching the problem space effectively. A mutation operator pool composed of three different mutation operators is built in this section. These operators will be implemented under the guidance of an adaptive selection mechanism at different stages of the search. The operators are as follows: H1: GA/derived from NSGA-II H2: DE/rand/2 H3: DE/current-to-best/1 derived from PSOwhere is the offspring and is the parent with being the upper bound on the parent component, is the lower bound, and is small variation which is calculated from a polynomial distribution [17].
Instead of operating directly on a solution space, the high level of the hyper-heuristics evaluates and grades the historical performance of the low-level operators during a given stage. A selection mechanism is devised to manage the calls of the low-level operators at different stages in the searching process. The underlying assumption is that each low-level operator has different performance in different regions of search space or at different stages of the whole searching process [15]. In fact, due to the function landscapes of the problems and the characteristics of the operators, it is difficult for the operator to keep stable optimization capability in different regions of the whole search area. It is common that some operators perform well in the early stage of iteration, but lack optimization capability in the later stage. Here, we introduce a sliding time window (STW) to dynamically evaluate the historical performance of the low-level operators. The length of the STW is much smaller than the maximum number of iteration. The uniformity of the nondominated solution and the dominated relationship between parent and offspring individuals are used as the evaluation indexes.
The specific steps of the evaluation and selection mechanism are as follows: Step 1 (evaluation): compute the spatial measure (uniformity) of nondominated solutions, . Step 2: calculate the number of the parent dominated by the offspring , and the number of the parent and the offspring does not dominate each other, . Step 3: calculate FIR score: Step 4: update the sliding window time, , and the historical information, . Step 5 (selection): collect the number of call frequency of each operator, . Step 6: determine the chosen probability based on historical information: Step 7: select the low-level operator using the roulette method.
3.2. Elitist Nondominated Sorting Mechanism
In different periods of iteration, the crowding distance between individuals is used as an index to guide the selection of parent individuals. In the early stage of iteration, the combination of nondominated sorting and crowding distance is used to select individuals. Firstly, the nondominated sorting of solutions is performed to obtain different Pareto frontiers of solutions. Then, the individuals are sorted in descending order of crowding distance in each Pareto frontier. Finally, 75% of the larger crowding distance in each Pareto frontier is selected as the paternal individuals until the total number of parent individuals is satisfied. At the later stage of the iteration, the population converges gradually and the density of solutions in the population is high. This paper combines nondominated sorting, crowding distance, and hypercube grid to select parent individuals. Individuals will be screened in two ways. One is nondominated sorting. The other is a combination of nondominated sorting, crowding distance, and hypercube grid. A parameter is introduced to decide which way to choose. First, individuals are screened using the same steps as in the early stage. Then, the selected individuals are put into the grid, and the grid with larger individual density is screened. In each selected grid, some individuals are randomly selected for deletion to guarantee the uniform distribution of individuals. By introducing crowding distance and hypercube grid, not only the diversity of the population is guaranteed but also the individuals are evenly distributed in the solution space. This helps to improve the diversity of the algorithm and maintain its search capability (Algorithm 1).
|
3.3. Constraint Handling
DEED involves a number of equality and inequality constraints. We cannot guarantee every generated solution will satisfy all constraints, especially real power balance and ramp rate constraints. Penalty-based methods are the most commonly used methods to deal with these complex constraints, but it has some shortcomings in the calculation efficiency. Furthermore, how to determine the punishment degree is also a challenge. In this section, we devised a heuristic process to deal with the constraints. Compared with [18], our method can avoid the situation that the output power of the adjacent hour violates the power balance constraints, especially at the later stage of the adjustment process. Handling of the SRRs in this chapter is referenced from Ref. [19]. Step 1: randomly select an hour . Step 2: calculate the corresponding feasible region of the adjustable output: Step 3: adjust the output of each generator at the current hour to satisfy (15) and (16) using (17): Step 4: tackle the SRR constraints: check if SRRs is satisfied and then go to Step 5; otherwise, calculate the violation of SRRs and set . Subtract from which are equal to their upper limits. Then, generate a new feasible individual at hour . Calculate and repeat the adjust process until the SRR constraints are satisfied at an hour. Step 5: tackle the demand-supply constraint by the following steps:(1)Calculate the system transmission loss using (5).(2)Calculate the difference between the power demand and the power grid output: Check if and then go to Step 2 with the next randomly selected hour; otherwise, go to the next step. Here, is a tolerance variable which decreases from early stages of the evolutionary process and is decreased gradually to a small value [20].(3)Generate a random sequence . Select as the first slack generator and adjust its power as
(4)Adjust using (17) to satisfy the feasible region limit. Update and select the next slack generator in turn according to . is a random value which lies between [0.2, 1]. Step 6: repeat the previous steps until a feasible solution is found.
3.4. Selection Hyper-Heuristic Algorithm for DEED
The detailed procedure of the proposed method for DEED is presented as follows: Step 1: initialize the setting parameters of the algorithm, such as population size, the maximum number of iteration, number of objective, and parameters of low-level operators. Step 2: generate the initial population, perform the nondominated sorting, and calculate the crowding distance. Step 3: select individuals to build up the parent population using the tournament selection, and then generate the initial offspring population by calling the low-level operators. Step 4: memorize the current iteration number and start the evolution process. Step 5: combine the parent population and the offspring population, conduct the nondominated sorting strategy, and update the crowding distance. Step 6: select the parent individuals using the proposed nondominated sorting strategy mentioned in the previous section, and then generate the offspring individuals. Step 7: grade and rank the low-level operators using the evaluation method presented in the previous section. Step 8: choose the low-level operator using the selection mechanism presented in the previous section. Step 9: memorize the optimal solution and check whether the current iteration is greater than the maximum iteration number. If not, continue the cycle; otherwise, output the global optima solution.
The flowchart of the proposed method is presented as follows.
4. Results and Discussion
For the experimental study, four test cases are studied to evaluate the performance of the proposed algorithm. The scheduling time of all cases is selected as one day with 24 hours, and the valve-point effect and the ramp rate constraints are taken into account. Case 1 and Case 2 are classic DEEDs with 5 thermal units and 10 thermal units without considering wind and PV generation, respectively. We use these two cases to verify the effectiveness of our proposed method by comparing them with the existing algorithms. Case 3 comes from an IEEE 57 bus system; three PV and two wind generators are incorporated into this system. Case 4 considers the ES unit on the basis of Case 3. We use this case to study the impact of ES on the system operation and its role in peak shaving and valley filling in the power supply.
Wind power output depends on wind speed, type of wind turbine, and parameters of its power performance curve. Wind speed is a random variable that needs to be predicted and has uncertainty. The previous study has proved that the wind speed at a given location closely follows a Weibull distribution over time [21]. The output power of wind turbines can be predicted by the PDF of wind speed and the wind speed-power curve. The detailed calculation steps and parameters can be found in [22]. The output power of the PV can be calculated by predicting solar radiation. However, the solar radiation reaching the Earth’s surface depends upon the climatic conditions of a location [23]. The calculation steps of power of PV are detailed in [24]. The parameters used for solving DEED in this study are shown in Table 1. For each case, 20 independent trials were conducted to compare the solution quality and convergence characteristics. All experiments were executed in the MATLAB 2016a computational environment on a Pentium-core 3.2 GHz personal computer with 8 GB of RAM.
4.1. Case Study
4.1.1. Case 1: Classic 5-Generator System
This test system consists of five thermal generators with considering the transmission loss. The generating system data and B coefficients to compute transmission line losses are adapted from [25]. The best cost and emission solutions are given in Tables 2 and 3, respectively. The compromise solution considering economic and emission is presented in Table 4. represents the dispatch hour, and denotes the output power allocated to the thermal generator. is the transmission line losses. and are the fuel cost and emission, respectively. The optimal cost and emission solutions are 44132 $/day and 17852 lb/day, respectively. The minimum cost and emission of the best compromise solution are 47503 $/day and 18066 lb/day, respectively. Table 5 gives the comparison among different methods, such as SA, MSL, EP, PS, PSO, MODE, and AGB-MOCDE. It is obvious that the proposed method can obtain solutions with lower total economic costs and emission than other reported methods.
4.1.2. Case 2: Classic 10-Generator System
This case consists of ten thermal generators with considering the transmission loss. The generating system data and B coefficients can be obtained in [26]. The best cost and emission solutions are given in Tables 6 and 7, respectively. The compromise solution considering economic and emission is presented in Table 8. The optimal cost and emission solutions are 2471526 $/day and 291988 lb/day, respectively. The minimum cost and emission of the best compromise solution are 2514885 $/day and 299585 lb/day, respectively. The simulation results of the proposed methods are compared in Table 5 with NSGA-II, EP, PSO, SQP, HMO-DE-PSO, RCGA, MADM, IBFA, and MODE. It is obvious that the proposed method can obtain solutions with lower total economic costs and emission than other reported methods.
To test the convergence characteristics of the proposed method, we conducted a comparison experiment in the same experiment situation for Case 2. Two typical methods, NSGA-II and MODE, are used to compare with SHHA. 20 independent trials were conducted to compare the convergence characteristics. Figure 1 shows the convergence comparison when aiming to optimize the fuel cost. The figure indicates that the SHHA converges faster than other methods. In addition, the statistical results of the 20 trails are given in Table 9. It can be seen that the SHHA outperforms the other methods regardless of minimum cost, average cost, and maximum cost. The SHHA has a lower standard deviation which also shows a better robustness (Figures 2 and 3).
4.1.3. Case 3: IEEE 57 Bus System (Seven Thermal Units, Three Wind Generators, and Two PV)
In this case, the proposed method is used to solve a DEED on an IEEE 57 bus system. Three PV and two wind generators are incorporated into this system. The data and information of this system can be obtained from [34]. The transmission loss coefficients were taken from [24]. The predictive wind power and PV power for the 24 hours of the next day are shown in Figure 4.
The optimal load allocation scheme considering the wind and solar power unit is given in Table 10. The cost and emission are 111203.2 $/day and 115221.2 lb/day, respectively. MODE and NSGA-II were used for comparison. It can be seen that the proposed method has better performance both in cost and emission. Moreover, we studied the impact of the introduction of RES on the system. Figure 5 shows that both cost and emission have improved. The cost and emission before introducing wind power and PV generators are 119305.6 $/day and 127416.7 lb/day, which means the introduction of RES will save about 8102.4 $ per day (approximately 2,957,376 $ per year) and reduce emissions 12195.5 lb per day (approximately 4,451,357 lb per year) (Table 11).
4.1.4. Case 4: Case 3 with Energy Storage
In this case, the ES unit is considered on the basis of Case 3. Relevant parameters of thermal units, wind power generators, PV generators, and transmission loss coefficient are consistent with Case 3. Parameters of ES can be found in Table 1. The initial value of SOC is set to 0.5, and the SOC at each hour can be accessed bywhere is the charge efficiency, is the discharge efficiency, and is the ES energy capacity. Table 12 lists the optimal compromise solution. The optimal compromise solution corresponding to the cost is 111317.6 $/day and the emission is 14565.4 lb/day. At the beginning of the day, the power load is in the low stage, at which time the power demand is less and the ES is in charge mode to store energy for future peak shaving. It is assumed that the initial SOC is 0.5 and the total charge power per hour should meet the upper limit. In the peak load stage, the ES is in the discharge mode to reduce the power generation of the grid. Finally, the ES state will be restored to 0.5 SOC at the end of the day. Figure 6 shows the effect of ES on load peak shaving and valley filling. Figure 6 reflects the impact of the introduction of ES on the Pareto optimal solution set. It can be seen that ES is beneficial to both cost and emission (Figure 7).
4.2. Effect of High-Level Strategy
In the hyper-heuristic method, the high-level strategy is used to evaluate each low-level operator and then determine the chosen probability of each operator for the following searching iterations. Figure 8 shows the change of the chosen probability of each low-level operator and it can be seen that the operators are adapting to the whole iteration. The operators show a significant difference in chosen probability at different iterations, which means that they have a different searching performance at different stages of the evolution process. Here, we only discuss the role of high-level strategy, without regard to the effect of the parameter setting of low-level operators.
4.3. Effect of Elitist Nondominated Sorting Strategy
Population diversity is an important indicator to evaluate convergence and search performance of algorithms, especially in dealing with optimization problems with nonlinear and multipeak characteristics. To demonstrate the effectiveness of the proposed nondominated sorting strategy, we recorded the results of individual selection with and without the elite sorting strategy. Figure 9 shows the comparison results, where the red points represent the distribution of individuals screened by using the elite strategy. We can find from the figure that the population has better performance in individual distribution. The elite strategy effectively keeps the population diversity.
4.4. Effect of Constraint Handling
Constraint handling presented in Section 3 aims to transform the infeasible individuals into high-quality feasible individuals. To demonstrate its effectiveness, we record the average cost with and without the handling procedure in the same experiment situation. A constraint handling method in [35] is used for comparison, and we have further improved it on the basis of this method. The difference in power demand is no longer added directly to the adjusted operator but multiplied with a random factor . This will make the adjustment of power balance constraints more balanced without losing the characteristics of random allocation. Table 13 shows that each case with the proposed constraint handling can find a better solution. When a feasible solution is obtained at the current hour, it may make it impossible for to meet the power balance and ramp rate constraints at the next hour. Failure times refer to the number of times of the situation that the constraints cannot be met. Lower failure times mean an effective constraint-handling strategy. We record the reduction percentage in the number of failures. It can be seen from Table 12 that the cost and emission are all improved and the number of failure times has been further reduced.
5. Conclusions
Building a low-carbon future has been an important and urgent task under Paris Global Agreement. Dynamic economic emission dispatch has long been a complex problem which aims to effectively save the fossil fuel cost, relief energy waste, and reduce environmental pollution. Moreover, as renewable energy is considered, the complexity of the DEED is further enhanced. In this paper, a mathematical model of the multiobjective dynamic economic and environmental load distribution problem including wind energy, solar power generation units, and ES is established. A multiobjective hyper-heuristic algorithm is proposed to solve the problem. Four test cases in 24-hour time period are tested to verify the effectiveness of the proposed algorithm.
The main contributions of this study can be summarized as follows: (1) we provide a new way of solving such coupled spatial-temporal scheduling problems in the power system. The hyper-heuristic algorithm is used to solve multiobjective DEED problems. Researchers can replace the low-level operators of our algorithm with their own methods to further improve the search capability. (2) For improving search capability and diversity distribution of Pareto front, an elitist nondominated sorting strategy is adopted to guide the screen of parent individuals. (3) An effective strategy is introduced to transform infeasible solutions towards feasible ones to enhance the quality of solutions and the convergence rate of the proposed algorithm.
In future work, we plan to further study the high level of the hyper-heuristics. The proposed algorithm is a selection hyper-heuristic algorithm, in which the low-level operators need to be pregiven. Further work will explore the generation hyper-heuristics to automatically designed operators. On the other hand, a renewable energy source such as electric vehicles will be integrated into the system dispatch to comprehensively study the power system.
Data Availability
The authors pledge to provide all the codes and the data underlying the findings of the study.
Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.
Acknowledgments
The authors would like to thank the National Natural Science Foundation of China (grant nos. 61773105, 61533007, 61873049, 61873053, 61703085, and 61374147) and the Fundamental Research Funds for the Central Universities (grant no. N182008004) for supporting this research work.