Multi-Objective Environmental Economic Dispatch of an Electricity System Considering Integrated Natural Gas Units and Variable Renewable Energy Sources

: This paper presents a multi-objective economic-environmental dispatch (MOEED) model for integrated thermal, natural gas, and renewable energy systems considering both pollutant emission levels and total fuel or generation cost aspects. Two cases are carried out with the IEEE 30-bus system by replacing thermal generation units into natural gas units to minimize the amount of toxin emission and fuel cost. Equality, inequality like active, reactive powers, prohibited operating zones (POZs) which represents poor operation in the generation cost function, and security constraints are considered as system constraints. Natural gas units (NGUs) are modeled in detail. Therefore, the ﬂow velocity of gas and pressure pipelines are also considered as system constraints. Multi-objective optimization algorithms, namely multi-objective Harris hawks optimization (MOHHO) and multi-objective ﬂower pollination algorithm (MOFPA) are employed to ﬁnd Pareto optimal solutions of fuel or generation cost and emission together. Furthermore, the technique for order preference by similarity to ideal solution (TOPSIS) is proposed to obtain the best value of Pareto optimal solutions. Three scenarios are investigated to validate the e ﬀ ectiveness of the proposed model applied to the IEEE 30-bus system with the integration of variable renewable energy sources (VRESs) and natural gas units. The results obtained from Scenario III with NGUs installed instead of two thermal units reveal that the economic dispatching approach presented in this work can greatly minimize emission levels as 0.421 t / h and achieve lower fuel cost as 796.35 $ / h. Finally, the results obtained show that the MOHHO outperforms the MOFPA in solving the MOEED problem.

* 1 denotes fuel cost, 2 denotes emission, 3 denotes power loss, and 4 denotes voltage stability. Table 2 presents some of the recent research works that investigated the presence of VRESs in the classical MOEED problem. Wang et al. [21] presented a multi-objective cross-entropy algorithm based on decomposition (MOCE/D) to solve the MOEED problem with wind/hydro/photovoltaic units incorporated into the studied power system, taking into account operational constraints and uncertainties of the VRESs considered, in which POZs limits were considered in the problem while the valve point effects were not included. Chen et al. [22] presented the multi-objective population extremal optimization (MOPEO) technique to minimize emissions and costs in the IEEE 30-bus system integrated with thermal, solar, and wind generation units. Bora et al. [23] presented the NSGA-II incorporated by a reinforcement learning method, called NSGA-RL for solving the MOEED problem. A formulation of six thermal generation units integrated with the wind energy system was presented to optimize fuel costs and emissions. The results explored that the NSGA-RL technique can solve multi-objective EED problems effectively. But it would have been better to use more than one renewable source like solar or hydropower. Biswas et al. [24] presented MOEA/D and summation based multi-objective differential evolution (SMODE) to minimize emissions and costs in the IEEE 30-bus system integrated with stochastic solar, wind, and hydropower generation units to solve the MOEED problem with a limited number of thermal plants. Yin et al. [25] presented a dynamic day-ahead stochastic scheduling of thermal/wind/photovoltaic (PV)/hydropower systems, but only the fuel cost was considered in the formulated optimization problem. Li et al. [26] presented a multi-objective moth-flame optimization (MOMFO) technique for solving dynamic EED of tradable green certificates-based hybrid renewable energy systems. However, the non-linear constraints arising from the power flow were not considered in the formulated optimization problem. Elattar [27] presented a modified shuffle frog leaping algorithm (MSFLA) to minimize fuel cost and emissions in a combined heat and power MOEED, taking into account the presence of wind and solar power. However, VRESs uncertainties were not considered. Also, arbitrary weightings have been used to turn the problem of multi-objective optimization into a single objective optimization problem. Joshi and Verma [28] presented a hybrid PSO and teaching learning-based optimization (TLBO-PSO) to minimize the total generation costs and emissions of a system integrated with thermal, solar, and wind generators. Different test cases were presented to investigate the impacts of using VRESs on the performance of the system. Also, VRESs uncertainties were not considered.
In this study, the standard IEEE 30-bus system is adapted with a limited number of thermal plants to engage solar PV, NG, wind, and small-hydro generation units. The stochastic nature of renewable sources like small-hydro power, wind, and solar are explored in detail, utilizing Gumbel, Weibull, and lognormal probability density functions (PDFs), respectively. Due to the uncertainty and intermittency of the VRESs, penalty cost for underestimation and reserve cost for overestimation is included in the proposed cost model. Then, a constrained multi-objective optimization problem is formulated to minimize the EED problem while complying with the power flow equality and inequality constraints, system voltage limits, NG constraints, prohibited operating zones (POZs) limits, and thermal capacity of lines. Furthermore, the sets of Pareto solutions for the multi-objective EED problem are found by using two multi-objective optimization (MOO) strategies: multi-objective Harris hawks optimization (MOHHO) and multi-objective flower pollination algorithm (MOFPA) to solve the EED problem formulated in this work. Recalling that most of the MOO techniques are usually adopted to deal with unconstrained optimization problems. Also, as the decision-maker may favor a single solution, hence, a ranking procedure named TOPSIS has been employed to obtain a single solution from the non-dominated solutions set of the problem under study. A TOPSIS metric has the advantages of consistency, simplicity, and comprehensibility in its calculation procedures.
The contributions of this work are briefly summarized as follows: • Formulation of the MOEED problem considering thermal, solar, wind, hydropower, and natural gas units.

•
Stochastic analysis of VRESs used has been presented using the most proper PDFs.

•
All system constraints like security, equality, and inequality power flow conditions, POZs limits, and NG constraints are considered in the formulated MOEED problem. • MOO techniques such as MOHHOA and MOFPA are employed to solve the MOEED problem. • A comparative analysis of the solutions obtained by the three optimization techniques is proposed. • TOPSIS is used to obtain the best compromise solution to the MOEED problem.
The rest of the paper is organized in different sections, in which configuration of the system studied is introduced in Section 2, a mathematical analysis of the EED issue incorporating VRESs and NG, formulation of the MOEED problem, system constraints, and the two optimization techniques used in work are introduced in Section 3, illustrations obtained are presented and explained in Section 4, and, lastly, conclusions and future studies are introduced in Section 5.

System Studied
The main goal of the power utilities is to schedule the output of the generators in order to fulfill the demand needs by reducing fuel costs without taking emissions into consideration. Currently, to comply with the national and global standards of environmental conservation, every nation has begun to help protect the atmosphere from pollution by different strategies and new approaches to decrease pollutants. Otherwise, they will be penalized. In addition to that, in electrical networks, minimization of the fuel cost of generation units plays a crucial role in the economic dispatch (ED) problem to meet the load requirements. On the other hand, the minimization of emissions (also referred to as the environmental objective) is considered to be one of the most significant aspect for alleviating the problems of climate change and environmental pollution. In the case of the conventional ED problem, only the economic objective is regarded and conceived as a single-objective optimization issue. However, because of the alleviation of global warming, more consideration has been devoted to controlling the emission of greenhouse gases. Therefore, it is of considerable significance to tackle the question of economic-emission dispatch (EED) by balancing the two goals at the same time in terms of economic and environmental aspects. The key enabler to do so is to integrate clean and economic VRESs into electrical networks to minimize environmental emissions. An adopted IEEE standard 30-bus system, that comprises thermal generators and renewable energy sources, is investigated in this work. As shown in Figure 1, the model comprises three thermal power generation units (TPGUs) installed on buses 1, 2, and 8, in addition to three different VRESs, namely wind, PV, and a hybrid power generation system of PV and small hydropower (PVSH) installed on buses 5, 11, and 13, respectively. The main parameters of the investigated system are given in Table 3.  Three different scenarios are tested to realize the minimization of both emission and fuel or generation cost. Generally, each scenario has three VRESs of wind, PV, and PV-small hydro units at buses 5, 11, and 13, respectively. It can be explained as follows: Scenario I: Using three TPGUs and three VRESs [24]. Scenario II: Replacing the fuel of TPGU at bus 1 into NGU. Scenario III: Replacing the fuel of TPGUs at buses 2 and 8 into NGUs.
Two optimization techniques of MOHHO and MOFPA are applied to the problem under study. The TOPSIS performance indicator is utilized to rank Pareto fronts (PFs) obtained from multiobjective optimization algorithms. The procedure of the proposed techniques is illustrated in Figure  2.  Three different scenarios are tested to realize the minimization of both emission and fuel or generation cost. Generally, each scenario has three VRESs of wind, PV, and PV-small hydro units at buses 5, 11, and 13, respectively. It can be explained as follows: Scenario I: Using three TPGUs and three VRESs [24]. Scenario II: Replacing the fuel of TPGU at bus 1 into NGU. Scenario III: Replacing the fuel of TPGUs at buses 2 and 8 into NGUs.
Two optimization techniques of MOHHO and MOFPA are applied to the problem under study. The TOPSIS performance indicator is utilized to rank Pareto fronts (PFs) obtained from multi-objective optimization algorithms. The procedure of the proposed techniques is illustrated in Figure 2.

Formulation of the Multi-Objective Function
In this work, three scenarios are investigated, in which the three VRESs are kept at buses 5, 11, and 13, in all the scenarios under focus. In Scenario 1, 3 TPGUs and 3 VRERs are considered. In Scenario 2, one TPGU, three VRESs, and two NGUs are considered, where the TPGU installed on buses 2 and 8 are replaced by natural gas units (NGUs), Furthermore, in Scenario 3, two TPGUs, three VRESs, and one NGU are considered, where the large TPGU installed on bus 1 is replaced by the NGU.

Objective Functions
The MOEED problem to be solved is a multi-objective problem, in which the first objective function represents the environmental impacts and the second one represents the economic aspects as expressed in Equation (1): where, and represent the total emissions and fuel costs to be minimized, respectively.

Total Fuel Costs
The total cost of the generated power is the summation of costs of the TPGUs, VRESs, and NGUs (if NGUs are considered), as expressed in Equation (2).
where, ( ) represents the total TPGUs cost, ( ) represents the total VRESs cost, and ( ) represents the total NGUs cost, where =0 with no NGUs considered and = 1 if NGUs are considered.

Fuel Cost Analysis of Thermal Power Generation Units (TPGUs)
The TPGUs cost, in dollars per megawatt-hour ($/MWh), should consider the flow of steam to the turbine blades and also the sudden changes in the valve's status. Steam in these plants is controlled by valves to run the turbine through a separate set of nozzles. That group of nozzles must be run at full output to achieve the best efficiency [29]. For the mandatory production, the valves are opened sequentially, which results in the discontinuity cost curve as shown in Figure 3. The cost function of TPGUs is given in Equation (3) [30].

Formulation of the Multi-Objective Function
In this work, three scenarios are investigated, in which the three VRESs are kept at buses 5, 11, and 13, in all the scenarios under focus. In Scenario 1, 3 TPGUs and 3 VRERs are considered. In Scenario 2, one TPGU, three VRESs, and two NGUs are considered, where the TPGU installed on buses 2 and 8 are replaced by natural gas units (NGUs), Furthermore, in Scenario 3, two TPGUs, three VRESs, and one NGU are considered, where the large TPGU installed on bus 1 is replaced by the NGU.

Objective Functions
The MOEED problem to be solved is a multi-objective problem, in which the first objective function represents the environmental impacts and the second one represents the economic aspects as expressed in Equation (1): where, E tot and C tot represent the total emissions and fuel costs to be minimized, respectively.

Total Fuel Costs
The total cost of the generated power is the summation of costs of the TPGUs, VRESs, and NGUs (if NGUs are considered), as expressed in Equation (2).
where, C tot (P TPGU ) represents the total TPGUs cost, C tot (P VRES ) represents the total VRESs cost, and C tot (P NGU ) represents the total NGUs cost, where µ = 0 with no NGUs considered and µ = 1 if NGUs are considered.

Fuel Cost Analysis of Thermal Power Generation Units (TPGUs)
The TPGUs cost, in dollars per megawatt-hour ($/MWh), should consider the flow of steam to the turbine blades and also the sudden changes in the valve's status. Steam in these plants is controlled by valves to run the turbine through a separate set of nozzles. That group of nozzles must be run at full output to achieve the best efficiency [29]. For the mandatory production, the valves are opened sequentially, which results in the discontinuity cost curve as shown in Figure 3. The cost function of TPGUs is given in Equation (3) [30].
where a T i , b T i and c T i are the cost parameters of the ith thermal generator unit (P TPGU i ). The two parameters d i and e i represent the valve point effect.P min TPGU i represents the minimum power of P TPGU i during the generator operation. These parameters are specified in Table 4. where , and are the cost parameters of the ith thermal generator unit ( ). The two parameters and represent the valve point effect.
represents the minimum power of during the generator operation. These parameters are specified in Table 4.   (4): However, each renewable source has a specific cost function. Also, the amount of underdelivered or over-delivered power might be calculated based on the PDFs of each renewable source. To cope with the intermittency nature of the VRESs, first, standby power generation units (SPGUs) might be installed when the generated power is less than the scheduled power. Second, energy storage (ES) units might be installed to reserve the extra-planned generated power.
Weibull, lognormal, and Gumbel distributions are applied to fit the random data of the wind speed, solar irradiance, and the flow rate of the small hydro unit, respectively, to formulate the cost expressions as illustrated in detail in the following subsections.
A. Cost calculation of wind plants  The VRESs cost, in dollars per megawatt-hour ($/MWh), is the sum of the total costs of WTs (C tot (P WT )), PVs (C tot (P PV )), and hybrid PV and SHP (C tot (P PVSH )) which can be expressed as in Equation (4): However, each renewable source has a specific cost function. Also, the amount of under-delivered or over-delivered power might be calculated based on the PDFs of each renewable source. To cope with the intermittency nature of the VRESs, first, standby power generation units (SPGUs) might be installed when the generated power is less than the scheduled power. Second, energy storage (ES) units might be installed to reserve the extra-planned generated power.
Weibull, lognormal, and Gumbel distributions are applied to fit the random data of the wind speed, solar irradiance, and the flow rate of the small hydro unit, respectively, to formulate the cost expressions as illustrated in detail in the following subsections.

A. Cost Calculation of Wind Plants
The cost equation of the WTs considers the direct investment costs, in addition to the costs of the standby and ES units.
The direct cost C d WT (P W sch ) of WTs represents the initial, operation, and maintenance costs, and is expressed in Equation (5).
where K d WT represents the direct cost parameter and P W sch represents the scheduled power of the WTs. when the actual power from the turbines is less than the scheduled power, the system involves possible standby units (reserve capacity) to maintain the demand requirements. The cost of the reserve capacity (C r WT ) is expressed in Equation (6).
K r WT represents the cost parameter of the standby units and P W act represents the actual delivered power from the WTs. Similarly, when P W act is greater than P W sch , the cost of the storage units will present as described in Equation (7): where P W r represents the rated wind power. f w (p w ) represents the PDF of the wind speed. In fact, the actual power of the standby and storage units depends on f w (p w ). The PDF of the wind energy system usually follows the well-known Weibull distribution to fit the random frequency of each wind speed level [31,32]. Figure 4a illustrates the Weibull PDF of the wind speed data by applying 8000 Monte-Carlo scenarios considering the scale and the shape parameters of the Weibull PDF, (denoted α and β) respectively. α and β are set to 9 and 2, respectively. The probability ( f v (v)) of the wind speed (v) is expressed as in Equation (8): The cost parameters of wind plants are illustrated in Appendix A. The yielded power of wind generators that depends on the wind speed is given by Equation (9): where, v in , v r , v out denote the cut-in speed, rated speed, and cut-out speed of the WTs, respectively. The wind power probability is expressed as in Equation (10): Finally, the total cost of the wind plant is given in Equation (11):

B. Cost Calculation of the Solar Plant
Likewise, the total cost equation of the solar plant has been built based on the same philosophy used to calculate the cost equations of the wind plants.
The direct cost C d PV (P PVsch ) of PVs represents the initial, operation, and maintenance costs, and is expressed in Equation (12).
where K d PV represents the direct cost coefficient of the PV system and P PVsch denotes the scheduled power of the PV system. When the actual power provided from the PV system (P PV act ) is less than P PVsch , then it is required to implement standby units (reserve capacity) to maintain the demand requirements. The cost of the reserve capacity (C r PV ) is expressed in Equation (13).
K r PV represents the cost coefficient of the standby units for the PV system. As mentioned before, the cost of the storage units (C s PV ) will present, if P PV act is greater than P PVsch and is described as in Equation (14): The delivered power from the standby and storage units depends on the PDF of the solar irradiance (G). The PDF of G is fitted via the lognormal distribution [33,34], as illustrated in Figure 4b. Equation (15) describes the probability of G at lognormal fit parameters set as µ = 5.6 and σ = 0.6; thus: The available power from the PV generation unit can be determined as in Equation (16): where the standard solar irradiance G std is equal to 1000 W/m2. During the operation irradiance R c is set as 120 W/m 2 . P PVr represents the rated output power of the PV units. The cost parameters of PV are illustrated in Appendix A. Finally, the total PV generation cost (C tot PV ) comprises the summation of the direct cost, standby unit cost, and the storage unit cost, and can be illustrated as in Equation (17):

C. Cost Calculation of the Photovoltaic and Small Hydropower (PVSH) Plant
Gumbel distribution [35] is applied to fit the river flow data, as shown in Figure 4c, in which the river flow rate probabilistic Q w that follows the Gumbel distribution with locational factor λ and the scaling factor γ is expressed as follows: The output power from the hydro plant P H (Q w ) mainly depends on Q w as given in Equation (19): where η w , ρ w , g w , and H w represent the hydro turbine efficiency, the water density, the gravity acceleration, and the effective pressure head, respectively. Estimated values of these coefficients shall be taken into account for the determination of P H (Q w ) are η W = 0.86; ρ W = 1000 kg/m 3 ; g W = 9.81 m/s 2 ; and H w = 26 m. In this work, the hydro plant at λ = 15 and γ = 1.2 is incorporated with a PV plant to enhance the hydro plant performance with the combined energy mixed-generation system. The cost equation of the PVSH system C tot PVSH is formulated in the same manner as the WT and PV cost equations, as in Equation (20); thus: where P PVSHsch and P PVSH act express the scheduled power and the actual power from the hybrid unit, respectively. C dPVSH (P PVSH ) represents the direct cost of the PVSH. C r PVSH and C S PVSH represent the costs of the standby and storage units, respectively. The cost parameters of PVSH are illustrated in Appendix A. where and express the scheduled power and the actual power from the hybrid unit, respectively.
( ) represents the direct cost of the PVSH. and represent the costs of the standby and storage units, respectively. The cost parameters of PVSH are illustrated in Appendix A.

Fuel Cost Analysis of Natural Gas Units (NGUs)
In this work, natural gas units (NGUs) are employed instead of thermal power units to attain both fuel cost and pollutant emission levels as low as possible. The total cost of NGUs consists of many expenditure parts; initial cost, operation, and maintenance cost, and fuel cost. Thus, the total cost of the NGUs ( ) is given in Equation (21) [36,37], thus: where , and are the cost parameters of the ith NGU ( ), and g denotes the coefficient of the initial and operating costs of NGUs. These parameters are specified in Table 5. To sum up, of the overall system is given in Equation (22). Fuel Cost Analysis of Natural Gas Units (NGUs) In this work, natural gas units (NGUs) are employed instead of thermal power units to attain both fuel cost and pollutant emission levels as low as possible. The total cost of NGUs consists of many expenditure parts; initial cost, operation, and maintenance cost, and fuel cost. Thus, the total cost of the NGUs ( C tot (P NGU )) is given in Equation (21) [36,37], thus: where a N i , b N i and c N i are the cost parameters of the ith NGU (P NGU i ), and g NGU i denotes the coefficient of the initial and operating costs of NGUs. These parameters are specified in Table 5. To sum up, C tot of the overall system is given in Equation (22).

Emission Levels
The emissions of thermal and natural gas units are only presented since the VRES has no emissions. The total emission (E tot ) due to these units is described in Equation (23): where E TPGUi and E NGUi are the total emissions of TPGU and NGU, respectively.

Emission Analysis of TPGUs
The total emission E tot (TPGU) is defined as the amount of emission of harmful gases negatively affecting the environment as SOx and NOx as a function in the generated output power. Emission in tones per hour (t/h) can be calculated as in Equation (24): where ϕ T i , ψ T i , ω T i ,τ T i and ξ T i are the parameters of emission levels related to the ith TPGU. These parameters are specified in Table 4.

Emission Analysis of NGUs
The total emissions due to NGUs is the same as the total emission due to TPGUs; hence, it is described as in Equation (25) [39]: (25) where ϕ N i , ψ N i , and ω N i are the parameters of emission levels related to the ith NGU. These parameters are specified in Table 5.

Constraints
The main constraints should be considered during the solving of the multi-objective function of any configuration that is summarized as follows:

Power Balance Constraints
Power balance limitations, active and reactive powers, are the summation of the power consumed from the total loads and power losses in the system network.

Active and Reactive Powers Limits
Equations (28) to (32) express the active power operational limits of all used generators; thermal, wind, solar, and natural gas units, respectively.
Also, the reactive power operational limits of all generator units are considered as in Equations (33) to (37): Prohibited Operating Zones Limits Because of some physical limitations of thermal generators like vibrations in a shaft bearing or failures in pumps, boilers, etc., POZs are allowed in certain operating regions, which, in turn, lead to a discontinuous operation in the thermal generation units. POZs can be described in Equation (38): where P minPOZ,j TPGU i and P maxPOZ,j TPGU i are the minimum and maximum boundaries in megawatt of the jth POZ of the ith thermal generator units.

Security Constraints
The generators' bus voltage security limits and the load bus voltages limits are illustrated in Equations (39) and (40), respectively. The branches' capacity limits are described in Equation (41): where V G i represents the voltage of the ith on generator bus and V L p represents the voltage of the pth on the load bus. N G , N L , and nl represent the number of generator buses, the number of load buses, and the branches number in the network, respectively.
In addition, the other two parameters of the network power loss (P loss ) and voltage quality indicator which called voltage deviation (VD) are also considered as network restrictions and can be calculated as in Equations (42) and (43): where G q(ij) denotes the transconductance of branch q connected to bus i and bus j. δ ij represents the phase difference between δ i and δ j of the buses i and j.

Natural Gas Limits
When it is planned to feed a power plant from the natural gas main network, an extension for the network should be executed. This extension has been serviced with a pressure reduction station for natural gas. Many operational conditions should be considered for each scenario such as the input pressure, flow rate, and input temperature. Usually, the natural gas pipeline carries a pressure higher than is required at the loads to ensure the reliability of the natural gas source. Therefore, a pressure regulator should be applied to keep the pressure within the proper limits. This regulator mainly consists of two chambers, one for the inlet (high) pressure and the other for the outlet (low) pressure. Its control is designed based on the closed-loop, where feedback is applied from the output to the pilot (control) unit to modify or correct its actions as shown in Figure 5.
where ( ) denotes the transconductance of branch q connected to bus i and bus j. represents the phase difference between and of the buses i and j.

Natural Gas Limits
When it is planned to feed a power plant from the natural gas main network, an extension for the network should be executed. This extension has been serviced with a pressure reduction station for natural gas. Many operational conditions should be considered for each scenario such as the input pressure, flow rate, and input temperature. Usually, the natural gas pipeline carries a pressure higher than is required at the loads to ensure the reliability of the natural gas source. Therefore, a pressure regulator should be applied to keep the pressure within the proper limits. This regulator mainly consists of two chambers, one for the inlet (high) pressure and the other for the outlet (low) pressure. Its control is designed based on the closed-loop, where feedback is applied from the output to the pilot (control) unit to modify or correct its actions as shown in Figure 5. In other words, the influence of the required NG volume on the remind loads of the NG supplying network should be studied [40]. The required volume from natural gas to generate a certain electrical power is described in Equation (44): where denotes the efficiency of the gas turbine, and represents the high heat value of natural gas. The length, diameter, material type, working pressure, and average flow rate of the natural gas ( , ) have the key role before supplying the natural gas into a gas turbine. Mathematically, it is described in Equation (45): In other words, the influence of the required NG volume ν NGU on the remind loads of the NG supplying network should be studied [40]. The required volume from natural gas to generate a certain electrical power is described in Equation (44): where η NGU i denotes the efficiency of the gas turbine, and HHV represents the high heat value of natural gas. The length, diameter, material type, working pressure, and average flow rate of the natural gas ( q xy,t ) have the key role before supplying the natural gas into a gas turbine. Mathematically, it is described in Equation (45): where C NG represents the constant-coefficient, T o represents the standard temperature (288 K), p o denotes the absolute pressure at atmospheric conditions (bar), p x , p y represent the absolute upstream (inlet) pressure and the absolute downstream (outlet) pressure, respectively. D P represents the internal diameter of the pipe in millimeters, S NG is the specific gravity of natural gas, Z represents the coefficient of difference between the gas in actual condition and in an ideal condition (that is, the compressor factor), and T av represents the mean temperature of the flowing gas. L P , f denote the pipeline length in meter and the hydraulic friction factor which is a range from 0.009 to 0.015 for corrugated Polyethylene pipes with smooth inner walls, respectively. If the flow velocity exceeds 20 m/s, the dust particles will change position, where the motion of the particles had a bad effect on the cooking appliances, pressure controllers, and it may lead to corrosion to the inner surface of the pipeline. Thus, the flow velocity (U NG ) in m/sec is expressed in Equation (47): The flow velocity is inversely proportional to the pipeline pressure in which both of the flow velocity U NG and pipeline pressure p p should be considered as system constraints [39,40] throughout the replacement of the thermal to gas generation units.
Equations (48) and (49) represent the equality constraints of line pack of pipeline x−y at hour t (L x−y,t ). Equation (50) represents the initial and final values of the line pack are equal. Also, Equations (51) to (54) represent the inequality constraints of flow velocity, pipeline pressure, the flow rate of NG suppliers, and the air compressor constraints, respectively.

Multi-Objective Optimization Techniques
Firstly, the offspring population is produced using the parent population, then a combination of the old and off-spring populations to form the total population is undertaken [41]. A non-dominated criterion is utilized to sort the total population. Secondly, the new population is then composed of diverse non-dominated fronts, in which the best non-dominated fronts are occupied. Then, the filling process continues with solutions of the second non-dominated front, then the third, and so on, as illustrated in Figure 6. The fronts that could not be accommodated are canceled. A niching (crowding) strategy is employed to choose members from the last front instead of discarding some members arbitrarily from the last front, which reside in the least crowded region in the front. The algorithm guarantees that the crowding strategy will be able to select a diverse set of solutions. Finally, the continuation of this algorithm will promise a better spread among the non-dominated solutions when the whole population converges to the Pareto-optimal front [43].
Crowded distance is the mean distance between two solutions along with each of the objectives on either side of a particular solution. The crowded distance calculation is illustrated in Figure 7 and the following steps are utilized to determine the crowded distance of every solution in the Fr set [44].
Step 1: Solutions are arranged in each objective domain.
Step 2: Crowded distances of the first solution and the last solution in the rank are selected as to infinity.
Step 3: For each of the other solutions, the crowded distance will be calculated in Equation (55): where M, i, and j represent the objectives number, the number of the solution, and the total number of solutions in the set Fr, respectively. and represent the m th objective functions of solution number (i+1) and (i-1) in the set Fr, respectively. and represent the minimum and maximum values of the m th objective function in the set Fr, respectively.  The fronts that could not be accommodated are canceled. A niching (crowding) strategy is employed to choose members from the last front instead of discarding some members arbitrarily from the last front, which reside in the least crowded region in the front. The algorithm guarantees that the crowding strategy will be able to select a diverse set of solutions. Finally, the continuation of this algorithm will promise a better spread among the non-dominated solutions when the whole population converges to the Pareto-optimal front [43].
Crowded distance is the mean distance between two solutions along with each of the objectives on either side of a particular solution. The crowded distance calculation is illustrated in Figure 7 and the following steps are utilized to determine the crowded distance of every solution in the F r set [44].
Step 1: Solutions are arranged in each objective domain.
Step 2: Crowded distances of the first solution and the last solution in the rank are selected as to infinity.
Step 3: For each of the other solutions, the crowded distance will be calculated in Equation (55): The fronts that could not be accommodated are canceled. A niching (crowding) strategy is employed to choose members from the last front instead of discarding some members arbitrarily from the last front, which reside in the least crowded region in the front. The algorithm guarantees that the crowding strategy will be able to select a diverse set of solutions. Finally, the continuation of this algorithm will promise a better spread among the non-dominated solutions when the whole population converges to the Pareto-optimal front [43].
Crowded distance is the mean distance between two solutions along with each of the objectives on either side of a particular solution. The crowded distance calculation is illustrated in Figure 7 and the following steps are utilized to determine the crowded distance of every solution in the Fr set [44].
Step 1: Solutions are arranged in each objective domain.
Step 2: Crowded distances of the first solution and the last solution in the rank are selected as to infinity.
Step 3: For each of the other solutions, the crowded distance will be calculated in Equation (55): where M, i, and j represent the objectives number, the number of the solution, and the total number of solutions in the set Fr, respectively. and represent the m th objective functions of solution number (i+1) and (i-1) in the set Fr, respectively. and represent the minimum and maximum values of the m th objective function in the set Fr, respectively.

Multi-Objective Flower Pollination Algorithm (MOFPA)
The flower pollination algorithm (FPA) settled by Yang [45] is a metaheuristic optimization method inspired by the nature-based flower pollination technique. Ultimately, the principal function of a flower is to replicate by pollination. Flower pollination is usually associating with pollen transmission and is also linked to pollinators including birds and insects. Indeed, some flowers and insects have a rather unique relationship with flower-pollinators, since certain flowers may attract only a certain type of insect or bird for efficient pollination. Pollination appears in two main forms: abiotic and biotic. Approximately 90% of flowering plants depend on the process of biotic pollination, in which pollinators transfer the pollen. Approximately 10% of pollination follows an abiotic process that needs no pollinators. Wind and diffusion aid in the process of pollination of these flowering plants.
Pollination can be divided into self-pollination and cross-pollination. Self-pollination is one flower's pollination from the pollen of the same flower. Cross-pollination is the pollination from the pollen of a flower of different plants. The goal of flower pollination is the existence of both the fittest and the optimal reproduction of plants in terms of both numbers and fittest. This can be known as plant species optimization method. All these variables and flower pollination processes produced the optimal reproduction of the flowering plants. The following four principles provide guidelines concerning the method of pollination used and the selection of step size [46]. (1) Global pollination is implemented when a uniform production arbitrary value rand ≤ ρ, illustrated in Equation (56), to sequentially change the location of the ith flower X i utilizing its spacing from the best flower X best . X t+1 where Lev(λ) denotes the Levy flight attitudes of the pollinators. This follows that the distribution of Levy reflects the intensity of the pollination. To control the size of the step, a scaling factor (F) is chosen.
(2) In the FPA technique, the local pollination will be performed to use a uniform distribution random value (rand i ) that lies from 0 to 1 to regulate the mutation in the ith flower.
To make the FPA able to solve the MOO problems, non-dominating sorting and crowding distance procedures are combined with the FPA to obtain a distributed PF. Initially, the population of the flowers is randomly created and the objective function is assessed. The new population is obtained by performing the local and global pollination steps of the basic FPA and the new solution is obtained by evaluating the objective function. To make the MOFPA algorithm robust no repository is used for saving the previous solution. The preceding solution and the new solutions are merged and sorted accordingly, and then it is truncated to the initial size to reduce computation time. These steps are repeated until the max iteration is attained. The flowchart of the MOFPA algorithm is provided in Figure 8. Moreover, the pseudo-code of MOFPA is illustrated as illustrates in Algorithm 1.

Multi-Objective Harris Hawks Optimization (MOHHO)
Heidari et al. [47] presented a technique of optimization called Harris hawks optimization (HHO). HHO is a nature-based optimization method inspired by the behavioral analysis of Harris Hawk birds. The spirit of the strategy is the collaboration between the hawks in the search for prey, in which Harris hawk's family strikes the prey from various angles to catch it by surprise. Apparently, the escape activity of the prey is related to the Harris hawk chase pattern. Birds are participating in the phase of the attack. While, the Harris hawks' leader reaches the expected target, records it, then falls out of control, and the next hawk starts hunting. This technique fatigues the target and ultimately ends in its capture. HHO, which is a global optimizer, will retain the equilibrium between the processes of development and discovery.
There are three steps to the HHO algorithm. The first step is the discovery function, which is described as follows: where X (t) , X (t + 1) , and X prey(t) represent the current location of a hawk, the location of the hawk in the following iteration t, and the location of the victim, respectively. r 1 , r 2 , r 3 , r 4 and q represent random values between (0,1). X a(t) and X rand(t) denote the average location of Harris Hawk and the random selection hawk among the population, respectively. In addition, LB and UB represent the lower and upper bounds, respectively. X a(t) can be formulated as follows: where X i(t) illustrates the location of each Harris hawk in iteration t and N denotes the Harris hawks' numbers. Diversification is the second process. The strength of the hawks is fishing and shooting. The prey energy can be expressed as: E o denotes the energy of the first point, T represents the maximum iterations number and E is the energy of escape. r 1 represents the random value from 0 to 1. At this point, when |E o | ≥ 1 diversification occurs, and when |E o | < 1 intensification happens.
The third phase is intensification, which is primarily aimed at enhancing local solutions from solutions previously obtained. This phase is a shocking attack by the hawks on the prey known in the previous phase. Based on the escape of the prey and the chasing of the hawks, 4 models were presented for the attack phase. The stages of the HHO technique are displayed in Figure 9. In addition, the method of implementing the proposed MOHHO is illustrated in Figure 10.
following iteration t, and the location of the victim, respectively. r1, r2, r3, r4 and q represent random values between (0,1). Xa(t) and Xrand(t) denote the average location of Harris Hawk and the random selection hawk among the population, respectively. In addition, LB and UB represent the lower and upper bounds, respectively. Xa(t) can be formulated as follows: where Xi(t) illustrates the location of each Harris hawk in iteration t and N denotes the Harris hawks' numbers. Diversification is the second process. The strength of the hawks is fishing and shooting. The prey energy can be expressed as: denotes the energy of the first point, T represents the maximum iterations number and E is the energy of escape.
represents the random value from 0 to 1. At this point, when | | ≥ 1 diversification occurs, and when | | < 1 intensification happens.
The third phase is intensification, which is primarily aimed at enhancing local solutions from solutions previously obtained. This phase is a shocking attack by the hawks on the prey known in the previous phase. Based on the escape of the prey and the chasing of the hawks, 4 models were presented for the attack phase. The stages of the HHO technique are displayed in Figure 9. In addition, the method of implementing the proposed MOHHO is illustrated in Figure 10.

A. Soft Besiege
The condition is available whenever ≥ 0 and | | ≥ 0, which is expressed in Equation (65): where Δx represents the gap between the actual position and the Hawk prey position in the next iteration. J represents the arbitrary jump power of the victim although it escapes, which corresponds to J = 2(1 − r5) and r5 denotes a random value between 0 and 1.

B. Hard Besiege
The condition of the hard besiege is available if r ≥ 0 and | |≺ 0. The victim is exhausted and does not have enough strength to run. This process shall be formulated as follows:

C. Soft Besiege with Progressive Rapid Dive
The condition of soft besiege with progressive quick dive is available if r ≺ 0 and | | ≥ 0. In this scenario, the victim has enough resources to flee. At this point, the hawk checks the next step to execute soft besiege measures, which can be articulated as follows:

A. Soft Besiege
The condition is available whenever r ≥ 0 and |E o | ≥ 0, which is expressed in Equation (65): where ∆x represents the gap between the actual position and the Hawk prey position in the next iteration. J represents the arbitrary jump power of the victim although it escapes, which corresponds to J = 2(1 − r 5 ) and r 5 denotes a random value between 0 and 1.

B. Hard Besiege
The condition of the hard besiege is available if r ≥ 0 and |E|≺ 0. The victim is exhausted and does not have enough strength to run. This process shall be formulated as follows:

C. Soft Besiege with Progressive Rapid Dive
The condition of soft besiege with progressive quick dive is available if r ≺ 0 and |E| ≥ 0. In this scenario, the victim has enough resources to flee. At this point, the hawk checks the next step to execute soft besiege measures, which can be articulated as follows: where D denotes the dimensional point and S denotes a vector by size 1 × D randomly and LF represents the function of the levy flight [19]. As a consequence, we have: 1.
Hard besiege with progressive quick dive.

2.
Hard besiege with progressive quick dive is suitable if r≺ 0 and |E| ≺ 0. In that event, the victim does not have sufficient energy to escape properly. This case is expressed in Equation (71): Moreover, the pseudo-code of MOHHO is presented as illustrates in Algorithm 2.
where D denotes the dimensional point and S denotes a vector by size 1×D randomly and LF represents the function of the levy flight [19]. As a consequence, we have: 1. Hard besiege with progressive quick dive.
2. Hard besiege with progressive quick dive is suitable if r≺ 0 and | |≺ 0. In that event, the victim does not have sufficient energy to escape properly. This case is expressed in Equation (71): Moreover, the pseudo-code of MOHHO is presented as illustrates in Algorithm 2.

A Technique for Order Preference by Similarity to Ideal Solution (TOPSIS)
To make a comparison between Pareto and optimization techniques to expedite and combine a range of alternatives, only one alternative should be favored through the decision-maker. Ranking procedures may be applied to get a set of non-dominating solutions to a unified solution. In this work, a ranked method named TOPSIS has been utilized for resolving this variation of multi-attribute decision making (MADM) [49]. TOPSIS tends to recognize that the best solution must involve the shortest length from the positive-ideal solution and the furthest length from the negative-ideal solution [50]. The primary goals of utilizing TOPSIS are the cohesive, comprehensible, and effortlessness of its calculations. In this method, the positive-ideal solution formed from all best attributes and negative-ideal solution formed from all worst attributes are obtained. TOPSIS works based on the calculation of Euclidean distance to the ideal alternative. TOPSIS method has been employed to rank the specified solutions of the Pareto front attained by the optimizations used. The main principle of TOPSIS relies on finding a solution that must have the shortest length from the positive ideal solution (v + ) and the furthest length from the negative ideal solution (v − ). In this work, the positive ideal solution (v + ij ) has the largest maximum index alteration and the smallest dispersal, and the negative ideal solution (v − ij ) is the opposite solution. The TOPSIS procedure is summarized by the following steps: Step 1: Determine a decision matrix S (composed of a set of non-dominant solutions of the Pareto solutions). The value G ij is an indication for the performance rating of the ith alternative concerning the jth function. Let, W = (w 1 , w 2 ) be the relative weight vector about the objectives, satisfying n j=1 w j = 1.
Step 2: Determine the normalized value Y ij by applying Equation (72), i.e., by normalizing the decision matrix: Step 3: Determine V ij using Equation (73) that denotes the weighted normalized decision matrix.
Step 4: Get v + ij and v − ij using Equations (74) and (75): Step 5: Using the n-dimensional Euclidean distance, determine the separation measures by Equations (76) and (77): Step 6: Determine the relative closeness (RC) or performance index to the ideal solution, which has been formulated as in Equation (78): Step 7: The preference order is to be ranked, so that the best compromise solution is considered as the solution with the greatest RC to the ideal solution.

Simulation Results and Comparative Analysis
In this section, three different scenarios are tested in order to realize the minimization of both emission and fuel or generation cost. Generally, each scenario has three VRESs of wind, PV, and PV-small hydro units at buses 5, 11, and 13, respectively. It can be explained as the following: Scenario I: Using three TPGUs and three VRESs [24]. Scenario II: Replacing the fuel of TPGU at bus 1 into NGU. Scenario III: Replacing the fuel of TPGUs at buses 2 and 8 into NGUs.
Two optimization techniques of MOHHO and MOFPA are applied to the problem under study. The TOPSIS performance indicator is utilized to rank Pareto fronts (PFs) obtained from multi-objective optimization algorithms.

First Scenario
In this scenario, a comparison of best Pareto fronts (PFs) obtained by the MOEA/D and SMODE was implemented as reported in [24] and illustrated in Table 6. It can be noted from that study, the diversity of SMODE is better than MOEA/D especially in the direction of cost objective. In other words, SMODE achieves fewer emission levels to the value of 0.4721 t/h, while MOEA/D achieves marginally lower fuel cost value 919.040 $/h. This event will be included in the next scenarios to obtain the PFs of minimum cost and emission alike.

Second Scenario
In this scenario, we replace the biggest thermal unit at bus 1 to NGU incorporating stochastic VRESs. Figure 11 displays the best PFs of MOHHO and MOFPA techniques with the TOPSIS indicator. Moreover, Table 7 illustrates the numerical simulation results of solutions from two optimization techniques utilized in this event with stochastic VRESs. indicator. Moreover, Table 7 illustrates the numerical simulation results of solutions from two optimization techniques utilized in this event with stochastic VRESs.   As seen from Table 7, detailed numerical results of control state variables for this scenario utilizing the optimization techniques of MOHHO and MOFPA, for best compromise solutions, are presented. Figure 11a illustrates the superiority of MOHHO over the MOFPA technique that achieves the conflicting objectives of minimum emissions and fuel costs. The best compromise solution is extracted from the PF with TOPSIS performance during all the runs of an algorithm, as shown in Figure 11b,c. Control parameters are all the generator active power (except swing generator TPGU1 linked to bus 1) and generator bus voltages. P PV , P PVSH , and P W are the scheduled power from PV, PV-Hydro, and wind generations, respectively. The tolerable limits of state variables are reported in [51]. However, the POZs that present poor operation in the generation cost curve for the TPGU2 linked to bus 2. The two POZs range between (35,45) MW and (60,65) MW.
Active power and reactive power of all generators are considered as system constraints to be realized by the algorithms. The solar and wind generation units are considered to be able to absorb and deliver reactive power of about 0.4 and 0.5 pu of rated capacity in line, respectively. Moreover, Table 7 also involves the calculation of power loss utilizing Equation (42) and accumulative voltage drop of load buses utilizing Equation (43).
The best cost objective functions attained by the two optimization techniques with the TOPSIS ranking are virtually the same and control variables are also comparable.
In the comparison to minimalization of emission levels, MOFPA achieves lesser emission levels to the value of 0.4521 t/h, while MOHHO achieves emission levels to the value of 0.5221 t/h. Therefore, MOFPA has better diversity PF than MOHHO in the direction of the cost function. In the comparison to total fuel costs, MOHHO achieves lower fuel cost value 798.18 $/h., while MOFPA achieves fuel cost to the value of 804.58 $/h.
In the comparison to convergence, MOHHO marginally better than MOFPA as PF of the former dominates some non-dominating alternatives as shown in Figure 11a.
We can summarize this scenario in the case of best PFs' comparison, the diversity of the PFs is found to be better in MOFPA than in MOHHO. On the other hand, convergence and uniform distribution of solutions are superior in MOHHO. Therefore, it is up to the network operators to choose which objective to exercise this on.

Third Scenario
In this scenario, we retain the biggest thermal unit at bus 1 and replace the thermal units at buses 2 and 8 into NGUs incorporating stochastic VRESs. Figure 12 displays the best PFs of MOHHO and MOFPA with the TOPSIS indicator. Moreover, Table 8 illustrates the numerical simulation results of solutions from two optimization techniques utilized in this event with stochastic VRESs. In the comparison to convergence, MOHHO marginally better than MOFPA as PF of the former dominates some non-dominating alternatives as shown in Figure 11a.
We can summarize this scenario in the case of best PFs' comparison, the diversity of the PFs is found to be better in MOFPA than in MOHHO. On the other hand, convergence and uniform distribution of solutions are superior in MOHHO. Therefore, it is up to the network operators to choose which objective to exercise this on.

Third Scenario
In this scenario, we retain the biggest thermal unit at bus 1 and replace the thermal units at buses 2 and 8 into NGUs incorporating stochastic VRESs. Figure 12 displays the best PFs of MOHHO and MOFPA with the TOPSIS indicator. Moreover, Table 8 illustrates the numerical simulation results of solutions from two optimization techniques utilized in this event with stochastic VRESs.    In the comparison to convergence, MOHHO marginally better than MOFPA as PF of the former dominates some non-dominating alternatives as shown in Figure 12a.
We can summarize this scenario in the case of best PFs' comparison, the diversity of the PFs is found better in MOFPA than in MOHHO. On the other hand, convergence and uniform distribution of solutions are superior in MOHHO. Therefore, it is up to the network operators to choose which objective to exercise this on.
In our proposed study, scenarios II and III are presented based on replacing the thermal units by natural gas units. In these scenarios, there are two points of view in the comparison, where we have a comparison between the algorithms for each scenario, and the other comparison is between the scenarios. It could be observed that the results obtained by the MOHHO technique are environmentally rejected due to high emission levels compared with the other techniques, while its cost objective might encourage the investment. By contrast, MOFPA has the lowest emission level and the highest cost objective. Consequently, the decision on the algorithm superiority might be impossible if the objective has a multidiscipline. This complexity could be resolved if we take a helicopter view of all the results. In other words, it is better to compare firstly between the two scenarios to determine which one feeds the governmental regulations and presents the benefits to the community. First, the results obtained due to old Scenario-I of the algorithms of MOEA/D and SMODE are easily excluded due to the highest levels of emissions and fuel costs compared to other proposed scenarios (scenarios II and III), as illustrated in Table 9. If we take a look at the results obtained due to scenario II, it is found that the results of MOHHO have the highest levels of emissions of 0.5221 t/h. By contrast, the results obtained due to scenario III revealed that MOFPA achieved the lowest levels of emissions of 0.421 t/h and MOHHO achieved the lowest fuel costs of 796.35 $/h. in addition the computational times due to Scenario III is lower than the computational times due to Scenario II.
Consequently, the results of Scenario III satisfy the system needs due to the lower emissions obtained from MOFPA and the lower fuel costs obtained from MOHHO. The question is which of two optimization techniques in Scenario III has presented the best compromise solution. We can say that the power system operator based on the given data can take the decision which optimization technique has the best solution. In general, Scenario III presents the best solution as the biggest thermal unit at bus 1 is remained and replace the thermal units at buses 2 and 8 into NGUs incorporating with stochastic VRESs Each profile in Figures 13 and 14 shows the load bus voltage profile of the two proposed scenarios for worst VD value among the VD values resulting from all non-dominating alternatives utilizing the two presented techniques. In Scenario II, the worst VD described by MOHHO is 0.56 pu whereas that by MOFPA is 0.79 pu. In Scenario III, the worst VD described by MOHHO is 0.5 pu while that by MOFPA is 0.64 pu. The change is because of the highest diversity of MOFPA which results in the technique achieving smaller emission or larger cost objectives. The study of the voltage profiles illustrates that the operational voltages of some buses are near or equal to the maximum allowed values. Consequently, these results also demonstrate the usefulness and effectiveness of a suitable handling constraint method for evolutionary techniques. The grouping of an evolutionary technique and an appropriate handling constraint method like SF can analytically lead the search procedure of an evolutionary algorithm in the direction of global feasible optima. utilizing the two presented techniques. In Scenario II, the worst VD described by MOHHO is 0.56 pu whereas that by MOFPA is 0.79 pu. In Scenario III, the worst VD described by MOHHO is 0.5 pu while that by MOFPA is 0.64 pu. The change is because of the highest diversity of MOFPA which results in the technique achieving smaller emission or larger cost objectives. The study of the voltage profiles illustrates that the operational voltages of some buses are near or equal to the maximum allowed values. Consequently, these results also demonstrate the usefulness and effectiveness of a suitable handling constraint method for evolutionary techniques. The grouping of an evolutionary technique and an appropriate handling constraint method like SF can analytically lead the search procedure of an evolutionary algorithm in the direction of global feasible optima.   utilizing the two presented techniques. In Scenario II, the worst VD described by MOHHO is 0.56 pu whereas that by MOFPA is 0.79 pu. In Scenario III, the worst VD described by MOHHO is 0.5 pu while that by MOFPA is 0.64 pu. The change is because of the highest diversity of MOFPA which results in the technique achieving smaller emission or larger cost objectives. The study of the voltage profiles illustrates that the operational voltages of some buses are near or equal to the maximum allowed values. Consequently, these results also demonstrate the usefulness and effectiveness of a suitable handling constraint method for evolutionary techniques. The grouping of an evolutionary technique and an appropriate handling constraint method like SF can analytically lead the search procedure of an evolutionary algorithm in the direction of global feasible optima.

Conclusions
This paper presents a multi-objective economic-environmental dispatch (MOEED) model for obtaining the best value of Pareto optimal solutions of an integrated IEEE 30-bus of thermal, natural gas, and variable renewable energy sources such as wind, PV, and PV-hydro considering both emission and total cost as a multi-objective function.
Three different scenarios are tested to realize the minimization of both emission and fuel or generation cost. Generally, each scenario has three VRESs of wind, PV, and PV-small hydro units at buses 5, 11, and 13, respectively. It is explained that the first scenario uses three TPGUs and three VRESs. The second scenario replaces the fuel of TPGU at bus 1 into NGU. The third scenario replaces the fuel of TPGUs at buses 2 and 8 into NGUs. The results obtained that achieve minimum emissions and total costs revealed the third scenario is the best compromise solution. MOHHO and MOFPA have been employed to get Pareto-optimal solutions concurrently. All system constraints in terms of equality, inequality, and natural gas limits have been satisfied. A comparative analysis has been implemented between the two optimization techniques to obtain the best values of the multi-objective EED pollutant emissions and fuel costs together. Moreover, the TOPSIS technique was implemented to enable the decision-maker to get the best alternative from the Pareto solutions with diverse preferences. The results obtained show that the MOHHO outperforms MOFPA towards attaining diversity, although the convergence was slightly better than the former.
The presented formulation on MOEED may be studied further utilizing other optimization techniques like the multi-objective zigzag search algorithm (MOZSA), multi-objective particle swarm optimization (MOPSO), etc., together with appropriate constraint handling techniques. Also, the dynamic EED problem with the consideration of variation in load demands over a time-period and generator ramping rate integrated with uncertainties of all the RES and all system limitations remains an issue for future studies. The scale factor of the wind turbine β The shape factor of the wind turbine C d PV The direct cost of the photovoltaic system C d PVSH The direct cost of the photovoltaic-small hydro system C d WT The direct cost of the wind turbine C NG The constant-coefficient of natural gas C r PV The reserve capacity cost of the photovoltaic system C r PVSH The reserve capacity cost of the photovoltaic-small hydro system C r WT The reserve capacity cost of the wind turbine C S PV The storage units cost of the photovoltaic system C S PVSH The storage units cost of the photovoltaic-small hydro system C S WT The storage units cost of the wind turbine C tot The total cost of the fuel or generation C tot PV The total cost of the photovoltaic generation unit C tot PVSH The total cost of the photovoltaic-small hydro generation unit C tot WT The total cost of the wind turbine generation unit C tot (P NGU ) The total cost of the natural gas unit C tot (P TPGU ) The total cost of the thermal power generation unit C tot (P VRES ) The total cost of the variable renewable energy sources δ ij The phase difference between the buses i and j d Pline The internal diameter of the pipe in millimeters η NGU The efficiency of the gas turbine η The efficiency of hydro turbine generator E tot The total emission f v (v) The probability of wind speed f Friction factor γ Scale parameter of the river g NGU i Initial and operation costs coefficient for ith natural gas units G Solar irradiance G std Standard solar irradiance G q(ij) The transconductance of branch q connected to bus i and bus j H W The effective pressure head for the water K d WT The direct cost parameter of the wind turbine K r WT The reserve capacity cost parameter of the wind turbine K S WT The storage unit cost parameter of the wind turbine λ Location parameter of the river L Pline Length of pipeline in meters N G Number of generator buses N L Number of load buses nl Number of branches in the network P 1 Absolute upstream (inlet) pressure P 2 Absolute downstream (outlet) pressure P b Absolute pressure P loss Network power loss P p Pipeline pressure P PV act The actual power of the photovoltaic system P PV r The rated power of the photovoltaic system P PV sch The scheduled power of the photovoltaic system P PVSH act The actual power of the photovoltaic small hydro system P PVSHsch The scheduled power of the photovoltaic-small hydro system P min

TPGU i
The minimum power of the ith thermal power generator unit P W act The actual power of the wind turbine P W r The rated power of the wind turbine P W sch The scheduled power of the wind turbine Q w River flow rate R c Operation irradiance ρ w Water density S Lp The branches' capacity limit S NG The specific gravity of natural gas T NGU The average temperature of the flowing gas in kelvin T s The standard temperature in kelvin U NG The flow velocity of the natural gas in m/sec v The wind speed V G i The voltage of the ith on generator bus v in Cut-in speed of the wind turbine V L p The voltage of the pth on load bus V NGU Volume on the remind loads of the natural gas v out Cut-out speed of the wind turbine v r The rated speed of the wind turbine Z Average compressibility factor of natural gas