Hybrid Harmony Search Algorithm Applied to the Optimal Coordination of Overcurrent Relays in Distribution Networks with Distributed Generation

: In recent years, distributed generation (DG) has become more common in modern distribution networks (DNs). The presence of these small-scale generation units within a DN brings new challenges to protection engineers, since short-circuit currents tend to increase; additionally, as with microgrids, modern DNs may feature several operational modes depending on their topology and the availability of DG. This paper presents a methodology for the optimal coordination of overcurrent relays (OCRs) in modern DNs with a high presence of DG. Given the fact that protection coordination is a non-linear and non-convex optimization problem, a hybrid harmony search and simulated annealing (HS-SA) approach was implemented for its solution and compared against other techniques, such as conventional HS, genetic algorithm (GA), particle swarm optimization (PSO) and hybrid PSO-HS. Several tests were performed on a DN, considering different operative scenarios as a function of the DG available within the network. A comparison with other works reported in the specialized literature was carried out, evidencing the applicability and effectiveness of the HS-SA technique in solving the optimal OCR coordination problem in modern DNs.


Introduction
The growth of distributed generation (DG) units within modern distribution networks (DNs) has opened up new possibilities in the diversification of energy supply, providing a better use of natural resources that contributes to sustainability. DG units are a promising alternative to meet the rising demand for energy, with a beneficial impact on electricity networks that enhances overall efficiency. Nonetheless, the inclusion of DG into traditional DNs leads to bi-directional power flows that alter the traditional load-generation dynamic interaction [1][2][3][4]. This interaction causes the problem of protection coordination to be more complex and challenging than in traditional DNs.
In current DNs, DG units are mainly related to solar and wind power resources; as a consequence, the traditional protection systems developed for DNs must contend with their inherent intermittency. Basically, DNs with the inclusion of DG units have become reconfigurable systems in which the level of current fault is variable; also, the inclusion of DG units converts load nodes into combined generation-load nodes, resulting in bi-directional power flows [5][6][7]. The inclusion of DG units has transformed the traditionally radial DNs into non-radial DNs with the advantage of making the system more flexible, but also with the drawback of the need for more complex protection coordination. Consequently, the coordination problem with the inclusion of DG that guarantees speed, selectivity and reliability has become one of the most challenging tasks when planning the operation of DNs [8].
Overcurrent relays (OCRs) are the preferred protection systems because of their performance, reliability, efficiency, simplicity, and facility of installation. However, their adequate integration and coordination in DNs when DG is incorporated is an emerging subject of research [1,2,9]. The choice and design of a protection coordination system are critical for the proper functioning and reliable operation of DNs [6]. Lately, researchers around the world have applied several optimization techniques to find a suitable protection coordination scheme that is characterized by minimum operation time and which, at the same time, guarantees reliability and safety [10,11].
In [12], the authors presented an impact index to evaluate the coordination of OCRs when DG units are integrated in DNs. The coordination problem of directional OCRs is also discussed in [13] for DNs with the integration of significant DG. The authors of [14] proposed an approach for the coordination of OCRs that involves DNs under unexpected faults. In [15], a multi-objective optimization strategy for the coordination of OCRs in DNs with DG was developed. In [16], the authors propose an online coordination protection system for OCRs that is adapted under topological variations of DNs. This scheme uses intelligent communication devices to obtain real-time information to update relay settings, simulating operation under dynamically changing conditions, and also including loss of loads, generations and lines. In [17], the authors also provided a coordination protection scheme considering topological variations of DNs. Basically, a support vector machine combined with an artificial neural network was used for DN state recognition. Nonetheless, the protection systems that can be adapted under topological variations are usually costly and also feature a complicated implementation [18].
The authors in [19] proposed the concept of non-standard features for the coordination of OCRs in DNs which are not included in IEEE or in IEC standards. Technological developments have provided the platform with the latest generation of digital OCRs that permit the implementation of new protection schemes. The constant expansion of DNs is making the coordination of protections an increasingly complex problem. Non-standard features then emerge as an attractive alternative to enhance the safety and reliability of DNs. In [20], the authors presented a protection coordination strategy considering photovoltaic units with different levels of penetration and locations along the DN; in general terms, they modified the characteristic curves of OCRs by varying the time multiplying setting (TMS) while keeping the pick up current fixed.
Microgrids and modern DNs feature different operational modes, which makes it more complex to achieve optimal settings of OCRs. For overcoming this issue, a constraint that involves the plug setting multiplier (PSM) was presented in [10]. The PSM is defined as the ratio between the fault current measured by the relay and the pickup current. The overcurrent relay tripping characteristic is restricted by the value of the PSM. Standard commercial relay characteristic curves are usually designed in a range in which the minimum value is 1.1 times the PSM, whereas the maximum value is 20 times this value. DG units in DNs significantly magnify the short-circuit level. Hence, for several faults, this defined maximum PSM value could be surpassed, impacting the protection sensitivity and leading to the loss of the coordination protection scheme. The new constraint presented in [10] addresses the above conditions by setting the maximum value of 20 times the PSM. In this case, the optimal coordination problem was solved using a genetic algorithm (GA). The authors in [11] implemented the constraint proposed in [10]; however, they used a non-standard curve to improve the efficiency of the OCR protection. This non-standard curve alters the region of the characteristic curve by considering a value of 100 times the PSM for its maximum value. In [21], the authors considered the upper limit of the PSM as a variable instead of a fixed parameter, as usually done in traditional approaches, and solved the coordination problem by means of a GA.
Smart and adaptive protection schemes are required for modern DNs. In this scenario, user-defined and non-standard features of OCRs should be implemented to meet the new challenges imposed by modern DNs and microgrids; however, they require specialized solution methods. For this reason, this paper proposes a hybrid harmony search and simulated annealing (HS-SA) algorithm for the solution of the coordination problem in modern DNs that integrate DG. Harmony search (HS) is a comparatively easy but efficient evolutionary algorithm. In the HS algorithm, a collection of plausible solutions is generated randomly. A new candidate solution is produced using all available solutions in the harmony memory (instead of only two as in the GA). If this new generated solution is better than the worst current solution in the harmony memory, the worst solution is substituted by this one. In this paper, the hybrid HS-SA algorithm proposed in [22] was implemented. SA introduces diversification in the search process and allows the HS to escape from locally optimal solutions. The efficiency and benefits of HS have been proven in several applications of electrical engineering, including optimal network reconfiguration [23] and harmonic elimination [24]. Additionally, hybrid versions of HS have been proposed in the domains of data mining [25], electricity price forecasting [26], job-shop scheduling [27] and constrained optimization [28], among others. As regards the protection coordination problem, HS was implemented in [29][30][31] to approach the optimal coordination of directional OCRs. The authors in [29,30] present applications for traditional power systems and interconnected power systems, respectively, while in [30], DG is also considered in the coordination problem. However, none of these previous works consider the protection coordination from the point of view of modern DNs. Instead, they solve the problem in power systems which operate at higher voltage levels; also, they do not consider different operative scenarios. Therefore, the main contribution of this paper is the application of a hybrid HS-SA algorithm to solve the optimal OCR coordination problem in modern DNs that include DG and feature several operative scenarios. The new constraint proposed by the authors in [21] was incorporated, showing better results than those reported in [10,11]. Furthermore, a comparison with the GA approach proposed in [21], conventional HS, particle swarm optimization (PSO) and hybrid PSO-HS was also carried out. The results on a benchmark network showed that the implemented HS-SA algorithm outperforms the aforementioned methods, allowing faster coordination times in all operative scenarios under consideration. In this sense, this paper opens a new filed of application for HS and its hybrid versions which consists of the optimal coordination of OCRs in modern DNs and microgrids.

Mathematical Formulation
The mathematical formulation for the optimal coordination of OCRs in distribution networks implemented in this paper is the one proposed in [21].

Objective Function
The objective function of the OCR coordination in electric distribution networks, labeled here as Z, consists of minimizing the time response of any fault. This is expressed in Equation (1) where t r f is the operation time of relay r when fault f takes place, while R and F are the sets of relays and faults, respectively.

Characteristic Curve
In this study, normal inverse characteristic curves are considered for all relays. This type of curve is indicated in Equation (2), where A and B are parameters and TMS r is defined as the time multiplying setting of relay r. Additionally, the ratio between the fault current I f r experienced by the relay and the pickup current ipickup r , labeled as PSM r f , is defined by Equation (3).

Coordination Criterion and Operation Time
The coordination criterion between the acting times of a main and back up relay is given by Equation (4). In this case, CTI stands for coordination time interval, which typically ranges between 0.2 and 0.3 s. t b f is the operation time of the backup relay and t r f is the operation time of the main relay. Equation (5) establishes the time limits of any relay r.

Limits on TSM and Pick Up Currents
Constraints given by (6) and (7) indicate the limits of the time multiplying setting TSM of each relay r and its corresponding pickup current ipickup r . In this case, TMS r is a decision variable for relay r.

Constraints Regarding PSM
The value of PSM rmax must be taken into account when computing the TMS for each relay. This represents the highest current level of the normally inverse curve programmed in the relay before reaching the definite time region. In this case, we consider a PSM rmax variable between values α and β as proposed in [21]. Equations (8) and (9) indicate the limits of PSM r f and PSM rmax , respectively.

Methodology
The model given by Equations (1)-(9) is nonlinear and non-convex; therefore, metaheuristic techniques are the best way to approach it. In this case, several techniques were implemented for comparative purposes, namely particle swarm optimization (PSO) and HS as well as hybrid HS-SA and PSO-HS algorithms.

Particle Swarm Optimization (PSO)
PSO is a stochastic optimization technique used to approach non-convex optimization problems. Within the PSO, the candidate solutions or particles move around the search space according to some rules regarding their position and velocity. Each particle's movement is influenced by its local best known position, and is also guided toward the best known positions in the search space [32]. Equations (10) and (11) are used for computing the velocity and position of a given particle, respectively. In this case, t is the number of the current iteration; x i and v i indicate, respectively, the ith particle's position and velocity; w(t) is the inertia weight; x gBest and x pBest i are the historically best position of the entire swarm and particle i, respectively; and c 1 and c 2 are the personal and global learnidng coefficients, respectively. Finally, r 1 and r 2 are uniformly distributed random numbers in the range [0,1]. Figure 1 depicts the flowchart of the implemented PSO. The algorithm begins with a random generation of candidate solutions or particles. Each particle is randomly located in the search space, bearing in mind the limits of the decision variables. As previously mentioned, every particle is associated with two vectors: one that indicates its position and another one that describes its velocity. These two vectors are updated for each particle in every iteration. In this updating process, the personal best and global best information are considered to update the new positions and velocity of the particles. The process stops once a given number of iterations has elapsed. In [33], the authors combine PSO with HS to approach high dimensional optimization problems. In this paper, a hybrid PSO-HS algorithm is also tested to solve the optimal coordination problem.

Start
Step 1 Read system data and define initial population

End
Step 2

Compute Fitness Function
Are the stopping criteria met? no

Yes
Step 3

Compute Global Best and Personal Best Solutions
Step 4 Update Velocity and Position

Harmony Search
HS is a population-based metaheuristic inspired on the musical process of searching for a perfect state of harmony. It was initially proposed in [34], and despite being relatively new, its effectiveness and applicability have been demonstrated in several areas, including continuous engineering optimization [35], vehicle routing [36], and water distribution networks [37], among others.
In HS, finding a perfect state of harmony is analogous to finding an optimal solution in an optimization problem. Each candidate solution within the HS is called a harmony. In this sense, some analogies can be derived between an optimization process and HS: • Every decision variable is analogous to a musical instrument; • An iteration of the HS is equivalent to the musicians improvising for new harmonies; • The objective function corresponds to the listeners' appreciation; • The optimal or quasi-optimal solution is analogous to the perfect harmony achieved among the musicians of the group.
The initialization of the HS takes place by creating an array called harmony memory (HM) that stores harmonies/individuals randomly generated and evaluated by their respective objective function. The flowchart of the HS algorithm is depicted in Figure 2.

Start
Step 1

HS parameters and optimization problem values initialization
Meet exit condition?

Yes no
Step 2

HM initialization
Step 3

New harmony improvisation
Step 4 Update HM It is better?
Step 5 Select the best harmony in HM no Yes Figure 2. Flowchart of the implemented HS algorithm.

Hybrid HS-SA Algorithm
HS has been successfully hybridized with other metaheuristic techniques to improve its performance. The main objective with hybridization is to combine the strengths of different techniques to find high quality solutions of complex problems. In [25], HS is combined with an artificial electric field algorithm for feature selection in high dimensional data sets. In [28], the authors hybridized a cultural algorithm with HS to solve the constrained optimization problem of diesel blending. In [26], HS was combined with a neural network for forecasting electricity prices. A comprehensive review of the development of the HS algorithm along with hybrid versions and applications can be consulted in [38].
In this paper, the hybrid HS-SA approach proposed in [22] was implemented. SA is a metaheuristic technique inspired by the annealing of certain metals proposed in [39]. This technique starts with a randomly selected solution S and an initial temperature T denoted as T 0 . In every iteration of the algorithm, a new solution S in the neighborhood of S is generated and the temperature is updated. A better solution is always accepted; nonetheless, inferior solutions are also accepted with a probability determined by its fitness and current temperature. Usually, low quality solutions might be accepted at the early stage of the algorithm to enforce diversification of the search, but as the temperature is updated, the probability of accepting inferior solutions decreases.
The implemented HS-SA algorithm is described in Algorithm 1. In this case, HMCR stands for harmony memory considering rate, which is the probability at which a new individual is generated from those already existing in the HM. PAR is the pitch adjusting rate, which is the probability value used for the generation of a new individual. HMS stands for the harmony memory size, which defines the number of individuals stored in the HM, and BW corresponds to the pitch bandwidth. As regards SA, T 0 is the initial temperature and α is the temperature reduction rate. if (∆Cost < 0 or rand < e −∆Cost/T ) then 19 Update HM by replacing the worst harmony by H 20 end if 21 Set T = αT 22 end while 23 Save the best harmony as the solution

Results
The proposed solution technique for the optimal coordination of OCRs was tested with the network depicted in Figure 3 (adapted from [21]), the parameters of which can be consulted in [40]. This system corresponds to an IEC test network used in several research papers [10,11,21,41,42]. Furthermore, for comparative purposes with other techniques, the same formulation of [21] was considered. Table 1 summarizes the four scenarios considered. These scenarios refer to different conditions under which the DN is expected to operate.
In scenario 1, the DG units are disconnected and the load is met using the main power supply; in contrast, in scenario 2, all DG units are connected and the load is fed using not only the DG plants but also the main power supply. In scenario 3, the DG3 and DG4 units are disconnected, so the load is supplied by the main power supply and also by the DG1 and DG2 units. In scenario 4, the DN is disconnected from the main power grid and the load is fed using the DG units only.  Figure 3. Distribution network with DG (adapted from [21]).
Five three-phase line failures were considered. In this case, fault 1 represents a fault on line DL-5; Fault 2 and Fault 3 are failures on lines DL-4 and DL-2, respectively; fault 5 represents a fault on the DL-1 line; and F5 represents a fault on the DL-3 line. The results obtained with the proposed approach were compared with those reported in [10,11]. Additionally, CTI was set to 0.3 s and IEC normal inverse curves were considered for all relays, which are labeled with numbers ranging from 1 to 15. Figure 3 shows the location of each relay. Switches CB-LOOP1 and CB-LOOP2 remained open for all scenarios under consideration.
Several tests were performed to adjust the HS and HS-SA parameters. Such tests were carried out by trial and error based on recommendations reported in [22,34]. The parameters that showed the best performance for both HS and HS-SA algorithms were: number of iterations 1000, harmony memory size 25, number of new harmonies 20, har-mony memory consideration rate 0.9, and pitch adjustment rate 0.1. Furthermore, the SA parameters consist of an initial temperature (T 0 ) of 0.025 and a temperature reduction rate (α) of 0.99, as recommended in [22]. As regards the parameters of the PSO and GA, these were adopted as suggested in [42] and limited to 1000 iterations.
As previously mentioned, PSO, HS, PSO-HS and HS-SA techniques were implemented in order to solve the coordination problem. The results obtained with these techniques are compared with those found in the literature [10,11,21]. Table 2 presents the total operating time (T(s)) obtained with each implemented technique. It is important to highlight that in all scenarios under analysis, HS-SA obtained adequate coordination and lower operating times than reported in [10,11,21]; additionally, HS-SA performed better HS, PSO and PSO-HS. Note that for all operative scenarios, HS, PSO and PSO-HS were also able to obtain a proper coordination with lower operational times than those reported in [10,11].
In conclusion, HS-SA remains an appealing option for solving the optimal coordination problem due to its ease of implementation and high performance.  7 show the convergence curves of the implemented algorithms, as well as the curve of the GA proposed in [21]. In Figure 5, GA presents the fastest convergence; nonetheless, the implemented hybrid HS-SA is the algorithm with the shortest operating time (best objective function). Figure 5 depicts the convergence curves for scenario 2. Note that in this case, HS is the fastest converging algorithm, while HS-SA is the algorithm with the best values of the objective function. In Figure 6, the performance of all algorithms is quite similar, while in Figure 7 it is evident that GA presents faster convergence but slightly lower objective function than HS-SA.

Results for Scenario 1
In scenario 1, the load is supplied from the main power grid and all DG units are offline. The five faults were simulated, the results of which were used as input data for the suggested coordination model. The obtained results are presented in Tables 3 and 4.  Table 3 shows the TMS and PSM imax values for each relay. Table 4 shows the operating times of the relays with the proposed method. Note that the information provided in Table 4 provides the operating times for the main and backup relays. It was observed that, in all cases, HS is able to ensure coordination between the main and backup relays.

Results for Scenario 2
In scenario 2, the load is supplied from the main power grid and all DG units. The coordination results obtained for this scenario are presented in Tables 5 and 6. Table 5 shows the TMS and PSM imax results obtained with the proposed model for each relay. Table 6 shows the operating times of the relays. Note that coordination between main and backup relays is guaranteed.

Results for Scenario 3
In scenario 3, DG3 and DG4 are disconnected, so the load of the distribution network can be supplied through the main power supply as well as through DG1 and DG2. The coordination results obtained for this scenario are presented in Tables 7 and 8. Table 7 shows the TMS and PSM imax results performed with the proposed model for each relay. Table 8 shows the operating times of the relays with the proposed method. In all cases, HS guaranteed coordination between the main and backup relays.

Results for Scenario 4
In scenario 4, the load of the distribution network can be supplied through all DG units. The coordination results obtained for this scenario are presented in Tables 9 and 10. Table 9 shows the TMS and PSM imax results performed with the proposed model for each relay. Table 10 shows the operating times of the relays with the proposed method. In all cases, there is coordination between the main and backup relays.

Conclusions
The complexity of today's distribution networks and the increasing presence of renewable DG require smarter and more adaptive protection schemes. In a context where distribution networks are expected to operate under different scenarios, traditional approaches for OCR coordination may not be reliable for certain distribution network conditions. In this paper, a novel application of HS and HS-SA is proposed regarding the optimal coordination of OCRs in distribution networks with DG. Several tests were performed on a test distribution network, which featured different types of DG units and several operational modes. A comparison of GA, PSO, HS, PSO-HS and HS-SA techniques was carried out. The implemented algorithms proved to be effective for solving the optimal coordination problem, presenting low operating times and quickly reaching convergence; however, HS-SA was the metaheuristic with the best performance in terms of quality of the objective function. The results evidenced the applicability and effectiveness of HS-SA to tackle the optimal coordination problem in modern DNs, achieving high-quality solutions. Other alternatives might be explored regarding HS approaches, such as online applications in real DNs which might require parallelization techniques to reduce computation time, clustering of operational scenarios and exploring other hybrid versions of HS. This opens a new field of applications for the HS technique regarding the safe operation of modern DNs and microgrids.