The Hybrid Harris Hawks Optimizer-Arithmetic Optimization Algorithm: A New Hybrid Algorithm for Sizing Optimization and Design of Microgrids

This paper presents a new hybrid metaheuristic algorithm, the hybrid Harris Hawks Optimizer-Arithmetic Optimization Algorithm (hHHO-AOA), as we have named it. It is proposed for sizing optimization and design of autonomous microgrids. The proposed hybrid algorithm has been developed based on operating the Harris Hawks Optimizer (HHO) and the Arithmetic Optimization Algorithm (AOA) in a uniquely cooperative manner. The developed algorithm is expected to increase the solution accuracy by increasing the solution diversity during an optimization process. The performance is verified with the evaluation metrics and the well-known statistical tests. According to the Friedman ranking test, the new algorithm performs 77.9% better than HHO and 78.6% better than AOA. Similarly, the performance checked with the Wilcoxon signed-rank test has revealed a significant superiority in solution accuracy compared to HHO and AOA alone. Later, the hybrid algorithm is tested on a microgrid that consists of a photovoltaic (PV) system, a wind turbine (WT) system, a battery energy storage system (BESS), diesel generators (DGs), and a commercial type load. For the optimal capacity planning of these components, a problem in which the loss of power supply probability (LPSP) and the cost of energy (COE) are defined as the objective function is formulated. The optimization done by the proposed algorithm has produced the lowest LPSP and the COE along with the highest rate of renewable fraction (RF). In conclusion, it is demonstrated that the new hHHO-AOA has proved itself in designing reliable, economical, and eco-friendly autonomous microgrids in the best optimal way.


INDEX TERMS
P BAT _max maximum power produced by the BESS (kW). P DG output power of diesel generator (kW). P DG r rated output power of diesel generator (kW). P o output power of converter (%). P r rated power of converter (%). P PV output power of PV array (kW). P PV _max maximum power produced by the PV system (kW). P PV r rated output power of PV array (kW). P WT output power of WT (kW). P WT _max maximum power produced by the WT system (kW). P WT r rated output power of WT (kW). PV photovoltaic. r 1 aoa random numbers inside (0,1). r 2 aoa random numbers inside (0,1). r 3 aoa random numbers inside (0,1). r 1 hho random numbers inside (0,1). r 2 hho random numbers inside (0,1). r 3 hho random numbers inside (0,1). r 4 hho random numbers inside (0,1). r 5 hho random numbers inside (0,1). S a random vector by size 1 × D. MOA (t) math optimizer accelerated function value in iteration t. O&MC WT operation and maintenance cost of WT system ($). O&MC WT INV operation and maintenance cost of WT converter ($). O&MC yearly comp operation and maintenance cost of a component in yearly ($). q random numbers inside (0,1). R j amplitude of j th micro-cycle of DOD history. RC replacement cost of system ($). RC BAT replacement cost of BESS ($). RC BC replacement cost of bidirectional converter ($). RC comp replacement cost of a component ($). RC DG replacement cost of DG system ($). RC PV replacement cost of PV system ($). RC PV INV replacement cost of PV converter ($). RC WT replacement cost of WT system ($). RC WT INV replacement cost of WT converter ($). RF renewable fraction (%). SV salvage cost of system ($). SV BAT salvage value/cost of BESS ($). SV BC salvage value/cost of bidirectional converter ($). SV comp salvage value/cost of a component ($). SV DG salvage value/cost of DG system ($). SV PV salvage value/cost of PV system ($). SV PV INV salvage value/cost of PV converter ($). SV WT salvage value/cost of WT system ($). VOLUME  upper bound value of the j th position. UC BAT unit cost of BESS ($). UC BC unit cost of bidirectional converter ($). UC DG unit cost of DG system ($). UC PV unit cost of PV system ($). UC PV INV unit cost of PV converter ($). UC WT unit cost of WT system ($). UC WT INV unit cost of WT converter ($). v random numbers inside (0,1). V wind speed (m/s). V cut−in cut-in wind speed (m/s). V cut−out cut-out wind speed (m/s). V r rated wind speed (m/s). w i i th weighting coefficient. w 1 1 st weighting coefficient. w 2 2 nd weighting coefficient. WT wind turbine. X a set of candidate solutions. X (t) current position vector of hawks in iteration t. X (t + 1) current position vector of hawks in the next iteration t. x i,j (t) j th position of the i th solution in iteration t. x i,j (t + 1) j th position of the i th solution in the next iteration t. x i (t + 1) current i th solution in the next iteration t. X m (t) average position of the current population of hawks in iteration t. X rand (t) a randomly selected hawk from the current population in iteration t. X rabbit (t) position of rabbit in iteration t. α aoa a sensitive parameter. α p temperature coefficient of maximum PV power (%/ • C). β a default constant set to 1.5. X (t) difference between the position vector of the rabbit and the current location in iteration t. MOA max maximum value of the math optimizer accelerated function. MOA min minimum value of the math optimizer accelerated function. MOP (t) math optimizer probability function value in iteration t. N AD number of autonomous days (days). Metaheuristic algorithms, which are inspired by the nature, are iteratively working stochastic algorithms. They are general-purpose but rather smart algorithms that scan the search space intelligently with techniques called the exploration and the exploitation phases. Metaheuristic algorithms can be applied to any problem with almost no restrictions, where the primary aim is to reach to the global optimum [1]- [3]. Nevertheless, researchers are putting dedicated and continuing efforts to improve the solution accuracy, ease of programming, faster convergence rate, flexibility, and so on when solving optimization problems. Combining two or more algorithms to solve the similar problems has been the preferred and an effective way of obtaining the greater overall performance. Therefore, the starting point of hybrid concept is to increase the performance of individually working algorithms in a hybrid structure and to solve the problems with better accuracy and faster convergence rate in a more efficient manner [4]- [6]. In metaheuristic algorithms, while progressing from the beginning of an optimization process towards the end of optimization, premature and/or slow convergence problems are often encountered because of the low number of solution diversity. Especially the premature convergence may prevent obtaining the global optimum, resulting in a poor accuracy of the solution. Here, the premature convergence is commonly defined as a solution that is acceptable since it is close to the optimal, but it is not the global optimum.
Solution accuracy is considered an extremely important issue when finding the performance of long-term optimization problems. For example, finding the total present cost of a microgrid design configured to have a lifetime of 20 years or more, and the economic and technical benefits over those long periods are quite important and desirable. Another issue is the slow convergence, which is also related with the quality of the solution. Stagnation and repetitive flat solutions (plateaus) that are typically seen on convergence graphs indicate slow convergence. Slow convergence and solution accuracy are related to each other, and both change as a function of solution diversity. While a higher number of solution diversity contributes to the better solution accuracy, it results in slow convergence. On the other hand, a lower number of solution diversity may yield faster convergence but poor accuracy. Therefore, there has to be a trade-off between the slow convergence and the solution accuracy if an individual algorithm is used to solve a particular problem. Typically, the solution diversity is quite high at the beginning of a solution process, but it starts to decline towards the end. This suggests that convergence can be faster but algorithm may not obtain the global optimum.
In order to find the global optimum without staying stuck at a local optimum, several things can be tried. Optimal planning of algorithms that facilitates smart scheduling of transitions between the exploration and exploitation phases and setting up a proper balance between these phases are among the things that can be done. Proper scheduling and balancing can help to solving the accuracy problem to some extent; however, the most effective technique is to supplement the most suitable variety during a solution process. For this reason, hybridization is considered as a powerful strategy that produces the desired diversity in the process of finding the global optimum without the risk of being stuck at the local optimum point [7]. Consequently, the first objective of this study is to propose a hybrid algorithm to increase the solution accuracy by keeping the solution diversity high all over the solution process.
The authors' second motivation is the sizing optimization/optimal capacity planning of an off-grid microgrid. The correct capacity planning of microgrid components for autonomous structures operating based on on-site production and on-site consumption principle is significantly important to achieve the desired economic, technical, and environmental benefits over the lifetime of the project. Therefore, the solution accuracy becomes the most important parameter during the optimization process. With this study, we want to demonstrate that the benefits can be maximized if a microgrid is implemented with optimally sized components. Optimal capacity planning, in other words, the sizing optimization of microgrids, has received significant attention and been studied intensively by many researchers. Based on our literature review, we have noticed that more can be done to increase the benefits of the microgrid concept in energy generation. This is somewhat important for microgrids that do not have a standard structure, whether operated off-grid or placed in a location/region with unique the meteorological conditions and energy consumption characteristics. Therefore, we are motivated by the fact that the performance of an algorithm that is proven effective can still be improved in the capacity planning of a microgrid by using it in a collaborative hybrid form.
As discussed above, it is expected that the hybrid algorithm will improve mainly the solution accuracy and the computational speed by refining the search capability. For this purpose, in this study, a new hybrid hHHO-AOA algorithm based on the HHO and the AOA has been developed, implemented, and tested for sizing optimization/optimal capacity planning of an off-grid microgrid. Next section provides discussion about the state-of-the-art algorithms that have appeared in the reviewed literature.

B. LITERATURE REVIEW
A detailed literature search has been done on metaheuristic algorithms used in the sizing optimization of microgrids. During our literature review, we have focused on three aspects of the state-of-the-art sizing optimization technologies, the individual algorithms, the hybrid algorithms excluding the HHO, and the hybrid algorithms including the HHO.
When the studies that prefer the individual algorithms in capacity planning are investigated, one can see that so many of them have been tried. We have identified 35 different algorithms used in sizing optimization. The list is quite long but it is worth mentioning the names of all the algorithms and the associated references dedicated to sizing optimization. These are the Harris hawks optimization algorithm, the firefly algorithm [8], the moth-flame optimization algorithm, the genetic algorithm, the grey wolf optimizer, the particle swarm optimization, the ant colony optimization, the salp swarm algorithm, and the dragonfly algorithm [9], the grasshopper optimization algorithm [10], the cuckoo search, the simulated annealing, the harmony search, the jaya algorithm, the flower pollination algorithm, the brainstorm optimization in objective space, and the simplified squirrel search algorithm [11], the antlion optimizer [12], the evolutionary particle swarm optimization [13], the water cycle algorithm [14], the smell agent optimization [15], the artificial ecosystem optimization [16], the accelerated particle swarm optimization algorithm, the generalized evolutionary walk algorithm and the bat algorithm [17], the krill herd algorithm [18], the equilibrium optimizer, the artificial electric field algorithm, the sooty tern optimization algorithm [19], the seagull optimization algorithm, the modified farmland fertility algorithm [20], the crow search algorithm [21], the whale optimization algorithm, the gravitational search algorithm [22], and finally the multi-verse optimization algorithm [23].
In addition to the individual algorithms, the studies that have preferred the metaheuristic algorithms in a hybrid form are also investigated in detail. The hybridized particle swarm optimization-ant colony optimization [24], the hybrid particle swarm optimization-grey wolf optimizer [25], and the hybrid particle swarm-gravitational search algorithm [26] are VOLUME 10, 2022 among the hybrid algorithms formed based on the particle swarm optimization. Similarly, based on the simulated annealing algorithm, the hybrid harmony search-simulated annealing method [27], the hybrid chaotic search/harmony search/simulated annealing algorithm [28], and the hybrid search algorithm and the simulated annealing algorithm [29] have been developed. In addition, the hybrid shuffled frog leaping and pattern search algorithm [30], the hybrid soccer league competition-pattern search optimization algorithm [31], and the hybrid nelder-mead and cuckoo search algorithm [32] have been identified.
In the third phase of our literature search, we have investigated the hybrid algorithms combined with the HHO algorithm, which is also the subject of this study. The HHO algorithm and the sine cosine algorithm are used in hybrid form in three studies. Kamboj et al. proposed a hybrid Harris hawks-sine cosine algorithm to free the HHO algorithm stuck in the local search space and speed up the global search process. The proposed algorithm was tested on the standard benchmark problems. With this hybrid algorithm, the HHO and the sine cosine algorithms were evaluated with the worst fitness value, the best fitness value, the average fitness value, and the standard deviation metrics. The proposed hybrid algorithm achieves superior results compared to the individual algorithms that it is compared with [33].
Fu et al. proposed the mutation sine cosine and Harris hawks optimization algorithm for fault detection in bearings. This hybrid algorithm is used to optimize the parameters of the support vector machine. With the proposed method and hybrid algorithm, superior diagnostic results have been obtained in the detection of faults in bearings and efficiency has been increased [34]. Hussain et al. proposed the hybrid sine-cosine Harris hawks optimization algorithm for future selection. In order to increase the exploitation performance of the HHO algorithm, the sine-cosine algorithm is integrated into the HHO algorithm. The accuracy rate is improved by increasing the convergence rate with the hybrid algorithm [35].
Elaziz et al. proposed a hybrid Harris hawks-moth-flame optimization algorithm to improve the exploration capability of the HHO. The proposed algorithm has been tested on the specified functions and engineering problems, and it was concluded that it is successful compared with the individual algorithms selected by the authors [36]. Li et al. used a combination of the elite evolutionary strategy and the HHO algorithm to improve the HHO algorithm's affinity for being stuck to local optima, and they proposed the hybrid elite evolutionary strategy Harris hawks optimization algorithm to solve this problem. Based on the test performed with the benchmark functions, it was shown that the convergence speed and optimization performance are improved compared to the individual algorithms [37].
Boa et al. studied the color image segmentation, and proposed the hybrid Harris hawks optimization-differential evolution algorithm. It is aimed to cope with the situation where color images contain too much information and the number of thresholds is high, and to perform the feature extraction process with the best performance. The results of the hybrid algorithm were evaluated in terms of standard deviation, the average fitness values, the structure similarity and the feature similarity index, and the peak signal to noise ratio. In addition, the individual algorithms and other selected algorithms were statistically compared with the Wilcoxon rank sum test. They obtained successful results with the hybrid algorithm [38]. Al-Wajih et al. proposed the hybrid binary gray wolf optimizer and the Harris hawks optimization algorithm to solve the problem of the gray wolf algorithm being stuck to local optimum point. When the results of the hybrid algorithm tested with the benchmark functions are compared with the binary gray wolf optimizer, it was seen that there is an improvement in accuracy, computation time, and the overall performance [39].
The literature search presented in the following paragraphs are dedicated to show the link between the state-of-theart HHO based hybrid algorithms and our proposed hybrid algorithm. In reference [33], researchers use the sine cosine algorithm in hybrid form with the HHO to enhance the exploration and exploitation phases, where they run the sine cosine algorithm after the HHO sequentially. The motivation in the hybrid approach in [34] is to update position in order to increase the diversity. They use the mutation operator to achieve a periodic updating. The HHO search agents that are positioned in the upper layer are updated by the sine cosine algorithm positioned in the lower layer. Basically, the objective has been to stay away from the local optimum point. Similarly, the study in [35] uses the sine cosine algorithm combined with the HHO to enhance the exploration phase. In contrast, they add the delta factor to enhance the exploitation process. The main objective in this study is to create a more dynamic exploitation phase. In [36], the moth-flame optimization algorithm is used in hybrid form with the HHO to increase the search capability of the HHO in exploration phase. The motivation here is to improve the overall performance. With the same motivation, the researchers in [37] use the elite evolution strategy with the HHO algorithm in order to improve the exploration capability of the HHO. In addition to the algorithms that prefer execution of individual algorithms sequentially in a certain order, the reference [38] suggests a strategy to run the differential evolution and the HHO in parallel form. According to this strategy, the number of populations is divided into half. Then, the HHO and differential evolution algorithms that have the same number of populations are run in parallel form. Later, the divided groups are combined again, and the best solution is sought by comparing the fitness values. Moreover, in [39], the binary grey wolf optimizer and the HHO are combined to avoid local optimum solutions and overcome premature convergence problems. Here, for obtaining the optimal solution, the authors make the HHO responsible for the exploration phase and the grey wolf optimizer responsible for the exploitation phase.
Based on the reviewed literature, it has been comprehended that the algorithms that perform well individually can do much better in a hybrid form. In brief, the major objectives of hybrid developments are avoiding the local optimum and increasing the solution accuracy and the convergence speed. It is also acknowledged that the superior results obtained by hybrid algorithms are not by chance alone, rather require more research and development. What we have noticed in those studies is that each is using unique update rules and operating schemes. Therefore, finding the right combination of algorithms, better update schemes, smart scheduling, and the determining of operating modes are all open to optimization within itself. Moreover, design of a superior hybrid structure requires the knowledge of pros and cons of each individual algorithm depending on the problem and the understanding of the best fit when combined [40]- [42]. We can conclude that our proposed hybrid structure is brand new, well organized, and configured with smart scheduling methods developed to solve advanced optimization problems such as off-grid microgrid designs.

C. CONTRIBUTIONS
The major contributions of this study are given below: • It is the first time that the HHO algorithm and the AOA are used in a hybrid form.
• A new hybrid metaheuristic algorithm, which is named as the hybrid Harris Hawks Optimizer-Arithmetic Optimization Algorithm (hHHO-AOA), is brought to the literature.
• For the best use of the proposed hybrid hHHO-AOA, a unique scenario is developed and presented.
• Performance tested with the Friedman ranking test method places the hHHO-AOA at the top with 77.9% better performance than the HHO and 78.6% better than the AOA alone.
• Performance tested with the Wilcoxon signed-rank test method also distinguishes the hHHO-AOA from the other individual algorithms in providing the highest solution accuracy in optimization.

D. ORGANIZATION OF THE ARTICLE
This article is organized into VIII Sections. In Section I, we state the problem, provide our motivation, describe the research gap, present a detailed literature review, and list the major contributions of the research. Section II explains the kinds of components used in the proposed microgrid structure and the mathematical model of each component for analysis purposes. Section III presents the detailed description of the problem. Section IV describes the proposed hHHO-AOA algorithm for the solution of the problem. In Section V, the results of the HHO, the AOA, and the hHHO-AOA algorithms are presented and their performances are compared against each other. In Section VI, the simulation results of the microgrid designed based on the hHHO-AOA algorithm is presented. Section VII presents the sensitivity analysis. Finally, Section VIII presents the conclusions and the future work.

II. COMPONENTS OF THE PROPOSED MICROGRID STRUCTURE
The structure of the proposed off-grid microgrid and its components is shown in Fig. 1. It consists of renewable energy sources, energy storage systems, backup energy sources, interfacing power electronic converters, an energy management system, and a real commercial type load. The following subsections explain the microgrid components in more detail.

A. CHARACTERISTICS OF THE LOAD AND METEOROLOGICAL DATA
As aforementioned, the selected load is an existing commercial load located inside the Gazi University Campus in Ankara, Turkey, at 39 • 46 46.6 N and 32 • 48 29.2 E latitude and longitude coordinates. It is the Technology Park Building hosting 128 companies at different sizes and running business in different industrial branches. The diversity is therefore resulting in a highly variable demand. In addition to the demand data, meteorological data had also been recorded by the meteorology station installed on the roof of the same building. The demand data, solar irradiance, wind speed, and temperature data belonging to this building are all given in Fig. 2-5, respectively. It is worth mentioning that data are recent and belong to 2018. We believe that the use of a yearlong data in a capacity planning is important in order to maximize the economic and technical benefits from a microgrid design. The yearlong data include all kinds of variations in consumption, from morning to evening, from weekdays to weekends, similarly variations in meteorological conditions, from daytime to nighttime, and from season to season. An investment can only be guaranteed to keep up with the project targets if a realistic input is used at the planning and designing stage. It is worth noting that the load and meteorological data presented in this Section will be used as the input data for the sizing optimization of the proposed microgrid. The following paragraphs explain the data in more details. In order to reflect the characteristics of the given data clearly, they are presented as the annual data, obtained by recording on an hourly basis, monthly average of the recoded data, and daily average of the recorded yearly data.
The load data seen in Fig. 2 have a minimum value of 29.54 kW, an average of 132.07 kW, and a maximum value of 349.62 kW. While the month with the lowest monthly average energy consumption is December with 109.41 kW, the month with the highest energy consumption is April with 150.87 kW. While the average daily energy consumption is the lowest at 06:00, it is the highest at 12:00. The solar irradiance data are seen in Fig. 3. It has an annual average of 0.18 kW/m 2 and 1.05 kW/m 2 maximum. Monthly average solar irradiance is the lowest in January, and the highest in June. The daily average is the highest at 13:00.
The second meteorological data, the wind speed, are given in Fig  with the lowest monthly average wind speed is October with 2.50 m/s, the month with the highest is March with 6.76 m/s. While the daily average wind speed is the lowest at 07:00, the fastest wind is at 17:00.
The last meteorological data, which is the temperature data, are given in Fig. 5. The temperature data have a minimum of −10.03 • C, an average of 13 • C, and a maximum of 34.70 • C. The month with the lowest monthly average temperature is January with 1.77 • C, while the highest month is August with 23.79 • C. While the lowest daily average temperature is at 07:00, the highest temperature is obtained at 16:00.

B. RENEWABLE ENERGY SOURCES
Because solar and wind energy are regional sources with high potential, they are chosen as the renewable energy sources of the proposed microgrid. The power obtained from a PV array is given in equation (1), as shown at the bottom of the page. The output power varies depending on the rated power of the selected panel, the derating factor, which is a parameter reflecting the actual environmental conditions, solar irradiance, ambient temperature, and the cell temperature of the PV modules [43], [44].
Similarly, the output power obtained from a WT is given in equation (2), as shown at the bottom of the page, which allows investigation of power in four regions. In region 1, power is zero since WT cannot work due to insufficient wind because the wind speed is below the cut-in wind speed. In region 2, WT can produce power since the wind speed is in the range of cut-in wind speed and nominal wind speed. The power obtained from WT in this region is determined by the available wind speed, cut-in wind speed, cut-out wind speed, and the rated power of the selected WT. The raise of the wind speed to the range of rated wind speed and cut-out wind speed corresponds to region 3 where the rated power is produced. In the last region, region 4, the wind speed exceeds the cutout wind speed, and WT is stopped due to safety reasons. The rated power of a WT given in equation (3) is determined by the maximum power coefficient, air density, the area covered by the rotor blades depending on the rotor diameter, and the rated wind speed [45], [46].

C. ENERGY STORAGE SYSTEM
In an off-grid microgrid, an energy storage system is used to store the surplus energy produced by the renewable energy sources to back it up when there is a need. A battery pack is used as the energy storage system. The BESS capacity is determined according to the load demand, the number of days that the microgrid is expected to operate autonomously, the depth of discharge, the converter efficiency, and the battery efficiency depending on the selected battery type and the model. The formula that gives the capacity of a battery-based energy storage system is given in equation (4) [47], [48].

D. BACKUP ENERGY SOURCE
In off-grid microgrids, backup energy sources in addition to the renewable energy sources and storage are important since they contribute to the uninterrupted operation. For this reason, DGs are used as the secondary energy source in this design. The formula that finds the fuel consumption of a DG in terms of the output power, the rated power of the selected DG, and the fuel consumption coefficients is given in equation (5) [49], [50].

E. INTERFACING CONVERTERS
Properly sized and configured power electronics converters act as the interfaces between the microgrid components (PV, WT, and BESS) and the AC bus. The efficiency calculation of these converters is given in equations (6)-(9) [51].

F. ENERGY MANAGEMENT SYSTEM
The final component of the proposed off-grid microgrid is the energy management system. This system sets up the communication between renewable energy sources, storage units, and the backup energy sources in order to meet the energy demand of the load at any time. This rule-based energy management system organizes the microgrid for the four states named as the storage state, the storage and dumped load state, the BESS and DG support state, and finally the DG support state.
• Storage State: Energy produced by the renewables is more than the demand, and the surplus energy is stored in the BESS.
• Storage and Dumped Load State: Excess energy produced by the renewables is stored in the BESS. However, when the BESS is charged to full capacity, the extra energy is wasted in the dump load.
• BESS and DG Support State: It is insufficient to meet the energy demand by renewable energy sources alone. In this case, energy is supplied to the load primarily by the BESS first, and when a low BESS capacity is detected, energy is supplied by the BESS and the DG together.
• DG Support State: It is insufficient to meet the energy demand by the renewables alone, and the BESS capacity is also at its minimum, then the energy is supplied to the load by the DG alone.

III. PROBLEM DEFINITION
The microgrid components and their mathematical models had been given in section II. In this section, the sizing optimization problem of the off-grid microgrid is defined and described. Our goal by defining the sizing optimization problem is to determine the optimal capacities of the components for the most economical and reliable microgrid design. In line with this goal, the objective function, decision variables, constraints, and sub-functions that make up this objective function are all determined. Following paragraphs explain these parameters in more detail.

A. OBJECTIVE FUNCTION, DECISION VARIABLES AND CONSTRAINTS
Multi-objective optimization problems are converted into single-objective optimization problem with the weighted sum method. The general representation of the scalarized objective function of the weight sum method is given in equations (10)-(11) [52]. Where, F WS represents the general-objective function. A weighting factor (w i ) is assigned to each objective function. Each function is then multiplied by the associated weighting factor as shown in (10). After adding all the elements, the objectives functions that were more than one are converted into a single scalar objective function with the weighted sum method [53].
In our sizing optimization problem, two objective functions, one being the economic and the other one being the reliability indicator, are used. The representation of the two individual objective functions with the weighted sum method is given in equation (12), and the scalarized objective function adapted for this problem is given in equation (13). When these two equations are examined, one can see that f 1 (x) is the first objective function, and represents the COE that is an economic indicator, f 2 (x) is the second objective function, and represents the LPSP that is a reliability indicator. x DV is the vector of the decision variables of the problem. Four decision variables, PV, WT, autonomous days, and DG numbers listed in equation (14) have been used in this problem. The constraints that determine the lower and upper limits of these decision variables are given in equations (15)- (18). Two weighting coefficients are used for our two objective functions. The sum of these weights must be equal to 1. Since both indicators have equal importance in the design of this microgrid, the weights are assigned a value of 0.5.
The mathematical models of the economic and reliability objective functions are discussed in detail in the following subsections separately.

B. FIRST OBJECTIVE FUNCTION -ECONOMIC
The first objective function based on the economic indicator is the cost of energy. The COE indicator allows comparison of superiority of projects that utilize different energy sources, capacity of components, technology, and lifetime. The COE, which is defined as the average cost per kWh of energy VOLUME 10, 2022 production, is given in equation (19). Total net present cost (TNPC) is calculated by subtracting all revenues obtained from the sum of all expenses incurred during the life of a project. As seen in equation (20), the total expenses of a system as being the investment cost, replacement cost, operation and maintenance cost, and fuel cost of DGs, the revenues, on the other hand, are the salvage value of the system. The TNPC is converted to annual cost when it is multiplied by the capital recovery factor (CRF) as given in equation (21) [54]- [56].
The investment cost (IC) includes the investment cost of PV, WT, BESS, DG, and the cost of power electronics (interfacing) converters. The formula that finds the total investment cost is given in equation (22). Equations (23)- (29) give the individual investment costs of all components listed in equation (22). When calculating the investment costs, the information containing the quantity, rated power, and the unit cost of the components are used for PV, WT, and DG; similarly, the maximum energy capacity and the unit cost of the components are taken into account for the BESS.
The capacity of the power electronics converters is decided based on the maximum power of the component they are interfacing. In this study, we added a proper safety margin above the maximum ratings when sizing the converters. After sizing of the converters for PV, WT, and the BESS, the investment cost of each converter is determined based on multiplying the particular rated power by the unit cost.
IC PV _INV = P PV _max · ℵ PV _INV · UC PV _INV (28) IC WT _INV = P WT _max · ℵ WT _INV · UC WT _INV (29) As in the investment cost, the replacement cost includes the individual replacement costs of PV, WT, BESS, DG, and the interfacing converters as shown in equation (30). Because PV and WT systems have long lifetimes, their lifetimes are generally considered equal to the project lifetime. However, especially in off-grid microgrids, BESSs and DGs are among the components to be renewed according to the usage conditions. The replacement costs of a BESS and a DG are given in equations (31)-(32), respectively. When calculating the replacement costs, the number of replacements of a component (N rep ) is obtained over the lifetime of the project (n p ) and the lifetime of the component (n comp ) as given in equation (33) [57]. Later, the cash flow in the year that a component was replaced must be carried over to present time. With the discount factor (f d ) given in equation (34), the monetary value of the component in the year it was renewed is carried over to today's value and therefore its present value is calculated [58].
N rep = n p n comp (33) The operation and maintenance cost, which is given in equation (35) The salvage value means the remaining life of a component when it is referenced to the project life. It is related to the renewal of life. The salvage cost of a system is given in equation (37) and the salvage value of a component is given in equation (38) [59], [60].
SV comp = RC comp · n rem_comp n comp (38) C. SECOND OBJECTIVE FUNCTION -RELIABILITY The second objective function based on the reliability indicator is the loss of power supply probability. As seen in equation (39), the LPSP is an index representing the ratio of the energy that cannot be supplied to the load. It takes a value in the range of 0 to 1 or 0% to 100%. 0% LPSP means that, the energy demand of the load is fully met, and reliability is at its maximum. If the LPSP is 100%, the system failed to energize the load and the system is unreliable [61], [62].

D. RENEWABLE FRACTION
The utilization of renewable energy sources in a microgrid is expressed as the renewable fraction (RF). The RF is given in equation (40), which is a measure of the energy generated from renewable energy sources in a system and transferred to the load [63].

IV. DEVELOPMENT OF HYBRID HHO-AOA ALGORITHM
In this section, we explain the development process of proposed hybrid algorithm, which we have named it as the hHHO-AOA algorithm. It is based on operating the HHO and the AOA algorithms in a unique way of cooperation. For a clear understanding of the development process, we initially treat each algorithm separately and explain their structures and the solution search processes individually. At the end of the Section IV, we will present the concept and the design of the hHHO-AOA algorithm with its scenarios.

A. HARRIS HAWKS OPTIMIZATION ALGORITHM
The HHO algorithm was developed by Ali Asghar Heidari et al. in 2019. It is a population-based metaheuristic algorithm inspired by the nature of the Harris' hawks. The concept of the HHO algorithm is based on the process of Harris' hawks catching a prey in a cooperative manner. These processes and their mathematical modeling divided into phases are explained below under three subsections: the exploration phase (phase I), the transition from exploration to exploitation phase, and the exploitation phase (phase II). The flowchart of the algorithm is shown in Fig. 6.

1) PHASE I -EXPLORATION
In the HHO algorithm, Harris' hawks are the candidate solutions used to hunt, in other words to reach the optimum solution. Harris' hawks have sharp eyes that help them to track and locate their prey immediately with precision. However, waiting, observing and tracking a prey in an environment like a desert can take hours. As shown in equation (41), during a hunting event, the process of the Harris' hawks perching in a random location and detecting their prey is modelled uniquely according to the status of q, which is either q < 0.5 or q ≥ 0.5. For the case q < 0.5, the Harris' hawks are considered to perch near the family members and the prey (rabbit), and for the case q ≥ 0.5, they are considered to perch on random tall trees within the range. X (t + 1) in equation (41) gives the next positions of hawks in the iteration process. In order to imitate the natural behavior of the hawks in its simplest form, the average position is used for the position update within the search range. The average position  of the current population of hawks in iteration t is given in equation (42).

2) TRANSITION FROM PHASE I TO PHASE II
The transition from the exploration phase to the exploitation phase is related to the energy of the prey chased by the Harris' hawks. The decrease in the energy of the prey escaping from the Harris' hawks over time is modeled by equation (43). The initial energy of the prey E 0 takes random values between -1 and 1 in each iteration. Increasing of E 0 from 0 to 1 means that the prey (rabbit) is getting stronger, and decreasing from 0 to -1 means that the rabbit is physically weaker. When evaluated during an iteration, the escape energy E shows a decreasing trend throughout the iteration.

3) PHASE II -EXPLOITATION
In the exploitation phase, surprise attacks (surprise pounce) are launched to the prey detected by Harris' hawks during the exploration phase. Like in the real life where the hunter chases the prey and the prey escapes, surprise pounce or seven kills strategies by the hawks are performed in different forms. The chance of escaping of the prey, which is constantly trying to escape from the hunter, is associated with the parameter r.
If r ≥ 0.5, it means that the prey could not escape during the surprise attack, and if r < 0.5, it means that the prey has escaped from the hunter successfully. The exploitation phase can be modeled by four strategies executed by the Harris' hawks, namely, the soft besiege, the hard besiege, the soft besiege with progressive rapid dives, and the hard besiege with progressive rapid dives. The following paragraphs explain each strategy.
a: r ≥ 0.5 AND E ≥ 0.5-SOFT BESIEGE When r ≥ 0.5 and |E| ≥ 0.5, we say that the rabbit has the energy to escape from the hawks. However, due to the soft besiege strategies of the Harris' hawks, the rabbit's allegedly misleading and deceptive moves and leaps do not work. He is soon caught by the hawks. The present position of the hawks, the difference between the position of the rabbit and the present location, and the random jump strength of the rabbit are given by equations (44)-(46), respectively.
In case of r ≥ 0.5 and |E| < 0.5, it is considered that the rabbit has too low energy to escape from the Harris' hawks. The rabbit is surrounded by the hard besiege strategy of the hawks and caught by a surprise attack. The present position of the hawks in the case of hard siege is given in equation (47).
The case of r < 0.5 and |E| ≥ 0.5 suggest that the rabbit has enough energy to escape from the Harris' hawks. At this stage, the strategy of the rabbits moving away from the place by zigzagging with deceptive movements is modeled with levy flight (LF). Hawks perform soft besiege before the surprise attack and make the evaluation and the decision of their next move based on equation (48). If the previous dive decisions were not good, they make rapid and irregular dives on the prey. Harris' hawks continue these dives until they catch the rabbit. This process is done using the LF pattern given in equation (49). The LF function is given in equation (50). The position update of the hawks for this case is performed by equation (51).

5-HARD BESIEGE WITH PROGRESSIVE RAPID DIVES
For the case r < 0.5 and |E| < 0.5, the rabbit does not have enough energy to escape from the Harris' hawks. Before a surprise attack, similar to soft besiege, hawks perform hard besiege to catch and kill the prey. For this to happen, it is aimed to reduce the distance between the prey and the average location of the hawks. The rules of this case are given in equations (52)-(54) [64].

B. ARITHMETIC OPTIMIZATION ALGORITHM
The AOA algorithm was proposed by Laith Abualigah et al. in 2021. It is a population-based metaheuristic algorithm inspired by the distribution behavior of arithmetic operators in mathematics consisting of addition, subtraction, multiplication, and division. With the AOA algorithm, the use of arithmetic operators in problem solving is carried out in three phases: the initialization, the exploration, and the exploitation phases. The flowchart of the AOA is given in Fig. 7. Following subsections explain the phases of the AOA in detail.

1) INITIALIZATION
In the first step of the AOA algorithm, the optimization process is started with the candidate solutions X given in equation (55). The optimization process, which is initiated with randomly generated candidate solutions, continues until the best result is obtained or the termination criterion is reached. Then, the choice of exploration and exploitation phases is judged with the math optimizer accelerated function (MOA) given in equation (56), which expresses a coefficient, and based on the rule given in equation (57). Here, r 1 aoa are the random numbers inside 0 and 1. When r 1 aoa ≥ MOA, it works in the exploration phase, and if r 1 aoa < MOA, it switches into the exploitation phase.

2) PHASE I -EXPLORATION
Among arithmetic operators, the multiplication and the division tend to take high-distributed values or decisions due to their high dispersion characteristics. Therefore, they cannot easily converge to the target. In the exploration phase, the multiplication and division operators are used to perform random exploration of different regions to find the best results. In case if r 2 aoa > 0.5, division operator, and in other cases the multiplication operator comes into play and the position update equations given in equation (58), as shown at the bottom of the next page, are used. The math optimizer probability function (MOP), which is a coefficient, is given in equation (59), as shown at the bottom of the next page.

3) PHASE II -EXPLOITATION
Among arithmetic operators, the addition and the subtraction tend to take high dense values or decisions due to their low dispersion characteristics. Therefore, unlike the multiplication and division operators, they can easily converge to the target. For this reason, multiplication and division are used in the random exploration of different regions, that is, in the exploration phase, while the addition and subtraction operators are used in the exploitation phase for in-depth exploration. In the exploitation phase, position update equations given in equation (60), as shown at the bottom of the page, are used. For the case when r 3 aoa > 0.5, the subtraction operator is used, and otherwise, the addition operator is employed [65].

C. DEVELOPMENT AND DESIGN OF PROPOSED HYBRID HHO-AOA ALGORITHM
Within the scope of this study, as aforementioned, we have developed a new hybrid algorithm by combining the HHO and the AOA algorithms with a unique strategy that achieves working cooperatively in an optimal way. By the proposed hybrid hHHO-AOA algorithm, as we have named it, we aim to solve the problems with better solution accuracy. By hybridization, we want to demonstrate that performance and accuracy of each algorithm can be improved beyond individual performances. The flowchart of the hHHO-AOA is given in Fig. 8. The following paragraphs explain each step of the flowchart in detail.
• First, the parameters of the HHO algorithm, which are the number of populations and iterations, are defined, and the random populations are initialized.
• Meantime, the parameters of the AOA algorithm, which are the number of solutions and iterations, are defined, and the random solution's positions are initialized.
• The HHO algorithm is run for the specified iteration time. Initial energy, jump strength, and escaping energy are all updated at each iteration. According to the x i,j (t + 1) = best x j − MOP · UB j − LB j µ + LB j r 3 aoa > 0.5 Subtraction best x j + MOP · UB j − LB j µ + LB j r 3 aoa ≤ 0.5 Addition Flowchart of hHHO-AOA. VOLUME 10, 2022 defined rules that constitute the basis of the HHO algorithm, exploration and exploitation phases are studied. Then the fitness values of the hawks are calculated and the best location is found.
• With the decrease of solution diversity, repeating flat results are obtained in the convergence graph. Since the diversity is high at the beginning of the solution process, the success rate in reaching to the optimal solution is also high. However, towards the end of the solution process, diversity decreases and the success rate in reaching to the optimal solution decreases. The most important point here is to provide the most appropriate variety and diversity in the optimization process. For this reason, in this study, a unique hybridization strategy is developed and applied to ensure the most appropriate and the richest diversity. The fitness values obtained by the HHO algorithm have been preferred and used to realize this objective. Therefore, when the number of repeating fitness values in the convergence graph reaches to a predetermined number, we make a transition from HHO to AOA. In Table 1, the list of transition scenarios from HHO to AOA are given. As shown in Table 1, for 20, 30, and 40 repeating fitness values of HHO, AOA is run for 3, 50, and 200 iterations.
• The AOA should be able to continue seeking solutions from where the HHO had left off. Moreover, in this transition process, the last best solution (location) found by the HHO should be assigned to the best solution value of the AOA algorithm. Therefore, the AOA algorithm will start looking for a new solution from where the HHO had left off.
• The AOA algorithm is run through the specified iteration period. The MOA and MOP are updated at each iteration. According to the defined rules that constitute the basis of the AOA algorithm, exploration and exploitation phases are studied. Then, the fitness values of the solutions are calculated and the best solution is found.
• When the number of iterations determined for the AOA algorithm is reached, hybrid algorithm assigns the best solution value to the best location value of the HHO algorithm and transfers the solution search process to the HHO algorithm. Thus, over the maximum iteration time of the HHO, the HHO and AOA algorithms work collaboratively.

V. RESULTS OF SIZING OPTIMIZATION OF THE OFF-GRID MICROGRID
Here, we present the detail performance analysis of the proposed hybrid algorithm via the sizing optimization of the autonomous microgrid structure described in Section II. Moreover, the performance analysis of the HHO and the AOA algorithms along with the proposed hybrid algorithm are presented. Therefore, computer specifications, definition of input parameters, results and comparison of sizing optimization based on evaluation metrics and statistical tests, and convergence curves are all given in this section.

A. COMPUTER SPECIFICATIONS
For the optimization problem, MATLAB software is used and algorithms are coded as M-files. The computer used to run the algorithms has an Intel(R) Core (TM) i7-4790 CPU@3.60GHz 24GB.

B. DEFINITION OF INPUT PARAMETERS
Before optimization, it is necessary to define both the microgrid components and the input parameters of the algorithms. While the technical and economic parameters of the components are effective in obtaining the most economical microgrid design, the optimization parameters are equally effective on the performance of the algorithms. The technical and economic parameters are given in Table 2, and the optimization parameters of the algorithms are given in Table 3. Moreover, before the optimization, we need to define the decision variables and their limits. The four decision variables have been defined for this study, which are also called the search space components, are listed in Table 4. For convenience, the component limits as shown in Table 4 are given by their actual units. However, in an optimization process, the unit of decision variables is generally taken as numbers such as the number of PV arrays, the number of WTs, the number of autonomous days (corresponding to certain BESS capacity), and the number of DGs. Once the optimization is completed and optimum numbers of components are determined, then the capacities are converted from numbers to actual units after they are multiplied by the rated values.

C. ASSESMENT OF SIZING OPTIMIZATION BASED ON EVALUATION METRICS AND STATISTICAL TESTS
The compared algorithms are set to 30 runs, each run with 500 iterations, and with 15 populations. All results have been recorded and a total of 15000 function evaluations are performed. Then, we have made a performance comparison based on two statistical methods and four evaluation metrics. Friedman ranking test and Wilcoxon signed-rank test have been used as the statistical methods. The mean, the standard deviation, the minimum, and the maximum are chosen as the evaluation metrics. Table 5 shows the optimization results of the tested algorithms. When the evaluation metrics results of the scalarized objective function and the rank numbers (from smallest to the largest) are examined, it can be seen the scenario-based operated hHHO-AOA algorithm is more successful than the HHO and the AOA individually. Especially the hHHO-AOA based on scenario 1, which is operating the HHO running 3 iterations of the AOA when the number of repeating fitness values is 20 (see Table 1), has the lowest values of the average, the minimum, and the maximum. When the algorithms are evaluated over the mean value only, the hybrid structure that is running the AOA algorithm with 3 iterations stands out. According to Table 5, the success sequence goes from 20, 30, and 40 repetitions; each sequence corresponds  to scenario numbers 4, 7, and 3, respectively. The similar success sequence is obtained for 20, 30, and 40 repetitions in 200 iterations and 50 iterations of the AOA corresponding to scenario numbers 6, 9, 2, 5, 8, 10, and 11, respectively, as given in Table 5.
The hHHO-AOA based on scenario 1, which has achieved the best mean result, has obtained 0.72% better results than the HHO algorithm, which is the last in the ranking list, and 0.80% better results than the AOA algorithm. In addition, the lowest standard deviation value is obtained again with the hybrid structure in which the 200 repetitions and 3 iterations of the AOA algorithm were in use. With the hHHO-AOA based on scenario 7, which has achieved the best standard deviation result, 60.17% more consistent results have been obtained compared to the HHO algorithm, which was placed at the bottom of the ranking, and 67.33% compared to the AOA algorithm. In brief, it can be concluded that the proposed hybrid algorithm performs much better than individual algorithms.
Computational times of algorithms are related to their nature and the computational complexity. The computational complexity of the HHO, AOA, and the hHHO-AOA algorithms are given in equations (61)-(63), respectively. In the HHO algorithm, hawks are updated after initialization and fitness value evaluation processes. The computational complexity involving these processes depends on the number of hawks, the maximum iteration time, and the size of the problem. In the AOA algorithm, similar to the HHO, the solutions are updated after initialization and fitness value evaluation processes. Again, the computational complexity is calculated depending on the number of solutions, the maximum iteration time, and the size of the problem [64], [65]. In the hybrid algorithm, the computational complexity is the sum of the computational times of the HHO working as the VOLUME 10, 2022 base algorithm and the AOA as the collaborator.
The running times of the algorithms are given in Table 6. It can be seen that the obtained data are suitable with the computational complexity. According to the results presented in Table 5, even though the AOA comes the last in finding the scalarized objective values, it is the algorithm with the lowest computation time. The hHHO-AOA based on scenario 1, which finds the best objective values, has the longest computation time. Obtaining of the best result in the longest time should be interpreted according to the uniqueness of the problem. The projected life of microgrids is 20 to 25 years. Such a long life suggests that a vise planning is very important to achieve the highest benefit from the investment at the end. In this case, capacity planning should focus on the most accurate result rather than the fastest calculated result. Therefore, in the sizing optimization, the result that is the best and the most accurate is referenced instead of the calculation time.
Friedman ranking test is a non-parametric statistical test method. With this test, which is equivalent to the variance analysis, groups can be evaluated. It is used to understand and compare the overall ranking of algorithms [66], [67]. Another non-parametric statistical method is the Wilcoxon signed-rank test. It is used to compare the performance of two algorithms that are selected to solve a particular problem in terms of statistical perspective. With this test, the algorithm that yields the higher statistical success is determined [68]. In this study, Friedman and Wilcoxon statistical tests both with 5% degree of significance and with 0.05 significant statistical value have been used to rank the 11 algorithms listed in Table 5 and to detect the major differences between their results. It is important to note that we have actually mentioned the names of only three algorithms throughout the paper. They are the proposed hybrid algorithm (the hHHO-AOA), the HHO alone, and the AOA alone. Since this paper is also devoted to demonstrating the development phases of the proposed hybrid algorithm, we decided to consider each case of the hybrid structure associated with a scenario as a unique algorithm. Therefore, the number of total algorithms evaluated has become effectively 11, and the 9 of them being the hybrid one with 9 distinct scenarios.  The results of the Friedman ranking test are given in Table 7. The hHHO-AOA based on scenario 1 is ranked the first and again the hHHO-AOA based on scenario 7 is ranked the second. The AOA algorithm took the 10th place among 11 algorithms and the HHO algorithm alone took the last place. The results in Table 7 reveal that the hHHO-AOA based on scenario 1 is 41.17% more successful from scenario 7, which is the second best, 47.36% from scenario 4, 57.74% from scenario 3, 59.18% from Scenario 6, 66.38% from scenario 9, 71.83% from scenario 5, 72.91% from scenario 2, 75.99% from scenario 8, 78.60% from the AOA, and 77.90% from the HHO. In conclusion, the Friedman test has clearly revealed that the proposed hybrid algorithm performs much better than the individual algorithms.
In addition to the successful tests made with the Friedman test method, Wilcoxon test is used to detect the major differences between the results of the algorithms and to tell whether they differ from each other significantly. The results of the Wilcoxon test are given in Table 8. In order to tell that a result is better than the one compared; the p-values must be less than 0.05. The statistical findings have revealed a significant difference between the hHHO-AOA based on scenario 1 in solving this problem compared to both the HHO and the AOA. Based on the Wilcoxon test, it can be concluded  again that the proposed hybrid algorithm yields better results than the individual algorithms.

D. CONVERGENCE CURVES
Based on statistical test results and overall evaluations, it has been concluded that the hHHO-AOA based on scenario 1 outperforms the HHO and the AOA, and the other hybrid structures described earlier. Fig. 9 shows the convergence curves expressing the process of finding the scalarized objective values of the hybrid algorithm based on scenario 1, the HHO, and the AOA. As it can be seen from the graph, the hybrid algorithm is more successful than the individual algorithms. As pointed out in previous Sections, the reason for this is that the collaborative hybrid algorithm works with higher diversity during optimization.
In conclusion, the benefits of obtaining high diversity during optimization by the developed strategy have been demonstrated in three ways: comparisons made with the evaluation criteria, the statistical tests, and with the convergence path that the hybrid algorithm follows in reaching the conclusion.

VI. SIMULATION OF OFF-GRID MICROGRID
The statistical evaluation of the algorithms, then the use of the evaluation metrics have proved that the hHHO-AOA is truly superior compared to other two algorithms. Therefore, we are motivated to form three microgrids based on the VOLUME 10, 2022  component capacities determined by the sizing optimization of each algorithm. In brief, this section presents the simulation of each microgrid and examines the performance of each design. In Table 9, the economic and reliability outputs of the system and the capacities of the components are given. When this table is examined in terms of component capacities, the proposed microgrid with hHHO-AOA has the highest PV, WT, and BESS capacity. This reflects positively on the reliability and utilization rate of renewable energy sources. With 6.51% LPSP and with 77.69% RF, the hHHO-AOA has produced a microgrid with the highest reliability and the highest clean energy utilization rate. The microgrid designed with the hHHO-AOA is more reliable with an LPSP value of 14.31% less than HHO and 12.51% less than AOA. In terms of the use of renewable resources, the hHHO-AOA algorithm uses renewable energy resources at a rate of 1.71% higher than HHO and 3.31% higher than AOA. When evaluated in terms of COE, the energy cost of the microgrid designed by the hHHO-AOA is 4.08% more expensive than HHO and 1.90% more expensive than AOA. There is a slight increase in the COE. However, we believe that 4.08% higher COE is justifiable compared to 14.31% lower LPSP. Creating a design that is highly reliable and producing results that are closer to the clean energy objectives is perfectly acceptable against the slight increase in COE. In brief, the microgrid design with the hHHO-AOA is considered satisfactory and successful in meeting the overall expectations.
The contribution of the component's capacities (that is, the distribution of the energy produced in the microgrid by the resources) to the energy production is given in Fig. 10. In the microgrid designed with the hHHO-AOA algorithm, 58% of energy is produced from PV, 24% from WT, and 18% from DG. In the design with the HHO algorithm, the energy production from PV is increased by 1% compared to the hybrid algorithm, and it becomes 59%. Accordingly, it is 22% from WT and 19% from DGs. In the design with the AOA algorithm, the energy production from PV is also 59%, same as the HHO case, 21% from WT, and 20% from DG. When Fig. 10 is evaluated in general, it can be seen that the highest rate of energy production from WT and the lowest DG usage are in the hHHO-AOA algorithm. This outcome supports the ecofriendly microgrid design and it is favorable.
The cost of the components that make up the individual microgrids are given in Table 10. According to Table 10, the largest expense item in the total cost is the investment costs. Therefore, the emerging investment costs of the hHHO-AOA, the HHO, and the AOA algorithms within the total cost becomes 58.61%, 56.53%, and 56.45%, respectively. Moreover, each microgrid infrastructure requires renewal of BESS and DG components throughout the project life. We can also see in Table 10 that the rate of the renewal costs within the total cost, according to above algorithms order, is 58.61%, 56.53% and 56.45%, respectively. Based the same order, the share allocated to the renewal of the DGs is 7.48% in the hHHO-AOA, 8.62% in the HHO, and 8.17% in the AOA. Here, it is seen that the hHHO-AOA algorithm has yielded the lowest DG replacement cost. This situation is related to the annual operating time of the DG, which directly affects the service life. When the annual usage rate, which is reflected in the usage percentage of the DG, is examined in terms of hours, the DG is used for 2621 hours per year in the hHHO-AOA algorithm, 2810 hours in the HHO algorithm, and 2807 hours in the AOA algorithm. From this, it can be concluded that the hybrid algorithm used 6.72% less DG than HHO and 6.62% less than AOA. It is an important advantage that the microgrid designed by the hHHO-AOA algorithm gives more priority on renewable energy sources and lowers the DG usage. Hence, the operating and maintenance costs within the total cost have a share of 3.96%, 3.82%, and 3.73% in the hHHO-AOA, the HHO, and the AOA algorithms, respectively. In the microgrid designed with the hHHO-AOA algorithm, the PV and BESS capacity, especially the WT, is larger than the microgrid components designed by the other two algorithms. The part of this issue that affects the operating and maintenance costs of these components and the power electronics converters that provides the required interface is reflected in Table 10. In addition to the operating and maintenance costs, the fuel cost of the DG, which depends on the usage rate, is also examined. The shares of the fuel cost within the total costs of the hHHO-AOA, HHO, and AOA algorithms are 29.48%, 31.92%, and 31.83%, respectively. The less DG usage in the hHHO-AOA algorithm means less fuel cost compared to other algorithms. The salvage cost is calculated as an income at the end of the project life, in addition to all costs incurred as expenses during the life of the project. The salvage values of the components are given in the Table 10. The 0.31% of the total cost in the hHHO-AOA, 0.17% in the HHO and the AOA algorithms are recovered as the salvage cost of the DG.
Two particular days in the year are chosen to examine the energy production and consumption of the microgrids formed by the algorithms. These days are the January 1, 2018 and August 23, 2018. Fig. 11, which is divided into five sections, shows the data for the day, January 1, 2018 and includes PV power, WT power, BESS's state of charge (SOC), and DG power in the listed order. In the fifth section of Fig. 11, we see the total power that is produced from renewable energy sources is subtracted from the load power. The negative part of the axis, that is, the part below the zero axis, is the surplus energy to be stored in the BESS. The positive part above the zero axis represents the energy demand that cannot be met from renewable energy sources. This demand therefore is supplied by the BESS and DG. Due to the winter season effect, the highest PV energy produced on this particular day has been realized as 239.19 kW, 237.66 kW, and 233.72 kW in spite of PV capacities that were determined as 499.22 kW by the hHHO-AOA, 496.01 kW by the HHO, and 487.80 kW by the AOA, respectively. While WT had a maximum generation capacity of 142.63 kW, 130.93 kW, and 125.63 kW according to the hHHO-AOA, HHO, and AOA algorithms, it has produced a maximum of 140.35 kW, 128.84 kW, and 123.62 kW of energy on the selected day, respectively. These results suggest that the energy produced from the PV system is decreased by approximately 52.08% due to the winter conditions; on the contrary, WTs are utilized close to their maximum levels.
Similarly, the status of energy production and consumption that took place on August 23, 2018 is given in Fig. 12. The major difference on this day compared to the selected winter day is that the BESS is charged to its maximum capacity (3rd graph in Fig. 12). However, this was only 17.69 kWh for the hHHO-AOA, 6.41 kWh for the HHO, and 1.67 kWh for the AOA on the selected winter day. As before, the positive part of the zero axis represents the energy demand that could not be met from the renewable energy sources. This demand is therefore met mostly from the BESS plus the DGs.
Due to the advantage of the summer season, PVs have produced close to their capacities. The maximum production capacity of PVs have been determined as 499. 22  Even though there is no definitive method to measure the degradation of batteries, it must be considered whenever the cost and lifetime of a system are evaluated. The degradation of batteries depends on many factors, but it is commonly evaluated in two categories. The first category is the calendar life, which is independent of battery operation. The second category is called the cycle life, in contrast, it is heavily dependent on the battery operation. The cycle life is affected by the operating temperature, the discharging rate, the charging rate, and the depth of discharge (DOD). According to reference [69], total degradation of battery life through, in other words the degradation of cycle-life is calculated by equation (64), where n c is the number of cycles, R j is the amplitude of micro-cycles within the DOD history, A and B are the empirical constants provided by battery manufacturers.
Based on our literature search, we have discovered that the rainflow counting algorithm can be used to count the number of cycles (n c ) and to obtain the R j parameter for a given DOD data [70]. It is first proposed by Matsuishi and Endo in 1968 as a counting technique, but it has evolved to a well-known cycle counting algorithm over the years [71]. Therefore, the DOD data belonging to the BESS component of our microgrid design, which are collected on hourly basis for over a year as shown in Fig. 13, are used as the input to the rainflow analysis. Based on the post processing performed by the rainflow algorithm, the number cycles that are subjected to the battery degradation is determined as 510 along with the R j data as shown in Fig. 14, which gives the amplitude of DOD versus the total cycle count.
In this microgrid design study, we have decided to evaluate battery degradation for three battery types: Lithium-Titanium, Lithium-LiFePO4, and Lithium-NMC. The A and B constants for each battery type are obtained from the datasheets as 16813 and -1.887 for Lithium-Titanium,  6372 and −1.433 for Lithium-LiFePO4, 3351 and −1.689 for Lithium-NMC [69]. Using all in equation (64), the annual degradation percentages of cycle-life for each battery is calculated as 1.1406% for Lithium-Titanium, 3.3198% for Lithium-LiFePO4, and 5.2098% for Lithium-NMC. Based on annual degradation percentages of cycle-life, the projected cycle-life degradation for 10-years is given in Fig. 15. As shown, at the end of 10 years, the usable capacity stands at 89.1617% for Lithium-Titanium, 71.3466% for Lithium-LiFePO4, and 58.5647% for Lithium-NMC.
Based on above analysis, we can conclude that degradation of batteries directly affects the investment, operation and replacement costs in a microgrid design. According to Fig. 15, when Lithium-NMC batteries with a low cycle-life are used to solve the optimization problem, the renewal period should be taken between 5-6 years. Here we assume that a battery can be used down to around 70% DOD before it is considered unusable. This number is typically 80% in electrical vehicles. Similarly, when Lithium-LiFePO4 batteries are used with the same system components, this period can be taken as 9-10 years. When Lithium-Titanium batteries with a high life cycle are preferred, then the renewal period becomes over 10 years. As it is seen, the life versus cost is a tradeoff. Therefore, the user's decision on the battery type becomes an important parameter in a microgrid optimization problem. Since Lithium-LiFePO4 is commonly used battery type in the industry with decent cost versus life performance compared to other types, we have considered using Lithium-LiFePO4 in our design and taken the life time of the battery as 10 years.

VII. SENSITIVITY ANALYSIS
After the evaluation of the microgrid simulation results, the performance of the hHHO-AOA algorithm in sizing optimization and optimal microgrid design is verified. From this point of view, in this section, a sensitivity analysis over the effect of changing of the capacities of the components on COE, LPSP and RF is made. The effect of the change of PV-WT, PV-BESS, PV-DG, WT-BESS, WT-DG and DG-BESS capacities on COE is given in Fig. 16. In the first graph, the effect of the change of PV and WT capacities on COE is given. The lowest COE is obtained for the WT power range of 100-200 kW against the maximum PV power. The change of PV and BESS capacities is given in the second graph. The most effective energy generation is realized with a PV power of over 400 kW and a BESS capacity from 200 kWh to 700 kWh. In the third graph, the change of PV and DG capacities is given. A PV power of more than approximately 350 kW and a maximum DG capacity of 185 kW result in an optimal energy production. Other graphs can be interpreted with similar point of view. In the microgrid designed with the hHHO-AOA algorithm, the COE is obtained as 0.2090 $/kWh with 499.22 kW PV, 142.63 kW WT, 644.04 kWh BESS, and 100 kW DG components. It is seen from the graphs that the most suitable component capacities have come together for the designed microgrid to obtain the most profitable COE.
In addition to the COE, the effect of changes in PV-WT, PV-BESS, PV-DG, WT-BESS, WT-DG, and DG-BESS capacities on the LPSP is given in Fig. 17. It is obvious that component capacities should be chosen approximately at  their maximum values for the lowest LPSP, which yields the most reliable microgrid design. However, this will result in the highest COE and the highest investment cost. Therefore, the LPSP alone is not a sufficient criterion. A compromise has to be made. It is this moment that optimization plays an important role in determining the most appropriate value for the LPSP and for the best trade off. In the microgrid designed with hHHO-AOA algorithm, the LPSP is obtained as 6.5064% with 499.22 kW PV, 142.63 kW WT, 644.04 kWh BESS, and 100 kW DG components. Based on these graphs, it is seen that the most suitable component capacities have come together to attain the best compromise, the most reasonable LPSP along with the most suitable COE.
In addition to the COE and the LPSP, the effect of changes in PV-BESS, PV-DG, WT-BESS, WT-DG, and DG-BESS capacities on the RF is given in Fig. 18. The choice of PV and WT capacities is the decisive factor of how well we can utilize the renewable energy sources. Higher the contribution of PVs and WTs in a microgrid as the energy producing components, the higher the RF will be. In the microgrid designed with the hHHO-AOA algorithm, an RF of 77.6879% is obtained with 499.22 kW PV, 142.63 kW WT, 644.04 kWh BESS, and 100 kW DG capacities. As a result of the sensitivity analysis, with the component capacities determined for the designed microgrid, the most suitable COE, the most reasonable LPSP, and the highest RF have been obtained.

VIII. CONCLUSION
In this paper, we have presented the design and development of a new hybrid metaheuristic algorithm, named as the hybrid Harris Hawks Optimizer-Arithmetic Optimization Algorithm (hHHO-AOA). The hybrid algorithm is expected to increase solution diversity and therefore the accuracy in sizing optimization of microgrid designs. The performance of the hHHO-AOA is assessed with the evaluation metrics and the statistical tests. According to the Friedman ranking test and the Wilcoxon signed-rank test, the hHHO-AOA has shown significant improvement in the solution accuracy of the sizing problem compared to the HHO and the AOA algorithms. Later, the developed hybrid algorithm is tested on a microgrid that consists of a PV power system, a WT power system, a BESS, DGs, and a commercial-type load. For the optimal sizing of these components, a problem that defines the LPSP and the COE as the objective function is formulated. The optimization done by the proposed hybrid algorithm has produced the lowest LPSP value (6.5064%) and the lowest COE (0.2090 $/kWh) along with the highest rate of RF (77.6879%). For consistency, we have taken the battery degradation into account in this design. The battery degradation is evaluated based on percent degradation of cycle-life using the rainflow algorithm. In conclusion, it is demonstrated that the microgrid that is designed by the hHHO-AOA is reliable, economical, and eco-friendly compared to other designs.
Our work will continue on autonomous microgrids, where local energy potentials are turned into advantages with renewable energy sources. In particular, our plan is to work on optimal operation of autonomous microgrids where the loads are considered in the complex category such as the commercialtype loads together with electric vehicles.