Fitness Distance Balance Based LSHADE Algorithm for Energy Hub Economic Dispatch Problem

This paper presents an improved version of Linear Population Size Reduction Success History Based Adaptive Differential Evolution (LSHADE) algorithm for solving global optimization problems. Fitness Distance Balance (FDB) selection method was used to redesign the mutation operator in order to enhance the search performance of the LSHADE algorithm. In order to test and validate the performance of the proposed algorithm, a comprehensive experimental study was carried out. For this purpose, it was tested on the CEC14 and CEC17 benchmark problems, consisting of different problem types and dimensions. Results of the FDB-LSHADE was compared to the performance of 8 other up-to-date and highly preferred metaheuristic search (MHS) algorithms. According to Friedman test results, the proposed FDBLSHADE algorithm ranked first among the all competing algorithms. Moreover, the proposed algorithm was used to solve single- and multi-objective energy hub economic dispatch (EHED) problems, which were a non-convex, a nonlinear, and high dimensional problems. To analyze the results of the proposed algorithm obtained from experimental studies, two non-parametric statistical methods, which are Wilcoxon and Friedman tests, were used. The simulation results of the proposed algorithm were compared to the results of the 8 other MHS algorithms. The results demonstrated that the FDB-LSHADE was a superior performance compared to other MHS algorithms for solving both benchmark and EHED problems.


I. INTRODUCTION
Meta-heuristic search (MHS) algorithms have been widely preferred in order to solve various optimization problems from past to present. In general, all of MHS algorithms have the same features from the point of search behavior: exploration and exploitation. In exploration phase, algorithms perform a random expansion search in the search space to increase the diversity of the solutions. On the contrary, exploitation performs local search in the neighborhood of a reference location in the search space for improving the solution quality. Therefore, the balance between exploration and exploitation is important so that the optimization algorithms do not get stuck in the local minimum.
In the literature, hundreds of MHS algorithms have been presented in the literature and they can be classified as four The associate editor coordinating the review of this manuscript and approving it for publication was Sun-Yuan Hsieh . categories [1]: (i) swarm intelligence [2]- [5], (ii) evolutionary [6], [7], (iii) human knowledge based algorithms [8], [9], and (iv) physics-based algorithms [10]- [12]. Differential Evolution (DE) algorithm presented by Storn and Price [6] is one of the most used and successful algorithm among evolutionary algorithms for the solving mathematical optimization problems. The DE algorithm mainly consists of four main steps: initialization, mutation, crossover, and selection. All three steps, mutation, crossover, and selection, play a vital role in the DE algorithm; however, the DE guides the population evolution through mutation operations based on difference vectors. Researchers have improved the performance of the DE algorithm by using different mutation strategies.
Various variants of the DE algorithms have been proposed until now. Among these variations, the Linear Population Size Reduction Success History Based Adaptive DE (LSHADE), which was the winner of the CEC14 Competition on Real-Parameter Single Objective optimization [13]. As a result of our research, the performances of the LSHADE algorithm and its variants were generally tested on CEC benchmark problems. However, we noticed that there are no studies on constrained engineering optimization problems with the LSHADE algorithm. The aim of this study was to diversify the search process and to improve the convergence performance of the LSHADE algorithm. For this purpose, a comprehensive study was carried out to redesign the mutation operator that guides the search process. Our motivation to achieve this was the Fitness Distance Balance (FDB) selection method proposed by Kahraman et al. in 2020 [14]. The FDB method was applied for the first time to improve the exploration and exploitation capabilities of the symbiotic organism algorithm (SOS). In another study, the FDB method was used to strengthen the diversity and balanced search capabilities of the SFS algorithm [45]. The FDB method was used to redesign the three candidate solutions used in the mutation strategy that guides the search process to improve the exploration capability of the AGDE algorithm, which is one of the versions of the DE algorithm [46].
In this study, we tried to improve the exploration ability of the LSHADE algorithm using the FDB method. Consequently, this redesigned algorithm was called FDB-LSHADE, which is an improved version of the LSHADE algorithm, has been developed. In the study, six variants of the FDB-LSHADE were proposed. In order to test and verify the performance of the proposed FDB-LSHADE variants, a comprehensive experimental study was carried out. For this purpose, CEC14 [15] and CEC17 [16] benchmark suites were used, where the benchmark test suites include different types of optimization problems from simple to complex. In order to show the performance of the proposed algorithm in this study, 8 recently up-to-date and often preferred MHS algorithms for optimization problems were used, consisting of moth-flame optimization algorithm (MFO) [17], salp swarm algorithm (SSA) [18], adaptive guided differential evolution algorithm (AGDE) [19], supply-demand-based optimization (SDO) [20], artificial electric field algorithm (AEFA) [21], artificial ecosystem-based optimization (AEO) [22], equilibrium optimizer (EO) [23], and marine predators algorithm (MPA) [24]. To analyzed data obtained from experimental studies, two non-parametric test methods, which are Wilcoxon and Friedman test, were used. The results of the analyses showed that FDB-LSHADE algorithm developed in this study is the most effective MHS algorithms compared to its competitors in two different benchmark suites.
In order to test and verify the performance of the proposed FDB-LSHADE on unconstrained optimization problems, the studies are explained above. Furthermore, a comprehensive experimental study has been conducted to test and validate the performance of the FDB-LSHADE on constrained real-world problems. Therefore, we chose the Energy Hub Economic Dispatch (EHED) problems for this. The following paragraph provides the information about the studies for this topic.
Energy Hub Economic Dispatch problem is based on the Economic Dispatch (ED) and Energy Hub (EH) concept.
ED is one of the most important and popular optimization problem in modern power systems. The goal of ED problem is to optimal planning of the power generation inputs in order to meet the load demand subject to equality and inequality constraints. For this reason, ED problem is a non-convex and nonlinear constrained optimization problem. The most used ED problem is the minimization of the total generation cost of electric power generation. In the past decade, heuristic optimization techniques have been used for solving ED problems [25]- [27].
Energy hub is defined as an integrated energy systems which includes the energy generation, conversion, and storage systems. That is, energy hub receives the input energies (i.e. electricity, natural gas, heat, and so on) at the input port and generates or converts them into another energy form by using technologies such as transformers, inverters, converters, heater exchangers, gas furnace, combined heat and power (CHP) unit, combined heat, cooling, and power unit (CHCP) unit, and etc. [28]- [30].
Many studies have been carried out about energy hub such as modelling, optimal power flow, and economic dispatch. In the study on energy hub modeling, Geidl and Andersson proposed a general matrix model for an energy hub system including various energy carriers such as electricity, natural gas, and heating. The aim of the study was to minimized the total cost of the system [31]. The same authors proposed an EH model that includes multiple energy carriers (MECs) such as electricity, natural gas, and heat, in order to solve the OPF problem. The objective was to minimize the total generation cost [32].
In the literature, MHS algorithms were used for solving energy hub optimization problems. Moeini-Aghtaie et al. focused on a method for decomposing the combined power flow studies with MECs into the classic optimal power flow (OPF) problem while retaining the significant benefits of simultaneous MEC analysis via the multi-agent genetic algorithm (MAGA). The goal was to keep the system's total energy cost as low as possible. The total costs of the system with and without energy hubs were calculated using the MAGA method [33]. In another study, Moeini-Aghtaie et al. used the MAGA method in order to solve a bi-level ED problem in energy hub in the presence of wind uncertainty. The inputs of the energy hub were electricity, natural gas, and wind energy. The presented model applied to 11-hub test system in a 24-h period. The results of MAGA showed that it has a better convergence behavior compared algorithm PSO and GA [34].
In the literature, the modified teaching-learning-based optimization algorithm (MTLBO) has been used to solve optimization problems in many studies. Shabanpour-Haghighi and Seifi proposed MLTBO in order to solve the OPF problem in MECs system. The problem was solved by the GA, PSO, and original TLBO and their results were compared to the MTLBO method. The results showed that the MTLBO algorithm provided better solutions than the other competing algorithms used in the study [35]. The same authors studied optimal operation of MEC system using the MTLBO algorithm. The aim of the study was to minimize the total operation cost and total emission of the system. The results obtained from the MTLBO algorithm were compared to the results of other techniques in the literature [36]. In another study using MLTBO algorithm, the same authors used it in order to solve the OPF problem for MEC system with an electrical, a natural gas, and a district heating subnetwork [37]. Shabanpour-Haghighi and Seifi used MTLBO algorithm in order to solve 24-h OPF problem of multi-career energy networks. The energy hubs used in the study have three inputs (i.e., electricity, natural gas, and district heating) and two outputs (i.e., electricity and heat). In order to show the superiority of the proposed algorithm, the results are compared with the other MHS algorithms [38].
In many studies presented in the literature, the time varying acceleration coefficient method (TVAC) has been used as a hybrid with MHS algorithms. In one of these studies, Beigvand et al. proposed time varying acceleration coefficient particle swarm optimization (TVAC-PSO) algorithm in order to solve Multiple Energy Carriers Economic Dispatch (MECED) problems. The results showed that the suggested TVAC-PSO algorithm reached a better convergence performance than the other presented algorithms such as GA, PSO, and DE [39]. The same authors introduced time varying acceleration coefficient gravitational search algorithm (TVAC-GSA) in order to solve Multiple Energy Carriers Optimal Power Flow (MECOPF) problems. The aim of the study was to minimize the energy cost and electrical transmission loss of electricity-gas network. In order to test and verify the TVAC-GSA algorithm, the results of the presented method were compared to the results obtained by GSA, PSO, GA, and DE. It showed that the proposed algorithm had better performance and fast convergence to solve the presented problems [40]. By developing TVAC-GSA algorithm, Beigvand et al. proposed self-adaptive learning with time varying acceleration coefficient-gravitational search algorithm (SAL-TVAC-GSA) in order to solve to solve single-and multi-objective Energy Hub Economic Dispatch (EHED) problems. They modelled the EHED problems and tried to determine the optimal operation strategy that minimizes energy cost and losses. For this motivation, a highcomplex, large-scale energy hub system consisting of 29 hub structures where electricity, gas, and heat were consumed and electricity, heat, cooling, and compressed air were produced. The results obtained by SAL-TVAC-GSA were compared to TVAC-GSA, enhanced GSA (EGSA), GSA, PSO, and GA in terms of quality solution and computational performance. It was seen that the performance of presented algorithm found better optimum solution than its competitors [41].
In this study, we aimed to solve the EHED problems using the LSHADE and FDB-LSHADE algorithms. Three objective function including total energy cost, total energy hub losses, and combination of energy cost and energy hub losses as a multi-objective problem were to be minimized within the equality and inequality constraints . The simulation results  obtained from the FDB-LSHADE for these problems were  compared to the results of the MHS algorithms, which were  GA, PSO, GSA, EGSA, TVAC-GSA, and SAL-TVAC-GSA, given in [41]. Moreover, the results of proposed algorithm were compared to 8 up-to-date MHS algorithms consisting of LSHADE, MFO, SSA, AGDE, SDO, AEO, EO, and MPA. The simulation results proved that FDB-LSHADE algorithm showed superior performance compared to the competing algorithms for the solving EHED problems.
The main contributions of the present study to the literature are as follows: • A powerful MHS algorithm called FDB-LSHADE was developed. The FDB selection method was used to redesign the mutation operator used in LSHADE algorithm in order to strengthen the exploration ability and balance search capabilities of the base LSHADE algorithm.
• The proposed FDB-LSHADE algorithm is potentially presented to the literature as one of the most powerful MHS algorithm. This paper contributes to the literature as one of the most comprehensive studies carried out to test and validate the performance of the proposed FDB-LSHADE algorithm using CEC14 and CEC17 benchmark suites. The results of the FDB-LSHADE algorithm were compared to 8 competing algorithms, which are up-to-date and most preferred MHS algorithms. In order to evaluate the results obtained from the experimental studies, statistical analysis methods were used. The results of them demonstrated that the proposed algorithm was shown a superior performance among its competitors.
• The proposed FDB-LSHADE algorithm was used for solving EHED problems, which is a non-convex, nonlinear, and high-dimensional constrained optimization problem. Three objective functions were used to evaluate the performance of the proposed algorithm. The superiority of the FDB-LSHADE algorithm has been proven, when compared with the results of other competing MHS algorithms and the results previously reported in the literature.
• The performance of the FDB-LSHADE algorithm has been verified both in benchmark problems and in a realword constrained engineering problem. Thus, we can say that the proposed method is a powerful algorithm and it could be used in different engineering problems by researchers. The remainder of this article is organized as follows. Section 2 gives the mathematical modelling of energy hubs and the formulation of the EHED problems. Section 3 consists of three sub-section introducing the LSHADE algorithm, the FDB selection method, and the proposed FDB-LSHADE algorithm. Section 4 presents the experimental study settings and the standards taken into account in conducting the experimental studies. Section 5 presents the performance of the algorithms and the statistical test results of the experimental studies. This section consists of three sub-sections. In the first sub-section, the performance of the variants of the FDB-LSHADE algorithm on the CEC14 and CEC17 benchmark problems was evaluated. The second sub-section presents the experimental and statistical analysis of the proposed FDB-LSHADE algorithm and the other MHS algorithms on benchmark problems. In the third sub-section, the EHED problems were solved by the FDB-LSHADE algorithms and competing MHS algorithms. In all sub-sections, statistical analysis methods were used to evaluate the data obtained from MHS algorithms. Finally, the conclusions are presented in Section 6.

II. FORMULATION OF THE PROBLEM
In this section, the energy hub model is explained in detail and the related equations of the energy hub structures given in the literature are listed. Moreover, the cost model of the energy sources is given. Then, the objective functions and constraints are presented.

A. MATHEMATICAL MODELING OF ENERGY HUB STRUCTURES
Energy hub is generally defined as an interface between different energy infrastructures and loads. It can be also defined as an integrated energy systems which includes energy generation, conversion, and storage systems. In other words, different forms of energy (i.e., electricity, natural gas, heat, etc.) consume at the input ports of the energy hub and supplies required energy services (i.e., electricity, heat, cooling, and compressed air, etc.) at the output ports for costumers [28]- [30], [42]. In the energy hub, using, i.e., transformers, CHP technology, heater exchangers, and other equipment, energy carriers are converted and conditioned. This equipment can be classified as: direct connections, convertors, and storage elements. With the direct connections, the energy carriers at the input ports are transmitted to the output ports without being converted to another form. By converter elements, energy is transformed into different forms. Storage devices are used to store the various types of energy and each type of energy is stored using a different technology.
The number of energy carriers at the input and output ports of the energy hub may vary. Based on this, energy conversions in the energy hub can be categorized as: single input and single output, single input and multiple outputs, multiple inputs and single output, and multiple inputs and multiple outputs. When the energy hub model has more than one input or output, a coupling matrix is defined. A single input and single output energy hub model can be described by Eq. (1). Here, E in α and E out β represents the input and output energy, respectively. C αβ describes the coupling factor that is the efficiency of the conversion element, and the subscript α and β describe the input and output of the energy carrier.
A multiple input and multiple output energy hub model can be expressed as: (2) where E in and E out are the input and output energy vector, respectively. C is the coupling matrix, the subscripts {α, β, . . . , ω} represent the energy carriers, and the elements of coupling matrix are the coupling factors. In the coupling matrix, each coupling factor is relevant to a specific input and output. In single input single output systems, the coupling factor is equal to the conversion element's efficiency; nevertheless, in multiple inputs and multiple outputs systems, the coupling factor is unequal to the conversion element's efficiency. An energy carrier at the input port can be split into numerous parts and used as a converter's input, or the output of the energy hub can be used as a conversion element's input.
In this case, a dispatch factor, describes how much energy will be distributed to each conversion element, is defined for the input or output energy. An special energy hub model created according to the above information is given in the Fig.1. In Fig. 1., the input energy carriers are electricity (e), natural gas (g), and heat (h) and the demands are electricity (e), heat (h), cooling (c), and compressed air (a). Moreover, two dispatch factor are used, in which v 1 is used for input energy and v 2 is used for output energy. The conversion elements used in the Fig. 1. are explained as below: • Transformer (T): It employs and supplies electricity both at the input and output port, respectively.
• Combined heat and power (CHP): This device produces heat by burning natural gas.
• Combined heat, cool, and power (CHCP): It transforms the natural gas into electricity, heat, and.
• Gas furnace (GF): This device produces heat by burning natural gas.
• Compressor (C): This device uses electricity to generate compressed air.
• Heater exchanger (HE): It consumes and provides heat at its input and output, respectively. In [41], 29 energy hub structures are given. The elements used in these structures, transformer, CHP, CHCP, gas furnace, and heater exchanger. Moreover, electricity, natural gas, and heat are used as input energy carriers. The energy hub structures used in this study are given in Fig. (2). The energy conversion expressions of these structures are given in Eqns. (3) - (31).
EH #1: The hub contains only transformer unit and is modelled as: where η T is the efficiency of transformer.  EH #2: It includes only CHP unit and is formulated as follows: E out where η CHP e and η CHP h are the CHP' s efficiencies related to electricity and heat. EH #3: Natural gas is converted into electricity, heat, and compressed air by CHCP unit as follows: where η CHCP e , η CHCP h , and η CHCP c are the CHCP's efficiencies related to electricity, heat, and cooling. EH #4: Only heat is generated by gas furnace employing natural gas as: where η GF is the efficiency of gas furnace.
EH #5: This hub consists of only one heater exchanger and is formulated as: where η HE is the efficiency of heater exchanger.
EH #6: Electricity and heat can be provided using electricity and heat through transformer and CHP as follows: EH #7: It produces electricity, heat, and compressed air employing only electricity as below: EH #8: Transformer and CHCP construct this hub and the energy conversion of it as follows: EH #9: This hub consists of CHP and gas furnace and can be modeled as: EH #10: In this hub, electricity and heat are produced by employing electricity, and heat as: EH #11: Two forms of energy (i.e., electricity and heat) are consumed by CHCP and heater exchanger in order to provide electricity, heat, and cooling as follows: EH #12: It includes CHCP and gas furnace, and produces electricity, heat, and cooling by consuming only electricity through the following formula: EH #13: Transformer and heater exchanger constructs this hub and can be formulated as: EH #14: Transformer and gas furnace convert electricity and natural gas into electricity and heat by the following expression: EH #15: Natural gas and heat is employed by gas furnace and heater exchanger, respectively and the energy conversion of this hub can be expressed as below: EH #16: This hub supplies the electricity, heat, cooling, and compressed air for consumers by employing only natural gas as below: (18) 66774 VOLUME 10, 2022 EH #17: Natural gas is used to produce electricity, heat, and compressed air through CHCP and compressor as follows: (19) EH #18: Three elements including transformer, CHCP, and heater exchanger is used in the structure of the hub in order to produce electricity, heat, and compressed air. The mathematical expression of the hub is as follows: EH #19: The hub has the same structure as the previous hub; however, CHCP is replaced with CHP in this hub. For customers, electricity and heat are generated and the energy conversion can be stated as: EH #20: In this hub, three types of energy are produced consuming electricity and heat by transformer, CHCP, and gas furnace as follows: EH #21: It consists of transformer, CHCP, and gas furnace in order to provide electricity, heat, and cooling via the following mathematical expression: EH #22: The hub structure is formed by transformer, CHCP, and compressor in order to supply electricity, heat, cooling, and compressed air as: EH #23: It is similar to previous hub structure, but CHCP is replaced with CHP. In this hub, electricity, heat, and compressed air are provided by the following mathematical formula: EH #24: In this hub, electricity and natural gas are converted into electricity, heat, cooling, and compressed air via transformer, CHCP, compressor, and gas furnace as (26), shown at the bottom of the page.
EH #25: This hub produces electricity, heat, and cooling from electricity and natural gas as (27), shown at the bottom of the page.
EH #26: Three forms of energy (i.e., electricity, heat, cooling, and compressed air) can be provided by using electricity, natural gas, and heat. The mathematical expression of the hub is as below: EH #27: In this hub, by employing electricity, natural gas, and heat, electricity, heat, and compressed air are generated via transformer, CHP, compressor, and VOLUME 10, 2022 FIGURE 2. The hub structures given in [41].
heater exchanger as follows: EH #28: It includes five types of converters which are transformer, CHCP, compressor, gas furnace, and heater exchanger and the energy conversion of it is formulated as (30), shown at the bottom of the next page.
EH #29: The structure of the hub is similar to previous hub; however, CHCP is replaced with CHP in this hub. The mathematical formulation of it can be stated as (31), shown at the bottom of the next page.

B. COST MODELS OF ENERGY SOURCES
The energy hub models consist of three input energy carriers (i.e., electricity, natural gas, and heat). The cost model of each energy carrier is calculated separately.

1) COST MODEL OF THERMAL GENERATING UNITS
In thermal generating units, the classical generation cost function is defined as a quadratic cost function in Eq. (3) that depends on the generated output active power [28], [41]. Here, a j , b j , and c j represent the cost coefficients of the j th thermal generating unit, P j shows the output power of j th thermal generating unit, and N is the total number of thermal generating units. In this study, the valve point effect has been taken into account when calculating the generation cost of the electrical energy carrier. The cost of the electrical energy carrier is calculated by Eq. (33). Here, d j,e and e j,e are the cost coefficients of the valve-point effect for electrical carrier of the j th source, E in j,e and E in,min j,e are respectively the energy production of j th source and the minimum value of energy production of j th source, CF e denotes the total energy cost of the electrical energy carrier, and N e is the total number of electrical energy carriers.

2) COST MODEL OF NATURAL GAS AND HEAT ENERGY CARRIERS
The energy production cost of natural gas and heat is calculated according to the classical generation cost of the thermal generating units given in Eq. (32) [41]. The cost models of the natural gas and heat are given in Eq. (34) and Eq. (35), respectively. Here, CF g and CF h represent the total energy cost of the natural gas and heat energy carrier, respectively, are the energy production of j th natural gas source and the energy production of j th heat source, respectively. N g and N h are the total number of the natural gas and heat energy sources.

C. OBJECTIVE FUNCTIONS
Objective functions evaluated in this study are minimization of total EH cost, minimization of total EH loss, and minimization of total EH cost and loss.

1) MINIMIZATION OF TOTAL ENERGY COST
The aim of this problem is to minimize the total energy cost of energy hub while satisfying various constraints [41]. The cost of all energy carriers (i.e., electricity, wind energy, solar energy, natural gas, and heat) used at the input port of energy hub structures is considered. Their costs are calculated separately and the total cost is the sum of the costs of all energy carriers given in Eq. (36)(37)(38)(39)(40). The cost function is: where OF 1 represents the objective function and N hub is the total number of energy hub structures.

2) MINIMIZATION OF TOTAL ENERGY HUB LOSS
The loss is defined as the difference between the input and output energy of a system. The objective is to minimize the energy losses for whole energy carriers in all energy hubs and it can be mathematically expressed as [41]: where OF 2 represents the objective function, E out m,n represents the output energy of the energy hubs for output energy carrier n.

3) MINIMIZATION OF TOTAL ENERGY COST AND HUB LOSS
The objective of this problem is minimizing the total energy hub cost and loss, simultaneously. It is a multi-objective problem and is converted into a single-objective problem given in Eq. (38) [41]. Here, OF 3 represents the objective function, E demand n is the demand value for output energy carrier n.

D. CONSTRAINTS OF THE PROBLEM 1) EQUALITY CONSTRAINTS
In this study, equality constraints of the objective function is described as the balance between the generated energy VOLUME 10, 2022 at the output of energy hubs and the demands. It can be mathematically explained as [41]:

2) INEQUALITY CONSTRAINTS
The inequality constraints of the problem can be classified as two sub-sections.

a: OPERATING CONSTRAINTS OF ENERGY HUBS
The minimum and maximum limits of the all energy carriers are expressed as follows [41]: All dispatch factors used for energy hubs should satisfy the following constraint [41]:

III. METHOD
This section is divided into three sub-sections to increase the clarity of the proposed FDB-LSHADE algorithm used in this study. First sub-section presents the explanation of the LSHADE algorithm. In the second sub-section, the FDB selection method is explained. In the last sub-section, the proposed FDB-LSHADE algorithm is explained.

A. OVERVIEW OF LINEAR POPULATION SIZE REDUCTION SUCCESS HISTORY BASED ADAPTIVE DE (LSHADE)
LSHADE is the extend version of SHADE [43] with Linear Population Size Reduction (LPSR) strategy. In LPSR, the population size is decreased continuously according to a linear function. Similar to the other evolutionary algorithms, LSHADE begins with a initial population P 0 . Each j th component (j = 1, 2, . . . , D) of the i th individuals (i = 1, 2, . . . , NP) in the P 0 is defined as [2,3]: where rand (0, 1) is the a uniformly distributed random number between [0,1)].
In LSHADE, a historical memory with H entries are maintained with control parameters F, CR, M F , M CR . F ∈ [0, 1] is scaling factor controlling the magnitude of the mutation operator and CR ∈ [0, 1] is the crossover rate. At the beginning, M CR,k and M F,k (k = 1, . . . , H ) are set to 0.5 and each generation G, F i and CR i for each x i are generated by randomly selecting an index r i from [1, H]. Then, F i and CR i are determined by Eq. (43) and (44), respectively. If the generated F i value is out of its boundary values, Eq. (43) is repeatedly until the value are produced within the limits. When generated CR i values is out of [0, 1], it is replaced with the limit value (0 or 1) which is closest the generated value [13].
In the process following the setting of the control parameters, a mutant vector v G i is generated according to Eq. (45) for each target vector x G i at generation G. For mutation process, DE/current-to-pbest/1/bin mutation strategy is used, which was proposed in JADE [44]. Here, r 1 and r 2 are randomly chosen indices from the population, x G pbest is the best individual vector which has the best fitness value in the population at generation G. F G i is the mutation scale factor which controls the magnitude of the difference vector [13].
For each dimension j, if the mutant vector element v G i is out of the upper and lower boundaries x min j , x max j , the correction given in Eq. (46) After generating the mutant vector v G i is crossed with the target vector x G i and the trial vector u G i is generated using the following scheme [13]: where rand is a uniformly distributed number in the interval (0,1), the crossover rate Cr ∈ [0, 1] controls the how many components are inherited from the mutant vector, j rand a uniformly distributed number in the interval [1, D]. After all of the trial vectors are generated, a selection process determines the survivors for the next generation. If and only if fitness function value of u G i is equal to or better than x G i , then u G i is set to x G+1 i . On the other hand, the old vector x G i is reserved. The selection scheme is for a minimization problem as [13]: In LSHADE, in order to maintain the diversity, external archive A is used. At generation G, target vectors x G i which is worse than the trial vectors u G i are preserved. At anytime, the size of archive exceeds the predefined maximum size, randomly selected individuals are removed in order to make space for newly inserted ones.
The last part of the LSHADE is Linear Population Size Reduction method. In LPSR, the population size decreases with a linear function given in Eq. (49). Here, maxFEs is the maximum number of fitness evaluation, FE is the current number of fitness evaluation, N init is the initial population size, and N min is the minimum number of population size. According to Eq. (49), the population size at generation 1 is N init and the population is N min at the end of the generation [13].
Finally, the pseudo code of LSHADE algorithm is shown in Algorithm 1 [13]. for  [14] in order to improve the search performance of the MHS algorithms. The selection method is important for the performance of MHS algorithms as it enables the selection of the solution candidates that will guide the search process in the most effective way and determine the direction of the search correctly. Selection methods used in the algorithms can be classified into three different categories: the greedy selection method, the random selection method, and the probabilistic selection method (roulette wheel and tournament method). The greedy selection method is the method in which the solution candidate with the highest fitness value is selected among the solution candidates in the population. In the random selection method, the solution candidates are randomly selected from the population. In the probabilistic selection method, the probability of selecting individuals in the population according to their fitness values is determined. For this reason, the individual with the highest fitness value is likely to be selected. The most common used probabilistic selection methods are roulette wheel and tournament methods [14], [45], [46].

Algorithm 1 The Pseudo Code of LSHADE Algorithm
Indeed, FDB is a greedy selection method. The most important difference of FDB from other selection methods is that the score values of the solution candidates are calculated and the selection process is carried out according to the score values. That is, the FDB scores are taken into account instead of the fitness values of the solution candidates. Besides, with this method, the selection of a solution candidate that is very close to the best solution in the population is prevented. Thus, it contributes to the solution of the problem of premature convergence, which is frequently encountered in the MHS process. In the FDB score calculation, the distance and fitness value of the solution candidate to the best solution (X best ), which has the best fitness value in the population, are taken into account. Accordingly, the FDB scores of the solution candidates in a population are explained below [14], [45], [46]: (i) The fitness value of each solution candidate in the population (P) is calculated. The dimension of the problem is n, and m is the number of solution candidates in the population. The i th candidate is denoted as The vector of population P and the fitness values F of the solution candidates is as follows: (ii) The Euclidean distance of each solution candidate X i from the X best is calculated by Eq. (52).
(iii) D P is the distance vector for the population P can be expressed as: (53) VOLUME 10, 2022 (iv) The other parameter used when calculating the FDB score is the fitness value. The fitness values of the solution candidates are represented by the vector F and are given in the Eq. (51). The distance and fitness values should be normalized so that they do not dominate each other in the score calculation. normF and normD P represent the normalized fitness and distance values of individuals, respectively. In the score calculation, normalized fitness and distance vectors are used. FDB score calculation of i th solution candidate is given in Eq. (54). The w is the weight coefficient which denotes the effect of fitness and distance values on FDB score and its value is changed in the range of [0, 1]. If the value of w is equal to 1, the fitness value effect is dominant in the score value. In this case, exploitation effect can be seen in the search process. On the contrary, if the value of w is equal to 0, the score value depends only on the distance value. Thus, the solution candidate farthest from P best is chosen in the population. Thanks to this, exploration capability for MHS algorithms is provided. As a result, w is a dynamic parameter and has the important function of balancing between exploration and exploitation.
The score vector (S P ), which represents the FDB score values of the solution candidates, is given in Eq. (55).
More detailed information about the FDB method can be obtained from the reference study [14].

C. PROPOSED ALGORITHM: FDB-LSHADESPACMA
In this study, the LSHADE algorithm was examined, which is a variant of the DE algorithm and is the winner of the CEC14 Competition on Real-Parameter Single Objective optimization. FDB is a selection method and provides efficient selection of reference positions that guide the algorithms in the MHS process. Therefore, in this study, a hybrid algorithm called FDB-LSHADE was proposed. The goal of the study was to improve the search capability and convergence performance of the LSHADE algorithm. For this purpose, the solution candidates in the mutation process were selected by the FDB method.
In Eq. (45), the mutation strategy of the LSHADE algorithm was given. In Eq. (45), x G pbest is the best individual vector which has the best fitness value in the population at generation G, x G i is the target vector at generation G, x G r1 and x G r2 are the randomly chosen vector, where r 1 and r 2 are randomly chosen indices from the population. Different combinations were created in selecting x G pbest , x G r1 , and x G r2 by using the FDB method. Among them, the best performance was obtained with variants when the FDB method was applied to instead of x G r1 vector and six different FDB-LSHADE variants were explained in this study. The x G FDB is the individual with the highest FDB score among the S P vectors in Eq. (55). The mathematical models of the FDB-LSHADE variants and are presented in the Table 1. When the FDB-LSHADE variants given in the Table 1, there are different values of the weight coefficient w. In the FDB method, the weight coefficient w is an important parameter used for the effects of the candidate's fitness and distance values on the score calculation given in Eq. (54) and its value is in the range of [0, 1]. In FDB-LSHADE variants, there are three different w value was considered. In Cases 1 and 4, w was taken as a constant value 0.5. In Case 2 and 5, w was changed with the function (exp (−FE/maxFEs)) and was changed with the function (1 − exp (−FE/maxFEs)) in Case 3 and 6. The pseudo code of the proposed FDB-LSHADE method is given in Algorithm 2. The FDB method is applied in the selection of given in line 14 of Algorithm 2. In the FDB-LSHADE algorithm, the mutation operator is applied using Eqns. (56) and (57) which FDB-LSHADE variant is used, as given in the line 15 of Algorithm 2.

IV. EXPERIMENTAL SETTINGS
A comprehensive experimental study was carried out in order to test and validate of the proposed FDB-LSHADE algorithm. The experimental study settings are given below: • The settings of the MHS algorithms given in their original articles were taken into account.
• The stopping criterion was used as the maximum number of fitness evaluations (maxFEs) to evaluate the all algorithms under equal conditions. Its value was set to 10000 * Dim and Dim is the dimension of a problem.
• In order to compare the performance of the algorithms, two mostly used benchmark test suites, CEC14 and CEC17, were used.
• Performances of the MHS algorithms change with the different dimension of the problem. CEC14 and CEC17 test functions were tested in 30, 50, and 100 dimensions.
• The MHS algorithms were run 51 times.
• Two mostly used statistical tests methods, Wilcoxon and Friedman, were applied to examine the performance of proposed methods. Wilcoxon test was conducted for a 5% level of significance.
• The experimental studies were implemented in MATLAB R R2016b and run on Intel R Core TM i5-CPU @ 2.90 GHz, 8 GB RAM, and x64-based processor. The information about the CEC14 and CEC17 test functions given in Table 2. In Table 2, number of test functions in benchmark suites, types and search range of the test functions, and dimension of the problems are given.
In order to test and verify the performance of the proposed FDB-LSHADE algorithm, 8 competing MHS algorithms and the base LSHADE algorithm were used in the experimental studies. The parameter settings of the MHS algorithms used in the study are given in Table 3. In determining the parameters of the algorithms, the settings given in their original articles were taken into account.

V. RESULTS AND ANALYSIS
This section consists of three sub-sections. In the first subsection, the performance of the FDB-LSHADE variants was analyzed. For this purpose, six different FDB-LSHADE variations and the LSHADE algorithm were compared in order to determine the most effective algorithm. In the second sub-section, experimental study was carried out to compare the performance of the proposed FDB-LSHADE and other MHS algorithms. In the last experimental study, the EHED problems were solved using the FDB-LSHADE algorithm. The results obtained from the proposed algorithm were compared to the results of the other studies and competing MHS algorithms.

A. DETERMINING THE BEST FDB-LSHADE ALGORITHM ON BENCHMARK SUITES
This sub-section provides the information about the experimental studies for determining the best FDB-LSHADE variant. Statistical analyses were performed among the LSHADE algorithm and six FDB-LSHADE variants (Case-1, Case-2, Case-3, Case-4, Case-5, and Case-6).

1) STATISTICAL ANALYSIS
In this section, firstly, the search performance of the LSHADE algorithm and the FDB-LSHADE variants were analyzed using Friedman test on CEC14 and CEC17 benchmark problems and for three dimensions (D = 30, 50, 100). Using two benchmark suites and three dimensions, six experiments were carried out. The results of Friedman test rankings are given in Table 4. According to Friedman test rankings given in Table 4, the FDB-LSHADE variants achieved the better score than the LSHADE algorithm in the six experiments. The last column of the Table 4 gives the mean ranking information obtained by the algorithms for the six experiments. The Case-2 variation performed more successful performance than the other variants of the FDB-LSHADE algorithm and the LSHADE algorithm and ranked first in mean rank value. The best Friedman score obtained for each experiment in Table 6 is indicated by bold text.
Besides Friedman test, the Wilcoxon test, which allows pairwise comparison between algorithms, was made between the LSHADE algorithm and each variant of the FDB-LSHADE algorithm. The results of the analysis obtained from a total of 58 different test functions are given in Table 5. ''Case-1 vs. LSHADE'' represents for the pairwise comparison between the LSHADE algorithm and the FDB-LSHADE variant. The sum of the numbers given in each cell refers to the total test function of its benchmark. In each cell, three different scores are given. For instance, in the third cell given in Table 5, a result of 19/9/1 is given. This experiment was conducted between the LSHADE and Case-5 in the CEC14 benchmark suite for D = 50. These scores showed that the LSHADE algorithm lost against Case-5 in 19 of 29 test functions, obtained similar results to Case-5 in 9 test functions, and performed better result in 1 problem.
According to the Wilcoxon pairwise test results given in Table 5, the performance of the FDB-LSHADE variants was superior to the LSHADE algorithm in the CEC14 and CEC17 benchmark suites. For example, the pairwise comparison results for D = 100 between the LSHADE and Case-5, the results of Case-5 and the LSHADE were 17/11/1 and 23/6/0, respectively. That is, Case-5 was shown better results in 40 of 58 problems. The pairwise comparison results clearly showed that the FDB based LSHADE variants were stronger than the LSHADE algorithm.

2) CONVERGENCE ANALYSIS
This subsection provides information on the convergence analysis of FDB-LSHADE and its variants. To demonstrate the search performance of the algorithms on different problem types (unimodal/multimodal/hybrid/composition), one of each problem type in CEC14 benchmark suite was selected. Box-plot graphics were drawn based on the error values found for F1 (unimodal), F5 (multimodal), F18 (hybrid), and F21 (composite) type problems of FDB-LSHADE and its variants as a results of 51 independent runs.
The box-plots graphics are shown in Fig. 3 were Unimodal type problem is used to show the exploitation ability of the algorithms. For the F1 unimodal type problem, the LSHADE and its variants converged steadily to the global optimum solution in low-dimensional search space. However, as the problem size increased, the error values obtained from the algorithms increased. For D = 50, Case-5 converged better than the competitor algorithms. For D = 100, Case-3, 4, 5, and 6 show the similar convergence performance among all algorithms.
Multimodal type problem has many local optimum trap. Compared to its competitors, for D = 30, all FDB variants showed the similar performance and good performance than the LSHADE. While Case-2 achieves better results compared to its competitors for D = 50, Case-3 surpassed all algorithms for D = 100.
For F18 hybrid type problem, in low-dimensional search space, for all dimensions, the algorithms performed different convergence performance. For D = 30, the error values of algorithms are quite low. However, when the size of the search space increases, Case-3 showed the steady performance among the competitors.
Composition type problems are very complex problems and used to test the balance between exploration and exploitation of algorithms. According to results of the LSHADE and its variants, for D = 30 dimension, the algorithms obtained almost the same error values. For D = 50, Case-2 performed stable performance than the competitors. In high-dimensional search space, Case-5 had a superior performance to their competitors.
Consequently, the effect of FDB method on the LSHADE for different types of the problems were evaluated. According to experimental studies, with the application of the FDB method, the exploitation, exploration, and balanced search capabilities of the LSHADE had improved. According to Friedman test and Wilcoxon signed-ranked test, Case-5 was the most successful algorithm among the FDB-LSHADE variations and will be referred to as FDB-LSHADE in the following sections.

B. COMPARISON OF FDB-LSHADE AND MHS ALGORITHMS
In this section, the performance of the FDB-LSHADE algorithm was tested, validated, and compared to competing MHS algorithms. Using CEC14 and CEC17 benchmark suites, the performance of the proposed FDB-LSHADE algorithm were compared to the base LSHADE algorithms and 8 of the powerful and up-to-date MHS algorithms. The parameter settings of the competing algorithms were same as the its original paper and detailed information was given in Section 4. In the following two sub-sections, the statistical analysis and convergence analysis of the algorithms are given.

1) STATISTICAL ANALYSIS
In order to evaluate the data obtained from the proposed FDB-LSHADE algorithm and competing MHS algorithms, Friedman test was used. Using three dimensions (D = 30, 50, 100) and two benchmark suites, six experiments were conducted for 58 problems. The Friedman rankings of the competing algorithms and proposed algorithm are given in Table 6. The best Friedman score obtained for each experiment in Table 6 is marked in bold. In the last column of the Table 6, the mean rank score value obtained by the algorithms for the six experiments are given. When the Table 6 was examined,  the FDB-LSHADE algorithm gave the best score value than the competing algorithms. Consequently, the FDB-LSHADE algorithm outperformed its competitors in all experiments and had the best mean rank value.
The Wilcoxon signed-rank test was used to conduct a pairwise comparison between the FDB-LSHADE and the competing algorithms. The pairwise comparison results between the FDB-LSHADE and its competitors are given in Table 7. The sum of the scores given in each cell in Table 7 show the number of test problems, where the FDB-LSHADE lost against its competitor, it performed similar performance with competing algorithm, and it outperformed the its competitor, respectively. When the pairwise results between the FDB-LSHADE and the LSHADE algorithms, in the first cell of the first row, the FDB-LSHADE performed similar result in 14 problems, and was superior in 15 problems. AGDE algorithm was ranked third in the Friedman test in terms of mean rank value. The pairwise comparison results between FDB-LSHADE and AGDE algorithms show that the FDB-LSHADE failed the 2 of the 58 problems for D = 30 dimensions, 1 of the 58 problems for D = 50 dimensions, and 2 of the 58 problems for D = 100 dimensions.
Consequently, the Friedman and Wilcoxon test results showed that the FDB-LSHADE algorithm was superior to its competing algorithms. According to experimental results, it was clearly shown that the performance of the MHS algorithm changed with the problem types and problem dimensions. With this motivation, in this study, it is aimed to propose a competitive algorithm to the literature.

2) CONVERGENCE ANALYSIS
In this sub-section, the convergence curves of the competing MHS algorithms and the proposed FDB-LSHADE algorithm are given. Four different types of problems were selected from test functions in CEC17 benchmark suite, which were the F3 (unimodal), F5 (multimodal), F12 (hybrid) and F21 (composite) types, to show the performance of the algorithms. The convergence curves were drawn for D = 30, 50, and 100 dimensions and are given in Fig. 4.
The convergence curves of unimodal-type F3 problem are shown in Fig. 4 (a), (b), (c). According to graphs for low dimensional problems (D = 30), the error values of the FDB-LSHADE, the base LSHADE, and AGDE algorithms decreased to 10 −15 . For D = 50, while the error values of the competing algorithms except the LSHADE were approximately 10 0 , the FDB-LSHADE algorithm had the best minimum error value with 10 −15 compared to competing algorithms. However, for D = 100, the only algorithm with the error value less than 10 −5 was the FDB-LSHADE algorithm. These results indicated that the exploitation ability of the FDB-LSHADE algorithm was superior to competing algorithms.
For multimodal-type F5 problem, the convergences curves of the algorithms are given in Fig. 4 (d), (e), (f ). Multimodaltype problems have many local solution traps. Therefore, the elimination of these traps depends on the exploration ability of the algorithms. When the Fig. 4 were examined, depending on the size of the problem, all algorithms except the FDB-LSHADE searched successfully until they reached a certain error value. For low-dimensional problem, the error values of  the both LSHADE and the FDB-LSHADE were lower than the 10 1 . However, the error values of the algorithms were greater than the 10 1 for D = 30, 50 and 10 2 for D = 100. That is, these competing MHS algorithms performed premature convergence in the multimodal problem space. On the contrary, the error value of the FDB-LSHADE algorithm were less than the 10 1 for D = 50 and 10 2 for D = 100. This results demonstrated that the exploration capability of the FDB-TLABC was stronger than its competitors.
In Fig. 4 (g), (h), (i), the convergence curves of the algorithms are given for the hybrid-type F12 problem. Hybrid-type problems are used to test the balanced search capabilities of the algorithms. When the convergence curves of the algorithms were analyzed, The FDB-LSHADE exhibited a performance relatively superior to its competitors.
In Fig. 4 (j), (k), (l), the convergence curves of the algorithms are shown for the composite-type F21 problem. Composite type problems are the highly complexity problem types and are more difficult to optimize than others. When the graphs were analyzed, the error values of the FDB-LSHADE algorithm were less than the its competitors.
In summary, the results of the examination of the convergence curves indicated that the FDB-LSHADE had powerful exploitation and exploration capabilities.

C. DETERMINING THE BEST METHOD FOR SOLVING THE EHED PROBLEMS
In this study, the proposed FDB-LSHADE, AGDE, MPA, AEFA, SDO, EO, AEO, SSA, and MFO algorithms were applied to the EHED problems. The test system is highly complex, non-linear and non-convex and high-dimensional with 103 control variables. In this test system, 76 of the control variables are belongs to the sources including 27 electrical sources, 34 gas sources, and 15 heat sources, and the remaining control variables are the dispatch factors of the energy hub structures. The demand values including the electricity, heat, cooling and compressed air from the system are 12 pu, 9 pu, 0.7 pu and 1.2 pu, respectively. The system data is given in Appendix.
Three objective functions were optimized in this study: (i) minimization of total energy cost, (ii) minimization of total energy hub loss, and (iii) minimization of total energy cost and hub loss. Three different cases considered as objective functions were analyzed and the simulation studies were identified according to the following test cases.
• Case-1: Optimizing the a cost function including the cost model of electricity, natural gas, and heat.
• Case-2: Optimizing the energy hub losses • Case-3: Optimizing a multi-objective problem including energy cost and hub loss, which is converted to a singleobjective problem VOLUME 10, 2022

1) CASE-1: MINIMIZATION OF TOTAL ENERGY COST
The minimization of the total cost value using the cost model of electricity, natural gas, and heat is studied in Case-1. The cost function given in Eq. (36) is considered as objective function to be minimized by proposed FDB-LSHADE and other MHS algorithms. The optimal solutions obtained by FDB-LSHADE are given in Table 8 for all cases.  the results of the Case-2 given in Table 9

4) LITERATURE COMPARISON, STATISTICAL ANALYSIS, AND CONVERGENCE ANALYSIS
The results of the proposed FDB-LSHADE algorithm for the EHED problems and of the TVAC-GSA and SAL-TVACGSA algorithms lately reported in the literature are given for Case-1, Case-2, and Case-3 as shown in Table 10. For Case-1, end for 17. for end if 27. end for 28.
If necessary, delete randomly selected individuals from the archive such that the archive size is |A|.

29.
Update memories M F and M CR 30.
Sort individuals in P based on their fitness values and delete worst N G − N G+1 members; 33.
Resize archive size |A| according to new |P| 34.
G + +; 35. end while the minimum cost value obtained by the FDB-LSHADE algorithm were 15451.9174 mu, which was 1.8758% and 1.7578% lower than the TVAC-GSA and SAL-TAV-GSA algorithms. The proposed algorithm obtained a feasible solution with 2.5524 pu for Case-2 was lower by 9.8618% and 9.1657% lower than the TVAC-GSA and SAL-TAV-GSA algorithms. For Case-3, the minimum objective value of the FDB-LSHADE algorithm was 16679.3993 mu, which was 8.3522% and 8.2019% lower than the TVAC-GSA and SAL-TAV-GSA algorithms.
In this study, three test cases were evaluated by the FDB-LSHADE, LSHADE, MPA, EO, AEFA, AEO, SDO, AGDE, SSA, and MFO algorithms. In Table 9, the minimum, mean, maximum, and standard deviation values of the 30 runs for each test case of all the optimization algorithms. According to these results, the advantage and efficiency of the FDB-LSHADE algorithm was demonstrated.  To evaluate the results of the FDB-LSHADE, LSHADE, MPA, EO, AEFA, AEO, SDO, AGDE, SSA, and MFO algorithms, Friedman test was applied on each test cases and the ranking results of the those algorithms were given in Table 11.  According to Friedman test results, the FDB-LSHADE was ranked first among the algorithms.
In this sub-section, the convergence curves of nine different optimization algorithms and the proposed FDBAGDE algorithm are shown in Figs. 5-7. The convergence curves given in Figs. 5 -7 clearly showed that the convergence performance of the FDB-LSHADE algorithm was better than the others in terms of convergence accuracy and speed. Consequently, the convergence curves indicated that the proposed method outperforms other algorithms in terms of best fitness value, convergence speed, and robustness for all test models.
The boxplot graphs of the algorithms provide easy understanding and interpretation of search performances. Boxplots graphs of the nine competing algorithms and the FDB-LSHADE algorithm for all case studies are given in Fig. 8. The examination of the graphs demonstrated that the minimum, maximum, mean and standard deviation margins of the FDB-LSHADE and LSHADE algorithms were reasonable; however, in all cases, the FDB-LSHADE was able to find better results than the LSHADE. When all five of the boxplots were examined, it was seen that the FDB-LSHADE demonstrated more stable and robust search than its competitors in all three cases.

VI. CONCLUSION
In this study, an improved version of the LSHADE algorithm called the FDB-LSHADE was presented. The use of the FDB selection method for redesigning of the mutation operator had enabled the algorithm to effectively establish the exploration and exploitation balance. The results of experimental studies carried out on CEC14 and CEC17 benchmark suites of different types showed that the proposed FDB-LSHADE algorithm had achieved a superior performance compared to its competitors.
The experimental studies in this study were carried out in three sections. In the first section, the six variants of the FDB-LSHADE algorithm were compared to the base LSHADE algorithm. The performance of these algorithms was tested on the four different types problem, which are unimodal, multimodal, hybrid, and composite, in the CEC14 and CEC17 test suites. The results obtained from the algorithms were evaluated using statistical test methods. The strongest FDB-LSHADE variant was determined according to comparing results. In the second section, the performance of the proposed FDB-LSHADE algorithm was compared to up-to-date and the most preferred MHS algorithms including the base LSHADE were compared. According to statistical test results, while the LSHADE algorithms ranked second, the proposed FDB-LSHADE ranked first in all dimensions of the test suites. These results show that the changes in the design of the LSHADE algorithm had been successful. Finally, the proposed FDB-LSHADE algorithm was applied to the solving EHED problems, which are non-convex, nonlinear, and high-dimensional. Three different EHED problems were solved by using the proposed algorithm and its competitors. According to the simulation results, it was proven by the statistical analysis methods that the FDB-LSHADE was superior to the LSHADE, MPA, EO, AEFA, AEO, SDO, AGDE, SSA, and MFO algorithms. Besides that, the FDB-LSHADE algorithm was found better results for EHED problems than the optimum solutions obtained from previously conducted studies.
Summing up, the FDB-LSHADE algorithm developed as a result of the comprehensive experimental study was presented to the literature as one of the most robust MHS methods, which can be used to solve unconstrained and constrained optimization problems.

APPENDIX
In this study, the energy hub system consists of 76 sources, where 27 electrical sources, 34 gas sources, and 15 heat sources. The cost coefficients and energy production limits of these sources are given in Table 12. Besides that, in the system, the efficiency coefficients of the hubs are given in Table 13.