Economic Emission Dispatch Considering Renewable Energy Resources—A Multi-Objective Cross Entropy Optimization Approach

: The conventional electrical power system economic dispatch (ED) often only pursues immediate economic beneﬁts but neglects the harmful environment impacts of gas emissions from thermal power plants. To address this shortfall, economic emission dispatch (EED) has drawn a lot of attention in recent years. With the increasing penetration of renewable generation, the intermittence and uncertainty of renewable energy such as solar power and wind power increase the difﬁculties of power system scheduling. To enhance the dispatch performance with signiﬁcant penetration of renewable energy, a modiﬁed multi-objective cross entropy algorithm (MMOCE) is proposed in this paper. To solve multi-objective optimization problems, a crowding–distance calculation technique and a novel external archive mechanism are introduced into the conventional cross entropy method. Additionally, the population updating process is simpliﬁed by introducing a self-adaptive parameter operator that substitutes the smoothing parameters, while the solution diversity and the adaptability in large scale systems are improved by introducing the crossover operator. Finally, a two-stage evolutionary mechanism further enhances the diversity and the rate of convergence. To verify the efﬁcacy of the proposed MMOCE, eight benchmark functions and three different test systems considering different mixes of renewable energy sources are employed. The dispatch results by the proposed MMOCE are compared with other multi-objective cross entropy algorithms and published heuristic methods, conﬁrming the superiority of the proposed MMOCE over other methods in all test systems.


Introduction
The economic dispatch (ED) is a fundamental issue in electrical power system scheduling that aims to maximize the economic profits by the optimal allocation of the output of each generator unit [1]. It is shown that more than 10% of the total energy consumption can be saved by means of economic dispatch [2]. However, in the last decades, environmental issues have drawn substantial attention worldwide, and much effort has been made to mitigate the negative effects of climate change and environment pollution [3]. As such, managing pernicious gases emissions has become an important consideration when it comes to the ED problem. The economic emission dispatch (EED) retains the original characteristics while incorporating the emission factors, resulting in a multi-objective optimization problem, which simultaneously minimizes the generating cost and the emission level to the lowest possible values.
In China, about 70% of electricity supply comes from coal-fired power stations; addressing climate change while meeting future energy needs will inevitably impose greater pressures on traditional energy sources [4]. Today, high penetration of renewable energy sources is expected for electric power systems, and many attempts to integrate renewable energy resources into the EED problems are reported [5]. As of 2015, more than 9% of electricity generation comes from renewable resources [6]. Renewable energy refers to the resources that can be used repeatedly to produce energy, such as solar energy, wind energy, biomass energy, geothermal energy, etc., [7]. It is also often referred to as alternative energy. Utilizing abundant renewable energy resources reduces the emission of harmful gases to a great extent. The majority of studies revealed that integrating renewable energy resources such as wind energy conversion systems (WECS) and PV solar systems into the power system is a promising technology and practically feasible. Nevertheless, the existing grid faces unpredictable and intermittent power supply from the resources, particularly from wind and photovoltaic solar. The electric power production from renewable energy resources such as wind speed and solar radiation can be very indeterminate, variable with time, and undispatchable with limited control, imposing major challenges on the scheduling of power systems [8]. Their intermittent power output imposes various difficulties on the scheduling, operation, and control of the power system networks. Therefore, handling renewable energy resources requires sophisticated planning and operation scheduling as well as state-of-the-art technologies.
As a research hotspot in the field of power system optimization, some progress has been reported on solving the EED problems incorporating renewable energy resources in recent years. In [9], the multi-objective economic emission dispatch problem is considered, which combines heat and wind power generation in a large micro-grid (MG), and IEEE 30-bus and 69-bus systems are used to represent a large MG. In [10], a solar-wind-thermal ELD problem integrating a pumped-storage hydraulic unit is considered. In [11], there is an analysis of the economic dispatch and unit commitment problem where five scenarios are considered to predict the power system operations by 2030. In [12], the EED model of an IEEE 30-bus system incorporating solar photovoltaic, wind and battery storage is built and multiple objectives are considered. The authors of [13] present an EED model that integrates thermal, natural gas, and renewable energy systems, considering both emission levels and generation cost. These reports describe the EED model considering renewable energy resources in detail, and various optimization strategies are utilized. Since there exist a number of EED problem types with renewable energy, most studies only consider one type of EED model, and the proposed methods may be difficult to apply to other kind of models. It is therefore meaningful to consider multiple types of EED models with different renewable energy sources. Under this consideration, this paper will investigate three different types of EED problems integrated renewable energy resources.
For the EED problem, a number of algorithms have been attempted so far. Generally speaking, these algorithms can be mainly grouped into two categories [14,15]. In the first category, all objectives are combined linearly so that the EED problem can be regarded as a single objective optimization problem [16][17][18]. Meanwhile, given the nonlinear and nonconvex nature of the dispatch problem, the classical deterministic techniques such as weighted sum method, −constraints, and linear programming may fail to solve the problem due to the sheer computational complexity of the problem, inclination to fall into local optima, and discontinuities and nonsmoothness in the process of handling functions. Hence, nature-inspired metaheuristic algorithms have attracted substantial interests. The drawbacks of this approach, however, are multi-faceted: First, only one solution is generated in a single run, so it is time-consuming. Second, the exact weight parameters in the objective function are rather difficult to optimize.
The second category of algorithms is to deal with the two objectives in the EED problem simultaneously. The Pareto front can be easily obtained with a single run so that the optimal solution can be chosen according to different principles. Many Pareto-based multi-objective algorithms have been adopted to solve the EED problem successfully, such as NSGAII, MOPSO, and SPEA. Qiao [19] employed a self-adaptive multi-objective differential evolution algorithm on a novel dynamic economic emission dispatching framework Sustainability 2021, 13, 5386 3 of 33 integrating both electric vehicles and wind farms and demonstrated that the proposed method can improve the results efficiently in different test systems based on 10-unit generators. Habibi [20] applied a multi-objective particle swarm optimization (MOPSO) algorithm by optimizing the objective function of the storage and additional cost to solve the EED problem incorporating stochastic wind in multi-area power systems. Jiang [21] used gravitational particle swarm optimization algorithm (GPSOA) to solve the EED problem considering wind power availability for the wind-thermal power system. Li [22] established a dynamic economic emission scheduling model incorporating wind power, solar power, and hydropower under tradable green certificate and proposed a new multiobjective moth-flame optimization (MOMFO) to resolve the model. Liao [23] presented a multi-objective optimization by learning automata (MOLA) and implemented it on the EED model such that the objective functions are formulated as a combination of cost, emission and voltage stability. Ghasemi [24] developed a new Honey Bee Mating Optimization (HBMO) with blended online learning mechanism for the EED model with uncertain wind power. However, these methods may still fail to obtain better results despite of their complexity or they have several parameters to pre-determine in the process of optimization that may greatly affect the optimization results. It is therefore worthwhile to develop a simple yet efficient algorithm to solve multiple objective optimization problems.
The cross entropy (CE) method proposed by Robinstein [25,26] is an innovative metaheuristic algorithm based on rare events resampling and the Kullback-Leibler distance minimization. The main idea of the CE method is to select elite samples from the best performing populations and generate new populations according to the distribution of elite samples. The superiority of CE lies in its diversified structure and its simple process of updating populations. Hence, the CE method has been successfully applied to a range of engineering problems [27][28][29] including the ED problem [30]. In [31], the CE method was firstly extended to multi-objective optimization and achieved good performance. In recent years, the multi-objective cross entropy method has been further improved and applied to a range of problems, such as finance [32], computer network [33], and engineering design [34]. For example, in [35], a cascaded algorithm combining the cross-entropy (CE) and Tabu search (TS) was successfully applied to building an effective boiler efficiency prediction model and improved the economic benefit and reliability of generating units. Dorini et al. [36] proposed a novel algorithm based on the noisy cross-entropy sensor locator (nCESL) to detect accidental and/or intentional contamination in water distribution systems. Sun et al. [37] improved the CE method and used it to solve the multi-objective energy routes problem in a WPTN system. Perelman et al. [38] extended the conventional CE method to multi-objective combinatorial optimization of water distribution systems design and demonstrated its better robustness comparing with NSGAII. Bekker et al. [39] adapted the CE method to multi-objective optimization and applied the proposed algorithm to a dynamic, stochastic problem and demonstrated its effectiveness and efficiency. However, it is worth mentioning that only limited papers reported to have applied the multi-objective CE methods to solving the EED problems so far. For instance, in [40], a cross entropy optimization based on decomposition was proposed to solve a multi-objective optimization problem for a model of a wind/hydro/thermal/photovoltaic power system with 10 generators. However, the performance of most CE methods tends to be less optimistic in solving relatively large-scale EED problems. In addition, they have a number of smoothing parameters, which may cause premature convergence and limit the algorithm performance. Therefore, it is necessary to improve the design of the multi-objective cross entropy algorithm with far fewer smoothing parameters and improve their applicability in solving complex EED problems.
As mentioned earlier, the conventional CE method has the advantage of a versatile structure and simple updating procedure. Nevertheless, the drawbacks of the CE method are also evident: they easily fall into local optimum and need to predetermine several constant parameters such as the smoothing parameters (different values of constant parameters may result in different optimization performance). While in the field of multi-Sustainability 2021, 13, 5386 4 of 33 objective optimization, the performances of most existing multi-objective CE methods are still yet to be satisfactory when dealing with large scale problems. In this paper, a modified multi-objective cross entropy algorithm (MMOCE) is presented to solve EED problems considering the presence of renewable energy sources. First, a self-adaptive operator substitutes the smoothing parameters in the updating formula for the mean value and the standard deviation. This new updating process is proposed to reduce the number of pre-defined constant parameters, which makes the algorithm simpler and more stable. Second, the crossover operator is employed in present population and the external archive to improve the global search mechanism and the scalability to large-scale test systems. Third, different parameter evolutionary mechanisms as stated in [41], which are implemented to calculate the mean value and the standard deviation of the present population according to population generation, are introduced into MMOCE to enhance the diversity of solutions and to further speed up the convergence. With these improvements, the algorithm efficiency is greatly enhanced. Eight benchmark functions and two performance indicators are used to verify the feasibility and superiority of the proposed MMOCE by comparing with other state-of-the-art algorithms such as NSGAII, MOPSO, MOEA/D, and PESA2. Additionally, in this paper, three different types of EED problems with different scales and combinations of renewable energy resources are introduced to testify the superiority of the proposed method. The first problem is a combined heat, emission, and economic dispatch (CHEED) problem integrating wind and solar power generations. The second problem is an IEEE 30-bus and six-generator system that incorporates stochastic wind and solar power. The third problem is a 40-unit combined emission economic dispatch problem with wind penetration. The results obtained from the above three examples confirm that the proposed method outperforms other competing algorithms.

Mathematical Model
This paper mainly focuses on three different environmental economic dispatch models considering renewable sources, including the combined heat, emission, and economic dispatch (CHEED), and applies the proposed method to a modified IEEE 30-bus and six-generator system with renewable energy and combined emission economic dispatch problems with wind penetration. The overall dispatching framework of the first model are shown in Figure 1.

Problem Formulation of CHEED
The system considered in this paper contains conventional generators, cogeneration units, and heat-only units. The feasible operating region (FOR) of the cogeneration units is illustrated in Figure 2 [42], which is enclosed by the boundary curve ABCDEF. For the CHEED problem, the power is derived from the thermal units and cogeneration units while the heat is derived from cogeneration units and heat-only units. Under the condition of ensuring the power balance of the system and satisfying the constraints, the output of each unit is reasonably allocated to achieve the goal of minimizing the cost of heat and power Sustainability 2021, 13, 5386 5 of 33 production while minimizing the emission level. Thus, CHEED can be mathematically elaborated as follows. Total cost is defined as: 2.
Emission is defined as:

Constraints
In this paper, the renewable energy is incorporated into the CHEED problem and plays the role of negative loads in order to decrease the demand load. Simultaneously, the electricity and heat production and capacity of each unit are subject to their respective constraints. The power transmission loss can be calculated by B-matrix coefficients as follows [43]: The overall dispatching framework of this model are shown in Figure 3. The cost of conventional thermal units primarily refers to fossil fuel cost. Meanwhile, in order to make the cost function more accurate and realistic, the valve effect, which is formulated as a sinusoidal function and added to the basic cost function, is also taken into consideration. Then, the ultimate cost function of convention thermal units can be expressed as follows in Equation (10).

Cost of Wind Energy
As one of the most common renewable energy resources, wind power is characterized by its intermittency and uncertainty. Provided that wind farms are unable to provide enough available scheduled power, an independent system operator (ISO) will take the responsibility of handling the deficit in order to maintain spinning reserves. This case is deemed as overestimation of renewable energy and the power generation cost will have to increase in order to maintain the spinning reserves. In contrast, if these farms provide excessive power, the ISO should accept the penalty. This case is called underestimation. In general, the total cost of wind power includes direct cost, reserve cost, and penalty cost [44].
The direct cost of wind power can be expressed as follows [45].
In addition, the situations of overestimation and underestimation may always exist and thus increase the overall costs. In summary, overestimation will lead to the reserve cost and underestimation will lead to the penalty cost. They can be described, respectively, as follows [45].
C RC,i (P sw,i − P wa,i ) = K RC,i (P sw,i − P wa,i ) = K RC,i P sw,i 0 In this paper, the Weibull probability density function (PDF) [46,47] is used to represent the wind speed distribution and Figure 4a illustrates the Weibull PDF of the wind speed data by employing 8000 Monte-Carlo scenarios. Thus, the probability of wind speed can be described as below, where c is the scale factor and k is a shape factor.
The mean of the Weibull distribution is formulated as: The power output of a wind turbine is approximately defined as follows: Next, wind power probabilities can be piecewise calculated by [48]:

Cost of Solar Photovoltaic Power
As the most popular renewable source, similar to wind power, the total cost of solar power also contains direct cost, reserve cost, and penalty cost. Here, the direct cost involved with solar power can be computed as: Reserve cost for the solar power plant is defined as: C SR,j P ss,j − P sa,j = K SR,j P ss,j − P sa,j = K SR,j × f s P sa,j < P ss,j × P ss,j − E P sa,j < P ss,j where f s P sa,j < P ss,j denotes the probability of event of solar power shortage associated with the scheduled power P ss,j . E P sa,j < P ss,j denotes the expectation of solar power below P ss,j . Penalty cost for underestimation is defined as [45]: C PS,j P sa,j − P ss,j = K PS,j P sa,j − P ss,j = K PS,j × f s P sa,j > P ss,j × E P sa,j > P ss,j − P ss,j Similarly, f s P sa,j > P ss,j denotes the probability of event of solar power shortage associated with the scheduled power P ss,j . E P sa,j > P ss,j denotes the expectation of the solar power above P ss,j .
In this paper, the Longnormal PDF is illustrated in Figure 4b and adopted to describe solar irradiance (G sr ). According to the Longnormal PDF with mean µ and standard deviation ς, the probability of solar irradiance can be calculated as follows [49]: The mean of Lognormal distribution can be expressed as: Based on the analysis from [45], the overestimation cost of Equation (21) is calculated as: C SR (P ss − P sa ) = K SR (P ss − P sa ) where P sn− represents the relevant power less than the scheduled power. f sn− denotes the relative frequency of occurrence of P sn− . N − stands for the number of pairs (P sn− , f sn− ) generated for the PDF. According to [45], the underestimation cost of Equation (22) can be computed as: where P sn+ represents the relevant power more than the scheduled power. f sn+ denotes the relative frequency of occurrence of P sn+ . N * is the number of pairs (P sn+ , f sn+ ) generated for the PDF.

Emission Function
It is generally understood that harmful gases (mainly SO x and NO x ) will be emitted into the atmosphere when electric power is generated from conventional fossil fuels. The functional relationship between pollutants emissions and the produced power from conventional thermal generators is defined as: In the past decade, the aggravated global warming and air pollution has driven many countries to commit to reducing carbon emissions [46]. The carbon tax (C tax ) emerged in order to penalize greenhouse gas emissions and encourage enterprise to invest renewable energy. Therefore, the cost of emission can be expressed as:

Objective Functions
Similar to conventional EED problems, the objective of this mathematical model is to optimize the output of each generating unit to meet the load demand and simultaneously reduce the generation cost and emission as much as possible under the condition of satisfying various constraints. The difference lies in that decision variables here consist of the produced actual power and generator bus voltages. The first objective function is formulated incorporating the cost of conventional thermal units, direct cost, reserve cost, and penalty cost as mentioned above. The second objective function principally takes environment factors into account and aims to minimize the total emission of conventional thermal generators. Moreover, the carbon tax is also incorporated into the objective function as discussed above. Thus, all the objective functions can be expressed as follows:

Constraints
There are equality and inequality constraints that must be met. Equality constraints refers to the balance between active and reactive power generated with load demand and transmission losses. Hence, equality constraints are formulated as follows: The inequality constraints cover the operation limits on equipment and security constraints on lines and load buses. Generator constraints: Q min ss,k ≤ Q ss,k ≤ Q max ss,k , k = 1, 2, . . . , N SG (38) Security constraints: Based on the above formula, Equations (33) and (34) represent the limitation of the active and reactive power generated by thermal generators. Equations (35) and (36) are the limitations of the active and reactive power generated by wind generators. Equations (37) and (38) formulate the limitation of the active and reactive power generated by solar PV. Equation (39) represents the constraints on voltage of generator buses. Equation (40) represents the limitation of voltage imposed on NL numbers of load buses. Equation (41) denotes the line capacity constraints with total nl numbers of lines in the network.

Combined Emission Economic Dispatch Problems with Wind Penetration
Analogous to the conventional EED problem, the combined emission economic dispatch problem with wind power penetration also aims to achieve optimal scheduling of power generators to minimize the fuel cost and emission generated by thermal generators while simultaneously satisfying all the equality and inequality constraints. The objective functions of this system are the same as defined in Equations (10) and (27). The only difference is that wind energy becomes a part of energy mix and supplies a portion of the power demand. Similarly to the CHEED model mentioned above, the wind energy is counted as negative loads. The power balance constraint is converted into: (42) where N is the number of thermal generators and P T i stands for the output power of i-th thermal generator. P wp and P D represent the outputs of the wind farm and the load demand of the system, respectively. P Loss is the transmission loss, which can be formulated by the B matrix coefficients: where B ij , B oi , and B oo are matrix coefficients of the transmission loss.

Overview of the Cross Entropy Method
As mentioned above, the cross entropy (CE) method was first proposed in 1997 by Robinstein [25,26] based on rare events resampling and the Kullback-Leibler distance minimization. The cross-entropy method was initially employed in mathematics for the calculation of low probability events and then extended to numerical optimization. It is an optimization method that updates new parameters by collecting information in each iteration. The implementation procedure of CE is briefly divided into two major parts: First, random samples are generated based on the probability density function (PDF) and obtain the fitness values according to objective functions. Second, the elite set is used to update the PDF and create new samples.
For simplicity, the critical process of the conventional CE method is presented here. It is assumed that every individual i solution has a D-dimensional position vector, which is formulated by X k i = x k 1 , x k 2 , . . . , x k D at iteration k. µ n , σ n are defined as the mean value and the standard deviation of the dimension n, respectively. The components of each individual solution are selected independently and specifically generated separately based on the N µ n , σ 2 n distribution. Then, the fitness values of all individuals are calculated and sorted in ascending order. For the minimization problems, ρ × N p best individuals are chosen as the elite set X elite to help guide the updating process of the mean value and the standard deviation where N p is the population size and ρ is a quantile. The population can be generated by Equation (44), where x j n stands for the value of the j-th individual in the Sustainability 2021, 13, 5386 11 of 33 n-th dimension and randn() is a random variable based on normal distribution function N(0, 1). The evolution of an individual to a new individual is conducted in terms of the mean and standard deviation of the elite set. The updating process are normally given as follows: However, it is far from being enough to merely rely on the elite set to update samples. Hence, smoothing parameters are introduced for the improvement of the searching capacity. In Equations (48) and (49), α and β are defined as smoothing parameters. µ(k + 1) and ς(k + 1) are employed to produce new samples, hoping to discover better solutions. Here, α and β are both constants (typically ranging from 0.8 and 0.99), corresponding to a fixed version of the CE method (FCE).
Moreover, it should be noted that FCE is able to show good performance in many cases but not in some other situations. Therefore, a dynamic version of the CE method (DCE) emerges. For µ, the same fixed parameter described in Equation (48) is used; the dynamic adjustment of smoothing parameters for standard deviation is expressed by: where β n denotes the value of β of the n-th decision variable in the k-th iteration. C is a constant between 0.8 and 0.99, and q is an integer between five and ten.
In the global search, the population is regenerated after the completion of the updating process of the mean value and the standard deviation. Then, the individuals in the new population are rearranged from largest to smallest according to fitness values and the elite set is selected from among them. The above steps are repeated for a specific number of iterations.

The Proposed Modified Multi-Objective Cross-Entropy Algorithm
Since most of the practical engineering optimization problems are multi-objective optimization problems, the CE method needs to be improved to adapt to such problems. In this paper, in order to enhance the adaptability of the algorithm in multi-objective optimization, the CE method is extended to a multi-objective algorithm due to the introduction of the framework of NSGA-II [50]. This method couples with a fast non-dominated sorting approach to choose the elite set and stores the Pareto solutions in the external archive. The modified multi-objective cross entropy algorithm (MMOCE) maintains the advantage of fast convergence. However, there is a risk of falling into local optimum in the process of global search, and the ability of the algorithm to solve large scale system is truly unsatisfactory. Hence, a crossover operator is integrated into the algorithm to increase the diversity of solutions and improve the evolutionary efficiency. Additionally, pre-determined parameters in the conventional CE method may influence the optimization results because of improper choice. Although DCE and FCE can be sometimes effective, the options of proper parameters are not so simple. To tackle this problem, a new version of the mean value and the standard deviation updating process is proposed. With no need for setting a constant parameter, the process of updating the mean value and the standard deviation is self-adaptive. This mechanism simplifies the algorithm and simultaneously ensure the convergence. Moreover, different parameter evolutionary mechanisms as stated in [41] are employed to this algorithm, which enhance the diversity of solutions and speed up the convergence rate. Analogous to [41], in order to facilitate the balance between exploitation and exploration searches, the whole iterative process is divided into the diversification stage and intensification stage. The diversification stage highlights the diversity of sampling points, while the intensification stage guarantees the rapid convergence. More details of the different parameter evolutionary mechanisms can be found in [41]. The specific introduction of the self-adaptive parameter updating process and crossover mechanism is given as follows.

Updating of Self-Adaptive Parameter
As mentioned earlier, the updating process of the mean value and the standard deviation involves the smoothing parameter. These constant parameters are often predetermined and application-dependent. To avoid this, a self-adaptive parameter updating process is introduced here, and the mean value and the standard deviation can be updated in a self-adaptive way. This self-adaptive updating process can be implemented as follows at iteration k: where rand represents a random value ranging from 0 and 1. Considering the widespread use of the self-adaptive thinking in real-life optimization, Equation (50) can be reformulated as Equation (52). From the listed formula above, we can see that the mean value and standard deviation are self-adaptive after initialization and avoid the trouble of setting parameters in advance. The test functions below can demonstrate the validity and effectiveness of the method.

Crossover Operator
To increase the diversity and expand the search range, the crossover operator is introduced. The genetic information of elite individuals in the external archive is introduced to guide the evolution. First, the random individuals X i and X j are chosen from the present population P and the external archive, respectively. Second, part of the code of the two individuals is exchanged, generating two new individuals. The above operation is repeated until a new population P new with N p individuals is generated, where N p is the size of the population. This operation can be described as follows: where x new i and x new j are the corresponding n-dimension encoding after the crossover operation. rand is a random number between 0 and 1; P c is the probability of crossover. Figure 5 presents the flowchart of the proposed algorithm.

Implementation of MMOCE for EED
This subsection presents the proposed MMOCE procedure for solving the EED problems with renewable energy in detail.

•
Step 1. Input data. The feasible range of control variables, the population size N p , the maximum number of function evaluations Max_FES, the maximum size of external archive Ar max , and the probability of crossover P c .

•
Step 2. Initialization. The initial population P is formed by normal distribution. It is assumed that the dimension n over ranges [X min , X max ], which denote the lower and upper limit, is assigned based on N µ n , σ 2 n , where µ n is a random number between the lower and upper limit and ς n is set to 10 · (X max −X min ), and FES = 0.

•
Step 3. The mean value µ and standard deviation ς are updated by Equations (51)- (53). The endpoint value method is used as the boundary constraints handling strategy for the population P.

•
Step 4. The objective function values are calculated and coupled with a fast nondominated sorting approach and crowding-distance metric to get the external archive.

•
Step 5. The crossover operator is carried out by Equation (54).

•
Step 6. The elite set X elite is selected from the external archive and the different parameter evolutionary mechanisms are applied to update the mean value µ and standard deviation ς. FES = FES + N p . • Step 7.
Step 2 to Step 6 are repeated. If the termination criterion is met, i.e., Max_FES is reached, the above procedure is stopped.

Experiments Based on Benchmark Functions
In order to verify the optimization performance of the proposed algorithm, eight medium-large scale test functions ZDT1-3, UF1-3, and WFG1-2 are selected. Meanwhile, five representative multi-objective evolutionary algorithms (NSGAII, MOEA/D, MOPSO, PESA2, MODE [51]) are also tested for comparison purposes. In addition, the DCE and FCE are both extended to multi-objective algorithms, denoted as DMOCE and FMOCE. The differences between the proposed algorithm and conventional CE variants are mainly reflected in the choice of smoothing parameters. Thus, four different multi-objective cross entropy optimization algorithms (DMOCE, FMOCE, MOO CEM [39], SMOCE [52]) are also implemented for comparison purposes. Moreover, this paper adopts the widely-used inverted generation distance and MaxSpread as performance indicators in multi-objective evolutionary algorithms. The inverted generation distance (IGD) calculates the Euclidean distance of the nearest Pareto optimal solution to the nearest Pareto optimal solution set, while MaxSpread (MS) is used to evaluate the coverage degree of the non-dominated solution obtained by the algorithm to Pareto theory optimal solution set. Their definitions are as follows: where Q denotes the Pareto optimal solution set obtained by the algorithm, |Q| is the number of non-dominated solutions in Q, Q * denotes the theoretical Pareto optimal solution set, |Q * | denotes the number of non-dominated solutions, and f where f i represents the value of the i-th objective function achieved by the algorithm corresponding to the solution and F i represents the objective value of the i-th objective of the theoretical optimal solution. If MS = 1, it means that the non-dominated solution found by the algorithm can completely cover the Pareto theoretical optimal solution.
In this paper, each algorithm runs 30 times independently to obtain the average values. The population size of each algorithm is set to 100, which is the same as [53][54][55], and the max function evaluations (FES) of each test function are 15,000 for all algorithms. The external archive set size is set to 100. The crossover probability of MMOCE is 0.9 and elite set size is set as 10. For NSGAII, the crossover probability is 0.9, the mutation probability is 0.1, and the mutation distribution index is 20. The mutation scale factors and crossover probability of MODE are 0.5 and 0.3, respectively. The probability of crossover and mutation of PESA2 is both 0.5 and the number of grids is 7. Finally, MOPSO uses the following parameters: the number of grids is set to 7, the probability of mutation is 0.1, both the individual and the global learning factors are 1 and 2 respectively, and the inertia weight is 0.5. All the simulations are run using MATLAB 2016a and on a 2.4 GHZ and 8 GM RAM computer in Windows 10.
Tables 1 and 2 are a summary of the influence of different constants in Equation (52) on the optimization performance of the proposed algorithm. IGD denotes the proximity between the non-dominated solution and the optimal Pareto front end, and MS can evaluate the coverage of the non-dominated solution to the optimal solution. It is evident that when the constant is 0.382, the proposed algorithm can achieve overall better optimization performance compared to other values. When the constant is 1.0, it is clear that the optimization performance of the algorithm deteriorates dramatically. The primary reason is the lack of perturbation, and the diversity of solutions decreases sharply and premature convergence occurs. Therefore, the diversity of solutions is of great significance to further improving the searching ability of the algorithm. Compared to other values, when the constant is set to 0.382, the eight test functions with different characteristics show better optimization performance.  Tables 3 and 4 present the IGD data of MMOCE and other algorithms. It can be observed from the tables that the MMOCE has better IGD characteristics than other algorithms for the function ZDT1-3 and WFG1-2, and it has the second best IGD for UF1-2 after MODE. Although not as good as other algorithms for UF3, it is obvious that the proposed algorithm is generally better than other algorithms on the IGD index. It implies that the non-dominated solutions obtained by the proposed algorithm is the closest to the theoretical optimal frontier.  Tables 5 and 6 compare the performance of MMO CE algorithm and the MS index with other algorithms. It can be seen from the tables that the MMOCE outperforms other algorithm on the MS index, which implies that the proposed algorithm has the maximum coverage of the non-dominated solution to the optimal solution. Meanwhile, it indicates that the MMOCE has the best diversity among all the algorithm. Additionally, Table 7 lists the simulation time of different algorithms on benchmark functions. It is obvious that MMOCE has certain advantages that benefit from its simple structure.

Simulation Results on Combined Heat, Emission, and Economic Dispatch (CHEED)
This section discusses the implementation of the proposed MMOCE algorithm on the combined heat, emission, and economic dispatch problem in [56]. In this work, it is assumed that wind and solar energy each account for 5% of the total load requirements in all test cases. In this section, two power systems are used to examine the superiority of this proposed algorithm in solving the CHEED problems considering renewable energy sources. The first power system is composed of four conventional thermal units, two cogeneration units, and one heat-only unit. P d and H d are set to 600 and 150 MW th , respectively. Note that this power system considers the power losses. The second power system consists of one conventional thermal unit, three cogeneration units, and one heat only unit, and two scenarios are considered. In the first scenario P d and H d are set to 300 MW th and 150 MW th , while in the second scenario, P d and H d are set to 250 MW and 175 MW th , respectively. The detail data of both power systems can be referred to in [42].
In order to confirm the feasibility and effectiveness of the proposed algorithm, the simulated results are employed to compare with one newly proposed MSFLA [56]. For coping with the given optimization model, the MSFLA introduces the price penalty factor (PPF) to combine the total cost and total emission into a single objective optimization model. The PPF in this paper is set to 35,863.1065 $/kg in the first system and 28,297.8499 $/kg in the second system. Thus, for convenient comparison, the best solution of the MMOCE (termed as S n ) is selected from the Pareto front according to minF = C t + ΛE t [56]. The maximum number of iterations is 100, and the FES for all the algorithms is set to 5000. Tables 8-10 list the detailed results of the two power systems derived from MMOCE. It can be observed that these results obtained by MMOCE are clearly more environmental and cost-efficient, which demonstrates the effectiveness and superiority of the method over MSFLA [56], GA [57], SFLA [58], TLBO [59], and ISFLA [60]. For the first system, the fuel cost and emission are 16  Moreover, to further verify the effectiveness and superiority of the proposed MMOCE, three other multi-objective cross entropy (DMOCE, FMOCE, and MOO CEM) methods and two classic meta-heuristic multi-objective algorithms (NSGAII and MOPSO) are also used for comparison. Figure 6 show the Pareto optimal solutions of other algorithms in contrast to the proposed method, respectively, in the first system. There are two explicit points: first, the Pareto optimal front is distinctly superior to other several multi-objective algorithms. Second, the Pareto solutions obtained by MMOCE are more evenly distributed, and the diversity is also excellent.  Figures 7 and 8 show the Pareto optimal solutions of different algorithms, and Tables 12 and 13 present all the compromise solutions of two scenarios in the second system of all the algorithms in comparison. Additionally, Figure 9 shows the comparison between the proposed MMOCE and other algorithms regrading to the total cost and emission for all scenarios. It can be clearly seen that overall the solutions obtained by MMOCE are better in convergence and diversity than those of other algorithms.
MMOCE and other algorithms regrading to the total cost and emission for all scenarios. It can be clearly seen that overall the solutions obtained by MMOCE are better in convergence and diversity than those of other algorithms.

Simulation Results on IEEE 30-Bus with Six-Generator System
In this section, three different cases of IEEE 30-bus with a six-generator system with renewable resources [61] are introduced, aiming to investigate the performance of the proposed MMOCE. In Case 1, the system incorporates stochastic wind and solar power. The Case 2 is a solar-free stochastic wind power system, while the third system is solar-only. Simultaneously, in order to verify the feasibility and validity of the proposed MMOCE, the CMOPEO-EED [61] method and CNSGAII-EED [50,62] method are selected as the competitors to participate in the optimization of three cases. Especially in Case 1, the results are also used to make comparison with that of the SHADE method in the literature [45]. The details of the three different cases of IEEE 30-bus with a six-generator system, renewable resources, and adjustable parameters of CMOPEO-EED and CNSGAII-EED can be found in [61]. The population size of the MMOCE is 60, and the maximum number of function evaluations is set to 24000, which is the same as that in the literature.

Case 1: System Incorporating Stochastic Wind and Solar
In this case study, the IEEE 30-bus system contains two wind farms and one solar PV array to generate electricity. In other words, this system considers the influence of stochastic wind and solar power at the same time. Firstly, the SHADE is used to make a comparison with the proposed MMOCE, as it treats the problem as a single-objective optimization problem. In order to highlight the outperformance of the MMOCE in contrast to the SHADE, the three optimal solutions (termed as S 1 ∼S 3 ) are extracted from the MMOCE on the basis of F = F 1 + C tax F 2 . Table 14 lists the results of the MMOCE and the SHADE in detail, and the possible range of each unit is also included. It is evident that the two Pareto solutions (i.e., S 1 , S 3 ) obtained by the MMOCE are found to be superior to the optimal solution of the SHADE. S 1 has a total cost of 809.060 $/h and an emission of 0.833 t/h, S 3 has a total cost of 809.655 $/h and an emission of 0.744 t/h, while the SHADE has a total cost of 810.346 $/h and an emission of 0.891 t/h. It is not difficult to discover that the proposed MMOCE yields lower total cost and emission than SHADE. The Pareto front of MMOCE is shown in Figure 10.  In this paper, to compare with other two multi-objective algorithm CMOPEO-EED and CNSGAII-EED proposed in the literature, three different solutions (termed as A0, B0, and C0 in Figure 10) of the MMOCE are chosen from the Pareto front. As stated in [61], because of the preference to cost, A0 can be named cost solution. Similarly, B0 can be named emission solution, and C0 can be named as compromise solution. Table 15 lists the results of the three solutions in detail. It is shown that compared with other two algorithms, the result obtained from the proposed MMOCE is better as it achieves a decrease in both the cost and emission simultaneously in terms of cost solution and compromise solution. Though for emission solution, the MMOCE fails to show distinct superiority comparing with CMOPEO-EED, the proposed MMOCE still exhibits excellent convergence as a whole. Furthermore, DMOCE, FMOCE, MOO CEM, NSGAII, and MOPSO are also tested for comparison purpose with the proposed MMOCE. Figure 11 compares the Pareto fronts between MMOCE and the above competitor algorithms. The compromise solutions of those algorithms are picked out and listed in Table 16. Apparently, the optimal solutions obtained by MMOCE unfold better convergence performance than MOO CEM, NSGAII, and MOPSO. In comparison to FMOCE and DMOCE, the best Pareto fronts are visually close to each other. However, by means of the compromise solutions, it is clearly shown that MMOCE produces the cost of 804.2061 $/h and emission of 0.4994 t/h, which is the lowest among these algorithms. In summary, the proposed algorithm MMOCE has demonstrated the superior performance in the optimization of this case.

Case 2: the system only incorporating stochastic wind
In this case study, the test system only contains two wind farms and does not include solar PV. For the sake of detailed comparison with two methods presented in literature [61], three different solutions of MMOCE are selected, and the simulation results are listed in Table 17. It is clear that the proposed MMOCE can find the lowest values of cost and emission, which means it can find better solutions than CMOPEO-EED and CNSGAII-EED. Figure 12 compares the Pareto fronts between MMOCE and competitor algorithms, and the compromise solutions are listed in Table 18. The data curves show that the proposed MMOCE converges closer to optimal solutions than other algorithms. Furthermore, the result obtained by MMOCE, with 834.5583 $/h in cost and 0.5896 t/h in emission, stands out from all algorithms.   In this case study, the test system only contains two wind farms and does not include solar PV. For the sake of detailed comparison with two methods presented in literature [61], three different solutions of MMOCE are selected, and the simulation results are listed in Table 17. It is clear that the proposed MMOCE can find the lowest values of cost and emission, which means it can find better solutions than CMOPEO-EED and CNSGAII-EED. Figure 12 compares the Pareto fronts between MMOCE and competitor algorithms, and the compromise solutions are listed in Table 18. The data curves show that the proposed MMOCE converges closer to optimal solutions than other algorithms. Furthermore, the result obtained by MMOCE, with 834.5583 $/h in cost and 0.5896 t/h in emission, stands out from all algorithms.

Case 3: The system only incorporating stochastic solar PV
The test system in this case contains only solar PV. Similar to the above two cases, the compromise solution of the proposed MMOCE is selected to compared with the algorithms in [64], and Table 19 lists the results clearly. Figure 13 depicts the Pareto front of MMOCE and other algorithms in Case 3. It is not hard to see that MMOCE still has the merits in convergence and effectiveness.
Furthermore, the comparison between MMOCE and other competitor algorithms is conducted to verify the superiority of MMOCE while the results are listed in Table 20. Though it does not demonstrate evident advantage over MOO CEM with 836.1867 $/h in The test system in this case contains only solar PV. Similar to the above two cases, the compromise solution of the proposed MMOCE is selected to compared with the algorithms in [63], and Table 19 lists the results clearly. Figure 13 depicts the Pareto front of MMOCE and other algorithms in Case 3. It is not hard to see that MMOCE still has the merits in convergence and effectiveness.

Combined emission economic dispatch problems with wind penetration
In this paper, the test system in this section consists of 40 thermal units and a wind conversion device. The characteristic data of the 40-unit system are given in [63]. The B matrix coefficients are taken from [64], since the transmission loss is taken into consideration. The wind power penetration is set to 10%, and P is 10,500MW. As the system grows in size, the optimization becomes more difficult. In this paper, several multi-objective algorithms (i.e., MMOCE, DMOCE, FMOCE, NSGAII, MOPSO, MOO CEM) are em- Furthermore, the comparison between MMOCE and other competitor algorithms is conducted to verify the superiority of MMOCE while the results are listed in Table 20. Though it does not demonstrate evident advantage over MOO CEM with 836.1867 $/h in cost and 0.7107 t/h in emission in terms of performance, a satisfactory result obtained by MMOCE is better than other algorithms.

Combined Emission Economic Dispatch Problems with Wind Penetration
In this paper, the test system in this section consists of 40 thermal units and a wind conversion device. The characteristic data of the 40-unit system are given in [64]. The B matrix coefficients are taken from [63], since the transmission loss is taken into consideration. The wind power penetration is set to 10%, and P d is 10,500MW. As the system grows in size, the optimization becomes more difficult. In this paper, several multi-objective algorithms (i.e., MMOCE, DMOCE, FMOCE, NSGAII, MOPSO, MOO CEM) are employed to optimize this system. Figure 14 presents the optimal Pareto solutions obtained by these multi-objective algorithms, and the compromise solution obtained by MMOCE is listed in Table 21. Table 22 shows the comparison performance between the compromise solutions obtained by these algorithms.  From Figure 14, it is shown that the proposed MMOCE is superior to all other competitor algorithms in terms of the convergence. Table 22 lists their compromise solutions, and it is clearly shown that the solution obtained by MMOCE represents a cost of 122,018.6133 $/h and emission of 207,075.3012 Kg/h, which is better than other solutions obtained by other algorithms. It implies that when it is applied to large scale optimization systems, MMOCE is still able to exhibit good convergence and has the potential to provide a competitive solution for the decision maker. Compared with other multi-objective cross entropy methods, the proposed MMOCE dramatically improves the global search ability for large scale systems and exhibits excellent competitiveness.
To examine the convergence performance of the proposed method, the convergence trends for the minimum cost and minimum emissions of MMOCE and other algorithms are shown in Figure 15. The max function evaluation is set to 20,000. According to Figure 15a, the minimum cost of MMOCE decreases quickly in the first 410 iterations, and it tends to be stable after 800. In Figure 15b, the minimum emission of MMOCE decreases quickly in the first 180 iterations. The emission tends to be stable after 200 iterations. Based on Figure 15a,b, MMOCE converges faster than other methods.  From Figure 14, it is shown that the proposed MMOCE is superior to all other competitor algorithms in terms of the convergence. Table 22 lists their compromise solutions, and it is clearly shown that the solution obtained by MMOCE represents a cost of 122,018.6133 $/h and emission of 207,075.3012 Kg/h, which is better than other solutions obtained by other algorithms. It implies that when it is applied to large scale optimization systems, MMOCE is still able to exhibit good convergence and has the potential to provide a competitive solution for the decision maker. Compared with other multi-objective cross entropy methods, the proposed MMOCE dramatically improves the global search ability for large scale systems and exhibits excellent competitiveness.
To examine the convergence performance of the proposed method, the convergence trends for the minimum cost and minimum emissions of MMOCE and other algorithms are shown in Figure 15. The max function evaluation is set to 20,000. According to Figure  15(a), the minimum cost of MMOCE decreases quickly in the first 410 iterations, and it tends to be stable after 800. In Figure 15

Conclusion
In this paper, a modified multi-objective cross entropy algorithm is proposed based on the conventional cross entropy method. The drawbacks of the conventional CE method have been effectively tackled, and the optimal performance is improved to a great extent.

Conclusions
In this paper, a modified multi-objective cross entropy algorithm is proposed based on the conventional cross entropy method. The drawbacks of the conventional CE method have been effectively tackled, and the optimal performance is improved to a great extent. To evaluate its performance, the proposed method is tested on eight well-known benchmark functions. The numerical results are compared with the different heuristic optimization algorithms, which confirms the performance of the proposed method in terms of the convergence and diversity. In addition, the proposed algorithm is used to optimize three different EED problems with renewable energy sources.
The first two models (i.e., CHEED system and the modified IEEE 30-bus and sixgenerator system) consider the availability of wind and solar power, and the comparisons with other optimization methods reveal that the proposed method significantly outperformed other methods in all cases, confirming that the proposed method is capable of enhancing the convergence and seeking the optimal global solution. The third 40-unit system uses wind power as negative loads, and the numerical results also exhibit better convergence characteristics for large-scale test systems. In summary, the test results confirm the feasibility and superiority of the proposed MMOCE method.
The future work will include further simplification of the parameters used in the conventional CE method and improvement of the scalability of the method for optimizing multi-objective engineering optimization problems. Additionally, in this paper, three different types of EED model are tested. However, the renewable sources of the first model and the third model are deterministic and are only regarded as negative loads. This handling approach for renewable sources is simple and neglects the stochastic nature of the renewable resources. The future work will explore EED model considering the stochastic characteristics of various intermittent energy sources.

Data Availability Statement:
No new data were created or analyzed in this study. Data sharing is not applicable to this article.

Acknowledgments:
There is no any support which is not covered by the author contribution or funding sections.

Conflicts of Interest:
The authors declare no conflict of interest.

Nomenclature
The following abbreviations are used in this manuscript: C total Total production cost E total Total emission C i (P i ) Production cost of thermal power unit i C j P j , H j Production cost of cogeneration unit j C k (H k ) Production cost of heat-only unit k E pi (P i ) Emission  Cost coefficients for thermal power unit i m i , n i , l i , x i , y i , z i Cost coefficients for cogeneration unit j α k , β k , γ k Cost coefficients for heat-only unit k δ i , i , ξ i Emission coefficients for thermal power unit i µ j Emission coefficients for cogeneration unit j σ k Emission coefficients for heat-only unit k N TU Total number of thermal power units P TU Power produced by generating unit j α j , β j , γ j Cost coefficients for thermal power unit j d j , e j Valve effect coefficients for thermal power unit j P TU min Minimum power of thermal power unit j h i Cost constant of wind power unit i P sw,i Scheduled wind power from unit i K RC,i Reverse cost coefficient of wind unit i P wa,i Actual available power of wind unit i f wp (P w,i ) Wind power probability function of wind unit i K PC,i Penalty cost coefficient of wind unit i P rw,i Rated output power of wind unit i V out , V r , V in Cut out wind speed, nominal wind speed, and cut in wind speed P rw Rated power v Current wind speed g i Cost constant of solar power unit j P ss,j Scheduled solar power from unit j K SR,j Reverse cost coefficient of solar power unit j P sa,j Actual available power of solar power unit j K PS,j Penalty cost coefficient of solar power unit j p i , q i , r i , s i , t i Emission coefficients for thermal power unit i