Next Article in Journal
B-Learning in Basic Vocational Training Students for the Development of the Module of Applied Sciences I
Previous Article in Journal
Two Forms of the Integral Representations of the Mittag-Leffler Function
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

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

by
Ahmed I. Omar
3,
Ziad M. Ali
1,2,*,
Mostafa Al-Gabalawy
4,
Shady H. E. Abdel Aleem
5 and
Mujahed Al-Dhaifallah
6
1
Electrical Engineering Department, College of Engineering, Prince Sattam bin Abdulaziz University, Wadi Addawaser 11991, Saudi Arabia
2
Electrical Engineering Department, Aswan Faculty of Engineering, Aswan University, Aswan 81542, Egypt
3
Electrical Power and Machines Engineering, The Higher Institute of Engineering at El-Shorouk City, El-Shorouk City 11837, Egypt
4
Pyramids Higher Institute for Engineering and Technology, Giza 12578, Egypt
5
15th of May Higher Institute of Engineering, Mathematical and Physical Sciences, Helwan 14531, Egypt
6
Systems Engineering Department, King Fahd University of Petroleum & Minerals, Dhahran 31261, Saudi Arabia
*
Author to whom correspondence should be addressed.
Mathematics 2020, 8(7), 1100; https://doi.org/10.3390/math8071100
Submission received: 19 May 2020 / Revised: 19 June 2020 / Accepted: 29 June 2020 / Published: 5 July 2020

Abstract

:
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 flow 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 flower pollination algorithm (MOFPA) are employed to find 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 effectiveness 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. Introduction

Recently, a significant increase in the ability of variable renewable energy sources (VRESs) capacity in modern electricity systems has been experienced in response to various environmental, technical, economic, social, and political issues [1,2]. In this regard, the environmental economic dispatch (EED) problem is one of the vital optimization issues in power systems planning and operation, particularly with the rising concerns of global warming and the environmental pollution caused by conventional fossil fuel-based power generation [3,4]. The main principle of the conventional EED problem in power generation is to plan the committed outputs of generating units to satisfy the load demand at the minimum operating cost while satisfying the system constraints as well as taking into account environmental issues to limit pollution caused by thermal power plants simultaneously. From a mathematical point of view [5,6], the EED problem is a non-convex problem that has conflicting objectives and non-linear constraints arising from power flow constraints and grid compliance conditions, valve points effects, and zones of prohibited operation of units, in which efficient algorithms to find the best compromise between environmental and economic requirements are required. Besides, the inclusion of VRESs in electricity systems further complicates the problem [7,8].
In Egypt, much attention has been paid to the integration of VRESs in the Egyptian grid, particularly wind and solar energy, because of the sustainable (economic-techno-environmental-social) added-values. At present, Egypt’s goal is that in 2022 it will be able to save 20% of Egypt’s total energy via renewable energy sources. Because of the large natural gas (NG) discoveries in Egypt, which contributed to meeting Egypt’s NG demand and reserves, the NG is now the enabler for a renewable energy future in Egypt [9]. Therefore, researchers and distribution network operators (DNOs) in the Egyptian electricity market are often asked to look for different alternatives to confirm that networks can handle this significant change in the Egyptian network safely and reliably [10].
In the literature, different mathematical-based or heuristic-based optimization methods for multi-objective optimization had been introduced in many studies to solve the non-convex EED problem in power systems with and without the integration of VRESs [11]. From the perspective of evolutionary and metaheuristic optimization techniques—an improved whale optimization algorithm (IWOA) [12], particle swarm optimization (PSO) [13], an interior search algorithm (ISA) [14], time-varying acceleration coefficient PSO (TVAC-PSO) [15], improved multi-objective moth-flame optimization (IMFO) [16], a combinatorial between the salp swarm algorithm (SSA) and PSO [17], dynamic population-based artificial bee colony (ABC-DP) [18], ABC [19], a multi-objective evolutionary algorithm based on decomposition (MOEA/D) [20]—two to three objective functions to minimize the cost, emissions and power loss are frequently considered in a multi-objective EED problem formulation, as presented in Table 1. As seen, a lot of research work has been done on the classical multi-objective optimal power, which considers only thermal generators. For instance, Medani et al. [12] proposed IWOA for solving the EED problem. In this approach, the authors considered the minimization of active power loss as the objective function but with no VRESs integrated. The results confirmed the robustness and efficacy of IWOA to get the optimal solution and reducing the power loss. Mason et al. [13] employed the PSO variants for solving dynamic EED problem, in which two objectives, cost, and emissions, were presented while considering hourly power demand uncertainty. The results obtained were compared with non-dominated sorting genetic algorithm (NSGA-II) and multi-agent reinforcement learning (MARL) to reveal the effectiveness of the different PSO variants used; however, no ranking method of the non-dominated solutions was employed to obtain the best solution. Karthik et al. [14] used ISA to obtain multi-objective EED (MOEED) problems of minimum fuel cost and emissions. Different optimization techniques were presented to analyze and compare the performance of the studied systems; however, valve points effects and prohibited operating zones (POZs) were not considered in the presented approaches. Nourianfar and Abdi [15] proposed an improved PSO called TVAC-PSO to solve the non-convex EED problem. The introduced technique was applied to many benchmark cases, in addition to a 48-unit combined heat and power (CHP) test system. A ranking procedure called technique for order preference by similarity to ideal solution (TOPSIS) had been utilized to find a single solution, best Pareto optimal front (POF), from the non-dominated solutions set of the EED problem. The results showed the capability of the proposed technique in solving the MOEED problem. Elsakaan et al. [16] proposed the IMFO algorithm to solve the common economic dispatch problem. However, the approach used did not consider some system constraints, such as POZs and system security. Besides, mathematically, arbitrary weightings have been used to turn the problem of multi-objective optimization into a single objective optimization problem. El Sehiemy et al. [17] presented a combinatorial optimization approach of the SSA and PSO swarm optimization for solving the MOEED problem. The presented approach was applied to single and multi-objective optimization techniques with many objectives such as minimization of generation costs, emissions, power loss, and maximization of the voltage stability of the systems studied. Ding et al. [18] and Liang et al. [19] used modified versions of the ABC to solve the MOEED problem with three objectives considered: generation costs, emissions, and power loss. In [18], the ABC-DP was presented with no ranking procedure employed to get the best solution, while, in [19], an improved ABC algorithm based on Pareto optimization was presented, in which the authors pointed out that a non-dominated solution could be selected based on the operator’s choice. Biswas et al. [20] presented the MOEA/D to minimize multiple objectives of emissions, costs, power losses, and voltage deviations in an IEEE 30-bus and IEEE 57-bus systems to solve the MOEED problem with a limited number of handling constraints. However, as seen from Table 1, no ranking procedure has been used in most of the approaches used to reach the best solution, nor have any VRESs been integrated into the systems studied.
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.

2. 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 multi-objective optimization algorithms. The procedure of the proposed techniques is illustrated in Figure 2.

3. 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.

3.1. 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):
Minimize   ( E E D ) = min ( C t o t , E t o t )
where, E t o t and C t o t represent the total emissions and fuel costs to be minimized, respectively.

3.1.1. 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).
C t o t = C t o t ( P T P G U ) + C t o t ( P V R E S s ) + μ × C t o t ( P N G U )
where, C t o t ( P T P G U ) represents the total TPGUs cost, C t o t ( P V R E S ) represents the total VRESs cost, and C t o t ( P N G U ) 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].
C t o t ( P T P G U ) = i = 1 N T P G U a T i + b T i P T P G U i + c T i P T P G U i 2 + | d i × sin ( e i × ( P T P G U i m i n P T P G U i ) ) |
where a T i , b T i and c T i are the cost parameters of the ith thermal generator unit ( P T P G U i ) . The two parameters d i and e i represent the valve point effect. P T P G U i m i n represents the minimum power of P T P G U i during the generator operation. These parameters are specified in Table 4.

Fuel Cost Analysis of Variable Renewable Energy Sources (VRESs)

The VRESs cost, in dollars per megawatt-hour ($/MWh), is the sum of the total costs of WTs ( C t o t ( P W T ) ) , PVs ( C t o t ( P P V ) ) , and hybrid PV and SHP ( C t o t ( P P V S H ) ) which can be expressed as in Equation (4):
C t o t ( P V R E S ) = C t o t ( P W T ) + C t o t ( P P V ) + C t o t ( P P V S H )
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 W T ( P W s c h ) of WTs represents the initial, operation, and maintenance costs, and is expressed in Equation (5).
C d W T ( P W s c h ) = K d W T × P W s c h
where K d W T represents the direct cost parameter and P W s c h 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 W T ) is expressed in Equation (6).
C r W T ( P W s c h P W a c t ) = K r W T ( P W s c h P W a c t ) = K r W T 0 P W a c t ( P W s c h p w ) f w ( p w ) d p w
K r W T represents the cost parameter of the standby units and P W a c t represents the actual delivered power from the WTs. Similarly, when P W a c t is greater than P W s c h , the cost of the storage units will present as described in Equation (7):
C s W T ( P W a c t P W s c h ) = K s W T ( P W a c t P W s c h ) = K s W T P W s c h P W r ( p w P W s c h ) f w ( p w ) d p w
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):
f v ( v ) = ( β α ) ( v α ) ( β 1 ) e ( v α ) β   for   0 < v <
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):
p w = { 0                                   v o u t v v i n P W r ( v v i n v r v i n )                                               v i n v v r P W r                                                                             v r v v o u t
where, v i n , v r , v o u t 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):
f w ( p w ) = β ( v r v i n ) α β × P W r [ v i n + p w P W r ( v r v i n ) ] β 1 × exp [ ( v i n + p w P W r ( v r v i n ) α ) β ]
Finally, the total cost of the wind plant is given in Equation (11):
C t o t W T = C d W T ( P W s c h ) + C r W T ( P W s c h P W a c t ) + C s W T ( P W a c t P W s c h )

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 P V ( P P V s c h ) of PVs represents the initial, operation, and maintenance costs, and is expressed in Equation (12).
C d P V ( P P V s c h ) = K d P V × P P V s c h
where K d P V represents the direct cost coefficient of the PV system and P P V s c h denotes the scheduled power of the PV system. When the actual power provided from the PV system ( P P V a c t ) is less than P P V s c h , then it is required to implement standby units (reserve capacity) to maintain the demand requirements. The cost of the reserve capacity ( C r P V ) is expressed in Equation (13).
C r P V ( P P V s c h P P V a c t ) = K r P V ( P P V s c h P P V a c t ) = K r P V ( P P V s c h p P V ) × f P V ( p P V )
K r P V represents the cost coefficient of the standby units for the PV system. As mentioned before, the cost of the storage units ( C s P V ) will present, if P P V a c t is greater than P P V s c h and is described as in Equation (14):
C s P V ( P P V a c t P P V s c h ) = K s P V ( P P V a c t P P V s c h ) = K s P V ( p P V P P V s c h ) × f P V ( p P V )
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:
f P V ( G ) = 1 G σ 2 π exp { ( l n G μ 2 ) 2 σ 2 } ,           G > 0
The available power from the PV generation unit can be determined as in Equation (16):
p P V ( G ) = { P P V r ( G 2 G s t d ) ,         0 < G < R c P P V r ( G G s t d ) ,                   G R c
where the standard solar irradiance G s t d is equal to 1000 W/m2. During the operation irradiance R c is set as 120 W/m2. P P V r 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 t o t P V ) comprises the summation of the direct cost, standby unit cost, and the storage unit cost, and can be illustrated as in Equation (17):
C t o t P V = C d P V ( P P V s c h ) + C r P V ( P P V s c h P P V a c t ) + C s P V ( P P V a c t P P V s c h )

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:
f Q ( Q w ) = 1 γ e x p ( Q w λ γ ) e x p [ e x p ( Q w λ γ ) ]
The output power from the hydro plant P H ( Q w ) mainly depends on Q w as given in Equation (19):
P H ( Q w ) = η w ρ w g w Q w H w
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 t o t P V S H ) is formulated in the same manner as the WT and PV cost equations, as in Equation (20); thus:
C t o t P V S H = C d P V S H ( P P V S H s c h ) + C r P V S H ( P P V S H s c h P P V S H a c t ) + C s P V S H ( P P V S H a c t P P V S H s c h )
where P P V S H s c h and P P V S H a c t express the scheduled power and the actual power from the hybrid unit, respectively. C d P V S H ( P P V S H ) represents the direct cost of the PVSH. C r P V S H and C S P V S H 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 (   C t o t ( P N G U ) ) is given in Equation (21) [36,37], thus:
  C t o t ( P N G U ) = i = 1 N N G U g N G U i × P N G U i + i = 1 N N G U a N i + b N i P N G U i + c N i P N G U i 2
where a N i , b N i and c N i are the cost parameters of the ith NGU ( P N G U i ) , and g N G U i denotes the coefficient of the initial and operating costs of NGUs. These parameters are specified in Table 5.
To sum up, C t o t of the overall system is given in Equation (22).
C t o t = i = 1 N T P G U a T i + b T i P T P G U i + c T i P T P G U i 2 + | d i × sin ( e i × ( P T P G U i m i n P T P G U i ) ) | + C d W T ( P W s c h ) + C r W T ( P W s c h P W a c t ) + C s W T ( P W a c t P W s c h ) + C d P V ( P P V s c h ) + C r P V ( P P V s c h P P V a c t ) + C s P V ( P P V a c t P P V s c h ) + C d P V S H ( P P V S H s c h ) + C r P V S H ( P P V S H s c h P P V S H a c t ) + C s P V S H ( P P V S H a c t P P V S H s c h ) + i = 1 N N G U g N G U i × P N G U i + i = 1 N N G U a N i + b N i P N G U i + c N i P N G U i 2

3.1.2. Emission Levels

The emissions of thermal and natural gas units are only presented since the VRES has no emissions. The total emission ( E t o t   ) due to these units is described in Equation (23):
E t o t = i = 1 N T P G U E T P G U i + i = 1 N N G U E N G U i
where E T P G U i and E N G U i are the total emissions of TPGU and NGU, respectively.

Emission Analysis of TPGUs

The total emission E t o t ( T P G U ) 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):
E t o t ( T P G U ) = i = 1 N T P G U [ φ T i + ( ψ T i   × P T P G U i ) + ( ω T i × P T P G U i 2 ) + τ T i × e ξ T i P T P G U i ]
where φ T i , ψ T i , ω T i , τ T i and ξ T i are the parameters of emission levels related to the ith T P G U . 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]:
E t o t ( N G U s ) = i = 1 N N G U [ φ N i + ( ψ N i × P N G U i ) + ( ω N i × P N G U i 2 ) ]
where φ N i , ψ N i , and ω N i are the parameters of emission levels related to the ith N G U . These parameters are specified in Table 5.

3.2. Constraints

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

3.2.1. 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.
P G U = P L i + P L o s s i
Q G U = Q L i + Q L o s s i

3.2.2. Rating Limitations

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.
P T P G U i m i n P T P G U i P T P G U i m a x           i N T P G U
P W T m i n P W T P W T m a x
P P V m i n P P V P P V m a x
P P V S H m i n P P V S H P P V S H m a x
P N G U m i n P N G U P N G U m a x
Also, the reactive power operational limits of all generator units are considered as in Equations (33) to (37):
Q T P G U i m i n Q T P G U i Q T P G U i m a x           i N T P G U
Q W T m i n Q W T Q W T m a x
Q P V m i n Q P V Q P V m a x
Q P V H m i n Q P V H Q P V H m a x
Q N G U m i n Q N G U Q N G U m a x

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):
P T P G U i m i n P O Z , j P O Z T P G U i j P T P G U i m a x P O Z , j
where P T P G U i m i n P O Z , j and P T P G U i m a x P O Z , j 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):
V G i m i n V G i V G i m a x             i N G
V L p m i n V L p V L p m a x           i N L
S L p S L p m a x                           i n l
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 n l 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 l o s s ) and voltage quality indicator which called voltage deviation ( V D ) are also considered as network restrictions and can be calculated as in Equations (42) and (43):
P l o s s = q = 1 n l G q ( i j ) [ V i 2 + V j 2 2 V i V j cos ( δ i j ) ]
V D = ( p = 1 N L | V L p 1 | )
where   G q ( i j ) denotes the transconductance of branch q connected to bus i and bus j. δ i j represents the phase difference between δ i and δ j of the buses i and j.

3.2.3. 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 V N G U 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):
V N G U i = 0.278 η N G U i H H V p N G U i
where η N G U i denotes the efficiency of the gas turbine, and H H V represents the high heat value of natural gas. The length, diameter, material type, working pressure, and average flow rate of the natural gas ( q ˜ x y , t ) have the key role before supplying the natural gas into a gas turbine. Mathematically, it is described in Equation (45):
q ˜ x y , t = ± C N G T O p o η p ( p x 2 p y 2 ) D p 5 S N G Z T a v L p f
q ˜ x y , t = q x y , t i n + q x y , t o u t 2
where C N G 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 N G 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 a v 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 N G ) in m/sec is expressed in Equation (47):
U N G = 353 × q N G U × p o D p 2 p p 2 3730   f   L p q ˜ x y , t 2 D p 5
The flow velocity is inversely proportional to the pipeline pressure in which both of the flow velocity U N G 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 xy 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.
L x y , t = G x y ( p x ,   t + p y , t 2 )
L x y , t = L x y , t 1 + q x y , t q y , t
t q x y , t i n = t q x y , t o u t
U N G 20   m / s e c
p p i p e m i n p p i p e p p i p e m a x
q p i p e m i n q p i p e q p i p e m a x
p x ,   t Γ C   p y ,     t

3.3. 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):
    d i = m = 1 M f m i + 1 f m i 1 f m m a x f m m i n , i [ 2 ,   j 1 ]
    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. f m i + 1 and f m i 1 represent the mth objective functions of solution number (i + 1) and (i − 1) in the set Fr, respectively. f m m i n and f m m a x represent the minimum and maximum values of the mth objective function in the set Fr, respectively.

3.3.1. 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].
Rule 1:
Global pollinators fly large distances for pollination, and are close to Levy flights in their movement.
Rule 2:
Local pollination is achieved using abiotic self-pollination.
Rule 3:
Local pollination occurs among the same flowers or flowers of the same species.
Rule 4:
A switch probability ( ρ ) lying in the interval (0,1) decides whether the pollination of a flower is local or global. The algorithm can be formulated as:
(1)
Global pollination is implemented when a uniform production arbitrary value r a n d ρ , illustrated in Equation (56), to sequentially change the location of the ith flower X i utilizing its spacing from the best flower X b e s t .
X i t + 1 = X i t + ɸ L e v ( λ ) × ( X b e s t X i t )
where L e v ( λ ) 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 ( ɸ ) is chosen.
l e v ( λ ) = r a n d 1 | r a n d 2 | 1 / λ × σ 1 ( λ ) σ 2 ( λ )
σ 1 ( λ ) = [ Γ ( 1 + λ ) × sin ( π λ 2 ) Γ ( 1 + λ 2 ) × λ × 2 ( λ 1 2 ) ] 1 / λ
σ 2 ( λ ) = 1
(2)
In the FPA technique, the local pollination will be performed to use a uniform distribution random value ( r a n d i ) that lies from 0 to 1 to regulate the mutation in the ith flower.
X i t + 1 = X i t + r a n d i × ( X j t X k t )
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.
Algorithm 1: Multi-objective flower pollination algorithm (MOFPA) pseudo code
Mathematics 08 01100 i001

3.3.2. 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:
X ( t + 1 ) = X r a n d ( t ) r 1 | X r a n d ( t ) 2 r 2 X ( t ) |                                                         q 0.5     X p r e y ( t ) X a ( t ) r 3 ( L B + r 4 ( U B L B ) )                                   q < 0.5  
where X(t), X(t + 1), and Xprey(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. 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:
X a ( t ) = 1 N i = 1 n X i ( t )
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:
E = 2 E o ( 1 1 T )
E o = 2 r 1 - 1
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.

A. Soft Besiege

The condition is available whenever r 0 and | E o | 0 , which is expressed in Equation (65):
X ( t + 1 ) = X ( t ) E J X p r e y ( t ) 2 X ( t )
X ( t ) = X p r e y ( t ) X ( t )
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 | E | ≺ 0. The victim is exhausted and does not have enough strength to run. This process shall be formulated as follows:
X ( t + 1 ) = X p r e y ( t ) E n X ( t )

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:
X ( t ) = X p r e y ( t ) E J X p r e y ( t ) 2 X ( t )
Z = Y + S + L F ( D )
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:
X ( t + 1 ) = { Y                                       f ( Y ) < F ( y ( t ) ) Z                                       f ( Z ) < F ( y ( t ) )    
  • Hard besiege with progressive quick dive.
  • 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):
    X ( t + 1 ) = { X p r e y ( t ) E J X p r e y ( t ) 2 X m ( t )                                 f ( Y ) < F ( y ( t ) )   Z = Y + S + L F ( D )                                                                                         f ( Z ) < F ( y ( t ) )
Moreover, the pseudo-code of MOHHO is presented as illustrates in Algorithm 2.
Algorithm 2: MOHHO pseudo code.
Mathematics 08 01100 i002

3.3.3. 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 i j + ) has the largest maximum index alteration and the smallest dispersal, and the negative ideal solution ( v i j ) 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 i j 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 j = 1 n w j = 1 .
Step 2:
Determine the normalized value Y i j by applying Equation (72), i.e., by normalizing the decision matrix:
Y i j = G i j i = 2 n G i j 2   i = 1 ,   2 ,   ,   n   &   j = 1 ,   2
Step 3:
Determine V i j using Equation (73) that denotes the weighted normalized decision matrix.
V i j = w j × Y i j     i = 1 ,   2 ,   ,   n   &   j = 1 ,   2
Step 4:
Get v i j + and v i j using Equations (74) and (75):
v i j + = { min ( V 11 , V n 1 ) , min ( V 12 , V n 2 ) }
v i j = { max ( V 11 , V n 1 ) , max ( V 12 , V n 2 ) }
Step 5:
Using the n-dimensional Euclidean distance, determine the separation measures by Equations (76) and (77):
S i j + = i = 2 n ( V i j v i j + ) 2     i = 1 ,   2 ,   ,   n   &   j = 1 ,   2
S i j = i = 2 n ( V i j v i j ) 2     i = 1 ,   2 ,   ,   n   &   j = 1 ,   2
Step 6:
Determine the relative closeness (RC) or performance index to the ideal solution, which has been formulated as in Equation (78):
C i j = S i j S i j + + S i j     i = 1 ,   2 ,   ,   n   &   j = 1 ,   2
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.

4. 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.

4.1. 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.

4.2. 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.
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. PPV, PPVSH, and PW 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.

4.3. 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.
Similarly, in the comparison to minimalization of emission levels, MOFPA achieves lesser emission levels to a value of 0. 0.421 t/h, while MOHHO gets emission levels to value 0.558 t/h. Therefore, MOFPA has better diversity PF than the MOHHO towards the cost function. In the comparison to total fuel costs, MOHHO achieves lower fuel cost valued at 796.35 $/h., while MOFPA gets fuel cost to value 807.89 $/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 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 Figure 13 and Figure 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.

5. 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.

Author Contributions

Conceptualization, A.I.O. and M.A.-G.; data curation, Z.M.A. and M.A.-G.; formal analysis, S.H.E.A.A.; methodology, A.I.O. and M.A.-G.; resources, M.A.-D. and Z.M.A.; software, A.I.O.; supervision, S.H.E.A.A.; validation, M.A.-G.; visualization, S.H.E.A.A.; writing—original draft, A.I.O. and M.A.-G.; writing—review and editing, S.H.E.A.A., Z.M.A. and M.A.-D. All authors together organized and refined the manuscript in the present form. All authors have read and agreed to the published version of the manuscript.

Funding

This research received no external funding.

Conflicts of Interest

The authors declare no conflict of interest.

Abbreviations

ABCArtificial bee colony
ABC-DPDynamic population-based artificial bee colony
CHPCombined heat and power
CSCuckoo search algorithm
DNOsDistribution network operators
EEDEnvironmental economic dispatch
ESEnergy storage
HHVHigh heat value of natural gas
IMFOImproved multi-objective moth-flame optimization
ISAInterior search algorithm
IWOAThe improved whale optimization algorithm
MADMMulti-attribute decision making
MARLMulti-agent reinforcement learning
MOCE/DMulti-objective cross-entropy algorithm based on decomposition
MOEEDMulti-objective environmental economic dispatch
MOEA/DDecomposition-based multi-objective evolutionary algorithm
MOFPAMulti-objective flower pollination algorithm
MOGOAMulti-objective grasshopper optimization algorithm
MOHHOMulti-objective Harris hawks optimization
MOMFOMulti-objective moth-flame optimization
MOOMulti-objective optimization
MOPEOMulti-objective population extremal optimization
MOPSOMulti-objective particle swarm optimization
MOSSAMulti-objective salp search algorithm
MOZSAMulti-objective zigzag search algorithm
NGNatural gas
NGUsNatural gas units
NSGANon-dominated sorting genetic algorithm
PDFsProbability density functions
PFsPareto fronts
POFPareto optimal front
POZsProhibited operating zones
RCRelative closeness to the ideal solution
SFFeasible solution
SMODESummation based multi-objective differential evolution
SSASalp swarm algorithm
TLBOTeaching learning-based optimization
TOPSISThe technique for order preference by similarity to an ideal solution
TPGUsThermal power generation units
TVACTime-varying acceleration coefficient
VDVoltage deviation
VRESsVariable renewable energy sources

Nomenclature

α The scale factor of the wind turbine
β The shape factor of the wind turbine
C d P V The direct cost of the photovoltaic system
C d P V S H The direct cost of the photovoltaic-small hydro system
C d W T The direct cost of the wind turbine
C N G The constant-coefficient of natural gas
C r P V The reserve capacity cost of the photovoltaic system
C r P V S H The reserve capacity cost of the photovoltaic-small hydro system
C r W T The reserve capacity cost of the wind turbine
C S P V The storage units cost of the photovoltaic system
C S P V S H The storage units cost of the photovoltaic-small hydro system
C S W T The storage units cost of the wind turbine
C t o t The total cost of the fuel or generation
C t o t P V The total cost of the photovoltaic generation unit
C t o t P V S H The total cost of the photovoltaic-small hydro generation unit
C t o t W T The total cost of the wind turbine generation unit
C t o t ( P N G U ) The total cost of the natural gas unit
C t o t ( P T P G U ) The total cost of the thermal power generation unit
C t o t ( P V R E S ) The total cost of the variable renewable energy sources
δ i j The phase difference between the buses i and j
d P l i n e The internal diameter of the pipe in millimeters
η N G U The efficiency of the gas turbine
η The efficiency of hydro turbine generator
E t o t The total emission
f v ( v ) The probability of wind speed
f Friction factor
γScale parameter of the river
g N G U i Initial and operation costs coefficient for ith natural gas units
GSolar irradiance
G s t d Standard solar irradiance
G q ( i j ) The transconductance of branch q connected to bus i and bus j
H W The effective pressure head for the water
K d W T The direct cost parameter of the wind turbine
K r W T The reserve capacity cost parameter of the wind turbine
K S W T The storage unit cost parameter of the wind turbine
λLocation parameter of the river
L P l i n e Length of pipeline in meters
N G Number of generator buses
N L Number of load buses
n l Number of branches in the network
P 1 Absolute upstream (inlet) pressure
P 2 Absolute downstream (outlet) pressure
P b Absolute pressure
P l o s s Network power loss
P p Pipeline pressure
P P V a c t The actual power of the photovoltaic system
P P V r The rated power of the photovoltaic system
P P V s c h The scheduled power of the photovoltaic system
P P V S H a c t The actual power of the photovoltaic small hydro system
P P V S H s c h The scheduled power of the photovoltaic-small hydro system
P T P G U i m i n The minimum power of the ith thermal power generator unit
P W a c t The actual power of the wind turbine
P W r The rated power of the wind turbine
P W s c h The scheduled power of the wind turbine
Q w River flow rate
R c Operation irradiance
ρ w Water density
S L p The branches’ capacity limit
S N G The specific gravity of natural gas
T N G U The average temperature of the flowing gas in kelvin
T s The standard temperature in kelvin
U N G 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 i n Cut-in speed of the wind turbine
V L p The voltage of the pth on load bus
V N G U Volume on the remind loads of the natural gas
v o u t Cut-out speed of the wind turbine
v r The rated speed of the wind turbine
Z Average compressibility factor of natural gas

Appendix A

Table A1. Direct, reserve, and standby cost parameters for stochastic VRESs.
Table A1. Direct, reserve, and standby cost parameters for stochastic VRESs.
Wind (Bus 5)Solar (Bus 11)Solar-Hydro (Bus 13)
Direct cost parameters ($/MW) K d W T = 1.7 K d P V = 1.6 K d P V S H = 1.5
Reserve cost parameters ($/MW) K r W T = 3 K r P V = 3 K r P V S H = 3
Penalty cost parameters ($/MW) K s W T = 1.4 K s P V = 1.4 K s P V S H = 1.4

References

  1. Deng, X.; Lv, T. Power system planning with increasing variable renewable energy: A review of optimization models. J. Clean. Prod. 2020, 246. [Google Scholar] [CrossRef]
  2. Yin, L.; Gao, Q.; Zhao, L.; Wang, T. Expandable deep learning for real-time economic generation dispatch and control of three-state energies based future smart grids. Energy 2020, 191, 116561. [Google Scholar] [CrossRef]
  3. Lin, Z.; Chen, H.; Wu, Q.; Li, W.; Li, M.; Ji, T. Mean-tracking model based stochastic economic dispatch for power systems with high penetration of wind power. Energy 2020, 193, 116826. [Google Scholar] [CrossRef] [Green Version]
  4. Wu, Y.; Wang, X.; Xu, Y.; Fu, Y. Multi-objective Differential-Based Brain Storm Optimization for Environmental Economic Dispatch Problem. In Adaptation, Learning, and Optimization; Cheng, S., Shi, Y., Eds.; Springer International Publishing: Cham, Germany, 2019; Volume 23, pp. 79–104. ISBN 978-3-030-15070-9. [Google Scholar]
  5. El-Sayed, W.T.; El-Saadany, E.F.; Zeineldin, H.H.; Al-Sumaiti, A.S. Fast initialization methods for the nonconvex economic dispatch problem. Energy 2020, 201, 117635. [Google Scholar] [CrossRef]
  6. Kaboli, S.H.A.; Alqallaf, A.K. Solving non-convex economic load dispatch problem via artificial cooperative search algorithm. Expert Syst. Appl. 2019, 128, 14–27. [Google Scholar] [CrossRef]
  7. Ismael, S.M.; Abdel Aleem, S.H.E.; Abdelaziz, A.Y.; Zobaa, A.F. State-of-the-art of hosting capacity in modern power systems with distributed generation. Renew. Energy 2019, 130, 1002–1020. [Google Scholar] [CrossRef]
  8. Rizk-Allah, R.M.; Abdel Mageed, H.M.; El-Sehiemy, R.A.; Abdel Aleem, S.H.E.; El Shahat, A. A new sine cosine optimization algorithm for solving combined non-convex economic and emission power dispatch problems. Int. J. Energy Convers. 2017, 5, 180. [Google Scholar] [CrossRef]
  9. Zhang, X.-P.; Ou, M.; Song, Y.; Li, X. Review of Middle East energy interconnection development. J. Mod. Power Syst. Clean Energy 2017, 5, 917–935. [Google Scholar] [CrossRef] [Green Version]
  10. Ismael, S.M.; Abdel Aleem, S.H.E.; Abdelaziz, A.Y.; Zobaa, A.F. Practical considerations for optimal conductor reinforcement and hosting capacity enhancement in radial distribution systems. IEEE Access 2018, 6, 27268–27277. [Google Scholar] [CrossRef]
  11. Qu, B.Y.; Zhu, Y.S.; Jiao, Y.C.; Wu, M.Y.; Suganthan, P.N.; Liang, J.J. A survey on multi-objective evolutionary algorithms for the solution of the environmental/economic dispatch problems. Swarm Evol. Comput. 2018, 38, 1–11. [Google Scholar] [CrossRef]
  12. Khaled, M.; Sayah, S.; Bekrar, A. Whale optimization algorithm based optimal reactive power dispatch: A case study of the Algerian power system. Electr. Power Syst. Res. 2018, 163, 696–705. [Google Scholar] [CrossRef]
  13. Mason, K.; Duggan, J.; Howley, E. Multi-objective dynamic economic emission dispatch using particle swarm optimisation variants. Neurocomputing 2017, 270, 188–197. [Google Scholar] [CrossRef]
  14. Karthik, N.; Parvathy, A.K.; Arul, R. Multi-objective economic emission dispatch using interior search algorithm. Int. Trans. Electr. Energy Syst. 2019, 29, e2683. [Google Scholar] [CrossRef] [Green Version]
  15. Nourianfar, H.; Abdi, H. Solving the multi-objective economic emission dispatch problems using Fast Non-Dominated Sorting TVAC-PSO combined with EMA. Appl. Soft Comput. 2019, 85, 105770. [Google Scholar] [CrossRef]
  16. Elsakaan, A.A.; El-Sehiemy, R.A.; Kaddah, S.S.; Elsaid, M.I. An enhanced moth-flame optimizer for solving non-smooth economic dispatch problems with emissions. Energy 2018, 157, 1063–1078. [Google Scholar] [CrossRef]
  17. El Sehiemy, R.A.; Selim, F.; Bentouati, B.; Abido, M.A. A novel multi-objective hybrid particle swarm and salp optimization algorithm for technical-economical-environmental operation in power systems. Energy 2020, 193, 116817. [Google Scholar] [CrossRef]
  18. Ding, M.; Chen, H.; Lin, N.; Jing, S.; Liu, F.; Liang, X.; Liu, W. Dynamic population artificial bee colony algorithm for multi-objective optimal power flow. Saudi J. Biol. Sci. 2017, 24, 703–710. [Google Scholar] [CrossRef]
  19. Liang, R.-H.; Wu, C.-Y.; Chen, Y.-T.; Tseng, W.-T. Multi-objective dynamic optimal power flow using improved artificial bee colony algorithm based on Pareto optimization. Int. Trans. Electr. Energy Syst. 2016, 26, 692–712. [Google Scholar] [CrossRef]
  20. Biswas, P.P.; Suganthan, P.N.; Mallipeddi, R.; Amaratunga, G.A.J. Multi-objective optimal power flow solutions using a constraint handling technique of evolutionary algorithms. Soft Comput. 2020, 24, 2999–3023. [Google Scholar] [CrossRef]
  21. Wang, G.; Zha, Y.; Wu, T.; Qiu, J.; Peng, J.; Xu, G. Cross entropy optimization based on decomposition for multi-objective economic emission dispatch considering renewable energy generation uncertainties. Energy 2020, 193, 116790. [Google Scholar] [CrossRef]
  22. Chen, M.-R.; Zeng, G.-Q.; Lu, K.-D. Constrained multi-objective population extremal optimization based economic-emission dispatch incorporating renewable energy resources. Renew. Energy 2019, 143, 277–294. [Google Scholar] [CrossRef]
  23. Bora, T.C.; Mariani, V.C.; dos Santos Coelho, L. Multi-objective optimization of the environmental-economic dispatch with reinforcement learning based on non-dominated sorting genetic algorithm. Appl. Therm. Eng. 2019, 146, 688–700. [Google Scholar] [CrossRef]
  24. Biswas, P.P.; Suganthan, P.N.; Qu, B.Y.; Amaratunga, G.A.J. Multiobjective economic-environmental power dispatch with stochastic wind-solar-small hydro power. Energy 2018, 150, 1039–1057. [Google Scholar] [CrossRef]
  25. Yin, Y.; Liu, T.; He, C. Day-ahead stochastic coordinated scheduling for thermal-hydro-wind-photovoltaic systems. Energy 2019, 187, 115944. [Google Scholar] [CrossRef]
  26. Li, X.; Wang, W.; Wang, H.; Wu, J.; Fan, X.; Xu, Q. Dynamic environmental economic dispatch of hybrid renewable energy systems based on tradable green certificates. Energy 2020, 193, 116699. [Google Scholar] [CrossRef]
  27. Elattar, E.E. Environmental economic dispatch with heat optimization in the presence of renewable energy based on modified shuffle frog leaping algorithm. Energy 2019, 171, 256–269. [Google Scholar] [CrossRef]
  28. Joshi, P.M.; Verma, H.K. An improved TLBO based economic dispatch of power generation through distributed energy resources considering environmental constraints. Sustain. Energy Grids Netw. 2019, 18, 100207. [Google Scholar] [CrossRef]
  29. Bai, W.; Eke, I.; Lee, K.Y. An improved artificial bee colony optimization algorithm based on orthogonal learning for optimal power flow problem. Control Eng. Pract. 2017, 61, 163–172. [Google Scholar] [CrossRef]
  30. Biswas, P.P.; Suganthan, P.N.; Amaratunga, G.A.J. Optimal power flow solutions incorporating stochastic wind and solar power. Energy Convers. Manag. 2017, 148, 1194–1207. [Google Scholar] [CrossRef]
  31. Hulio, Z.H.; Jiang, W.; Rehman, S. Techno—Economic assessment of wind power potential of Hawke’s Bay using Weibull parameter: A review. Energy Strateg. Rev. 2019, 26, 100375. [Google Scholar] [CrossRef]
  32. Panda, A.; Tripathy, M. Security constrained optimal power flow solution of wind-thermal generation system using modified bacteria foraging algorithm. Energy 2015, 93, 816–827. [Google Scholar] [CrossRef]
  33. Zhai, R.; Liu, H.; Chen, Y.; Wu, H.; Yang, Y. The daily and annual technical-economic analysis of the thermal storage PV-CSP system in two dispatch strategies. Energy Convers. Manag. 2017, 154, 56–67. [Google Scholar] [CrossRef]
  34. Tan, Q.; Mei, S.; Dai, M.; Zhou, L.; Wei, Y.; Ju, L. A multi-objective optimization dispatching and adaptability analysis model for wind-PV-thermal-coordinated operations considering comprehensive forecasting error distribution. J. Clean. Prod. 2020, 256, 120407. [Google Scholar] [CrossRef]
  35. Gómez, Y.M.; Bolfarine, H.; Gómez, H.W. Gumbel distribution with heavy tails and applications to environmental data. Math. Comput. Simul. 2019, 157, 115–129. [Google Scholar] [CrossRef]
  36. Tso, W.W.; Demirhan, C.D.; Floudas, C.A.; Pistikopoulos, E.N. Multi-scale energy systems engineering for optimal natural gas utilization. Catal. Today 2019. [Google Scholar] [CrossRef]
  37. Chen, S.; Wei, Z.; Sun, G.; Wang, D.; Zhang, Y.; Ma, Z. Stochastic look-ahead dispatch for coupled electricity and natural-gas networks. Electr. Power Syst. Res. 2018, 164, 159–166. [Google Scholar] [CrossRef]
  38. Algabalawy, M.A.; Abdelaziz, A.Y.; Mekhamer, S.F.; Abdel Aleem, S.H.E. Considerations on optimal design of hybrid power generation systems using whale and sine cosine optimization algorithms. J. Electr. Syst. Inf. Technol. 2018, 5, 312–325. [Google Scholar] [CrossRef]
  39. Li, G.; Zhang, R.; Jiang, T.; Chen, H.; Bai, L.; Li, X. Security-constrained bi-level economic dispatch model for integrated natural gas and electricity systems considering wind power and power-to-gas process. Appl. Energy 2017, 194, 696–704. [Google Scholar] [CrossRef]
  40. Avalos, R.; Fitzgerald, T.; Rucker, R.R. Measuring the effects of natural gas pipeline constraints on regional pricing and market integration. Energy Econ. 2016, 60, 217–231. [Google Scholar] [CrossRef] [Green Version]
  41. Abul’Wafa, A.R. Optimization of economic/emission load dispatch for hybrid generating systems using controlled Elitist NSGA-II. Electr. Power Syst. Res. 2013, 105, 142–151. [Google Scholar] [CrossRef]
  42. Dhanalakshmi, S.; Kannan, S.; Mahadevan, K.; Baskar, S. Application of modified NSGA-II algorithm to Combined Economic and Emission Dispatch problem. Int. J. Electr. Power Energy Syst. 2011, 33, 992–1002. [Google Scholar] [CrossRef]
  43. Zhao, F.; Yuan, J.; Wang, N. Dynamic economic dispatch model of microgrid containing energy storage components based on a variant of NSGA-II Algorithm. Energies 2019, 12, 871. [Google Scholar] [CrossRef] [Green Version]
  44. Basu, M. Combined heat and power economic emission dispatch using nondominated sorting genetic algorithm-II. Int. J. Electr. Power Energy Syst. 2013, 53, 135–141. [Google Scholar] [CrossRef]
  45. Yang, X.-S.; Karamanoglu, M.; He, X. Flower pollination algorithm: A novel approach for multiobjective optimization. Eng. Optim. 2014, 46, 1222–1237. [Google Scholar] [CrossRef] [Green Version]
  46. Abdelaziz, A.Y.; Ali, E.S.; Abd Elazim, S.M. Flower pollination algorithm to solve combined economic and emission dispatch problems. Eng. Sci. Technol. Int. J. 2016, 19, 980–990. [Google Scholar] [CrossRef] [Green Version]
  47. Heidari, A.A.; Mirjalili, S.; Faris, H.; Aljarah, I.; Mafarja, M.; Chen, H. Harris hawks optimization: Algorithm and applications. Future Gener. Comput. Syst. 2019, 97, 849–872. [Google Scholar] [CrossRef]
  48. Aleem, S.H.E.A.; Zobaa, A.F.; Balci, M.E.; Ismael, S.M. Harmonic overloading minimization of frequency-dependent components in harmonics polluted distribution systems using harris hawks optimization algorithm. IEEE Access 2019, 7, 100824–100837. [Google Scholar] [CrossRef]
  49. Lin, Y.-K.; Chang, P.-C.; Yeng, L.C.-L.; Huang, S.-F. Bi-objective optimization for a multistate job-shop production network using NSGA-II and TOPSIS. J. Manuf. Syst. 2019, 52, 43–54. [Google Scholar] [CrossRef]
  50. Deb, M.; Debbarma, B.; Majumder, A.; Banerjee, R. Performance-emission optimization of a diesel-hydrogen dual fuel operation: A NSGA II coupled TOPSIS MADM approach. Energy 2016, 117, 281–290. [Google Scholar] [CrossRef]
  51. Zhang, H.; Lu, Z.; Hu, W.; Wang, Y.; Dong, L.; Zhang, J. Coordinated optimal operation of hydro-wind-solar integrated systems. Appl. Energy 2019, 242, 883–896. [Google Scholar] [CrossRef]
Figure 1. Adopted IEEE 30-bus system studied.
Figure 1. Adopted IEEE 30-bus system studied.
Mathematics 08 01100 g001
Figure 2. The procedure of the proposed techniques.
Figure 2. The procedure of the proposed techniques.
Mathematics 08 01100 g002
Figure 3. Fuel cost function with and without the impact of valve point.
Figure 3. Fuel cost function with and without the impact of valve point.
Mathematics 08 01100 g003
Figure 4. Probability density functions (PDFs): (a) Weibull distribution of wind speed, (b) Lognormal distribution of solar irradiance, and (c) Gumbel distribution of river flow rate.
Figure 4. Probability density functions (PDFs): (a) Weibull distribution of wind speed, (b) Lognormal distribution of solar irradiance, and (c) Gumbel distribution of river flow rate.
Mathematics 08 01100 g004aMathematics 08 01100 g004b
Figure 5. The flow of natural gas in pipelines.
Figure 5. The flow of natural gas in pipelines.
Mathematics 08 01100 g005
Figure 6. Pareto front (PF) plot [42].
Figure 6. Pareto front (PF) plot [42].
Mathematics 08 01100 g006
Figure 7. Calculation of crowded distance.
Figure 7. Calculation of crowded distance.
Mathematics 08 01100 g007
Figure 8. Flowchart of multi-objective flower pollination algorithm (MOFPA) [46].
Figure 8. Flowchart of multi-objective flower pollination algorithm (MOFPA) [46].
Mathematics 08 01100 g008
Figure 9. Different steps of the multi-objective Harris hawks optimization (MOHHO) technique [48].
Figure 9. Different steps of the multi-objective Harris hawks optimization (MOHHO) technique [48].
Mathematics 08 01100 g009
Figure 10. The method of implementing the proposed MOHHO technique.
Figure 10. The method of implementing the proposed MOHHO technique.
Mathematics 08 01100 g010
Figure 11. Obtained PFs by using MOHHO and MOFPA techniques. (a) Obtained PF using the two optimization techniques. (b) Obtained PF using the MOHHO technique with TOPSIS. (c) Obtained PF using the MOFPA technique with TOPSIS.
Figure 11. Obtained PFs by using MOHHO and MOFPA techniques. (a) Obtained PF using the two optimization techniques. (b) Obtained PF using the MOHHO technique with TOPSIS. (c) Obtained PF using the MOFPA technique with TOPSIS.
Mathematics 08 01100 g011
Figure 12. Obtained PFs by using MOHHO and MOFPA techniques with TOPSIS. (a) Obtained PF using the two optimization techniques. (b) Obtained PF using the MOHHO technique with TOPSIS. (c) Obtained PF using the MOFPA technique with TOPSIS.
Figure 12. Obtained PFs by using MOHHO and MOFPA techniques with TOPSIS. (a) Obtained PF using the two optimization techniques. (b) Obtained PF using the MOHHO technique with TOPSIS. (c) Obtained PF using the MOFPA technique with TOPSIS.
Mathematics 08 01100 g012aMathematics 08 01100 g012b
Figure 13. Load bus voltage profiles for worst VD values utilizing the two optimizations in Scenario II.
Figure 13. Load bus voltage profiles for worst VD values utilizing the two optimizations in Scenario II.
Mathematics 08 01100 g013
Figure 14. Load bus voltage profiles for worst VD values utilizing the two optimizations in Scenario III.
Figure 14. Load bus voltage profiles for worst VD values utilizing the two optimizations in Scenario III.
Mathematics 08 01100 g014
Table 1. Summary of multi-objective economic-environmental dispatch (MOEED) of thermal plants with no variable renewable energy sources (VRESs) considered.
Table 1. Summary of multi-objective economic-environmental dispatch (MOEED) of thermal plants with no variable renewable energy sources (VRESs) considered.
Ref.YearSystems UsedOptimization TechniquesComparative AnalysisObjective Functions *Ranking MethodComments
1234
[12]2018IEEE 14-bus, IEEE 30-bus and Algerian 114-busIWOAPSO
PSO-TVAC
Different cases are examined. No ranking method is used to get the best solution
[13]2017Hypothetical systemPSONSGA-II
MARL
The aims are to test different variants of PSO. Many constraints are not considered.
[14]20193, 10, 20 and 40 generating units and IEEE 30-busISAGA
FPA
CS
ABC
POZ and valve point effects constraints are not considered.
[15]201948-unit CHPTVAC-PSOFCMSpinning reserve requirements, ramp rate limits, valve points effects, and multiple fuel units are considered as additional constraints. TOPSIS is utilized to get the best optimal solutions.
[16]20186, 40 and 80 generating unitsIMFOFPA
GA
PSO
Multi-objective optimization is transformed into a single objective, which made the problem simple. No ranking method is utilized. Some constraints are not considered.
[17]2017IEEE 30-busABC-DPNSGA-II
MOABC
Only equality and inequality constraints are considered.
[18]2016IEEE 30-bus
IEEE 118-bus
ABCIABC
GA
DE
Operators can select one of the non-dominated solutions according to the situation. Some constraints are not considered.
[19]2020IEEE 30-bus
IEEE 57-bus
IEEE 118-bus
SSA-PSOABC
GA
DE
WOA
Transformer tapping limits, voltage magnitude limits of load buses, and power flow limits of transmission lines are considered as inequality constraints. No ranking method is employed to get the best solution.
[20]2019IEEE 30-bus
IEEE 57-bus
MOEA/D MTLBO
MGBICA
MOICA
POZ and valve point effects constraints are not considered. No ranking method is employed to obtain the best solution.
* 1 denotes fuel cost, 2 denotes emission, 3 denotes power loss, and 4 denotes voltage stability.
Table 2. Summary of MOEED of thermal plants integrated with VRESs.
Table 2. Summary of MOEED of thermal plants integrated with VRESs.
Ref.YearTest SystemOptimization TechniquesComparative AnalysisObjective Functions *Ranking IndexComments
1234
[21]2020IEEE 30
IEEE 118
MOCE/DPSO
NSGA-II
Valve point effects are not considered. No ranking method is employed to get the best solution.
[22]2019IEEE 30MOPEODE
NSGA-II
Other sources such as small-hydro power and electric vehicles are not considered in the problem formulation.
[23]2019IEEE 30NSGA-RLNSGA-IIGenerational distance and spread of evaluation performance matrices are used to get the best optimal solution.
[24]2018IEEE 30MOEA/D
SMODE
Stochastic natures of VRESs are considered. The obtained results were not compared to other algorithms to evaluate its performance.
[25]20196 bus power systemCopula functionSingle objective optimization is applied. No ranking method is used. Some constraints are not considered.
[26]2020IEEE 39MOMFONSGA-II
MOPSO
Multi-objective functions are not considered. Only equality and inequality constraints are considered.
[27]2019Different test systemsMSFLASLFA
GA
TLBO
Single objective optimization is applied. Various constraints are included. Valve point effects are not considered.
[28]2019Different systemsTLBO-PSOGA
CTLBO
Multi-objective functions are investigated. Various constraints are not considered such as valve point effects.
* 1 denotes fuel cost, 2 denotes emission, 3 denotes power loss, and 4 denotes voltage stability.
Table 3. The model specifications [24].
Table 3. The model specifications [24].
ItemQuantitySpecifications
Generators63 TPGUs and 3 VRESs
TPGUs3Bus 1 (slack), bus 2, and bus 8
Wind turbine (WT)25Bus 5, 75 MW
Photovoltaic array (PV)1Bus 11, 50 MW
Hybrid PV and small-hydro (PVSH)1Bus 13, 45 + 5 MW
Active load demand-283.4 MW
Reactive load demand-126.2 MVAr
Number of PQ buses2424 load buses
Minimum load voltage allowed-0.95 pu
Maximum load voltage allowed-1.10 pu
Table 4. Emission and cost parameters of the TPGUs [30].
Table 4. Emission and cost parameters of the TPGUs [30].
Emission Parameters
GeneratorBus φ T ψ T ω T τ T ξ T
(t/h)(t/pu. MWh)(t/pu. MW2h)(t/h)(pu. MW−1)
TPGU110.04091−0.055540.06490.00026.667
TPGU220.02543−0.060470.056380.00053.333
TPGU380.05326−0.03550.03380.0022
Cost Parameters
GeneratorBus a T b T c T d T e T
($/h)($/MWh)($/MW2h)($/h)(MW−1)
TPGU113020.00375180.037
TPGU22251.750.0175160.038
TPGU38203.250.00834120.045
Table 5. Emission and cost parameters of the NGUs [38].
Table 5. Emission and cost parameters of the NGUs [38].
Emission Parameters
GeneratorBus φ N ψ N ω N
(t/h)(t/pu. MWh)(t/pu. MW2h)
NGU110.02091−0.075540.04490
NGU220.02543−0.050470.03638
NGU380.03326−0.055500.01380
Cost Parameters
GeneratorBus a N b N c N
($/h)($/MWh)($/MW2h)
NGU11141.060.00175
NGU22151.050.0105
NGU38171.250.02434
Table 6. Detailed numerical results of the optimization techniques for Scenario I.
Table 6. Detailed numerical results of the optimization techniques for Scenario I.
State VariablesMin.Max.MOEA/DSMODE
PNGU1 (MW)50140117.118111.91
PTPGU2 (MW)20806565
PTPGU3 (MW)103518.40323.555
Pw (MW)07555.44754.058
Ppv (MW)05017.64918.436
Ppvh (MW)05015.32615.755
Q1 (MVAr)−501402.1282.788
Q2 (MVAr)−206021.4134.504
Q5 (MVAr)−157037.72736.169
Q8 (MVAr)−306027.10220.376
Q11 (MVAr)−203024.91123.41
Q13 (MVAr)−202520.32815.62
V1 (pu)0.961.101.0761.0761
V2 (pu)0.961.101.06481.0662
V5 (pu)0.961.101.04441.0362
V8 (pu)0.961.101.04021.0362
V11 (pu)0.961.101.08781.0778
V13 (pu)0.961.101.06021.0432
Ploss (MW) 5.54295.3148
VD (pu) 0.45300.4215
Total cost ($/h) 919.040927.049
Emission (t/h) 0.62210.4721
Table 7. Detailed numerical results of the optimization techniques for Scenario II.
Table 7. Detailed numerical results of the optimization techniques for Scenario II.
State VariablesMin.Max.MOHHOMOFPA
PTPGU1 (MW)5014072.8464.86
PTPGU2 (MW)208080.0074.20
PTPGU3 (MW)103535.0034.97
Pw (MW)07533.7431.35
Ppv (MW)05039.7339.60
Ppvh (MW)05028.4337.83
Q1 (MVAr)−5014027.31−29.19
Q2 (MVAr)−206016.3367.74
Q5 (MVAr)−157030.8864.07
Q8 (MVAr)−306051.1828.01
Q11 (MVAr)−20302.4511.82
Q13 (MVAr)−20253.45−8.32
V1 (pu)0.961.101.091.06
V2 (pu)0.961.101.081.10
V5 (pu)0.961.101.051.08
V8 (pu)0.961.101.051.03
V11 (pu)0.961.101.011.03
V13 (pu)0.961.101.010.97
Ploss (MW) 4.504.89
VD (pu) 0.560.79
Wgencost 102.39124.88
Sgencost 147.24150.24
Shgencost 97.38143.72
mass flow (m3) 63.5571.36401
Emission of P2 and P3 (t/h) 0.06910.066868
Emission of P1 (t/h) 0.01540.005905
P1_cost 9.5310.7046
Total cost ($/h) 798.18804.58
Emission (t/h) 0.52210.4521
Fuelvlvcost 340.82327.40
Table 8. Detailed numerical results of the optimization techniques for scenario III.
Table 8. Detailed numerical results of the optimization techniques for scenario III.
State VariablesMinMaxMOHHOMOFPA
PTPGU1 (MW)50140121.4852.81
PNGU2 (MW)208080.0080.00
PTPGU3 (MW)103535.0043.19
Pw (MW)07569.5652.29
Ppv (MW)05023.9332.95
Ppvh (MW)05018.5528.58
Q1 (MVAr)−5014030.12−5.72
Q2 (MVAr)−2060−5.866.15
Q5 (MVAr)−154017.8536.58
Q8 (MVAr)−303548.5549.17
Q11 (MVAr)−20255.5217.46
Q13 (MVAr)−202545.6826.71
V1 (pu)0.961.101.051.00
V2 (pu)0.961.101.051.00
V5 (pu)0.961.101.000.99
V8 (pu)0.961.101.011.01
V11 (pu)0.961.101.021.06
V13 (pu)0.961.101.081.05
Ploss (MW) 7.214.04
VD (pu) 0.500.64
Wgencost 99.53193.20
Sgencost 32.9051.60
Shgencost 27.8939.78
mass flow_02 172.83172.83
mass flow_08 75.6293.31
Total cost ($/h) 796.35807.89
Emission (t/h) 0.5580.421
Fuelvlvcost 336.87147.93
Table 9. Comparative analysis between optimization techniques for three scenarios.
Table 9. Comparative analysis between optimization techniques for three scenarios.
ScenariosScenario IScenario IIScenario III
AlgorithmsMOEA/DSMODEMOHHOMOFPAMOHHOMOFPA
No. of iterations200200200200200200
No. of population200200200200200200
Control parameters δ = 0.9
Mutation factor = 0.5
Crossover rate = 0.9
Mutation factor = 0.5
Crossover rate = 0.9
r   1 = [0,1]
r 5   = [0,1]
ρ = 0.5
λ = 1.5
r a n d i = [0,1]
r   1 = [0,1]
r 5   = [0,1]
ρ   = 0.5
λ   = 1.5
r a n d i   = [0,1]
Computation time (min)45.1445.644.2314.084.1613.34
Total cost ($/h)919.040927.049798.18804.58796.35807.89
Emission (t/h)0.62210.47210.52210.45210.5580.421
* The value in boldface represents the best value achieved by the proposed algorithms.

Share and Cite

MDPI and ACS Style

Omar, A.I.; Ali, Z.M.; Al-Gabalawy, M.; Abdel Aleem, S.H.E.; Al-Dhaifallah, M. Multi-Objective Environmental Economic Dispatch of an Electricity System Considering Integrated Natural Gas Units and Variable Renewable Energy Sources. Mathematics 2020, 8, 1100. https://doi.org/10.3390/math8071100

AMA Style

Omar AI, Ali ZM, Al-Gabalawy M, Abdel Aleem SHE, Al-Dhaifallah M. Multi-Objective Environmental Economic Dispatch of an Electricity System Considering Integrated Natural Gas Units and Variable Renewable Energy Sources. Mathematics. 2020; 8(7):1100. https://doi.org/10.3390/math8071100

Chicago/Turabian Style

Omar, Ahmed I., Ziad M. Ali, Mostafa Al-Gabalawy, Shady H. E. Abdel Aleem, and Mujahed Al-Dhaifallah. 2020. "Multi-Objective Environmental Economic Dispatch of an Electricity System Considering Integrated Natural Gas Units and Variable Renewable Energy Sources" Mathematics 8, no. 7: 1100. https://doi.org/10.3390/math8071100

APA Style

Omar, A. I., Ali, Z. M., Al-Gabalawy, M., Abdel Aleem, S. H. E., & Al-Dhaifallah, M. (2020). Multi-Objective Environmental Economic Dispatch of an Electricity System Considering Integrated Natural Gas Units and Variable Renewable Energy Sources. Mathematics, 8(7), 1100. https://doi.org/10.3390/math8071100

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop