A Metaheuristic Harris Hawk Optimization Approach for Coordinated Control of Energy Management in Distributed Generation Based Microgrids

Cost management of microgrids represents a real challenge since the power generation of microgrids is usually composed of different renewable and non-renewable sources. Additionally, it is always desired to make a connection between the microgrid and national grid to secure the load demand and to fit the regulations of liberated energy markets. Because of all these reasons, it is essential to develop a smart energy management unit to control different energy resources within the microgrid to achieve minimum operation costs. This paper presents a proposal for a smart unit for the cost management and operation of multi-source based microgrids. The proposed unit utilizes the Harris hawk optimization (HHO) algorithm which is used to optimize the cost of operation based on current load demand, energy prices and generation capacities. The proposed unit is tested on a microgrid with different energy resources using MATLAB while applying different operation scenarios. All simulation results show that the proposed unit succeeds in operating the microgrid at minimum cost. Obtained results are compared with other optimization algorithms and the proposed Harris hawk algorithm gives superior performance.


Introduction
The demand for renewable energy as a resource for electrical power is being increased rapidly. This is because of several reasons, including that they don't have negative impacts on environment and have low operation costs when compared to classical power generation methods. However, there are a lot of limitations and challenges that must be surmounted to harness the energy from renewable resources effectively [1]. One of these challenges of renewable energy resources is that they must be installed in limited locations. Moreover, their output power depends on the current environmental and weather conditions, which means that they may not be able to supply the demand at certain times. All these limitations have directed researchers to a new concept of electric distribution, which is the microgrid (MG) [2,3]. A microgrid is a distribution network that supplies certain loads from different distributed energy resources (DERs), such as photovoltaic (PV), wind turbines (WT), diesel engine generators (DEGs), fuel cells (FC) and battery storage [4]. Additionally, a microgrid usually has a connection with the unified grid. This connection is established to enable the power exchange between the grids to maximize the load security. Furthermore, the microgrid can sell excess power to unified grids [5].
To get the maximum benefits of microgrids, it is important to have an energy management system (EMS) to be responsible for the operation and control of different DERs and power exchange between the microgrid [6,7] and the grid as shown in Figure 1.
To get the maximum benefits of microgrids, it is important to have an energy management system (EMS) to be responsible for the operation and control of different DERs and power exchange between the microgrid [6,7] and the grid as shown in Figure 1. The main role of the EMS is to decide the amount of energy that needs to be supplied by the DERs within the microgrid and to organize the power exchange between the unified grid and the microgrids. The decisions of the EMS are based on the ability of DERs and the current prices of energy. This means that the EMS must take accurate decisions to balance between increasing local production and maximizing economic revenue. Consequently, energy management of microgrids forms a highly constrained optimization problem and this optimization problem is usually considered an offline one [8].
From the above-mentioned explanation, it is clear that the EMS's objective function depends on providing various types of data. Some of these data can be easily determined, such as the generation capacities of the grid and microturbines. However, most of them depend on forecasts and estimations such as loading limits, wind speed, solar irradiance, energy prices, etc. As a result, the researchers have followed two techniques in optimizing the EMS's objective function [9,10].
The first is the deterministic approach: in this approach, all forecasted data are taken into consideration as they are assuming that they have very high accuracy. The other approach is the probabilistic approach in which all forecasted data are put in the form of variables. Each variable has a probability that resembles the accuracy of the forecast.
Researchers have applied many optimization algorithms to solve the EMS problem. One of the famous techniques is linear programming (LP) which considers the power The main role of the EMS is to decide the amount of energy that needs to be supplied by the DERs within the microgrid and to organize the power exchange between the unified grid and the microgrids. The decisions of the EMS are based on the ability of DERs and the current prices of energy. This means that the EMS must take accurate decisions to balance between increasing local production and maximizing economic revenue. Consequently, energy management of microgrids forms a highly constrained optimization problem and this optimization problem is usually considered an offline one [8].
From the above-mentioned explanation, it is clear that the EMS's objective function depends on providing various types of data. Some of these data can be easily determined, such as the generation capacities of the grid and microturbines. However, most of them depend on forecasts and estimations such as loading limits, wind speed, solar irradiance, energy prices, etc. As a result, the researchers have followed two techniques in optimizing the EMS's objective function [9,10].
The first is the deterministic approach: in this approach, all forecasted data are taken into consideration as they are assuming that they have very high accuracy. The other approach is the probabilistic approach in which all forecasted data are put in the form of variables. Each variable has a probability that resembles the accuracy of the forecast.
Researchers have applied many optimization algorithms to solve the EMS problem. One of the famous techniques is linear programming (LP) which considers the power balance and generation limits of the distributed generation units. However, this method suffers from one main disadvantage, which is a high computational burden [11][12][13].
One of the popular optimization techniques that have been incorporated in solving the EMS problem is dynamic programming (DP). The main advantage of this method is its ability to divide the EMS into smaller sub-problems. This means that the optimization technique will solve sequential problems instead of solving one sophisticated problem, which helps in reaching the optimum solution with quick and accurate performance. However, the operation of this algorithm is considered difficult since it includes a high number of recursive functions [14].
In [15], a genetic algorithm is used in optimizing the EMS's objective function. This method has better computational burden when compared to LP. However, still the computational complexity exists.
Multi agent (MA) optimization has gained great researchers attraction in solving the EMS optimization problem. Although this method has relatively high accuracy considering many constrains, it suffers also from high computational complexity and time [16][17][18][19].
Particle swarm optimization (PSO) has attracted many researchers as well to use it in the optimization of the EMS problem. This method showed better results compared to previously mentioned methods in terms of accuracy and computational burden [15,20].
In [21], authors proposed the artificial bee colony (ABC) to solve the EMS problem. Although this method is considered quite simple and has robust performance, the formulation of the EMS problem to adapt between it and the algorithm is considered complex as this technique is suitable for objective functions that do not have high number of parameters.
This technique has directed many researchers to use various meta-heuristic optimization algorithms such as grey wolf optimization (GWO) [22], evolutionary algorithms (EA), etc.
Another research direction led by many researchers depends on coordinated control. In [23] a coordinated controller for a microgrid based on a predictive model and voltage control is presented. The main advantage of this method is that it needs less effort in the tuning of the controller. Authors of [24] also proposed a multi-layer coordinated control technique for the energy management of microgrids based on forecasting customer loading. Furthermore, another method of coordinated control is presented in [25]. This method depends on bi-level stochastic modelling.
Most of papers mentioned in literature have suffered from drawbacks (particularly in the modelling of the system) that can be summarized in the following points [11,12,15,21,22,[26][27][28][29]: • Power losses are not taken into considered in some papers. • Battery charging and discharging intervals are not considered in some papers.

•
High computational complexity in some of the above indicated references.
In this paper, the Harris hawk optimization (HHO) is used in solving the EMS optimization problem in a grid connected microgrids. The algorithm is applied to a highly constrained objective function which takes into consideration many important constraints such as power balance, generation power capacities, spinning reserve and energy storage limits (including charging and discharging rates). The proposed algorithm shows an improved performance compared to other algorithms used previously in literature in terms of speed of convergence and minimization of cost of operation when compared head-to-head and applied to same test bench. The proposed HHO algorithm may form a good optimization technique for solving other optimization problems in the power and energy fields.
The rest of this paper is organized as follows: Section 2 shows the problem formulation of the EMS optimization problem, while Section 3 explains the HHO algorithm and its application to the proposed EMS optimization algorithm. In Section 4 the results of applying the proposed EMS optimization algorithm to a real microgrid with two different operating scenarios are presented. In Sections 5 and 6, respectively, the discussion and conclusion are presented for the evaluation of the effectiveness of the proposed EMS algorithm.

EMS Problem Formulation
The EMS problem is defined as minimizing the operational cost of the microgrid by means of adjusting the generation parameters of different DERs and the utility grid. The optimization problem also considers open market prices for energy received from the grid and DERs.

Objective Function
As mentioned before, the objective function of the EMS is to minimize the operation cost of the microgrid as shown in Equation (1): where P t Gi is the generated active power of the generator i at time t. P t Grid is the power exchanged between the utility grid and the microgrid at time t. B Gi is the cost per kW of the active power generated by the ith DER. MP t is the market price of active power exchanged between the utility grid and the local microgrid.
NT is the total hours included in the objective function which is usually a 24 h period and NG is the total number of distributed generation units including storage units.

Constraints of the Objective Function
To simulate the real response of a microgrid against the decisions of the EMS, it is essential to put many constraints that are mandatory for smooth operation of the microgrid.
The first constraint is that all active power generated by all DERs and the grid must equal the load active power as per Equation (2).
where P t L D is the value of the load active power at time t, ND is total number of load levels. Additionally, it is important to include the generation limits for each DER and the utility grid as illustrated in Equations (3) and (4). These are mandatory to guarantee stable operation of the microgrid.
where P t Gi,min and P t Gi,max are the generating limits of the generator i at time t while P t Grid,min and P t Grid,max are the generating limits of the utility grid at time t. To keep the system running with high reliability, it is essential to include a portion of the microgrid generation capacity as spinning reserve, especially if most of the DERs are from renewable resources that may suffer from output instability. Equation (5) describes how the criteria of spinning reserves are met.
where R t is the scheduled spinning reserve at time t.
One of the important constraints for energy storage elements is the limitation of charge and discharge cycles, especially the time interval needed for each cycle. The amount of energy stored in the battery is calculated based on the following rule: W ess,min ≤ W ess,t ≤ W ess,max P charge,t ≤ P charge,max ; P discharge,t ≤ P discharge,max where W ess is the amount of energy stored in the battery. P charge is the allowed rate of charge while P discharge is the allowed rate of discharge. ∆t is a predefined time period in which the rate of charge or discharge is allowed.

Calculation of Energy Cost and Bids for the Microgrid Elements
According to the EMS objective function mentioned in Equation (1), the main role of the EMS is to reduce the energy cost of the microgrid. Because of this, it is very important to calculate the energy cost of each source used in the microgrid and the energy price of the utility grid.
Regarding the fuel cell and the microturbine, the bid price is calculated according to (8): where B G is the bid price in €/h and C f uel is the price of fuel in €/kW h. P G is the electrical output power in kW. η G is the efficiency of the generation unit. C inv is the payback amount of the cost of investment of the generation unit which can be get from (9) as follows: where P Gnom is the nominal generation power of the generation unit. AP is the annual energy production in kW h/kW. AC is the annual investment cost which is based on the following rule: where i is the interest rate, n is the depreciation period and IC is the installation cost.

Proposed HHO-Based EMS
Harris hawk optimization is a metaheuristic optimization which was proposed in [30], and its mathematical model is inspired from the cooperative behavior of Harris hawks in hunting, chasing and besieging their victims. The HHO is based on population optimization without having any gradients which gives it competitive edge over other techniques in terms of speed of conversion.
The HHO consists of two main phases: exploration and exploitation. Additionally, there is a transition phase though which the algorithm is switched from exploration to exploitation.
In the exploration phase, Harris hawks start to search randomly for victims as per the following equation: where X(t + 1) is the location of the hawks in the iteration (t + 1), X rabbit (t) is the location of the rabbit (the victim), r 1 to r 4 and q are random numbers that can vary between 0 and 1, X rand (t) represents a hawk which is chosen randomly and X m is the average location of the current population of hawks which can be calculated from Equation (12): where X i (t) indicates the position of each hawk at iteration t while N represents the total number of hawks.
As mentioned above, after finishing the exploration stage, there is a transient stage before moving to the exploitation stage. At this transient stage, it is necessary to model the energy of the rabbit as per Equation (13): where E is the escaping energy of the rabbit, T is the maximum number of iterations and E 0 is the initial state of the rabbit energy. The value of E 0 is varying between −1 and 1 based on the physical fitness of the victim. When E 0 goes towards −1, this means that the victim is losing its energy and vice versa.
According to the behavior of rabbits, the relation between the rabbit energy and the time is inversely proportional. This means that as long as t increases, the E is decreased. Based on E, Harris hawks also decide to either search different areas to detect the location of the rabbit when |E| ≥ 1 or move forward to the exploitation phase.
In the exploitation phase, there are two behaviors that need to be modelled. The first is the soft besiege in which the rabbit energy is still high and can run fast; in this condition Harris hawks try to softly follow and put it under surveillance until it starts to get exhausted. The second behavior is the hard besiege; the prey in this behavior is tired and has no sufficient energy to escape. As a result, the Harris hawks in this mode form closed circles to make a sudden attack. Figure 2 shows Harris hawk attack patterns. To mathematically model the two behaviors, let r be the percentage of the successful escape of the rabbit. If |E| ≥ 0.5 and r ≥ 0.5, this means that the rabbit has relatively high escaping energy and at the same time the chance of successful escape is higher than 50%. This means that the Harris hawks will perform a soft besiege and will update their location according to Equation (14): where ∆X(t) is the position difference between the rabbit and the hawks. This value can be calculated as follows: where J is a random number that represents the jump strength that can be met from Equation (16) as follows: where r 5 is a random number that varies between 0 and 1. If |E| ≥ 0.5 and r < 0.5, this means that the rabbit has high energy. However, the chance of successfully escaping is not big. In this case, the Harris hawks will perform a soft besiege but with progressive and rapid dives. The next movement of the hawks will be updated according to: The hawks then will compare the current position with the previous dive to evaluate which is better. If the previous dive is better, the hawks will use it. If not, the hawks will then apply new dive using the levy flight (LF) equation: where D is the problem dimension and S is a random vector with size of 1 × D. The LF function can be calculated according to (19) as noted in [31]: where u and v are random numbers that vary between 0 and 1. β is a constant value of 1.5. σ is calculated using: The Harris hawks will then evaluate the positions Y and Z and then select the best position based on (21): where only the better location of either Y or Z will be used to update the hawk's position. If |E| < 0.5 and r ≥ 0.5, this means that the rabbit has relatively low energy but it has a moderate chance of a successful escape. In this condition, the hawks will perform a hard besiege and will update their equation based on Equation (22): If |E| < 0.5 and r < 0.5, this means that the victim has low energy and also has a low chance to escape. In this situation, the hawks will also perform a hard besiege but with progressive rapid dives at which the next position of the hawks will be updated using Equation (21). Z will be calculated from Equation (18) and Y will be calculated using Equation (23) as follows: The proposed EMS uses the HHO in evaluating the optimum values for generated powers from different DGs within the microgrid. The HHO in this optimization problem will consider the DG-generated powers as the location of the Harris hawks and the prey will represent the energy cost. Figure 3 shows a flowchart of the proposed HHO algorithm while Figure 4 shows a flowchart for the proposed EMS algorithm. Start the random search process equation (11) Calculate the average location Xm as per equation (12) Calculate E as per equation (13) E > 0.5? r > 0.5? r > 0.5?

No
No Yes Yes X(t+1)=Y update X(t+1) as per equation (14) Evaluate next move mean , Y as per equation (17) Calculate the levy flight LF(D) as per equation (19) Find the required dives based on LF(D) , Z as per equation (18) Y < x(t)?
X(t+1)=Z update X(t+1) as per equation (22   Calculate the spinning reserve according to equation (4) Check that generated power equals load power as per equation (12) Apply energy storage limits of battery as per equation (6) Calculate the energy bid price using equation (8)- Optimize the objective function of equation (1) using HHO by applying equation (14)

Simulation Results
To validate the effectiveness of the proposed algorithm, a microgrid is simulated using MATLAB. The output power of its DG units and grid is set by the proposed EMS. The simulated microgrid consists of PV, WT, DEG, FC, battery storage and a connection to the utility grid. This benchmark has been used by previous research [31] and has been chosen to be simulated in this paper in order to be able to compare obtained results with previous research results available in literature. Table 1 shows the minimum and maximum output power of each element of the microgrid as well as coefficients of bid functions, where the negative sign indicates direction of power flow. Table 2

Simulation Results
To validate the effectiveness of the proposed algorithm, a microgrid is simulated using MATLAB. The output power of its DG units and grid is set by the proposed EMS. The simulated microgrid consists of PV, WT, DEG, FC, battery storage and a connection to the utility grid. This benchmark has been used by previous research [31] and has been chosen to be simulated in this paper in order to be able to compare obtained results with previous research results available in literature. Table 1 shows the minimum and maximum output power of each element of the microgrid as well as coefficients of bid functions, where the negative sign indicates direction of power flow. Table 2 presents sample forecast values for the load demand and DG units' powers along with market price. Two case studies are simulated to evaluate the effectiveness of the proposed HHO-based EMS on different operating conditions. Cost is expressed in units of euro cents (€ ct).  In this case, the PV and WT will supply their maximum available generating power according to the hourly forecast. The proposed EMS controls and sets the output power of the other DG units to decrease the generating cost and according to their corresponding power limits and constraints.
The results of the proposed algorithm are tabulated in Table 3. When observing the behavior of the proposed EMS during the first and last eight hours of the day, results show that during these periods the proposed EMS depends on the FC and the utility grid for generating the largest portion of output power. This occurs because the energy price during these periods is low which means that generating electricity from the utility grid is cheaper than the DEGs. This is not the case at midday, when the EMS mainly depends on the PV, WT and DEG units as the energy price of the utility grid is expensive.

Case Study 2
In this case, the output power of all DG units including PV and WT are controlled by the EMS within their constrains and power limits. The results of this case, which are shown in Table 4, show a huge reduction in the total cost of energy compared to that of case 1. The reason for this is that the EMS fully controls DG units aiming to reduce the total cost. This leads to higher dependency on the utility grid during low energy bids.
Results of both cases show the HHO algorithm functions properly during different modes of operation for the EMS. Figure 5, which indicates the conversion curve sample of the HHO algorithm, clearly shows the effectiveness of the proposed method as the conversion rate is fast during both cases. Results of both cases show the HHO algorithm functions properly during different modes of operation for the EMS. Figure 5, which indicates the conversion curve sample of the HHO algorithm, clearly shows the effectiveness of the proposed method as the conversion rate is fast during both cases.

Discussion and Analysis
In order to validate the obtained results confidently, it is important to compare the results of HHO algorithm in both cases with earlier research results that used the exact microgrid simulation model but with different optimization algorithms, such as fast surrogate-assisted particle swarm optimization (FSAPSO) and several modified PSO algorithms [31]. Table 5 shows the obtained results of the HHO algorithm along with results of other previously proposed algorithms in literature. Number of objective function evaluations (OFE) carried out till best solution is reached is indicated where it is equal to 1x number of agents x iterations. Additionally, the number of algorithm iterations needed till best solution is reached is indicated in the table along the trials number. Although there is a lack in information in [31], it can be noted that the proposed HHO algorithm successfully reaches lower operating costs with lower iterations. It is suggested as future work to take into consideration other constraints such as PV lifetime and power degradation, and battery state of health and power degradation, where some modifications can be made to power constraints, mentioned earlier in Table 1, accordingly, to give more realistic results especially for long term simulations.

Discussion and Analysis
In order to validate the obtained results confidently, it is important to compare the results of HHO algorithm in both cases with earlier research results that used the exact microgrid simulation model but with different optimization algorithms, such as fast surrogate-assisted particle swarm optimization (FSAPSO) and several modified PSO algorithms [31]. Table 5 shows the obtained results of the HHO algorithm along with results of other previously proposed algorithms in literature. Number of objective function evaluations (OFE) carried out till best solution is reached is indicated where it is equal to 1× number of agents × iterations. Additionally, the number of algorithm iterations needed till best solution is reached is indicated in the table along the trials number. Although there is a lack in information in [31], it can be noted that the proposed HHO algorithm successfully reaches lower operating costs with lower iterations. It is suggested as future work to take into consideration other constraints such as PV lifetime and power degradation, and battery state of health and power degradation, where some modifications can be made to power constraints, mentioned earlier in Table 1, accordingly, to give more realistic results especially for long term simulations.

Conclusions
Power generation cost optimization is considered one of the main problems facing modern microgrids having different energy resources and interconnections. In this paper, the generation cost optimization problem with all its constraints is explained in detail. A new optimization technique, Harris hawk optimization, is presented thoroughly as a proposed solution to the energy management optimization problem. The proposed algorithm is tested using MATLAB simulation on a benchmark microgrid and results are compared with previous research. The Harris hawk algorithm succeeded in achieving much-improved results as it produced the lowest generation cost in comparison with other widely used optimization algorithms. Obtained results indicate that the Harris hawk algorithm should be tested further on other applications in the field of power and energy management as it may form a strong and reliable computation algorithm.

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

Nomenclature
The following symbols are used in this manuscript: min x F minimum operation cost B Gi cost per kW of active power generated by i P t Gi generated active power of i t time P t Grid power exchange between utility grid and microgrid MP t market price of active power exchanged NT total hours NG total number of distributed generations ND total number of load levels MG Microgrid P t L D load active power at time t P t Gi,min minimum generation limit of generator i at time t P t Gi,max maximum generation limit of generator i at time t P t Grid,min minimum generation limit of utility grid at time t P t Grid,max maximum generation limit of utility grid at time t R t cheduled spinning reserve at time t W ess,t energy stored in battery W ess,min minimum energy to be stored in battery W ess,max maximum energy to be stored in battery P charge rate of charge P discharge rate of discharge P charge,max max charge P discharge,max max discharge ∆t predefined time for charge or discharge η charge charge efficiency η discharge discharge efficiency B G bid price C f uel price of fuel P G output power η G generation unit efficiency C inv payback amount P Gnom nominal power generation B G bid price AP annual energy production AC annual investment cost IC installation cost X(t + 1) location of hawk at iteration t+1 X rabbit (t) location of rabbit X rand (t) random choice of hawk X m average location of current population of hawks X i (t) position of hawk i at iteration t N total number of hawks E escaping energy of rabbit E 0 initial energy of rabbit T maximum number if iterations r percentage of successful rabbit escape ∆X(t) position difference between rabbit and hawks J jump strength r 1 to r 5 , q random numbers between 0 and 1 LF levy flight Y, Z next movement update D problem dimension S random vector with size 1 × D u, v levy flight random value between 0 and 1 β constant of value 1.5