Optimized Adaptive Overcurrent Protection Using Hybridized Nature-Inspired Algorithm and Clustering in Microgrids

Microgrids have been popularized in the past decade because of their ability to add distributed generation into the classic distribution systems. Protection problems are among several other problems that need solutions in order to extract the full capability of these novel networks. This research follows the branches of two major solutions, namely adaptive protection and protection optimization. Adaptive protection implementation with a novel concept of clustering is considered, and protection setting optimization is done using a novel hybrid nature-inspired algorithm. Adaptive protection is utilized to cope with the topology variations, while optimization techniques are used to calculate the protection settings that provide faster fault clearances in coordination with backup protection. A modified IEEE 14 bus system is used as the test system. Validation was done for the fault clearing performance. The selected algorithm was effective than most other algorithms, and the clustering approach for adaptive overcurrent protection was able to reduce the number of relays’ setting groups.


Microgrid Protection Challenges
The microgrids have accumulated a considerable amount of interest within the past years and turn out to be an indispensable asset in the power industry. The capability to incorporate renewable power generation into the distribution system is one of the vital reasons for microgrids recognition. A broad array of distributed generation (DG) technologies such as sustainable micro-turbine generation, including wind generation, photovoltaic generation, and energy storage systems, make the microgrid more feasible in both islanded and grid-connected modes [1].
There are several technical disputes to be confronted when trying to access the entire potential of microgrids, and protection is one of those challenging areas [2]. Many solutions were presented, influenced by the improvement of protection techniques. Microgrids containing DG can cause variations in short circuit levels, which is one of the main reasons for these protection challenges. The ability of microgrids to operate islanded of the main utility can make these variations more drastic. Also, the existence of different operating topologies makes the protection selectivity more complicated.
The existence of DG and different operating topologies mainly causes the technical challenges that microgrids face. Protection challenges can be identified as a major area that requires specialized solutions in order to manage efficient running microgrids. A few of the main protection challenges can be identified as follows.

Methodology and Paper Structure
This paper proposes a novel protection scheme for microgrids inspired by adaptive protection and optimization techniques. Section 2 describes the introduction of directional overcurrent protection and different variations of adaptive protection implementations. Furthermore, the clustering methods are proposed to group similar topologies together to reduce the complexity of the adaptive settings. Then the protection coordination using optimization algorithms is introduced, and initial stages of problem formulation are presented. In Section 3, the system model used in the simulation is presented. Initial data for the optimization problem are also presented along with the objective function formulation and constraint formulation. Then the selected algorithm and topology clustering method is introduced. In Section 4, simulation results are presented and discussed. Section 5 elaborates on the conclusion and other prominent research achievements, followed by recommendations for future improvements.

Overcurrent Protection
Protection of the distribution level networks, including microgrids, are usually handled through four types of relays. These four types can be identified as overcurrent, distance, voltage based, and differential relays. This research is mainly focused on distribution line protection. Overcurrent and distance are the most commonly used types for line protection. Out of these two, the overcurrent solution is more straightforward and cost-efficient. Even though it is the plainest option, it is more challenging to achieve the protection of novel distribution networks with varying topologies and embedded generation.
Because of the meshed structure and the presence of DG in these networks, directional over current relays (DOCR) are utilized more commonly. This directionality function is extremely important in a meshed system to correctly discriminate the origin of faults. Relays will not operate for fault currents in the "reverse" direction.
The operation of an overcurrent relay is governed by its characteristic curve. The inverse time curve is the widely used type where the operating time is increasing with the magnitude of the fault current. The operating time of a relay can be expressed by Equation (1), with several relay parameters and fault current according to IEC 60255 standard [25]. These parameters are listed in Tables 1 and 2.

Adaptive Overcurrent Protection
In this paper, adaptive overcurrent protection is implemented with an offline calculation of protection settings. These pre-calculated relay settings are applied according to the operating mode. A set of profound network topologies must be identified, and the respective protection settings for each state are recalled through an event table.
Operating topologies can be identified by considering single contingencies. Each line, generator, and transformer can be disconnected one at a time to obtain new topologies. While it is possible to identify more meaningful topologies by using double contingency, the obtained topologies are sufficient for validating the proposed design.

Clustering of Topologies
In convention, there is the same number of protection setting groups as the number of identified operating topologies. By clustering a similar set of topologies together, they can be operated using the same protection setting group. This effectively reduces the number of protection setting groups to be stored in the earlier mentioned event table.
Clustering is a phenomenon used in data science to group similar data objects. It is primarily a search function that is used to identify similar and dissimilar objects and group them accordingly [26]. Data clustering is vastly utilized in the fields of data mining and pattern recognition. There are several methods to develop clusters, such as hierarchical methods, partition methods, grid-based methods, model-based methods, etc. The k-means clustering algorithm is a widely used method used for data clustering. The main objective of this algorithm is to minimize the distance from each data point to the center of its cluster. Figure 1 illustrates the basic flowchart of the k-means algorithm.
Energies 2020, 13, x FOR PEER REVIEW 5 of 23 In convention, there is the same number of protection setting groups as the number of identified operating topologies. By clustering a similar set of topologies together, they can be operated using the same protection setting group. This effectively reduces the number of protection setting groups to be stored in the earlier mentioned event table.
Clustering is a phenomenon used in data science to group similar data objects. It is primarily a search function that is used to identify similar and dissimilar objects and group them accordingly [26]. Data clustering is vastly utilized in the fields of data mining and pattern recognition. There are several methods to develop clusters, such as hierarchical methods, partition methods, grid-based methods, model-based methods, etc. The k-means clustering algorithm is a widely used method used for data clustering. The main objective of this algorithm is to minimize the distance from each data point to the center of its cluster. Figure 1 illustrates the basic flowchart of the k-means algorithm. Initially, the data points, number of desired clusters, maximum iterations, and stating points of centers are determined. Then the distances from every single data point to the closest center is calculated from Equation (2). Dij is the distance between datapoint x and initial center point c, where d is the number of dimensions in the selected space. After calculating the distances, the average of data in the whole cluster is calculated by Equation (3). Here m is the dimension, j is the cluster no, n is the number of data points in the cluster.
After that, the new averages are assigned as centers, and the process is repeated until the desired iterations, or there is no considerable change in centers between iterations. Initially, the data points, number of desired clusters, maximum iterations, and stating points of centers are determined. Then the distances from every single data point to the closest center is calculated from Equation (2). D ij is the distance between datapoint x and initial center point c, where d is the number of dimensions in the selected space. After calculating the distances, the average of data in the whole cluster is calculated by Equation (3). Here m is the dimension, j is the cluster no, n is the number of data points in the cluster.

Protection Coordination with Optimization
Energies 2020, 13, 3324 6 of 23 After that, the new averages are assigned as centers, and the process is repeated until the desired iterations, or there is no considerable change in centers between iterations.

Protection Coordination with Optimization
Overcurrent relay coordination is done by converting it into an optimization problem. As per the conventional structure of an optimization problem, it includes an objective function, a set of constraints, and the algorithm that is used to generate the solution. The minimization approach is used as the problem is handling time quantities, and operational times should be minimum for an optimal solution.

Different Objective Functions
The sum of relay operating times for a defined fault is the most widespread objective function (OF) utilized by the literature Equation (4) [27][28][29]. The probability coefficient (w) is used to represent each fault occurrence probability.
The probability coefficient w is typically set as 1, which forces all faults to exhibit an equal rate of occurring. The reflecting three-phase to ground fault for each primary protection zone can either be a near-bus fault Equation (5) [30,31], or both near-bus and far-bus faults Equation (6) [32,33].
A fault occurring at the beginning of a line near the relay is considered a near-end fault. This near end fault has the maximum short circuit current level for the considered relay. A far-end fault provides the minimum short circuit value for the considered relay while occurring at the furthest end of the line. Operating time for each near-bus fault and far-bus fault can be denoted as t (i,near) and t (j,far) , respectively. When both near-bus and far-bus faults are included in the OF, the optimization is expected to become more generalized.
The third type of common OF is the sum of squares of the relay operating times, as in Equation (7) [34,35]. This is mainly used to minimize the backup relay operating times while prioritizing the primary operating times with the squares. The coordination time interval (CTI) constraints are included within the OF second summation term. Here ∆t j is for the relay coordination as in Equation (8). The first term can be used alone for primary relay settings with coordination term included in the constraints as in Equation (9).
It was found that these different objective functions do not make a considerable impact on the final results [36]. Therefore, the popular choice of using near-bus fault Equation (5) was used in the proposed scheme.

Optimization Algorithms and Constraint Handling
A large number of optimization algorithms have been used in relay coordination throughout the years. Different optimization algorithms can produce different optimized datasets as output values. Better the algorithm performance, the more optimal the results will be. To identify the degree of optimization for a problem, its OF value can be used as a performance indicator. Hybrid NIA can be identified as a novel type of algorithms that combines two or more primary algorithms. They tend to perform well in obtaining the best optimal solution when compared with other algorithms. The proposed scheme utilizes a hybrid algorithm to get a better solution set.
The algorithms are used in both linear and non-linear problem formulation. In linear scenarios, only a single variable is considered in the optimization, usually the time dial setting (TDS). In these types of linear cases, the other variables, namely pick-up setting (PS) and the curve type are pre-selected using either conventional methods or previous experiences [37]. In non-linear formulation, multiple variables can be included in the optimization. Usually, both TDS and PS are optimized simultaneously [38]. The proposed scheme uses a non-linear optimization formulation.
Multi-objective optimization (MOO) is required when there are multiple objectives to be improved. The relay coordination problem is a multi-objective problem as there are multiple relays to be considered. When performing a MOO simultaneously, compromises are being made in order to achieve the optimal set of solutions. This set of solutions is called the Pareto front. In attaining the required Pareto front solution, there are two popular methods as multi-objective decision-making (MODM) practices and multi-objective optimization techniques.
When using the MODM method, the multi-objective problem is converted into a single objective model, usually through the weighted sum method [39]. Different weights are assigned according to the importance of the objective and added together to formulate a single objective function. At each run of the optimization algorithm, it will provide a solution from its respective Pareto optimal solution set.
The multi-objective optimization technique utilizes a particular multi-objective version of the metaheuristic algorithm. Most of the novel optimization algorithms have their own multi-objective version proposed by the literature. This method is not widely used because of its long execution times and more extensive standard deviation of the solution vector. They usually exhibit weaker algorithm performance when compared with their single objective counterpart [40]. The proposed scheme uses MODM with a weighted sum method.
When solving an optimization problem, there can be constraints that limit the feasible solution space. The constraints should not be violated for a feasible solution. There are two main types of constraints as equality constraints and inequality constraints.
For the above optimization of function Z Equation (4), there are q number of inequality constraints g(x) Equation (10) and m number of equality constraints h(x) Equation (11). Optimization algorithms handle constraints in diverse ways. Some algorithms have inherited constraint handling, where it only requires defining the equality and inequality constraint functions. For example, the GA can perform constrained optimization without any modification of the usual input syntax. Some algorithms cannot process constrained optimization problems and need to use modifications in input data to obtain a feasible solution set, i.e., PSO. The selected algorithm has inherited constraint handling capability, which enables fewer constraint violations than using the penalty method.

Problem Formulation and Simulation
The modified IEEE 14 bus system was considered as the microgrid under testing. Different topologies are identified, and fault current and load flow data for each topology has been obtained. The optimization problem is then formulated and solved in order to obtain protection settings for each operating topology.

Modified IEEE 14-Bus System
The 33 kV Distribution section of the IEEE 14 bus system was used as the test model to verify the protection solution. As shown in Figure 2, the 33 kV section is connected with a high voltage section through transformers from bus 6 and bus 9. When disconnected from these links, the distribution section acts as an islanded microgrid. There is only one synchronous generator connected to bus 6 in the original 33 kV section (Table 3). To increase the DG penetration of the system, two similar generators were added to bus 13 and bus 9. Network line parameters and general load parameters are listed in Tables 4 and 5 (Per-unit values are calculated on a 100 MVA base). Complete system parameters can be found in [41].
Energies 2020, 13, x FOR PEER REVIEW 8 of 23 The 33 kV Distribution section of the IEEE 14 bus system was used as the test model to verify the protection solution. As shown in Figure 2, the 33 kV section is connected with a high voltage section through transformers from bus 6 and bus 9. When disconnected from these links, the distribution section acts as an islanded microgrid. There is only one synchronous generator connected to bus 6 in the original 33 kV section (Table 3). To increase the DG penetration of the system, two similar generators were added to bus 13 and bus 9. Network line parameters and general load parameters are listed in Tables       Each line is equipped with two DOCRs near each bus, and their operating directions are pointing away from the busses. For a fault in the line, the two relays will operate as primary protection, and respective backup relay can be identified from the adjacent lines while considering the directionality.

Single Element Contingencies
In order to represent the variation of topologies, single element contingencies for the modified 14 bus system were considered. New topologies are obtained by disconnecting each line, generator, and utility connection one at a time (Table 6). Eight topologies are obtained by disconnecting one line at a time. Ninth topology is the one with everything connected. Three more topologies can be obtained by disconnecting each generator and utility connection. Additional topology was obtained by disconnecting all three DGs to incorporate a larger variation in short circuit current levels. Line 13-6 10 Gen-bus 13 4 Line 13-14 11 Gen-bus 9 5 Line 14-9 12 Gen-bus 6 6 Line 9-10 13 Utility 7 Line 10-11 14 All 3 DGs

Load Flow and Short Circuit Analysis
The above modified 14 bus system is simulated in DigSILENT Powerfactory software (DIgSILENT GmbH, Gomaringen, Germany). Load flow analysis is done for each topology to obtain current flowing through each relay at its nominal operation. During the topology 13, load curtailing was done to balance the generation and demand. Every load value was adjusted by a scale factor of 0.3 only in topology 13. Load flow values are used as lower limits when calculating the pickup settings of the relays. These values are indicated in Table 7. Current transformer (CT) ratios are also indicated in the table, and all current values in the table are in amperes. Throughout this paper, R is used to denote DOCRs followed by the relevant relay number (e.g., R2, DOCR number 2), and T is used to denote the topology number (e.g., T4, topology number 4 as in Table 6). Table 7. Load flow analysis for normal operation.  T01  T02  T03  T04  T05  T06  T07  T08  T09  T10  T11  T12  T13  T14   1  -103  226  75  132  94  108  112  102  132  106  98  11  Problem formulation for optimization requires short circuit current values through each relay. Short circuit analysis is done for the test system model, according to IEC60909 [42] because of the presence of distributed generation [43]. Three-phase to ground fault is considered as the maximum short circuit value. Near-bus faults are taken as the maximum fault current value while the far-bus fault current is considered as the minimum.
For the modified 14 bus system, short circuit analysis was done for all selected topologies. All primary backup relay pairs were considered, and both near-bus and far-bus 3 ph fault values were recorded. Table 8 presents the fault analysis for topology no 1.

Optimization Algorithm Implementation
Optimization algorithm implementation can be separated into two areas as objective function formulation and constraint formulation. Relay operating time is the principal consideration of the optimization, and the minimization approach is used. Short circuit data and load flow data are used in this formulation while the minimization problem is solved to find the best set of values for the decision variables.

Objective Function Formulation
Objective function formulation can be considered the fundamental stage of the optimization process. As described in Section 2.3.1, the summation of the relay operating times is used to form the objective function. The final goal is to minimize this summation Equation (12), subject to different constraints.
From the Equation (13) for the operating time of a relay, TDS and PS are considered variables. They are called decision variables, and the final solution is to find the optimal set of values for these decision variables. As for the curve type, the IEC standard inverse curve is selected.
Summation of all relay operating times for maximum short circuit currents is taken. Fault current values are substituted from the short circuit analysis, and the time-current curve type was taken as the standard inverse. For minimization, we get a function of TDS and PS Equation (14).

Constraint Formulation
There are two types of constraints applied in the optimization process to achieve an accurate set of solutions. They are decision variable bounds and functional constraints.
Decision variable bounds are applied to the decision variables according to their physical limitations. TDS values are limited between 0.05 and 1.1 according to relay manufacturer limitations, and PS values are bound between overload currents and minimum fault currents [28,30,44] (Equations (15)- (19) Functional constraints are used to keep the coordination between primary and backup relays. The operation time difference between a primary relay and a backup relay is set to be greater than the CTI Equation (20). CTI value was selected as 0.2 s. Relay operating times are also limited to achieve safe fault clearances as in Equation (21). 0.05 ≤ t i,op ≤ 1.00 (21) This system has a total of 16 relays; therefore, with two variables per each relay makes the total decision variables to 32. An objective function is formed by taking the summation of maximum fault operating times of each relay. Both decision variable bounds and functional constraints are formulated with the load flow data and short circuit data from Section 3.2. Primary relay operating time limits are also included in the main constraint function.

WCMFO Algorithm
Hybrid water cycle-moth flame optimization algorithm (WCMFO) is a combination of two primary algorithms, namely the water cycle optimization algorithm and the moth flame optimization algorithm. This algorithm was initially proposed by Khalilpourazari et al. [45].
Water cycle algorithm (WCA) is a novel optimization algorithm inspired by the flow of water streams toward the sea. Figure 3a illustrates a simple schematic of the algorithm composition. This algorithm includes evaporation and raining conditions to provide necessary randomization between iterations. This feature enables the optimization results to avoid local bests and get closer to the global solution.
Functional constraints are used to keep the coordination between primary and backup relays. The operation time difference between a primary relay and a backup relay is set to be greater than the CTI Equation (20). CTI value was selected as 0.2 s. Relay operating times are also limited to achieve safe fault clearances as in Equation (21).
0.05 ≤ t i,op ≤ 1.00 (13) This system has a total of 16 relays; therefore, with two variables per each relay makes the total decision variables to 32. An objective function is formed by taking the summation of maximum fault operating times of each relay. Both decision variable bounds and functional constraints are formulated with the load flow data and short circuit data from Section 3.2. Primary relay operating time limits are also included in the main constraint function.

WCMFO Algorithm
Hybrid water cycle-moth flame optimization algorithm (WCMFO) is a combination of two primary algorithms, namely the water cycle optimization algorithm and the moth flame optimization algorithm. This algorithm was initially proposed by Khalilpourazari et al. [45].
Water cycle algorithm (WCA) is a novel optimization algorithm inspired by the flow of water streams toward the sea. Figure 3a illustrates a simple schematic of the algorithm composition. This algorithm includes evaporation and raining conditions to provide necessary randomization between iterations. This feature enables the optimization results to avoid local bests and get closer to the global solution.
(a) (b) Moth flame optimization (MFO) algorithm is inspired by the flying patterns of moths around a light source. Figure 3b illustrates the fly patterns used for the algorithm. Each moth represents a solution and flies around their corresponding flame. With each algorithm iteration, the number of flames is reduced to obtain the final solution.
A combination of these algorithms makes a high-performance metaheuristic algorithm. WCA has improved solution space exploration ability constructed on its streams and river formation. MFO has a great ability in exploitation because of the localized formation of flames. WCMFO is designed to combine these advantages of the underlying algorithms. WCA is the base algorithm of the A combination of these algorithms makes a high-performance metaheuristic algorithm. WCA has improved solution space exploration ability constructed on its streams and river formation. MFO has a great ability in exploitation because of the localized formation of flames. WCMFO is designed to combine these advantages of the underlying algorithms. WCA is the base algorithm of the WCMFO hybrid algorithm. The stream positions are then updated according to the spiral movement from the MFO algorithm. WCMFO has an improved raining process to increase the randomization using levy flight.

Clustering for Adaptive Protection
In the modified IEEE 14 bus system model, optimization was done using the WCMFO algorithm for different operating topologies. Instead of calculating a separate set for every topology, clustering was used to identify similar topologies and they were given the same set of protection settings.
As described in Section 3.1, there are 14 topologies under consideration. Most industrial overcurrent relays have facilities to save at least four setting groups. Here the available 14 topologies are divided into four clusters according to their overall short circuit levels. The number of clusters to be created will depend on the relay's ability and deviation of the topologies.
The deviation of the short circuit current for each relay at every topology is considered the clustering factor. The mean value of each primary relay's maximum fault current was calculated using Equation (22). Here, Ri is the relay number, and I f, Ri, Tx is the maximum fault current (3-ph near-bus) of relay Ri for the operating topology Tx. N non-zero is the population size without including zero-current scenarios (e.g., for R1, N non-zero is 13 without including the T1 scenario where R1 is disconnected). These mean values are presented in Table 9. The sum of deviation of currents for each topology was calculated using Equation (23). Here, Tx is the topology number, and I f, Ri, Tx is the same as above. Deviation of each relay's maximum fault current from its calculated mean (M Ri ) was aggregated for each topology, obtaining the summation of deviations for a given topology. These values are presented in Table 10.  Figure 4 provides the summation of deviations for each topology. Here we can identify four distinct groups. After deciding on the number of clusters, each cluster can be finalized using the k-means algorithm. Topology no 14 is separated into its own cluster, and topologies 1, 2, 4, 5, 6, 7, 8, 9 are grouped into another cluster. The remaining two clusters are selected as topologies 3, 11, and topologies 10, 12, 13. Table 11 presents the grouping of topologies into clusters.   The objective function was formed with the maximum currents of the set of topologies belonging to the relevant cluster. When imposing constraints, some topologies had similar constraints, which had to be put only once effectively reducing the number of constraints. Minimum and maximum operating times constraints were also imposed for the minimum and maximum fault values considering the whole cluster. Decision variable bounds were set considering all the values ranging throughout the cluster.

Results and Discussion
Optimization is run on MATLAB software to obtain the best data set for the decision variables. These obtained relay settings are then put into the Equations, and the validity of operating time for the faults are observed. The settings are further applied to the DigSILENT simulation model, and relay time-current curves are also obtained.
Several popular algorithms were used to compare the performance of the WCMFO algorithm. These algorithms and corresponding references are presented in Table 12   The objective function was formed with the maximum currents of the set of topologies belonging to the relevant cluster. When imposing constraints, some topologies had similar constraints, which had to be put only once effectively reducing the number of constraints. Minimum and maximum operating times constraints were also imposed for the minimum and maximum fault values considering the whole cluster. Decision variable bounds were set considering all the values ranging throughout the cluster.

Results and Discussion
Optimization is run on MATLAB software to obtain the best data set for the decision variables. These obtained relay settings are then put into the Equations, and the validity of operating time for the faults are observed. The settings are further applied to the DigSILENT simulation model, and relay time-current curves are also obtained.
Several popular algorithms were used to compare the performance of the WCMFO algorithm. These algorithms and corresponding references are presented in Table 12. The cluster 3 data were used for the comparison calculation. Obtained protection settings and OF values of each algorithm are depicted in Table 13. It is evident that WCMFO has the lowest OF value of all comparing algorithms. OF value can be identified as the total operating time of relays, and the net gain in operating time is compared in Table 14. This is about a 54% decrease in total relay operating time from the worst solution which is the PSO data set. It is a 5.3% decrease from the second-best solution obtained through GWO.   When comparing the computational burden of each algorithm, PSOGSA and HHO took the least central processing unit (CPU) time to converge into a solution within 500 and 1000 iterations, respectively. PSO algorithm took 1103 s of CPU time to converge from 8000 iterations. GWO and WOA were executed with 10,000 iteration limits and took 2105 s and 1530 s, respectively. WCMFO WCMFO was used to obtain the protection settings for the rest of the clusters. The full list of relay settings is presented in Table 15. Here PS values are presented as CT secondary current values. Topology 1 relay operating times were obtained by using the settings given by cluster 1 solution. Full relay operating times for each maximum 3 ph fault for every relay is presented in Table 16. Currents through both primary relay and backup relay are listed under the fault current column. Both primary and backup relay operating times are included under the operating time column, and CTIs are calculated. Primary relay operating times does not exceed the operating limits of 1 s or 0.05 s. CTI values also do not fall below 0.2 s minimum limit.  1  5  -----1  14  -----2  4  -----3  1  0  0  ---4 Table 17 presents the fault currents and operating times for the relay pair 4 and 7 throughout the thirteen topologies. Relay 4 is the primary relay for faults between bus 12 and 13, and the corresponding backup relay is relay 7. Fault currents seen by the primary relay from both near-end and far-end faults are considered, and both operating times are within the safe limits. Figure 5a presents the Relay 4 fault current, and Figure 5b depicts the corresponding relay operating times for each topology.  Table 17 presents the fault currents and operating times for the relay pair 4 and 7 throughout the thirteen topologies. Relay 4 is the primary relay for faults between bus 12 and 13, and the corresponding backup relay is relay 7. Fault currents seen by the primary relay from both near-end and far-end faults are considered, and both operating times are within the safe limits. Figure 5a presents the Relay 4 fault current, and Figure 5b depicts the corresponding relay operating times for each topology.  Figure 6a presents fault currents seen by the backup relay R7, and Figure 6b illustrates the relevant operating times for every topology. Even though the backup relay operating times were not limited to 1 s constraints, they are within acceptable limits and follow perfect coordination. Table 18 depicts the CTI values for both near-end faults and far-end fault. Figure 6a presents fault currents seen by the backup relay R7, and Figure 6b illustrates the relevant operating times for every topology. Even though the backup relay operating times were not limited to 1 s constraints, they are within acceptable limits and follow perfect coordination. Table 18   The following graphs (Figures 7-10) were obtained by feeding the optimized settings back to the DigSILENT model. Three-phase to ground fault was made at the middle of the line (50%-line length), and relay operation curves were obtained. Relay current values are represented in the x-axis in amperes, and the y-axis is time in seconds with log scaling. The red curve represents the primary relay time-current curve, and the respective current seen by the primary relay is represented in the red vertical line. Operating time for this primary operation is at the intersection point. Green curves are backup relay curves, and when there is more than one backup relay, the blue curve also represents the backup curve. If the fault current seen by the backup relay is different from the fault current seen by the primary relay, they are represented in vertical lines of relevant color, and their intersection points are the relevant operating times. Figure 7a presents the R1, R5, and R14 characteristics and currents through each relay when a fault occurs within the primary protection area of R1 (the line from bus 12 to bus 6). A fault is created at 50% of line length, and the considered topology is T2. Here R5 and R14 operate as backup protection according to their own characteristics and fault currents. Figure 7b presents the curves for R4, R7, and R16 and currents through each relay when a fault occurs within the primary protection area of R4 (the line from bus 12 to bus 13). A fault is created at the midpoint of the line as previously, and the considered topology is T9. Here R7 and R16 operate as backup protection according to their own characteristics and fault currents. Similarly, Figure 8a shows the same fault scenario for R7 and R9, where the topology is T3, and the primary relay is R7. Figure 8b shows the fault scenario for R15,  The following graphs (Figures 7-10) were obtained by feeding the optimized settings back to the DigSILENT model. Three-phase to ground fault was made at the middle of the line (50%-line length), and relay operation curves were obtained. Relay current values are represented in the x-axis in amperes, and the y-axis is time in seconds with log scaling. The red curve represents the primary relay time-current curve, and the respective current seen by the primary relay is represented in the red vertical line. Operating time for this primary operation is at the intersection point. Green curves are backup relay curves, and when there is more than one backup relay, the blue curve also represents the backup curve. If the fault current seen by the backup relay is different from the fault current seen by the primary relay, they are represented in vertical lines of relevant color, and their intersection points are the relevant operating times. Figure 7a presents the R1, R5, and R14 characteristics and currents through each relay when a fault occurs within the primary protection area of R1 (the line from bus 12 to bus 6). A fault is created at 50% of line length, and the considered topology is T2. Here R5 and R14 operate as backup protection according to their own characteristics and fault currents. Figure 7b presents the curves for R4, R7, and R16 and currents through each relay when a fault occurs within the primary protection area of R4 (the line from bus 12 to bus 13). A fault is created at the midpoint of the line as previously, and the considered topology is T9. Here R7 and R16 operate as backup protection according to their own characteristics and fault currents. Similarly, Figure 8a shows the same fault scenario for R7 and R9, where the topology is T3, and the primary relay is R7. Figure 8b shows the fault scenario for R15, R5, and R2, where the topology is T11, and the primary relay is R15. Similarly, Figures 9 and 10 present  characteristic curves and similar fault scenarios for T10, T13, and T14. Energies 2020, 13, x FOR PEER REVIEW 19 of 23 R5, and R2, where the topology is T11, and the primary relay is R15. Similarly, Figures 9 and 10 present characteristic curves and similar fault scenarios for T10, T13, and T14. It is evident that the settings group for each cluster is valid for all topologies belonging to that cluster. This clustering approach vastly reduces the number of adaptive stages while grouping similar topologies.  R5, and R2, where the topology is T11, and the primary relay is R15. Similarly, Figures 9 and 10 present characteristic curves and similar fault scenarios for T10, T13, and T14. It is evident that the settings group for each cluster is valid for all topologies belonging to that cluster. This clustering approach vastly reduces the number of adaptive stages while grouping similar topologies. R5, and R2, where the topology is T11, and the primary relay is R15. Similarly, Figures 9 and 10 present characteristic curves and similar fault scenarios for T10, T13, and T14. It is evident that the settings group for each cluster is valid for all topologies belonging to that cluster. This clustering approach vastly reduces the number of adaptive stages while grouping similar topologies.

Conclusions
This paper focused on two principal areas of microgrid protection as the existence of multiple operating modes and the difficulty of relay coordination due to the multi-sourced and multi-loop network architecture. Identified solutions were the use of adaptive protection and using a computational intelligence-based relay setting optimization method. The clustering method was used in handling the different adaptive scenarios. Network element single contingencies were used to define operating topologies. Clustering was used to group similar topologies together in terms of the short circuit current. It proved to be effective in reducing the number of effective protection settings groups needed for the correct operation of the protection system.
A novel hybrid algorithm was used in obtaining the optimized protection settings. WCMFO algorithm confirmed its superiority among the other algorithms, such as GWO, WOA, HHO, PSOGSA, and PSO. As these are some of the better-performing novel algorithms, it is safe to assume WCMFO is better than most of the algorithms. However, there is always a possibility of a betterperforming algorithm with a different hybrid combination. The proposed method has a 50% to 5% decrease in the total relay operating time among the results obtained through other algorithms. Final relay operating times obtained by the settings were well within the initialized constraints. Coordination times did not go below the 0.2 s margin, and the fault operating times did not exceed 1 s for the primary protection. This research presents the successful implementation of adaptive protection using metaheuristics and clustering. This method limited itself to a single type of standard time-current curve, and the selected topologies were based on single element contingencies. Exploring these areas can be identified as future research directions for a flexible protection solution.   It is evident that the settings group for each cluster is valid for all topologies belonging to that cluster. This clustering approach vastly reduces the number of adaptive stages while grouping similar topologies.

Conclusions
This paper focused on two principal areas of microgrid protection as the existence of multiple operating modes and the difficulty of relay coordination due to the multi-sourced and multi-loop network architecture. Identified solutions were the use of adaptive protection and using a computational intelligence-based relay setting optimization method. The clustering method was used in handling the different adaptive scenarios. Network element single contingencies were used to define operating topologies. Clustering was used to group similar topologies together in terms of the short circuit current. It proved to be effective in reducing the number of effective protection settings groups needed for the correct operation of the protection system.
A novel hybrid algorithm was used in obtaining the optimized protection settings. WCMFO algorithm confirmed its superiority among the other algorithms, such as GWO, WOA, HHO, PSOGSA, and PSO. As these are some of the better-performing novel algorithms, it is safe to assume WCMFO is better than most of the algorithms. However, there is always a possibility of a better-performing algorithm with a different hybrid combination. The proposed method has a 50% to 5% decrease in the total relay operating time among the results obtained through other algorithms. Final relay operating times obtained by the settings were well within the initialized constraints. Coordination times did not go below the 0.2 s margin, and the fault operating times did not exceed 1 s for the primary protection. This research presents the successful implementation of adaptive protection using metaheuristics and clustering. This method limited itself to a single type of standard time-current curve, and the selected topologies were based on single element contingencies. Exploring these areas can be identified as future research directions for a flexible protection solution.