Renewable Energy Micro-Grid Interfacing: Economic and Environmental Issues

: This paper presents a study on the technical, economic, and environmental aspects of renewable energy resources-based distributed generation units (DGs). These units are connected to the medium-voltage network to create a new structure called a microgrid (MG). Renewable energies, especially wind and solar, are the most important generation units among DGs. The stochastic behavior of renewable resources increases the need to ﬁnd the optimum operation of the MG. The optimal operation of a typical MG aims to simultaneously minimize the operational costs and the accompanied emission pollutants over a daily scheduling horizon. Several renewable DGs are investigated in the MG, consisting of biomass generators (BGs), wind turbines (WTs), and photovoltaics (PV). For the proposed operating strategy of the MG, a recent equilibrium optimization (EO) technique is developed and is inspired by the mass balance models for a control volume that are used to estimate their dynamic and equilibrium states. The uncertainties of wind speed and solar irradiation are considered via the Weibull and Beta-probability density functions (PDF) with different states of mean and standard deviation for each hour, respectively. Based on the developed EO, the hourly output powers of the PV, WT, and BGs are optimized, as are the associated power factors of the BGs. The proposed MG operating strategy based on the developed EO is tested on the IEEE 33-bus system and the practical large-scale 141-bus system of AES-Venezuela in the metropolitan area of Caracas. The simulation results demonstrate the signiﬁcant beneﬁts of the optimal operation of a typical MG using the developed EO by minimizing the operational costs and emissions while preserving the penetration level of the DGs by 60%. Additionally, the voltage proﬁle of the MG operation for each hour is highly enhanced where the minimum voltage at each hour is corrected within the permissible limit of [0.95–1.05] Pu. Moreover, the active power losses per hour are greatly reduced. WT output power considering the uncertainties in wind speed.


Introduction
Due to the continuous increase in power demand and rapid depletion of fossil fuels, researchers all over the world have no other option but to look for alternative energy sources by utilizing small-scale distributed power generation (DG) and energy storage systems (ESS) [1]. However, due to the inherent intermittency and volatility of renewable energy sources (RESs), its large-scale integration into the power system will increase the regulation burden and affect the security operation of the main grid. The microgrid is defined as smallscale controllable electrical distribution systems, which have the important advantage of operating either islanded or interconnected to the main grid. Microgrids (MGs) are usually composed of Distributed Energy Resources (DERs), ESSs, and controllable loads [2]. The DERs are based on conventional resources, such as diesel generators and RESs such as photovoltaics (PV) and wind turbines (WT) [3][4][5][6]. There are two modes of operation for MGs. In a grid-connected operation, the MGs draw/supply power from/to the grid based on load and generation conditions with regard to the market prices. On the other hand, it will be disconnected from the grid to provide electricity to associated critical loads in the event of faults [7].
In the first mode of operation, the MGs should be operated economically and reliably, where a supervisory control and data acquisition (SCADA) system is activated to monitor, control, and dispatch all DERs to guarantee the economic and secure operation of the MGs. For this targeted operation of the MGs, various conventional and artificial intelligent programming techniques were applied for dispatching the DERs [8]. In [9][10][11][12], the operation of distribution systems was optimally controlled via DERs commitment, Capacitor Banks (CBs) switching, SVC, and reconfiguration using the jellyfish search algorithm and manta ray foraging optimization algorithm, respectively. In both studies, the wasted energy of power losses was minimized considering the daily load variations, but the uncertainties of the DERs were not taken into account. In [13], the optimal operation and energy management method for a hybrid MG including photovoltaics, wind turbines, a pump as a turbine system, and a diesel generator was introduced, with a study on day-ahead scheduling. The optimal energy management minimizes the fuel cost of diesel generators, the daily operating cost, as well as the balance between the generation and load for both warm and cold days using an imperialist competitive algorithm. In [14], a hybrid ant lion optimization with a bat algorithm was utilized for the power management of the MG considering a droop controller strategy. The main target of this droop controller was stabilized by the MG by minimizing the errors of real and reactive power under a power shortage and power maximum. In [15], an optimal energy scheduling mechanism was presented in multi-MGs in order to minimize the total operational costs of their committed DERs. This study provided different DERs types and their associated uncertainties in a multi-MG system, but the linked lines and their losses were completely ignored. In [16], a Tabu search algorithm was applied for the design of the MG system components with minimum investment, operation, and emission costs. This study utilized the Monte Carlo simulation to deal with the uncertainties due to load forecasting and the random outages of the units. However, the uncertainties due to the intermittent sources of WTs and PVs were not taken into consideration where their outputs were directly evaluated from the daily wind velocity and solar irradiance, respectively. In [17], particle swarm optimization (PSO) was dedicated to minimize the MG operational costs considering the variations in loadings, DGs, and requirements of stable grid operational constraints.
Despite that, only the fuel costs of the committed DGs were handled where the quadratic cost models were utilized for the fuel cells and the micro turbines. In [18], a manta ray foraging optimizer was developed to optimally solve the economic dispatch problem with wind power inclusion considering the valve point effects of the generators, while the wind power effects were ignored in [19,20]. However, the environmental impacts of this operation were not considered. In [21], an artificial ecosystem-based optimizer was applied considering the demand side management for minimizing the techno-economic evaluations of hybrid energy systems. In [22], an optimal operational strategy was presented for MGs including hydrogen storage to integrate RESs and decrease the emissions. In [23], in off-grid MG simulation and tests, an adaptable energy management strategy relying on a mixed energy systems has been provided to preserve the stable operating condition of the off-grid MG and lengthen the lifetime of batteries.
In [24], a harmony search algorithm was combined with differential evolution for the optimal operation of MGs. Added to that, an optimal operation of MGs was considered with WTs, PVs, battery energy storage (BES) systems, electric vehicles (EVs), and demand response for minimizing the total operating costs [25]. In this study, Sparse Nonlinear OPTimizer (SNOPT) solver was utilized using Generalized Algebraic Modeling Systems (GAMS) software. In both studies [24,25], the outputs of WTs and PVs were directly evaluated from the hourly wind velocity and solar irradiance, respectively. Therefore, the uncertainties due to the intermittent sources require an effective handling procedure. In [26,27] a moth flame optimization algorithm was employed for the optimal operation of a hybrid energy system including WTs, PVs, gas turbines, and energy storage.
In [28], an optimal MG operation strategy was presented to minimize the fuel costs and the produced emissions. In this study, the DER uncertainties were considered via the probability distributions and confidence. However, the application was limited to a very small number of sources where linear programming and quadratic programming (QP) were utilized as solvers. In [29], RES was allocated and optimized in a distribution system using Mixed Integer-Linear Programming (MILP) and was solved by FICO ® XPRESS optimization software. In [30], the Energy Storage (ES) and RES were integrated using MILPmethod.
Recently, a new effective optimization algorithm of Equilibrium optimizer (EO) was presented [31]. EO provides strong exploratory and exploitative search mechanisms to adjust solutions at random, assisting in local minima avoidance in the optimization process, which is a common drawback of many optimization algorithms [31]. EO was efficiently used for solving the optimal power flow (OPF) problem in the AC power systems [32] and hybrid AC/DC power grids [33]. In [34,35], EO was utilized for handling the economic and the combined economic environmental dispatch problems, respectively, considering the power constraints, effects of the valve point, transmission losses, and ramp rate limits. In [36], an adaptive EO was developed for an optimal allocation procedure of biomass DGs to enhance the performance of the distribution systems and to reduce the related environmental emissions. In [37], EO was used to deal with the energy management optimization (EMO) in the MG considering the variations of WTs, PVs, and load demand for cost minimization and voltage magnitude improvements. In [38], an improved EO integrated for with optimal allocation of multiple PV units with batteries has been described. In [39], EO was utilized for the EMO considering the energy storage devices and the emissions from the associated DERs. In [40], EO was employed for estimating the undefined parameters for the lithium-ion batteries. In [41], EO was employed for identifying the prediction of oil breakdown voltage considering the barrier impact.
In this paper, an optimal operating strategy based on EO is developed for the technoeconomic and environmental optimization scheme for MGs with multiple RESs. For this target, two objective functions are represented for minimizing the generation costs and minimizing the emissions of environmental pollution caused by them. The proposed algorithm is tested on two systems in order to verify its effectiveness and efficiency. The two systems are IEEE 33-bus and a practical large-scale 141-bus system of AES-Venezuela in the metropolitan area of Caracas. Added to that, EO is compared with other recent algorithms of DE and RAO algorithms. The main contributions of this paper can be summarized as follows:

•
Proposing an operational optimization problem of MG incorporating PV, WT, and BG with consideration of the uncertainties of PV and wind as RESs.

•
Simultaneously minimizing the operational costs and the accompanied pollutants.

•
Hourly load variations over a 24-h scheduling horizon are handled. • EO technique is employed with higher performance compared to DE and RAO.

•
The validity of the proposed methodology is validated on a large-scale 141-bus real distribution system.
Four additional rest sections of the current paper are organized thusly: Section 2 shows the mathematical formulation of the MG operation. The EO algorithm is developed in Section 3. Section 4 reveals the simulation results on the two small and large scale tested distribution systems. Section 5 concludes the paper findings.

Problem Formulation
DGs based on renewable energies are the main technologies that greatly affect the operation of the MG units. Because of their stochastic nature, optimal MG operation has to effectively handle their types and associated uncertainties in hourly basis. In this section, an optimization scheme is presented for the techno-economic and environmental operation of MGs that are established with several DGs based on renewable energies, including PV/WT/BG considering hourly variable load demand.

Types of DGs
Four types of the DG resources are classified based on their ability to deliver active and reactive powers [42][43][44]. In the first type, the DG units supply active power only. In Type-2, the DG units provide reactive power only. The DG units in Type 3 inject real and reactive power. The 4th type of DG units injects the active power and consumes reactive power. In this work, the WT model is considered Type-4. The related reactive power to be consumed by these units is obtained by [45].
where P DG,i and Q DG,i are the active and reactive powers of the DG unit i, respectively.

Photovoltaic DGs
The small-scale PV DGs are the most common renewable sources in the MGs. PV modules are usually modeled using single, double, and triple diode-equivalent circuits [46,47]. Unfortunately, the produced power is intermittent and varied in each hour with high uncertainty levels due to their dependency on the solar irradiances. Therefore, they should be treated with an effective way for each hour where the uncertainty of solar irradiance is modeled by the Beta Probability Density Function (PDF). Consequently, in each hour, different states of the solar irradiance are considered to generate its Beta-PDF. For each state of solar irradiance s, the output power from the PV module, P pv0 (s), can be expressed as follows [48,49]: Electronics 2022, 11, 815 5 of 27 P pv (t) = 1 0 P pv0 (s)ρ(s)ds (11) where N is the module number; s is the solar irradiance kW/m 2 ; K i and K v are current and voltage temperature coefficients (A/ • C and V/ • C), respectively; T cy and T A are cell and ambient temperatures ( • C), respectively; N OT is the nominal operating temperature of the cell in • C; FF is the fill factor; V oc and I sc are the open circuit voltage (V) and short circuit current (A), respectively; V MPP and I MPP are the voltage and current at the maximum power point, respectively; f b (s) is the Beta-PDF of s; α and β are the parameters of the Beta-PDF; µ and σ are the mean and standard deviation of the random variable s, respectively; s 1 and s 2 are the solar irradiance limits of state (s); and ρ(s) is the probability of the solar irradiance state (s) during any specific hour. Figure 1 describes in detail the calculation of the total average power at each hour including the uncertainties of solar irradiance. As shown, the related uncertainties are modeled using Beta-PDF where the day is split into 24-h periods, each of which is 1 hr. From the collected historical data, the mean and standard deviation of the hourly solar irradiance of the day is estimated. For each hour, different states of the solar irradiance among the Beta-PDF with equal steps are taken. In this study, each hour has 20 states for solar irradiance with a step of 0.05 kW/m 2 . Accordingly, the PV output power is obtained for each state using Equation (3). Besides, the probability of the solar irradiance state (s) is estimated using Equation (10). Thus, the average output power of the PV module at any specific hour can be obtained using Equation (11) [48]. This study considers that a PV unit is associated with the type of converter that can deliver active power only (i.e., unity power factor) as the standard IEEE 1547 [50]. (11) where N is the module number ; s is the solar irradiance kW/m 2 ; Ki and Kv are current and voltage temperature coefficients (A/°C and V/°C), respectively; Tcy and TA are cell and ambient temperatures (°C), respectively; NOT is the nominal operating temperature of the cell in °C; FF is the fill factor; Voc and Isc are the open circuit voltage (V) and short circuit current (A), respectively; VMPP and IMPP are the voltage and current at the maximum power point, respectively; fb(s) is the Beta-PDF of s; α and β are the parameters of the Beta-PDF; µ and σ are the mean and standard deviation of the random variable s, respectively; s1 and s2 are the solar irradiance limits of state (s); and ρ(s) is the probability of the solar irradiance state (s) during any specific hour. Figure 1 describes in detail the calculation of the total average power at each hour including the uncertainties of solar irradiance. As shown, the related uncertainties are modeled using Beta-PDF where the day is split into 24-h periods, each of which is 1 hr. From the collected historical data, the mean and standard deviation of the hourly solar irradiance of the day is estimated. For each hour, different states of the solar irradiance among the Beta-PDF with equal steps are taken. In this study, each hour has 20 states for solar irradiance with a step of 0.05 kW/m 2 . Accordingly, the PV output power is obtained for each state using Equation (3). Besides, the probability of the solar irradiance state (s) is estimated using Equation (10). Thus, the average output power of the PV module at any specific hour can be obtained using Equation (11) [48]. This study considers that a PV unit is associated with the type of converter that can deliver active power only (i.e., unity power factor) as the standard IEEE 1547 [50].

Wind Turbine DGs
The WT-DGs are renewable sources, as well, where their produced powers are intermittent and varied in each hour with high uncertainty levels due to the dependency on wind speed. Generally, the power output from WT is calculated as Equation (12) where v co and v ci are the cut-out and cut-in wind speeds of Wind Turbines (WT), respectively; v is the average wind speed of the hour; v r is nominal operating wind speed of WT; and P r is the maximum power generated by the WT. For the WT-DGs, the uncertainty of wind speed is modeled by the Weibull-PDF. The Weibull-PDF f (v) is formulated as in Equations (13)- (15) [52,53]: where K and C are the shape and scale indexes of the Weibull-PDF; v m is the mean wind speed; and σ is the standard deviation. Based on that, in each hour, different states of the wind speeds are generated. Accordingly, the output power of the WT-DG is obtained for each state using Equation (12). Besides, the probability of the wind speed for each state is estimated using Equation (16). (16) where ρ(v) is the probability of wind speed in each state whereas vw 1 and vw 2 are the regarding limits of the wind speed. Thus, the average output power of the WT-DG for each hour is calculated using Equation (17) [48]. Figure 2 describes in detail the calculation of the total average power at each hour including the uncertainties of wind speed, where each hour has 20 states for wind speed with a step of 5% of the maximum wind speed.

Biomass DGs
The biomass generators (BGs) are firm generation DGs. The resulted BGs power is constant with no associated uncertainties in its rated value. The BG active and reactive powers at bus i can be expressed as in [36]: where Q BG,i and P BG,i are the reactive and active powers associated with BG unit i, respectively. PF BG,i refers to the lag/lead operating biomass power factor, and a i takes positive or negative signs for the BG when supplying and consuming reactive powers, respectively.

Biomass DGs
The biomass generators (BGs) are firm generation DGs. The resulted BGs power is constant with no associated uncertainties in its rated value. The BG active and reactive powers at bus i can be expressed as in [36]: where QBG,i and PBG,i are the reactive and active powers associated with BG unit i, respectively. PFBG,i refers to the lag/lead operating biomass power factor, and ai takes positive or negative signs for the BG when supplying and consuming reactive powers, respectively.

Techno-Economic and Environmental Optimization Scheme of MGs
The optimal operation of a typical MG is investigated to reduce the economical operational costs and the accompanied pollutants simultaneously over a 24-h scheduling horizon.

Objective Function
Nowadays, the environmental emission reduction comes with great importance in order to reduce emission levels [54]. In this work, the goal is to optimize two conflicting cost functions simultaneously, which are the operating cost and the emissions of the WT, PV, BG, and grid. The objective function is modeled as follows: where OCt and Et are the total operational and emission costs in the MG for each hour (t), which are functions of the output powers of the renewable sources and grid; PBG is the output of BGs; Pw is the scheduled output of WT units; Ppv is the scheduled output of PV units; Pgrid is the grid output power; and H is the price penalty factor ($/Ton Co2). The price

Techno-Economic and Environmental Optimization Scheme of MGs
The optimal operation of a typical MG is investigated to reduce the economical operational costs and the accompanied pollutants simultaneously over a 24-h scheduling horizon.

Objective Function
Nowadays, the environmental emission reduction comes with great importance in order to reduce emission levels [54]. In this work, the goal is to optimize two conflicting cost functions simultaneously, which are the operating cost and the emissions of the WT, PV, BG, and grid. The objective function is modeled as follows: Min OC t (P BDG , P w , P pv , P grid ) + H. E t (P BDG , P w , P pv , P grid ) (20) where OC t and E t are the total operational and emission costs in the MG for each hour (t), which are functions of the output powers of the renewable sources and grid; P BG is the output of BGs; P w is the scheduled output of WT units; P pv is the scheduled output of PV units; P grid is the grid output power; and H is the price penalty factor ($/Ton Co 2 ). The price penalty factor was used [55], which is the ratio between the maximum fuel cost and the maximum emission of the associated generator The total operational cost function is expressed in Equation (21) as The total emission cost function can be expressed in Equation (22) as follows: where N w is the number of WT units; N pv is the number of PV units; N BG is the number of BG units; C pw is the operational costs of WT ($); C pv is the operational costs of PV units ($); C BG is the cost coefficients of the operation and maintenance of BG; C grid is the cost coefficients related to grid; and E w , E pv , E BG , and E grid are the emission coefficients related to WT, PV, BG, and the grid, respectively.

Operational Costs
First, for the WT-DGs, the operational cost is the summation of two parts. The first part refers to the cost of the actual output power of the WT. The second part is the penalty cost of unavailable wind output power that refers to the deviation between the available and scheduled wind powers. At time t, the total operational costs of each WT (i) can be represented, as in [56], as where C pw,i,t is the total operating cost of the wind, c w is the operational cost of the wind output ($/kW), c pw is the penalty costs for under-estimation power of the WT ($/kW), P wav,i,t is the available ith WT power at time t (kW), and c rw is the reserve cost for overestimation wind turbine power ($/kW). Second, for the PV units, the operational costs for each unit (i) at each time (t) are described as in [56]: where C pv,i,t is the total operating cost of the PV, c pv is the operational cost of the PV units ($/kW), c ppv and c rpv are, respectively, the penalty costs for under-estimation and overestimated powers ($/kW), and P pv_av,i,t is the available power that can be produced (kW). The first parts of Equations (23) and (24) represent the cost based on the actual PV and WT generation power, while the second part represents the costs due to under-estimation or over-estimation. The under-estimation costs are related to penalties due to the generated PV, and WT output power is less than the scheduled one, whereas the over-estimation costs are related to the reserve costs in case the scheduled output power is greater than the available generated power for both PV and WT.

Constraints
For the optimal operating strategy in MGs, different equality and inequality constraints have to be maintained as presented in Equations (25)- (31). Equations (25) and (26) represent the active and reactive power balance constraints mentioned as P s and Q s that are the supplied active and reactive power by the main feeder and the active and reactive power losses denoted by P Tloss and Q Tloss in the MG, respectively [57,58]. Equations (27)- (29) present the allowable limits of the output powers of different RERs. Equation (30) represents the output power for BG; PV and WT units must follow their operating limits. The power factor (PF) limits for each BG are maintained as in [59]. The bus voltage (V) limits should be in the range of [0.95-1.05] for each bus j [59], as presented in Equation (31). Equation (32) preserves the thermal capacity of the branches below their maximum thermal capacity for each branch [60,61]. Equation (33) gives the penetration bounds of the total renewable sources capacity in the system (KP) as [62,63]:

Equilibrium Optimization for Optimal Operation Strategy in MGs EO Algorithm
The dynamic balance of mass on the control volume is a key source of inspiration for the EO technique. The following three steps can be used to explain its mathematical model: Step 1: Initialization: In the starting of the optimization process, the EO randomly generates the population. The initial concentrations are calculated using uniform random initialization based on the particle number, population, and dimensions, as in Equation (34) [64]: where C i,initial refers to the initial vector for each particle (i), C i,min and C i,max are the lower and upper bounds of the control variables, rand i refers to a random distributed vector in [0,1], and N i is the number of particles. In this step, the fitness of the initialized particles is estimated, and the best scores are used to find the nominee solutions.
Step 2: Equilibrium pool and candidates: In this step, the EO finds the particle's equilibrium state. The algorithm reaches a near-optimal solution at its equilibrium state. It assigns the best four particles in the population at equilibrium candidates and the fifth one containing the average of the previous best four particles. The pool of equilibrium (Ceq, pool) that helps the EO features in their exploitation and exploration operations is expressed by these five equilibrium candidates.
Step 3: Updating the concentration: The evaluation process for updating each concentration vector (C) is carried out as where C eq is a randomly generated vector (C eq, pool ) from the pool of equilibrium; λ is a random vector [0,1]; G is the generation rate; V stands for the volume unit, which is equal to one [31]; and F is an exponential term that helps the EO algorithm in achieving a balance between the exploration and extraction phases. It can be determined as follows: where t 0 is the initial start time, and the time (t) depends on the number of iterations (Iter) as follows: where Iter and Max_iter are the initial and maximum iteration numbers, respectively; and a 2 is a constant value equal to 1 that is used to monitor exploitation potential [31]. The following formula can be considered to boost the developed technique's exploration and exploitation abilities.
where a 1 is a constant value of 2 that is used to control exploration ability [31], r is a random vector in the range of 0 to 1, and the term (sign (r − 0.5)) affects exploration and exploitation directions. The generation rate (G) is calculated as follows: where G cp is the control parameter of the generation rate that is used to update the EO technique, G P is the generation probability that equals 0.5 [31], and r 1 and r 2 are random numbers in the range [0,1]. Figure 3 shows the EO based optimal operation procedure of MG in the tested distribution systems.

Simulation Results
The EO is applied to determine the optimal operation of MG to achieve technical and economic benefits with respecting the associated operational and emission costs for two distribution systems. For the optimal operation, the control variables are the active power resulted from different RERs, including BGs, wind, and PV, and the power factor of BGs.

Simulation Results
The EO is applied to determine the optimal operation of MG to achieve technical and economic benefits with respecting the associated operational and emission costs for two distribution systems. For the optimal operation, the control variables are the active power resulted from different RERs, including BGs, wind, and PV, and the power factor of BGs.

Test Systems
Two test systems are considered to apply the proposed operation strategy. The 1st test system is the 12.66 kV IEEE 33-bus distribution system. It contains 33-node and 32 branches whilst its total active and reactive demand is 3715 kW + 2300 kVAr [65]. Figure 4 displays the first MG configuration where this system includes 1.1 MW-WT connected at bus 3, and five PV generators with ratings of 20, 25, 30, 40, and 50 kW are connected at buses 13, 17, 20, 27, and 33, respectively. In addition to that, four BGs with ratings of 400, 500, 700, and 1000 kVA are installed at buses 15, 25, 9, and 31, respectively [36].  The 12.47 kV AES-Venezuela in the metropolitan area of Caracas that is modeled b 141-node is considered the large-scale test system. The total peak load is 12.19 MW 6.2894 MVAR [66,67]. Figure 5 displays the second MG configuration where the 141-bu system is modified to include one WT with a capacity of 3 MW connected at bus 50, an five PV generators with ratings of 40, 50, 60, 80, and 100 kW are connected at buses 60, 70 80, 90, and 100, respectively. Added to that, four BGs with ratings of 650, 2500, 3500, an 2000 kVA are installed at buses 109, 16, 78, and 63, respectively [36]. For both studie systems, the load profile is taken into consideration where Figure 6 describes this hourl loading as a percentage of the peak load [68]. The 12.47 kV AES-Venezuela in the metropolitan area of Caracas that is modeled by 141-node is considered the large-scale test system. The total peak load is 12.19 MW + 6.2894 MVAR [66,67]. Figure 5 displays the second MG configuration where the 141-bus system is modified to include one WT with a capacity of 3 MW connected at bus 50, and five PV generators with ratings of 40, 50, 60, 80, and 100 kW are connected at buses 60, 70, 80, 90, and 100, respectively. Added to that, four BGs with ratings of 650, 2500, 3500, and 2000 kVA are installed at buses 109, 16, 78, and 63, respectively [36]. For both studied systems, the load profile is taken into consideration where Figure 6 describes this hourly loading as a percentage of the peak load [68].
6.2894 MVAR [66,67]. Figure 5 displays the second MG configuration where the 141-bus system is modified to include one WT with a capacity of 3 MW connected at bus 50, and five PV generators with ratings of 40, 50, 60, 80, and 100 kW are connected at buses 60, 70, 80, 90, and 100, respectively. Added to that, four BGs with ratings of 650, 2500, 3500, and 2000 kVA are installed at buses 109, 16, 78, and 63, respectively [36]. For both studied systems, the load profile is taken into consideration where Figure 6 describes this hourly loading as a percentage of the peak load [68].  Also, the maximum iterations number equals 400 for the tested systems. The particles number is considered to be 50 for IEEE 33-bus and 100 for the large scale 141-bus system. The numbers of iterations and particles are considered the same for all competitive algorithms for a fair assessment of their performance. Additionally, the allowable number of the RERs units, BGs, wind, and PV are assumed as 4, 1, and 5, respectively. The maximum penetration level is limited by 60% of the load demand [69]. Table 1 shows the data of operational costs and the emission for the shared energy sources in the microgrid.

Cases Studied
In this paper, three different cases studied are performed for each system, and these cases are described as follows: Case 1: this is considered the basic system condition without adding DGs units. Case 2: in this case, the MG is operated without considering the optimal control settings of BGs, WT, and PV with daily load variations.
Case 3: this case develops the optimal operation of the MG with optimal control set- Also, the maximum iterations number equals 400 for the tested systems. The particles number is considered to be 50 for IEEE 33-bus and 100 for the large scale 141-bus system. The numbers of iterations and particles are considered the same for all competitive algorithms for a fair assessment of their performance. Additionally, the allowable number of the RERs units, BGs, wind, and PV are assumed as 4, 1, and 5, respectively. The maximum penetration level is limited by 60% of the load demand [69]. Table 1 shows the data of operational costs and the emission for the shared energy sources in the microgrid. PV E pv = 0.045 [70] c pv = 6 [71] WT E w = 0.016 [72] c w = 7 [71] Price penalty factor H = 10 $/Ton Co 2 [73]

Cases Studied
In this paper, three different cases studied are performed for each system, and these cases are described as follows: Case 1: this is considered the basic system condition without adding DGs units. Case 2: in this case, the MG is operated without considering the optimal control settings of BGs, WT, and PV with daily load variations.
Case 3: this case develops the optimal operation of the MG with optimal control settings by using the EO algorithm for BGs, WT, and PV, simultaneously with considering the daily load variations. In the second case, MG is operated without optimal control on DGs whereas the WTs and PVs generate its available power, and the BGs generate its full capacity at a unity power factor. For the small-scale PV DGs, the Beta-PDF is used to model the uncertainty of the solar irradiance and the evaluation of PV output power is performed as described in Figure 1 considering the uncertainties in solar irradiance. From the collected historical data, the mean and standard deviations of the hourly solar irradiance of the day are estimated. Based on that, the Beta-PDF is generated for each hour. Consequently, 20 different states of solar irradiance among the Beta-PDF are taken with equal step sizes of 0.05 kW/m 2 for each hour. Figure 7 displays their produced powers from each PV source at each hour. There is no power output to be generated from the PV sources in the first and last five hours in the morning and evening since there are no solar irradiances at these hours. Additionally, the highest output of each PV source occurred at 12 pm, with 55.2, 55.3%, 55.4%, 55.2% and 55.3% from their installed capacities for the PV sources at buses 13, 17, 20, 27, and 33, respectively.  For the WT-DGs, the Weibull-PDF is used to model the uncertainty of the wind speed and the evaluation of WT output power, considering that the uncertainties in wind speed is performed as described in Figure 2. For each hour, the mean and standard deviations of the hourly wind speed of the day are estimated, and the associated Weibull-PDF is generated. Similarly, 20 states of the wind speed, for each hour, are considered where the step is 1.25 m/s. Figure 8 displays the hourly active and reactive powers that are produced from the WT source at bus 3. The type 4 distributed generation, WT-DG, is used in this For the WT-DGs, the Weibull-PDF is used to model the uncertainty of the wind speed and the evaluation of WT output power, considering that the uncertainties in wind speed is performed as described in Figure 2. For each hour, the mean and standard deviations of the hourly wind speed of the day are estimated, and the associated Weibull-PDF is generated. Similarly, 20 states of the wind speed, for each hour, are considered where the step is 1.25 m/s. Figure 8 displays the hourly active and reactive powers that are produced from the WT source at bus 3. The type 4 distributed generation, WT-DG, is used in this study to deliver the active power and to consume reactive power using the induction generators at a fixed speed. At hour 18, the WT-DG produces the highest active power with 692.1 kW and consumes the highest reactive power with 69.16 kVAr. Thus, it is operated with 62.92% from its full capacity. At hour 7, the WT-DG is operated with 36.53% from its full capacity, where it produces the least active power with 401.8 kW and consumes the least reactive power with 56.45 kVAr. For the WT-DGs, the Weibull-PDF is used to model the uncertainty of the wind speed and the evaluation of WT output power, considering that the uncertainties in wind speed is performed as described in Figure 2. For each hour, the mean and standard deviations of the hourly wind speed of the day are estimated, and the associated Weibull-PDF is generated. Similarly, 20 states of the wind speed, for each hour, are considered where the step is 1.25 m/s. Figure 8 displays the hourly active and reactive powers that are produced from the WT source at bus 3. The type 4 distributed generation, WT-DG, is used in this study to deliver the active power and to consume reactive power using the induction generators at a fixed speed. At hour 18, the WT-DG produces the highest active power with 692.1 kW and consumes the highest reactive power with 69.16 kVAr. Thus, it is operated with 62.92% from its full capacity. At hour 7, the WT-DG is operated with 36.53% from its full capacity, where it produces the least active power with 401.8 kW and consumes the least reactive power with 56.45 kVAr.  On the other side, the BGs are operated at full capacity without optimal control. By running the load flow algorithm for these hourly circumstances, Figure 9 shows the total active power of each type of DG and the summation of loads and the hourly power losses. Then, Figure 9 concludes the following findings:

•
During nine hours (1-8 and 24), the sign of the grid power is negative. At these hours, the DGs generated power in the MG is greater than the loads and the MG power losses so that the excess power is back to the grid bus. Consequently, there is surplus active power to be generated from the DGs in the MG and injected to the grid bus. On the other side, the BGs are operated at full capacity without optimal control. By running the load flow algorithm for these hourly circumstances, Figure 9 shows the total active power of each type of DG and the summation of loads and the hourly power losses. Then, Figure 9 concludes the following findings:  In the third case, the EO algorithm is applied to operate the economical solution of MG considering the emission concerns by simultaneously minimizing the hourly operational costs and the accompanied pollutants. Added to that, the assessment of EO perfor-

Case 3: MG with Optimal Control on BGs, WT, and PV Simultaneously with Daily Load Variation
In the third case, the EO algorithm is applied to operate the economical solution of MG considering the emission concerns by simultaneously minimizing the hourly operational costs and the accompanied pollutants. Added to that, the assessment of EO performance is carried out compared with other optimization algorithms such as DE [74] and RAO [75] techniques at different hours. The convergence comparisons between the developed EO, DE, and RAO algorithms at hours 6 and 12 are described in Figure 10 where the associated outcomes of them are tabulated in Tables 2 and 3 at hours 6 and 12, respectively. At hour 6 (  Based on the developed EO, the output powers of the PV, WT, and BGs are optimized besides the associated power factors of the BGs. Figure 11a shows the percentage apparent power of the BGs at buses 15, 25, 9, and 31 for each hour, whereas Figure 11b displays the hourly optimized value of the power factor. From Figure 11a,b:

•
The first BG at bus 31 is operating within range [40:100] % of its capacity to inject apparent power demand to the MG during the day. Similarly, the second BG at bus 15 is operating within range [60:100] % of its capacity during the day excluding hour 6. • Both BG1 and BG2 at buses 31 and 15 are approximately operating at a power factor of 0.85, which has the minimum power factor to be considered, and so their injected reactive power is directly proportional to their scheduled output power. • Both BG3 and BG4 at buses 25 and 9 produce very small apparent power at the first five hours, respectively, which has the smallest loadings. Then, their output powers increase through the loading increases, whereas BG3 at bus 25 is operating with 100% of its capacity at hour 6. Added to that, BG4 at bus 9 is operating with 100% of its capacity at hours 12, 14, 17, and 18. • BG3 at bus 25 is operating at a varied power factor through the day. BG4 at bus 9 is operating at a varied power factor at the first nine hours with low loading values, whilst its power factor is nearly fixed at 0.85 for the next hours.
(a) Percentage apparent power   Based on the developed EO, the output powers of the PV, WT, and BGs are optimized besides the associated power factors of the BGs. Figure 11a shows the percentage apparent power of the BGs at buses 15, 25, 9, and 31 for each hour, whereas Figure 11b displays the hourly optimized value of the power factor. From Figure 11a,b:

•
The first BG at bus 31 is operating within range [40:100] % of its capacity to inject apparent power demand to the MG during the day. Similarly, the second BG at bus 15 is operating within range [60:100] % of its capacity during the day excluding hour 6. • Both BG1 and BG2 at buses 31 and 15 are approximately operating at a power factor of 0.85, which has the minimum power factor to be considered, and so their injected reactive power is directly proportional to their scheduled output power. • Both BG3 and BG4 at buses 25 and 9 produce very small apparent power at the first five hours, respectively, which has the smallest loadings. Then, their output powers increase through the loading increases, whereas BG3 at bus 25 is operating with 100% of its capacity at hour 6. Added to that, BG4 at bus 9 is operating with 100% of its capacity at hours 12, 14, 17, and 18. • BG3 at bus 25 is operating at a varied power factor through the day. BG4 at bus 9 is operating at a varied power factor at the first nine hours with low loading values, whilst its power factor is nearly fixed at 0.85 for the next hours.
five hours, respectively, which has the smallest loadings. Then, their output powers increase through the loading increases, whereas BG3 at bus 25 is operating with 100% of its capacity at hour 6. Added to that, BG4 at bus 9 is operating with 100% of its capacity at hours 12, 14, 17, and 18. • BG3 at bus 25 is operating at a varied power factor through the day. BG4 at bus 9 is operating at a varied power factor at the first nine hours with low loading values, whilst its power factor is nearly fixed at 0.85 for the next hours.
(a) Percentage apparent power (b) Power factor Figure 11. Hourly Percentage apparent power and power factor for each BG for the first MG system. Figure 11. Hourly Percentage apparent power and power factor for each BG for the first MG system.
To study the hourly power balance operation of the MG, Figure 12 illustrates the total active power for each type of DG and the summation of loads and hourly power losses. This figure leads to the following findings:

•
The total generation equals to the loads and power losses at each hour, where the summations of the output powers from the BGs are varied, which is optimally controlled by the developed EO algorithm at each hour.

•
The developed EO algorithm optimizes the output powers from WT and PV sources at their available outputs to make use of their environmentally friendly characteristics.

•
During the day, there are no excess powers in the MG back to the grid bus. The MG was supplied with the active power from the grid bus. The maximum supplied power is achieved at hour 12 with 1520 kW, whilst the minimum supplied power is reached at hour 5 with 843.97 kW. • As a result, the supplied power represents 40.58% and 40.34% referred to the loading and power losses at hours 12 and 5, respectively. Therefore, the total produces power from the DGs with 59.42% and 59.66% referred to the loading and power losses at hours 12 and 5, respectively. • Nevertheless, the penetration of the total DGs power production in the MG does not exceed the penetration limit of 60%.
To assess the voltage quality of the MG operation for each hour, Figure 13 displays the minimum voltage at each hour for the three cases studied. It is shown that the minimum voltage is corrected in cases 2 and 3 where, when the voltages at all hours are below the permissible limit of 0.95 Pu in the case (case 1), the highest minimum voltage occurs at bus 18 at hour 5, whereas the least minimum voltage occurs at the same bus at hour 15. At both hours, the voltage profiles at all MG buses are described in Figure 14. The voltage profile at each bus is improved at light loading at hour 5. At this loading hour, the minimum voltage at bus 18 is corrected from 0.9495 to 0.9877 Pu, which exceeds the minimum limit of 0.95 Pu in case 3; consequently, this improvement represents 4.02%. Additionally, the voltage profile at each bus is improved at peak loading at hour 15, where the minimum voltage at bus 18 is improved with 7.822% and corrects the voltage from 0.9038 to 0.9745 Pu, which exceeds the minimum limit of 0.95 Pu in case 3.
reached at hour 5 with 843.97 kW.
• As a result, the supplied power represents 40.58% and 40.34% referred to the loading and power losses at hours 12 and 5, respectively. Therefore, the total produces power from the DGs with 59.42% and 59.66% referred to the loading and power losses at hours 12 and 5, respectively. • Nevertheless, the penetration of the total DGs power production in the MG does not exceed the penetration limit of 60%. To assess the voltage quality of the MG operation for each hour, Figure 13 displays the minimum voltage at each hour for the three cases studied. It is shown that the minimum voltage is corrected in cases 2 and 3 where, when the voltages at all hours are below the permissible limit of 0.95 Pu in the case (case 1), the highest minimum voltage occurs at bus 18 at hour 5, whereas the least minimum voltage occurs at the same bus at hour 15. At both hours, the voltage profiles at all MG buses are described in Figure 14. The voltage profile at each bus is improved at light loading at hour 5. At this loading hour, the minimum voltage at bus 18 is corrected from 0.9495 to 0.9877 Pu, which exceeds the minimum limit of 0.95 Pu in case 3; consequently, this improvement represents 4.02%. Additionally, the voltage profile at each bus is improved at peak loading at hour 15, where the minimum voltage at bus 18 is improved with 7.822% and corrects the voltage from 0.9038 to 0.9745 Pu, which exceeds the minimum limit of 0.95 Pu in case 3. In addition, Figure 15 shows the hourly active power losses for the three cases that are greatly reduced from case 1 to cases 2 and 3, whereas the optimal operating strategy based on the developed EO algorithm in case 3 provides the minimum power losses at In addition, Figure 15 shows the hourly active power losses for the three cases that are greatly reduced from case 1 to cases 2 and 3, whereas the optimal operating strategy based on the developed EO algorithm in case 3 provides the minimum power losses at each hour through the day. Compared to case 1, the percentages of the reduction in power losses that are achieved by case 3 reach 81.9, 79.6, 73.7, 80.6, 80.9, 74.2, 71.8, 75.5, 84.9, 83.7, In addition, Figure 15 shows the hourly active power losses for the three cases that are greatly reduced from case 1 to cases 2 and 3, whereas the optimal operating strategy based on the developed EO algorithm in case 3 provides the minimum power losses at each hour through the day. Compared to case 1, the percentages of the reduction in power losses that are achieved by case 3 reach 81.9, 79.6, 73.7, 80.6, 80.9, 74.2, 71.8, 75.5, 84.9, 83.7, 85.3,  85.2, 85.5, 85.2, 84.2, 84.8, 82.1, 81.1, 82.4, 83.1, 83.7, 84.5, 84.3, and 79.8% for hours 1-24, respectively. Compared to case 2, the percentages of the reduction in power losses that are achieved by case 3 reach 41.9, 48.5, 48.6, 61.7, 61.8, 48.6, 31.6, 19.5, 22.3, 18.7, 19.7, 19.4, 19.9, 19.4, 18.3, 19.4, 16.7, 15.7, 17.6, 20.3, 19.1, 19.7, 21.3, and 27.6% for hours 1-24, respectively. Additionally, the hourly operational costs and the associated emissions of the distributed energy sources in the MG for the cases studied are depicted in Figure 16. As shown, the operational costs and the associated emissions in the MG at each hour are greatly reduced from case 1 to cases 2 and 3. Despite that case 2 provides the least operational costs and the emissions in the MG in comparison to case 3, the penetration limit of the total output of the DGs in the MG exceeds the limits of the 60% penetration ratio, as detailed in Figure 9. In Figure 9

Simulation Results of a large-Scale 141-Bus Test System
In this subsection, the proposed solution methodology is applied on a large-scale 141- Additionally, the hourly operational costs and the associated emissions of the distributed energy sources in the MG for the cases studied are depicted in Figure 16. As shown, the operational costs and the associated emissions in the MG at each hour are greatly reduced from case 1 to cases 2 and 3. Despite that case 2 provides the least operational costs and the emissions in the MG in comparison to case 3, the penetration limit of the total output of the DGs in the MG exceeds the limits of the 60% penetration ratio, as detailed in Figure 9. In Figure 9, the penetration level in case 2 exceeds 100% at nine hours, where it reaches 129. 96  Additionally, the hourly operational costs and the associated emissions of the distributed energy sources in the MG for the cases studied are depicted in Figure 16. As shown, the operational costs and the associated emissions in the MG at each hour are greatly reduced from case 1 to cases 2 and 3. Despite that case 2 provides the least operational costs and the emissions in the MG in comparison to case 3, the penetration limit of the total output of the DGs in the MG exceeds the limits of the 60% penetration ratio, as detailed in Figure 9. In Figure 9, the penetration level in case 2 exceeds 100% at nine hours, where it reaches 129. 96

Simulation Results of a large-Scale 141-Bus Test System
In this subsection, the proposed solution methodology is applied on a large-scale 141bus system. Three cases studied are investigated. The system's initial power loss (without DGs) was 603.821 kW, where bus 87 has the minimum voltage that is 0.9294 Pu. In the 2nd Figure 16. Hourly optimal operational costs and the associated emissions in the MG for the first MG.

Simulation Results of a Large-Scale 141-Bus Test System
In this subsection, the proposed solution methodology is applied on a large-scale 141-bus system. Three cases studied are investigated. The system's initial power loss (without DGs) was 603.821 kW, where bus 87 has the minimum voltage that is 0.9294 Pu. In the 2nd case, the resulted uncertainties of the solar irradiance and wind speed are modeled by Beta-PDF and Weibull-PDF for the PV and WT sources, respectively. Therefore, the PV output power is evaluated at buses 60, 70, 80, 90, and 100. The active and reactive powers are hourly produced from the WT source at bus 50, as displayed in Figure 17.  On the other side, the BGs are operated at their full capacity without optimal control. By running the load flow algorithm for these hourly circumstances, Figure 18 shows the total active power of each type of DG and the summation of loads and losses at each hour. During 9 h (1-8 and 24), there is surplus active power to be generated from the DGs in the MG and injected to the grid bus, where the generated power from the DGs in the MG is greater than the loads and the power losses so that the excess power is back to the grid bus. In the other 15 h, the MG absorbs active power from the grid bus. In the third case, the developed EO algorithm is applied to optimally operate the MG in order to minimize the operational costs and the accompanied pollutants simultaneously over a 24-h scheduling horizon. Based on the developed EO, the output powers of the PV, WT, and BGs are optimized besides the associated power factors of the BGs. Figure  19a shows the percentage apparent power of the BGs at buses 109, 16, 78, and 63 for each hour, whereas Figure 19b displays the hourly optimized value of the power factor. On the other side, the BGs are operated at their full capacity without optimal control. By running the load flow algorithm for these hourly circumstances, Figure 18 shows the total active power of each type of DG and the summation of loads and losses at each hour. During 9 h (1-8 and 24), there is surplus active power to be generated from the DGs in the MG and injected to the grid bus, where the generated power from the DGs in the MG is greater than the loads and the power losses so that the excess power is back to the grid bus. In the other 15 h, the MG absorbs active power from the grid bus.  On the other side, the BGs are operated at their full capacity without optimal control. By running the load flow algorithm for these hourly circumstances, Figure 18 shows the total active power of each type of DG and the summation of loads and losses at each hour. During 9 h (1-8 and 24), there is surplus active power to be generated from the DGs in the MG and injected to the grid bus, where the generated power from the DGs in the MG is greater than the loads and the power losses so that the excess power is back to the grid bus. In the other 15 h, the MG absorbs active power from the grid bus. In the third case, the developed EO algorithm is applied to optimally operate the MG in order to minimize the operational costs and the accompanied pollutants simultaneously over a 24-h scheduling horizon. Based on the developed EO, the output powers of the PV, WT, and BGs are optimized besides the associated power factors of the BGs. Figure  19a shows the percentage apparent power of the BGs at buses 109, 16, 78, and 63 for each hour, whereas Figure 19b displays the hourly optimized value of the power factor. In the third case, the developed EO algorithm is applied to optimally operate the MG in order to minimize the operational costs and the accompanied pollutants simultaneously over a 24-h scheduling horizon. Based on the developed EO, the output powers of the PV, WT, and BGs are optimized besides the associated power factors of the BGs. Figure 19a shows the percentage apparent power of the BGs at buses 109, 16, 78, and 63 for each hour, whereas Figure 19b displays the hourly optimized value of the power factor.  To describe the power balance operation of the MG for each hour, Figure 20 illustrates the total active power of each type of DG and the summation of loads and power losses at each hour. At each hour, the generated grid power is always higher than 40% of summation loads and power losses whereas the penetration level of 60% is preserved by means of the developed EO algorithm. For each hour, Figure 21 displays the minimum voltage for the three-cases studied. For hours 8-23, the minimum voltages are corrected in cases 2 and 3. In addition, in cases 2 and 3, the minimum voltages at all hours are above the permissible limit of 0.95 Pu. In the initial case (case 1), the highest minimum voltage occurs at bus 87 at hour 5, whereas the least minimum voltage occurs at the same bus at hour 15. At both hours, the voltage profile at all MG buses is described in Figure 22. However, the voltage profile at each bus is improved at light loading at hour 5. At this loading hour, bus 87 has the minimum voltage level in case 3. It is declined from 0.9616 to 0.9914 Pu. Additionally, the voltage profile at each bus is corrected at peak loading at hour 15, where the bus 87 has the minimum voltage level in case 3 that is declined from 0.9294 to 0.9772 Pu which exceeds the minimum limit of 0.95 Pu; consequently, this improvement represents 5.14%. To describe the power balance operation of the MG for each hour, Figure 20 illustrates the total active power of each type of DG and the summation of loads and power losses at each hour. At each hour, the generated grid power is always higher than 40% of summation loads and power losses whereas the penetration level of 60% is preserved by means of the developed EO algorithm.  To describe the power balance operation of the MG for each hour, Figure 20 illustrates the total active power of each type of DG and the summation of loads and power losses at each hour. At each hour, the generated grid power is always higher than 40% of summation loads and power losses whereas the penetration level of 60% is preserved by means of the developed EO algorithm. For each hour, Figure 21 displays the minimum voltage for the three-cases studied. For hours 8-23, the minimum voltages are corrected in cases 2 and 3. In addition, in cases 2 and 3, the minimum voltages at all hours are above the permissible limit of 0.95 Pu. In the initial case (case 1), the highest minimum voltage occurs at bus 87 at hour 5, whereas the least minimum voltage occurs at the same bus at hour 15. At both hours, the voltage profile at all MG buses is described in Figure 22. However, the voltage profile at each bus is improved at light loading at hour 5. At this loading hour, bus 87 has the minimum voltage level in case 3. It is declined from 0.9616 to 0.9914 Pu. Additionally, the voltage profile at each bus is corrected at peak loading at hour 15, where the bus 87 has the minimum voltage level in case 3 that is declined from 0.9294 to 0.9772 Pu which exceeds the minimum limit of 0.95 Pu; consequently, this improvement represents 5.14%. For each hour, Figure 21 displays the minimum voltage for the three-cases studied. For hours 8-23, the minimum voltages are corrected in cases 2 and 3. In addition, in cases 2 and 3, the minimum voltages at all hours are above the permissible limit of 0.95 Pu. In the initial case (case 1), the highest minimum voltage occurs at bus 87 at hour 5, whereas the least minimum voltage occurs at the same bus at hour 15. At both hours, the voltage profile at all MG buses is described in Figure 22. However, the voltage profile at each bus is improved at light loading at hour 5. At this loading hour, bus 87 has the minimum voltage level in case 3. It is declined from 0.9616 to 0.9914 Pu. Additionally, the voltage profile at each bus is corrected at peak loading at hour 15, where the bus 87 has the minimum voltage level in case 3 that is declined from 0.9294 to 0.9772 Pu which exceeds the minimum limit of 0.95 Pu; consequently, this improvement represents 5.14%.  Moreover, the active power losses at each hour for the three cases studied are depicted in Figure 23. The power losses at each hour are greatly reduced from case 1 to cases 2 and 3, whereas the optimal operating strategy based on the developed EO algorithm in case 3 provides the minimum power losses at each hour through the day. Compared to case 1, the percentages of the reduction in power losses that are achieved by case 3 reached 78. 9    Moreover, the active power losses at each hour for the three cases studied are depicted in Figure 23. The power losses at each hour are greatly reduced from case 1 to cases 2 and 3, whereas the optimal operating strategy based on the developed EO algorithm in case 3 provides the minimum power losses at each hour through the day. Compared to case 1, the percentages of the reduction in power losses that are achieved by case 3 reached 78. 9  Moreover, the active power losses at each hour for the three cases studied are depicted in Figure 23. The power losses at each hour are greatly reduced from case 1 to cases 2 and 3, whereas the optimal operating strategy based on the developed EO algorithm in case 3 provides the minimum power losses at each hour through the day. Compared to case 1, the percentages of the reduction in power losses that are achieved by case 3 reached 78. 9  Additionally, the operational costs and the associated emissions of the distributed energy sources in the MG at each hour for the three cases studied are depicted in Figure  24. However, the operational costs and the associated emissions in the MG at each hour are greatly reduced from case 1 to cases 3 and 2. Despite case 2 providing the least operational costs and the emissions in the MG in comparison to case 3, the penetration limit of the total output of the DGs in the MG exceeds the limits of the 60% penetration ratio, as described by Figure 18.   Additionally, the operational costs and the associated emissions of the distributed energy sources in the MG at each hour for the three cases studied are depicted in Figure 24. However, the operational costs and the associated emissions in the MG at each hour are greatly reduced from case 1 to cases 3 and 2. Despite case 2 providing the least operational costs and the emissions in the MG in comparison to case 3, the penetration limit of the total output of the DGs in the MG exceeds the limits of the 60% penetration ratio, as described by Figure 18. Additionally, the operational costs and the associated emissions of the distributed energy sources in the MG at each hour for the three cases studied are depicted in Figure  24. However, the operational costs and the associated emissions in the MG at each hour are greatly reduced from case 1 to cases 3 and 2. Despite case 2 providing the least operational costs and the emissions in the MG in comparison to case 3, the penetration limit of the total output of the DGs in the MG exceeds the limits of the 60% penetration ratio, as described by Figure 18.

•
The optimal operating strategy based on the developed equilibrium optimization (EO) technique has been presented and is proposed in order to minimize the operational costs and the accompanied pollutants simultaneously over a 24-h scheduling horizon. • Several renewable DGs have been investigated in the understudy MG, consisting of biomass generators (BG), wind turbines (WTs), and photovoltaics (PV).

•
The uncertainties of wind speed and solar irradiation have been considered via Weibull and Beta-PDF with different states of mean and standard deviation at each hour, respectively.

•
The developed EO methodology has been tested on the IEEE 33-bus system and the practical large-scale 141-bus system of AES-Venezuela in the metropolitan area of Caracas to verify its effectiveness and efficiency. • Added to that, the optimal operation of MG using the developed EO has been compared with other recent algorithms of DE and RAO algorithms. • Based on the developed EO, the hourly output powers of the PV, WTs, and BGs are optimized, as are the associated power factors of the BGs. • Significant benefits of the optimal operation of a typical MG using the developed EO have been demonstrated by minimizing the operational costs and emissions while preserving the penetration level of the DGs by 60%.

•
The voltage profiles of the MG operation for every hour have been greatly improved, with the minimum voltage adjusted within the permitted limit of [0.95-1.05] p.u.

•
The active power losses per hour have been significantly decreased.