A Framework for Profit Maximization in a Grid-Connected Microgrid With Hybrid Resources Using a Novel Rule Base-BAT Algorithm

In this paper, an energy management system (EMS) is proposed for optimal operation of a microgrid (MG). Dispersed photovoltaic arrays (PVs) and wind turbine generators (WTs) as renewable energy sources (RES) supply a major part of the network demanded energy. Also, an energy storage system (ESS), a micro-turbine unit (MT), and a fuel cell unit (FC) are integrated. The uncertainty and stochastic nature of the network load and RES data are treated via probabilistic modeling and scenario-selection approach. The predicted day-ahead data of the most diverse hourly scenarios are entered into the proposed EMS to determine the active and reactive power (P-Q) participations of local distributed resources. Likewise, it specifies the discharging/charging power and state of the ESS in addition to the exchanged active/reactive power amounts with the main network. The main goal is to maximize the profit of the secondary grid while satisfying all technical constraints. In the proposed EMS, the day-ahead energy management is developed as a comprehensive optimization problem. Moreover, the paper proposes novel modifications to improve the BAT optimization technique. The optimization problem of the energy management in the microgrid is implemented using a new integrated rule base–improved BAT method. Furthermore, the proposed EMS competence is proven by comparing its performance to recent literature.


PDF
Probability density function SOC State of charge MPPT Maximum power point tracking MG PRO,t The profit of MG at the t th hour TOC t Operation cost of the MG at the t th hour TMC t Maintenance cost of the MG at the t th hour s Index for scenario number N s Total number of selected scenarios ρ s Occurrence probability of the s th scenario.

SP t
The hourly magnitude of sold active power (kWh) λ P,buy,t The hourly buying price of active power of the market ($/kWh) SQ t The hourly magnitude of sold reactive power (kVArh) λ Q,buy,t The hourly buying price of reactive power of the market ($/kVArh) P DGr,i,y The rated power of the i th DG (kVA) of the y th DG type MC DG,i,y The i th DG maintenance cost ($/kVA) of the y th DG type n DG,y Number of DGs of the y th type Y The number of installed DG types n ESS Number of ESSs MC ESS,k , k th ESS maintenance cost ($/kVA) S ESSr,k Rated capacity of the k th ESS (kVA) MC OLTC The OLTC maintenance cost ($/kVAh) S OLTC,t The OLTC power flow (kVAh) at hour t C grid,t The hourly operating costs at time t for main grid C d dg,t The hourly operating costs at time t for the d th DG unit C k ESS,t The hourly operating costs at time t for the k th ESS unit n dg Number of dispatchable DG units n ESS Number of ESS units C Pgrid,t Hourly active power cost of the main grid at time t C Qgrid,t Hourly reactive power cost of the main grid at time t λ P,sell,t Price of active power sold by the MG to the main grid at time t λ Q,sell,t Price of reactive power sold by the MG to the main grid at time t P grid,t Active power of the main grid at time t Q grid,t Reactive power of the main grid at time t C dg,t The operating cost of a dispatchable DG unit at time t λ P,t Hourly generation costs of active power at time t λ Q,t Hourly generation costs of reactive power at time t C ESS,t The hourly operating cost of the ESS at time t P ESS,t The average hourly ESS power at time t η ch The charging efficiency η dis Discharging efficiency of the ESS P imp, max Maximum imported main grid active power from the MG P exp, max Maximum exported main grid active power to the MG P dg,t The active power of a DG unit at time t P dg, max Maximum active power limit of DG unit P dg, min Minimum active power limit of DG unit Q dg,t The output reactive power of a DG unit at time t Q dg, max Maximum reactive power limit of DG unit Q dg, min Minimum reactive power limit of DG unit P max,Ch Maximum charging active power of ESS P max, Dis Maximum discharging active power of ESS

SOC t
The hourly charge state of ESS at the t th hour SOC min The minimum SOC t of ESS SOC max The maximum SOC t of ESS V i i th bus voltage V min Allowable minimum voltage of the bus V max Allowable maximum voltage of the bus S ij The apparent power flow between node i and node j S max Maximum power flow between node i and node j T OLTC,t Tap position of OLTC at t th hour T min Minimum tap position of OLTC Velocity and position of bat w b,it Velocity weight factor at the it th iteration rand (·) Normally distributed random number r b,it Bat pulse rate at the it th iteration A b,it Loudness level of the b th bat at the it th iteration Value of the objective function F at the position x new b,it+1 α a , α r and α w Loudness weight factor, pulse rate weight factor, and velocity damping factor, respectively F new New value of the objective function F old Old value of the objective function

I. INTRODUCTION
Nowadays, microgrids (MGs) are acquiring more widespread as they can lower the electrical energy supply cost for consumers, improve ecosystem, and mend energy efficiency. MG is defined as a distribution network incorporating several kinds of distributed energy resources, mostly renewables, to provide powers of loads in a grid-tied mode or isolated mode [1]- [3]. To maximize the economic and technical benefits of MG, the integrated distributed generators (DGs) and energy storage systems (ESS) must have a coordinated control. Therefore, the MG needs an energy management system (EMS) [4]. EMS is the supervisory monitoring and control system of the MG that properly runs its generation and transmission facilities to accomplish its techno-economic objectives [5]. The adopted assessment algorithms and control strategies in the EMS define its aptness. The most common methods utilized in EMS of MG involve particle swarm optimization (PSO), artificial intelligence approaches, BAT algorithm, hyper-spherical search method and mathematical programming [6]- [15]. A central optimization solver collects all the system data and then processes them to obtain the solution and issue the MG operation decisions [3], [5]. However, papers [10]- [15] ignored the associated uncertainties in the forecasts of renewable energy sources (RES) output power and load power. Reference [6] and [7] have the same shortage. On the other hand, a PSObased probabilistic day-ahead EMS considering uncertainties in wind, solar and load power is presented in [8], [9]. Similarly, many other optimization methods are employed to design EMS for MG under uncertainty [16]- [26]. Reference [16] describes a second order cone program approach for the same problem. The authors in [17] present a non-dominated sorting genetic algorithm-based stochastic day-ahead EMS considering uncertainties in RES and load power. Cuckoo optimization algorithm and a grey wolf optimization algorithm are used in [18] and [19], respectively, for the EMS design problem. Reference [20] studies a nonlinear programming-based day-ahead EMS considering the forecasting error of photovoltaic generator (PV) output-power. A mixed-integer programming algorithm-based day-ahead EMS is presented in [21]. It accommodates uncertainties in PV and wind turbine (WT) output-power. However, some of works in [16]- [26] overlooked uncertainty of load power and the use of non-renewable dispatchable resources such as micro-turbine (MT) and fuel cell (FC). In other some, the developed EMS is only applied to a simple one-bus MG.
Different approaches are employed to handle the random-nature and uncertainty in input variables, such as the load and RES forecasts. In [8], authors used a number of probabilistic scenarios generated by the point estimate method. In [16], authors presented a set of stochastic scenarios which are generated from probabilistic forecasts. In [9], [17]- [19], [22]- [24], authors produced probabilistic scenarios using the Monte-Carlo method. Wind speed and solar radiation are modeled as Weibull and Beta probability distribution, respectively. The load is modeled as a random variable with normal probability distribution. Reference [21] behaves similarly but assumed normal probability distributions for both solar radiation and wind speed. Besides, uncertainty in PV power-output is modeled using the K-Means clustering method to produce a number of possible scenarios in [25]. Authors in [26] got a number of random scenarios using the Latin hypercube sampling method. Wind speed and solar radiation are described using Log-normal probability distributions. Reference [16] considers the correlation between load and other variables like temperature and wind speed. So, random input variables are treated as dependent. On the other hand, most of the published literature in MG analysis assumes the input random variables to be independent.
Moreover, all the aforementioned research works [3]- [26], except [6], consider only active power shares of the dispatchable DGs. Also, only the active power exchange among the MG DGs and the main grid is concerned. No attention is paid neither to the optimal reactive power dispatch among DGs nor to the exchange of reactive power between the main grid and the MG. So, this can amplify the probabilities of enlarged power loss, voltage deviations, and line congestion. Operation cost of the MG cannot be optimum and hence the MG profit gets compromised. As a distinct work, the EMS reported in [6] coped both active and reactive power (P-Q) shares among local DGs. Nonetheless, optimal utilization of commonly installed equipment like capacitor banks and on-load tap changer of the transformer (OLTC) connecting the MG to the main grid is not accommodated into the EMS. This can limit the MG profit on one side. It can endanger its security on the other side by giving higher chances to the occurrence of line congestions and voltage violations.
In this paper, the optimal day-ahead operation of a hybridenergy grid-tided MG is catered. The MG comprises various types of DGs and battery ESS. The DGs include PV units and WT units dispersed over the MG area. These RESs are assumed to be non-dispatchable DGs where maximum available energy must be extracted from them whenever possible. Besides, the MG has FC and MT units as dispatchable DGs. The MG is tied to the nearby main grid through a limitedcapacity transformer with OLTC. Based on the next day forecasts of metrological data, MG load variation, and time-ofuse electricity prices, the MG operation should be optimized by maximizing its profit satisfying all technical constraints. This is accomplished by determining the P-Q shares of all DGs, the charging/discharging state and power setting of the ESS, the amount of power exchange with main grid, and OLTC tap position for each hour in the next day. A comprehensive formulation of the constrained probabilistic MG mixed-integer optimization problem is developed. The paper takes the uncertainty in next-day forecasts of solar radiation, wind speed, and the MG load into account. An appropriate probabilistic model is defined for each variable. A largeenough number of possible input variables combinations (scenarios) are generated randomly by MCS. Then, a proper scenario-reduction technique is used to limit the required computation time keeping an adequate level of diversity in the selected scenarios. The optimal management of both P-Q shares of all DGs and ESS is concerned. The proposed design of the EMS is proved to be more efficient and faster in comparison to other recent methods reported in literature.
The main contributions of this paper can be summarized as follows: Uncertainty in input data is considered. Novel modifications are made to the BAT technique to enhance its searching capability avoiding trapping in local or sub-optimal solutions and shorten the search time.
An adaptive method is presented to identify the modified BAT algorithm parameters.
A new rule base-assisted improved BAT method is proposed as the solution algorithm for the MG optimization problem.
Optimal locations and sizes of all DGs and ESS are determined as a preliminary stage. A method for post-processing the optimization results is proposed to further improve the MG performance by intuitive readjustment of OLTC. In the rest of this paper, the problem statement is presented in Section II. Modeling of random input variables and scenarios selection are discussed in Section III. The proposed EMS design is described in Section IV. The case study system is introduced in Section V. Performance analysis and comparative evaluation are provided in Section VI. Finally, main conclusions are given in Section VII.

II. PROBLEM STATEMENT
In this paper, the day-head P-Q setting of each dispatchable DG, main grid, and ESS are determined to maximize the MG profit. The profit of MG for the t th hour (MG PRO,t ) depends on the outputted P-Q of the main grid, ESS, and dispatchable DGs. It is noticed that the reactive power supplied by the main grid is an explicit function of the tap position adjustment of the OLTC. MG PRO,t is expressed as: Random variation of load and RES are considered in (1). Sufficient diverse scenarios are selected to reflect the stochastic nature of these input variables as discussed next in Section III. The first two items in (1) mean the total hourly revenue/income from the selling P-Q to the consumers' loads, respectively. The next two terms signify the operation and maintenance costs of the MG, respectively. s and N s are an index and total number of selected scenarios, respectively. ρ s is occurrence probability of the s th scenario. SP t and SQ t are the hourly magnitudes of P-Q of the MG load at time t. λ P,buy,t and λ Q,buy,t are the hourly P-Q purchasing prices from the market at time t.
The total hourly maintenance cost of DGs and ESS is expressed as follows: where, P DGr,i,y , MC DG,i,y , n DG,y are the rated power of the i th DG (kVA) of the y th DG type, the i th DG maintenance cost ($/kVA) of the y th DG type, and number of DGs of the y th type, respectively. Y is the number of installed DG types. n ESS , MC ESS,k , S ESSr,k are number of ESSs, the k th ESS maintenance cost ($/kVA), and the rated capacity of the k th ESS (kVA), respectively. MC OLTC is the OLTC maintenance cost ($/kVAh) and S OLTC,t is the OLTC power flow (kVAh) at hour t.
The MG operating cost at the t th hour (TOC t ) is estimated by (3).
where, C grid,t , C d dg,t C k ESS,t are the hourly operating costs at time t for main grid, the d th DG unit, and the k th ESS unit respectively. n dg and n ESS are the number of dispatchable DG units and ESS units, respectively.
C grid,t is computed as [9], [12]: where, C Pgrid,t and C Qgrid,t are the hourly costs of P-Q of the main grid at time t. λ P,sell,t , λ Q,sell,t are the prices of P-Q purchased by the main grid from the MG at time t. P grid,t , Q grid,t are P-Q of the main grid at time t. In this work, λ P,sell,t , λ Q,sell,t are taken as 90% of λ P,buy,t , λ Q,buy,t , respectively [12]. The operating cost of a dispatchable DG unit at time t (C dg,t ) is given as [6], [12], [27]: where, λ P,t and λ Q,t are the hourly generation costs of P-Q at time t. P dg,t and Q dg,t are the mean generated P-Q at time t, respectively. The hourly operating cost of the ESS at time t (C ESS,t ) is estimated as [27]: where P ESS,t is the average hourly ESS power at time t, η ch is the charging efficiency, and η dis is the discharging efficiency of the ESS.
To summarize and clarify more the structure of equation (1), the MG profit has two components. The first is the profit due to selling power to the MG loads. It is calculated as the difference between the selling revenue minus the MG DGs and ESS cost excluding the main grid revenue/cost. This profit component is slightly affected by selling price as the loads buy the power at the power buying rate of the market. The second profit component is profit due to selling power to the main grid. It is computed as the revenue of selling power to the main grid using equations (4)-(6).

The objective function of the tackled operation problem of MG is
Subject to the following set of technical constraints: (1) Main grid constraints P imp,max , P exp,max are maximum imported and exported main grid active power from/to the MG, respectively.
(2) DG constraints where, P dg,t is the output active power of a DG unit at time t. P dg,max , P dg,min are maximum and minimum active power limits of this DG unit, respectively. Q dg,t is the output reactive power of a DG unit at time t. Q dg,max , Q dg,min are maximum and minimum reactive power limits of this DG unit, respectively.
(3) ESS constraints [12] SOC min ≤ SOC t ≤ SOC max (13) O Charging mode (P ESS,t > 0) O Discharging mode (P ESS,t < 0) where P max,Ch and P max,Dis are maximum charging and discharging active power of ESS, respectively. SOC t is the hourly charge state of ESS at the t th hour. t is a time step. SOC min and SOC max are the minimum and maximum SOC t of ESS, respectively. (4) Bus voltage constraints where V i is the i th bus voltage, V min and V max are the allowable minimum and maximum voltage of the bus. (5) Power flow constraints of line and transformer where S ij is the apparent power flow between node i and node j, and S max is the maximum power flow between node i and node j.
where T OLTC,t is the tap position of OLTC at t th hour, T min and T max are the minimum and maximum tap position of OLTC. N mov is the total number of daily movements of OLTC, N mov,max is the maximum number of daily movement steps.

III. PROBABILISTIC LOAD AND RENEWABLE GENERATION MODELS
To represent the prediction errors of load and RES, scenariobased technique is used [28]. At first, Monte Carlo simulation (MCS) is employed to generate abundant scenarios that simulate the uncertainty of load and RES based on probability distribution functions (PDFs) [29]. Hourly forecasts of load, solar irradiance and wind speed over the scheduling time should be obtained and used as mean values [17]. Furthermore, the standard deviations of these variables are calculated from their historical data [30].

A. LOAD MODELING
Load is modeled as a normally-distributed random variable [29].

B. WIND SPEED MODELING
Wind speed is represented with Weibull distribution [24]. The PDF of wind speed (v w ) is given as [24]: where, c and k are scale factor and shape factor, respectively. The Weibull parameters c and k are determined as [31]: µ v and σ v are mean and standard deviation of wind speed, respectively.

C. SOLAR IRRADIANCE MODELING
The solar irradiance (S i ) is represented as Beta-distributed random variable [18]. Its PDF is expressed as [18]: where, α and β are parameters of the Beta PDF. They are estimated as follows: µ s and σ s are mean and standard deviation of solar irradiance, respectively. Wind speed, solar irradiance and load are assumed to be independent variables [32]. Though the simulation accuracy gets better by increasing the number of scenarios, the enormous number of scenarios is obstructing in the probabilistic optimization techniques [21]. Consequently, after generation of a large enough number of possible scenarios for random input variables, an appropriate scenario reduction method must be applied to minimize the number of analyzed scenarios [24]. The scenario reduction algorithm computes the distance between each pair of scenarios and builds the distance array. Then, it finds the minimum allowed distance value from the distance array and removes scenarios with the lowest distance. After that, it updates the distance array for the remaining scenarios. The reduction stops when the number of remaining scenarios comes to a predefined accepted value [28], [29]. In this paper, a scenario reduction technique using backward reduction algorithm is used [28].

D. GENERATED POWER BY WT AND PV UNITS
The wind speed and solar irradiance values from the selected scenarios after the reduction algorithm are converted into the corresponding output powers. This is done by using the WT output characteristic curve and PV module characteristics, respectively as provided in [33], [34].

IV. PROPOSED DESIGN OF THE EMS A. PROPOSED MODIFICATION OF BAT TECHNIQUE
BAT is a recent optimization method basically developed in 2010 [35]. It is used for solving the different optimization problems [12]. The conventional BAT technique is presented with enough details in [35]. The BAT algorithm may converge to local solutions due to the lack of diversity of the bats in the search space. Also, conventional BAT technique suffers from slow convergence [12]. So, the conventional BAT algorithm is modified to overcome these drawbacks. Reference [12] changed only the velocity updating manner of the conventional BAT technique. To adjust the bat position, the authors in [36] advised a position updating way that depends on a sigmoid transfer function. This scheme has a demerit as the bats positions will remain unchanged when their velocity values continues to increase [37]. The authors analyze the effect of loudness and pulse rate values on the convergence speed and accuracy of the BAT algorithm in [37]. It is found that the convergence speed and precision of BAT technique improve by increasing the pulse rate and decreasing loudness of bats. Accordingly, a modification of BAT algorithm is made in [38], [39]. It depends on the adaptation of loudness or position equation of bats during the optimization process to accelerate the convergence speed. But no mechanism is given to adapt pulse rate and velocity of bats during the optimization process. Reference [40] replaced the loudness and pulse rate equations of BAT algorithm with the chaotic maps equations to heighten the convergence rate. Nevertheless, this work did not upgrade velocity update of bats during the optimization process. Reference [41] reformed the frequency equation only to improve the BAT algorithm performance. However, this may only increase velocity of bats during the optimization process. The obtained solution accuracy remains unchanged. To overcome the aforementioned limitations, this paper proposes a new modified BAT algorithm as explained below.
1. Initialize the bat population, each bat b generates an initial position vector x b,0 (initial vector of control variables to be determined) and a velocity (v b,0 ) randomly. 2. Calculate the objective function using (1) for each bat, and assign it as the best value of solution of bat b (F b,best ). Next, identify the global best value of solution (G best ) for all the population of bats. 3. Update the velocity (v b ) and the position (x b ) of each bat as proposed below: where, f b,min and f b,max are arbitrarily set to 0 and 2, respectively. w b,it is velocity weight factor at the it th iteration. rand(·) is a normally distributed random number.
Herein, the weight factor w b,it is added to the classic bat velocity equation in [35]. This weight factor, if intuitively selected, provides the balance between exploration and exploitation process during optimization. Thus, it controls the convergence speed to the global optimal solution avoiding the local optimal solution. 4. Generate an updated position (solution vector) for every bat locally by random walk based on pulse rate at the it th iteration (r b,it ) as follows.
x new b,it+1 = Gbest it + for rand(·) > r b,it (31) where, the factor is a normally distributed random number to limit the step sizes of random walks. 5. Check the continuity condition that is: If (rand(·) < A b,it and F x new b,it+1 ≥ F b,best ). Where, A b,it is the loudness level of the b th bat at the it th iteration and is the value of the objective function F at the position x new b,it+1 . If this continuity condition is satisfied, adjust the present position and optimal local solution of bat b as follows: Then, decrease A b,it , increase r b,it and decrease w b,it as follows: where, α a , α r and α w are loudness weight factor, pulse rate weight factor, and velocity damping factor, respectively. Here, a new update equation of the bat velocity weight factor w b,it is introduced. Also, the update equation of pulse rate r b,it is modified in this research VOLUME 8, 2020 to change linearly. An exponential change is assumed in the previous studies. These proposed modifications improve the BAT algorithm performance as demonstrated in later sections. If the above continuity condition is not satisfied, go to step 6. 6. Obtain the updated G best and the winner bat by comparing the updated F b,best of all bats. 7. Update the loudness and pulse rate of other bats to be the same as these of the winner bat b that got the latest best solution. This is a new step proposed in this work. It provides the rest of bats with enhanced parameters to reach a better global solution. 8. Increase the number of iterations by 1 and repeat steps 2 to 7. If the number of iterations is more than the preset maximum number of iterations, then stop the algorithm.

B. SOLUTION ALGORITHM OF THE MG OPERATION PROBLEM
The ESS will charge mainly from the main grid during the off-peak load periods at the low buying prices of the market. ESS will also be charged from the surplus renewable energy production from WT and PV units, if it is not fully charged. The ESS will be discharged during the on-peak load periods (i.e., high energy price of market). One charge/discharge cycle in the same day is permitted for the ESS. This means that the ESS is not permitted to be charged again in the same day if it had a discharging period following an earlier charging period in the same day. If this constraint is not considered, the ESS units economic feasibility becomes unsure [23]. A proposed combination of the new modified BAT method described above and rule base system is deployed to solve the operation problem of MG presented in Section II. The control vector of the investigated problem comprises the hourly values of P-Q of the dispatchable DGs (FC and MT), the ESS discharging/charging power, the operation mode of the ESS, the main grid power, and OLTC tap position. The solution algorithm is outlined in the following steps: 1. For the day hour number t, initial data of the problem are defined. The hourly probabilistic scenarios of the random input variables, solar radiation, wind speed, and the load data, are selected as in Section III. The technical and economic data for each DG unit, hourly P-Q selling/buying prices of the market, the BAT technique parameters, and the upper and lower limits of each variable are specified. 2. Generate a population of bats with a number of NPOP.
For each bat b, choose a random initial position vector x b,0 (control vector of the problem) and velocity v b,0 . Set x b,0 as the best position of bat b (F b,best ). 3. Every bat determines the discharge/charge state of the ESS based on the MG operator experience using the given below rule base: i. If the hourly active power buying price from the market (λ P,buy,t ) is below the price of ESS charging (λ ESS,Ch ) (i.e., 40% of the utmost value of active power buying price [9], [42]), then the bat defines ESS mode as a charging state and generates a positive random number of P ESS within the ESS charging limits, P ESS ∈ {P max,Ch , 0}, for the t th hour of the day. ii. If λ P,buy,t is more than λ ESS,Ch , and SOC of ESS at time t is below SOC max , and the expected sum of generated powers from all RES is above the expected sum of powers of loads (P load ), then the bat defines the ESS mode as a charging state. Then, it adjusts the power value of ESS P ESS as the difference power between the sum of generated powers from all RES and the sum of powers of loads at the t th hour of the day. If the calculated P ESS exceeds the ESS charging power limit P max,Ch , it is set at P max,Ch . iii. If the hourly active power selling price to the market (λ P,sell,t ) is higher than the price of ESS discharging (λ ESS,Dis ) (i.e., 70% of the utmost value of active power selling price [9], [42]), and SOC t is more than SOC min , then the bat defines ESS mode as a discharging state. Then, it generates a negative random number for P ESS within the ESS discharging power limits, P ESS ∈ {0, P max,Dis }, at the hour t of the day. 4. For each selected scenario of the input random variables, conduct a power flow analysis for the MG system using Newton-Raphson method. Then, evaluate the objective function in (1). 5. If the position vector x b of a bat fulfills all the constraints (10)-(21) for every selected scenario, go to step 6. Otherwise, this bat is infeasible and is excluded by imposing a high penalty [35]. 6. Compare a bat's objective function value with its individual F b,best . If the objective function value is higher than the current F b,best , set this value as the new F b,best . Then, identify G best . 7. Update the position and velocity of each bat as in (28)- (30). Then, generate a new position for each bat locally using random walk as in (31)-(32). 8. If the continuity condition is satisfied, modify the bat position and best individual solution as in (33)- (34). Also, update the loudness level, pulse rate, and velocity weight factor as in (35)- (37). Otherwise, move to step 9. 9. Obtain the updated G best and the winner bat by comparing the updated individual local best solutions of all bats. 10. Update the loudness and pulse rate of other bats to be the same as these of the winner bat that got the latest best solution. 11. If the highest number of iterations is reached, then go to step 12. Otherwise, increase the iterations counter by 1 and move to step 4. 12. Identify G best and the corresponding position vector x b .
These give the optimal solution of the problem for the day hour number t. 13. The optimization algorithm changes the hourly tap position of OLTC, keeping the total daily movements limit, to minimize the MG power loss and thereby increase MG profit. Thus, as the proposed algorithm works on hourly basis, it may consume the allowable daily OLTC movements in the early hours of the day and freeze the tap position for the rest of the day hours. This can influence the overall daily techno-economic performance of the MG. This is because the attained increase in the MG profit in some hour due to an avoidable tap movement can compromise the MG profit in the upcoming hours when the total tap movements is over. Accordingly, the hour-by-hour OLTC tap position profile obtained in the preceding step is re-checked to assure its necessity as follows: a. If the tap position of OLTC at this hour number t (T t ) differs from its position at the preceding hour number t − 1, (T t−1 ), Increase/decrease the current OLTC tap position by 1 to revert to T t−1 (if possible). Else, move to step 14. b. Run power flow for each selected scenario, and check all the MG constraints (10)- (21). If the constraints are satisfied, then save this new tap position (T t,n ). Otherwise, move to step 14. c. If T t,n equals T t−1 , then T t,n = T t−1 and go to step 14. Else, increase/decrease T t,n by 1, and repeat step 13.b above. The rule in 13.a checks if it is possible to make no change in tap position between hour t-1 and hour t. If this causes no violation of any operation constraint according to step 13.b, the procedure is stopped. In contrary, if this assumption is not permitted as it violates operation constraints, one tries to lessen the change in tap movements between hour t-1 and hour t as in step 13.c. 14. Solve the optimization problem for the following hour t+1 of the day and move to step 1 Fig.1 illustrates the solution procedures using the proposed algorithm.

C. TUNING OF THE MODIFIED BAT ALGORITHM PARAMETERS
The values of the weight factors of the BAT algorithm parameters should fit the nature of the problem under study. In this work, a scheme is proposed to determine proper values of the weight factors to make the BAT algorithm-based solution of the studied MG operation problem more accurate and faster. This scheme is as follows: 1. For an arbitrarily chosen day hour, run the solution algorithm explained in Section IV. B for the steps from 1 to 12 using the initial BAT algorithm parameters listed in Table 1.

If the obtained new value of the objective function F new
is greater than or equal to its old value F old , update the BAT algorithm weight factors as given below. Otherwise, go to step 4.
3. Go to step 1. 4. Save the values of the BAT algorithm weight factors. The adaptation procedure of the proposed BAT technique is illustrated in Fig.2. The finally determined BAT parameters are provided in Table 1. The values of the obtained weight  factors are noted to be unchanged when other day hours are taken.

V. SIMULATION SETUP AND DGS SIZING
The MG model used to evaluate the proposed EMS is a modified IEEE 33-bus distribution network shown in Fig.3 [43]. The total system load is 5.83 MW and 2.83 MVAr. The solution algorithm described in Section IV. B is coded in MATLAB and executed by a personal computer with 2.3GHz processor and 4GB RAM. The maximum exporting and importing power of the main grid are 5 MW and -3 MW, respectively.
The minimum and maximum tap position of OLTC are ±12. Table 2 provides operating costs for FC and MT [6]. Table 3 shows the ESS specifications [6], [8]. Technical maintenance costs per year for MT, FC, ESS, PV, and WT are 25,25,20,20, and 75 $/kVA, respectively [6], [12]. The hourly mean and standard deviation of solar irradiance, wind speed and load are depicted in Table 4 for a day in summer season as reported in [30]. The MCS technique is employed  to generate 2000 random scenarios for each day hour as discussed in Section III. The backward reduction algorithm is used to reduce the number of scenarios to 30 by selecting the most distinct ones [17], [28]. Fig.4 shows the probability distribution of solar irradiance, wind speed and load power at the time 08:00 A.M for this day. The daily variations of the MG feeders' loads are shown in Fig.5. Feeder1 includes loads on buses 1-6. Feeder2 comprises loads on buses 19-22. Feeder3 includes loads on buses 23-25. Feeder4 includes loads on buses 26-33. Feeder5 includes loads on buses 7-18. Fig.6 provides the forecasted mean hourly total load of MG, total active power output of PV units, and the total active power output of WT units for the same day. Fig.7 provides the forecasted hourly P-Q purchasing prices of market for the same day [9], [32]. The ESS initial state of charge is assumed as 25%.
To Determine the appropriate sizes and locations of the ESS and the DGs in the MG, the method described in [44] is applied. It is based on optimization-problem formulation to minimize active power losses and improving system volt-VOLUME 8, 2020  age level in the MG. Annual load growth of 3%, and 25years life time of the project are assumed. The mean metrological and load variation data given above in this section are used. This DGs sizing problem is separately coded in MATLAB and solved using PSO algorithm. The selected maximum iterations, population size, social parameter and cognitive parameter of the PSO algorithm are 100, 50, 2 and 2, respectively.
The obtained optimal ratings of renewable DGs, FC, MT, and ESS are specified as indicated in Table 5. The maximum ESS capacity got by the above optimization problem is 20MWh. However, this value may not produce the maximum MG profit as targeted in the operation problem of MG formed in Section II. So, a minimum ESS capacity is selected as 4MWh. Then, the optimal ESS size is identified by direct search between 4 and 20MWh in 1.5MWh step such that the profit of MG is maximized. An optimal ESS size of 16.5MWh  is determined. Therefore, four PV units are installed at buses 12, 19, 23, and 29. Four WT units are installed at buses 17, 21, 25, and 32. The ESS, FC, and MT units are installed at buses 5, 9, and 26, respectively. Specifications of PV module and WT are obtained from [31]. All PV and WT units operate at unity power factor and maximum power point tracking (MPPT) control mode. Meanwhile, FC, MT, and ESS operate at PQ control mode [12]. The maximum capacity for each DG unit is shown in Table 5.

A. SIMULATION RESULTS OF THE EMS
The solution algorithm described in Section IV. B, C is applied. Fig.8 depicts the determined hourly mean P-Q shares of all dispatchable DGs in the MG. It is observed that ESS will charge from the main grid and begin charging state from nearly 02:00 to the early morning at 07:00 a.m., during low buying price of energy and off-peak load time, as depicted in Fig.8. The daily MG profit is found as 5594 $. However,  the daily profit of MG will be reduced to 4898 $ if the energy storage system is not installed. Keeping in mind that the annual cost of the ESS is 175604 $/year [8], [12], an extra net annual profit of 78436 $/year is obtained due to the ESS installation in the MG. Thus, the pay-back period of the ESS is about two years. This is economically very encouraging as the expected ESS life span is more than 20 years according to the applied ESS charging/discharging strategy. Fig.9 provides the hourly change of the OLTC tap position using the adaptive operation imposed by step 13 in the solution algorithm in Section IV. B. The daily number of OLTC tap movement steps is maintained below the nine-movement steps limit as dictated by constraint (21) even this adaptive operation is halted. It is worthy to mention that the referred adaptive operation of OLTC improves the performance of the proposed algorithm and enlarges the maximum MG profit. This is because, without this adaptation procedure, the OLTC will consume all the daily allowable tap movements during the first hours of the day as shown in Fig.9. Therefore, after t = 10:00 a.m., the OLTC cannot make any further tap movements. So, during the discharging period of ESS from   [6], PSO-based method in [8], and BAT-based method in [12]; (a) mean value (p.u), (b) voltage standard deviation (p.u). 04:00 p.m. to 10:00 p.m., some buses voltages rise above the 1.05 p.u limit for some of the selected input variables scenarios. This forces the MT and FC power shares to the levels that prevent bus voltage violations not the levels that can further maximize the MG profit. Moreover, the rulebase system explained in step 3 of the solution algorithm in Section IV. B for controlling the ESS charge/discharge process obviously enhances the performance of the proposed EMS. If this rule-base is excluded, the daily profit of MG will be reduced to 5562 $.
Also, the maximum stored energy level in the ESS would be curtailed to 76.5% before beginning the state of discharge degrading its economic utilization. Besides, the computation time of the proposed solution algorithm on a commercial PC will increase to 2.67 hours instead of 2.31 hours. Thus, this rule-base system assists the proposed algorithm to converge faster and globally maximize the objective function over the  FIGURE 11. SOC of ESS using the proposed EMS, PSO-based method in [6], PSO-based method in [8], and BAT-based method in [12]. search space increasing the MG profit. Table 6 summarizes the proposed EMS performance. Table 7 holds a comparison to evaluate the performance of the proposed MG EMS versus other recent optimization-based approaches analyzed in [6], [8] and [12]. It should be noted that the proposed EMS achieves the best economic performance. From the technical aspect, it maintains the daily voltage profiles of all the MG buses within the allowable limits (between 0.95 and 1.05p.u) for all scenarios. Whereas, voltage violations arise by the methods in [8] and [12] as shown in Fig.10. Moreover, the proposed EMS maintains the lines power flow under the maximum permissible capacity. Convergence pattern at the 3rd hour by the proposed EMS, PSO-based method in [6], PSO-based method in [8], BAT-based method in [12], and the proposed EMS with the BAT modifications in [40].

B. COMPARISON TO OTHER METHODS
In addition, Fig.11 shows the evolution of ESS state of charge using the proposed EMS and recent methods presented in [6], [8] and [12]. It is observed that the input variables uncertainty clearly affects the operation of the ESS, in particular under the use of methods reported in [8] and [12]. The different scenarios considered to reflect uncertainties in input variables cause voltage violations problems that obstruct the optimal utilization of the ESS and decrease the profit of the MG. So, it is clear that the proposed EMS is more efficient and robust to manage MG under different situations and scenarios of uncertainties. Fig.12 portrays the convergence pattern of different algorithms during the solution of the energy management problem to maximize the MG profit at the 3 rd hour of the day. Evidently, the proposed EMS has the best convergence rate and accuracy. Besides, the convergence behavior of the proposed EMS using the chaotic maps' equations in [40] instead of equation (35)-(36) is also traced in Fig.12.
It is observed that this replacement would lower the speed and accuracy of convergence of the proposed BAT algorithm. Furthermore, the replacement will attenuate the daily MG profit to 5580 $. So, it can deduce that the proposed EMS is more effective, more swift, and more accurate than the approaches presented in [6], [8], [12], and [40].

VII. CONCLUSION
The paper proposes an EMS for optimal operation of a gridconnected MG. Uncertainties in RES and load power forecasts are analyzed by appropriate probabilistic models. The MG includes hybrid renewable DGs, other dispatchable DGs, and battery bank ESS. It is desired to maximize the profit of the MG meeting all technical limitations. The constraints include equipment capacity limits, nodes voltages, and line power flows. The problem is mathematically formulated as a stochastic constrained nonlinear optimization. Thousands of scenarios representing random variations of metrological and load data are generated. Then, the most diverse scenarios are assessed using a scenario reduction technique. New modifications are done to improve the BAT optimization solver. Moreover, a method for tuning the BAT algorithm weight factors is presented. Based on the MG operator experience, a rule-base is formed to regulate the working of the ESS. Then, a combination of the rule-base and the newly improved BAT technique are adopted to solve the optimization problem. A post-processing mechanism of the optimization outcomes is carried out to further refine the results by tuning the OLTC tap position. Thus, the optimal day-ahead P-Q dispatch of each dispatchable distributed generator, the discharging/charging power and state of the ESS, and the exchanged active/reactive power amount with the main grid are decided. The proposed EMS design maximizes the MG profit in a shorter computation time without breaking any constraint for all possible scenarios. The apparently superior performance of the proposed EMS is proven by comparison to other recent approaches. It is always found to be faster and more profitable to the MG.