# Proactive and Automatic Underfrequency Load Shedding via PMUs and Particle
Filters
††thanks: This work was partially funded by the U.S. Department of Energy
under Contract DE-AC05-76RL01830.

Gian Paramo  Electrical and Computer Engineering
University of Florida
Gainesville, FL, USA
<EMAIL_ADDRESS>Arturo Bretas  Distributed Systems Group
Pacifc Northwest National Laboratory
Richland, WA, USA
<EMAIL_ADDRESS>Newton Bretas  Electrical and Computer Engineering
University of Sao Paulo
Sao Carlos, SP, Brazil
<EMAIL_ADDRESS>
###### Abstract

Underfrequency (UF) load shedding schemes are traditionally implemented in two
ways: One approach is based on manual load shedding, with system operators
requesting loads to be shed ahead of anticipated stressful operating
conditions. Manual load shedding is usually done through phone calls. The
second method is automatic load shedding via underfrequency relays. Using
static settings, these schemes can be designed to operate in stages and drop
previously identified loads. The main limitation of traditional load shedding
schemes is that they are reactive and leave little room for optimized
corrective actions. This work presents a proactive and automatic
underfrequency load shedding solution for power systems. Measurements are
captured via phasor measurement units (PMUs) at relatively low sampling rates
of 30 Hz. These measurements are then processed by particle filters who
predict the future state of the system’s frequency. Based on these predictions
excess load is determined and shed. Comparative case studies are performed in
simulated environments. Easy-to-implement models, without hard-to-derive
parameters, highlight potential aspects for real-life implementation.

###### Index Terms:

underfrequency load shedding, particle filters, phasor measurement units,
power system protection.

††publicationid: pubid: 978-1-6654-8537-1/22/$31.00 © 2022 IEEE

##  I Introduction

Currently, the majority of automatic underfrequency load shedding solutions
are built on decentralized architectures, and operate at the feeder level. UF
relays are installed at feeders considered non-critical by power system
operators [1]. Once the frequency seen by the relay drops below a threshold,
the relay issues a trip signal and the entire circuit is disconnected. Small
degrees of selectivity and coordination can be achieved by applying time
delays or multiple pick-up settings [2]. Despite these improvements, these
solutions, referred to in this work as traditional UF schemes, take a drastic
and unforgiving approach: if the relay sees the system’s frequency drop, the
customers on that circuit are disconnected. Traditional UF schemes are
reactive, as corrective actions occur only after a disturbance has been
observed. UF load shedding strategies of this type suffer from two commonly
observed issues: delayed response and overshedding [2].

While traditional UF schemes are by far the most widely utilized solution,
more intricate techniques leveraging the processing power of digital relays
have been suggested [1]. Semi-adaptive techniques, those considering time-
derivatives, and fully adaptive techniques have been proposed [3, 4, 5]. While
these offer improved performance compared to traditional schemes, the problems
of delayed response and overshedding are still present. Details regarding
frequency stability and mitigation methods can be found in [6, 7].

Clearly, corrective actions in the face of UF events could be orchestrated
better; however, a well choreographed response requires time. Precious seconds
can be gained by switching from the traditional reactive approach into a
proactive scheme [8].

Predictive schemes are not exactly new, one such approach was suggested in [9]
over two decades ago. The method presented in [9] has the same foundation as
many of the techniques proposed today: Measurements are collected via PMUs,
data is processed as a time series, and short term predictions are made. An
approach similar to [9] is presented in [10]. However, [10] focuses on steady
state UF mitigation. This technique makes decisions based on the predicted
steady state value of frequency after a disturbance is detected. More
recently, an elegant solution was presented in [11]. Measurements taken via
PMUs are processed by a prediction algorithm based on simple polynomial curve
fitting.

One significant limitation of the techniques found in literature is that they
cannot be supported by technology that is currently available or already in
service. For instance, these approaches normally rely on PMU measurements
collected at high sampling rates (in some cases over 100 Hz as in [9]),
meanwhile, the PMU sampling rate of modern digital relays is only 30 Hz[12].
This is the equipment demographic this method aims to exploit.

Another limitation found in contemporary UF frameworks is the use of non-
adaptive physics-based models. Any model carries an implicit modelling
uncertainty, and this problem is compounded when the state of the physical
system changes.

In light of the limitations of traditional and contemporary UF solutions, this
work presents a proactive framework for automatic UF control where model
uncertainty is continuously updated through particle filters [13]. This work
found that particle filters offer an exciting degree of accuracy and
robustness with relatively modest requirements in terms of hardware and
expertise.

A different approach to UF mitigation is introduced in [14]. Deviations in
frequency are compensated in real-time by actuating DERs. The results of [14]
are promising and are used as a benchmark in the third case study presented in
this work.

The rest of the paper is structured as follows: Section II provides a
mathematical introduction into the particle filter and the equations used in
system modeling. Section III showcases three case studies and discusses the
findings. Concluding remarks are given in Section IV.

##  II Theoretical Background

###  II-A Particle Filter

In the last decade filtering techniques, such as the Kalman filter (KF) have
gained attention from researchers in the area of power systems. This has been
driven by advancements in hardware, computing power, and the need to establish
a framework for dynamic state estimation [15]. While KF based techniques have
provided encouraging results, some limitations of the KF have been observed.
In particular, a common assumption that data points follow a Gaussian
distribution has been called into question in [16]. A lesser known filtering
technique, the particle filter (PF), solves this limitation by making multiple
predictions for each state being tracked, with each prediction having a
different probability. At each time step predictions and corrections,
equivalent to those in the KF are performed considering data from previous
time steps, an underlying system model, and the predictions made. The result
is a filter that is more flexible than the KF as illustrated in Figure 1.

Figure 1: UKF estimate (top) vs PF estimate (bottom).

When working with the PF, the optimal solution in Bayesian form is calculated
as a sum of samples. These samples are referred to as particles, and each one
is assigned a weight.

$p(x_{0:k}|z_{1:k})\approx\Sigma^{N_{s}}_{i=1}w^{i}_{k}\delta(x_{0:k}-x^{i}_{0:k})$
(1)

A set containing $N_{s}$ samples and weights can be expressed as
$\\{w^{i}_{k},x^{i}_{0:k}\\}^{N_{s}}_{i=1}$. The weight, or importance of each
sample $x^{i}_{0:k}$ is represented by $w^{i}_{k}$. The sum of all weights
$\Sigma^{N_{s}}_{i=1}w^{i}_{k}=1$. Samples thought to be of higher accuracy
are given a higher weight relative to samples of lower accuracy. Finally, the
Dirac delta function is represented by $\delta(\cdot)$. Weights are computed
via (2):

$w^{i}_{k}\>\alpha\>w^{i}_{k-1}\frac{p(z_{k}|x^{i}_{k})p(x^{i}_{k}|x^{i}_{k-1})}{q(x^{i}_{k}|x^{i}_{k-1},z_{k})}$
(2)

Without knowledge of a posterior, but given density $q$, the relationships in
(2) can be used to assign particle weights. This equation can be thought of as
a ratio between posterior and importance density for each particle. With
further manipulations it can be shown that the corresponding posterior can be
expressed as:

$p(x_{k}|z_{1:k})\approx\Sigma^{N_{s}}_{i=1}w^{i}_{k}\delta(x_{k}-x^{i}_{k})$
(3)

A key takeaway from (3) is that as the number of particles is increased, the
solution moves closer to the real values.

Resampling in the PF plays a similar role to the correction step in the KF.
Resampling attempts to correct imbalances in weight assignments that might
skew the overall performance of the filter.

Algorithm 1 Particle filter with resampling

0: $\\{x^{i}_{k-1},w^{i}_{k-1}\\}^{N_{s}}_{i=1},z_{k}$

0: $\\{x^{i}_{k},w^{i}_{k}\\}^{N_{s}}_{i=1}$ $w_{sum=0}$

1: for $i=1,...,N_{s}$ do

2: draw sample $x^{i}_{k}\approx q(x^{i}_{k}|x^{i}_{k-1},z_{k})$ assign weight
$w^{i}_{k}$ using (2)$w_{sum}=w_{sum}+w^{i}_{k}$

3: end for

4: for $i=1,...,N_{s}$ do

5: $w^{i}_{k}=w^{i}/w_{sum}$

6: end for

7: Resample $N_{s}$ particles with replacement

8: for $i=1,...,N_{s}$ do

9: $w^{i}_{k}=1/N_{s}$

10: end for

A complete derivation of the PF was omitted due space constraints but it can
be found in [13]. A comparison among several types of Bayesian filters can be
found in [17].

The PF can be considered a non-Gaussian extension of the KF. The trade-off for
this increase in flexibility is a modest increase in complexity. For this
reason, conservative processing delays are included in the case studies
presented in Section III.

###  II-B System Frequency

A combination of hard and soft thresholds were use in this work. The hard
thresholds are those seen in traditional UF schemes. Action is taken after
frequency drops below a set point. For the soft thresholds, actions are taken
based on the rate of change of frequency:

$R=\frac{f_{2}-f_{1}}{dt}$ (4)

were $R$ is the average rate of change in frequency, $f_{1}$ is the initial
frequency, while $f_{2}$ is the frequency at the end of the time window $dt$.
These soft thresholds have delays corresponding to the magnitude of the
difference between the two frequency measurements. These values are presented
as a lookup table in [18], with large values of $R$ (2.33 to 15 Hz/sec) having
a delay of only 3 cycles, while smaller values of $R$ (0.33 to 0.37 Hz/sec)
have a corresponding delay as large as 21 cycles.

###  II-C Prediction Problem

In their original formulation, the KF and the PF make predictions one time
step into the future ($k+1$). A simple way to extend the horizon of these
predictions is to feed artificial data points (ADPs) into the filter. The
number of ADPs required can be found by considering the sampling frequency and
the number of seconds into the future one wishes to predict:

$N_{ADP}=\frac{t_{p}}{f_{s}}$ (5)

Here $N_{ADP}$ is the number of artificial measurements required, $t_{p}$ is
the number of seconds into future, and $f_{s}$ is the sampling frequency of
the measuring device. In order to emulate the dynamic nature of the system,
the first and second deviates are calculated for each of the last ten time
steps, corresponding to the last ten data points before a prediction is made.
The derivatives are then averaged, before being used to systematically adjust
the last data point received to produce a vector of ADPs. This is a sequential
process based on the following equation:

$\begin{split}ADP_{i}&=ADP_{i-1}+t_{s}f^{\prime}+t^{2}_{s}f^{\prime\prime}\\\
\end{split}$ (6)

Where $ADP_{i-1}$, is the previous ADP. The average first derivative is
represented by $f^{\prime}$, while the average second derivative is
represented by $f^{\prime\prime}$. Finally, $t_{s}$ represents the time window
of the derivatives.

Algorithm 2 ADP Generation

0:

0: Initialisation: $ADP_{i-1}$ = Last measurement $f^{\prime}$ = Average first
derivative in last 10 measurements $f^{\prime\prime}$ = Average second
derivative in last 10 measurements $N_{ADP}$ = Number of ADPs required

1: for $i=1$ to $N_{ADP}$ do

2: $ADP_{i}=ADP_{i-1}+t_{s}f^{\prime}$
$f^{\prime}=f^{\prime}+t_{s}f^{\prime\prime}$

3: end for

The PF then processes the ADPs and iterates through predictions and correction
steps to generate future estimates as illustrated in Figure 2.

Figure 2: PF predictions with varying degrees of curvature.

The result is an algorithm that is able to capture the dynamics of the event
using a single curve, unlike the models presented in [10, 9, 11].

###  II-D Decision Making Problem

After a prediction is made, equations derived from the swing equation are used
to calculate the power imbalance [1]:

$L=\frac{R_{p}H(1-\frac{f^{2}_{p}}{f^{2}_{1}})}{p(f_{p}-f{1})}$ (7)

$L$ represents the load excess factor, $H$ is the inertia factor, $p$
represents the power factor, and $R_{p}$ is the predicted average rate of
change in frequency found with (4) using the current frequency measurement
$f_{1}$, and the predicted frequency measurement $f_{p}$. In this work,
frequency is predicted three seconds into the future; therefore, $f_{p}$ is
the predicted frequency value three seconds after $f_{1}$. A visual overview
of this framework is illustrated in Figure 3.

Figure 3: Conceptual overview of the solution.

##  III Case Studies

Three case studies are presented in this section. Each one utilizes the
proposed algorithm to mitigate frequency deviations in slightly different
ways. The simulations were performed in a reduced-order system model. This is
an accepted assumption in the study power system frequency stability [19]. All
related works mentioned in Section I also made this assumption. In order to
test the limits of this solution, a system with low inertia was chosen. In
addition to the low inertia of the system, parameters in the speed controller
of the generator were modified to decrease its performance. This produces a
system where highly dynamic frequency deviations can be observed. These
constraints would severely hinder the performance of the approaches suggested
in [10], and in [11]. The model used in the case studies and all its
parameters can be found in [20]. Gaussian noise with a variance of 0.025 was
added to the measurements.

###  III-A Case Study I: Single Stage Load Shed

In this scenario a single load shedding stage is used. Frequency deviations
start at 1.5 s, as illustrated in Figure 4. Thresholds are exceeded and a
prediction is made at 2.5 s.

Figure 4: Frequency deviation and prediction.

The frequency rate of change $R$ is calculated using (4). $f_{1}$ is the
frequency value estimated at 2.5 s, while $f_{2}$ is the predicted frequency
at 5.5 s. $R$ is then used in (7) to estimate the load excess factor $L$.

As depicted in Figure 5, the prediction in this case is not perfect; however,
the algorithm still manages to bring the system frequency back to a level
where generator governors can correct the deviation, as illustrated in Figure
6.

Figure 5: Discrepancy between predictions and actual values.

Corrective action is taken at 3.5 s, a full second after the thresholds were
exceeded.

Figure 6: Frequency at the end of the mitigation process.

The ability to take corrective actions early made up for a less than perfect
prediction where the load excess factor was underestimated. The accuracy of
the predictions and calculations can be improved by adjusting the thresholds
and switching to a multi-stage scheme as shown in Case Study II.

###  III-B Case Study II: Multi Stage Load Shed

In the second scenario, multiple load-shedding stages are used. This test case
highlights the ability of the solution to adapt and compensate for
inaccuracies made during the prediction step. In this case a prediction is
made at 3 s, as shown in Figure 7.

Figure 7: Initial Prediction.

The same process seen in Case Study I is applied, except that this time only
half of the calculated excess load is dropped, while a new prediction is made
one second after the first one.

Figure 8: Prediction and actual values.

The first load shedding action takes place at 4 s. Meanwhile, the prediction
is updated using data received after the initial prediction was made. This
correction is illustrated in Figure 9.

Figure 9: Corrected prediction.

With a new prediction, load-shedding action is once again executed at 5 s.
Both load-shedding stages can be observed in Figure 10.

Figure 10: Frequency at the end of the mitigation process.

As illustrated in Figure 10, the first load-shedding stage stops the decline
in frequency, while the second one sends it back to its normal range. A
subsequent prediction is made at 5 s but no corrective action is taken as the
algorithm predicts the frequency will be returning to normal levels.

###  III-C Case Study III: Real Time UF Compensation with Distributed Energy
Resources

The goal of this final test case is to highlight the flexibility of the
solution and use it to drive Distributed Energy Resources (DERs) in real-time.
This test was run under a similar set of assumptions as those made by [14].
Once again, in order to test the algorithm under demanding conditions, faster
frequency deviations than those seen in [14] were generated. Most importantly,
the total delay time involved in the processing of data and actuation of DERs
was increased to 500 ms; up from the 40 ms time delay used in [14]. That’s a
response time over ten times slower.

Frequency deviations start at 1 s, with a significant loss in generation at
1.5 s. As illustrated in Figure 11, DER actuation take place 0.5 s after
deviations in system frequency.

Figure 11: Real-Time UF mitigation via DERs.

As before, measurements are made via PMUs at 30 Hz. The power mismatch is
calculated continuously in this test case per (4)-(7) in the form of a
controller. As shown in Figure 11 above, the frequency and power mismatch
follow virtually the same trend but in opposite directions. When frequency
deviates from the 60 Hz reference, a corresponding current output is seen from
the DERs based on the power mismatch calculated. Despite a 0.5 second delay
before DER actuation, the system successfully mitigates the frequency
deviations.

##  IV Conclusion

A proactive and automatic underfrequency load shedding solution was presented
in this work. The solution leverages the particle filter along with existing
technology to deliver a real-time and predictive UF mitigation scheme. Several
limitations of existing techniques were addressed and improved upon,
particularly in regards to real-life implementation. Test results indicate
that regardless of the decision making strategy in place, load shedding or
DERs actuation, frequency stability is considerably improved by this technique
compared to traditional UF schemes. Results from the case studies also
highlight the robustness and the flexibility of the particle filter.

## References

  * [1] S. H. Horowitz and A. G. Phadke, Power System Relaying, 4th, Wiley, West Sussex, United Kingdom, 2014.
  * [2] L. Sigrist, L. Rouco, and F. M. Echavarren, “A review of the state of the art of UFLS schemes for isolated power systems,” in International Journal of Electrical Power & Energy Systems, vol 99, pp. 525–539, 2018, doi: 10.1016/j.ijepes.2018.01.052.
  * [3] P. M. Anderson and M. Mirheydar, ”An adaptive method for setting underfrequency load shedding relays,” in IEEE Transactions on Power Systems, vol. 7, no. 2, pp. 647-655, May 1992, doi: 10.1109/59.141770.
  * [4] M. S. Pasand and H. Seyedi, ”New Centralized Adaptive Under Frequency Load Shedding Algorithms,” 2007 Large Engineering Systems Conference on Power Engineering, 2007, pp. 44-48, doi: 10.1109/LESCPE.2007.4437350.
  * [5] S. Abdelwahid, A. Babiker, A. Eltom and G. Kobet, ”Hardware Implementation of an Automatic Adaptive Centralized Underfrequency Load Shedding Scheme,” in IEEE Transactions on Power Delivery, vol. 29, no. 6, pp. 2664-2673, Dec. 2014, doi: 10.1109/TPWRD.2014.2331495.
  * [6] P. Kundur, Power System Stability and Control, 2nd; McGraw-Hill, New York, NY, USA, 1994.
  * [7] P. Kundur et al., ”Definition and classification of power system stability IEEE/CIGRE joint task force on stability terms and definitions,” in IEEE Transactions on Power Systems, vol. 19, no. 3, pp. 1387-1401, Aug. 2004, doi: 10.1109/TPWRS.2004.825981.
  * [8] G. Paramo, A. Bretas, and S. Meyn, “Research Trends and Applications of PMUs,” Energies, vol. 15, no. 15, p. 5329, Jul. 2022, doi: 10.3390/en15155329.
  * [9] N. G. Bretas and A. G. Phadke, ”Real time instability prediction through adaptive time series coefficients,” IEEE Power Engineering Society. 1999 Winter Meeting (Cat. No.99CH36233), 1999, pp. 731-736 vol.1, doi: 10.1109/PESW.1999.747547.
  * [10] M. Larsson and C. Rehtanz, ”Predictive frequency stability control based on wide-area phasor measurements,” IEEE Power Engineering Society Summer Meeting,, 2002, pp. 233-238 vol.1, doi: 10.1109/PESS.2002.1043222.
  * [11] U. Rudez and R. Mihalic, ”WAMS-Based Underfrequency Load Shedding With Short-Term Frequency Prediction,” in IEEE Transactions on Power Delivery, vol. 31, no. 4, pp. 1912-1920, Aug. 2016, doi: 10.1109/TPWRD.2015.2503734.
  * [12] SEL-421 https://selinc.com/products/421 (accessed on 17 Feb 2022).
  * [13] J. Elfring, E. Torta, and R. van de Molengraft, “Particle Filters: A Hands-On Tutorial,” Sensors, vol. 21, no. 2, p. 438, Jan. 2021, doi: 10.3390/s21020438.
  * [14] S. Pulendran and J. E. Tate, ”Energy Storage System Control for Prevention of Transient Under-Frequency Load Shedding,” in IEEE Transactions on Smart Grid, vol. 8, no. 2, pp. 927-936, March 2017, doi: 10.1109/TSG.2015.2476963.
  * [15] J. Zhao et al., ”Power System Dynamic State Estimation: Motivations, Definitions, Methodologies, and Future Work,” in IEEE Transactions on Power Systems, vol. 34, no. 4, pp. 3188-3198, July 2019, doi: 10.1109/TPWRS.2019.2894769.
  * [16] S. Wang, J. Zhao, Z. Huang and R. Diao, ”Assessing Gaussian Assumption of PMU Measurement Error Using Field Data,” in IEEE Transactions on Power Delivery, vol. 33, no. 6, pp. 3233-3236, Dec. 2018, doi: 10.1109/TPWRD.2017.2762927.
  * [17] M. S. Arulampalam, S. Maskell, N. Gordon and T. Clapp, ”A tutorial on particle filters for online nonlinear/non-Gaussian Bayesian tracking,” in IEEE Transactions on Signal Processing, vol. 50, no. 2, pp. 174-188, Feb. 2002, doi: 10.1109/78.978374.
  * [18] SEL-751 https://selinc.com/products/751 (accessed on 1 July 2022).
  * [19] E. Farantatos, R. Huang, G. J. Cokkinides and A. P. Meliopoulos, ”A Predictive Generator Out-of-Step Protection and Transient Stability Monitoring Scheme Enabled by a Distributed Dynamic State Estimator,” in IEEE Transactions on Power Delivery, vol. 31, no. 4, pp. 1826-1835, Aug. 2016, doi: 10.1109/TPWRD.2015.2512268.
  * [20] G. Sybille, ”Simplified Synchronous Machine - Speed Regulation”, MathWorks, https://www.mathworks.com/help/sps/ug/simplified-synchronous-machine-speed-regulation.html (accessed on 1 Sep 2022).