Event-Based Under-Frequency Load Shedding Scheme in a Standalone Power System

Abstract

Under-frequency load shedding (UFLS) prevents a power grid from a blackout when a severe contingency occurs. UFLS schemes can be classified into two categories—event-based and response-driven. A response-driven scheme utilizes 81L relays with pre-determined settings while an event-based scheme develops a pre-specified look-up table. In this work, an event-based UFLS scheme is presented for use in an offshore standalone power grid with renewables to avoid cascading outages due to low frequency protection of wind power generators and photovoltaic arrays. Possible “N-1” and “N-2” forced outages for peak and off-peak load scenarios in summer and winter are investigated. For each forced outage event, the total shed load is minimized and the frequency nadir is maximized using particle swarm optimization (PSO). In order to reduce the computation time, initialization and parallel computing are implemented using MATLAB/Simulink because all forced outage events and all particles in PSO are mutually independent. A standalone 38-bus power grid with two wind turbines of 2 × 2 MW and photovoltaics of 7.563 MW was studied. For each event, the proposed method generally obtains a result with a smaller shed load and a smaller overshoot frequency than the utility and existing methods. These simulation results verify that the proposed method is practically applicable in a standalone power system with penetration of renewables.

1. Introduction

A power system must generate adequate power, such that the total power loads, losses, and required spinning reserves are met under normal conditions to maintain the system frequency at the rated value (60 or 50 Hz). If an outage of one or more large power generators occurs in the power grid, then a power mismatch between supply and demand is likely to occur, leading to a rapid decline in system frequency. A standalone or autonomous power system may experience a severe frequency drop if major power generation sources are lost on account of the small inertias of synchronous machines. Electric utility companies generally develop various protection schemes, such as under-frequency load shedding (UFLS), to deal with such severe events [1]. UFLS schemes can be classified into two categories—event-based and response-driven. When appropriate loads are shed, the declining frequency can be restored gradually to its nominal value.

Response-driven UFLS schemes use wide area measurement systems (WAMS) [2,3,4,5,6] or local real-time measurements [7,8,9,10,11,12,13,14,15,16,17,18,19] to activate a special protection system. Specifically, the prediction of the frequency nadir [2,3], the estimation of the rate of change of frequency (RoCoF) [4,5], and the optimization [6] were explored in WAMS-based UFLS schemes. Local measurement-based response-driven schemes incorporate digital relays, such as 81L and RoCoF relays, in which multi-stage frequencies and the time delay are set. Regarding UFLS schemes in a bulk power system using local measurements, minimization of shed load and maximum of frequency nadir in [7], the inertia time constant and generation deficiency in [8], and the dynamics of governors and load damping in [9] were implemented. Flexible random and fixed load priorities in [10], wind and photovoltaic (PV) measurements in [11], a closed-form-based shed load solution in [12], uncertainty of generation loss in [13], minimization the total amount of disconnected load using harmonic search in [14], and usage of a linear function for the determination of each individual load shed in [15] were investigated. In relation to response-driven UFLS schemes that have been developed for standalone/isolated power systems, a fuzzy rule base [16], the measured frequency deviation [17], an estimated/updated droop coefficient [18] and time-continuous load flow [19] were used to estimate the deficit power generation to prioritize the shed loads.

The event-based UFLS scheme uses a pre-determined look-up table to activate the shedding of low-priority loads as soon as a forced outage of one or more power generators occurs [20,21,22,23,24,25,26,27,28]. More specifically, two-loop integrated optimization was used to coordinate load shedding and corrective line switching under constraints of transient frequency/voltage deviation and overload capacity of the transmission lines [20]. A pre-specified UFLS scheme and post load shedding strategies were presented by considering the event type in [21]. The RoCoF and inertia constant were used to estimate the total power imbalance and a machine learning-based method was proposed to classify forced outage events by severity [22]. Shed loads were minimized subject to constraints on transient angle stability and voltage/frequency security [23]. Feedforward artificial neural networks with various PV penetration levels, total system load demands, and the frequency decay rates as inputs were proposed to determine the minimal shed load [24]. An online method for predicting and correcting potential load shedding by redistributing load dispatching through an artificial neural network was presented [25]. A fuzzy logic-based UFLS scheme was proposed to shed loads if a small hydro generator outage occurs [26,27]. A UFLS look-up table was developed using frequency information, RoCoF, customers’ willingness to pay, and loads histories in [28].

The works that were cited above [20,21,22,23,24,25,26,27,28] have at least one of the following limitations.

• (a). Although forced outage events were classified or identified, no optimization was implemented to mitigate the shed load [21,22,24,26,27,28], leading to large shed loads.

• (b). System responses, such as voltages and frequency, were used as event classifiers/identifiers [22]; however, load levels and the statuses of breakers at generators are the most intuitive signals for classification or identification.

• (c). Event-based UFLS schemes require many off-line time-domain simulations to design a look-up table or artificial neural networks [20,21,22,24,25,28]. Parallel computation is required to reduce the computational burden [23].

This paper proposes an event-based UFLS scheme for a realistic standalone power system. The contributions of this work are summarized as follows.

• (a). Particle swarm optimization (PSO) is employed to minimize the shed loads and maximize the frequency nadir without estimating power deficiency. This can be achieved by a co-simulation between a MATLAB code and Simulink time-domain simulation. Because the number of unknown variables is limited, the number of population size in PSO can be set at a small value to alleviate computational burden.

• (b). Only the status of breakers at the generators and the net power load, excluding renewable power generation, are utilized to identify the forced outage without the use of an extra detection or classification algorithm to reduce computation time. The aforementioned information is generally available in a dispatch center; matching data in a look-up table is faster than running complex detection or classification algorithms.

• (c). Parallel computation and initialization are implemented in the MATLAB/Simulink to reduce the computational burden. Specifically, the result of a normal base case scenario serves as an initial condition for the corresponding contingency cases to greatly reduce the time-domain simulation. Moreover, all contingencies and all particle studies in PSO are mutually independent; thus, they can be studied in parallel using multi-core in a CPU of a personal computer.

The rest of this paper is organized as follows. Section 2 describes and formulates the problem. Section 3 presents the proposed method for solving the problem, which involves PSO, parallel computation, and initialization. Section 4 presents the results of a simulation of a standalone offshore power system. Section 5 draws conclusions and provides directions for future work.

2. Description and Formulation of Problem

2.1. Difference between Event-Based and Response-Driven UFLS Schemes

Operation and dispatch in a standalone power system become increasingly challenging as the penetration of renewables increases. Specifically, the inertias of the power system become small owing to inverter control in distributed generations; the protection scheme must accommodate bi-directional power flows; and severe forced outages are likely to lead to cascading distributed generation outages.

The differences between event-based and response-driven UFLS schemes are summarized as follows.

• An event-based UFLS scheme has a single stage in which loads are shed immediately after a pre-specified outage occurs while a response-driven UFLS scheme comprises multiple stages with various frequency thresholds and time delays.

• An event-based UFLS scheme generally requires a communication system to identify or detect pre-specified events while a response-driven UFLS scheme needs local frequency measurements at load buses.

• An event-based UFLS scheme is implemented using a look-up table in an energy management system (EMS) while a response-driven UFLS scheme is realized using local digital relays, such as the 81L relay.

Further important considerations associated with event-based and response-driven UFLS schemes are as follows.

• The common communication system that is needed in the event-based UFLS scheme is ready in utilities; thus, no additional work is required. However, the latency that arises from the communication system must be reduced as much as possible. The triggering time for load shedding is typically 0.2–0.3 s, which includes local fault detection in protection relays and switch devices, and the communication latency from the control center to the breakers [23]. Thus, additional detection or classification algorithm must be avoided.

• Once a response-driven UFLS scheme is implemented, a latent cascading unanticipated outage at a renewable generation resource is most likely to happen. For example, the under-frequency protection of a Vestas V80 wind-turbine generator is set to 58.5 Hz. Before the breakers are open in the first stage (say 57.4 Hz) to shed loads, the wind-turbine generator will automatically trip, causing a more severe frequency nadir. The event-based UFLS scheme can avoid such an unanticipated wind-turbine outage.

• The frequency threshold in the first stage of a response-driven UFLS scheme relies on the power quality requirement in the power system. However, an event-based UFLS scheme activates pre-determined shed loads independently of frequency threshold. Restated, the moment of load shedding in an event-based UFLS scheme is earlier than that in the response-driven one. Accordingly, the frequency nadir that is obtained by an event-based UFLS scheme generally exceeds that by a response-driven UFLS scheme.

• The peak-load scenario is generally considered in developing a response-driven UFLS scheme; the load and power generation pattern in other scenarios are assumed proportional to those in the peak-load scenario to avoid many studies on multiple scenarios. Thus, only the peak-load scenario is studied in the response-driven UFLS scheme. However, this assumption may not hold in a power system with a high penetration of photovoltaic power and wind power because the power generation from renewables is extremely uncertain. To deal with this situation, multi-scenario considering stochastic (probabilistic) models shall be explored to develop the response-driven UFLS scheme.

2.2. Studied Standalone Power System

The studied Kinmen offshore standalone power system (24.26° N/118.20° E), operated by Taiwan Power Company, has 38 buses. The area of Kinmen is approximately 153.01 km² (https://en.wikipedia.org/wiki/Kinmen, accessed on 8 September 2021). The peak load in summer and the off-peak load in winter are estimated to be 66.62 MW and 22.8 MW, respectively. Figure 1 presents a one-line diagram of this standalone power system. The major diesel generators are located at buses 1~8 with capacities of 7.91~8.25 MW each. Buses 25 and 27 have minor diesel generators with a capacity of 3.168~3.512 MW each. Two doubly fed induction generators (DFIGs, Vestas V80, 2 MW × 2) are installed at buses 36 and 38. The total installed photovoltaic power is 7.563 MW, allocated at buses 14, 19, 29, and 34.

Table 1 shows the scheduled operation of generators at buses 1~8, 25, and 27, which was developed by the electric utility company according to various net load levels, excluding renewable power generation. Isochronous speed control (open loop) and droop control (closed loop) modes are implemented in the diesel generators at buses 1~8. The isochronous speed control maintains the speed of a synchronous machine at its nominal speed despite any small variation/imbalance of power in the power system. As a rule of thumb, when multiple diesel engines are connected at the same bus in parallel, the power system must be able to maintain a constant speed as long as the isochronous engine is capable of accommodating any disturbance. Restated, the power system is capable of maintaining a constant speed as long as the remaining isochronous engines can accommodate any power deficit that is caused by any N-1 or N-2 generator outage. Specifically, at least three on-line generators with isochronous engines are needed to stabilize the power system at buses 1~8. Generally, a droop governor lowers the speed from 3% to 5% of the reference speed [29]. The model of diesel engine governor can be found in the Appendix A.

2.3. Problem Formulation of a Single Event

The aforementioned problem in an offshore standalone power system can be formulated as the following multi-objective problem.

(1)$\mathrm{Min} F_1\left(\rho \right)$

(2)$\mathrm{Min}-F_2\left(x,y, \rho \right)$

(3)$\mathrm{s}.\mathrm{t}. {\rho }^{min}\le F_1\left(\rho \right)\le {\rho }^{max}$

(4)$\dot{x}=\mathrm{H}\left(x,y, \rho \right)$

(5)$y=\mathrm{G}\left(x,y, \rho \right)$

Specifically, $F_1\left(\rho \right)$ is the total shed load. ${\rho }^{max}$ and ${\rho }^{min}$ are set at 40% and 5% of the total system load, respectively. Generally, ${\rho }^{max}$ is set between 25% and 55% of the total system load [14,15,30]; the mean (40%) of 25% and 55% is adopted herein. ${\rho }^{min}$ is approximately set as the percentage of the largest single feeder load among all feeders. $F_2\left(x,y, \rho \right)$ is the frequency nadir; the minimum of the negative frequency nadir in Equation (2) is equivalent to its maximum. The symbol x denotes the vector of the dynamic state variables, such as the angles and velocities of synchronous machines. The symbol y is the vector of algebraic state variables, such as bus voltages. The total shed load in Equation (3) must be less than a given value, such as 40% of the total load. Equation (4) denotes a system of dynamic equations, while Equation (5) is the system static algebraic equation, involving an N-1 or N-2 contingency. Equations (1)–(5) constitute a multi-objective nonlinear programming problem that incorporates dynamic equations.

3. Proposed Method

3.1. Pareto Optimum

A multi-objective problem, such as Equations (1)–(5), has many Pareto optimums [31,32]. According to the concept of the least upper bound, Equations (1)–(5) are equivalently reformulated as [31,32]:

Min Z (6)

subject to

(7)$w_1\left[F_1\left(\rho \right)-F_1^{*}\right]\le Z$

(8)$w_2\left[F_2\left(x,y, \rho \right)-F_2^{*}\right]\le Z$

and Equations (3)–(6).

Where Z is the least upper bound in Equations (1) and (2). The least upper bound implies that the minimum of Z in Equation (6) is maximum limit in Equations (7) and (8). $F_1^{*}$ is set at $\frac{{\rho }^{max}}{4};$ $F_2^{*}$ is set at 60 Hz. The weighting factors $w_1$ and $w_2$ can be evaluated using Equation (9).

(9)$w_{\kappa }=\frac{F_{\kappa }^{max}-F_{\kappa }^{*}}{{{\displaystyle \sum }}_{\kappa =1}^2\left(F_{\kappa }^{max}-F_{\kappa }^{min}\right)}$

$F_1^{max}$ and $F_1^{min}$ are set at 40% and 5% of the total system load, respectively. $F_2^{max}$ and $F_2^{min}$ are set at 60.5 Hz and 57.4 Hz, respectively.

3.2. Particle Swarm Optimization

PSO is used to solve Equations (3)–(9) in this paper for two reasons. (a) PSO searches for the global optimum and outperforms traditional evolutionary algorithms (such as genetic algorithms) [33] and (b) PSO does not have to compute the gradient of implicit $F_2$, which is obtained herein by time-domain simulation using Simulink.

PSO uses the behavior of biological swarms to conduct a global search (exploration) and a local search (exploitation) to approach a target (global optimum). PSO involves a population of individuals (particles) that represent possible solutions and move in an $\mathsf{\Omega }$-dimensional search space. The velocity $V_p^t$ and position $X_p^t$ of particle p are updated using Equations (10) and (11), respectively [34,35]:

(10)$V_p^{t+1}={\omega }^t\times V_p^t+C_1\times r_1\times \left(p_{best}^t-X_p^t\right)+C_2\times r_2\times \left(g_{best}^t-X_p^t\right)$

(11)$X_p^{t+1}=X_p^t+V_p^{t+1}$

$V_p^{t+1}$ is the (t + 1)-th updated amount ($\Delta X_p^t$) to $X_p^t.$ In this paper, the number of dimensions $\mathsf{\Omega }$ is 2 and the 2-dimensional variable $X_p^{}$ includes $F_1\left(\rho \right)$ and Z only. Because the number of unknown variables is few, the number of population size is set at 10. The inertia weight ${\omega }^t$ allows the previously updated features to be reused in the next (t + 1) iteration. The random numbers $r_1$ and $r_2$ are between 0 and unity. Theoretically, learning factors $C_1$ and $C_2$ are positive constants constrained within the range from 0 to 2 such that $C_1$+ $C_2$≤ 4 to ensure the convergence of PSO. The values of $C_1$ and $C_2$ are set at 0.5 initially and then increased to unity at the maximum iteration herein. $p_{best}^t$ and $g_{best}^t$ are the best position of a particle in the t-th iteration and the best known position that is sought by any particle in the swarm so far, respectively. The initial shed loads are “30% plus a random number within [0, 0.5]” of the total system load.

The original PSO was developed for unconstrained optimization. However, Equations (3)–(9) constitute a single objective with nonlinear and linear constraints. To deal with these constraints, some treatments concerning the violation of constraints have to be implemented; otherwise, the algorithmic steps of PSO cannot converge at the maximum iteration and an infeasible result, which is not an anticipated solution for operators, is obtained. For example, if $F_1\left(\rho \right)>{\rho }^{max}$, this indicates that shed loads exceed the maximum allowed amount. The operators have to re-evaluate the scheduled on-line generators to enhance the system inertias. Table 2 shows these treatments. Specifically, if either Equation (7) or (8) is violated, then the least upper bound Z is set at the maximum of $w_1\left[F_1\left(\rho \right)-F_1^{*}\right]$ and $w_2\left[F_2\left(x,y, \rho \right)-F_2^{*}\right]$ to satisfy both Equations (7) and (8). Figure 2 plots relations among the MATLAB/Simulink modules that involve PSO-based optimization, the time-domain simulation (Simulink), and treatment of the violation of the constraints.

3.3. Initialization and Parallel Computing

Because PSO produces many particles, which are independent of one another, in each iteration, the time-domain simulation to obtain the individual $F_2\left(x,y, \rho \right)$ can be run in parallel. Moreover, all scenarios that involve different base cases and contingencies in the look-up tables are mutually independent. Time-domain simulation is known to require a very long CPU time, so proper initialization of system states for each particle is essential. In this work, the result of a normal base case scenario without generator outages is used as an initial condition for any corresponding N-1 and N-2 outages.

3.3.1. Initialization

To save a proper initial condition, the MATLAB command “set_param” is used to control parameters in Simulink; “SaveFinalState” provides an option for the final state of results while “SaveStateName” specifies the corresponding file name. The MATLAB command “sim” may be used to run Simulink; “SaveOutput” enables its output to be sent back to MATLAB while the option “StopTime” specifies the total simulation time. After Simulink is run, the option “SaveFinalState” is used to set to “off” to terminate the process.

To apply the aforementioned initial condition, the MATLAB command “set_param” can again be used to retrieve the saved results; “LoadInitialState” provides a means of retrieving the initial condition while “InitialState” specifies the corresponding file name. When a scenario with major generator outages is studied, the same MATLAB command “set_param” with the option “StopTime” can be used to set the total simulation time.

3.3.2. Parallel Computing

MATLAB provides a parallel computing function using multiple cores in a CPU of a personal computer. The MATLAB command “parpool” activates the available cores in a personal computer such that the developed parallel algorithm can be carried out independently in separate cores. The traditional command “for” to run a loop is replaced by “parfor”.

Parallel computing is very likely to lead to data concurrency errors because a single Simulink file produces different data saved in a file at the same time. To solve this problem, temporary directories may be opened for various cores at the beginning of the running of the MATLAB code. These temporary directories are removed after the corresponding Simulink task has been completed.

3.4. Block Diagram of Proposed Method

Figure 3 shows the block diagram of the proposed method. Block A identifies the N-1 and N-2 generator outages that may be likely to cause low system frequencies. The peak, off-peak and in-between base cases in the summer and winter seasons are identified. Different solar and wind power generations are estimated in block C. All information in blocks A, B, and C is considered to develop two sets of operating base case conditions for summer and winter, respectively, in block D. Each intact base case in block B is run using Simulink to obtain an initial condition for its corresponding contingencies in block E.

Each contingency that is associated with an N-1 or N-2 generator outage can be studied independently. In block F, Equations (10) and (11) are used to obtain each particle consisting of $F_1\left(\rho \right)$ and the least upper bound Z in PSO; in block G, Simulink is used to study the transient phenomenon in parallel, taking advantage of multi-core in a CPU. Finally, in block H, two look-up tables that are associated with shed loads for all scenarios that are related to summer and winter are established.

Figure 4 shows the flow chart of the proposed event-based UFLS method when a single scenario is considered. As shown in Figure 4, the net load level in a single scenario has to be estimated first to determine the control (droop or isochronous) mode and status (on or off) of each generator, which are indicated in Table 1. The corresponding contingencies for this scenario will be identified by checking their transient phenomena. The developed MATLAB/Simulink code considering inequality constraints (see Figure 2 and Table 2) can be then run in parallel. Finally, essential information about all contingencies can be output comprising the status of each generator, shed loads allocated at buses and net load level.

4. Simulation Results

4.1. Base-Case Scenarios and Contingencies

Different base-case scenarios in summer and winter were investigated herein. Table 3 and Table 4 show the load levels and their corresponding renewable power generations. Once the net load levels are known, power generations from the diesel generators are scheduled (see the second columns in Table 3 and Table 4). With these scheduled power generations, the corresponding diesel generators, as shown in Table 1, are set on-line to generate adequate power to balance the system load, considering isochronous speed control (open loop) and droop control (closed loop) modes. Table 5 and Table 6 show the on-line diesel generators for different base-case scenarios.

After the on-line diesel generators are specified, the definite N-1 and N-2 generator outages shall be studied to evaluate the frequency nadir. A total of 33 and 27 contingencies in summer and winter, respectively, shall be studied by running Simulink time-domain simulations. If the system frequency of a contingency is lower than 57.4 Hz, then the corresponding event-based UFLS scheme will be further studied. The predetermined shed load will be shed in the event-based UFLS scheme upon detection of the breaker status (on or off) of each diesel generator as well as the estimated net system load only without detecting the system frequency or other states.

Table 7 and Table 8 show the frequency nadirs of individual contingencies with N-1 and N-2 generator outages in the summer and winter seasons, respectively. The frequency nadirs that are caused by the generator outages at buses 25 and 27 are obviously higher than those at buses 1–8. Since proper isochronous and droop controls are implemented, as shown in Table 5 and Table 6, the system frequencies that arise from all contingencies are stable, as shown in Table 7 and Table 8. However, the event-based scheme is still required to prevent decline of the system frequency below 57.4 Hz. Fifteen out of 33 and 19 out of 27 contingencies require further bi-objective optimization by solving Equations (3)–(8) for the summer and winter seasons, respectively. Obviously, the conditions in the winter are worse than those in the summer.

4.2. Optimization of UFLS

The aforementioned 15 and 19 contingencies in the summer and winter seasons, respectively, are further investigated to minimize the total shed load and maximize the frequency nadir by PSO. Table 9 and Table 10 show the scenarios, contingencies, frequency nadirs (Hz) and the shed load (MW) for the summer and winter, respectively. The worst lowest frequency is 58.41 Hz, which greatly exceeds 57.4 Hz, after the shedding of 9.45 MW in response to an N-2 generator outage at buses 5 and 6 in summer. The worst lowest frequency is 57.61 Hz, which also exceeds 57.4 Hz, after the shedding of 8.54 MW for the same outage in winter. This contingency can be regarded as the most severe generator outage in both summer and winter when the off-peak load level is considered. Figure 5 displays the data flow, which consists of the net system load level and the status of the breaker of each generator, in a look-up table that is used to identify the specific contingency.

4.3. Comparison of Results Obtained Using Electric Utility Method, Singh’s Method and Proposed Method

The simulation results that are obtained using the proposed method are compared with those obtained using the original method of electric utility (Taiwan power company) and the Singh’s method in [22]. The UFLS scheme of the electric utility is a response-driven scheme, which covers four stages with frequency thresholds of 57.4, 57.0, 56.5 and 56 Hz in that order. The shed loads in the feeders at each stage are prioritized according to locations, such as residential area, hospitals, military bases, transportation hub, and a district office. The method in [22] is an event-based UFLS scheme, which used the RoCoF and inertia constant to estimate the total power imbalance. Then the amount of total shed load is equal to the estimated total power imbalance minus the total spinning reserve.

The event-based UFLS scheme relies on a wide area measurement system. Additional delays that are caused by PMUs, data serializing output, the communication system, and the phasor data concentrator, are around 270 ms [23], which is the estimated latency to activate the shed load following the forced outage in the following simulation.

According to experiences of the electric utility, all non-critical loads are allocated at four stages to be shed in the response-driven UFLS scheme; specifically, the percentages of four stages to shed loads are 18.06%, 24.35%, 21.34%, and 26.81%, respectively. Only 9.44% of total load is critical and remains. Obviously, the response-driven UFLS scheme is improperly designed due to this excessive shed loads. Accordingly, four base case scenarios s1, s5, w1, and w5, associated with the most severe N-2 contingency at buses 5 and 6, serve as examples in this section. As shown in Table 11, more loads in off-peak load scenarios (s1 and w1) are shed than those in peak load scenarios (s5 and w5) using the utility scheme because off-peak load scenarios have low inertias of synchronous machines. Moreover, more than one stage is required in these off-peak load scenarios when the utility response-driven scheme is used. The shed loads in scenarios s1, s5, and w1 obtained by the proposed method are smaller than those obtained by the utility scheme. Only the shed load of scenario w5 that was obtained by the proposed method (10.81 MW) exceeds that (9.40 MW) obtained by the utility scheme. Table 11 also shows the frequency nadirs and overshoot frequencies that are obtained by the proposed method and the utility scheme for these four contingencies.

Figure 6 and Figure 7 plot the corresponding transient responses that are obtained using the Simulink time-domain simulation. The overall frequency responses that are obtained by the proposed method are better than those obtained using the utility scheme because (i) the load is shed by the proposed method before it is shed by the utility scheme, and (ii) the shed load and frequency nadir are obtained by PSO in the proposed method. As shown in Table 11, the proposed method sheds smaller loads than the method in [22], such that overshoot frequencies obtained by the proposed method are smaller (see Figure 6 and Figure 7). Detailed comparisons among three methods are discussed below.

For the off-peak load in summer, two stage loads of 13.31 MW are shed if the utility scheme is implemented, leading to a very low frequency of 56.92 Hz and a large overshoot frequency of 63.35 Hz, as shown in Figure 6a. Moreover, the method in [22] obtained a shed load of 11.49 MW and an overshoot frequency of 61.08 Hz. When the proposed method is used, a smaller load of 9.45 MW is shed, leading to a much more moderate frequency response.

For the peak load in summer, one stage load of 11.32 MW is shed if the utility scheme is implemented; the method in [22] obtained a shed load of 12.12 MW; the proposed method, in contrast, sheds 10.33 MW. Although 11.32 and 10. 33 MW are very close, the load is shed by the proposed method before it is shed using the utility scheme, as presented in Figure 6b, so that frequency nadir (59.04 Hz) associated with the proposed method is higher than that (57.29 Hz) of the utility scheme. The moment that the method in [22] sheds the load is the same as that adopted by the proposed method. The frequency response obtained by the method in [22] becomes stable earlier than that obtained by the proposed method because more load is shed by the method in [22].

The case of the off-peak load in winter is the most severe among these four cases. The utility scheme sheds three stage loads of 13.62 MW, leading to the lowest frequency of 56.28 Hz and a very large overshoot of 62.19 Hz. It can be found that a very large overshoot of 61.35 Hz is obtained by the method in [22], too, due to a large shed load of 11.26 MW. However, the proposed method sheds a small load of 8.54 MW and yields a smooth frequency response, as shown in Figure 7a.

The utility scheme sheds one stage load of 9.40 MW if the peak load in winter is considered. The proposed method sheds a slightly higher load of 10.81 MW. The allover frequency response, which is between 58.94 and 60.29 Hz obtained by the proposed method, is much smoother than that obtained by the utility scheme, as shown in Figure 7b. The method in [22] shed the largest load among these three methods.

4.4. Discussions

According to the simulation results in Section 4.1, Section 4.2 and Section 4.3, four discussions are given below.

• (1). PSO is able to obtain smaller shed loads than those obtained by the method in [22]. PSO is also attains a smaller shed load in scenarios s1, s5, and w1 associating with the most severe N-2 contingency at buses 5 and 6 than the utility method. It can be found that the frequency nadirs obtained by the proposed method are higher than those by the utility method; the frequency nadirs obtained by the proposed method are very close to those by the method in [22]. All these evidences imply that PSO implemented in the proposed method is crucial.

• (2). According to Table 5 and Table 6, in additional to the net load levels in summer or winter, the number of status of breakers at generators, which are required in the look-up table, are few: (i) 8 at buses 1–8 and 3 at buses 25 and 27 in summer, (ii) 7 at buses 1~8 and 2 at buses 25 and 27. Accordingly, it is easy to implement the event-based UFLS scheme in a standalone power system without the use of a complex detection or classification algorithm.

• (3). When a MATLAB/Simulink time-domain simulation is carried out, it generally needs simulation time of 60 s to approach to a steady state due to a default initial condition. Thus, any occurrence of a forced outage may be set at 60 s and it takes more time, say 25 s (see Figure 6 and Figure 7), to approach to a new steady state. Thus, the optimal result of a normal condition obtained by PSO at 60 s may serve as an initial condition for all corresponding forced outages with the same net load level. A PSO considering 10 particles generally needs 20 iterations to attain an optimal solution. For the summer season, 15 contingencies require PSO studies (see Table 9). If all initial conditions of scenarios s1~s5 are ready, then total period of time-domain simulation is estimated to be 25 s × 10 × 20 × 15 (75,000 s = 20.83 h). If a 4-core CPU in a computer is used for parallel computing, then the period can be reduced to 18,750 s. However, if initial conditions and parallel computing are not implemented, then total period of the time-domain simulation is about 85 s × 10 × 20 × 15 (255,000 s = 70.83 h).

• (4). The proposed method can be applied to large-scale power systems with high penetration of renewable energies. Rather than diesel generators with isochronous control modes in a standalone power grid, combined cycle generators are used in a large power system to compensate the power deficit quickly. The average mode-based time-domain simulation tool (such as PSS/E) can be used to calculate the effective values of all state variables to save computational time. This time-domain simulation tool can be integrated with a Python-based PSO. Moreover, the net load levels, the number of status of breakers at generators and shed loads allocated at buses are required to build up a look-up table.

5. Conclusions

A novel particle swarm optimization based on the least-upper bound method is proposed for the development of an event-based under-frequency load shedding scheme in this paper. Initialization and parallel computation are performed to reduce the computational burden because all contingencies and all particles in PSO are mutually independent. Only the status of breakers at the diesel generators and the net power load, excluding the renewable power generation, are utilized to detect a forced outage without extra detection or classification algorithm. The net load level determines the on-line synchronous machines in the power system. A realistic 38-bus standalone power system is used to validate the proposed method, considering scheduled on-line generators with droop and isochronous control modes. According to the simulation results, identifying droop and isochronous control modes is essential in developing the event-based UFLS scheme. The power system can maintain a power balance as long as the remaining isochronous engines can accommodate any power deficit caused by any N-1 or N-2 generator outage. This observation was not addressed in the previous works. Besides, the off-peak load level is generally severe, owing to its low inertias once an N-2 forced outage occurs. The proposed event-based load shedding scheme outperforms the current utility scheme and an existing method, which is based on estimation of power imbalance. This implies that the optimization is crucial to minimize the amount of shed loads and maximize the frequency nadir; excessive shed loads obtained by the utility scheme and the existing method also result in large overshoot frequencies. Future works will differentiate the total system load and renewable power generation as input data for an event-based UFLS scheme. Moreover, a hybrid event-based and response-driven UFLS scheme will be also explored by considering severe forced outages and high penetration of renewables.

Author Contributions

Conceptualization, Y.-Y.H.; methodology, Y.-Y.H. and C.-Y.H.; software, C.-Y.H.; writing—original draft preparation, Y.-Y.H.; writing—review and editing, Y.-Y.H.; visualization, C.-Y.H.; project administration, Y.-Y.H.; funding acquisition, Y.-Y.H. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by Ministry of Science and Technology (MOST) in Taiwan, grant number MOST 110-3116-F-008-001 and the APC was funded by MOST in Taiwan.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Data sharing is not applicable to this article.

Conflicts of Interest

The authors declare no conflict of interest.

Nomenclature

Abbreviations:
CPU	central processing unit
EMS	energy management system
ICON	isochronous control
MBASE	MW base
PELEC	electric real power
PMECH	mechanical power
PMU	phase measurement unit
PSO	particle swarm optimization
PV	photovoltaic
RoCoF	rate of change of frequency (Hz/s)
SBASE	apparent power base
UFLS	under-frequency load shedding
WAMS	wide-area measurement system
Indices:
p	index for a particle, p = 1, 2, …, P (P = 10)
s1~s5	scenario indices for summer
t	iteration index in PSO
w1~w5	scenario indices for winter
Parameters:
$C_1$ and $C_2$	positive learning factors between 0 and 2
$F_k^{*}$	the expected value for $F_{\kappa }$, $\mathsf{\kappa }$ = 1, 2
$F_{\kappa }^{max}$, $F_{\kappa }^{min}$	the maximum and minimum values of $F_{\kappa }$
$r_1$, $r_2$	random numbers within [0, 1] in PSO
$w_{\kappa }$	the weighting factor for $F_{\kappa }$, $\mathsf{\kappa }$ = 1, 2
ρ	total shed load, the same as $F_1\left(\rho \right)$
${\omega }^t$	the inertia weight in the t-th iteration for PSO
Variables:
$g_{best}^t$	the best position of a particle’s neighbors so far in the t-th iteration for PSO
$p_{best}^t$	the best position of a particle in the t-th iteration for PSO
$V_p$	$\mathsf{\Omega }$-dimensional velocity vector of particle p
x	the vector of the dynamic state variables, such as the angles and velocities of synchronous machines
$X_p$	$\mathsf{\Omega }$-dimensional position vector of particle p
y	the vector of algebraic state variables, such as bus voltages
Z	the least upper bound of the multiple objective values
Functions:
$F_1$	a scale function for the total shed load
$F_2$	the frequency nadir

Appendix A

The governors of diesel generators strongly affect system frequency. Figure A1 shows a governor model. The value of droop is 0.03–0.05 for a close loop operation (ICON = 0). When ICON = 1, the diesel generator is under isochronous speed control (open loop). Table A1 shows parameters of different generators.

Figure A1

Mathematical model of diesel engine governor [36].

[Figure omitted. See PDF]

Table A1

Governor Parameters.

Parameters	Buses 1–4	Buses 5–8	Buses 25, 27
T1	0	0	0
T2	0	0	0
T3	0	0	0
K	1.2	1.5	0.1
T4	0.2	0.8	0
T5	0.17	0.17	1
T6	0	0	0
TD	0.055	0.045	0.055
TMAX	1	1	1
TMIN	0.06	0	0.06
Droop	0.05	N/A	0.05
TE	0.15	N/A	0.15

Footnotes

Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Figures and Tables

Figure 1. Studied 38-bus standalone power system.

Figure 2. Relation blocks of proposed method.

Figure 3. Block diagrams of proposed method.

Figure 4. Flow chart of proposed event-based UFLS method.

Figure 5. Data flow used in look-up table.

Figure 6. Frequency responses obtained by proposed method, utility scheme, and method in [22] in summer: (a) off-peak load; (b) peak load.

Figure 7. Frequency responses obtained by proposed method, utility scheme, and method in [22] in winter: (a) off-peak load; (b) peak load.

Table 1

Scheduled operation of generators at buses 1~8, 25, and 27.

Net Load Level	Buses 25 and 27	Buses 1~4	Buses 5~8
	droop mode	droop mode	isochronous mode
	scheduled power × quantity	scheduled power × quantity
–20.4 MW	0	5.1 MW × 4
20.401–23.0 MW	2.6 MW × 1	5.1 MW × 4
23.001–31.1 MW	2.6 MW × 1	5.7 MW × 5
31.101–41.6 MW	2.6 MW × 1	6.5 MW × 6
41.601–52.1 MW	2.6 MW × 2	6.7 MW × 7
52.101–62.2 MW	2.6 MW × 3	6.8 MW × 8
62.201–65.6 MW	2.6 MW × 4	6.9 MW × 8
65.601–68.2 MW	2.6 MW × 5	6.9 MW × 8

Table 2

Treatments of violations of constraints.

Condition	Treatment
$F_1\left(\rho \right)<{\rho }^{min}$	$F_1\left(\rho \right)={\rho }^{min}$
$F_1\left(\rho \right)>{\rho }^{max}$	$F_1\left(\rho \right)={\rho }^{max}$
Violation of Equation (7) or (8)	Z = maximum of $w_1\left[F_1\left(\rho \right)-F_1^{*}\right]$ and $w_2\left[F_2\left(x,y, \rho \right)-F_2^{*}\right]$

Table 3

Conditions (in MW) of base case in summer.

Scenario	Total Power of Diesel Generators	Load Level	Wind Power	PV Power	NetLoad Level
s1	31.101–41.60	33.36	4.0	0.0	29.36
s2	41.601–52.10	46.85	3.0	1.9	41.95
s3	41.601–52.10	57.15	2.0	3.8	51.35
s4	52.101–62.20	63.90	1.0	5.7	57.20
s5	52.101–62.20	66.62	0.0	7.6	59.02

Table 4

Conditions (in MW) of base case in winter.

Scenario	Total Power of Diesel Generators	Load Level	Wind Power	PV Power	NetLoad Level
w1	–20.04	22.80	4.0	0.0	18.80
w2	20.401–23.00	27.05	3.0	1.9	22.15
w3	23.001–31.10	36.35	2.0	3.8	30.55
w4	31.101–41.60	46.85	1.0	5.7	40.15
w5	41.601–52.10	55.53	0.0	7.6	47.93

Table 5

On-line diesel generators along with their MW generation and control mode in summer.

Scenario	Buses 1–4with Droop Mode	Buses 5–8 with Isochronous Mode	Buses 25 and 27with Droop Mode
s1	5.7 × 1 (MW)	5.7 × 4 (MW)	2.6 × 1 (MW)
s2	6.7 × 3 (MW)	6.7 × 4 (MW)	2.6 × 2 (MW)
s3	6.7 × 3 (MW)	6.7 × 4 (MW)	2.6 × 2 (MW)
s4	6.8 × 4 (MW)	6.8 × 4 (MW)	2.6 × 3 (MW)
s5	6.8 × 4 (MW)	6.8 × 4 (MW)	2.6 × 3 (MW)

Table 6

On-line diesel generators along with their MW generation and control mode in winter.

Scenario	Buses 1–4with Droop Mode	Buses 5–8 with Isochronous Mode	Buses 25 and 27with Droop Mode
w1	5.1 × 1 (MW)	5.1 × 3 (MW)	0 (MW)
w2	5.1 × 1 (MW)	5.1 × 3 (MW)	2.6 × 1 (MW)
w3	5.7 × 1 (MW)	5.7 × 4 (MW)	2.6 × 1 (MW)
w4	6.5 × 2 (MW)	6.5 × 4 (MW)	2.6 × 1 (MW)
w5	6.7 × 3 (MW)	6.7 × 4 (MW)	2.6 × 2 (MW)

Table 7

Individual frequency nadir (Hz) caused by various generator outages for each scenario in summer.

Scenario	N-1 atBus 1	N-1 atBus 5	N-2 atBuses 1 and 2	N-2 atBuses 5 and 6	N-1 atBus 25	N-2 atBus 25	N-1 at Bus 1; N-1 at Bus 5
s1	57.29	56.98	-	44.99	58.88	-	51.41
s2	57.50	58.15	55.28	55.92	59.05	58.17	55.64
s3	57.63	57.41	55.36	51.06	59.10	58.34	53.03
s4	57.66	57.92	55.72	55.76	59.16	58.32	55.70
s5	57.82	57.78	53.89	50.85	59.17	58.33	52.39

Table 8

Individual frequency nadir (Hz) caused by various generator outages for each scenario in winter.

Scenario	N-1 atBus 1	N-1 atBus 5	N-2 atBuses 1 and 2	N-2 atBuses 5 and 6	N-1 atBus 25	N-2 atBus 25	N-1 at Bus 1; N-1 at Bus 5
w1	56.58	55.94	-	44.99	-	-	45.15
w2	56.91	56.56	-	44.99	58.54	-	54.82
w3	57.39	56.98	-	44.99	58.92	-	58.92
w4	57.38	57.12	55.90	54.68	59.09	-	55.68
w5	57.57	57.62	54.46	50.40	59.08	58.27	52.46

Table 9

Frequency nadir and shed load for each contingency in summer.

Scenarios	Contingencies	Frequency Nadir (Hz)	Shed Load (MW)
s1	N-1 outage at bus 1	59.52	5.08
	N-1 outage at bus 5	59.41	6.62
	N-2 outage at buses 5 and 6	58.41	9.45
	N-1 at bus 1; N-1 at bus 5	58.48	9.19
s2	N-2 outage at buses 1 and 2	58.89	10.83
	N-2 outage at buses 5 and 6	59.21	8.38
	N-1 at bus 1; N-1 at bus 5	59.13	9.40
s3	N-2 outage at buses 1 and 2	59.18	11.46
	N-2 outage at buses 5 and 6	58.88	12.80
	N-1 at bus 1; N-1 at bus 5	59.05	11.84
s4	N-2 outage at buses 1 and 2	59.22	11.82
	N-2 outage at buses 5 and 6	59.17	12.10
	N-1 at bus 1; N-1 at bus 5	59.41	10.86
s5	N-2 outage at buses 1 and 2	59.41	11.32
	N-2 outage at buses 5 and 6	59.04	10.33
	N-1 at bus 1; N-1 at bus 5	59.26	12.06

Table 10

Frequency nadir and shed load for each contingency in winter.

Scenarios	Contingencies	Frequency Nadir (Hz)	Shed Load (MW)
w1	N-1 outage at bus 1	59.30	4.22
	N-1 outage at bus 5	59.24	4.61
	N-2 outage at buses 5 and 6	57.61	8.54
	N-1 at bus 1; N-1 at bus 5	57.77	8.10
w2	N-1 outage at bus 1	59.34	4.61
	N-1 outage at bus 5	59.28	4.67
	N-2 outage at buses 5 and 6	58.01	8.85
	N-1 at bus 1; N-1 at bus 5	58.16	8.38
w3	N-1 outage at bus 1	59.06	10.18
	N-1 outage at bus 5	58.98	9.89
	N-2 outage at buses 5 and 6	58.38	10.29
	N-1 at bus 1; N-1 at bus 5	58.55	9.57
w4	N-1 outage at bus 1	59.58	6.18
	N-1 outage at bus 5	59.43	7.16
	N-2 outage at buses 1 and 2	58.79	11.04
	N-2 outage at buses 5 and 6	58.66	11.73
	N-1 at bus 1; N-1 at bus 5	58.77	11.15
w5	N-2 outage at buses 1 and 2	59.04	11.46
	N-2 outage at buses 5 and 6	58.94	10.81
	N-1 at bus 1; N-1 at bus 5	59.12	11.10

Table 11

Comparisons of results (N-2 contingency at buses 5 and 6) obtained by electric utility scheme, method in [22] and proposed method.

Scenario	Methods	Shed Load(MW)	FrequencyNadir (Hz)	OvershootFrequency (Hz)	Number ofStages
s1	proposed	9.45	58.41	60.65	-
	utility	13.31	56.92	63.35	2
	[22]	11.49	58.60	61.08	-
s5	proposed	10.33	59.04	60.20	-
	utility	11.32	57.29	60.58	1
	[22]	12.12	59.08	60.37	-
w1	proposed	8.54	57.61	60.28	-
	utility	13.62	56.28	62.19	3
	[22]	11.26	57.98	61.35
w5	proposed	10.81	58.94	60.29	-
	utility	9.40	57.19	60.32	1
	[22]	11.56	58.95	60.44	-