Equilibrium solutions for microscopic stochastic systems in population dynamics

1.   Introduction

Integrating distributed generators into power systems has emerged as a transformative force in today's dynamic and ever-evolving energy landscape. Distributed generators, encompassing a wide range of renewable and non-renewable sources such as solar panels, wind turbines, diesel and microturbines, are increasingly prevalent in modern power networks. While these DGs promise cleaner and more sustainable energy generation, their integration presents unique challenges, particularly concerning the coordination of power system relays. Power relays play a critical role in safeguarding power systems by detecting and isolating faults, thereby ensuring the continuous delivery of electricity to consumers. However, traditional relay coordination practices face new complexities as DG penetration grows. This article delves into the intricate issue of power system relay coordination and explores the profound influence of DGs on this essential aspect of grid operation.

Power system relay coordination is a well-established discipline aimed at designing relay protection schemes that ensure the selective operation of protective devices during fault conditions. Traditionally, this coordination has focused on centralized generation sources, allowing for precise adjustment of relay settings based on a predetermined hierarchy. In this classical approach, generators are typically large, centralized units whose behavior can be predicted with a high degree of accuracy. Relay coordination is achieved by adjusting settings such as time-current characteristics, fault-clearing times, and coordination margins to ensure that the nearest relay responds to a fault while isolating it, leaving the rest of the system intact.

The integration of DGs challenges the traditional relay coordination paradigm in several ways:

Intermittency and variability: DGs, especially renewable sources like solar and wind, are inherently variable and intermittent. Their power output depends on weather conditions, time of day, and other factors, making it difficult to predict and coordinate relay settings accurately.

Bi-directional power flow: Unlike conventional generators, DGs can inject power into the grid and draw power from it. This bi-directional power flow introduces the potential for reverse power flows during fault conditions, which can complicate relay coordination strategies ^[1].

Islanding operation: DGs can create isolated sections or "islands" within the power grid during disturbances or outages. These islands pose unique challenges for relay coordination, as traditional settings may lead to unnecessary tripping of distributed generators.

Protection coordination time scales: DGs often operate on shorter time scales than centralized generators, responding rapidly to changes in grid conditions. Relay coordination must now consider these fast-acting devices while protecting against grid faults.

To deal with these challenges, power system engineers are developing innovative solutions considering the generators' distributed nature. Adaptive relay coordination schemes, which use real-time data and advanced algorithms, are gaining prominence. These schemes allow for dynamic adjustments of relay settings based on the actual operating conditions of the grid, including the presence and behavior of DGs. Adaptive relay coordination promises enhanced reliability and resilience in power systems, ensuring that protective devices respond appropriately to faults while minimizing unnecessary tripping of distributed generators. Moreover, these strategies can improve the overall efficiency and stability of the grid, making it more accommodating to the variability of renewable energy sources.

As mentioned, DG can impact the protection coordination of over-current relays on distribution systems. ^[2] proposed a "protection coordination index" (PCI), which can serve as an effective measure when planning the protection of meshed distribution systems with DGs. Then, a two-phase non-linear programming (NLP) optimization method is applied to determine the PCI by calculating variations in the maximum DG penetration level with changes in the protection coordination time interval. ^[3,4] resents the multi-agent theory to deal with protection coordination. In ^[3], the effects of DG grid connection have been considered, and the Java Agent Development Framework platform open software is used to adjust the protection coordination curve. ^[4] considers the issue of feeder reorganization when DG is connected to the grid and uses the signals measured by the phasor measuring unit (PMU) to dynamically adjust the relay settings by the central control system in the power station. The protection coordination problem was regarded as an optimization model in traditional mathematics, and the mixed integer nonlinear programming (MINLP) and mixed integer linear programming (MILP) were used to decide the TMS setting for fault current variation of the DGs integration in ^[5,6]. ^[7] uses the interior-point (IP) algorithm to solve the optimal adaptive protection coordination strategy under high penetration of green generation and network reconfiguration. Many contributions to the application of metaheuristic intelligence-based random search techniques applied to the adaptive coordination decision are represented here by document ^[8], which uses firefly optimization theory (FA) and evolutionary programming method (EP) to simultaneously consider feeder reorganization, DGs integration, DGs capacity planning and optimal adaptive protection coordination strategy. ^[9] utilizes the genetic algorithm (GA) to deal with the optimal DG placement problem to maximize the penetration level of DG in MV distribution networks without changing the original relay protection parameter. GA is used to discover the optimal sizes and locations of DG. Coincidentally, ^[10] also worked on the optimal placement issue of DGs and using optimization methodology to classify the cases, either coordination holds or coordination lost, by considering DG location changing and DC capacity variation. ^[11] presents an approach for protection coordination in microgrids that incorporates non-standard characteristic features of directional over-current relays (DOCRs), and the coordination model corresponds to an MINLP problem. There are four famous metaheuristic techniques, i.e., particle swarm optimization (PSO), GA, teaching-learning based optimization (TLBO) algorithm and shuffled frog leaping algorithm (SFLA), implemented for solving the optimal coordination problem. ^[12] proposes the use of adaptive protection, using voltage and current measurement of the system to overcome the challenges of overcurrent protection in distribution systems with DGs. The trip characteristics of OCRs are updated by detecting system operating states (grid-connected or island) and the faulted section.

References ^{[13,14,15,16]} apply artificial neural networks (ANN), using known training samples to train ANN weight values and construct a decision-making neural network system. When the system status changes–such as DG ON/OFF–the line after switching or line reorganization, the change in line current is used as the input value of the system and the corresponding relay setting adjustment value is output through the network decision-making process. In ^[13], an adaptive distance relaying method was presented. The multi-layered perceptron (MLP), also known as the error backpropagation network, was applied to estimate the actual power system condition and to calculate the appropriate tripping impedance under varying power system conditions. ^[14] and ^[15] apply ANN based on the Levenberg-Marquardt algorithm to model the DOCR characteristics where the trained ANN model can compute TMS, the pickup current setting (PCS) and the operation time of each DOCR in terms of changes caused by DG integration in the radial distribution system and looped distribution system respectively.

In this paper, the relay protection coordination issue was regarded as an optimization model by minimizing the objective function of the accumulated sum of all relays' operating times in the distribution power system subject to the TMS limits digital OCR and CTI restriction of relay operation time. The optimal protection coordination model is applied to a power system by considering several conditions of DG incorporation. According to DGs' various operating conditions, each OCR's TMS setting value is investigated. First, the author proposes the concept of tracking paths. Tracking paths can be divided into two categories: the main feeder path and the branch line path. The trace path will be represented as a set of relay numbers in order according to the topology of feeder connectivity. Then, the TMS of each relay can be used to calculate the OT of the relay, and the cumulative OT of the relay on each tracking path will be used to determine the fitness of the TMS setting combination. Furthermore, the relay protection coordination problem can be modeled to minimize the accumulated summation of all primary and backup relay OT subject to the CTI limitation.

Second, as the author's earlier contribution, an adaptive evolution mechanism with a self-adjusting scheme was proposed in ^[17]. The self-adjusting scheme was combined successfully with the PSO algorithm to improve the shortcomings of premature convergence. Then the modified particle swarm optimization algorithm with adaptive self-cognition and society operation scheme was proposed in this paper. The MPSO-ASSOS determines TMS for each relay on the tracking path, subject to CTI and OT limitation. A 16-bus test MV system with five DGs will be applied to test the adaptive protection coordination approach proposed in this paper. The simulation results show that the proposed method can quickly and stably obtain the solution with maximum fitness (i.e., minimum overall OT) and overcome the impact on protection coordination settings of the system with DG incorporation.

2.   Impact of DERs on RDNs protection coordination

As mentioned in the previous section, the power injection of RESs can affect the protection coordination of over-current relays (OCR) for the RDNs; there is a simple radial feeder, with six buses, five feeders, four DGs and five OCRs with paired circuit breakers (Figure 1).

The circuit breakers, CB 18, CB 19, CB 20, and CB 21, switch and simulate whether the DG is connected to the power distribution system. The fault analysis and protection coordination calculation of the demonstrated radial feeder was performed by using ETAP software. The parameter setting of each relay of the mentioned RDN in Figure 1 includes PCS and TMS, and these are assumed to carry out protection coordination calculation according to formal design steps in advance, as shown in Table 1. From Table 1, it can be found that the TMS of the relay is designed to be arranged in a sequential order of value, and the larger the value of the TMS, the longer the operation time of the relay. Therefore, it can think of an OCR with a larger TMS as the backup protection of an OCR with a smaller TMS, which is called primary protection; for instance, because Relay 8 is the primary protection of the line section 7–8, Relay 7 will be the backup protection of the Relay 8. Similarly, Relay 6 is the backup protection of Relay 7, and Relay 5 serves as the backup relay of Relay 6.

Table 2 shows the fault current distribution of the three-phase grounding fault at Bus 6 for the condition without DG injection and with DG1 integration. It can be seen from Figure 1 that only DG1 is incorporated into the end of this feeder. Table 2 shows that the source of the fault current flowing through Relay 5 and Relay 6 is the power source of this feeder, and due to the integration of DG1, it will generate an opposite-direction load flow from the power source. For this reason, the value of fault current through Relay 7 and Relay 8 when DG1 is integrated will be higher than if DG1 is not incorporated. Then, the higher fault current will affect the protection coordination mechanism of the original design. The protection coordination impact of DG1 integration for the fault that occurred at Bus 6 is shown in Figure 2. The operation time of Relay 5 and Relay 6 will not change due to the integration of DG1, which are 0.894 seconds and 0.646 seconds, respectively, as shown at points B and A. However, after DG1 is incorporated, the fault current value flowing through Relay 7 and Relay 8 becomes higher, so that the operation time of Relay 7 and Relay 8 can be calculated as 1.08 seconds and 0.557 seconds respectively, as shown at points D and C., according to the computed operation time, the order of operation sequence of Relay 5–8 will be changed to Relay 8, Relay 6, Relay 7 and Relay 5, which has caused the original protection coordination out of order. Simply put, when Bus 6 fails, Relay 6 should trip the circuit breaker, but due to the incorporation of DG1, Relay 8 will act instead and the protection coordination is lost.

A further example illustrates the situation when the power connection point of DG is not located at the end of the feeder. Considering the feeder in Figure 1, DG1 and DG2 are connected at Bus 8 and Bus 6, respectively, to the supply power at the same time. The resulting fault current distribution when the three-phase grounding fault occurred at Bus 5 is shown in Table 3. The fault current flowing to Bus 5 comes from DG1, DG2 and the power company, and it is not difficult to find that only the fault current flowing through Relay 5 is provided by the power company, regardless of whether the feeder is incorporated into DG or not. When DG1 and DG2 are combined, the fault current of Relay 6, Relay 7 and Relay 8 is increasing due to the contribution of DG1 and DG2. The protection coordination impact of two DGs integration for the fault that occurred at Bus 5 is shown in Figure 3. The operation time of Relay 5 remains the same value with the integration of DG1 and DG2, which are 0.854 seconds, as shown at point F. in Figure 3, however, after two DGs are incorporated, the fault current value flowing through Relay 6 to Relay 8 becomes higher, so that the operation time of three relays can be found as 1.03, 1.11 and 0.537 seconds, respectively, and are marked at points E, H and G in Figure 3. The order of operation sequence of Relay 5 to Relay 8 will be changed to Relay 8, Relay 5, Relay 6 and Relay 7, which has also caused the original protection coordination to be out of order. From the results of the above two case studies, it can be concluded that the existing relay settings need to be reviewed and updated to deal with the effect of DER integration to reduce the occurrence of relay miscoordination.

3.   Protection coordination optimization formulation

The relay protection coordination problem can be formulated as a nonlinear programming problem for minimizing all relay operating time summation. A relay's operation time depends on its PCS and TMS for a given fault current through the relay. Obtaining suitable PCS and TMS parameter values for all the relays installed in the distribution system and appropriately coordinating all the relays' operations is challenging ^[18,19]. To simplify the protection coordination optimization problem, this paper will assume that the PCS was given and has undergone a reasonably rigorous design and calculation process. Then, the setting of TMS can be obtained by formulating the OCRs coordination issue as an optimization problem with the objective function expressed as

where

Obj_F : Object function of calculating the total summation of operation time for each specified relay.

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

$ i $ : The relay identifier.

$ j $ : The fault location identifier.

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

FN : The total number of fault locations is considered.

TN : The total number of coordination tracking paths.

SN : The total relay number of the set of tracking paths.

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

$ {PCS}_{i} $ : The given 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 path.

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

Equation 1 is an objective function to minimize the summation of the operation times of the primary and backup relays when faults occur at different pre-setting positions. According to different fault locations, TN different protection coordination tracking paths can be determined and these paths may partially overlap, so the relay operation time on the overlapping paths will not be accumulated repeatedly. The normal inverse time-current characteristic function of the IEC standard ^[20], as shown in Eq 2, was only considered in this paper for evaluating the relay's operation time, and the fault current in Eq 2 can be simulated previously with several power engineering software, such as ETAP, PSS/E, etc. The parameters, $ {\alpha }_{i} $ and $ {\beta }_{i} $, of the standard normal time-current curve are set to the value of 0.14 and 0.02 in this paper, respectively.

Coordination time interval (CTI) is the operation time difference between the primary and backup relays. When the CTI of all the primary and backup paired relays in the power system is within a limit value or an allowable value range, it represents that protection relays are well-coordinated. The CTI evaluation and its limitations are clearly described in Eqs 3 and 4. Then, the selection of the TMS value will depend on the brand and model of the selected relay, and there are also restrictions on its upper limit and lower limit, as shown in Eq 5.

where

$ {CTI}_{i, j}^{k} $ : Object function of calculating the total summation of operation time for each specified relay.

$ {t}_{{i}_{p}, j}^{k} $ : Operation time of primary relay ip when the j^th fault occurs in the kth coordination tracking path.

$ {t}_{{i}_{b}, j}^{k} $ : Operation time of backup relay ip when the j^th fault occurs in the kth coordination tracking path.

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

$ {CTI}_{}^{max} $ : The upper CTI bound is set to 0.24 seconds in this paper.

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

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

4.   Solution algorithm-modified particle swarm optimization algorithm with adaptive self-cognition and society operation scheme (MPSO-ASSOS)

In the procedure of the PSO, a given objective function is optimized by simulating the natural behavior of birds (particles) flocking. Each particle (pbest) knows its best value and position so far; this information is the analogy to the personal experiences of each particle. Moreover, each particle knows the group's best value (gbest) among pbests so far. This information is an analogy to the knowledge of how the other particles around them have performed ^[21,22,23].

4.1.   Generation of particles and TMS initialization

4.1.1.   Simple self-cognition and society operation scheme (SSSOS)

In the original model of PSO, the velocity v will control the direction and distance of particle movement and can be classified into three operating modes as follows ^[21,22,23], and for the sake of distinction, it is called simple self-cognition and society operation scheme (SSSOS) here

1) Self-cognition operation mode

2) Society operation mode

3) Hybrid operation mode

where

rand₁ : The uniform random number in (0, 1) for the self-cognition process.

rand₂ : The uniform random number in (0, 1) for the society process.

g : The current iteration numbers.

$ {v}_{i}^{\left(g\right)} $ : The velocity of particle i at iteration g.

C₁, C₂ : Weighting factors.

$ {p}_{i}^{\left(g\right)} $ : Current position of particle i at iteration g.

$ pbes{t}_{i}^{\left(g\right)} $ : pbest of particle i at iteration g.

$ gbes{t}_{i}^{\left(g\right)} $ : gbest of particle i at iteration g.

A certain velocity, which gradually gets close to pbest and gbest, can be calculated using the above equation. The particle position of the next iteration (searching point within the solution space) with a unit time can be modified by the following equation:

4.1.2.   Adaptive self-cognition and society operation scheme (ASSOS)

The PSO algorithm prefers to use the hybrid operation to generate the next generation of new solution combinations in the SSSOS. A higher weighting factor, C₁, allows exploring solution territory around the local best solution and controlling the probability of introducing new particles. The weighting factor C₂ controls the velocity toward the area of the global best solution space. Then, the hybrid operation has the functions of the above two operations to generate the next generation of new solution combinations. If C₁ is too high, the solution might settle at a local optimum. On the contrary, a lower rate could generate too many possibilities of uncertainty. The offspring particle loses its resemblance to the parent particles, and the algorithm will not learn from the past and could become unstable. Choosing a suitable weighting factor of self-cognition and society operation for PSO is a dilemma. This article proposes an adaptive self-cognition and society operation scheme below to avoid such difficulty. The ASSOS scheme was illustrated as follows,

1) Assume that the operation mechanism of each particle generated in the last iteration process is recorded. There are three operation mechanisms: self-cognition, society and hybrid.

2) Select one particle orderly for producing a next-generation new particle according to

a) If rand₁ < P₁^(g) and rand₂ < P₂^(g): neither self-cognition nor society process be executed;

b) If rand₁ ≥ P₁^(g) and rand₂ < P₂^(g): execute self-cognition process only;

c) If rand₁ < P₁^(g) and rand₂ ≥ P₂^(g): execute society process only;

d) If rand₁ ≥ P₁^(g) and rand₂ ≥ P₂^(g): hybrid process be executed.

where

P₁^(g) : control parameter of self-cognition with initial value P₁⁽⁰⁾ = 0.5 and 0$ \leqslant $ P₁ $ \leqslant $ 1.

P₂^(g) : control parameter of society process with initial value P₂⁽⁰⁾ = 0.5 and 0 $ \leqslant $ P₂ $ \leqslant $1.

The next-generation new particles will be generated until all parent particles are processed. Figure 4 shows the initial relationship between self-cognition, society, and hybrid operation in PSO, which can be performed to generate next-generation particles in equal initial probability. If the random number of self-cognition(rand₁) or society(rand₂) processes is less than the corresponding control parameter, the related procedures will not be implemented. Self-cognition operation will play a more important role than that in PSO since the mutation can explore new regions. Suppose the search is very close to the local or global optimum. In that case, self-cognition may need to become dominant, especially in the absence of the critical suitable particles in a generation. Since all the procedures are randomly operated, there is no telling which is better.

3) A competition mechanism is thus implemented in the search process according to the fitness score. For instance, if the best current particle comes from a hybrid process, there is more likelihood for this procedure to generate a better solution for the next iteration. The area of the hybrid procedure must get bigger by reducing P₁^(g) and P₂^(g) to expand the probability, as shown in Figure 5.

If the best fitness of the iteration g-1 is greater than the generation g, i.e., $ {Fitness}_{min}^{\left(g-1\right)} > {Fitness}_{min}^{\left(g\right)} $ comes from hybrid procedure, both of control parameters will increase. We have

K₁ and K₂ are the regulating factors, and in general, K₁ < K₂. g_max is the maximum iteration number. Figure 5 shows the decrease in the probability of hybrid operation area. On the contrary, there is a greater likelihood for the other two procedures to generate better particles. Both control parameters must decrease to subjoin the probability, as shown in Figure 6.

If the best solution remains the same, hybrid operation needs to be held back to recover the related area. If the best fitness of generation g-1 is less than that of generation g, i.e., $ {Fitness}_{min}^{\left(g-1\right)}\le {Fitness}_{min}^{\left(g\right)} $ comes from the hybrid operation, the control parameters will decrease as the following manner,

4) If the best fitness of the iteration g-1 is greater than the iteration g, (i.e., $ {Fitness}_{min}^{\left(g-1\right)} > {Fitness}_{min}^{\left(g\right)} $) comes from only self-cognition procedure, then the control parameter will increase by using Eq 10. Conversely, if $ {Fitness}_{min}^{\left(g-1\right)}\le {Fitness}_{min}^{\left(g\right)} $ comes from self-cognition, then the control parameters will decrease by using Eq 12. In this situation, the control parameter P₂ is fixed. The probability variation of self-cognition is shown illustrated in Figures 7 and 8.

5) If the best fitness of the iteration g-1 is greater than the iteration g, i.e., $ {Fitness}_{min}^{\left(g-1\right)} > {Fitness}_{min}^{\left(g\right)} $, comes from only society procedure, then the control parameter will increase by using Eq 11. Conversely, if $ {Fitness}_{min}^{\left(g-1\right)}\le {Fitness}_{min}^{\left(g\right)} $ comes from society operation, then the control parameters will decrease by using Eq 13. In this situation, the control parameter P₁ is fixed. The probability variation is shown illustrated in Figures 9 and 10.

4.2.   Procedure of particles generation and TMS initialization

Particle generation is the first procedure of the PSO, and it has been mentioned in the previous section that the relay protection coordination problem can be regarded as a nonlinear programming problem. Therefore, in this paper, the MPSO-ASSOS will be used to combine the concept of nonlinear programming and apply it to the TMS parameters setting for each relay. First, the protection coordination tracking path in the system must be determined. The relay numbers on the path must be arranged in sequence according to the order of the main relay and the backup relay, and the set of relay numbers can be obtained; for instance, two tracking paths can be determined in Figure 1 as Path1: { 5 6 7 8 } and Path2: { 7 9 }. Path 1 can be called as the main feeder, and Path 2 is the branch separated from the main feeder. The TMS of relays on the main feeder and branch lines need to be designed in different manuals as explained in the following sections.

4.2.1.   TMS estimation of relays on the main feeder

As a main feeder path set example, Path: { x y z }, the last relay number z can be chosen as the starting point for protection coordination and its operation time can be set to the minimum value allowed by power regulations, such as 0.2 seconds. This relay can be regarded as the primary relay due to the fault current and PCS are given. Then, its TMS can be calculated using Eq 2 in the following way,

$ {t}_{z}^{k} $ is the presetting operation time of the starting point relay. Then, the $ {TMS}_{y} $ of the up-stream (backup) relay y will be estimated as

where UPB_y is the upper bound of TMS_y, it can be obtained by adding the maximum CTI allowable value to the operation time of the previous relay. Similarly, LOWB_y is the lower bound can be calculated by adding the minimum CTI allowable value to the operating time of the previous relay. The operation time $ {t}_{y}^{k} $ can be found by

Since relay x is the next up-stream (backup) protection of relay y, $ {TMS}_{x} $ can be calculated repeat according to Eqs 15–18, as long as the subscript y in these equations are changed to x.

For instance, the tracking path mentioned earlier in this section, Path1: { 5 6 7 8 }, if the operation time of relay 8, $ {t}_{8}^{k} $, is set to 0.2 sec., the TMS₈ can be obtained by using Eq 2 and Eq 14. Then, the TMS₇ could be generated randomly within upper and lower bound by

with

The operation time $ {t}_{7}^{k} $ can be found by

As mentioned above, the TMS of remaining rely in the path can be obtained by repeating Eqs 15–18.

4.2.2.   TMS estimation of relays on the branch line

After designing the relay parameters of the main feeder, the relay TMS of branch lines that branched from the main feeder can estimated next. Let us consider the continuation of the tracking path in section 4.2.1. and assume that the set of branch path is { y w }. This path means relay w is located on the branch line that is branching from main feeder after relay y. The TMS and operation time of relay w can be appraised by using Eqs 15 and 18 respectively, as long as the subscript y change to w and the upper and lower bound should be estimated as

After all the TMS setting values of all relays in the system be estimated, they can be used as a combination of possible solutions (a particle) and used as the basis for calculating the next-generation solution in PSO and MPSO-ASSOS.

4.3.   Fitness evaluation

The fitness score of each particle is obtained by calculating the objective function in Eq 1 and considering the inequivalent constraints mentioned in Eq 4 and Eq 5. If one or more variables of the particle violate their limits, the corresponding particle will be punished by multiply a punishing factor leading to a lower fitness value. The fitness function of the coordination problem is defined as

where

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

$ {PUNF}^{CTI} $ : The total summation of CTI punishing factor of all relays of each particle.

$ \rho, \varepsilon, \mu $ : Proportional constants.

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

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

$ {P}_{1} $ : TMS punishing value of relay i is set to more than 10³ when the TMS value is against the boundary.

$ {P}_{2} $ : CTI punishing value of relay i is also set to more than 10³.

5.   Simulation result and discussion

The proposed MPSO-ASSOS algorithm is simulated on 16-bus radial distribution systems. The proposed MPSO-ASSOS algorithm was applied to determine the optimal protection coordination issue on a 16-bus radial distribution system with DGs. The test system structure is shown in Figure 11; it has one main transformer, sixteen load buses and five DGs integration. The main substation consists of a 161 kV/22.8 kV, 50 MVA, Wye-Delta connection transformer. The DGs and load data are demonstrated in Table 4 and Table 5, respectively. The optimal protection coordination problem solved by MPSO-ASSOS was coded by MATLAB 2023 software. All programs were executed on a personal computer with Intel Core i7-10510U 2.3 GHz CPU and 16.0 GB RAM.

5.1.   Fault current calculation

The corresponding fault current, which flow through the specified relay under different fault, can be calculated from the ETAP software. In the simulation test of this paper, when no DG is integrated into the system, the ETAP is first used to calculate the fault current value flowing through each relay under conditions of different fault locations as the initial value of the relay's TMS setting. This test system will consider five DG-integrated operation scenarios, so there are eight operating conditions, including the above-mentioned initial operating states, as shown in Table 6. The fault current variation when the fault occurred at Bus 8, Bus 10 and Bus 14 for each DG operation condition are demonstrated in Figures 12, 13 and 14, respectively.

5.2.   Determination of coordination tracking path

According to the structure of the test distribution system shown in Figure 11, the coordination tracking path can be investigated and arranged as shown in Table 7. The Table lists 15 tracking paths. In addition to being used as the initial TMS solutions for estimating each relay, any two adjacent relays listed in each path also describe the connectivity relationship between the primary and backup relay. Therefore, the difference in operating time between the two relays is the CTI value.

5.3.   Total relay operation time minimization

The proposed MPSO-ASSOS and PSO original models will be used to solve the protection coordinating optimization problem considering the parallel operation of DGs. These two approaches were applied to investigate the required adjustments to the TMS parameters of existing protection relays with the same 16-bus MV system. In this section, the total relay operation time mentioned in Eq 1 is the objective function that calculates the total operation time summation according to the list of relay No. as shown in Table 7. Table 8 shows the total operation time of relays in the coordination tracking paths, which are evaluated from the optimal TMS setting value optimized by PSO and MPSO-ASSOS under different numbers of DGs integrated conditions, as shown in Table 6. It can be found from the data that as the number of DGs incorporated increases, the distribution of fault current changes. The reverse fault current causes the original fault current to decrease. The TMS of the power relay must be revised downwards, which shortens the operation time of the power relay so that the total operation time is reduced. Two heuristic search methods are tested using the same fitness function and constraint sets. Then, the convergence solutions were shown from the average value of each method after executing the corresponding program 30 times, with each time having 160,000 particles for 100 iterations. It is obvious that the proposed MPSO-ASSOS possesses a better search ability to obtain the suitable TMS parameter combination for the shorter total relay operating time than the PSO.

5.4.   Optimal TMS setting and corresponding CTI variation

In this study, all the relays are assumed to be numerical OCR, which considers continuous values of TMS and PCS. As mentioned earlier, the PCS of each relay was prespecified in the work. The optimum TMS settings of OCRs solved by the proposed MPSO-ASSOS are used to protect the network in each given condition, as shown in Table 9. These settings can provide proper decisions for TMS setting of protection coordination in any possible DGs integrated conditions.

The operating times of primary-backup relay pairs with their CTI are shown in Figures 15, 16, 17, 18 and 19 for different DG quantities mentioned in Table 6. Similar results of other DG combinations are not demonstrated because of the brevity. 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. It is to be noted that relay pairs 9 and 12 present the miscoordination situation in Figures 18 and 19 because of the high value of fault current flow through the backup relay, resulting in the CTI greater than 0.24 except for the above situations, and the remainder of these figure show that the CTI is between 0.2 and 0.24 sec. for all the relay pairs.

5.5.   Convergence test of the MPSO-ASSOS

The convergence of the MPSO-ASSOS, PSO and GA of searching for the best solution for optimal TMS setting is presented in Figure 20.

These metaheuristic search methods were tested based on the same fitness function of Eq 21 and constraint set, and then the convergence diagrams were drawn from the average value of each method after executing the corresponding program 50 times, with each time having total 500 chromosomes or particles populations for 100 iterations. Figure 20 demonstrates that the proposed MPSO-ASSOS can obtain the best and stable solution about in the 36th generation (or iteration) and that it performs better than the other two algorithms in finding the maximum fitness value. Moreover, it obtains the best TMS combination of each relay and the approximate lowest total accumulated operation time.

6.   Conclusions

This article proposes an approach to the adaptive protection coordination problem considering the MV distribution system with high DG incorporation. The protection coordination model was formulated as the optimization problem model. This paper included two main contributions: First, a concept of tracking paths is proposed, the distribution system has divided into main feeder path and branch line path. Each tracking path set contains the relay number installed on the transmission line topology. Then, the TMS of each relay can be used to calculate the OT of the relay to determine the fitness of the TMS setting combination. Second, an MPSO- ASSOS algorithm with adaptive self-cognition and society operation scheme was proposed, and it is utilized to determine TMS for each relay on the tracking path subject to CTI and TMS limitation. A protection coordination problem of a 16-bus MV system with five DGs and sixteen OCRs is applied to test the proposed method. The simulation results show that the proposed MPSO-ASSOS reduces the overall OT and overcomes the impact on protection coordination settings of the system with different conditions of the DGs integration.

Although the proposed algorithm can successfully handle the problem of relay protection coordination, according to the test of this article, there may still be situations in the system that cannot be coordinated due to the inherent design nature. This study recommends replacing all the existing electromechanical OCRs with modern numerical types. A central protection agent can apply the optimum TMS setting of the numerical OCRs to coordinate relays through telecommunication or wireless technology links.

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 NL1120262 is appreciated.

Conflict of interest

The authors declare no conflict of interest.

Author contributions

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