Probabilistic Modeling and Equilibrium Optimizer Solving for Energy Management of Renewable Micro-Grids Incorporating Storage Devices

: Recently, micro-grids (MGs) have had a great impact on power system issues due to their clear environmental and economic advantages. This paper proposes an equilibrium optimizer (EO) technique for solving the energy management problem of MGs incorporating energy storage devices concerning the emissions from renewable energy sources (RES) of MGs. Because of the imprecision and uncertainties related to the RESs, market prices, and forecast load demand, the optimization problem is described in a probabilistic manner using a 2m + 1 point estimation approach. Then, the EO approach is utilized for solving the probabilistic energy management (EM) problem. The EM problem is described according to the market policy on the basis of minimizing the total operating cost and emission from RESs through optimal settings of the power generated from distributed generators (DGs) and grids connected under the condition of satisfying the operational constraints of the system. The proposed EO is evaluated based on a grid-connected MG that includes energy storage devices. Moreover, to prove the effectiveness of the EO, it is compared with other recently meta-heuristic techniques. The simulation results show acceptable robustness of the EO for solving the EM problem as compared to other techniques.


Introduction
The exponential increase in universal power demand has had a major impact on the rapid depletion of fossil fuels, which has resulted in increasing gas emissions from conventional generators. To solve this problem, the world has tended to encourage obtaining energy from renewable energy resources (RESs) such as wind turbines (WT), photovoltaics (PV), biomass, fuel cell (FC), micro-turbines (MT), hydropower, tides, and energy storage (ES) devices [1,2]. It is worth noting that micro-grids (MGs) systems have an important vital role in solving different energy-related problems such as reducing the minimization of the operating costs, transmission losses, environmental pollution, and peak load as well as improving voltage regulation and participating in frequency regulation. The MGs systems consist of different technologies of non-dispatchable and dispatchable renewable energy sources (RES), ES devices, and dispatchable loads. Further, MGs enable the active sharing of customers by giving them real-time control and access to the information [3][4][5].
Considering all mentioned benefits, more and more extra small MGs have been deployed over all the world such as facilities, military bases, hospitals, and industries in recent ten years [6]. With the high penetration of energy from RESs in the existing power systems, it also has the important impact of reducing energy losses, pollutant emissions, and enhancing power system performance. However, the improper utilization of RESs leads to problems in the operation of the system [7,8].

•
Formulate the optimization problem of EM incorporating ES devices with consideration of the emissions from RESs via converting the multi-objective function problem into a coefficient single objective function using a price penalty factor and weighting factors by handling the operational constraints.

•
Model the uncertainties of RESs with ES devices, demand load, and market prices, then apply (2m + 1) point estimate method. • Propose EO for solving the optimal problem of EM as an efficient technique and compare with other new techniques that are newly employed here for solving EM problem. • Implement the proposed technique for solving the EM problem based on deterministic and probabilistic EM problems with emission. • Investigate the effectiveness and applicability of the EO when compared with other recent optimization techniques through different scenarios.

The Mathematical Modeling of the EM Problem
A model of low-voltage (LV) MG connected to grid consists of different types of RESs, storage device, and loads as shown in Figure 1. The MG is connected to the utility (MV Energies 2021, 14, 1373 4 of 24 distribution network) by the MV/LV transformer. The MG central controller (MGCC) acts as an interface between the power system and the MG. The MGCC is existing downstream of utility. The functions of MGCC are designed to monitor the real and reactive power produced from DG units, maximize the MG's value, and optimize the MG operations through sending control signal to MGs to evaluate the level of DG's production. The local controllers (LCs) are located close to the DGs. The functions of LCs are assigned to control the MG's components through adapting each type of DG units, controllable loads, and storage devices. Additionally, the LCs are used to control both the voltage and frequency in an islanding operating mode of MG. The optimization proceedings applied in the MGCC are based on the implemented market policy in MG operation. Most policies of operating MGs are the MGCC goal to serve the demand load of the MG through maximizing utilization of its local production. From the consumers' point of view, the MGCC aims to minimize the total operational cost of the MG with considering DG bids, market price, and demand load through optimal adjustment of the generated power from DG units while satisfying the operational constraints [10].
(MV distribution network) by the MV/LV transformer. The MG central controller (MGCC) acts as an interface between the power system and the MG. The MGCC is existing downstream of utility. The functions of MGCC are designed to monitor the real and reactive power produced from DG units, maximize the MG's value, and optimize the MG operations through sending control signal to MGs to evaluate the level of DG's production. The local controllers (LCs) are located close to the DGs. The functions of LCs are assigned to control the MG's components through adapting each type of DG units, controllable loads, and storage devices. Additionally, the LCs are used to control both the voltage and frequency in an islanding operating mode of MG. The optimization proceedings applied in the MGCC are based on the implemented market policy in MG operation. Most policies of operating MGs are the MGCC goal to serve the demand load of the MG through maximizing utilization of its local production. From the consumers' point of view, the MGCC aims to minimize the total operational cost of the MG with considering DG bids, market price, and demand load through optimal adjustment of the generated power from DG units while satisfying the operational constraints [10].

The Operating Cost Function
The total operating cost of the MG system contains several DG units bids and the power exchange cost between the main grid and MG. The control variables of the optimization problem consist of the real power from the DG units including the storage units, whilst the dependent variable is the real power that is sold to or bought from the main grid. Thus, the mathematical formulation of the EM problem may be described as follows [7]:

The Operating Cost Function
The total operating cost of the MG system contains several DG units bids and the power exchange cost between the main grid and MG. The control variables of the optimization problem consist of the real power from the DG units including the storage units, whilst the dependent variable is the real power that is sold to or bought from the main grid. Thus, the mathematical formulation of the EM problem may be described as follows [7]: where R indicates the optimal output power of different DG units, main grid, and their binary status. It is defined by: R t = P G1,t , P G2,t , . . . , P GNG,t , P S1,t , P S2,t , . . . , P SNS,t , P grid,t , u 1,t , u 2,t , . . . , u NG+NS,t From (1), it is seen that the objective function includes the total operating cost of the DG units along with their costs of startup/shutdown in the first term. The second term represents the operating cost of the ES devices with their costs of startup/shut-down. The third term is the exchange power cost between the main grid and MG.

The Pollutant Emission
The minimization of the pollutant emission includes the amount of emissions from the main grid, DG units, and ES devices in MG. There are different pollutant emissions; however, the considered important pollutants which are used to evaluate the total emission are CO 2 , SO 2 , and NO x . Minimizing these emissions is considered as the second objective function that can be formulated as follows [30]: Not only does the objective function include the amount of emission from the main grid, but it also considers the emissions produced from the DG units as well as the ES devices as stated in (4). The total amount of pollutants emissions are E grid,t , E Gi,t , and E Sk,t , that can be described as follows [30]: The total generated power from DG units, ES devices, and the main grid must supply the total load demand of the MG at each period t. However, the active power losses are neglected in the MG. Accordingly, the constraints of the power balance equation can be expressed as follows [31]:

Constraints of Power Generation
The energy generated from each DG unit, ES devices, and the main grid should be limited in the range of its defined capacity.

P min
Gi,t ≤ P Gi,t ≤ P max P min Sk,t ≤ P Sk,t ≤ P max P min grid,t ≤ P grid,t ≤ P max grid,t

Spinning Reserve
It is necessary to maintain the power system reliability due to the continuous fluctuations in RESs. However, for achieving the spinning reserve, the following expression should be employed [10].

Constraints of Energy Storage Devices
Concerning ES devices, the amount of stored energy and the limits on the charge and discharge rate of these ES devices during each time t interval with the constraints are defined as [10]: E es,t = E es,t−1 +η chr P chr ∆t − (1/η dischr )P dischr ∆t

The Optimization Problem
The total demand load is supplied from the sources of MG by the local production, as available. Therefore, the output power from the main grid P grid,t is set as a dependent variable for enforcing the power balance equation constraint. However, the power from the main grid is determined from (8) as follows: The constraints can be handled by converting it into unconstrained one using the penalty function technique. However, the penalty factors with the constraints are involved into the objective function as quadratic penalty terms. Further, the goal of the penalty factors that are chosen as positive large values is to reduce the impractical solutions [32]. Therefore, P grid,t is added as a quadratic penalty term to the objective function in (1). Consequently, the new extended objective function by the penalty factor becomes: MP grid,t P grid,t +γ P grid,t − P lim grid,t 2 (18) In this paper, the above two objective functions: the minimization of both cost and pollution emission of MGs, are integrated. The multi-objective function is converted into a single objective function using the coefficients and price penalty terms. Moreover, the price penalty function for each generating unit is defined as the ratio between the operating cost and the pollution emission. There are various types of price penalty techniques as discussed in [33]. Here, the max/max price penalty factor is utilized. Consequently, the developed objective function can be described as follows:

Probabilistic EM of the MG
The probabilistic formulation problem requires representation of the input random variables, statistically. Solving the EM problem based on the probabilistic approach has two key steps. In the first step, the input random variables are characterized statistically. The source of the uncertainty for the dependent variables in the system arise from the uncertainty of the independent input random variables. Furthermore, due to the random nature of input variables such as the variations of demand load, output power from DG units, and market prices, the obtained EM results such as the total operating cost of the system should be considered as a random variable. Evaluation of the statistical characteristics for the EM result is the second step in the probabilistic EM approach. However, in the case of one or more uncertain input variables, the EM problem turns into a probabilistic problem.

The Statistical Characterization of Input Random Variables
Evaluating the output power of WT mainly depends on two factors: the power curve of WT and the wind speed at a specific location, as described in [34] as follows: Likewise, the output power of the PV depends on both the solar irradiance and the ambient air temperature of the PV. The output power of PV can be computed as described in [35] as follows: where, the temperature of PV module can be determined as follows:

Modeling of Wind Speed
At a specific location, the probability density function (PDF) of the wind speed, as well as the output power of WT, can be described by a Weibull distribution as in [34] as follows: For the Weibull distribution, the cumulative density function (CDF) is: For calculating the wind speed, the CDF as well as its inverse has been used as follows: The parameters of the Weibull distribution can be calculated as follows: where the gamma function Γx is described as: 3.1.2. Modeling of Solar Irradiance, Market-Price, and Load Demand In this paper, it was suggested that the solar irradiance, market price, and demand of load have a normal distribution function. Therefore, the corresponding PDF of any variable z i can be formulated as follows [10,35]: Additionally, the CDF for the normal distribution can be computed by: For determining the variable z i , the CDF and its inverse has been used as follows: where, the error function and its inverse are defined as follows:

Statistical Evaluation of the Output
In order to formulate the probabilistic problem, it is necessary to characterize the input random variables in a statistical nature as well as the method for evaluating the output variables statistically. The probabilistic EM can be formulated mathematically as follows: There are several statistical approaches for handling the output variables as introduced in some literature, in which they are classified into three categories: analytical methods, approximation approaches, and the Monte Carlo simulation (MCS). The MCS randomly generates many uncertain input variables, then uses these values to solve a deterministic problem. The simulation process is repeated many times (hundreds to thousands) to attain an acceptable accuracy for the statistical characteristics of the output results. The major problem of using MCS is the large number of simulations that are required to achieve convergence.
The point estimation approach is considered one of the approximation approaches that present a convergent description of the statistical output random variables. This method, like MCS, employs a deterministic approach for solving a probabilistic problem but over a much lower number of the simulations that is implemented in MCS. Moreover, this method guarantees a considerable reduction in the calculation's efforts compared to the MCS.
The basic idea of the point estimation approach is to focus on the information collected from few forecasted values of the m input variables on K = 2 points for each interest input variable that are called central moments. Using these two central moments for each input, the output variables are calculated in the form of statistical moments 2m + 1 times for selected uncertain input variables [10]. Depending on the statistical moments, the PDF of output random variables of interest can be approximated based on the Gram-Charlier series approach [36]. The following steps of the 2m + 1 point estimation method can be summarized as follows: Step 1: Define the number of selected uncertain inputs.
Step 2: Set the output variable vector of the jth moment E(V j )= 0.
Step 4: Set the uncertain parameter z r .
Step 5: Calculate the two standard locations of the variable z r : Step 6: Define the two locations of z r : Step 7: Run the deterministic EM based on the optimization algorithm for both locations z r .
Step 10: Repeat Steps 5-9 for r = r + 1 until the selected random inputs are taken.
Step 11: Compute the deterministic EM based on the following vector: Step 12: Determine the weight factor for the EM problem obtained in step 11: Step 13: Update E(V j ) Step 14: After obtaining the output statistical moments, the µ and σ can be computed: Depending on the values of µ and σ along with the method of Gram-Charlier as in [36], PDF and CDF are calculated for the output random variable.

The Equilibrium Optimizer Algorithm Overview
The EO algorithm is inspired by the analytical solution of a well-mixed dynamic equilibrium in a control volume. The balance equation of mass is utilized to describe the nonreactive constituent concentration on a control volume. It is known that the balance equation plays an important role in providing principles of physics to maintain the mass produced entering and leaving a control volume. Like other several recent optimization algorithms, at the beginning of the optimization procedure of EO, a random population is generated in the dimensional of search space for the optimization problem, where the particle acts as a solution and the concentration acts as the position of the particle [28].
According to the number of particles, the population and the dimensions of search space and the uniform random initialization of the initial concentrations are described as follows: At the beginning of the search process, the concentration degree along with the equilibrium state is unknown. For providing a search style of particles, only the equilibrium candidate is calculated. However, for estimating the equilibrium candidate solutions, the fitness function of all particles is evaluated, then the best scores are saved and sorted. The end convergence state is called the equilibrium state which is called the global optimal solution.
Based on various experiments of optimization problems, the EO algorithm allocates the four best particle solutions in the generated population at equilibrium candidate solutions during the implementation of the whole optimization procedure in addition to another particle having a concentration equal to the mean value of the previous proper four best solutions. However, these five candidate solutions are stored in a vector, which is called an equilibrium pool.
It worth noting that the four particle candidates make the EO algorithm capable of improving the exploration, while others support the exploitation process [27]. C e,pool = {C e,1 , C e,2 , C e,3 , C e,4 , C e,mean } For each iteration, the particle concentration is updated using a random style between the selected candidates that are selected with the same probability. However, in the first iteration, the first selected particle will update all concentrations based on C. Then, in the second iteration, the concentrations of the same selected particle are updated based on C e,mean . The updating process is repeated for each particle until reaching the process end. The exponential parameter term U contributes to the concentration updating rule for achieving the appropriate balance between exploitation and exploration, where it is computed as follows: The value of t is reduced while the iteration number increased.
However, the component sign(r 1 − 0.5) has an impact on the directions of both exploration and exploitation. The following formula should be implemented for improving the capability of both the exploration and exploitation of the EO algorithm to guarantee convergence. To improve the exploitation phase, the generation rate is developed as the most important step in the procedure of the EO algorithm for providing the optimal solution.
where the updating rule of the EO algorithm is as follows: The second term will be set to find an optimal position through globally searching the space upon its discovering a position. The third term will contribute to making the solution more accurate. The added memory is utilized to store the procedures and assist the individual particle to follow the path of its coordinates in the dimension of search space and notify it about its better score. The fittest value of every particle in the current iteration is compared to an identical value in the previous iteration and will be stored if it results in an improved value. This procedure assists the ability of exploitation but may raise the chance of falling in local minima if the technique does not achieve benefit from the ability of global exploration. More details of the EO algorithm are reported in [27]. Figure 2 shows the flowchart of the EO.

The Simulation Results
To prove the effectiveness of the EO algorithm for solving both deterministic and probabilistic EM problems with the consideration of the emission from RESs as well as ES devices in the MG system, the EO is investigated through the MG test system shown in Figure 1 device such as a NiMH-battery [30,31]. The MG is connected to the main grid via a power transformer.
The power is exchanged between the main grid and MG according to the decisions from the control center of MG for trading energy. Further, all MG sources only produce active power at unity power factor. Moreover, the demand loads inside MG during the day consist of one feeder for supplying light commercial loads, a second feeder for serving the residential region, and the third feeder for supplying a small industrial load. However, approximately the total demand load is equivalent to 1695 kWh per day. The study is carried out for the next day (24 h) at different loads. The data of the MG system are detailed and taken from [7,10,24]. Table 1 displays the minimum and maximum generation, the bid coefficients, and the emissions of MGs [8,30]. The forecasted demand load bid coefficients of the main grid and the expected output power of PV and WT over one day are shown in Figure 3 [10].

Case 1: The Optimization of EM without Emission
In this case, the EO algorithm is implemented based on the first scenario for solving the EM deterministic optimization problem to minimize the total operating cost without considering the pullout emissions. The obtained results of the EO algorithm are compared with those results obtained from reported algorithms in [30,31]. All MG sources operate at each hour during the day, which means that the operator should buy at least the minimum output power generated from these sources. The initial charging of the battery is infinitive.
The simulation results obtained using the EO algorithm for Scenario 1 are tabulated in Table 2. It is seen that because the PEM-FC unit has a lower bid compared to the other units, the PEM-FC operates at each hour of the day and delivers its maximum output power. In contrast, MT units which have a higher value of bid deliver their minimum value of output power for 13 h per day.
The results also show that, during the light load periods, a large part of the load is supplied by the main grid and PEM-FC because their bids are low. In addition, the battery In this study, the important parameters of the EO algorithm are chosen before its applications for solving the EM optimization problem. These parameters are computed using experimental tests. In addition, all simulation evaluations have been executed using MATLAB 2016b on a 2.9-GHz i7 PC with 16-GB of RAM. The EO is applied to compute EM in both deterministic and probabilistic manners according to three considered scenarios: • Scenario 1: It is supposed that all MG sources operate over the examined interval. The PV and WT are represented to deliver the forecasted maximum output power during each hour. In contrast, the main grid, PEMFC, MT, and battery operate according to their output power limits for achieving the constraints. • Scenario 2: All MG sources are operated at their output power limits while satisfying operational constraints. • Scenario 3: All MG sources except the main grid act as in Scenario 2, while the main grid is represented as an unconstrained source. Therefore, the energy will be exchanged without any limitations.
These scenarios are employed in the following cases studies.

Case 1: The Optimization of EM without Emission
In this case, the EO algorithm is implemented based on the first scenario for solving the EM deterministic optimization problem to minimize the total operating cost without considering the pullout emissions. The obtained results of the EO algorithm are compared with those results obtained from reported algorithms in [30,31]. All MG sources operate at each hour during the day, which means that the operator should buy at least the minimum output power generated from these sources. The initial charging of the battery is infinitive.
The simulation results obtained using the EO algorithm for Scenario 1 are tabulated in Table 2. It is seen that because the PEM-FC unit has a lower bid compared to the other units, the PEM-FC operates at each hour of the day and delivers its maximum output power. In contrast, MT units which have a higher value of bid deliver their minimum value of output power for 13 h per day. The results also show that, during the light load periods, a large part of the load is supplied by the main grid and PEM-FC because their bids are low. In addition, the battery was charged during this period. In contrast, in the periods of large load, the bid of the main grid is expensive, which required increasing the output power of MG. Additionally, the battery operates in discharge mode, and the excess energy from MG is exported to the main grid.
Further, 20 trial runs were performed for each algorithm for the purpose of comparing results and the simulation results as shown in Table 3. The results validate the effectiveness of the EO algorithm with the best value of the objective function when compared with the other techniques. Moreover, the convergence characteristic curve of the EO is shown in Figure 4, which demonstrates the robustness and fast convergence characteristic of the EO technique.

Case 2: The Optimization of EM with Emission Consideration
In this case, the deterministic EM problem is solved for minimizing the total objective function that includes minimizing emissions. The optimization problem is solved according to the mentioned previous three scenarios, and the results for each scenario are tabulated in Tables 4-6, respectively. The best results are achieved based on the EO algorithm; all operational constraints are also satisfied. In addition, Figure 5 shows the convergence characteristic of the EO algorithm for all scenarios in this case that validate the fast convergence characteristic of the EO method for all scenarios.  By inspecting the numerical results, during the period of light load, a large amount of load is supplied from the main grid and battery because these units offer both low bids and low emission as compared with other sources over these intervals. In the interval of heavy load, the output power from MTs and FCs is increased due to increasing bids through these intervals. Therefore, the main grid buys the energy from MG that assists in decreasing the cost of MG. From the results of the second scenario, the WT and PV are modeled for operating within their output power limits. The objective function is reduced by 2.74% when compared to the first scenario. Furthermore, when the main grid acts as an unconstrained source in the third scenario, the objective function is reduced by 9.46% and 6.91% as compared to the first and second scenarios, respectively. Due to the WT and PV having higher bids when compared with other MG sources, these sources have low contributions for supplying the demand load in the third scenario.
To prove the search capability and examine the superiority of the EO algorithm for solving the optimization problem in this case, the analysis is implemented for other recently published techniques to compare with the EO method. The compared methods are the artificial electric field algorithm (AEFA) [37], the bonobo optimizer (BO) [38], the modified moth flame optimization technique (MMFO) [39], and ESSA [8]. The results after 30 independent runs of EO are compared with AEFA, BO, MMFO, and ESSA. The statistical results are shown in Tables 7-9 for all considered scenarios, which reflect the satisfying performance of EO over other techniques that have the best value of the optimized objective function. It is seen that the EO has the best values of SD when compared with other techniques, which confirm the premium robustness of the EO technique.
Finally, due to the importance time for calculating this type of EM problem related to the computation of operating costs, Table 10 displays the elapsed time of implementing all compared algorithms. It is seen that the EO takes a shorter time as compared to other algorithms for all three suggested scenarios.

Case 3: The Probabilistic EM Problem
In this case, for obtaining the optimal solution of the EM problem in the probabilistic frame, the point estimation method (2m + 1) is evaluated along with the EO algorithm and other compared techniques. Moreover, the output power generated from PV and WT, market price, and demand load level are adopted as uncorrelated random input variables. Utilizing proper PDF modeling for hourly data of solar irradiance and wind speed as described in Sections 3.1.1 and 3.1.2, the output power of PV and WT units are computed. Additionally, the market price and demand load follow the same PDFs as given in Section 3.1.2. It is supposed that the output power of the PV and WT units, market price, and demand load have a normal distribution at an STD of 5% [10]. After that, the 2m + 1 method is employed, and the concentrations (locations) of the input random variables are calculated during the implementation of 2m + 1 for solving the probabilistic EM problem as displayed in Table 11. The EO algorithm and other compared techniques are applied for each case data results obtained from (2m + 1) as presented in Table 11. Finally, the statistical moments, mean values, and STD of the EM output random variables (the operating cost of MG) are determined.
Only the second scenario is considered in this case, where all MG sources are represented within their power limits. The predictable output power generated from MG sources and the corresponding cost based on EO algorithm are obtained using the mean values of the random input variables as displayed in Table 11 (the last four columns). The best simulation results of the EO algorithm are shown in Table 12, which shows that all operational constraints are met. However, the convergence curve of EO for this case is shown in Figure 6, in which the results indicate to the fast convergence characteristic of EO. It is seen that, during light load periods, the produced energy is supplied from the main grid and battery, causing these sources to have lower market prices and lower emissions when compared with other sources over these periods. Additionally, the contributions of PV and WT are lower than other sources because these units have a higher market price.
Using the Gram-Charlier expression for series expansion of the obtained results based on solving EM using the EO along with the two locations of each input random variables, the PDFs and CDFs are obtained for the optimal solutions of the probabilistic EM as shown in Figure 7. It is clear that there is no difference between the continuous random variables and the normal distributions as seen in Figure 7, which also shows that the EO algorithm gives the best results in the case of solving the EM in a probabilistic manner.  Using the Gram-Charlier expression for series expansion of the obtained results based on solving EM using the EO along with the two locations of each input random variables, the PDFs and CDFs are obtained for the optimal solutions of the probabilistic EM as shown in Figure 7. It is clear that there is no difference between the continuous random variables and the normal distributions as seen in Figure 7, which also shows that the EO algorithm gives the best results in the case of solving the EM in a probabilistic manner. The superiority of the EO is inspected here over the stated previous techniques in the previous subsection. Table 13 shows the statistical results of 30 independent runs of the The superiority of the EO is inspected here over the stated previous techniques in the previous subsection. Table 13 shows the statistical results of 30 independent runs of the EO algorithm compared with AFEA, BO, MMFO, and ESSA. The results of the proposed EO for solving probabilistic EM show that the EO has better effectiveness and performance in comparison with other techniques. Additionally, the EO has premium robustness due to having a lower value of SD compared with other techniques. Moreover, the superiority of the EO is inspected here over the stated previous techniques in the previous subsection.
By analyzing the results for Cases 2 and 3 based on Scenario 2, the total objective function is increased in Case 3 due to considering the impact of uncertainty when solving the optimization problem. This proves the impact of considering the uncertainty of input variables of the EM problem and the decision making of operators because the probabilistic analysis of EM gives more accurate and trustworthy results.

Conclusions
An efficient EO algorithm was suggested in this paper with successful implementation to solve the EM optimization problem for a typical MG system. Furthermore, the probabilistic framework using the 2m + 1 method was proposed to model the uncertainties of input random variables. The proposed EO was utilized for finding the optimal solution of EM incorporating ES devices, considering the emissions of MG sources based on converting the multi-objective function using the penalty factor price to the single objective function to minimize the total operating costs (fuel cost, startup and shutdown cost, as well as emission cost). The efficiency and advantages of the EO have been proven using a standard MG test system operating in grid-connected mode with different operating scenarios of the deterministic and probabilistic EM problem. Solving EM based on a probabilistic approach provided useful realistic decision making for operators of the MG system and assisted in finding the impact of the uncertainties for the input random variables on the statistical indicators which describe the state of the MG system. It was confirmed that the total operating cost for Scenario 2 based on probabilistic EM had higher value than deterministic EM by 6.53 EUR (0.346%) per day. The results confirm smooth and reliable EO convergence without any oscillation in the response. This confirmed the superiority of the EO over other compared optimization algorithms. In this context, the EO could be a useful decision-making tool for operators in the MG control center.

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

F 1 (R)
The objective function of operating cost F 2 (R) The objective function of pollutant emissions R Vector of the control variables T Total number of time periods in hours NG, NS Total number of DG units and energy storage devices, respectively P Gi,t , P Sk,t Active power produced from both DG unit i and energy storage device k at specified time t, respectively The skewness and kurtosis of input random variable z r , respectively µ z r , σ z r Mean and standard deviation of z r , respectively W r,l Weight factors of z r C k d The initial concentration vector of kth particles in dth dimension C min d , C max d The minimum and maximum values for search space dimension dth, respectively rand k d , ν, r 1 , r 2 , r 3 Random vector numbers between 0 and 1 m The number of population particles itr The current iteration number maxitr The maximum iterations number h 1 , h 2 Constant values that used to control both exploitation and exploration capabilities, were set to 1 and 2 respectively δ Denotes a decay constant GP Represents generation rate probability P WT The output power of WT P nom The nominal output power ν ci Cut-in wind speed of WT ν nom Nominal wind speed of WT ν co Cut-out wind speed of WT ν The wind speed

P PV
The output power of PV P STC The maximum power of PV module at standard test conditions (STC)

I S
The solar irradiance on the surface of PV module B Temperature coefficient of the PV module T c Temperature of the PV module T a The ambient air temperature T NOCT The nominal operating cell temperature (C) of the PV module K shape parameter of the Weibull distribution C Scale parameter of the Weibull distribution R uniformly distributed random numbers on [0, 1] ν m The mean wind speed µ , σ The Mean and standard deviation, respectively