Deep belief improved bidirectional LSTM for multivariate time series forecasting

1.   Introduction

Incorporating distributed generators (DGs) or distributed energy resources (DERs) into power systems has become a significant driving force in today's dynamic and evolving energy landscape. These generators, including renewable and non-renewable sources such as solar panels, wind turbines, diesel engines, and microturbines, are increasingly common in modern power networks. While they promise cleaner, more sustainable energy generation, their integration introduces unique challenges, particularly in relay coordination. Power relays play a crucial role in maintaining power system stability by detecting and isolating faults and ensuring uninterrupted electricity supply to consumers. However, traditional relay coordination practices face new complexities as DG penetration increases. We explore the intricate protection coordination issue of overcurrent relays (OCRs) and examines the substantial impact of DGs on this critical aspect of grid operation.

The discipline of power system relay coordination is a well-established practice aimed at crafting protective schemes that guarantee the selective operation of protective devices during fault occurrences. Historically, this coordination has centered on centralized generation sources, enabling meticulous calibration of relay settings according to a predefined hierarchy. Historically, relay coordination has focused on centralized generation sources, allowing for precise calibration of relay settings according to a predefined hierarchy. In this traditional framework, generators are typically large, centralized units with predictable operations. Relay coordination involves adjusting parameters such as time-current characteristics, fault-clearing times, and coordination margins to ensure that the nearest relay responds quickly to a fault, isolating it while maintaining the stability of the overall system. The integration of DGs presents several challenges to the traditional relay coordination paradigm:

1. Bidirectional power flow ^[1,2]:

DGs are typically integrated into feeders' middle or end sections in medium-voltage (MV) and low-voltage (LV) power systems. This integration causes the power flow provided by DGs to be opposite to the original flow supplied by the utility company. Such a situation will also occur during system faults.

2. Fault current variability ^[2,3]:

The presence of distributed generation can alter fault current levels, complicating the coordination of OCRs. This variability can result in situations where the primary protection may not operate as intended, leading to potential equipment damage or safety hazards. Depending on the location of the fault, the fault current may increase at some points compared to the original design, while at other places, the fault current may be offset and reduced. This change can cause malfunctions in existing protection schemes and instability in the network.

3. Intermittency and variability of power generation ^[1,2]:

DGs exhibit inherent variability and intermittency, particularly renewable sources like solar and wind. DGs can cause voltage variations, including overvoltage, particularly when generation exceeds local load demand. Their power output depends on weather conditions and time of day, making it difficult to accurately predict and coordinate relay settings.

4. Protection coordination ^[1,2,3]:

The traditional overcurrent relay coordination schemes may need to be more effective in the presence of DG. The existing protection schemes may not adequately respond to the changes introduced by fault current variability. It necessitates new parameterization curves for overcurrent relays (ANSI 50/51 function) to ensure proper protection coordination. Therefore, we need new coordination strategies to maintain reliability and security in the protection system.

To address these challenges, power system engineers are developing innovative solutions that account for the distributed nature of generators. Adaptive relay coordination schemes, which leverage real-time data and advanced algorithms, are becoming increasingly important. These schemes enable dynamic adjustments to relay settings based on the grid's actual operating conditions, including the presence and behavior of DGs. Adaptive relay coordination enhances reliability and resilience in power systems by ensuring that protective devices respond appropriately to faults while minimizing unnecessary DG tripping. Additionally, these strategies improve the overall efficiency and stability of the grid, making it more adaptable to the variability of renewable energy sources.

As mentioned, DG integration can significantly affect the coordination of overcurrent relays in distribution systems. Reference ^[4] introduced a 'protection coordination index' (PCI) as an effective metric for planning the protection of meshed distribution systems with DGs. They employed a two-phase nonlinear programming (NLP) optimization method to calculate the PCI, considering variations in maximum DG penetration levels as the protection coordination time interval (CTI) changes. In ^[5] and ^[6], the multi-agent theory was applied to address protection coordination challenges. In ^[5], the impact of DG grid connections was considered using the Java Agent Development Framework platform to adjust the protection coordination curve. Moreover, the researchers in ^[6] addressed feeder reorganization issues due to DG connections, dynamically adjusting relay settings using signals from phasor measuring units (PMUs) via a central control system at the power station.

The protection coordination problem has traditionally been modeled as an optimization problem, utilizing mixed-integer nonlinear programming (MINLP) and mixed-integer linear programming (MILP) to determine relay settings in response to fault current variations from DG integration. As shown in ^[7], the authors present an integrated approach that considers both distribution system reconfiguration and protection coordination, which are formulated as a mixed integer nonlinear programming (MINLP) problem and a nonlinear programming (NLP) problem, respectively. The researchers in ^[8] addressed the concept of clustering the various network topologies into a limited set is proposed where an optimal setting group is determined for each network topology cluster. Due to the complexity and non-linearity, the setting group protection coordination problem was reformulated as a MILP problem. In ^[9], the interior-point (IP) algorithm was applied to develop optimal adaptive protection coordination strategies for scenarios involving high penetration of green generation and network reconfiguration. Reference ^[10] proposed an approach for protection coordination in microgrids incorporating non-standard characteristic features of directional over-current relays (DOCRs), modeled as a MINLP problem. In reference ^[11], the authors proposed an adaptive protection scheme that updates relay trip characteristics based on the system's operating condition, whether grid-connected or islanded. Local information, such as DG status and fault currents, determines the system's mode, and the relays adjust their settings accordingly. The overcurrent relay (OCR) trip characteristics are updated based on the system's operational state and faulted section detection. This adaptive system relies solely on local data, avoiding the complexity and cost of communication between relays. The system ensures faster fault clearance, particularly in islanded mode when fault currents are lower. Simulations demonstrate that the adaptive system clears faults faster in islanded mode, preventing unnecessary generator disconnections and minimizing load losses.

Metaheuristic intelligence-based random search and evolutionary process techniques have been widely applied to adaptive coordination decisions, as documented in references ^[12,13], where firefly optimization (FA) was used to develop optimal adaptive protection coordination strategies. The researchers in ^[14] employed a genetic algorithm (GA) to solve the optimal DG placement problem, maximizing DG penetration in medium-voltage (MV) distribution networks without altering original relay protection parameters. GA was also used to determine the optimal sizes and locations of DGs. Similarly, reference ^[15] addressed optimal DG placement, using optimization methodologies to classify cases as either 'coordination holds' or 'coordination lost', considering DG location and capacity variation changes. The researhers in ^[16] introduced a new methodology for optimizing overcurrent relay (OCR) settings using the Gorilla Troops Optimizer (GTO). This approach aimed to enhance the coordination of protective devices by factoring in transformer phase shifts and DG integration to minimize outages during faults. The methodology was tested on a real distribution network in a desert environment. GTO effectively solved the relay coordination problem while accounting for practical constraints like transformer phase shifts and fault conditions. Various metaheuristic techniques, including GTO ^[17,18], the water cycle algorithm (WCA) ^[17], particle swarm optimization (PSO), GA, teaching-learning-based optimization (TLBO), and the shuffled frog leaping algorithm (SFLA) ^[18], have been implemented to solve the optimal coordination problem.

In this paper, we address the protection coordination challenge using an optimization model. The authors introduce the concept of tracking routes, categorized into main feeder and branch line routes, represented as relay numbers based on the feeder connectivity topology. The objective is to minimize the cumulated operating time (COT) of all relays on each tracking route, using this as an index to evaluate the suitability of relay coordination. Since different combinations of time multiplier setting (TMS) and pickup current setting (PCS) result in varying relay operation time (OT), the focus is on determining the optimal parameter combinations while adhering to the TMS and PCS limits of numerical OCRs. The relay protection coordination problem is then formulated to incorporate CTI constraints for each primary and backup relay pair along every tracking route. The optimal protection coordination model will be simulated under various DG integration scenarios at different penetration levels.

The immune algorithm (IA) draws inspiration from the operation of the human immune system, our fundamental biological defense against viruses and pathogens ^[19]. This intricate defense mechanism involves gene combinations to combat invading antigens. Leveraging this heuristic algorithm, IA exhibits several advantages and characteristics that often result in superior performance compared to many other algorithms ^[20]. In their earlier work ^[21,22], the authors delved into an adaptive evolution mechanism featuring a self-adjusting scheme, successfully combined with the other metaheuristic intelligence-based algorithm to address premature convergence issues. This article presents a refined immune algorithm (RIA) with an auto-tuning reproductive mechanism (ATRM). The authors also implemented the tabu search mechanism ^[23] to avoid searching the local optima solutions, violating constraints, and solutions visited. This algorithm determines TMS and PCS for each OCR along the tracking routes while considering CTI constraints and limitations of the OCR parameters. Then, a 16-bus actual distribution network and the IEEE 37 Bus test system with DGs are employed for simulation to evaluate the effectiveness of the proposed adaptive protection coordination approach. The results demonstrate that the proposed method can rapidly and consistently identify solutions with maximum fitness, minimizing COT. Moreover, it effectively mitigates the impact on protection coordination settings caused by incorporating DGs into the system. All the TMS and PCS solution combinations were verified using the renowned power system simulation software ETAP, where the time-current curves (TCC) of each tracking route were plotted, and the protection coordination sequence was validated. The proposed RIA-ATRM was used to evaluate the optimal parameters of OCRs of the mentioned distribution system with distributed energy resources. The results demonstrate that the RIA-ATRM can successfully deal with the optimal coordination issue and effectively reduce the total COT on the overall tracking routes, further mitigating the impact of DERs on protection coordination.

2.   Relay setting and optimal formulation of protection coordination

The parameter setting of the OCRs is an important measure to ensure that the circuit breaker can cut off the power supply at the pre-scheduling time in the event of an overload or short circuit fault. When a short circuit fault occurs, the fault current will be generated more significantly than the system's regular operation. To maintain the stability of system operation, the OCRs must trip accurately and timely to effectively isolate the faulty area, protect equipment, and reduce the impact on the typical operating region. Two parameters, TMS and PCS, can be used to set up or tune the time-current curve (TCC), which describes the relationship between the OCR's operating time, pick-up current, and fault current. According to the definition in the IEC 60255-3 standard ^[24], the curve equation is shown in Eq (1). The start-up current of relays, $IP$, can be calculated by multiplying the CT ratio by the corresponding PCS as shown in Eq (2).

where

$t$: Relay OT.

$TMS$: Time multiplier settings.

$PCS$: Pick-up current settings.

$IP$: Pick-up or start-up current of the relay.

$If$: Fault current through the relay.

$CT\_R$: CT ratio setting of the corresponding relay.

$\alpha , \beta$: Characteristic parameters of the standard time-current function.

The curve can generally be classified into short-time, long-time, inverse time, very inverse time, and extremely inverse time, with the corresponding characteristic parameters provided in references ^[24,25]. The coordination of adjacent relays primarily depends on the operating time difference between the primary and backup relay, known as the mentioned CTI. All the relays are appropriately coordinated when the CTI for all primary and backup relays in the system is within a specified limit or allowable range. The CTI value is usually set between 0.2~0.5 sec. ^[24], the calculation method of CTI is

where

$CTI$: Coordination time interval (sec.).

${t}_{b}$: OT of backup relay (sec.).

${t}_{p}$: OT of primary relay (sec.).

The relay protection coordination problem can be framed as a nonlinear programming problem aimed at minimizing the total operating time of all relays. The operating time of a relay is contingent upon its TMS and PCS settings for a specific fault current. Determining appropriate PCS and TMS values for all relays of the distribution system and effectively coordinating their operations presents a significant challenge ^[26,27]. In this paper, the settings of TMS and PCS can be obtained by formulating the OCRs coordination issues as an optimization problem in the objective function expressed as

where

Obj_F: Object function of calculating the COT for all specified tracking routes.

${W}_{i, j}^{k}$: The weighting factor of relay i when the j^th fault occurs in the k^th coordination tracking route condition.

$i$: The relay identifier.

$j$: The fault location identifier.

$k$: The coordination tracking route identifier (relay connectivity).

FN: The total number of fault locations is considered.

TP: The total number of coordination tracking routes.

RN: The total relay number of the set of tracking routes.

${TMS}_{i}$: TMS parameter of relay i.

${PCS}_{i}:$ PCS parameter of relay i.

${If}_{i, j}^{k}$: Fault current through the relay i when the j^th fault occurs in the k^th coordination tracking route.

${IP}_{i}$: Pickup or Start-up current of relay i.

${\alpha }_{i}, {\beta }_{i}$: Parameters of the standard normal time-current characteristic function ^[24] for relay i.

${CT\_R}_{i}$: CT ratio corresponds to relay i.

Eq (4) serves as the objective function, aiming to minimize the COT for both primary and backup relays on the tracking routes across fault positions. Based on different fault locations, a set of protection coordination tracking routes can be determined, which may partially overlap. As a result, the relay operation time on overlapping paths is accumulated repeatedly. To evaluate the relay's operation time, consider only the normal inverse time-current characteristic function outlined in the IEC standard ^[25], represented by Eq (1). The fault current in Eq (1) can be simulated using various power engineering software such as ETAP, PSS/E, etc. Here, the parameters α_i and β_i of the standard normal time-current curve are set to 0.14 and 0.02, respectively.

According to Eq (3) it indicates that the protection relays are effectively coordinated. The evaluation of CTI and its constraints are explicitly outlined in Eqs (7) and (8). Subsequently, the selection of the TMS and PCS value hinges on the brand and model of the chosen relay. Additionally, there are limitations on its upper and lower bounds, as depicted in Eqs (9) and (10).

where

${CTI}_{x, j}^{k}$: CTI of the x^th pair of primary/backup relay when the j^th fault occurs in the k^th coordination tracking route.

${t}_{{i}_{p}, j}^{k}$: OT of primary relay i_p when the j^th fault occurs in the k^th coordination tracking route.

${t}_{{i}_{b}, j}^{k}$: OT of backup relay i_b when the j^th fault occurs in the k^th coordination tracking route.

${CTI}_{}^{min}$: The lower CTI bound is set to 0.2 seconds.

${CTI}_{}^{max}$: The upper CTI bound is set to 0.35 seconds.

${TMS}_{}^{min}$: The lower bound of TMS is set to 0.05.

${TMS}_{}^{max}$: The upper bound of TMS is set to 1.0.

${PCS}_{}^{min}$: The lower bound of PCS is set to 0.05.

${PCS}_{}^{max}$: The upper bound of PCS is set to 5.0.

3.   Optimal algorithm—a refined immune algorithm with an auto-tuning reproductive mechanism

A refined IA with an auto-tuning reproductive mechanism (RIA-ATRM) was developed to enhance the performance of IA, as discussed in the following sections. Furthermore, the crossover and the mutation mechanism were refined by a competition and auto-selection scheme according to the quality of the solutions of the two generations before and after to avoid prematurity, and a competition mechanism was implemented to determine the choice of either one or both automatically. With the advantages of both heuristic ideals and AI, RIA-ATRM supersedes the original ideals threefold: The complicated problem can be solved, it can achieve better performance than IA, and it is more likely to reach a global optimum solution than other heuristic methods.

3.1.   Immune system

Figure 1 illustrates the immune system's primary function: identifying and eradicating foreign viruses or antigens. This adaptive biological system comprises four crucial subsystems: The antibody recomposition system, restraint system, memory system, and lymphatic system. The lymphatic system plays a pivotal role in identifying invading antigens or viruses by cross-referencing them with the pattern memory system. If a matching pattern is found, lymphocytes facilitate the production of corresponding antibodies through clonal proliferation within the antibody recomposition system, aiming to eliminate the invader. Conversely, if the antigen presents a novel pattern, it is stored in the memory system, prompting the recomposition system to generate new lymphocytes tailored to combat this new threat. Genetic evolution operations such as crossover, mutation, and reproduction occur within the recomposition system.

The lymphatic system comprises two major cell types: T-lymphocytes, produced by the thymus, and B-lymphocytes, produced by bone marrow. Newly formed lymphocytes undergo scrutiny by the restraint system to ensure compatibility with biological constraints. Once approved, they are memorized and undergo clonal proliferation to combat the external antigen. The immune system exhibits four key characteristics: self/non-self-recognition, specificity, diversity, and memory. Self/non-self-recognition distinguishes between internal and external antigens. Specificity allows for the classification of different antigens. Diversity enables unique responses to various antigens. Memory ensures immunization upon subsequent encounters with the same virus.

IA operates as a stochastic search algorithm inspired by natural selection and genetics. Antigens represent the problem's objective functions in this framework, while antibodies symbolize feasible solutions. The genetic structure closely resembles that of the GA, encompassing crossover, mutation, and reproduction stages for antibody recomposition. The following subsections include a breakdown of the algorithm.

3.2.   Encoding of genes and chromosomes

The coding scheme of chromosomes of the antibody is illustrated in Figure 2, where it was divided into two categories of chromosomes—one is for the determination of the setting parameter of TMS, and the other is for PCS. These two parameters of each relay will be multiplied by 100 and then converted into an 8-bit binary code, each as the gene, and then a chromosome was assembled by several genes' binary strings. Each chromosome indicates a combination of each relay's TMS or PCS settings. If the IA or RIA-ATRM search is terminated, each gene will then be decoded and divided by 100 separately. For instance, in Figure 2, BIN_TMS_i(PN) and BIN_PCS_i(PN) are the binary strings of TMS and PCS of the ith relay, and PN indicates the population index—from 1 to the maximum population number. The entire encoding calculation process can be completed by Eqs (11)–(14).

where

D2B: The decimal to binary conversion subroutine.

${TMS}_{i}$: TMS random value for the i^th relay.

${PCS}_{i}$: PCS random value for the i^th relay.

${TMS}_{i}^{max}, {TMS}_{i}^{min}$: The upper and lower limit of the TMS setting of the i^th relay.

${PCS}_{i}^{max}, {PCS}_{i}^{min}$: The upper and lower limit of the PCS setting of the i^th relay.

$bit$: The number of bits in a chromosome.

${resol}_{x}$: The resolution of the split between the upper and lower limit value.

3.3.   Affinity and diversity evaluation

The immune system produces various antibodies based on the affinity recognition between antigens and antibodies, or between two antibodies. In immune affinity, there are two categories. The first pertains to the affinity between antigens and antibodies, indicating the strength of their interaction ^[19,20]. The second category concerns the affinity between two antibodies, reflecting their similarity. From Information theory, entropy can be applied to measure the diversity of antibodies. It can be computed as

where N is the number of antibodies, and P_ik is the probability of the i^th allele coming out of the k^th allele. For example, if all alleles at the k^th antibody are the same, ${E}^{k}\left(N\right) = 0$. Thus, the total diversity of the k^th antibody is

where M is the number of genes of the k^th antibody.

Furthermore, two affinity forms must be taken into account in the evaluation procedure. One is the affinity between two antibodies, j^th and k^th, which can be calculated below

where E(2) is the total diversity of the j^th and k^th antibody only. Note that if two antibodies are of the same, $({Aff}_{b}{)}^{jk}$ is equal to 1. Moreover, $({Aff}_{b}{)}^{jk}$ is set between zero and one. The other one shown in the expression below is applied to investigate the affinity between antibodies and antigens.

where ${{Obj}_{F}}_{k}$ is the objective function value for the k^th antibody. The affinity score of each antibody is obtained by calculating the objective function. If one or more variables violate their limits, the corresponding chromosome will be put into the tabu list ^[23] to avoid generating the same infeasible solution again.

3.4.   Fitness evaluation

The affinity score of every chromosome is derived by evaluating the objective function described in Eq (4) while taking into account both equivalent and inequivalent constraints outlined in Eqs (7)–(10). If one or more variables of the antibody violate their limits, the antibody will be punished by multiplying a punishing factor, leading to a lower fitness value. The affinity function in Eq (18) of the coordination problem must be modified as

where

${PUNF}^{TMS}$: The total summation of the TMS punishing factor of all relays of each antibody.

${PUNF}^{PCS}$: The total summation of the PCS punishing factor of all relays of each antibody.

${PUNF}^{CTI}$: The total summation of CTI punishing factor of all relay pairs of each antibody.

$\rho , \epsilon , \mu , \phi$: Proportional constants.

${PUNF}_{i}^{TMS}$: TMS punishing factor of relay i.

${PUNF}_{i}^{PCS}$: PCS punishing factor of relay i.

${PUNF}_{j}^{CTI}$: CTI punishing factor of relay pairs j.

${P}_{1}{, P}_{2}{, P}_{3}$: The punishing number of relay i is set to more than 10³ when the TMS/PCS/CTI value is against the boundary.

3.5.   Production of the offspring

In IA, the offspring are the new chromosomes generated from crossover and mutation processes ^[19,20]. These two operations are called simple crossover and mutation schemes (SCM). To avoid the prematurity of the convergence of the IA, an auto-selective crossover and mutation mechanism is proposed instead of SCM, which is mentioned as ATRM and is stated in detail in this section. These two schemes are described as follows:

1. SCM scheme

The crossover operation selects two parents randomly for exchanging chromosomes, governed by a crossover rate (P_C) in a uniform probability distribution. The position of a gene within a chromosome is termed a loci, and the crossover point is chosen randomly from these loci. If either or both offspring are deemed infeasible, another mate is selected for crossover. Conversely, using a uniform probability distribution, the mutation process randomly picks one parent with a mutation rate (P_M), with loci randomly selected for mutation. If the resulting offspring is infeasible, another parent is chosen until a feasible solution is attained. In the computation of IA, both P_C and P_M are fixed values ranging between 0.0 and 1.0, where the sum of P_C and P_M equals 1.0. A detailed description of IA can be found in Refs. ^[19,20] and will not be reiterated here.

2. ATRM scheme

In the IA search process, crossover typically precedes mutation. Within IA, a higher P_C facilitates exploration of the solution space around the parent solution. Moreover, the P_M governs the likelihood of introducing new genes, thereby exploring fresh solution territories. If P_M is too low, the solution might converge prematurely to a local optimum. Conversely, a high P_M could introduce excessive uncertainty. When offspring deviate significantly from their parents, the algorithm may fail to leverage past learning, potentially leading to instability. Hence, selecting appropriate crossover and mutation rates poses a dilemma in IA optimization. The ATRM scheme was proposed to avoid such a difficulty, and it is illustrated as follows:

Ⅰ. The operation process encompasses three categories: crossover, mutation, and hybrid (crossover and mutation). Then, the operation mechanism of each antibody generated in the previous iteration has been documented.

Ⅱ. Select two parent antibodies randomly to produce offspring according to:

(a) If rand_C < P_C^(g) and rand_M < P_M^(g): neither crossover nor mutation process be executed;

(b) If rand_C ≥ P_C^(g) and rand_M < P_M^(g): execute the crossover process only;

(c) If rand_C < P_C^(g) and rand_M ≥ P_M^(g): execute the mutation process only;

(d) If rand_C ≥ P_C^(g) and rand_M ≥ P_M^(g): hybrid process be executed.

where

rand_C: The uniform random number in (0, 1) for crossover.

rand_M: The uniform random number in (0, 1) for mutation.

g: The current generation numbers.

P_C^(g): The control parameter crossover process with initial value P_C⁽⁰⁾ = 0.5 and 0 ≤ P_C ≤ 1.

P_M^(g): The control parameter mutation process with initial value P_M⁽⁰⁾ = 0.5 and 0 ≤ P_M ≤ 1.

Offspring generation continues until all parents are processed. Figure 3 illustrates the initial relationship among three operational procedures, each with an equal probability of generating offspring. The related procedures are not executed if the randomly generated crossover (rand_C) or mutation(rand_M) numbers are below their control parameters. Mutation is more significant than IA and GA, as it is better at exploring new regions. When the search nears a local or global optimum, mutation may need to take precedence, especially if crucial beneficial genes are absent in a generation. Since all procedures are random operators, determining the superiority of one over the others is impossible.

Ⅲ. A competition mechanism is implemented in the search process according to the fitness score. For instance, if the best antibody in the present generation comes from the hybrid process, there is more likelihood that this procedure will generate better offspring for the next generation. The area of the hybrid procedure must be expanded by reducing P_C^(g) and P_M^(g) to increase the probability, as shown in Figure 4.

Speak more precisely, if the best affinity of generation g-1 is less than that of generation g, i.e., ${Aff}_{max}^{\left(g-1\right)}\le {Aff}_{max}^{\left(g\right)}$, comes from the hybrid operation, the control parameters will decrease using the following manner

where K₁ and K₂ are the regulating factors, and in general, K₁ < K₂. g_max is the maximum iteration number. Figure 4 shows the variation in the probability of the crossover and mutation areas. On the contrary, if the best affinity of generation g-1 is greater than that of generation g, i.e., ${Aff}_{max}^{\left(g-1\right)} > {Aff}_{max}^{\left(g\right)}$, comes from the hybrid procedure, both control parameters should increase

The other two procedures are more likely to produce superior offspring, requiring both control parameters to increase to reduce this likelihood, as depicted in Figure 5. If the best solution remains unchanged, crossover and mutation operations must also restrain themselves to regain the associated territory.

It is worth noting that there is no restriction that P_M plus P_C must be equal to 1.0 in the ATRM mechanism, and it is not difficult to find that these two variables will operate independently from Eqs (23)–(26).

Ⅳ. If the best affinity of generation g-1 is less than that of present generation g, i.e., ${Aff}_{max}^{\left(g-1\right)}\le {Aff}_{max}^{\left(g\right)}$, comes from only the crossover procedure, the control parameter will decrease using Eq (23). Conversely, if ${Aff}_{max}^{\left(g-1\right)} > {Aff}_{max}^{\left(g\right)}$, comes from the crossover, the control parameters will increase by employing Eq (25). In this situation, the control parameter P_M is fixed. The probability variation of the crossover is illustrated in Figures 6 and 7.

Ⅴ. If the best affinity of generation g-1 is less than that of generation g, i.e., ${Aff}_{max}^{\left(g-1\right)}\le {Aff}_{max}^{\left(g\right)}$, which comes only from the mutation procedure, the control parameter will decrease by employing Eq (24). If ${Aff}_{max}^{\left(g-1\right)} > {Aff}_{max}^{\left(g\right)}$ the control parameters will increase by employing Eq (26). In this situation, the control parameter P_C is fixed. The probability variation is illustrated in Figures 8 and 9, separately.

3.6.   Antigen recognition

In this paper, we utilize a tabu list ^[23] as the central logic for the judgment of the antigen recognition system, including:

Ⅰ. Visited solutions, excluding the best solution in the current generation.

Ⅱ. Any previously encountered local optima.

Ⅲ. Antibodies that violate constraints.

Ⅳ. Solution spaces that do not meet the bargaining condition.

3.7.   RIA-ATRM procedure for optimal protection coordination

The RIA-ATRM process of generating new antibodies with the best affinity is continued until the affinity values are optimized or the maximum generation number is reached. The flow chart of the RIA-ATRM process is shown in Figure 10 and described as below:

Step 1. Data Input and Collection;

Read power system data, including feeders, DGs, relays, loads, and fault analysis data.I Initialize the RIA-ATRM program setting.

Step 2. Generate TMS and PCS and antibodies Coding;

Randomly generate the TMS and PCS combination for each relay, and then the code converts these values to the binary string for assembling antibodies.

Step 3. Fitness evaluation;

Calculate operation time, CTI, and objective function according to Eqs (1), (3), and (4), respectively. The code also evaluated each antibody string's fitness, affinity, and diversity to judge adaptation.

Step 4. ATRM process;

To produce new TMS and PCS combinations(antibodies) under the ATRM procedure. If the recomposition system makes a new antibody, the new antibody must be recorded in the memory cell(Tabu list).

Step 5. Did the convergence occur? Or is the maximum generation number reached?

Yes: Go to Step 6; No: Go to Step 3;

Step 6. Program Termination. List the optimal TMS and PCS combinations.

4.   Simulation result and discussion

The proposed RIA-ATRM is implemented to address the optimal coordination problem of a 16-bus actual distribution network and IEEE 37-bus radial distribution systems with varying DG penetration rates (PR). The structures of the test systems are shown in Figures 11 and 12, while the DG data are presented in Tables 1 and 2. The network topology of the 16-bus system is shown in Figure 11 as Case 1, and the system includes a main transformer with a rated capacity of 50 MVA, 161 kV/22.8 kV, and six distributed generators. The location, capacity, and generation power of each DG are shown in Table 1. The rightmost column provides two generation values, representing the power output for 60% PR and 80% PR cases, respectively. Likewise, Figure 12 demonstrates the IEEE 37-bus distribution system with five DGs integration as Case 2, and Table 2 also describes the operation data of the DGs for different load penetrations of 60% and 80%. The optimal protection coordination problem solved by RIA-ATRM was coded by MATLAB 2024 software. All programs were executed on a personal computer with Intel Core i7-10510U 2.3 GHz CPU and 16.0 GB RAM.

4.1.   Impact of DERs on distribution systems protection coordination

In this section, we first conducted a protection coordination design for the test systems without considering DG integration. After DG integration, protection coordination violations in the system are observed under different DG penetration rates. We utilize ETAP to solve the fault current passing through every relay when the fault occurred at the bus with high penetration conditions. In Case 1, when DG's PR is 60% and a fault occurs at Bus 3, the fault current is concentrated toward the fault location because the DG was placed at the feeder terminal. This fault current increases the fault current through Relay 4 from 2.459 kA to 4.788 kA, as shown in Figure 13, causing the OT of Relay 4 to be 1.42 seconds earlier than expected. The DGs' fault current contribution leads to an ill-coordinated condition, altering the relay tripping sequence to Relay 3 → Relay 2 → Relay 4 → Relay 1, as shown in the time-current curves (TCC) in Figure 14.

Similarly, in the same case, when PR rises to 80%, a fault also occurs at Bus 3. The variation of fault current solved by ETAP is shown in Figure 15. The fault currents for Relay 4 and Relay 5 are more than twice as large compared to the scenario without DG integration, even when PR is 60%. Therefore, it results that the relay tripping sequence changes to Relay 3 → Relay 4 → Relay 2 → Relay 5 → Relay 1, as shown in the TCC of Figure 16. The integration of DG impacts the protection coordination design by introducing additional fault current, which disrupts the original protection coordination sequence. It is essential to retune the relevant protection relay settings to ensure the protection coordination mechanism of the distribution system remains stable and reliable under these new operating conditions.

In Case 2, the fault current passing through Relay 12 increases from 167 A to 928 A and 1117 A, as shown in Figure 17, when a three-phase short-circuit fault occurs at bus 702 in the IEEE 37-bus system, with PR at 60% and 80%, respectively. Additionally, the fault current flowing through other relays—Relay 1, Relay 2, and Relay 3 shows a slight decrease as PR increases.

These fault current variations cause the OT of Relay 12 to be 1.91 seconds, which is the maximum OT value allowed by the preset parameters. The DGs' fault current contribution also results in ill-coordinated conditions, altering the relay tripping sequence to Relay 3 → Relay 2 → Relay 12 → Relay 1, as shown in the TCC in Figures 18 and 19.

The analysis reveals that a fault occurs in the power system; DGs and utility generators contribute fault current, which converges toward the fault location. Depending on the fault position, the fault current flowing through each relay differs from the original pick-up current values set during the initial design (without DG integration). Generally, fault current values on the main feeders decrease while those near the DGs increase. To address this phenomenon, it is necessary to adjust relay settings to meet protection coordination requirements.

4.2.   Tracking route and Primary/Backup relay pair arrangement

As discussed earlier in this paper, the tracking route should be predefined based on the feeder topology to optimize the OCR parameters. The route is a set of relay numbers that are listed sequentially along the route. In each tracking route, the primary and backup protection relays are defined based on their upstream and downstream relationship within the route. For example, as shown in Figure 11 and Table 3, Relay 10 in Routes 5 and 6 is the backup protection relay for Relay 11 and 12. Relay 14 is the backup protection relay for Relay 15 and 16 in Routes 9 and 10. The 16-bus system (Case 1) was set up with ten tracking routes; each route can separate many primary and backup relay pairs. Then, the CTI value represents the difference in OT between two adjacent relays. The tracking routes and primary/backup relay pair arrangement for Case 2, which contains 26 routes, were also listed in Table 4 according to the network topology of Figure 12. The relay number of the second column of these tables are related to calculating the COT of Eq (4) for each case. Each combination of the third column of tables was used to evaluate the CTI of the primary/backup relay pair according to Eq (3) or Eq (7).

4.3.   Determination of optimal TMS and PCS settings using RIA-ATRM

The proposed RIA-ATRM was used to optimize the TMS and PCS parameters of each relay to correct the previously mentioned protection coordination violations. Tables 5 and 6 show a comparison of TMS and IP between the original design and the RIA-ATRM solution for Case 1 and Case 2, respectively. In both tables, the 'original design' represents the relay parameters from the initial design of the distribution system without considering DG integration. We apply the RIA-ATRM method to determine the optimal relay parameters for distribution systems under PR of 60% and 80%, with the results displayed separately in the two tables. It should be noted that the IP values shown in the tables are obtained by multiplying the optimal PCS, as determined by RIA-ATRM, by the relay's associated CT ratio, as calculated using Eq (2) or Eq (6).

The operating times of primary-backup relay pairs with their CTI are shown in Figures 20, 21, 22, and 24 for different DG penetrations. Figures 20 and 21, respectively, show the relay operating time and the corresponding CTI between relay pairs for the 16-bus system in Case 1 at 60% PR and 80% PR. The symbols, tb and tp, inside the figure are the operation time of the backup relay and primary relay, respectively. CTI can be evaluated by subtracting tb by tp. These values are calculated using the optimized relay parameters obtained through the RIA-ATRM method (refer to Table 5) with COT minimization and complied with relevant constraints, thereby correcting the violations of protection coordination. As mentioned in Section 2, the CTI constraint is set between 0.2 and 0.35, ensuring the smooth operation of the protection coordination mechanism; therefore, the figure clearly shows that the CTI value for each primary/backup relay pair is maintained at around 0.2.

This optimization model can similarly be applied to solve the optimal protection coordination problem for Case 2 with multiple DGs integrated. As shown in the bar charts in Figures 22 and 23, the OT for each relay and the CTI for each relay pair are derived using the IP and TMS parameters from Table 6. It can also be observed that the CTI values for all relay pairs remain within the specified limits.

4.4.   TCC verification of optimal relay parameter settings

In this section, the optimal relay parameters for Case 2 from Table 6 are utilized to simulate and generate TCC using ETAP software and are shown in Figures 24 and 25. The optimal relay parameters in Table 6 were obtained using the RIA-ATRM method to deal with the ill-coordination issues caused by DG integration, as shown in Figures 18 and 19.

The fault currents passing through each relay in Figures 24 and 25 are consistent with those in Figures 18 and 19; however, due to adjustments made to the relay parameters, the operation times of the relays will be affected. Therefore, there will be slight differences in the OT of the same relay between Figures 24 and 18 under the same fault current conditions, which can also be observed in a comparison between Figures 25 and 19. On the other hand, these two figures also show that the CTI between each relay is within 0.2 to 0.35 seconds, complying with the CTI requirements for protection coordination outlined in this paper. Additionally, the TCC for Relay 12 is no longer displayed in the figure because, after optimizing the parameters using the RIA-ATRM method, the fault current back-fed from DG1 and DG2 will no longer cause a malfunction in Relay 12 when the fault occurred at Bus 702.

4.5.   Convergence and robustness test of the RIA-ATRM

In this paper, we use the GA, PSO, IA, and the proposed RIA-ATRM to solve the protection coordination problem, considering the integration of the DGs to verify the performance and superiority of the proposed algorithm. These four approaches are applied to investigate the required adjustments to the TMS and PCS parameters of existing protection relays with the same test system. In this section, the COT mentioned in Eq (3) is the objective function that calculates the accumulative operation time summation according to the list of relays on the tracking routes, as shown in Tables 3 and 4.

The reverse fault current induces a decrease in the original fault current. Consequently, the TMS and PCS of the power relay need downward revision to shorten its operation time, thereby reducing the COT; conversely, if the fault current increases, it will lead to a longer operating time for the relay at the end of the feeder. Four heuristic search methods employing the same fitness function and constraint sets are evaluated. The convergence solutions are presented as the average values of each technique, derived from thirty executions of the corresponding program. Each execution involves 200,000 antibodies or chromosomes over 100 iterations. Table 7 shows the COT comparison for two study cases in two different PR. Compared to the other methods, the proposed RIA-ATRM exhibits superior search ability to find optimal TMS and PCS combinations to curtail the COT.

Additionally, we statistically analyze the optimal PCS and TMS settings obtained after running the RIA-ATRM program thirty times. We present these values as box-and-whisker plots in Figures 26 and 27. The plots show that the data from all thirty optimal solutions are highly concentrated, with no dispersion or divergence, indicating the robust stability of the proposed RIA-ATRM method in finding the optimal convergent solution for protection coordination parameters. The relay parameters in Figures 26 and 27 represent the optimal solutions obtained for Case 1 under the 80% PR condition, and then these parameters can be further used to calculate OT for each relay. Figure 28 shows the box-and-whisker plot of the relay OTs derived from the 30 sets of optimal relay parameters in Figures 26 and 27. However, for representation, we display only the OTs for the relays on tracking routes No. 2 and No. 3 as listed in Table 3.

We also applied the RIA-ATRM method to Case 2 under the 60% PR condition, performing thirty execution tests. The resulting thirty sets of optimal parameters and the calculated relay OT values are each presented as box-and-whisker plots, shown in Figures 29, 30, 31, and 32. Figures 29 and 30 also demonstrate that all thirty optimal solutions data are highly concentrated without dispersion and indicate the stability of the proposed RIA-ATRM algorithm. Figures 31 and 32 are box-and-whisker plots of relay OTs calculated from the thirty sets of optimal relay parameters shown in Figures 29 and 30. Here, only the OTs for the relays on tracking routes 4, 5, and 21, as listed in Table 4, as representative examples.

The convergence analysis of GA, PSO, IA, and RIA-ATRM in searching for the optimal TMS and PCS settings is illustrated in Figure 33. These metaheuristic search methods are evaluated using the same fitness function and constraint set. Convergence plots are generated from the average results of each technique over fifty program executions, each utilizing a total population of 200,000 chromosomes or particles across 100 iterations. Figure 33 highlights that the proposed RIA-ATRM achieves its best and most stable solution around the 37^th generation (or iteration), outperforming other methods in maximizing the fitness value. Additionally, it identifies the optimal TMS and PCS combinations for each relay, resulting in the approximate minimum total accumulated operation time.

5.   Conclusions

We present an RIA-ATRM approach to address the adaptive protection coordination challenge in MV distribution systems with high DG penetration. The protection coordination model is formulated as an optimization problem, offering two major contributions. First, we introduce the concept of tracking routes, dividing the distribution system into main feeder paths and branch line paths. Each tracking route set comprises the relay numbers installed on the transmission line topology. This allows for the calculation of relay operating times based on each relay's TMS and PCS, enabling the assessment of the suitability of each TMS and PCS setting combination. Second, we propose an RIA-ATRM algorithm featuring adaptive crossover and mutation operations to determine TMS and PCS settings for each relay on the tracking route, subject to CTI limitations. Simulation investigates a protection coordination problem of the 16-bus MV system and the IEEE 37 bus system to validate the proposed method. The results demonstrate that the proposed RIA-ATRM effectively reduces the COT and mitigates the impact on protection coordination settings under various conditions of DG integration.

Although the proposed algorithm can successfully handle the problem of relay protection coordination, according to the test of this article, there may be situations in the system that cannot be coordinated due to the inherent design nature. To address this issue, we suggest moderately loosening the set CTI limitation values while optimizing relay setting parameters. This adjustment facilitates algorithm convergence while maintaining protection coordination among relays.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in creating this article.

Acknowledgments

The author would like to thank the Financial support given to this work by the National Atomic Research Institute, Taiwan, R.O.C. under contract number NL1130113 is appreciated.

Conflict of interest

The authors declare no conflicts of interest.

Author contributions

T.-S. Z. provided project ideas and performed system models, data analysis, model simulations, article writing, and editing. Y.-D. L. and J.-L. J provided the research funding support and related experiences. All authors have read and agreed to the published version of the manuscript.