Calculation of Distance Protection Settings in Mutually Coupled Transmission Lines: A Comparative Analysis

Abstract

Protection of mutually coupled transmission lines poses special challenges to protection engineers, since the general principles outlining the expected performance of the protection elements are based on simplified versions of the issues affecting them. Distance protection elements are among the most commonly accepted methods for protecting mutually coupled transmission lines, however, the calculation of settings and the validation of their performance is a complex task, involving lots of simulations and analysis of results to come up with an acceptable solution. Though various methods have been proposed to deal with the problem of calculating the optimal settings of distance protection elements, there is still room for improvement, since these methods have been formulated for the general case considering simple transmission lines. In this paper, a variation of one of these methods is introduced to enhance the computing times and reliability of the solutions for the specific application of quadrilateral distance protection to mutually coupled transmission lines, through the characterization of the expected performance of the operating characteristic, and the formulation of the optimization problem and the solution method. The obtained results show the improvements in computing times and quality of the solution provided by the proposed algorithm.

1. Introduction

Parallel transmission lines have been extensively used in transmission systems all over the world, given their economic advantages in increasing the power transfer capacity of an existing transmission link between two substations over other alternatives, like the use of series compensation or building a separate transmission line, and taking advantage of almost the same right-of-way, thus reducing the environmental impact. The different operation status of one of the circuits in comparison with its parallel, and the mutual coupling between them, brings a series of challenges to the performance of the protection system.

Mutual coupling between transmission lines is the effect caused by the magnetic field through the flow of current in one of the coupled lines into the others when they are in close proximity to each other, which finds its most common application in parallel transmission lines. Positive and negative sequence mutual impedances are negligible [1], but the zero-sequence mutual impedance is usually significant, and therefore, it takes part in the sequence networks representing ground faults and cannot be ignored.

This paper addresses the problem of determining the optimal reach settings for the distance protection quadrilateral characteristic applied to mutually coupled lines, which is one of the most commonly used methods to protect these elements [1]. The optimization problem raises from the difficulty related to the variation of the apparent impedance, as calculated by the distance relay, under the different operating status of the mutually-coupled transmission lines, which includes three main scenarios: (1) both transmission lines are in service, (2) the mutually-coupled transmission line is out of service and grounded at both ends, and (3) the mutually-coupled transmission line is out of service and isolated from the ground at both line ends.

As indicated in [2], it will not be possible to ensure proper performance of the distance protection element for all the possible operating conditions of the mutually-coupled transmission lines, especially when the typical setting criteria is applied for the calculation of the reach settings and the earth compensation factor is calculated according to Equation (1), since such equation is obtained from the assumption that the protected transmission line has no mutual coupling [3], and therefore, the effects of the zero-sequence mutual impedance are not considered:

$$k_0=\frac{Z_0-Z_1}{3Z_1},$$

(1)

The different alternatives for coping with these issues are described in [1,2,4,5], and they consist of the use of different earth compensation factors for the different protection zones of the distance protection element, for both applications with a single settings group or several, or in the variation of the reach settings related to these zones, in order to avoid overreach (loss of selectivity) or underreach (loss of sensitivity). A third option, as described in [2,6] is to compensate the zero-sequence current induced by the mutually coupled circuit by directly measuring it, and taking it to the relay in the protected line for the calculation of a second compensation factor, called the mutual coupling compensation factor; however, this alternative has limited use, due to reliability issues in the protection system related to the addition of a failure mode, additional wiring and cost, and because this alternative is effective only in transmission lines with the same infeed at both ends, as described in [2].

The alternative involving the use of different earth compensation factors to different protection zones is usually applied to prioritize selectivity for underreaching zones (Zone 1 in most applications) and sensitivity for overreaching zones, since such zones work as backup protection elements for other protection schemes in the proximity of the protected line. The main advantage of this method is that all modern numerical protection relays include the option to provide different earth compensation factor settings for the different protection zones. Its main disadvantage is the complexity involved in the calculation of the different earth compensation factors, which involves certain assumptions not representing all the different conditions to which the distance protection will be subjected to, as described in [2,4,5].

Finally, the alternative involving the setting of the zone reaches for the different operating conditions of the mutually-coupled transmission lines will be the method used for the optimization problem to be described and resolved in this paper. The main advantage of such approach lies in the fact that the reach settings are adapted according to the actual performance of the transmission line, based on a single earth compensation factor, calculated from the normal operating conditions of the mutually-coupled transmission lines, according to the method described in [2,4], which is clearly the most probable operating scenario. The main disadvantage of this method is that not all distance protection relays currently available have the possibility to set different reaches only applicable for the phase-to-ground impedance loops and so, the phase-to-phase loops and related operation characteristic settings will also be affected by the settings required to cope with a problem only affecting the phase-to-ground impedance loops; however, most state of the art numeric relays include different reach settings for the phase and ground operation characteristics.

The main contribution of this paper, when compared to previous works used as reference [2,7,8,9,10,11,12], is the simplification of some of the premises of the optimization problem formulation, based on a comparative analysis with the techniques currently employed to determine the settings of the distance protection relays, according to its application to mutually-coupled transmission lines, which involve a series of additional problems and special considerations, not included in the reference works. On the other hand, the optimization techniques described in previous works [7,8,9,10,11,12] do not include details on the methodology to set its parameters according to the natural bounds and application limits of distance protection relays to transmission lines, which causes higher computing times and the addition of data that should not affect the actual performance of the distance protection relays. The solution method proposed in [10] uses the Monte Carlo method to estimate system model parameters, thus improving the accuracy of the calculations related to the impedance loops; however, its main disadvantage is related to the excessive number of fault conditions required to provide enough reliability to the proposed settings, with the high computation cost related. Other works [11,12] deal with the optimization of the Zone 1 reach settings, in order to improve sensitivity over the traditional setting criteria based on engineering experience, as described in [1]; however, the authors of [11] assume the voltage phasors of the power sources of the equivalent fault circuit to be known parameters, which is very difficult to accomplish by using the measured values available for the actual distance protection relay. The method presented in [12] is limited in both scope and application, since it only improves the sensitivity of Zone 1 reach, and it is formulated for series compensated transmission lines. Finally, the approaches proposed in [7,8,9], based on probability distribution functions related to the fault conditions simulated have the inconvenient that these functions are difficult to estimate; however, for this particular case, the proposed method will eliminate some of the probability distribution functions according to performance criteria expected from the distance relay under specific fault conditions, as otherwise, the protection system could perform undesirably for events with low probability of occurrence, included in the preference objective function.

The proposed method to solve the optimization problem related to the determination of the reach settings of the distance protection is thus based on the formulation outlined in [7], improving the technical basis for the formulation of the optimization problem, which will include a short-circuit sensitivity analysis for establishing the actual performance of the relay to be set under the multiple fault and system conditions affecting the measured impedance loops, and including the probability of occurrence of every single fault simulated in an specialized software tool, in order to determine the preference function related to the optimization problem to be solved. In addition, the solution method will be flexible enough to accommodate to special performance requirements of the distance protection elements, which will depend on the evaluation of the fault conditions obtained from the short-circuit sensitivity analysis.

The methodology and proposed optimization algorithm (genetic algorithm) are entirely based on the main factors that actually affect the performance of these protection schemes, therefore, the results obtained for the distance protection reach settings have a wider and stronger base for ensuring its proper performance in terms of reliability when compared to previous works [7,8,9,10,11,12]. The proposed optimization algorithm, successfully used in [13] in comparison with conventional calculation methods, will prove to be a compromise between the in-depth short-circuit sensitivity analysis, obtained through simulation, and the computing time required to obtain a good solution.

2. Methodology

Distance protection elements based on quadrilateral operation characteristics, as the one shown in Figure 1, are the most widely used in the protection schemes of transmission lines, given their good fault resistance coverage and higher immunity to load when compared to other operating characteristics, according to [1], which also discusses the applications and limitations of other commonly used characteristics. The reach settings of every zone of the quadrilateral characteristic are composed of a resistive setting and a reactive setting, which are commonly calculated according to widely known criteria, presenting very minor variations among the technical literature and the practices used by utility companies all over the world. The main principles and applications of these setting criteria are outlined in [1,5,6]. Of course, the application of these setting criteria does not ensure proper performance of the distance protection elements for all the possible fault conditions, being subjected to different phenomena affecting mainly selectivity and sensitivity, like the infeed effect, which can affect both, and varies according to the system pre-fault conditions and fault resistance [1,5,6].

2.1. Expected General Performance of the Distance Protection Zones

The reach settings calculated for every zone of a quadrilateral operation characteristic, as that shown in Figure 1, must be based on a trade-off between selectivity and sensitivity. Selectivity is the property of any protection system for which all its protection devices operate in a coordinated fashion, so that the outage area can be reduced when a faulted component of the power system is isolated from the system, thus providing a certain degree of security to the protection system. Sensitivity refers to the minimum operating quantities required by the protection element to detect a fault or abnormal operating condition in the protected component.

Zone 1 of the distance protection element shown in Figure 1 is usually called “underreach zone”, since it is usually set to operate with no time delay when a measured impedance loop is inside this zone. Due to the combination of measurement errors from the current and voltage transformers [14,15,16], the measured impedance loops can have a tolerance as high as ±15% of the actual impedance from the relay location to the line section where a fault occurs. As a result, Zone 1 reactive reach setting is usually set below 85% of the positive sequence impedance of the protected line [1], so that faults in adjacent elements (transformers, compensation devices or other lines) near the remote end of the protected line can be tripped by the protection elements of such elements, thus ensuring a certain degree of selectivity.

Since Zone 1 cannot cover all the faults occurring in the entire protected line, Zone 2 is aimed to cover all these faults. So, for the same reason exposed previously, related to the maximum deviation in the measured impedance loops, a minimum reactive setting of 120% of the positive sequence impedance of the protected line is required to fulfill the condition of covering all the faults applied in the protected line. The former is the minimum sensitivity reach of the distance protection Zone 2 characteristic. The upper limits for the reactive reach setting of Zone 2 are aimed to avoid selectivity issues with adjacent protection elements. The first limit is intended to avoid overlapping with another distance protection element in Zone 2, since this would require to coordinate both elements by reducing the speed of the one with higher reach, as presented in Figure 2. In order to avoid this condition, a commonly applied criterion is to limit the Zone 2 reactive reach to the sum of the positive sequence impedance of the protected line and 50% of the positive sequence impedance of the shortest adjacent line [1]. The second limit is related to a general constraint of the distance protection elements, for which such elements must not operate for faults occurring in voltage levels different from that of the protected line, due to the difficulty to coordinate the definite-time operating characteristic of the distance protection elements with inverse time overcurrent elements, given the reduced speed of the later for high impedance faults. As a result, the reactive reach of Zone 2 is usually limited to the sum of the positive sequence impedance of the protected line and 80% of the equivalent positive sequence impedance of the transformers connected to the remote substation. The Zone 2 time delay is usually set from 250 ms to 400 ms.

Zone 3 is usually set to provide time-delayed remote backup protection in case the main protection of the remote substation fails to operate; therefore, it is typically set to 120% of the impedance measured by the relay for a fault at the remote end of the adjacent line with highest impedance [1]. Due to the same limitation described for Zone 2, for which distance protection elements must not detect faults in other voltage levels, the reach must not be higher than 80% of the equivalent positive sequence impedance of the transformers connected to the remote substation.

Zone 4 is a reverse zone, with the highest time delay amongst all the protection zones (over 1.5 s), and serves both as a delayed backup for faults near the protected line at the local substation, and as a blocking element to the communication assisted trip schemes. Its setting is given depending on the application, and it usually does not require further performance evaluation, given its high delay time and short coverage.

Finally, regarding the resistive reach setting, it must ensure that the distance protection does not trip for heavy power flow, and therefore, the reach can be as high as that corresponding to the real part of the minimum load impedance seen by the distance protection under such conditions, as presented in Equation (2), considering the maximum expected power factor angle. Equation (3) can be used to calculate the minimum load resistance based on the former constraints, assuming that the minimum load impedance seen by the distance protection element is related to extreme loading conditions, which according to the NERC (North American Electric Reliability Corporation) guidelines described in [1], can be obtained considering the minimum voltage, $V_{L-L,min}$, and the maximum current flow, $I_{max,load}$, during extreme loading conditions, and considering the composite error of the measured impedance from the combination of measurement errors of the current and voltage transformers [14,15,16], ${\epsilon }_z$, and the maximum power factor angle, ${\phi }_{pf}$, expected through the protected line:

$$Z_{min,load}=\frac{V_{L-L,min}}{\sqrt{3}·I_{max,load}}$$

(2)

$$R_{min,load}={\epsilon }_z·cos\left({\phi }_{pf}\right)·Z_{min,load}$$

(3)

All the terms in Equation (3) depend on the application, country regulations and utility company’s operation philosophies, therefore, the minimum load impedance, used to limit the resistive setting reach of all the zones of the operation characteristic, is a parameter for the calculation of settings of the distance protection element.

In some cases, the resistive reach of a particular zone, not complying with the previously mentioned performance criteria, especially the selectivity criteria, must be reduced in order to shape the operating characteristic of the distance protection element to the specific requirements of a given application.

2.2. Performance of the Distance Protection for Mutually Coupled Transmission Lines

For the particular application of quadrilateral distance protection for mutually coupled transmission lines, the variation in the impedance seen by the distance relay for the different operating conditions of the protected lines makes more difficult estimating and validating the reach settings of the distance protection.

Figure 3 presents the different operating conditions of mutually coupled transmission lines. There must be emphasized that there is no need to share terminals for the mutually coupled lines, since the mutual coupling would still represent the same issues for the application of distance protection. A detailed explanation of the model used to represent the operating conditions of mutually coupled transmission lines is presented in [17,18], which correspond to the technical references of DIgSILENT PowerFactory, the software used to carry out the short-circuit sensitivity analysis and the calculation of the impedance seen by the distance protection relay model used.

According to [1], the most recommended alternatives to cope with the problems related to mutually coupled transmission lines involve the modification of the zone reach settings or the residual compensation factor, so that the operating characteristic of the distance protection element can have an adequate balance between selectivity and sensitivity. Modern numerical protection relays include the option to use different earth compensation factors for the specific selectivity and sensitivity requirements of the operating characteristic zones.

In terms of reliability, the method used to set the distance protection element should result in the use of a single setting group, which will perform correctly for all earth faults under the different operating conditions presented in Figure 3. Addressing the problem of determining the earth compensation factors for every zone by using currently available algorithms has the difficulty of stablishing bounds for the magnitude and angle of the earth compensation factors, since this phasor cannot be directly used to address the performance of the distance protection element for a given set of fault conditions.

The main advantage of modifying the reach settings related to the zones of the operating characteristic of a given quadrilateral distance protection element is that the reach settings can be used to evaluate the performance of every zone for a given set of faults for which the measured impedance will not change its value, as it would if the earth compensation factor is to be modified. Therefore, a set of optimization and metaheuristic algorithms can be applied to the set of faults resulting from the simulation of relevant fault conditions required to address the performance of each zone of the operating characteristic, as previously described, and stablish the optimal zone reach settings.

The issue still remaining before formulating the optimization problem for determining the optimal zone reach settings is the correct earth compensation factor to be used, given that the general formulation for obtaining it, as presented in Equation (1), will affect the accuracy of the distance protection element when it comes to the calculation of the phase-to-ground impedance loops, since this corresponds to the operating condition presented in Figure 3c, for which the effects of the zero-sequence mutual impedance are not considered, thus affecting the performance of the zone for which the optimal reach settings are to be determined. The proposed solution for this inconvenience is to use the earth compensation factor calculated from the normal operating condition of the mutually coupled transmission lines, represented in Figure 3a, which correspond to the operating condition with the highest probability of occurrence, given the fact that the other two conditions do not have a duration higher than a few minutes, for line outages caused by permanent faults, as represented by the condition shown in Figure 3c; or higher than a few hours or even a few days, for maintenance works on one of the mutually coupled transmission line, as represented by the operating conditions shown in Figure 3b,c.

In order to calculate the earth compensation factor under normal operating conditions, a single line to ground fault is to be simulated at the remote bus of the protected line, with zero fault resistance, using the condition shown in Figure 3a, for which the distance protection element should ideally measure a phase-to-ground impedance, $Z_{\phi -G}$, equivalent to the positive sequence impedance of the protected line, according to Equation (4). Different equations proposed for the calculation of the impedance at the distance relay location, regarding the different operating status of the mutually coupled transmission lines, are described in [2]; however, due to the limitations of the classic approaches to deal with the performance of distance protection elements in these applications, as described in [1], the proposed solution will not be based on variations to the calculated impedance used by most modern numerical protection relays including the distance protection function, as presented in Equation (4):

$$Z_{\phi -G}=\frac{V_{\phi -G}}{I_{\phi }+k_0·I_N}$$

(4)

Since the value of $Z_{\phi -G}$ corresponds to the positive sequence impedance of the protected line, and the phase voltage $V_{\phi -G}$, phase current $I_{\phi }$ and residual current $I_N$, are available from the simulated fault, the earth compensation factor $k_0$ can be obtained from Equation (5):

$$k_0=\frac{1}{I_N}·\left(\frac{V_{\phi -G}}{Z_{\phi -G}}-I_{\phi }\right)$$

(5)

Equation (5) is the equivalent of Equation (1) for the case of mutually coupled transmission lines, since it is obtained using the same premises applied to define the later, as detailed in [3]. It is important to consider that the measured impedance loops of the proposed method are to be obtained directly from the advanced protection relay models available in DIgSILENT PowerFactory, and the method can be applicable to any other given software including these models, thus making its application highly flexible, since it will not consider the specific case of any given distance protection element including variations to the way the phase-to-ground impedance loops are calculated, according to Equation (1). With a general earth compensation factor defined, what remains is to formulate the optimization problem for the calculation of the zone reach settings.

2.3. Formulation of the Optimization Problem to Calculate the Zone Reach Settings

As previously mentioned, the formulation of the optimization problem to calculate the zone reach settings will follow the approach outlined in [7], with some modifications required to increase reliability over the obtained solution, based on practical considerations for the performance of every zone of the quadrilateral distance protection element.

The optimization problem will be oriented to obtain a trade-off between the property of being selective and the property of being sensitive, for every zone, based on representative fault locations, organized in a hierarchical fashion for every zone, and entirely based on the results from a detailed short-circuit sensitivity analysis. For example, Figure 4a shows the operating characteristic of Zone 1, defined using the setting criteria previously described, and faults expected to be covered by Zone 1, this is, faults simulated at a maximum of 85% of the protected line length, given the performance criteria previously described for this zone; therefore, all those faults outside Zone 1 represent a loss of sensitivity.

Figure 4b shows the same Zone 1 operating characteristic as in Figure 4a, and faults expected to be outside the reach of Zone 1, this is, faults simulated beyond 85% of the protected line length, or in adjacent lines. Accordance to the performance criteria discussed previously, all the impedances entering in Zone 1 represents a loss of selectivity. This approach for addressing the problem based on the selectivity and sensitivity performance will be applied to every zone of the forward-looking quadrilateral distance protection element.

The decision variables of the problem are the zone reaches (R and X couples), with the solution being obtained per zone. The objective function will be a preference function ($M$), corresponding to the sum of the probability of presenting loss of selectivity $p\left(T\right)$, and loss of sensitivity $p\left(S\right)$, as presented in Equation (6), where $k$ represents the weight assigned to the performance properties of selectivity and sensitivity in the objective function. Therefore, the optimization problem with the preference function proposed in Equation (5) will correspond to a minimization type:

$$Minimize M=k·p\left(T\right)+\left(1-k\right)·p\left(S\right)$$

(6)

Figure 5 shows the sets representing the distribution of faults included in the short-circuit sensitivity analysis, and its classification according to the factors affecting the performance of the reach settings for Zone 1, which in turn can be generalized for all the other zones. G1 represents the set of calculated impedances obtained from the short-circuit sensitivity analysis, and which are relevant for the calculation of Zone 1 reach settings, this is, faults along the protected line and adjacent lines connected to the remote bus. F1 is a subset of G1, and stores all the calculated impedances, which are expected to be outside of Zone 1 reach (beyond 85% of the protected line). D1 is a subset of G1, and stores all the calculated impedances which are expected to be within Zone 1 reach (below 85% of the protected line).

As the impedances seen by the distance protection element are obtained directly from the distance protection relay model included in the network model under analysis, all the simulated fault conditions which do not pass the starting and direction determination blocks in the distance protection element are discarded, which will improve the computing times of the proposed algorithm, since this deletes all the irrelevant simulation cases for determining the settings of the operation characteristic.

For every apparent impedance stored in subsets D1 and F1, a probability of occurrence is calculated, considering only the following categories, and corresponding distribution:

(1)

Fault type p(Tf), for which the following probability distribution is assumed [19]: single-phase to ground faults (80%), double-phase to ground faults (12%) and three-phase faults (8%).

(2)

Operative status of the mutually-coupled circuit p(Eo), assuming an empirical distribution: both circuits in service (90%), coupled circuit out of service and grounded at both ends (6%), coupled circuit out of service and isolated from ground (4%).

According to this, the probability related to every apparent impedance ($Z_{ap}$) obtained from the short-circuit sensitivity analysis is calculated by using Equation (7), which includes the probability functions depending on the factors which affect the most the performance of the zone being calculated:

$$p\left(Z_{ap}\right)=p\left(T_f\right)·p\left(E_o\right)$$

(7)

As it can be observed in [7], other conditions of the simulated faults were considered within the set of probability functions, like the pre-fault load flow, which should not be considered to give weight to the apparent impedance to be part of the objective function, since the distribution of generation and demand varies during all day, and the distance protection element should perform accurately for all of those conditions, as they will be exposed to these conditions always. Within the short-circuit sensitivity study, only the most variable pre-fault load flow conditions should be considered, since distance relays only present major variations in their operation when subjected to extreme conditions, specially directionality and starting. The same reasoning can be applied to the probability function related to the faulted element, since they are not based on actual acceptable data, which is difficult to get, and again, the distance protection element should perform correctly, no matter the faulted element (line, bus, transformer, etc.), which at the same time, has no deep effect on the expected performance of the zone being set, as discussed before. Regarding the probability function related to the fault resistance, as described in [7], it has no basis, since this is pretty much related to the fault location and the variations to which the fault resistance is subjected to, due to the ever changing conditions of the power system and those subjected to the conditions of the soil and environment. As a result, the distance protection must ensure proper performance for all the fault resistances covered by the operation characteristic, and therefore, this should not be included as a probability function. Finally, regarding the probability function related to the measurement error of the impedance, this was already included in the performance criteria for every zone, which plays a key role in the adjustment of the solution method to obtain feasible solutions.

For every solution obtained for the zone reach settings, the apparent impedances ($Z_{ap}$) will be stored in subsets T1 ($T1\in F1$) and S1 ($S1\in D1$), according to the comparison of every $Z_{ap}$ with the zone reach settings of every solution available. The value of $p\left(T1\right)$ will be obtained as the sum of all the $p\left(Z_{ap}\right)$ related to the apparent impedances stored in T1, whereas the value of $p\left(S1\right)$ will be obtained as the sum of all the $p\left(Z_{ap}\right)$ related to the apparent impedances stored in S1. With the calculated probability functions, the objective function for zone 1 is described in Equation (8):

$$M1=k·p\left(T1\right)+\left(1-k\right)·p\left(S1\right)$$

(8)

The constraints of the formulated optimization problem correspond to the limits described in the performance criteria for every zone, as discussed earlier. Therefore, the constraints of the optimization problem involving the calculation of the reach settings for Zone 1, specifically the resistive reach (R1) and reactive reach (X1), are obtained from Equations (9) an (10), which are in turn based on the performance criteria previously stablished for Zone 1.

$$0.2·XL1\le R1\le R_{min,load}$$

(9)

$$0.2·XL1\le X1\le 0.85·XL1$$

(10)

Constant XL1 is the positive sequence reactance of the protected transmission line, and $R_{min,load}$ is the minimum load resistance, calculated from Equation (3). A similar formulation will be followed for Zone 2 and Zone 3.

2.4. Sample Transmission System and Methodology for Obtaining the Input Data of the Solution Strategy

As mentioned early, the typical setting criteria, defined in Section 2.1, will be used to define the initial reach settings of the zones of the quadrilateral distance protection element. Next, by means of a short-circuit sensitivity analysis, performed through the simulation of the given network model and protection relay in DIgSILENT PowerFactory software (2018 version), a set of considerable apparent impedances will be obtained, as seen by the distance protection relay to be set. The short-circuit sensitivity analysis will cover multiple types, locations and conditions of the simulated faults in the sample transmission system, in order to characterize the required performance of the distance protection element for a given application.

Three fault types will be considered: single-phase to ground, double-phase to ground and three-phase, in order to capture the infeed effect on the calculated impedance loops at the relay location, considering the influence of the sequence quantities on the performance of the operating characteristic. The fault impedance simulated will be purely resistive, as commonly assumed in protection coordination studies, and given the actual nature of this impedance, predominantly resistive, as seen by the distance protection relay. The set of values for the fault resistance will be user defined, as indicated in the pseudocode of

Table A1. Detailed pseudocode for the solution method algorithm.

Solution Method Algorithm
1:	Setting of the parameters related to the algorithm
2:	Variable declaration:
3:	Sets: auxiliary sets for faulted elements: sLines, sBars
4:	Integers: auxiliary counters, vector sizes, set and matrices sizes, integer operators of the GA
5:	Objects: short-circuit command, distance relay measured impedance, faulted line, faulted bar
6:	Doubles: auxiliary operators for probability counters, reach settings bounds
7:	Definition of the fault impedance vector
8:	Definition of the fault location vector (for faulted lines only)
9:	Edition of matrices G1, D1 and F1
10:	Start row counters for G1, D1 and F1
11:	Part 1: simulation and classification of faults
12:	Initialize counter for status of the mutually coupled line
13:	while status of the mutually coupled line <= 3:
14:	if status = both lines are in service:
15:	Energize mutually coupled line
16:	iStatus = 1
17:	else if status = mutually coupled line out of service and grounded at both line ends:
18:	Isolate and ground mutually coupled line
19:	iStatus = 2
20:	else if status = mutually coupled line out of service and ungrounded at both line ends:
21:	Isolate mutually coupled line
22:	iStatus = 3
23:	Initialize counter for fault type
24:	while fault type <= 3:
25:	if fault type = single-phase to ground:
26:	Set fault type in the short-circuit command to single-phase to ground
27:	iFaultType = 1
28:	else if fault type = double-phase to ground:
29:	Set fault type in the short-circuit command to double-phase to ground
30:	iFaultType = 2
31:	else if fault type = three-phase:
32:	Set fault type in the short-circuit command to three-phase
33:	iFaultType = 3
34:	Set the short-circuit method to be used
35:	Set the faulted object (line or bus bar)
36:	Initialize counter for fault location
37:	while Fault location available at the fault location vector:
38:	Set the fault location according to the value of the counter for fault location and its corresponding position in the fault location vector
39:	Initialize counter for fault impedance
40:	while Fault impedance available at the fault impedance vector:
41:	Set the fault impedance according to the value of the counter for fault impedance and its corresponding position in the fault location vector
42:	Execute short-circuit simulation
43:	if Distance protection starts, detects the fault in the forward direction and the resistive component of the apparent impedance measured from the fault is lower than the resistive component related to the minimum load transfer impedance:
44:	if fault type = single-phase to ground or three-phase:
45:	Store the apparent impedance of loop A to ground in set G1
46:	else if fault type = double-phase to ground:
47:	Store the apparent impedance of loop B-C in set G1
48:	if fault location is inside the protected line and until 85% of the total distance from the relay location:
49:	if fault type = single-phase to ground or three-phase:
50:	Store the apparent impedance of loop A to ground in set D1
51:	Calculate p(Zap) according to fault type and status of the mutually coupled line
52:	else if fault type = double-phase to ground:
53:	Store the apparent impedance of loop B-C in set D1
54:	Calculate p(Zap) according to fault type and status of the mutually coupled line
55:	else if fault location is outside the protected line or inside and above 85% of the total distance from the relay location:
56:	if fault type = single-phase to ground or three-phase:
57:	Store the apparent impedance of loop A to ground in set D1
58:	Calculate p(Zap) according to fault type and status of the mutually coupled line
59:	else if fault type = double-phase to ground:
60:	Store the apparent impedance of loop B-C in set D1
61:	Calculate p(Zap) according to fault type and status of the mutually coupled line
62:	Part 2: Genetic Algorithm (GA) implementation
63:	Part 2.1: Stablishing the initial population
64:	Creation of the initial population matrix
65:	Initialize auxiliary counter for the population matrix
66:	Definition of the upper and lower limits for the reaches (R and X) of Zone 1:
67:	RLowerLimit = 0.2*XL1 (XL1: positive sequence reactance of the protected line)
68:	XLowerLimit = 0.2*XL1
69:	RUpperLimit = 0.45*ZminLoad
70:	XupperLimit = 0.85*XL1
71:	Calculation of steps in R and X for the calculation of the initial population
72:	DeltaR = (RUpperLimit − RLowerLimit)/25
73:	DeltaX = (XUpperLimit − XLowerLimit)/25
74:	while Population <= Population Size:
75:	Creation of the initial population individuals based on an uniform distribution of solutions between the upper and lower limits for the reaches of Zone 1 (Rreach, Xreach)
76:	while Population <= Population Size:
77:	Edition of matrices S1 from D1 and T1 from F1
78:	Start row counters for S1 and T1
79:	Initialize probability accumulators p(S1) and p(T1) at zero
80:	for all (Rap,Xap) entries in D1:
81:	if {(Rap ∈ D1) > Rreach} or {(Xap ∈ D1) > Xreach}:
82:	Store Rap, Xap and p(Zap) in S1
83:	p(S1) = p(S1) + p(Zap)
84:	for all (Rap,Xap) entries in F1:
85:	if {(Rap ∈ F1) < Rreach} and {(Xap ∈ F1) < Xreach}:
86:	Store Rap, Xap and p(Zap) in T1
87:	p(T1) = p(T1) + p(Zap)
88:	Calculate the preference function for the current solution in the initial population as:
89:	M1 = k*p(T1) + (1 − k)*p(S1)
90:	Part 2.2: Natural Selection
91:	Creation of the new generation population matrix
92:	Initialize auxiliary counter for the new generation population matrix
93:	Initialize the new generation population matrix with the initial population matrix values
94:	Initialize counter of the number of iterations for the GA
95:	while below the maximum number of iterations:
96:	Sort next generation population matrix according to the objective function
97:	Natural Selection: select the 25 individuals with the lowest value of the preference function (selection rate of 2/3)
98:	Creation and initialization of the auxiliary population matrix
99:	Store the 25 individuals from the natural selection in the auxiliary population matrix
100:	Part 2.3: Mating
101:	Creation and initialization of the crossover matrix
102:	Initialize counter for number of mates
103:	while Number of Mates <= 25:
104:	Random selection of three candidates to father from the subset of the natural selection
105:	Selection of the father among the three candidates by tournament (candidate with the lowest objective function)
106:	Store selected father in the crossover matrix
107:	Random selection of three candidates to mother from the subset of the natural selection
108:	Selection of the mother among the three candidates by tournament (candidate with the lowest objective function)
109:	if the selected mother is the same as the selected father:
110:	Change the mother by another candidate among the three in the tournament
111:	Store selected mother in the crossover matrix
112:	Part 2.4: Crossover
113:	Initialize counters for number of crossings and number of children
114:	while Number of Crossings <= 25:
115:	Get the mother and father from the crossover matrix
116:	Get the values for R and X from both the mother and father stored in the auxiliary population matrix (the 25 chosen from the natural selection)
117:	if Crossover Rate = 0.5:
118:	Random selection of the variable (R or X) to be crossed
119:	if the selected value to be crossed is R:
120:	β = Random{0,1}, with β ∈ ℜ (real number)
121:	α = Random{0,1}, with α ∈ Ζ (integer number)
122:	Perform crossover for the first child:
123:	Rchild1 = β*Rfather + (1 − β)*Rmother
124:	Xchild1 = α*Xfather +(1 − α)*Xmother
125:	Store the first child in the auxiliary population matrix
126:	Increase counter for number of children by 1
127:	β = Random{0,1}, with β ∈ ℜ (real number)
128:	α = Random{0,1}, with α ∈ Ζ (integer number)
129:	Perform crossover for the second child:
130:	Rchild2 = β*Rfather + (1 − β)*Rmother
131:	Xchild2 = α*Xfather +(1 − α)*Xmother
132:	Store the second child in the auxiliary population matrix
133:	Increase counter for number of children by 1
134:	else if the selected value to be crossed is X:
135:	β = Random{0,1}, with β ∈ ℜ (real number)
136:	α = Random{0,1}, with α ∈ Ζ (integer number)
137:	Perform crossover for the first child:
138:	Rchild1 = α*Rfather +(1 − α)*Rmother
139:	Xchild1 = β*Xfather + (1 − β)*Xmother
140:	Store the first child in the auxiliary population matrix
141:	Increase counter for number of children by 1
142:	β = Random{0,1}, with β ∈ ℜ (real number)
143:	α = Random{0,1}, with α ∈ Ζ (integer number)
144:	Perform crossover for the second child:
145:	Rchild2 = α*Rfather +(1 − α)*Rmother
146:	Xchild2 = β*Xfather + (1 − β)*Xmother
147:	Store the second child in the auxiliary population matrix
148:	Increase counter for number of children by 1
149:	else if Crossover Rate = 1:
150:	Random selection of the variable (R or X) to be crossed
151:	β = Random{0,1}, with β ∈ ℜ (real number)
152:	Perform crossover for the first child:
153:	Rchild1 = β*Rfather + (1 − β)*Rmother
154:	Xchild1 = β*Xfather + (1 − β)*Xmother
155:	Store the first child in the auxiliary population matrix
156:	Increase counter for number of children by 1
157:	β = Random{0,1}, with β ∈ ℜ (real number)
158:	Perform crossover for the second child:
159:	Rchild2 = β*Rfather + (1 − β)*Rmother
160:	Xchild2 = β*Xfather + (1 − β)*Xmother
161:	Store the second child in the auxiliary population matrix
162:	Increase counter for number of children by 1
163:	Increase counter for number of crossings by 1
164:	Part 2.5: Mutation
165:	Set the auxiliary population matrix row counter at 26, so that mutation is only applied to the 50 children obtained from the crossover
166:	while Population matrix row counter <= Population Size (75):
167:	Edition of matrices S1 from D1 and T1 from F1
168:	Start row counters for S1 and T1
169:	Initialize probability accumulators p(S1) and p(T1) at zero
170:	for all (Rap,Xap) entries in D1:
171:	if {(Rap ∈ D1) > Rreach} or {(Xap ∈ D1) > Xreach}:
172:	Store Rap, Xap and p(Zap) in S1
173:	p(S1) = p(S1) + p(Zap)
174:	for all (Rap,Xap) entries in F1:
175:	if {(Rap ∈ F1) < Rreach} and {(Xap ∈ F1) < Xreach}:
176:	Store Rap, Xap and p(Zap) in T1
177:	p(T1) = p(T1) + p(Zap)
178:	Calculate the preference function for the current solution in the auxiliary population as:
179:	M1 = k*p(T1) + (1-k)*p(S1)
180:	Update objective function in the auxiliary population matrix
181:	Sort the auxiliary population matrix according to the value of the objective function
182:	Calculate the Number of Mutations:
183:	#Mutarions = MutationRate*(PopulationSize − ElitePopulation)*#Variables
184:	Initialize mutation counter
185:	while mutation counter <= #Mutations:
186:	Random selection of the cell of the auxiliary population matrix to be mutated, avoiding to include in the mutation operator those individuals in the elite population:
187:	iAuxRow = Random(ElitePopulation + 1, PopulationSize), with iAux Row ∈ Ζ (integer)
188:	iAuxColumn = Random(1,2), with iAuxColumn ∈ Ζ (integer)
189:	if iAuxColumn = 1:
190:	dAux = Random(RLowerLimit,RUpperLimit), with dAux ∈ ℜ (real number)
191:	Change the cell of the auxiliary population matrix indicated by iAuxRow and iAuxColumn to the value obtained for dAux
192:	else:
193:	dAux = Random(XLowerLimit,XUpperLimit), with dAux ∈ ℜ (real number)
194:	Change the cell of the auxiliary population matrix indicated by iAuxRow and iAuxColumn to the value obtained for dAux
195:	Part 2.6: Setting of the new generation population matrix
196:	Assign all the final values in the auxiliary population matrix to the new population matrix
197:	Increase by 1 the counter of the number of iterations for the GA
198:	Part 3: Local Search implementation
199:	Take the best solution from the final population of the GA
200:	Define neighborhoods for R and X with steps according to the minimum combined measurement error of the instrument transformers
201:	Calculation of the objective function for all the solutions in the neighborhood
202:	Select the best solution for changes in R and X reaches (separately)
203:	if the solution obtained improves the best solution from the final population of the GA
204:	Repeat until it is not possible to improve the previous solution
205:	Report the final solution

, in order to adapt the algorithm to the specifics of a given application, and considering the maximum resistive reach setting, as described in Section 2.1. This is important, because fault impedances beyond a certain value will necessarily represent a loss of sensitivity for the algorithm, but since they will lie beyond the maximum practical resistive reach setting, it makes no sense including such high values; therefore, it is important to properly set the maximum fault impedance according to the expected performance of the distance protection relay for a specific application case.

Another important issue to be addressed for the short-circuit sensitivity analysis is the fault location, which is composed of a configurable set of faulted elements and a variable set of locations depending on the element type to be considered (transmission line, power transformer or bus bar). The set of elements is configured according to the operating characteristic zone for which the settings are to be calculated. So, for Zone 1, only the protected line and the forward adjacent lines are to be considered, according to the expected performance of this zone, as described in Section 2.1 and Section 2.2. For the other zones, a different set of elements, including power transformers and transmission lines beyond the immediately adjacent to the protected line, must be considered.

Regarding the operating scenarios considered for all the simulated fault conditions, which represent the generation-demand balance for a set of specific conditions, only those affecting the dispatched generators in the influence zone of the protected line should be considered in order to stablish the extreme conditions for evaluating the infeed effect and polarization of the distance protection relay to be set.

According to [6] and [1], the infeed effect is dependent on the short-circuit current contribution from the line ends, and therefore, at least two extreme generation dispatch scenarios should be considered: the one increasing at a maximum the short-circuit current contribution from the line end where the distance protection relay to be set is located, and the one reducing this contribution to the minimum with respect to the opposite line end. A combination of the scenarios affecting the short-circuit current contribution at both line ends is recommended in order to consider all the possible extreme conditions affecting the infeed effect, and thus the performance of the distance protection relay.

For establishing the most representative operating scenarios affecting the polarization of distance, since according to [6] and [1] the pre-fault load transfer through the protected line has a direct impact on the performance of the cross-polarization and voltage memory polarization techniques, it is very common to find that the same operating scenarios selected for maximum and minimum short-circuit current contribution at both line ends matches those scenarios for maximum and minimum load flow through the protected line. However, there must be ensured that all of the selected scenarios represent the most severe conditions affecting the performance of the distance protection relay for the main factors described (infeed at both line ends during the fault and load flow prior to the fault), since the solution method proposed will consider both the polarization and starting blocks of the distance protection relay advanced models to be used, in order to characterize the expected performance of the operating characteristic for the protected line in a given transmission network. This is why the network model is of utmost importance for the application of the proposed solution, and it requires a high level of detail for representing the different operating conditions the power system is usually is subjected to. These models are usually developed and constantly updated by the transmission system operator or by the protected transmission line owner, depending on the regulations applicable by the related country.

The sample transmission system, presented in Figure 6, was built from an existing network model, reduced to the portion of a network which will give the required topology to test the proposed solution strategy of the formulated minimization problem. The parameters of all the elements used in the sample system are presented in Appendix B, from

Table A2. Electrical parameters for the generation transformers used in the sample power system.

Transformer	Rated Power [MVA]	Nominal Frequency [Hz]	Rated Voltage [kV]		Vector Group			Short-Circuit Voltage [%]
			HV-Side	LV-Side	HV-Side	LV-Side	Phase Shift
Generation Transformers	220	60	230	16.5	YN	D	1	9.99

,

Table A3. Electrical parameters for the network transformers used in the sample power system.

Transformer	Rated Power [MVA]			Rated Voltage [kV]			Vector Group	Short-Circuit Voltage [%]
	HV-Side	MV-Side	LV-Side	HV-Side	MV-Side	LV-Side		HV-MV	MV-LV	LV-HV
XFMR1	100	70	30	220	110	34.5	YN0yn0d1	6.93	6	10.899
XFMR2	100	70	30	220	110	34.5	YN0yn0d1	6.93	6	10.899
TR_3W	180	180	60	220	110	44.6	YN0yn0d1	9.31	9.15	13.06

,

Table A4. Electrical parameters for the transmission lines used in the sample power system.

Line	Length [km]	Rated Current [kA]	Positive-Sequence Impedance [Ohms/km]		Zero-Sequence Impedance [Ohms/km]		Mutual Zero-Sequence Impedance [Ohms/km]		Susceptance [uS/km]
			R1	X1	R0	X0	R0M	X0M	B1	B0
LT_1	84	0.755	0.048831	0.35786	0.29136	1.0086	0.24072	0.5563	4.6529	2.643
LT_2	84	0.755	0.048831	0.35786	0.29136	1.0086	0.24072	0.5563	4.6529	2.643
Power Transformer Connections	0.01	0.971	0.0536	0.4998	0.2632	1.1175	--	--	3.4209	2.1723

,

Table A5. Electrical parameters for the generators used in the sample power system.

Generator	Rated Power [MVA]	Rated Voltage [kV]	Rated Power Factor	Connection	xd″ [p.u.]	xd′ [p.u.]	x0 [p.u.]	x2 [p.u.]	xd [p.u.]	xq [p.u.]
Generators	224	18	0.85	YN	0.216	0.265	0.111	0.217	1.781	1.74

,

Table A6. Electrical parameters for the network equivalents used in the sample power system.

Network Equivalent	Angle [°]	Voltage Setpoint [p.u.]	Maximum Short-Circuit					Minimum Short-Circuit
			Short-Circuit Current [kA]	R/X	Z2/Z1	X0/X1	R0/X0	Short-Circuit Current [kA]	R/X	Z2/Z1	X0/X1	R0/X0
RTS 110 kV	−5.13	1.03	17.658	0.271	1.015	1.87	0.202	13.572	0.263	1.026	1.96	0.26
ST 220 kV	−3.64	1.07	14.5	0.106	1.021	1.895	0.202	10.68	0.13	1.012	1.94	0.193

and

Table A7. Electrical parameters for the demand equivalents used in the sample power system.

Demand Equivalent	Rated Power [MVA]	Power Factor	Leading/Lagging
Demand STS	50	0.9	lagging
Demand TS 110 kV	110	0.92	lagging

.

2.5. Solution Strategy

Figure 7 shows the flowchart of the proposed solution strategy for the minimization problem given by Equations (7)–(10). As it can be seen, it is composed by the short-circuit sensitivity analysis, previously described, from which the apparent impedances of all the fault conditions included in the analysis were obtained. Then, from the block calculating the initial settings, based on the typical criteria described in Section 2.1, the initial population for the genetic algorithm proposed as the main solution strategy for the problem is obtained.

As mentioned in Section 1, a genetic algorithm (GA) has successfully been used in [13] to solve a similar problem, in comparison with conventional calculation methods. A general description of the GA is presented in [20]. The GA is inspired by Darwin’s principle of natural selection and the survival of the fittest, and has been successfully used to solve different optimization problems, recently achieving substantial progresses [21,22,23,24]. The GA was selected as the main solution strategy for the formulated optimization problem due to its advantages over other optimization algorithms, as its ease to converge and its effectiveness in obtaining a global view of the search space, as a result of its diversification and intensification technique [25]. Finally, the parameters of the GA, as it will be presented later, allow for a good integration with the nature of the formulated problem.

Upon meeting the stopping criteria of the GA, the solution strategy includes an optional local search component, whose main goal is to improve the solutions obtained by the GA. As a result, the solution strategy consists of a hybrid algorithm, in which the local search algorithm is optional, as the GA has proven to give high-quality solutions by itself, and given that the practical solutions to be included in the search space for the local search algorithm (neighborhood of alternative solutions) are conditioned by the tolerance defined by the combined error of the instrument transformers, for most of the application cases, it is expected to get good enough solutions by only applying the GA.

Appendix A includes the detailed pseudocode of the algorithm used to implement the solution method.

2.6. Setting Parameters for the Proposed Genetic Algorithm

The initial population of the GA is established based on the degree of minimum precision associated with the instrument transformers, of 1% for a current transformer of class 5P, and 3% for a voltage transformer of class 3P, which results in a composite error of maximum 4%, which will be taken as the interval size for the variation of ranges between the reach setting limits given by the aforementioned constraints of the problem (Equations (9) and (10)). This is done in order to guarantee an acceptable degree of diversity in the initial population, to ensure that the search space delimited by these constraints is adequately explored.

Given the two main factors affecting the performance of the operation characteristic of the distance protection element, selectivity and sensitivity, which correspond to the operative state of the line and the infeed effect [1,6], which may require the reduction of the zone reaches in X or in R, or in both, and that due to the size of the selected interval, which represents 25 steps, the population size is set to 75, in such a way that 25 solutions are considered in the initial population, with X1 in the maximum possible value and R1 variable, 25 with R1 in the maximum possible value and X1 variable, and 25 with R1 and X1 variating proportionally to their respective setting intervals, according to the combined error of the instrument transformers. A graphical summary of the method for stablishing the initial population is presented in Figure 8.

This otherwise non-related criteria for selecting the population size is intended to provide enough diversity to the population, so that the GA is not locked in a local minimum, and to improve the convergence, since it allows a more efficient search of the solution space. Therefore, the condition imposed on the population size is more related to achieve the previously described goal, and upon the sensitivity analysis performed for the application problem, any other number that helps achieving such goal can be used.

As stated in [26], given the nature of the problem, where fast convergence to a local minimum is always a risk, an otherwise large population size of 75, for a continuous GA, has proven to be effective in ensuring the diversity component that improves the quality of the solutions, according to the parameter setting process to be discussed below, and at the same time mitigating the risk of fast convergence to a local minimum.

For the natural selection stage, in each iteration, the 25 individuals with the best fitness (lowest value of the objective function) are selected (selection rate of 2/3), from which 50 offspring are obtained by crossing. In this way, the size of the population remains constant at 75.

To choose the pairs of parents, three candidates are taken, at random, from the subset resulting from the selection stage, and the one with the best objective function is chosen by tournament. In this way, the search is intensified in the most promising regions of the search space. This assignment of couples includes a control that prevents selecting the same individual as both father and mother. The crossover matrix is thus built up, of 25 rows for all the selected mates, and two columns, one for the position of the father solutions in the population matrix and another for the mother solutions.

Since there are only two decision variables, the crossover rate can be 0.5 or 1.0 (selectable). A crossover rate of 0.5 implies that only one of the two variables is used to cross between the father and the mother, thus obtaining the value of such variable in the child; with the variable that will not cross remaining as a random selection of the value of the father or mother, which will then be transmitted directly to the child. This is done by randomly selecting an integer number α between 0 and 1 and applying Equation (11). To perform the crossover, a real random number β between 0 and 1 is generated, and the linear combination is performed as presented in Equation (12). Setting the crossover rate at 1.0 implies that both variables will take part in the linear combination described in Equation (12). The detailed process for natural selection and crossover can be better understood in the pseudocode presented in Table A1 of Appendix A, and in the examples presented in Figure 9:

$$Va{r}_{child}=\alpha ·Va{r}_{father}+\left(1-\alpha \right)·Va{r}_{mother}$$

(11)

$$Va{r}_{child}=\beta ·Va{r}_{father}+\left(1-\beta \right)·Va{r}_{mother}$$

(12)

This crossover method allows, for the formulated problem, to guarantee the feasibility of the population, by keeping the reaches in R and X always within the borders delimited by the established constraints. After obtaining the offspring solutions, and before the mutation process, it is necessary to calculate the objective function for them, in such a way that if a new elite solution is obtained among the offspring, this will not be undone by the mutation operator.

For the process of mutation of the population formed by the parents of natural selection and the offspring of the crossover process, the mutation rate is left as an adjustable parameter, taking a value of 20% as a reference, and calculating the number of mutations based on [26], as presented in Equation (13):

$$\#Mutations=MutationRate\times \left(PopulationSize-ElitePopulation\right)\times \#Variables,$$

(13)

Before performing the mutations, a control was included to avoid mutating the best solutions (elite solutions), which will form a subset of solutions of adjustable size, and which allow to integrate a memory component that reinforces the component of intensification in the search. Then, the mutation operator consists on selecting random positions of the population matrix obtained from the crossover operator, as many as mutations are required, calculated from Equation (13), and in each of these positions, which represent a setting reach in R or in X, a random value is selected between the limits defined by the constraints of the problem, in such a way that feasibility is maintained in the solutions that will be part of the new population. This can be better understood in the pseudocode presented in Table A1 of Appendix A, and in the examples presented in Figure 10.

Finally, the aforementioned process is repeated until the stopping criteria is met, which has been established based on a maximum number of iterations, finding high quality solutions with at least 100 iterations, according to the sensitivity analysis performed in Section 3.4.

Regarding the local search algorithm, this will be optional, depending on the quality of the solution obtained from the GA, and the search strategy will be of maximum improvement. The neighborhood will be defined by a tolerance band of ±15% (corresponding to the maximum measurement error of the impedance, as discussed previously) with respect to the solution to be improved, with a movement that will divert the initial solution towards the limits of the neighborhood with an accuracy of 1%.

3. Tests and Results

3.1. Calculation of the Earth Compensation Factor (k0)

For the sample transmission system shown in Figure 6, the common earth compensation factor is calculated according to Equation (5) for a single-line-to-ground fault of zero fault impedance at bus SE_ST 220 kV, considering both mutually coupled lines, LT_1 and LT_2, in service. The resulting phasor quantities required for the calculation of the earth compensation factor are presented in

Table 1. Phasor results for the calculation of the earth compensation factor.

${\mathit{V}}_{\mathit{\phi }-\mathit{G}} \left[\mathbf{V}\right]$	${\mathit{I}}_{\mathit{\phi }} \left[\mathbf{A}\right]$	${\mathit{I}}_{\mathit{N}} \left[\mathbf{A}\right]$	${\mathit{Z}}_{\mathit{\phi }-\mathit{G}} \left[\mathsf{\Omega }\right]$
$82416\angle 1.475°$	$1322\angle -66.24°$	$1240\angle -78.88°$	$30.34\angle 82.23°$

.

After applying Equation (5) with the quantities presented in Table 1, the earth compensation factor results in $k_0=1.189\angle -{14.866}^{°}$. This will then be used to compensate the phase to ground impedance loops for ground faults, according to Equation (4).

3.2. Calculation of the Initial Settings for Zone 1

The settings calculated using the typical setting criteria described in Section 2.1 represent the start point for any protection engineer dealing with the problem of calculating the zone reach settings for mutually coupled transmission lines. The initial proposed settings for Zone 1 of the quadrilateral distance protection element to be set for a Siemens 7SA87 relay protecting line LT_1 from substation SE, of the sample transmission system shown in Figure 6, are presented in

Table 2. Initial zone 1 reach settings.

R1 [Ω]	X1 [Ω]
58.23	25.55

. The measured fault impedances stored in sets D1 and F1 are presented in Figure 11 for single-phase-to-ground and three-phase faults; and in Figure 12 for double-phase-to-ground faults. All the faults correspond to all those simulated as part of the short-circuit sensitivity analysis.

3.3. Calculation of the Reach Settings for Zone 1 Using the Proposed Method

As it can be seen, especially for most of the faults stored in set F1, the problem represents a major challenge to be solved, specially without considering the probability of occurrence of the faults simulated in the short-circuit sensitivity study.

The results presented in Figure 13 show the contrast between the obtained solution for the resistive and reactive reaches of Zone 1 and the faults included in the short-circuit sensitivity analysis, with a preference factor of 0.5 in the objective function (equal sensitivity and selectivity), and the probability distribution previously discussed. The parameters for the GA are also according to the values previously discussed (population size fixed at 75, crossover rate equal to 1, mutation rate equal to 0.2, elite population of five individuals, and a maximum number of iterations equal to 200). The solution obtained and presented in Figure 13 represents significant reductions in both the resistive (37.91 Ω) and reactive (17.03 Ω) reaches, when compared to the initial reach settings presented in Table 2. With the initial settings, the preference function for the cumulative probability of presenting a non-selective or non-sensitive operation during the set of faults included in the short-circuit sensitivity study is 17.72, whereas for the solution obtained from the proposed method, it decreases to 5.02, which represents a serious improvement in the performance of Zone 1.

However, some system operators and transmission companies give the highest priority to selectivity, sacrificing sensitivity for high-impedance faults, arguing that the impact of a non-selective trip during high power transfer through the protected line could yield to the violation of reliability indicators dictated by the local regulation or stability issues in the area of influence of the mutually coupled transmission lines. Since high-impedance ground faults represents lower risks to the protected elements and the stability of the generators near the area of influence of the protected line, the common philosophy is to let these faults to be detected by a backup directional ground overcurrent element, much more sensitive than the distance protection for these fault.

As a result, a common practice is to select the reach settings for every zone, in order to reduce as much as possible, the risk of obtaining a non-selective operation of the distance protection element. Figure 14 shows the solution obtained by assigning the same probability to the three operative states of the mutually coupled line and with a preference factor of 0.95 (solution 95% selective and 5% sensitive). The resistive reach obtained is 23.98 Ω, and the reactive reach is 16.69 Ω, with an objective function of 1.30. The significant reduction in the objective function is related to the fact that the sensitivity only represents 5% of the preference function of the algorithm.

For the results presented in Figure 14, there are still faults causing non-selective operations in the operating characteristic of Zone 1. This is due to the probability distribution among the three fault types considered. For the most conservative approach in terms of selectivity, the same probability can be assigned to all the fault types considered, and using a preference factor of 1 in the objective function. This flexibility allows that the proposed algorithm can be applied no matter the protection philosophy used to evaluate the proper performance of the distance protection element.

3.4. Sensitivity Analysis for the Mutation Rate and Number of Iterations of the Genetic Algorithm

The results presented

Table 3. Performance tests of the GA with a preference factor k = 0.8 and a mutation rate of 0.2.

Number of Iterations	Solution		Objective Function	Population Deviation from Solution	Number of Solutions with the Same Value of the Objective Function
	R Reach	X Reach
2	21.391	16.411	6.167	6.215	1
5	22.972	16.485	5.674	3.529	1
10	22.887	16.662	5.453	4.001	1
20	24.024	16.644	5.301	5.471	6
50	23.846	16.685	5.328	5.592	1
100	23.477	16.688	5.275	6.482	2
200	24.168	16.690	5.195	8.811	1
500	24.073	16.663	5.301	5.715	3
1000	24.132	16.691	5.195	6.502	1

,

Table 4. Performance tests of the GA with a preference factor k = 0.8 and a mutation rate of 0.1.

Number of Iterations	Solution		Objective Function	Population Deviation from Solution	Number of Solutions with the Same Value of the Objective Function
	R Reach	X Reach
2	20.639	16.534	6.007	4.385	1
5	22.279	16.537	5.701	3.119	1
10	23.973	16.741	5.355	1.790	1
20	56.269	16.679	5.248	2.202	35
50	22.699	16.744	5.315	6.042	1
100	34.750	22.006	5.248	6.358	8
200	24.285	16.658	5.301	7.549	1
500	23.981	16.686	5.248	8.128	1
1000	41.950	23.094	5.195	8.187	19

,

Table 5. Performance tests of the GA with a preference factor k = 0.8 and a mutation rate of 0.05.

Number of Iterations	Solution		Objective Function	Population Deviation from Solution	Number of Solutions with the Same Value of the Objective Function
	R Reach	X Reach
2	27.315	16.711	5.874	5.156	1
5	26.496	16.173	6.007	2.714	1
10	25.782	16.369	5.688	0.363	2
20	25.846	16.005	6.021	0.313	2
50	31.028	16.972	5.195	2.136	25
100	56.130	19.481	5.248	6.524	25
200	24.081	16.690	5.195	7.069	1
500	52.913	23.834	5.248	6.561	25
1000	54.210	19.057	5.195	5.085	25

,

Table 6. Performance tests of the GA with a preference factor k = 0.95 and a mutation rate of 0.2.

Number of Iterations	Solution		Objective Function	Population Deviation from Solution	Number of Solutions with the Same Value of the Objective Function
	R Reach	X Reach
2	31.079	15.391	1.864	1.734	1
5	22.488	16.344	1.502	2.656	1
10	25.393	16.161	1.479	6.294	1
20	23.852	16.683	1.332	4.251	2
50	24.289	16.615	1.339	8.209	1
100	23.080	16.688	1.339	5.960	1
200	23.999	16.675	1.325	2.793	2
500	24.025	16.646	1.325	5.183	1
1000	24.126	16.661	1.325	9.230	5

,

Table 7. Performance tests of the GA with a preference factor k = 0.95 and a mutation rate of 0.1.

Number of Iterations	Solution		Objective Function	Population Deviation from Solution	Number of Solutions with the Same Value of the Objective Function
	R Reach	X Reach
2	22.581	16.298	1.502	2.792	1
5	19.581	16.836	1.485	1.558	1
10	25.814	16.288	1.395	3.491	1
20	26.025	16.024	1.505	2.056	1
50	24.116	16.585	1.339	6.078	1
100	24.231	16.680	1.312	4.998	2
200	23.381	16.690	1.319	8.667	2
500	38.312	21.141	1.312	2.462	25
1000	23.960	16.676	1.312	3.896	1

and

Table 8. Performance tests of the GA with a preference factor k = 0.95 and a mutation rate of 0.1.

Number of Iterations	Solution		Objective Function	Population Deviation from Solution	Number of Solutions with the Same Value of the Objective Function
	R Reach	X Reach
2	31.079	15.391	1.865	1.503	1
5	33.962	14.652	1.841	3.904	1
10	27.184	15.827	1.532	0.309	1
20	23.206	16.660	1.345	2.544	1
50	33.298	15.742	1.319	0.964	25
100	50.784	23.095	1.299	3.040	25
200	49.696	16.688	1.299	6.092	9
500	24.288	16.692	1.299	6.890	4
1000	45.948	23.014	1.299	6.800	25

correspond to a sensitivity analysis aimed to determine the impact of the GA parameters related to the number of iterations and the mutation rate. The same cases covered before, in Section 3.3, are now presented, with each table showing the results obtained with different values for the number of iteration parameters, for a given preference factor and mutation rate. The population deviation from the final solution (PDFS) has been obtained from Equation (14), where $M_i$ and $M^{\ast }$ represent the value of the preference function for every single solution available at the final population matrix and for the selected solution with the lowest value, respectively; and $N$ corresponds to the population size, which has been fixed at 75:

$$PDFS=\sqrt{\frac{{{\displaystyle \sum }}_{i=1}^N{\left(M_i-M^{\ast }\right)}^2}{N}}$$

(14)

The value of PDFS shows how diversified the final population is, and if the population size and mutation rate parameters are accomplishing the goal of keeping the GA from stagnating in local minimum solutions within the search space, as defined by the problem constraints.

From Table 3, Table 4 and Table 5, the results for the preference factor of 0.8 and variable mutation rate, according to the values recommended in [26], are presented. From Table 6, Table 7 and Table 8, a preference factor of 0.95 and variable mutation rate have been used. The two values chosen for the preference factor are intended to guide the solutions to the most commonly used practices among protection engineers, which gives a higher priority to selectivity rather than sensitivity, when stablishing the reliability balance of the distance protection element.

The results show, in general, that the number of iterations required to obtain high quality solutions is at least 100, and that the mutation rate affects the diversity of the population, for which the highest value of 0.2 is recommended, in order to avoid any risk of compromising the quality of the solution obtained from the GA, caused by a fast convergence into a local minimum.

Another important subject to be discussed, based on the results, is related to the different solutions obtained with the same value in the objective function, but very different values for the reach settings of the distance protection operating characteristic. Since the formulated problem corresponds to a multi-objective one, it is expected that more than a single solution satisfies the balance between selectivity and sensitivity, which is the main goal of the proposed algorithm.

The different solutions with the same value of the objective function represent different alternatives to achieve such balance, which has, as its main parameter, the preference factor k included in Equation (6). This factor, as discussed previously, is higher than 0.5, for which selectivity is prioritized over sensitivity in the objective function, according to the most commonly used practices being used by protection engineers, according to the guidelines included in [1] and similar references. Since the objective function prioritizes selectivity, the best way to select the proper reach settings, among those solutions with the same value in the objective function, is according to the maximum reach possible, which can be estimated by multiplying the values of the resistive and reactive reaches, and selecting the solution with the maximum reach among them, in order to increase sensitivity as much as possible without affecting the selectivity.

This is a very simple, yet practical and valid alternative for choosing the best possible solution, without the use of more complex algorithms, like NSGA-2, since it is entirely based on the protection philosophy and expected performance, as presented in Section 2.2.

It is important to perform the short-circuit sensitivity analysis as detailed as possible, including different operating scenarios, variating the dispatch of the generators in or near the influence zone of the protected line, and adding critical contingencies that could change the performance of the infeed effect in the distance protection relays. This is of utmost importance, since the short-circuit sensitivity analysis is the element used to characterize the proper performance of the distance protection element.

4. Conclusions

In this paper, an alternative method to solve the optimization problem related to the determination of the reach settings of the quadrilateral distance protection element of mutually coupled transmission lines is proposed, improving the technical basis for previous formulations of the optimization problem, which included a short-circuit sensitivity analysis for establishing the actual performance of the relay to be set under the multiple fault and system conditions affecting the measured impedance loops, and including the probability of occurrence of every single fault simulated in an specialized software tool, in order to determine the preference function related to the optimization problem to be solved. In addition, the proposed solution method is flexible enough to accommodate to special performance requirements of the distance protection elements, which will depend on the evaluation of the fault conditions obtained from the short-circuit sensitivity analysis. As the main component of the solution strategy, a GA was considered, taking advantage of the relation of the formulated problem with the multiple operators used by this algorithm to ensure obtaining high-quality solutions.

The solutions obtained with the solution strategy implemented were compared in performance with the general procedure for setting and validating the quadrilateral distance protection elements, resulting in a significate improvement of the performance of such elements.

It was found that the local search component of maximum improvement was effective only when the number of simulated faults was very low, which makes sense, since the greater space between impedances necessarily represents a smaller impact on the objective function evaluated in the individuals of the population of the GA. For this reason, and given that the developed algorithm allows defining the number of faults simulated on the same analyzed network, it was finally decided to make the application of the local search component optional, and consequently, it is recommended that for the effective application of the algorithm, it is necessary to reduce as much as possible the distance between fault impedances and fault locations, in such a way that the proposed objective function captures more faithfully the behavior of the solutions obtained in each iteration, significantly improving the performance of the GA.

The proposed algorithm can be used to determine the reach settings of any protection zone, under any operating scenario, and the computation times are considerably low compared to the time it takes to perform this type of analysis in practice, and can be reduced even more if the user has some experience in the analysis of this type of problems, by allowing the adjustment of the amount and characteristics of the simulations necessary to obtain a high-quality solution.

Since the solutions obtained will always leave faults undetected due to the balance reached between selectivity and sensitivity, which depends on the preference factor used in the objective function of the formulated optimization problem, for a complete validation of the protection scheme of the transmission lines, it is important to continue with the research work in order to include the operation of the backup protection elements for the set of undetected faults, considering a basic evaluation of sensitivity and selectivity with the adjacent protection elements. This would allow detecting critical faults which could yield to undesired operations due to loss of selectivity or loss of sensitivity of the entire protection scheme of a given transmission line in any system.