#### Abstract

The power system continuously deals with frequency fluctuations. When a power disturbance occurs, the transmission system operators rely on the underfrequency load shedding (UFLS) scheme to address severe underfrequency (UF) events to maintain the frequency at the permissible level and prevent blackouts. Defining settings of a conventional UFLS scheme is a very complex problem due to the nonlinear nature of the frequency response, and the size of the problem is vast because of the number of UF-relays spread on the power system. The under or over total load disconnection produced by the wrong setting of the UF-relays can create secondary frequency events or even a total blackout. This paper introduces a novel method to compute an optimally parametrized conventional UFLS scheme in specific given operating conditions by formulating it as a constrained problem and using the Improved Harmony Search (IHS) algorithm to solve it. Since there is no previous knowledge of using IHS to solve the UFLS scheme, a numerical parameter sensitivity analysis (PSA) is developed to tune the parameter of the IHS algorithm. The IEEE 39-bus system was modelled in DIgSILENT® PowerFactory™ and used as a test system. The optimally parametrized conventional UFLS methodology presented in this paper reveals superior results against the conventional UFLS scheme, and the suitability of using the IHS algorithm is confirmed.

#### 1. Introduction

The underfrequency load shedding (UFLS) is an emergency frequency control used in the power system to deal with underfrequency (UF) events caused by a deficit of power generation or an excess of power demand [1]. It is intended to re-establish the balance between power demand and power generation by disconnecting (shedding) an appropriate load amount to stop the frequency from declining [2]. The UFLS scheme uses multistage load shedding (LS) UF-relays designed to operate on the instantaneous value of system frequency [3]. The UF-relays operate any time the frequency drops below a frequency setpoint, and LS is done over several stages at the exact location where the frequency sensing is performed. The classical problem of specifying an automatic UFLS scheme is defining the fixed settings of the UF-relays [4]. The procedure to implement the conventional automatic UFLS scheme [5] based on local UF-relays consists of first identifying the worst possible power imbalance (commonly caused by the sudden disconnection of the most significant power generation unit) in the power system. Then, the total amount of LS that guarantees the frequency recovers above a minimum permissible value is estimated. After the total amount of LS is estimated, the fixed settings of the UF-relays are determined by using a trial-error procedure based on the utility company operators’ experience [3]. The process of setting the UF-relays using the conventional method becomes a complicated and easily mistaken task due to a large number of UF-relays installed in the power system. The wrong setting the UF-relays can cause two undesirable consequences:(i)An excessive LS, or overshedding, at the early frequency response stage, might result in overfrequency conditions or loss of the service continuity, negatively affecting the revenue of the power utility.(ii)An underestimated shedding, or undershedding, in the initial stages might result in failing to arrest the frequency decline, which may, in turn, lead to further loss of generation on underfrequency or even system-wide blackout.

A detailed survey of power system blackouts and cascading events are addressed in reference [6]. The semi-adaptative strategy, based on the rate of change of frequency, activates the UF relays when a frequency deviation exceeds a predefined threshold value. The amount of LS is in proportion to the frequency deviation [7]. Moreover, adaptative schemes involves the following: (i) event-based approach takes measurements of the tripped generator capacity using wide-area measurement systems (WAMS) and calculates the amount of LS [8], (ii) response-based method uses the swing equation to calculate the amount of load to be disconnected [9], and (iii) frequency prediction approaches use the predicted minimum frequency or the distance to the frequency threshold to activate the UF-relays [10]. Although these techniques offer advantages against the conventional method, they still can inaccurately estimate the amount of LS due to noise and/or delay in measurements, wrong disturbance estimation, and predictions errors, leading to an over/under LS.

In recent years, computational algorithms have been taking relevance to compute the optimal parameters of the conventional UFLS scheme. Several methodologies have relied on implementing artificial intelligence techniques and the use of WAMS [11] or have proposed appliances at the industrial level [12]. Moreover, other techniques formulate the UFLS as an adaptative scheme by applying artificial neural network (ANN) methods to compute the optimal parameters of UF-relays by determining the worst power imbalance [13]. In reference [14], an emergency LS model-based is proposed to minimize the total LS of the power system using trajectory sensitivity technique. Meanwhile, references [15, 16] use the particle swarm optimization (PSO) algorithm to calculate the optimal load disconnection considering the islanding operation conditions and multiple objective functions. Lastly, in reference [17], a model based on WAMS is employed to compute the optimal parameters of the conventional UFLS scheme using the imperialist competitive algorithm. For further information about the UFLS schemes, the authors of reference [18] have carried out a detailed and comprehensive survey of UFLS schemes available in the literature. After an extensive literature review, two main disadvantages have been raised: (i) Several methodologies such as in references [11, 13–16] only compute one optimal parameter of UF-relays (block size of LS) and the remaining parameters (number of LS stages, frequency setpoint, and time delay) are fixed values; and (ii) only one set of optimal parameters is computed, and this set of optimal parameters is the same for all UF-relays [11–13, 17]. These two drawbacks can produce an overestimation of the total amount of LS because the optimization problem is simplified, and therefore, the searching space is reduced. The metaheuristic algorithms adapted to solve the UFLS problem mentioned above offer satisfactory results; however, these algorithms are designed to solve unconstrained optimization problems. Consequently, dealing with the constrained UFLS problem using those methods implies additional work to turn it into an unconstrained problem, that is, selecting a suitable technique and adjusting its parameters [19]. Therefore, in this paper, the Improved Harmony Search (IHS) metaheuristic algorithm is used as an alternative to solve the UFLS problem since it is designed to handle constrained/unconstrained problems. The IHS algorithm is based on the harmony search (HS) algorithm proposed by Geem et al. [20]. It is focused on emulating the musical composition process to solve constraint/unconstraint optimization problems.

The purpose of this research paper is to calculate the UFLS scheme settings by formulating it as a constrained optimization problem, considering three parameters of UF-relays, which directly influence the frequency response to minimize the total amount of LS and guarantee the frequency of the power system. The proposed approach is designed for planning purposes using off-line simulations and classical UF-relays; no communication is needed. It allows the system operator to update the UFLS settings as many times as required. The main contributions of this paper are listed below:(1)The problem of the UFLS scheme is formulated as an optimization problem in a way that a set of optimal parameters (block size of LS, frequency setpoint, and the time delay) is computed for each UF-relay installed in the power system to minimize the total amount of LS.(2)The IHS algorithm has adapted to solve the optimal UFLS problem without turning it into an unconstrained problem.(3)A dedicated numerical parameter sensitivity analysis (PSA) is used to identify the appropriate parameters of the IHS algorithm to solve the optimal UFLS problem.(4)The IEEE 39-bus system was modelled in DIgSILENT® PowerFactory™. The IEEE 39-bus system model includes the dynamic of the governor and the automatic voltage regulator, and it was equipped with UF-relays provided in the DIgSILENT® PowerFactor™ library.(5)The proposed optimal UFLS approach is compared against the conventional UFLS scheme and solved by PSO algorithm, demonstrating superior performance.

This research paper is organized as follows: Section 2 provides a brief depiction of the conventional UFLS scheme and the methodology proposed in this paper, where Section 2.1 describes the formulation of the optimal UFLS approach; Section 2.2 is dedicated to present foundations of the IHS algorithm; Section 2.3 depicts the developed PSA procedure to tune the IHS parameters, and Section 2.4 describes the joint simulation outline designed to solve the optimal UFLS approach. Moreover, Section 3 provides the top results, which describe the test system used to assess the optimal UFLS approach, the numerical results of the PSA procedure (Section 3.1), and the optimal UFLS approach (Section 3.2). Finally, Section 4 presents the principal remarks and conclusions.

#### 2. The Optimally Parametrized Conventional UFLS Approach

The system frequency fluctuates second by second due to changes in production or demand for energy. Immediately after a severe disturbance produced by an outage of single/multiple generation units, the system frequency will decline. Depending on the generator prime-mover and spinning reserve, the frequency will eventually return to its desired value, which happens only if the loss of generation is lower than the spinning reserve [21]. However, UF-relays may initiate LS to stop frequency decline if the system frequency drops too low. The implementation of the conventional UFLS scheme relies basically on setting the values of four parameters at each individual UF-relay (N_{FR} = number of UF-relays): (i) the number of steps (*N*_{s}) or stages in which the UF-relay will shed the locally connected load, (ii) block size of LS (as a function of the total locally connected load) that will be disconnected at each stage at the *k*-th UF-relay (∆*P*_{s1,k}, ∆*P*_{s2,k}, …, ∆*P*_{sNs,k}), *k* = 1, 2, …, N_{FR}, (iii) the frequency setpoint (*f*_{s1,k}*, f*_{s2,k}*, …, f*_{sNs,k}) at which the load must be shed in each stage at the *k*-th UF-relay, *k* = 1, 2, …, N_{FR}, and (iv) the time delay (*t*_{ds1,k}*, t*_{ds2,k,}*…,t*_{dsNs,k}) between activating the consecutive stages (see Figure 1(a)).

**(a)**

**(b)**

Mathematically speaking, the k-th UF-relay is characterized by three vectors:where Δ*P*_{sj,k} is the block size of LS of the *j*-th stage at the *k*-th UF-relay, *f*_{sj,k} defines the frequency setpoint of the *j*-th stage at the *k*-th UF-relay, and *T*_{dj,k} is the time delay of the *j*-th stage at the *k*-th UF-relay.

The calculation of the parameters for each individual UF-relay using the conventional methodology becomes extremely difficult as the number of UF-relays installed in the power system increases. Considering a power system with N_{FR} UF-relays, the total number of parameters (*N*_{s}, Δ**P**_{k,}**F**_{k}, and **T**_{k}) to be calculated is (3*N*_{s} + 1) × N_{FR}. It is essential to disconnect the proper amount of load to ensure the frequency will recover into its permissible values and guarantee the security and economic operation of the power system [5, 22].

In addition, to compute the UF-relay parameters, the UFLS scheme must limit the maximum frequency deviation (Δf_{max}) to offer protection to larger turbine-generator units from continuous UF operation [23, 24] (i.e., 49.2 Hz at 50 Hz rated frequency). Moreover, the UFLS must limit the depth of the frequency response (*f*_{nadir}) for massive load levels and consider the system frequency will settle at a level depending upon the initial system overload and load reduction (Δ*P*_{L}) per frequency reduction (Δ*f*) [MW/Hz]. The idea is to provide time to other controllers to act and/or the system operator to survey the emergency condition and manually initiate additional LS and/or increase generation (such as peaking unit start-up). Furthermore, at the time of calculating the UF-relay parameters, it must be provided with an adequate boundary between the first frequency setpoint and the nominal frequency (*f*_{0}) to avoid the tripping on severe but nonemergency frequency situations. Therefore, it is required the UFLS scheme to be formulated as an optimization problem in order to calculate the optimal setting of UF-relays, ensure the frequency recovery, and at the same time avoid under/overload disconnections.

##### 2.1. Formulation of UFLS Scheme as an Optimization Problem

Nowadays, the widespread use of WAMS provides flexibility to transmission system operators for monitoring and controlling the power system. However, the use of WAMS also implies several challenges to deal with, such as data delay, data loss, and noise [25]. These challenges can affect the performance of the UFLS scheme if they are not proper to address. Thus, in this paper, the formulation of the optimal UFLS problem attempts to exploit the local measurements by including in addition to the ∆P other UF-relay parameters such as frequency setpoint (*f*_{s}) on the optimization problem. The location of the UF event directly impacts the local measurements taken by the UF-relay measurements.

###### 2.1.1. Decision Variables

The general problem of adjusting the UFLS scheme is defining the settings of each UF-relay, where (3*N*_{s} + 1) parameters are used at each UF-relay: *N*_{s}, Δ**P**_{k,}**F**_{k}, and **T**_{k}. Those parameters might be used as the decision-maker controls in the optimization formulation. The most common approach focuses on minimizing the total LS using ∆**P** as a decision variable. However, other approaches consider N_{s}, **F**, the frequency interval between two subsequent frequency setpoints (∆**F**), **T,** and even the best location of the UF-relays as decision variables. Therefore, in this paper, the optimally parametrized conventional UFLS problem is formulated considering ∆**P**, *f*_{s1}, ∆**F,** and **T** in each UF-relay as decision variables, and the vector of decision variables (**x**) containing the parameters of N_{FR} UF-relays is written as follows:where Δ**P**_{k} and **T**_{k} are defined in (1) and

Note that number of decision variables is reduced by one at each UF-relay (graphical interpretation is shown in Figure 1(b)) when the frequency interval between two subsequent frequency setpoints (∆**F**) is used. The decision variables are bounded by considering upper and lower limits as follows:

The superscripts *L* and U indicate the lower and upper bounds of the block size of LS (∆P), frequency setpoint (*f*_{s1}), and the frequency interval between two subsequent frequency setpoints (∆**F**), respectively.

The total number of decision variables of the optimal UFLS problem is *N* = 3N_{FR}N_{s}. Since the vector of decision variables has only considered the first frequency setpoint (*f*_{s1}) for the *k*-th UF-relay, the remainder frequency setpoints of the *k*-th UF-relay are calculated as follows (see Figure 1(b)):

###### 2.1.2. Objective Function

The purpose of solving the optimal UFLS problem is to determine the minimum amount of LS () required to recover the steady-state frequency into its operating values. Therefore, the objective function is formulated as the sum of the percentage of LS percentage of each activated stage in each UF-relay multiplied by the active power of its respective load. The objective function is defined as follows:where *P*_{load,k} is the total active power of the load in the *k*-th UF-relay, ∆*P*_{si,k} is the block size of LS (% of *P*_{load,k}) at the *i*-th activated stage of the *k*-th UF-relay, and *n*_{s} is the number of activated stages during the UF event (note that *n*_{s} ≤ *N*_{s}).

###### 2.1.3. Inequality Constraints

The main goal of solving the optimally parametrized conventional UFLS problem is finding the minimum LS () and ensure specific frequency requirements which are different in each country/utility company [2, 26]. In this paper, the frequency requirements are formulated as inequality constraints defined in function of the frequency setpoint of each stage (*f*_{s(i),k}) in the UF-relays and the steady-state frequency (*f*_{ss}).

###### 2.1.4. Steady-State Frequency (*f*_{ss})

After the action of the UF-relays, the steady-state frequency (*f*_{ss}) must remain inside its operative values to maintain the correct operation of the synchronous generators. Therefore, *f*_{ss} must be inside specific limits, and it is formulated as an inequality constraint as follows:where *f*_{0} is the nominal frequency of the power system and ∆*f*_{ss} is the permissible steady-state frequency deviation. The steady-state frequency requirement is considered in the optimal UFLS problem by evaluating that *f*_{ss} is above its minimum limit (*f*_{0}−∆*f*_{ss}). Therefore, a set of N_{FR} inequality constraints are collected in the vector **G***α*:

###### 2.1.5. Frequency Setpoint (*f*s)

The frequency setpoint in each stage of the *k*-th UF-relay (*f*_{s(i),k}, *i* = 1,…, *N*_{s}) must fulfil the technical requirements of the utility companies in which *f* shall be inside specific permissible values, that is, *f*^{L} ≤ *f*_{sik} ≤ *f*^{U} where *f*^{L} and *f*^{U} are the lower and upper limits of the frequency setpoint, respectively. Since the optimization problem is formulated to obtain the first frequency setpoint (*f*_{s1}) of each UF-relay, and it is already bounded, it is necessary to define the inequality constraint for the lower bound. A set of (N_{FR}N_{s}) equations are collected in the vector **G**_{β}:where the vector **J** is the vector of ones, that us, every element in that vector is equal to one. The vector of inequality constraints **G**(**x**) is used to consolidate all the previously defined inequality constraints as follows:and the UFLS optimization problem has a total number of inequality constraints *N*_{ineq} = *N*_{FR} (*N*_{s} + 1).

The formulation of the optimal UFLS problem is made in a way that it can easily add or remove parameters of the vector of decision variables. Therefore, according to user requirements, this proposed formulation can be adapted to the number of UF-relays and the number of stages at each UF-relay. The size of the vector of decision variables depends directly on the number of UF-relays to be set. Therefore, the complexity of the optimal UFLS problem increases as the number of UF-relays to be set rises. Moreover, since the optimally parametrized conventional UFLS problem is formulated as a constrained optimization problem, the authors have decided to adopt the IHS algorithm, which is a metaheuristic algorithm able to solve constrained optimization problems.

##### 2.2. Improved Harmony Search (IHS) Algorithm

The Improved Harmony Search (IHS) is a metaheuristic algorithm that follows the musical composition process of the IHS algorithm, the decision variables of the optimization problem represent the musicians' improvised pitches, and each solution vector is a harmony [20]. The procedure of the IHS is based on five main steps [27]:

###### 2.2.1. Initialize the IHS Parameters

The IHS algorithm has five main parameters that required to be initialized: (i) the harmony memory size (HMS) defines the size of the harmony memory matrix (HM) where the harmony vectors are stored; (ii) the harmony memory considering rate (HMCR) specifies the probability of whether new decision variables are selected from HM or created randomly; (iii) the minimum and maximum pitch adjustment rates (PAR_{min} and PAR_{max}) denote the probability of adjusting the decision variables; and (iv) the minimum and maximum bandwidth (*bw*_{min} and *bw*_{max}), and (v) the number of improvisations (NI) is the number of times the objective function is assessed.

###### 2.2.2. Initialize the HM

The IHS algorithm uses a set of harmony vectors created randomly to initialize the HM in the first improvisation.

###### 2.2.3. Improvisation of the New Harmony

The IHS algorithm uses HMCR to determine if the new harmony is selected from HM or is created randomly. The new harmony is chosen from HM using a harmony memory selection procedure if the probability is HMCR. On the other hand, if the probability is 1−HMCR, the new harmony is randomly created. When the new harmony is chosen from HM, it is adjusted considering the probability PAR; otherwise, the new harmony is not adjusted.

###### 2.2.4. Update the HM

The IHS algorithm updates the *HM* in each improvisation; it compares the worst harmony stored in the *HM* against the new harmony. If the new harmony is better than the worst harmony stored in the HM, the algorithm replaces the worst harmony with the new harmony. Otherwise, the HM does not change.

###### 2.2.5. Stopping Criteria

The criterion to stop the optimization process is based on the number of improvisations (NI). When the current improvisation is greater than *NI*, the IHS algorithm ends. On the other hand, if the current improvisation is smaller than *NI*, Step 2 and Step 3 are repeated until *NI* is reached.

As a summary, Figure 2 shows the flowchart of the IHS algorithm with the main steps involved.

##### 2.3. Tuning the Parameters of the IHS Algorithm

The performance of the IHS algorithm is impacted by the adjusting of its parameters (HMS, HMCR, PAR_{min}, PAR_{max}, *bw*_{min}, and *bw*_{max}). The HMS and HMCR parameters play an important role in the performance of the IHS algorithm. If the HMCR is very small, it produces a slow convergence of the objective function. Otherwise, if the HMCR is extremely large, it leads to converging at the local optimum point [28]. However, the worst-case scenario is a mistaken selection of HMCR that might lead to a potentially incorrect solution [29]. Meanwhile, a high value for the HMS parameter deteriorates the IHS algorithm’s performance; however, it must be high enough to maintain the diversity in the harmonies [28]. The values of the IHS algorithm parameters change depending on the optimization problem to be solved. Unfortunately, there is no sufficient evidence about the fine-tuning of the IHS algorithm parameters to solve the UFLS problem. Due to PAR and *bw* are already automatically adjusted over a range of values, a numerical parameter sensitivity analysis (PSA) is proposed to determine the best values of HMCR and HMS for solving the optimally parametrized conventional UFLS problem.

The PSA proposed in this paper focuses on obtaining a combination for HMS and HMCR, which produces the best performance of the IHS algorithm. The author decided to start with a default set of parameters of the IHS algorithms reported in the literature that allows the solution of an optimization problem. The default parameters of the IHS algorithm are taken from reference [30] (HMS = 1.0, HMCR = 0.80, PAR_{min} = 0.35, PAR_{max} = 0.99, *bw*_{min} = 1 × 10^{−5}, and bw_{max} = 1.0) as a base. Since HMS and HMCR are fixed values, a procedure to carry out the PSA is defined, and it is presented as follows:

*Step 1. *The range of values that the HMS and HMCR can take in the PSA is defined: HMS∈[ HMS_{min}, HMS_{max}] and HMCR∈[ HMCR_{min}, HMCR_{max}].

Each of those intervals is discretized over a specific number of equally speciated steps N_{HMS} and N_{HMCR} for HMS and HMCR, respectively.(a)The number of steps to define the HMS values can be any integer number, but this paper uses the size of the harmony (*N*) as an element of parametrization of HMS. The minimum limit is the default parameter (HMS_{0} = 1), HMS_{min} = 1, and the maximum limit is equal to the size of the harmony, HMS_{max} = *N.* The HMS interval is discretized over seven steps (N_{HMS} = 7). Therefore, the values to be evaluated are HMS∊ [1, (1/6) *N*, (2/6) *N*, (3/6) *N*, (4/6) *N*, (5/6) *N*, *N*].(b)The literature indicates that the value of the HMCR is typically selected between 0.7 and 0.95 for solving engineering problems [29]. Therefore, the limits of HMCR are selected as HMCR_{min} = 0.70 and HMCR_{max} = 0.95 with N_{HMCR} = 6 steps and the size of step (∆N_{HMCR}) is 5%, that is, ∆N_{HMCR} = 0.05. The values to be evaluated are HMCR ∊ HMCR_{min}[1, 1 + ∆N_{HMCR}, 1 + 2∆N_{HMCR}, 1 + 3∆N_{HMCR}, 1 + 4∆N_{HMCR}, 1 + 5∆N_{HMCR}].The values of the remaining parameters of the IHS algorithm are the same as the default values.

*Step 2. *Execute the IHS algorithm to evaluate all combinations of HMS and HMCR and obtain the objective function values () for a given number of improvisations.

*Step 3. *Identify the minimum value of the objective function () and compute the variation of the objective function through all combinations of HMS and HMCR concerning with as follows:

*Step 4. *The objective function final value can variate inside a range that does not produce an incorrect setting in the UF-relays, that is, the variation must be less than the lower bound of the block size of LS (∆*P*^{L}) in the UF-relay installed in the minimum load. Therefore, the maximum allowable variation of the objective functions is as follows:Select all pairs (HMS and HMCR) that produce a variation in the objective function () less than the maximum allowable variation ():

*Step 5. *Obtain the number of improvisations that the IHS algorithm requires to reach the maximum allowable value of the objective function, , for all .

*Step 6. *Select the pair (HMS and HMCR) that requires the minimum number of improvisations and has the minimum value of the objective function.

Figure 3 shows the flowchart of the procedure followed to tune the parameters of the IHS algorithm.

##### 2.4. Join Simulation Outline

The authors have developed a dedicated joint simulation outline for this research paper, and it is depicted in Figure 4. DIgSILENT® PowerFactory™ is used to perform time-domain simulations and produce time-series and discrete events. The high-level, general-purpose programming language Python is used to perform the optimization and to automatize simulation during the optimization purpose. The IHS is implemented in Python, and the optimization is performed in a close loop with DIgSILENT® PowerFactory™. For each improvisation, the IHS algorithm gives a vector of decision variables that contains the settings of the UF-relays. These settings are taken and put into UF-relays using the interface DIgSILENT® PowerFactory™-Python. Then, DIgSILENT® PowerFactory™ performs a time-domain simulation involving a series of discrete events (an outage of a synchronous generator and the activation of the UF-relays). The time-series results of the frequency and trigger signal of each stage of UF-relays are taken and sent back to Python and are used to evaluate the objective function defined in (5) and the constraints vector defined in (9). The optimization process ends when the *NI* is reached.

#### 3. Numerical Results

A classical simplified version of the New England power system in the USA, known as the IEEE 39-bus system, is used to illustrate the proposed approach. The test system consists of ten generators, 19 loads, 34 transmission lines, 12 transformers, and 39 buses; the full description of the model and the parameters is given in reference [31]; the total active power demand is *P*_{d} = 6097.10 MW. The dynamic model of the test system, including generation unit controllers, for example, governors and automatic voltage regulators (AVRs), is implemented in DIgSILENT® PowerFactory™ 2020. The UFLS scheme is based on a multistep UF-relay function ANSI 81 provided in the DIgSILENT® PowerFactory™ Global Library.

As the test system is representative of the New England area, the authors decided to use the South-eastern Electric Reliability Council (SERC) technical requirements to define the operational needs of the test system. The SERC UFLS Standard (PRC-006-SERC-02) establishes specific requirements for the UF-relay settings: (i) The steady-state frequency (*f*_{s}) shall be between 59.5 Hz and 60.5 Hz. The frequency setpoint (*f*_{s}) of the UF-relays shall be no lower than 58.4 Hz and not higher than 59.5 Hz, the interval between two subsequent frequency setpoints (∆*f*_{s}) shall be at least 0.2 Hz but no greater than 0.5 Hz, and the time delay (*t*_{d}) shall be at least six cycles (0.1 seconds) [32].

Starting from a steady-state operation, a frequency disturbance is inserted in the test system; it consists of the sudden disconnection of the generation unit G06 causing a 10.58% decrease of the total active power, and the total imbalance is Δ*P*_{T} = *P*_{G06} = 650 MW. The UFLS scheme is based on three UF-relays installed in the IEEE-39 bus system; the relays have been installed near to system frequency disturbance since UFLS should be performed as close as possible to the area deficient in the active power production [5]. The rationale behind implementing just 3 UF-relays is because this study does not intend to find the optimal location of the relays, just the optimal settings of the UF-relays available in an area of the power system. Nevertheless, the proposed methodology is readily expandable to larger power systems with a larger number of relays. Therefore, Relay 1, Relay 2, and Relay 3 are placed at buses 16, 21, and 23, respectively. The single-line diagram of the IEEE-39 bus system and the location of the UF-relays is depicted in Figure 5.

The number of UF-relays is N_{FR} = 3, the number of LS steps in each UF-relays is *N*_{S} = 3, therefore, the number of the decision variable is *N* = 2N_{FR} × *N*_{S} = 18, and the number of inequality constraints is *N*_{ineq} = N_{FR}(*N*_{S} + 1) = 12. The boundaries of the decision variables are defined following the requirements in the SERC UFLS Standard: ∆*P*^{L} = 1%, ∆*P*^{U} = 50%, *f*^{L} = 59.3 Hz, *f*^{U} = 59.5 Hz, ∆*f*^{L} = 0.2 Hz, and ∆*f*^{U} = 0.5 Hz. Furthermore, the time delay of all UF-relays is set constant and equal to *T*_{d} = 0.2 seconds. The nominal frequency is *f*_{0} = 60 Hz, and the permissible steady-state frequency deviation is *∆f*_{ss} = ±0.5 Hz.

##### 3.1. Tuning HMS and HMCR Parameters Using PSA

The fixed values for the IHS algorithm are *NI* = 500, PAR_{min} = 0.35, PAR_{max} = 0.99, *bw*_{min} = 1 × 10^{−5}, and *bw*_{max} = 1. The range of values of HMS and HMCR to perform the PSA is HMS ∊ [1, 3, 6, 9, 12, 15, 18] and HMCR∊[0.70, 0.75, 0.80, 0.85, 0.90, 0.95]. These values are obtained using Step 1 of the PSA procedure. The IHS algorithm is run N_{HMS} × N_{HMCR} = 42 times to evaluate the combination of the different sizes of the harmony memory with the probabilities of choosing the new decision variable from the harmony memory. The results of the evaluation of the HMS and HMCR are discussed below.

Figure 6 presents the convergence curves obtained during the PSA procedure, and evaluation of the objective function represents the total LS in MW during the UF event simulated in each improvisation through the optimization process. The objective function value in the first improvisation is different even if the HMS is the same; this situation is caused by the randomness included in the optimization method. At improvisation 500, for all pairs (HMS and HMCR), the objective function values vary between 288.5251 and 300.3221 MW, and the minimum value of the objective function is produced when HMS = 3 and HMCR = 0.85; therefore, (3,0.85) = 288.5251 MW.

**(a)**

**(b)**

**(c)**

**(d)**

**(e)**

**(f)**

An essential indicator of the performance of the optimization processes is the CPU burden; in this paper, it is obtained by measuring the simulation time (time elapsed during the optimization). Figure 7 presents the simulation time that is taken to execute the IHS algorithm over 500 improvisations for all pairs (HMCR and HMS). It can be observed that as the HMCR value approaches HMCR_{max}, the simulation time increases slightly; more precisely, it grows on average 1.08 minutes. Meanwhile, the simulation time increases on average 8.15 minutes as HMS increases. The correct tuning of the HMS will enhance the overall simulation time of the IHS algorithm for solving the optimal UFLS problem. Even though HMCR does not significantly impact the simulation time, it is essential to identify the best values of HMS and HMCR since it will allow solving the optimal UFLS problem in the shortest time. In addition, the optimal configuration of UF-relays is achieved without having under/over LS scenarios.

The maximum allowable variation of the objective function is computed using (11), as described in Step 4. The minimum load is in Relay 3 with 247.5 MW and ∆*P*^{L} = 1%. Therefore, the maximum permissible deviation of the objective function is . The value ensures selecting the best pair (HMS and HMCR) without causing a stage activation due to the variation of the objective function results since these variations are less than the lower bound of the block size of LS (∆*P*^{L}).

Following the PSA procedure, Table 1 shows the number of improvisations required to reach a value less or equal to the maximum allowable value of the objective function, which is computed using Step 5 and is = 291.0 MW.

The minimum number of improvisations that the IHS algorithm requires to solve the optimal UFLS problem is 92 for pair (1.0, 0.70), but this pair does not produce the minimum value of the objective function. Moreover, the literature advises using a value of HMS high enough to preserve the diversity in the harmonies [28]. On the other hand, the pair (18.0, 0.70) produces a value of the objective function lower that pair (1.0, 0.70) and requires 119 improvisations; however, the HMS is high, and it causes a significant increase in the simulation time. Therefore, the pair selected for the best setting of HMS and HMCR is (3, 0.85); this pair produces the objective function's minimum value and requires a relatively small number of improvisations compared to the total number of iterations (*NI* = 500) used in PSA. Results of the numeric PSA used to tune HMS and HMCR to solve the optimal UFLS problem rise the following remarks:(1)The precision of the HMCR value does not affect the simulation time considerably; however, if this parameter is wrong tuned, it can cause that the objective function converges slowly or converge to the local optimum point. Therefore, based on the PSA results, it is recommended to use a value of HMCR between 0.85 and 0.90 to solve the optimal UFLS problem.(2)It is recommended that the HMS be chosen between one-sixth and one-third of the size of the vector of decision variables since between those values, the IHS algorithm requires a small number of improvisations to reach the optimal solution and allows maintaining the diversity in the harmonies (decision variables). Moreover, the simulation time does not rise significantly regarding harmony memory size equal to one (HMS = 1.0).

For the solution of the optimal UFLS problem and based on the PSA carried out, the parameters of the IHS algorithm are tuned as HMS = 3.0, HMCR = 0.85, and *NI* = 250. The remainder parameter values are the same as the default values (PAR_{min} = 0.35, PAR_{max} = 0.99, *bw*_{min} = 1 × 10^{−5}, and *bw*_{max} = 1) taken from reference [30]. The IHS algorithm is executed 30 independent times, and its fitness statistics are compared against the ones obtained by executing the IHS algorithm using the default parameter values presented in the literature (HMS = 1, HMCR = 0.80). The fitness statistics of both scenarios are presented in Table 2. The graphical representation of the fitness statistics in each improvisation is shown in Figure 8.

**(a)**

**(b)**

The performance of the IHS algorithm is enhanced using the PSA to set HMS and HMCR instead of applying the default values. From the results depicted in Table 2, the dispersion of the fitness values is reduced by 31.2% (see the gray area in Figure 8), and the difference between the best and worst values is decreased by 16.88%. Moreover, using default values causes the IHS algorithm to require a greater number of improvisations to find the optimum solution. This is because the best value obtained using default parameters is greater than the one obtained using the tuned parameters in the same 250 improvisations.

##### 3.2. Assessing the Performance of the Proposed Optimally Parametrized Conventional UFLS Approach

The objective function is evaluated 250 times, and its convergence curve is presented in Figure 9. It shows that the IHS algorithm reaches the optimal solution around improvisation 141, and the value of the objective function is 287.7161 MW. The objective function value represents the total optimum amount of active power to be disconnected by the action of the UF-relays, that is, .

The solution vector of the decision variables (**x**) contains the settings of the block size of LS (∆*P*) for all stages of the UF-relays, the frequency setpoint of the first stage (*f*_{s1}) of each UF-relay, and the frequency interval between two subsequent frequency setpoints (∆*f*). From this result, the frequency setpoints of the remaining stages are computed using (4). The optimal settings of the UF-relays are presented in Table 3. The frequency setpoint (*f*_{s}) of each UF-relay satisfies the inequality constraint defined in (8), which establish that *f*_{s} must be inside given predefined limits, that is, 58.4 Hz ≤ *f*_{s} *≤* 59.5 Hz.

The optimal settings of the UF-relays are evaluated by performing the dynamic simulation of the IEEE-39 bus system when the sudden disconnection of generator G06 occurs. The frequency response is shown in Figure 10. After the action of UF-relays, the minimum frequency is 59.01 Hz, and it is reached at 8.705 seconds. The steady-state frequency is *f*_{ss} = 59.508 Hz, and this value is inside the operational limits, which has a permissible steady-state frequency deviation of *∆f*_{ss} = ±0.5 Hz, that is, 59.5 Hz ≤ *f*_{ss}≤60.5 Hz. Therefore, the steady-state frequency satisfies the inequality constraints defined in (7) to ensure the security of the power system.

The total amount of LS is Δ*P*_{LS} = 287.716 MW, and this value represents 44.26% of the total active imbalance (Δ*P*_{T} = *P*_{G06} = 650 MW) and is the amount of active power that the power system requires to recover the frequency into its predefined operational limits which are 59.5 Hz ≤ *f*_{ss} ≤ 60.5 Hz. The dynamic simulation results show that the setting of the UF-relays obtained from the optimization is correctly computed and satisfies the technical and operational requirements written as a set of inequality constraints (*fs* of each stage and *f*_{ss}). Moreover, the optimal settings of the UF-relays prevented the unnecessary disconnection of 362.2839 MW, which represent 5.94% of the total active power demand (*P*_{d} = 6097.10 MW).

The superiority of the proposed formulation to determinate the settings of the UFLS scheme is demonstrated by comparing it against the conventional UFLS scheme. Furthermore, to validate the suitability of the IHS algorithm to solve the proposed formulation for the optimal UFLS problem, its results are compared against those obtained by solving the problem using the PSO metaheuristic algorithm.

The parameters of the UF-relays for the conventional UFLS scheme are assumed the same on the three UF-relays; the block size of LS (Δ*P*) is 25%, 15% and 10%; and the frequency setpoint (*f*_{s}) is 59.3 Hz, 59.1 Hz, and 58.9 Hz, for stage 1, stage 2, and stage 3, respectively. This configuration has been taken according to the criteria recommended in reference [5]. Meanwhile, the proposed optimal UFLS problem was solved using the PSO algorithm, as explained in reference [15]. The optimal parameters of UF-relays are depicted in Table 4; the number of generations was the same as the one used to solve the IHS algorithm (250 generations).

The IEEE-39 bus system is excited by disconnecting the generator G06 producing a UF event. The frequency response indicators (*f*_{min}, *t*_{min}, and *f*_{ss}) [33] and the total amount of LS () are used to illustrate and assess the performance of the UFLS scheme using (i) the conventional UFLS approach, (ii) optimal UFLS approach solved by the PSO algorithm, and (iii) optimal UFLS approach solved by the IHS algorithm. The frequency response is depicted in Figure 11, and Table 5 summarizes the frequency response indicators as well as the total amount of load disconnected (**)**.

The main purpose of formulating the UFLS scheme as an optimization problem is obtaining the optimal parameters of the UF-relays that produce minimum load disconnection, at the same time arrest the frequency decaying and settle it into permissive operational values. From this perspective, the results presented in Table 5 certify the superiority of the optimal UFLS approach solved by the IHS algorithm because it prevents the unnecessary disconnection of 52.48 MW and 37.13 MW against the conventional UFLS approach and optimal UFLS approach solved by the OPS algorithm, respectively. Moreover, it ensures the frequency recovers inside its operational limits (59.5 Hz ≤ *f*_{ss} ≤ 60.5 Hz).

The suitability of the IHS algorithm for solving the optimal UFLS approach is demonstrated since for the same number of generations (improvisations), the PSO algorithm cannot find the optimum solution and the simulation time increases 89% concerning the simulation time taken by the IHS algorithm.

#### 4. Conclusions

An optimally parametrized conventional UFLS approach based on the block size of load shedding and the frequency setpoint was developed. This approach allows calculating the setting of the UFLS scheme, which produces the minimum load disconnection. The proposed numerical parametric sensitivity analysis procedure allows correctly determining the best values of HMS and HMCR to solve the optimal UFLS problem. The performance of the IHS algorithm is enhanced, creating a precedent of the HMS and HMCR values that should be used to solve optimal UFLS problem. Using the IHS algorithm to solve the optimal UFLS problem provides advantages, such as formulating a UFLS scheme as a constrained optimization problem since this algorithm can solve constrained optimization problems. Moreover, it shows superior performance in terms of simulation time and finding the optimal solution in contrast to the PSO algorithm. It demonstrated the superiority of the optimally parametrized conventional UFLS approach since it prevents excessive or insufficient disconnection of the load. Furthermore, the optimal setting obtained from the optimization fulfils the technical and operational requirements of the utility company and ensures the steady-state frequency recovering inside its allowable limits.

#### Data Availability

The simulation results used to support the findings of this study are available from the corresponding author upon request.

#### Conflicts of Interest

The authors declare that they have no conflicts of interest.

#### Acknowledgments

The research visit of J. M. Riquelme-Dominguez to the DIgEnSys-Lab has been supported by the research project funded by the Spanish National Research Agency/Agencia Estatal de Investigacion, grant nos. PID2019-108966RB-I00/AEI/10.13039/501100011033. Martha N. Acosta acknowledges the support of the University of South-Eastern Norway and Universidad Autónoma de Nuevo León and would like to thank the financial support provided by CONACYT (Mexico). F. Gonzalez-Longatt acknowledges DIgSILENT GmbH for supporting his research. This research has been funded by the Spanish Ministry of Economy and Competitiveness under the project TED2021-131311B-C21. Also by the Andalusian Regional Government under the project PYC20 RE 078 US.