A Robust Adaptive Overcurrent Relay Coordination Scheme for Wind-Farm-Integrated Power Systems Based on Forecasting the Wind Dynamics for Smart Energy Systems

: Conventional protection schemes in the distribution system are liable to su ﬀ er from high penetration of renewable energy source-based distributed generation (RES-DG). The characteristics of RES-DG, such as wind turbine generators (WTGs), are stochastic due to the intermittent behavior of wind dynamics (WD). It can ﬂuctuate the fault current level, which in turn creates the overcurrent relay coordination (ORC) problem. In this paper, the e ﬀ ects of WD such as wind speed and direction on the short-circuit current contribution from a WTG is investigated, and a robust adaptive overcurrent relay coordination scheme is proposed by forecasting the WD. The seasonal autoregression integrated moving average (SARIMA) and artiﬁcial neuro-fuzzy inference system (ANFIS) are implemented for forecasting periodic and nonperiodic WD, respectively, and the fault current level is calculated in advance. Furthermore, the ORC problem is optimized using hybrid Harris hawks optimization and linear programming (HHO–LP) to minimize the operating times of relays. The proposed algorithm is tested on the modiﬁed IEEE-8 bus system with wind farms, and the overcurrent relay (OCR) miscoordination caused by WD is eliminated. To further prove the e ﬀ ectiveness of the algorithm, it is also tested in a typical wind-farm-integrated substation. Compared to conventional protection schemes, the results of the proposed scheme were found to be promising in fault isolation with a remarkable reduction in the total operation time of relays and zero miscoordination. hybrid ANFIS–seasonal autoregression integrated moving average (SARIMA) forecasting algorithm coupled with the hybrid Harris hawks optimization and linear programming (HHO–LP) optimization algorithm.


Introduction
The integration of renewable energy source-based distributed generation (RES-DG), such as wind turbine generators (WTGs), in power systems is continuously increasing due to extensive technical developments, as well as clean and low-cost energy production [1,2]. The wind power share of the predicted FCL, the ORC problem is formulated as an NLP problem and is tackled by hybrid Harris hawks optimization and linear programming (HHO-LP). The main purpose of employing HH-LP is to achieve a remarkable reduction in the total operating time of relays by optimizing the relay settings to reduce the level of equipment damage and nuisance-related WTG disconnection, and to improve the reliability.
The key contributions of this paper are as follows: • The delay in updating the relay settings and the coordination with other relays can cause the malfunctioning of OCRs. A considerable delay time is evaded when updating the relay settings by predicting the wind speed and FCL variation in advance.

•
The hybrid ANFIS-SARIMA is devised for predicting periodic and nonperiodic wind series. • An efficient optimization model HHO-LP is established for the existing constraints. • A significant reduction in the overall tripping time of relays is achieved.

•
The is no record of miscoordination or limit violation.
The rest of the paper is organized as follows: Section 2 describes the problem formulation and objective functions for relay coordination. Section 3 addresses the proposed methodology. The investigated test systems and the obtained results are reflected in Section 4. Finally, Section 5 concludes the manuscript.

Objective Function
The fault current in a WTG WF varies with the variation in wind speed. This variation can adversely affect the ORC. Thus, the FCL needs to be determined in advance by forecasting the WD, which is coupled with an optimization algorithm (OA) for the overcurrent relay coordination problem. The OA reduces the total operation time of relays. The operation time of backup protection is of vital importance in protection coordination; thus, the time of operation for primary and backup relays is concurrently minimized.
where F OT is the overall operation time of relays for fault isolation, t p if and t b ijf represent the operation time of the i-th primary relay and j-th backup relay, respectively, for a fault at location f, and F, I, and J denote the set of fault points, the total primary relays, and the backup relays, respectively.
Normally, the OCR follows inverse time characteristics as follows: where t is the operation time of the relay, I is the fault current, A denotes the constant for relay TCCs, and B denotes the inverse time type. As standard definitions are considered in Equation (2), A, B, and C are assumed constant. The terms t p if and t b ijf can be expressed from Equation (2) as follows: Appl. Sci. 2020, 10, 6318 4 of 25 TMS p i and TMS b j are the time multiplier settings of the i-th primary relay and j-th backup relay, respectively. Similarly, I P i and I P j are the pickup currents of the primary and backup relays.

OCR Coordination Constraints
There should be enough time difference between the operating times of primary and backup relays, termed the coordination time interval (CTI), for the security of protection coordination.
The optimization variables TMS and I p are assessed within the lower and upper bounds given as I min p,i ≤ I p,i ≤ I max p,i .
The value of I p should be more than the maximum load current and less than the minimum short-circuit current value. During temporary faults, the relays should not operate and, therefore, there is a minimum limit for the operation time of relays.

EEMD
Empirical mode decomposition (EMD) is effective in extracting the characteristic information from an original wind speed series, which can be decomposed into a set of intrinsic mode functions (IMFs). The IMFs indicate the oscillatory mode of the original wind speed series. EMD is a self-adaptive time series processing method, which can be perfectly used for complicated processing [40]. The main drawback of EMD is its mode mixing problem. To resolve the mode mixing problem, the EEMD method was proposed in [41]. The procedure of EEMD is described as follows: (1) Initialize the number of ensembles (M) and the amplitude of the added white noise; set i = 1.
(2) Add a white noise series to the original wind speed series x(t).
x i (t) = x(t) + n i (t), (9) where n i (t) denotes the i-th added white noise series, and x i (t) denotes the series with the added white noise. (3) Decompose the series xi(t) into J IMFs cij(t) (j = 1, 2, . . . , J) using the EMD method, where cij(t) is the j-th IMF after the i-th trial, and J is the number of IMFs.

ANFIS
ANFIS is a multilayer feed forward network, which integrates the merits of neural networks and fuzzy inference systems [31]. In this paper, ANFIS with type-3 reasoning mechanisms was applied. The typical ANFIS with type-3 reasoning mechanisms consists of five layers, which are shown in Figure 1, the detailed descriptions of which are given in [31]. The functions of each layer are given below.
Layer 3: All the input variables are normalized in this layer, and the output of this layer is calculated as where O3,i is the output of Layer 3, and Wi is the incentive strength of rule i. Layer 4: The following node function is applied in this layer: where {pi,qi,ri} is the parameter set of the nodes. Layer 5: The single node in this layer summarizes all incoming series as follows: Figure 1. The architecture of the adaptive neural network-based fuzzy inference system (ANFIS) with type-3 reasoning mechanisms.

SARIMA
SARIMA is the most popular method for periodic time series prediction, which is described as follows: where x or y denotes the wind speed series, O 1,i is the membership degree of fuzzy set {A 1 , A 2 } or {B 1 , B 2 }, and µ(x) or µ(y) is the membership function.
The following membership function is utilized: where µ Ai (x) is the Gaussian function, and c i and σ i are the mean and standard deviation of the membership function, respectively. Layer 2: This layer is the operation layer. Layer 3: All the input variables are normalized in this layer, and the output of this layer is calculated as where O 3,i is the output of Layer 3, and W i is the incentive strength of rule i. Layer 4: The following node function is applied in this layer: where {pi,qi,ri} is the parameter set of the nodes. Layer 5: The single node in this layer summarizes all incoming series as follows: Appl. Sci. 2020, 10, 6318 6 of 25

SARIMA
SARIMA is the most popular method for periodic time series prediction, which is described as follows: where F(B) and U(B s ) denote nonperiodic and periodic autoregressive polynomials, respectively, Q(B) and V(B s ) denote nonperiodic and periodic moving average polynomials, respectively, Z t denotes the wind speed series, et represents the white noise series, d is the level of integration, D is the level of periodic integration, s is the order of periodicity, and B is the back-shift operator. More details about SARIMA can been found in [42].

Hybrid ANFIS-SARIMA Model
A composite algorithm comprising ensemble empirical mode decomposition (EEMD), ANFIS, and SARIMA can be utilized for short-term wind speed forecasting [43,44]. Due to the stochastic nature of WD, deep insight into the actual wind speed series (WSS) is paramount to get accurate results. The original WSS contains periodic and nonperiodic series. Thus, a composite method with the capability of modeling periodic and nonperiodic series is an efficient choice for wind-speed forecasting. The EEMD decomposes the original WSS into a set of intrinsic mode functions (IMFs) of periodic and nonperiodic series. SARIMA forecasts the periodic WSS and ANFIS forecasts the nonperiodic WSS. The artificial neural network-based ANFIS-SARIMA constitutes better nonlinear ability, adaptivity, associative learning ability, and fault tolerance. To predict the WSS data of the n-th interval, the data of the (n − 1)-th interval are also taken into account for better results. The inputs taken into account are local time, temperature, pressure, and humidity, whereas the outputs are wind direction, wind speed, and WT output power. The step-by-step procedure from WD forecasting to FCL calculation is described below.
(i) The WSS is decomposed into IMFs, and one residual series is given as where S(t) is the wind speed series, I i (t) represents the IMFs, and R n (t) is the residual series. (ii) The periodic and nonperiodic series of I i (t) and R n (t) are defined as P j (t) and N i (t), respectively.
Thus, the original wind speed series can be given as where N i (t) is the nonperiodic WSS, and P j (t) is the periodic WSS. (iii) For P j (t), the SARIMA model is implemented and the results are defined asP j (t), Whereas, for N i (t) and R n (t), the ANFIS model is implemented and the results are defined asN i (t) andR n (t).
The sum of results of ANFIS-SARIMA is the forecasted wind speed given aŝ (iv) On the basis of the predicted wind speed in Equation (20), the wind power can be expressed in the form of wind power flux or kinetic energy flux given as Appl. Sci. 2020, 10, 6318 7 of 25 where ρ (t) is the density of air, P (t) is the atmospheric pressure, R s is the specific gas density, and T (t) is the atmospheric temperature. A is the rotor swept area, C p is the coefficient of maximum power, and S(t) is the forecasted wind speed. If the wind hits the turbine at an angle ϕ (t) , as shown in Figure 2, then the azimuthal angle variation in the airflow can be considered as cosϕ (t) , and Equation (21) can be written as Appl. Sci. 2020, 10, x FOR PEER REVIEW 7 of 26 Usually, the quantity of interest is the temporal average of the power. In order to derive an expression for the temporal average of the power, we use Reynold's decomposition [45].
S S s and ϕ ϕ ϕ where ̅ (t) and (t) are the temporal means of the wind speed and wind direction, respectively, while s'(t) and φ'(t) are perturbations or fluctuations about their respective means. Hereafter, for simplicity the notation (t) is removed from all terms. Substituting Equation (24) into Equation (23) and performing Taylor's expansion and neglecting higher-order terms [46], we have where is the variance of wind speed and is the variance of direction.
(v) The squirrel-cage induction generator (SCIG) and doubly fed induction generator (DFIG) are used almost exclusively in the energy conversion stage of the induction generator wind power system. In this study, SCIG was used. The most commonly used system topology is an SCIG directly connected to the power grid, as shown in Figure 3. This topology implies a constant frequency and voltage of the SCIG that establishes a fixed-speed operation. In such a system, the SCIG relies on the grid (or capacitor bank) to provide reactive power, which is necessary to build electromagnetic excitation for the rotary field. The generating mode of the SCIG is triggered by driven torque, which acts opposite to the generator speed within the super-synchronous speed operation region. Due to the absence of a power electronics interface, such a system can only serve the grid support applications, wherein just limited control (pitch-angle control) can be applied. Usually, the quantity of interest is the temporal average of the power. In order to derive an expression for the temporal average of the power, we use Reynold's decomposition [45].
where S (t) and ϕ (t) are the temporal means of the wind speed and wind direction, respectively, while s' (t) and ϕ' (t) are perturbations or fluctuations about their respective means. Hereafter, for simplicity the notation (t) is removed from all terms. Substituting Equation (24) into Equation (23) and performing Taylor's expansion and neglecting higher-order terms [46], we have where σ 2 u is the variance of wind speed and σ 2 ϕ is the variance of direction.
(v) The squirrel-cage induction generator (SCIG) and doubly fed induction generator (DFIG) are used almost exclusively in the energy conversion stage of the induction generator wind power system. In this study, SCIG was used. The most commonly used system topology is an SCIG directly connected to the power grid, as shown in Figure 3. This topology implies a constant frequency and voltage of the SCIG that establishes a fixed-speed operation. In such a system, the SCIG relies on the grid (or capacitor bank) to provide reactive power, which is necessary to build electromagnetic excitation for the rotary field. The generating mode of the SCIG is triggered by driven torque, which acts opposite to the generator speed within the super-synchronous speed operation region. Due to the absence of a power electronics interface, such a system can only serve the grid support applications, wherein just limited control (pitch-angle control) can be applied.
SCIG relies on the grid (or capacitor bank) to provide reactive power, which is necessary to build electromagnetic excitation for the rotary field. The generating mode of the SCIG is triggered by driven torque, which acts opposite to the generator speed within the super-synchronous speed operation region. Due to the absence of a power electronics interface, such a system can only serve the grid support applications, wherein just limited control (pitch-angle control) can be applied.  Now that the mechanical power of the rotor is known, it must be determined how much of this power is transferred to the electrical grid. A simplified overview of energy transfer from wind to the electrical grid is shown in Figure 4. Now that the mechanical power of the rotor is known, it must be determined how much of this power is transferred to the electrical grid. A simplified overview of energy transfer from wind to the electrical grid is shown in Figure 4.
Here, the efficiency for the gearbox is 0.95, that for the generator is 0.97, and that for the power electronics device is 0.98 [47]. the current from the WTG in the wind farm can be calculated as where cos φ is the power factor, and VL is the line voltage. (vii) The fault current from a three-phase fault in a squirrel-cage induction machine is calculated using the network shown in Figure 5 [48]. The short-circuit current value at t = 0 is given as where E′ is the voltage behind transient reactance X′ , Rs is the stator resistance, Xs and Xr are the leakage reactance of the stator and rotor, respectively, Xm is the magnetizing reactance, and RL and XL are the resistance and reactance of the line connecting the WTG and grid. Substituting Equation (29) into Equation (28), the short-circuit current for a particular instance can be given as All the terms in Equation (31)   (vi) The electrical power transferred to the grid is given as P e = η gb η gn η p (P WT ). (26) Here, the efficiency for the gearbox is 0.95, that for the generator is 0.97, and that for the power electronics device is 0.98 [47]. the current from the WTG in the wind farm can be calculated as where cosφ is the power factor, and V L is the line voltage. (vii) The fault current from a three-phase fault in a squirrel-cage induction machine is calculated using the network shown in Figure 5 [48]. The short-circuit current value at t = 0 is given as where E is the voltage behind transient reactance X , R s is the stator resistance, X s and X r are the leakage reactance of the stator and rotor, respectively, X m is the magnetizing reactance, and R L Appl. Sci. 2020, 10, 6318 9 of 25 and X L are the resistance and reactance of the line connecting the WTG and grid. Substituting Equation (29) into Equation (28), the short-circuit current for a particular instance can be given as All the terms in Equation (31) are constant for an SCIG; thus, the fault current depends upon the value of I WTG calculated in Equation (16), which depends upon the wind speed. The total fault current from a wind farm is the sum of fault currents from all the WTGs.
using the network shown in Figure 5 [48]. The short-circuit current value at t = 0 is given as where E′ is the voltage behind transient reactance X′ , Rs is the stator resistance, Xs and Xr are the leakage reactance of the stator and rotor, respectively, Xm is the magnetizing reactance, and RL and XL are the resistance and reactance of the line connecting the WTG and grid. Substituting Equation (29) into Equation (28), the short-circuit current for a particular instance can be given as All the terms in Equation (31) are constant for an SCIG; thus, the fault current depends upon the value of IWTG calculated in Equation (16), which depends upon the wind speed. The total fault current from a wind farm is the sum of fault currents from all the WTGs.

Hybrid HHO-LP Optimization Algorithm
The simultaneous optimization of TMS and I p makes the ORC a nonlinear problem (NLP). A hybrid HHO-LP is proposed to solve this NLP by converting it into a linear programming (LP) problem. The basic technique involves the decomposition of the ORC problem into two subproblems. In the first subproblem, a random value is assigned to I p within its limits. This is only for the first iteration. Later on, its value is updated by the HHO. This converts the NLP into an LP. HHO calls the second subproblem in each iteration, which optimizes the TMS variable by using the standard LP method. This process continues until the convergence of the solution to an optimal value. Detailed descriptions of HHO and LP are given below.

Harris Hawks Optimization
HHO is a population-based algorithm proposed in [49], which mimics the hunting behavior of Harris hawks. It comprises exploration and exploitation. During the exploration phase, the position of hawks is updated on the basis of switches (ε) in attacking tactics as follows: where P (t) and P (t+1) are the position vectors of hawks at t and t + 1 iterations, respectively, P rand(t) is a randomly selected hawk position from the current population, P avg(t) is the average position of hawks, Pbest(t) is the prey position, LL and UL are the lower and upper limits of the position variables, and r 1 , r 2 , r 3 , and r 4 re random values selected from the range [0, 1]. The next stage is the transition from exploration to exploitation, which is executed by the change in different exploitative expressions, which depends on the escaping energy (E) of prey, as given in Equation (35).
where t is the current iteration, t max is the total number of iterations, and E 0 and E are the initial and current escape energies of prey randomly selected taken from [−1, 1]. During the attack of hawks in the exploitation phase, the prey has (r) probability of escaping. On the basis of the escape energy and the escape probability of prey, the hawks can espouse one of the four strategies tabulated in Table 1.
Hard siege (HS) Hard siege with progressive rapid dives (HSPRD) The positions of hawks are updated during SS and SSPRD using Equations (36) and (37), respectively.
where J is the random escape power of prey, ∆P (t) is the difference between the position vectors of the prey and hawk, and r 5 is a randomly selected number within the range [0, 1]. The positions of hawks are updated during HS and HSPRD using Equations (38) and (39), respectively.
where C and R are the values of current movements and rapid dives, LF is levy flight, and P m(t) is the average position of hawks.

Linear Programming
To convert the NLP into a linear one, the value of I p is fixed as extracted from HHO. The linear programming subproblem is called repeatedly by HHO to compute the value of TMS and the fitness of each hawk corresponding to the I p. A penalty relative to the severity of violation is added to the fitness value of each hawk if it violates the inequality coordination constraint. The complete algorithm of the proposed methodology for forecasting wind coupled with HHO-LP optimization algorithm is shown in Figure 6. In this study, short-term wind forecasting is used. The predicted wind speed, wind direction, and metrological variables are compared with the actual ones, and the predicted data are modified if the error exceeds the limits. On the basis of the predicted wind dynamics, the fault current is calculated in advance, and HHO-LP optimizes the ORC problem. The algorithm continuously checks if there is any change in wind data during a specified time interval; then, the fault current is calculated accordingly, and the relay settings are updated. The total time from forecasting to the upgrading of relay settings is given as where ∆T F is the time taken by ANFIS-SARIMA to calculate the predicted wind power and fault current, ∆T HHO-LP is the time consumed by HHO-LP to compute the optimum values of relay variables, and ∆T t is the time required to transfer the values of I p and TMS to the relay.

Case Studies
The impact of wind speed, wind direction, and metrological conditions on the variation in FCL in a wind-farm-integrated power system was observed by using the wind-speed prediction model, i.e., hybrid ANFIS-SARIMA. To increase the reliability of the power system, relay settings are updated continuously according to the predicted variation in FCL caused by wind speed. HHO-LP is incorporated to find the optimal values of TMS and Ip for a minimization of the total operation time of primary and backup relays for the fastest fault isolation.
The short-term interval (5 min) wind-speed data for the year 2019 collected from the southern parts of Pakistan, Jamshoro city, Sindh province, were used to train the ANFIS-SARIMA model. The latitude and longitude of this geographical location are 25.4007 and 68.2662, respectively. The predicted and actual wind speed, wind direction, and metrological conditions for one day in all four seasons are shown in Figures 7-9, respectively. It can be observed that the prediction of wind direction is harder than that of wind speed, as it requires extensive computation time.
It can be seen that the wind speed is highest during the spring season, resulting in a maximum chance of disturbing the protection coordination. The efficiency of the hybrid ANFIS-SARIMA model in terms of error compared with the state-of-the-art techniques reported in the literature such as fine

Case Studies
The impact of wind speed, wind direction, and metrological conditions on the variation in FCL in a wind-farm-integrated power system was observed by using the wind-speed prediction model, i.e., hybrid ANFIS-SARIMA. To increase the reliability of the power system, relay settings are updated continuously according to the predicted variation in FCL caused by wind speed. HHO-LP is incorporated to find the optimal values of TMS and I p for a minimization of the total operation time of primary and backup relays for the fastest fault isolation.
The short-term interval (5 min) wind-speed data for the year 2019 collected from the southern parts of Pakistan, Jamshoro city, Sindh province, were used to train the ANFIS-SARIMA model. The latitude and longitude of this geographical location are 25.4007 and 68.2662, respectively. The predicted and actual wind speed, wind direction, and metrological conditions for one day in all four seasons are shown in Figures 7-9, respectively. It can be observed that the prediction of wind direction is harder than that of wind speed, as it requires extensive computation time.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 12 of 26 neural network-radial basis function neural network(CNN-RBFNN) [52] and stacked extreme learning machine(SELM) [53] is reflected in Table 2.     neural network-radial basis function neural network(CNN-RBFNN) [52] and stacked extreme learning machine(SELM) [53] is reflected in Table 2.    neural network-radial basis function neural network(CNN-RBFNN) [52] and stacked extreme learning machine(SELM) [53] is reflected in Table 2.    It can be seen that the wind speed is highest during the spring season, resulting in a maximum chance of disturbing the protection coordination. The efficiency of the hybrid ANFIS-SARIMA model in terms of error compared with the state-of-the-art techniques reported in the literature such as fine tuning support vector machines(FTSVM) [50], modified persistence model(MPM) [51], convolutional neural network-radial basis function neural network(CNN-RBFNN) [52] and stacked extreme learning machine(SELM) [53] is reflected in Table 2.

Test System Specification
The IEEE-8 bus system and a real wind-farm-integrated substation named the Jhimpir power substation in Jamshoro, Pakistan, were considered testbeds for an evaluation of the performance of the proposed approach.
The ANFIS-SARIMA model and HHO-LP optimization technique were implemented in a registered version of Matlab-2020(R2002a). The fault analysis was carried out on the Electrical Analysis Transient Program (ETAP) software. The operating system had a processor Intel(R) Core™ i5-4210 central processing unit (CPU) at 1.7 GHz-2.4 GHz with 4 GB of random-access memory (RAM).

IEEE-8 Bus System
The standard IEEE-8 bus system is a meshed distribution system with multiple sources [54]. It consists of eight buses and seven line segments. Three-phase faults were simulated at the midpoint of each line segment. There were a total of 14 directional overcurrent relays (DOCR). These relays comprised 20 primary-backup relay pairs, as given in Table 3. To study the effect of wind speed on FCL variation, two wind farms were integrated at bus-3 and bus-6. The details of WTGs in the wind farms are given in Appendix A Table A1. A one-line diagram of the IEEE-8 bus system is given in Figure 10.

Test System Specification
The IEEE-8 bus system and a real wind-farm-integrated substation named the Jhimpir power substation in Jamshoro, Pakistan, were considered testbeds for an evaluation of the performance of the proposed approach.
The ANFIS-SARIMA model and HHO-LP optimization technique were implemented in a registered version of Matlab-2020(R2002a). The fault analysis was carried out on the Electrical Analysis Transient Program (ETAP) software. The operating system had a processor Intel(R) Core™ i5-4210 central processing unit (CPU) at 1.7 GHz-2.4 GHz with 4 GB of random-access memory (RAM).

IEEE-8 Bus System
The standard IEEE-8 bus system is a meshed distribution system with multiple sources [54]. It consists of eight buses and seven line segments. Three-phase faults were simulated at the midpoint of each line segment. There were a total of 14 directional overcurrent relays (DOCR). These relays comprised 20 primary-backup relay pairs, as given in Table 3. To study the effect of wind speed on FCL variation, two wind farms were integrated at bus-3 and bus-6. The details of WTGs in the wind farms are given in Appendix A Table A1. A one-line diagram of the IEEE-8 bus system is given in Figure 10.  GE Multilin750/760 numerical relays were employed here. The CT ratio for all relays was 400:1. The minimum and maximum limits of TMS were 0.05 and 1.1, respectively, whereas the limits for Ip were 1.1 × ILoad and 1.5 × ILoad. The CTI can take a value between 0.2 and 0.5 (REFF). However, for analysis of the IEEE-8 bus system, it was taken as 0.3. The minimum and maximum operation times of the relay were kept as 0.1 and 2.5, respectively, to assure the reliability of the proposed protection GE Multilin750/760 numerical relays were employed here. The CT ratio for all relays was 400:1. The minimum and maximum limits of TMS were 0.05 and 1.1, respectively, whereas the limits for I p were 1.1 × I Load and 1.5 × I Load . The CTI can take a value between 0.2 and 0.5 (REFF). However, for analysis of the IEEE-8 bus system, it was taken as 0.3. The minimum and maximum operation times of the relay were kept as 0.1 and 2.5, respectively, to assure the reliability of the proposed protection scheme.
Two case studies were carried out. In the first case, the wind speed and wind direction were taken as 10 m/s and 5 • , respectively, whereas, in the second case, the wind speed was 20 m/s and the wind direction was 10 • . Figure 11 depicts the variation in FCL due to wind speed for all 20 relay pairs for both cases. It is reflected that the fault current increased as the wind speed increased, as described in Equation (16). If the conventional relay settings are not updated according to wind-speed variation after WTG integration, then the variation in FCL due to WS variation can cause miscoordination in primary-backup relay pairs. Thus, using the predicted FCL by ANFIS-SARIMA, the relay settings were optimally updated by hybrid HHO-LP. The relay settings for a wind speed of 10 m/s obtained using conventional, PSO-LP [55], and HHO-LP methods are given in Table 4. Table 5 provides the time of operation of primary-backup relays for the first case. It can be seen that, in some cases, relay pairs violated the CTI limit. For example, in relay pair 6 in the first case, the backup relay R10 and primary relay R9 could not maintain a CTI of 0.3, thereby disturbing the protection coordination. On the other hand, in relay pair 10, both primary and backup relays operated at the same time. Figure 12a shows the characteristic curve of relay pair 6 taken during the simulation in ETAP, with conventional relay settings during the first case. The relay settings were updated optimally by HHO-LP on the basis of the predicted fault current, and the relay characteristics for the same relay pair 6 with the updated relay settings are shown in Figure 12b. It can be seen that the CTI of 0.3 was maintained between primary and backup relays for this pair. Two case studies were carried out. In the first case, the wind speed and wind direction were taken as 10 m/s and 5°, respectively, whereas, in the second case, the wind speed was 20 m/s and the wind direction was 10°. Figure 11 depicts the variation in FCL due to wind speed for all 20 relay pairs for both cases. It is reflected that the fault current increased as the wind speed increased, as described in Equation (16). If the conventional relay settings are not updated according to wind-speed variation after WTG integration, then the variation in FCL due to WS variation can cause miscoordination in primary-backup relay pairs. Thus, using the predicted FCL by ANFIS-SARIMA, the relay settings were optimally updated by hybrid HHO-LP. The relay settings for a wind speed of 10 m/s obtained using conventional, PSO-LP [55], and HHO-LP methods are given in Table 4. Table 5 provides the time of operation of primary-backup relays for the first case. It can be seen that, in some cases, relay pairs violated the CTI limit. For example, in relay pair 6 in the first case, the backup relay R10 and primary relay R9 could not maintain a CTI of 0.3, thereby disturbing the protection coordination. On the other hand, in relay pair 10, both primary and backup relays operated at the same time. Figure  12a shows the characteristic curve of relay pair 6 taken during the simulation in ETAP, with conventional relay settings during the first case. The relay settings were updated optimally by HHO-LP on the basis of the predicted fault current, and the relay characteristics for the same relay pair 6 with the updated relay settings are shown in Figure 12b. It can be seen that the CTI of 0.3 was maintained between primary and backup relays for this pair.       In the first case, with the conventional relay settings, the total operation time of primary and backup relays was 28.716 s and 43.541 s, respectively. With the relay settings obtained using PSO-LP [55], the overall operation time for the same primary and backup relays was 26.470 s and 35.568 s, respectively. Lastly, the time achieved upon optimizing the relay settings with the proposed HHO-LP was 23.93 s and 32.434 s, respectively, which is much lower than that obtained using the conventional and PSO-LP methods. In the second case, with a 20 m/s wind speed, HHO-LP reduced the total operation time of primary-backup relays by 19.86% compared to conventional settings and by 11.522% compared to PSO-LP. No miscoordination was reported with HHO-LP in all 20 relay pairs. The optimal relay settings and relay operating time for the second case are shown in Figures 13 and 14, respectively. The CTI between primary and backup relay pairs with conventional settings and those obtained using PSO-LP [55] and HHO-LP for both cases is reflected in Figure 15. It can be seen that, with the proposed approach, the CTI of 0.3 was maintained in all relay pairs, and the results were satisfactory. The proposed algorithm was also verified by changing the location and size of WTGs, and it is shown in Table 6 that the proposed HHO-LP worked well in all conditions. The overall performance in terms of an improvement in the reduction of time in all cases is given in Table 7.

. Jhimpir Wind-Farm-Integrated Substation
The Jhimpir power substation is a typical wind-farm-integrated substation in the Jamshoro city of Sindh province in Pakistan [56]. The geographical coordinates of this location are latitude = 24.4769 and longitude = 67.9240, and the hub height is 80 m. This wind farm consists of 31 wind turbine generators, each of capacity 1.6 MW. A doubly fed induction type generator (GE 1.6 MW/103) is installed in the WF. The cut-in, cut-out, and rated wind speed are 3.5 m/s, 25 m/s, and 12 m/s, respectively. The generated voltage is 0.690 kV, which is stepped up to 22 kV. Then, the 22 kV lines are connected to 132 kV lines with a step-up transformer of 22 kV/132 kV, 100 MVA. A one-line diagram of the wind-farm-integrated substation is given in Figure 16. This system has 36 overcurrent relays and there are 35 primary-backup relays pairs, as given in Table 8. The predicted and measured wind speeds for one day in all four seasons are already reflected in Figure 7. It can be seen that wind speed is high during the spring season; thus, it can have a high impact on the variation in fault current

Jhimpir Wind-Farm-Integrated Substation
The Jhimpir power substation is a typical wind-farm-integrated substation in the Jamshoro city of Sindh province in Pakistan [56]. The geographical coordinates of this location are latitude = 24.4769 and longitude = 67.9240, and the hub height is 80 m. This wind farm consists of 31 wind turbine generators, each of capacity 1.6 MW. A doubly fed induction type generator (GE 1.6 MW/103) is installed in the WF. The cut-in, cut-out, and rated wind speed are 3.5 m/s, 25 m/s, and 12 m/s, respectively. The generated voltage is 0.690 kV, which is stepped up to 22 kV. Then, the 22 kV lines are connected to 132 kV lines with a step-up transformer of 22 kV/132 kV, 100 MVA. A one-line diagram of the wind-farm-integrated substation is given in Figure 16. This system has 36 overcurrent relays and there are 35 primary-backup relays pairs, as given in Table 8. The predicted and measured wind speeds for one day in all four seasons are already reflected in Figure 7. It can be seen that wind speed is high during the spring season; thus, it can have a high impact on the variation in fault current level. Therefore, to validate the proposed algorithm in the Jhimpir wind-farm-integrated substation, one hour was selected from the day of the spring season, i.e., 10:00 a.m. to 11:00 a.m. on 15 July of 2019.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 18 of 26 the relay settings were determined earlier on the basis of the predicted fault current level. If the difference between the predicted and actual fault current due to wind-speed variation was greater than 2%, the actual values were then updated, and the HHO-LP algorithm was implemented during the current interval to optimize the relay settings on the basis of the actual FCL.    The minimum and maximum limits of TMS were 0.5 and 1.1, whereas, for I p , these limits were 1.1 × I Load and 1.5 × I Load respectively. The CTI, in this case was set to 0.3. The proposed algorithm was implemented, and the relay settings were updated after every 5 min. The TMS and I p for all relays using PSO-LP and HHO-LP updated at each interval are given in Table 9. The operation time of primary-backup relay pairs and the CTI for the 12 intervals of one hour are reflected in Figure 17. The results show that the proposed algorithm reduced the overall operation time of relays and also maintained the CTI between primary-backup relay pairs, which ensured the reliability and security of the power system. The overall performance in terms of an improvement in the reduction of time in all cases is given in Table 10. The computation time was reduced in the proposed algorithm because the relay settings were determined earlier on the basis of the predicted fault current level. If the difference between the predicted and actual fault current due to wind-speed variation was greater than 2%, the actual values were then updated, and the HHO-LP algorithm was implemented during the current interval to optimize the relay settings on the basis of the actual FCL.   TMS I p  TMS I p  TMS I p  TMS I p  TMS I p  TMS I p  TMS I p  TMS I P  TMS I P  TMS I P  TMS I

Conclusions
The disturbance in protection coordination caused by variations in the fault current level due to wind-speed variation in a wind-farm-integrated power system was investigated in this paper. The proposed algorithm consists of two steps. In the first step, the wind speed is predicted by forecasting it with hybrid ANFIS-SARIMA for a time interval of five minutes, and the fault current level is calculated in advance. Then, to reduce the overall operating times of relays in the second step, the hybrid HHO-LP is proposed for optimal relay coordination by optimizing the relay settings on the basis of the predicted fault current level. These optimal settings are transferred to relays at the start of each interval. If the difference between the actual and predicted value is more than 2%, then the HHO-LP is again run during the interval, and the relay settings are updated, which is a very rare case. The proposed algorithm was tested on a modified IEEE-8 bus system with WTGs and a local wind-farm-integrated substation. The results show that the proposed algorithm provided the optimal relay settings following the variation in fault current level due to wind-speed variation. No miscoordination was seen, and the proper CTI was maintained in all primary-backup relay pairs in all cases for both test benches with a considerable reduction in the overall operation time of relays, which shows the effectiveness of the proposed algorithm.

Conclusions
The disturbance in protection coordination caused by variations in the fault current level due to wind-speed variation in a wind-farm-integrated power system was investigated in this paper. The proposed algorithm consists of two steps. In the first step, the wind speed is predicted by forecasting it with hybrid ANFIS-SARIMA for a time interval of five minutes, and the fault current level is calculated in advance. Then, to reduce the overall operating times of relays in the second step, the hybrid HHO-LP is proposed for optimal relay coordination by optimizing the relay settings on the basis of the predicted fault current level. These optimal settings are transferred to relays at the start of each interval. If the difference between the actual and predicted value is more than 2%, then the HHO-LP is again run during the interval, and the relay settings are updated, which is a very rare case. The proposed algorithm was tested on a modified IEEE-8 bus system with WTGs and a local wind-farm-integrated substation. The results show that the proposed algorithm provided the optimal relay settings following the variation in fault current level due to wind-speed variation. No miscoordination was seen, and the proper CTI was maintained in all primary-backup relay pairs in all cases for both test benches with a considerable reduction in the overall operation time of relays, which shows the effectiveness of the proposed algorithm.
Funding: This research is funded by the Deanship of Scientific Research at King Saud University through research group No (RG-1438-089) and in part by the National Natural Science Foundation of China (NNSFC) through grant number 51707034. The authors also thank the Electrical Engineering Department, Southeast University and Lucheng Hong for the support given in conducting the research.

Conflicts of Interest:
The authors declare no conflict of interest.