Instantaneous Frequency Identification in Microgrids Through Adaptive Data Analysis

(1)

Erlend WestadInstantaneous Frequency Identification in Microgrids Through Adaptive Data Analysis NTNU Norwegian University of Science and Technology Faculty of Information Technology and Electrical Engineering Department of Electric Power Engineering

Master ’s thesis

Erlend Westad

Instantaneous Frequency Identification in Microgrids Through Adaptive Data Analysis

Master’s thesis in Energy and Environmental Engineering Supervisor: Olav Bjarte Fosso and Marta Molinas

June 2020

(2)

(3)

Erlend Westad

Instantaneous Frequency

Identification in Microgrids Through Adaptive Data Analysis

Master’s thesis in Energy and Environmental Engineering Supervisor: Olav Bjarte Fosso and Marta Molinas

June 2020

Norwegian University of Science and Technology

Faculty of Information Technology and Electrical Engineering

Department of Electric Power Engineering

(4)

(5)

Abstract

With increased integration of power electronic equipment and distributed energy generation, the modern power system has become more prone to harmonic pollution. As isolated systems categorized by low inertia such as microgrids are more common, the presence of non-linear distortion is becoming an increasing problem. The most common methods used for surveillance in the power system are as of now based on average value calculations. They are thus not suited for the rising non-linearity caused by the harmonics in the evolving system. For monitoring and control equipment to keep operating with adequate performance, there is a need for alternative surveillance methods capable of instantaneous frequency and amplitude detection.

This thesis aims to explore the use of adaptive data analysis as an alternative surveillance method for harmonic detection in the power system. Empirical Mode Decomposi- tion (EMD) and its real-time extension, Online EMD, have been used in conjunction with Hilbert-Transform (HT) and Fast Fourier Transform (FFT) for extensive instantaneous frequency and amplitude identification. First, a comparison of the decomposition quality of EMD and Online EMD is made using synthetic signals, looking at different stoppage criteria and their reaction to disturbances. The methods have then been applied for harmonic detection in an actual voltage and current measurement taken from a wind turbine generator. The measurements have been decomposed by EMD, utilizing FFT and a recreated boundary condition map to construct appropriate masking signals to avoid mode mixing.

The use of masking signals to avoid mode-mixing has then been adapted and used in the Online EMD for the same measurements. As adequately pure IMFs are extracted, they are removed from the original signal, and the new resulting signal is treated as new input for further decomposition. As all relevant modes are extracted, they are analyzed by HT and FFT, and the results from the EMD and Online EMD are compared.

Essentially, the methods proved to be powerful tools for harmonic detection but requires techniques to handle mode mixing when applied to complex signals. The use of masking signals is shown to be a highly effective mode mixing separation technique and is seen to also be effective when adapted for use with the Online EMD.

(6)

(7)

Sammendrag

Den økende integreringen av kraftelektroniske komponenter har sammen med mer dis- tribuert generering av energi ført til at kraftsystemet er mer usatt for harmonisk foruren- sning. Isolerte systemer typisk kategorisert av lavt treghetsmoment er spesielt usatt, og har gjort ikke-lineære harmonisk støy til et økende problem. Dagens monitoreringsutstyr benytter seg i stor grad av metoder som baserer seg p˚a gjennomsnittlige verdier og utreg- ninger, og er dermed ikke egnet for ˚a benyttes ved større grad av ikke-linearitet. Det trengs dermed en alternativ metode med mulighet for ˚a beregne momentant frekvens og amplitude slik at kontrollutstyr skal kunne fungere med nødvendig nøyaktighet og presisjon.

Denne oppgaven utforsker bruken av adaptiv data analyse som en alternativ metode for detektering av harmoniske komponenter i kraftsystemet. Empirical Mode Decomposition (EMD) og sanntids versjonen Online EMD har blitt brukt sammen med Hilbert-Transform (HT) og Fast Fourier Transform (FFT) for ˚a se nærmere p˚a muligheten for momentan frekvens og amplitude beregning. En sammenlikning av dekomposisjon kvaliteten til EMD og Online EMD er gjort for forskjellige stop-kriterier og forstyrrelser. Deretter har metodene blitt brukt for ˚a indentifisere det harmoniske innholdet i en virkelig strøm og spenningsm˚aling henter fra en vind turbin. EMD ble brukt til ˚a dekomponere m˚alingene, og FFT ble brukt sammen med et grensetilstands-kart for ˚a konstruere maskerings signal for ˚a forhindre frekvens miksing. Bruken av maskerings signal ble s˚a modifisert til bruk sammen med Online EMD, for de samme m˚alingene. Tilfredsstillende IMFer ble deretter fjernet fra det originale signalet og det resulterende signalet ble analysert p˚a nytt. Etter at alle relevante komponenter har blitt fjernet fra signalet blir de videre analysert ved hjelp av HT og FFT, og resultatene fra EMD og Online EMD blir sammenliknet.

Metodene brukt i denne oppgaven viste seg ˚a være gode verktøy for ˚a finne det harmoniske innholdet i et signal, men krever teknikker som kan behandle mode-miksing ved mer kom- plekse signaler. Maskerings-signal var veldig effektiv seperasjons metode, og den viste seg og ogs˚a fungere godt sammen med Online EMD.

(8)

(9)

Preface

This is the master’s thesis to conclude the Master of Science degree in Energy and Envi- ronmental Engineering at the Norwegian University of Science and Technology. The work was carried out at the Department of Electric Power Engineering during spring 2020.

I want to express my sincere gratitude to my supervisor Professor Olav Bjarte Fosso, for excellent guidance and for providing the opportunity to work on such an interesting topic.

I also want to thank my co-supervisor, Professor Marta Molinas. Her insightful lectures sparked my interest for the world of signal analysis.

Trondheim, June 4, 2020 Erlend Westad

(10)

(11)

List of Figures

2.1 Waveforms of linear and non linear loads . . . 7

2.2 FFT applied to a linear and stationary signal . . . 11

2.3 FFT applied to a signal with varying amplitude . . . 11

2.4 FFT applied to a signal with increasing frequency . . . 11

2.5 FFT applied to a signal with increasing frequency and varying amplitude 12 2.6 FFT applied two superimposed signals . . . 12

2.7 A signal with varying amplitude and its envelope, calculated using the magnitude of the Hilbert Transform . . . 13

2.8 The complex phasor representation, instantaneous frequency and instantaneous amplitude resulting from the HT applied to a linear and stationary signal . . . 15

2.9 The complex phasor representation, instantaneous frequency and instantaneous amplitude resulting from the HT applied to a signal with varying amplitude . . . 15

2.10 The complex phasor representation, instantaneous frequency and instantaneous amplitude resulting from the HT applied to a signal with increasing frequency . . . 15

2.11 The complex phasor representation, instantaneous frequency and instantaneous amplitude resulting from the HT applied to a signal with increasing frequency and varying amplitude . . . 16

2.12 The complex phasor representation, instantaneous frequency and instantaneous amplitude resulting from the HT applied to a signal with offset . . . 16

2.13 The complex phasor representation, instantaneous frequency and instantaneous amplitude resulting from the HT applied to two superimposed signals . . . 16

2.14 Flow diagram of the EMD algorithm . . . 18

2.15 Example of the mean of a signal calculated by use of local envelopes based on local extrema. . . 19

2.16 EMD applied to a signal composite of three different frequencies. . . 19

2.17 Example of improper and proper boundary conditions. . . 21

(14)

they get closer(right), an interpretation in terms of a single tone modulated in amplitude is favored . . . 22 2.19 The mode-mixing criteria visualized . . . 23 2.20 Boundary condition map for mode-mixing in EMD. Red area represent

mode-mixing, blue area represent clean separation of two signals. . . 24 2.21 EMD decomposition of synthetic signalx(t) =cos(2π10t) +cos(2π5t).

The series has been decomposed into four equal segments each with a separate EMD indicated by the color scheme. After adding the resulting IMFs, there are clear boundary errors at t = 1.25, 2.5 and 3.75. . . 25 2.22 Schema of the SEMD algorithm . . . 27 2.23 Error in resulting IMFs given signal as in Figure 2.21 . . . 27 2.24 Overview of Online EMD sliding window and stitching procedure. In this

case, the l = 10 extrema,Mⁱgiven by the blue signal, represents the fastest oscillation and is extracted using the classical EMD. The window function φ(t)is plotted in green, and the resulting weighted IMF,M¯, stitched to the previous uncovered IMFs is red. Note that herel⁰=l. . . 29 2.25 Execution time of Online EMD (l = 20) and classical EMD (both using

Rilling stopping criterion) with a white noise signal and a sinusoid with a trend . . . 31 3.1 Signal consisting of three different amplitudes and modes of oscillations. . 34 3.2 EMD with stoppage criteria SD = 0.01. Overall and rms error is calculated

for each extracted IMF. . . 35 3.3 EMD with stoppage criteriaα= 0.001 θ1= 0.001 θ2= 0.01. Overall

and rms error is calculated for each extracted IMF. . . 35 3.4 mal.png . . . 36 3.5 EMD with stoppage criteria FIX = 10. Overall and rms error is calculated

for each extracted IMF. . . 36 3.6 Average root-mean-square error for all IMFs as the widow size increases

for the Online EMD. . . 36 3.7 Signal consisting of three different amplitudes and modes of oscillations

with a added Gaussian white noise component, matched to scale by the MATLAB awgn(’measured’) command. . . 37 3.8 The first 8 IMFs from using the stoppage criteriaSD= 0.01. . . 38 3.9 FFT of the first 8 IMFs using the stoppage criteriaSD= 0.01. . . 38 3.10 EMD with stoppage criteriaSD = 0.01. Overall and rms error is calcu-

lated for each extracted IMF. . . 39 3.11 The first 8 IMFs from using the stoppage criteria α = 0.001 θ₁ =

0.001 θ₂= 0.01. . . 40 3.12 FFT of the first 8 IMFs using the stoppage criteriaα = 0.001 θ₁ =

0.001 θ2= 0.01. . . 40 3.13 EMD with stoppage criteriaα= 0.001 θ1= 0.001 θ2= 0.01. Overall

and rms error is calculated for each extracted IMF. . . 41 3.14 The first 7 IMFs from using the stoppage criteria FIX = 10. . . 42

(15)

3.15 FFT of the first 7 IMFs using the stoppage criteria FIX = 10. . . 42

3.16 EMD with stoppage criteria FIX = 10. Overall and rms error is calculated for each extracted IMF . . . 43

3.17 Online EMD with FIX = 10 as stoppage criteria. The red parts represents the section of the IMFs not yet completed due to lacking number of extrema. 44 3.18 Average root-mean-square error for all IMFs as the widow size increases for the Online EMD. . . 44

4.1 Current measurement from a WTG at Hundhammarhjellet, taken the 8^th of October 2014 from 23:34:02 to 23:34:06 with a sampling frequency of 25600 Hz. . . 48

4.2 Voltage measurement from a WTG at Hundhammarhjellet, taken the 8^th of October 2014 from 23:34:02 to 23:34:06 with a sampling frequency of 25600 Hz. . . 48

4.3 Closer look at the current and voltage measurement. . . 48

4.4 Current measurementf(t). . . 49

4.5 FFT of the dominant frequency, and a close-up of the other frequency components present inf(t). . . 49

4.6 Resulting IMFs from EMD applied tof(t). . . 50

4.7 Closer look at the FFT of IMF 1. . . 50

4.8 Resulting IMFs from EMD applied tof(t)after adding masking signal f_m1(t). . . 51

4.9 FFT of the IMFs from EMD applied tof(t)after adding masking signal fm1(t). . . 51

4.10 Closer look at the FFT of IMF 1 after adding masking signalfm1(t). . . . 52

4.11 Resulting signalf₁(t), after removing 50 Hz component. . . 52

4.12 Resulting IMFs from EMD applied tof1(t). . . 53

4.13 FFT of the IMFs from EMD applied tof₁(t). . . 53

4.14 Closer look at FFT of IMF 1 of signalf1(t). . . 54

4.15 Resulting signalf2(t), after removing 50 Hz and 950 Hz component. . . . 54

4.16 Resulting IMFs from EMD appliedf₂(t). . . 55

4.17 FFT of the IMFs from EMD applied tof2(t). . . 55

4.18 Resulting IMFs from EMD appliedf₂(t)after adding masking signalf_m2(t). 56 4.19 FFT of the IMFs from EMD applied tof2(t)after adding masking signal fm2(t). . . 56

4.20 Voltage measurementg(t) . . . 58

4.21 FFT of the dominant frequency, and a close-up of the other frequency components present ing(t). . . 58

4.22 Resulting IMFs from EMD applied tog(t). . . 59

4.23 FFT of the IMFs from EMD applied tog(t). . . 59

4.24 FFT of IMF 1 after applying EMD . . . 60

4.25 Resulting IMFs from EMD applied tog(t)after adding masking signal gm1(t). . . 61

4.26 FFT of the IMFs from EMD applied tog(t)after adding masking signal gm1(t). . . 61

(16)

4.28 Resulting signalg1(t), after removing 50 Hz component. . . 62

4.29 FFT ofg₁(t)with and without boundary condition errors . . . 63

4.30 Resulting IMFs from EMD applied tog1(t). . . 64

4.31 FFT of the IMFs from EMD applied tog₁(t). . . 64

4.32 Resulting IMFs from EMD applied tog1(t)after adding masking signal gm2(t). . . 65

4.33 FFT of the IMFs from EMD applied tog₁(t)after adding masking signal g_m2(t). . . 65

4.34 Resulting signalg2(t), after removing 50 and 950 Hz component. . . 66

4.35 Resulting IMFs from EMD applied tog₂(t). . . 67

4.36 FFT of the IMFs from EMD applied tog2(t). . . 67

4.37 Resulting IMFs from EMD applied tog₂(t)after adding masking signal g_m3(t). . . 68

4.38 FFT of the IMFs from EMD applied tog2(t)after adding masking signal gm3(t). . . 68

4.39 Resulting signalg₃(t). . . 69

4.40 Online EMD applied toL1.From left to right, the total data points are 662, 2650, and 10600. The black part of the signal is related to the first 2650 data points. The blue part is related to the parts 2650 to 5300, and teal is related to the parts from 5300 to 10600. Red indicates the parts of the IMFs still incomplete . . . 70

4.41 Online EMD applied tof(t). In red are parts of the IMFs still incomplete. 71 4.42 FFT of the resulting IMFs off(t)given by Online EMD. . . 71

4.43 IMFs after Masking signalf_m1(t)is applied to signalf(t)for the Online EMD . . . 72

4.44 FFT of the IMFs after masking signalfm1(t)is applied to signalf(t)for the Online EMD . . . 72

4.45 Resulting IMFs after applying Online EMD tof1(t). . . 73

4.46 FFT of resulting IMFs after applying Online EMD tof1(t) . . . 74

4.47 Online EMD applied toL12. From left to right, the total data points number 5120, 10240, and 20480. The black part of the signal is related to the first 5120 data points. The blue part is related to the parts 5120 to 10240, and teal is related to the parts from 10240 to 20480. Red indicates the parts of the IMFs still incomplete. . . 74

4.48 The resulting IMFs after applying Online EMD toL12 . . . 75

4.49 FFTs of the resulting IMFs after applying Online EMD tog(t) . . . 75

4.50 Resulting IMFs after applying masking signalg_m1(t)to the Online EMD 76 4.51 FFT of resulting IMFs after applying masking signalgm1(t)to the Online EMD . . . 76

4.52 Resulting IMFs after applying Online EMD tog₁(t). . . 77

4.53 FFT of resulting IMFs after applying Online EMD tog1(t) . . . 77

4.54 Resulting IMFs after adding masking signalgm2tog1(t) . . . 78

4.55 FFT of resulting IMFs after adding masking signalgm2tog1(t) . . . 78

(17)

4.56 Final resulting IMFs for current measurementL1using classical EMD. . 79

4.57 FFT of final resulting IMFs for current measurementL1 using classical EMD. . . 80

4.58 IF of final resulting IMFs for current measurementL1using classical EMD. 80 4.59 IA of final resulting IMFs for current measurementL1using classical EMD. 81 4.60 Final resulting IMFs for current measurementL1using Online EMD. . . 82

4.61 FFT of final resulting IMFs for current measurementL1 using classical EMD. . . 82

4.62 IF of final resulting IMFs for current measurementL1using Online EMD. 83 4.63 IA of final resulting IMFs for current measurementL1using Online EMD. 83 4.64 Final resulting IMFs for voltage measurementL12using classical EMD. . 84

4.65 Close-up of the resulting IMFs for voltage measurementL12using classical EMD. . . 84

4.66 FFT of final resulting IMFs for voltage measurementL12using classical EMD. . . 85

4.67 IF of final resulting IMFs for voltage measurementL12using classical EMD. . . 85

4.68 IA of final resulting IMFs for voltage measurementL12 using classical EMD. . . 86

4.69 Final resulting IMFs for voltage measurementL12using Online EMD. . . 86

4.70 Closer look at the final resulting IMFs for voltage measurementL12using Online EMD. . . 87

4.71 Closer look at the final resulting IMFs for voltage measurementL12using Online EMD. . . 87

4.72 IF of final resulting IMFs for voltage measurementL12using Online EMD. 88 4.73 IA of final resulting IMFs for voltage measurementL12using Online EMD. 88 4.74 The frequency component of the HT for IMF 1 of signalL1compared to the frequency component of the HT of a signal given bycos(2π950t) + 0.07cos(2π1050t). . . 89

4.75 Closer look a the amplitude component of the HT for IMF 4 of signalL1 90 A.1 IF of the IMFs from the EMD applied tof(t). . . 101

A.2 IA of the IMFs from the EMD applied tof(t). . . 102

A.3 IF of IMFs after adding masking signalf_m1(t)tof(t). . . 102

A.4 IA of IMFs after adding masking signalf_m1(t)tof(t). . . 103

A.5 IF of the IMFs from the EMD applied tof₁(t). . . 103

A.6 IA of the IMFs from the EMD applied tof₁(t). . . 104

A.7 IF of the IMFs from the EMD applied tof₂(t). . . 104

A.8 IA of the IMFs from the EMD applied tof₂(t). . . 105

A.9 IF of IMFs after adding masking signalf_m2(t)tof₂(t). . . 105

A.10 IA of IMFs after adding masking signalf_m2(t)tof₂(t). . . 106

A.11 Closer look at the IF of resulting IMFs off(t). . . 106

A.12 Closer look at the IA of resulting IMFs off(t). . . 107

A.13 Recreation of original signalf(t)from the resulting IMFs. . . 107

A.14 IF of the IMFs from the EMD applied tog(t). . . 108

A.15 IA of the IMFs from the EMD applied tog(t). . . 108

(18)

A.18 IF of the IMFs from the EMD applied tog₁(t). . . 110

A.19 IA of the IMFs from the EMD applied tog₁(t). . . 110

A.20 IF of IMFs after adding masking signalgm2(t)tog1(t). . . 111

A.21 IA of IMFs after adding masking signalgm2(t)tog1(t). . . 111

A.22 IF of the IMFs from the EMD applied tog2(t). . . 112

A.23 IA of the IMFs from the EMD applied tog2(t). . . 112

A.24 IF of IMFs after adding masking signalgm3(t)tog2(t). . . 113

A.25 IA of IMFs after adding masking signalgm3(t)tog2(t). . . 113

A.26 Closer look at the IF of resulting IMFs ofg(t). . . 114

A.27 Closer look at the IA of resulting IMFs ofg(t). . . 114

A.28 Recreation of original signalg(t)from the resulting IMFs. . . 115

(19)

Abbreviations

DER Distributed Energy Resources DFT Discrete Fourier Transform EMD Empirical Mode Decomposition FFT Fast Fourier Transform

HHT Hilbert-Huang Transform HT Hilbert Transform IA Instantaneous Amplitude IF Instantaneous Frequency IMF Intrinsic Mode Function PCC Point of Common Coupling RMS Root-Mean-Square

STFT Short-Time Fourier Transform THD Total Harmonic Distortion WT Wavelet Transform WTG Wind Turbine Generator

(20)

(21)

Chapter 1 Introduction

1.1 Background and Motivation

This thesis is a continuation of a specialization project done prior to this work, thus the background and motivation is here presented as an adapted and revised version from [1].

For a long time, the power system has been characterized by excellent power and frequency quality, ensured by the high inertia rotating masses of the generators in the system.

As such, the frequency was not likely to deviate far from the fundamental frequency, usually keeping a steady 50 Hz or 60 Hz. However, as the electrical grid is evolving, new components, such as power electronics converters, rectifiers, and adjustable speed drives, now interact at a larger scale. Combined with the increase of non-linear loads, the electrical grid is more prone to harmonic pollution as harmonic injection and non-linear distortions are becoming common.

Power systems with no connection to a stiff grid, typically categorized by low inertia, are especially susceptible to harmonic pollution. Isolated microgrids such as wind farms, marine vessels, and solar microgrids are only some examples of such systems as large frequency deviations caused by harmonics can cause significant damage. Furthermore, as a result of more advanced inverter control, conditions with varying load demand may lead to time-varying oscillations in such power systems [2]. To satisfy the operational conditions of the low inertia microgrids, sophisticated control and monitoring methods are needed.

Monitoring methods linked to frequency-related issues have been relying on Fourier based analysis, and the concept ofinstantaneous frequencyhas yet to be used on a larger scale for monitoring or control purposes in the electric power system [3].

The standard methods used for data acquisition in the power system today assumes a constant fundamental frequency, relying on average value calculations and only consid- ered constant harmonics that are integer multiples of the fundamental frequency [4, 5, 6].

Looking at the IEC Standard 61000-4-7, it is apparent that the currently used methods are

(22)

inaccurate when dealing with environments of time-varying angular frequency [7]. For control and monitoring equipment to keep adequate control, measurement, and estimation techniques relying on average values should transition to rely on instantaneous frequency and amplitude [2][8]. However, using mathematical formulations to describe physical behavior in large scale systems can be challenging. Controllers relying on mathematical models and approximations may lead to complex control schemes whose operations might be challenging to understand, and implementation might be costly due to their complex nature. As opposed to this high-fidelity modeling, this thesis explores the use of adaptive data analysis as an alternative method for frequency and harmonic monitoring in the power system. Using methods capable of handling both non-linear and non-stationary signals, available measurements are used to characterize the grid under investigation. Relying on instantaneous values, the control systems used in microgrids should be able to perform faster and more accurate actions, resulting in a more secure and stable power grid.

1.2 Objective and Scope of Work

This thesis’s objective is to explore adaptive data analysis as an alternative to present methods used for monitoring and control in the power system. This is done by assessing the potential of EMD and online EMD when used in conjunction with HT and FFT. Thus, the main methods are thorough review of existing literature and further implementation of the methods in MATLAB. Further investigation is done by constructing synthetic signals to compare decomposition quality, and by applying the detection methods to an actual voltage and current measurement to identify the harmonic content. Techniques for mode-mixing separation in the EMD are also explored, mainly focusing on masking signals constructed by the help of frequency identification through FFT and a boundary condition map. This technique is then adapted and investigated when used with the Online EMD.

1.3 Structure of Thesis

• Chapter 1: Introduction- Here the background and motivation for this thesis is presented. Also included are the Objective and Scope of work, and the structure of the thesis.

• Chapter 2: State-of-the-art in Methods for Frequency Identifications in Mi- crogrids- This chapter presents a more thorough theoretical background for the relevant techniques used in this thesis. A brief introduction to the concepts of microgrids and harmonics is also given.

• Chapter 3: Comparison of Classical and Online EMD- In this chapter a comparison of the decomposition quality of classical and Online EMD is made by various synthetic signals.

• Chapter 4: Frequency and Amplitude Identification Through HHT and Online EMD- In this chapter a real voltage and current measurement is analysed by EMD, Online EMD, HT and FFT.

(23)

1.3 Structure of Thesis

• Chapter 5: Conclusion and Future Work- Conclusion of the thesis and recom- mendations for future work.

• Appendix A:- Contains more graphs from the step-by-step decomposition of sig- nalsL1andL12, and a recreation of the original signals from the resulting IMFs. A closer look at the IF and IA plot for the resulting IMFs is also included.

(24)

(25)

Chapter 2 State-of-the-art in Methods for Frequency Identification in

Microgrids

The state of art were reviewed, and an identification of some of the relevant background material were carried out in the project preceding this thesis [1]. The relevant sections have thus been revised and adapted to be presented in this masters thesis.

2.1 Concept of Microgrids and Harmonics

2.1.1 Microgrids

The definition of a microgrid given by the US department of energy is [9]: ”a group of interconnected loads and distributed energy resources with clearly defined electrical boundaries that act as a single controllable entity for the grid and can connect and dis- connect from the grid to enable it to operate in both grid-connected or island modes”.

Being a combination of storage devices, distributed energy sources, and flexible loads, the microgrid is an entirely localized energy system capable of operating both connected and disconnected to the traditional macro grid. Disconnecting from the grid is often called

“islanding,” while the interconnection point is referred to as Point of Common Coupling (PCC) [10].

An essential aspect of the microgrid is local generation capability, thus integrating local renewable energy sources such as solar, wind, small hydro, geothermal, waste-to-energy, and combined heat and power is done to ensure reliable generation of power when not connected to the main grid. Less distance between generation and load leads to lower transmission losses, resulting in a more efficient power system. By integrating the local generation, the microgrid can operate as a power source for the main grid in times of out-

(26)

age or failure. This not only ensures a more secure and reliable main grid but results in more reliable and affordable energy for both urban and rural communities.

Contrary to the high inertia of the traditional power system, lower inertia is a significant concern for the microgrid as frequency, and speed deviation is inversely proportional to the inertia of the system. The high inertia rotating masses, i.e., synchronous generators in the traditional power grid, serve as storage for additional energy in the system during disturbances, thus ensuring constant voltage and frequency levels [11]. Furthermore, an increased load on the system causes a mismatch between the mechanical and electrical power, thus slowing the rotor using the kinetic energy to compensate for the difference in supply and demand. However, the low inertia of the microgrid makes such control schemes ineffective when dealing with system transients. The microgrid thus needs sophisticated frequency control and surveillance, as well as efficient energy storage to deal with transients when operating independently from the main grid.

2.1.2 Harmonics

This section provides a short introduction to the concept of harmonics. The information presented here is mainly from [12], where more thorough information on harmonics and its impact on the power system can be found. Here harmonics refers to the components of a waveform that is an integer multiple of that waveform’s fundamental frequency [12].

When mentioned in the power system, harmonics usually refer to the harmonic content found in either the current or voltage of the system. Thus, for a power system with a fundamental frequency of 50 Hz, the third harmonic implies the presence of a 150 Hz component and the fourth harmonic a 200 Hz component. Thus, harmonic components are referred to by order and is given by:

fh=h·fn (2.1)

Wheref_n is the fundamental frequency andhis the harmonic order. Such that the fifth order harmonicf₅is given by

f5= 5·50Hz= 250Hz (2.2)

Non-integer multiples of the fundamental frequency may also appear, though less common, this phenomenon is referred to as inter-harmonics. Harmonic voltages and currents are sources of many problems within the power system. These problems include over- heating of equipment and cabling, reduced energy efficiency, and reduced functionality due to loss of electromagnetic compatibility [12]. Furthermore, the harmonic current from installations can flow back into the network propagating as voltage harmonics, causing the supply waveform to be distorted, leading to increased network losses and lower reliabil- ity for several types of equipment [12]. The presence of harmonic current in the electric supply system is not a new phenomenon. The phenomenon was initially produced by the mercury arc rectifiers when converting AC to DC for railway electrification, by industrial DC variable speed drives, and by direct half-wave rectification used in radio and television sets [12]. However, the modern power system relying more on power electronic components and other passive equipment interacting at a larger scale has caused the number of

(27)

2.1 Concept of Microgrids and Harmonics

units causing harmonics to increase[12].

Total Harmonic Distortion (THD) is used to measure the degree of harmonics present in the power system. THD is defined as the ratio between the harmonics and the fundamental frequency present in the signal. Let the current with harmonic presence be given as:

is(t) =is1(t) +X

h6=1

ish(t) (2.3)

Whereis1is the fundamental component andish(t)the harmonic component. The THD of the current in 2.3 referreing to RMS value is given as [13]:

%T HD= 100·

pi²_s−i²_s1 is1

= 100· s

X

h6=1

(ish

is1

)² (2.4)

In an ideal power system, the current and voltage waveforms are pure sinusoids with no distortions. Realistically the system will contain some harmonics. When a non-linear re- lation between the applied voltage and current flowing in the load occurs, non-sinusoidal currents (and voltages) are generated. This should not be confused with a phase shift between the waveforms, as usually seen in simple circuits containing capacitance or in- ductance, as the current flowing will still be proportional to the voltage. However, given a circuit consisting of a full-wave rectifier and a capacitor, the supply voltage must exceed that stored on the reservoir’s capacitor for the current to flow. Thus, the current only flows when the voltage is close to the peak of the voltage sine wave [12]. This produces a cyclic waveform as seen in Figure 2.1(b).

(a)Current waveform in a linear load (b)Current waveform in a non linear load Figure 2.1:Waveforms of linear and non linear loads [12]

Though this waveform is a cropped version of the waveform given a linear load as seen in Figure 2.1(a), this current waveform can be deconstructed into several sine waves that together result in what can be seen in Figure 2.1(b). All the additional sine waves contributing to the flat shape will be the resulting harmonics in the wave. This makes it apparent

(28)

that the distorted current waveform, as seen in Figure 2.1(b) is the sum of the fundamental harmonic and a percentage of some higher-order harmonics.

(29)

2.2 Fourier Transform

The Fourier transform (FT), developed by Jean B.J Fourier in 1807, is regarded as a powerful tool in spectrum analysis. However, the method did not gain its current popularity before the development of the Fast Fourier Transform (FFT) in 1965 [14]. Based on the Discrete Fourier Transform (DTF), the FFT is a method for efficiently computing the finite FT ofNcomplex data points inN log2N operations [15]. Defined in Equation (2.5), the FT results in the Fourier expansion seen in Equation (2.6), a trigonometric function that is complete, convergent, orthogonal and unique.

F{g(t)}=G(f) = Z ∞

−∞

g(t)e^{−2iπf t}dt (2.5)

This allows every signal to be decomposed into a combination of sinusoidal wavesx(t), each with constant amplitudeajand frequencyωj[16].

x(t) =X

j

aje^iω^j^t (2.6)

The discrete implementation of the FT is given in Equation (2.7).

A_r=

N−1

X

k=1

X_ke^−2πirk^N r= 0, ...., N−1 (2.7)

WhereAris therth cofficient of the DFT,j=√

−1andXkdenotes thekth sample of the time series consisting of N samples.

The FFT iteratively calculates coefficients, combining progressively larger weighted sums of data samples to produce the resulting DFT coefficients seen in Equation 2.7 [17]. By leveraging the use of symmetries in the transformation, the FFT is a more efficient implementation of the DFT, effectively reducing the number of computing operations needed fromO(N²)toO(N logN). However, it is also important to note that as the FFT is an integral transformation that imparts a domain change from time to frequency, it is limited by the uncertainty principle. Thus, the determination of frequency is time interval dependent, imposing a strong limitation when the demand for both frequency and time resolution is very high[16].

2.2.1 Fourier Transform in the Power System

The FFT is as of today the standard method used for harmonic measuring and monitoring in the power system [4, 5, 6]. When applying Fourier series for harmonic measurement in the electrical power grid, it should be implemented as given by IEEE 519 and IEC 61000-4-7 [6, 4, 5] as seen in Equation (2.8) to (2.15).

f(t) =c0+

∞

X

m=1

cmsin(m

Nω1t+ϕm) (2.8)

(30)

cm=p

a²_m+b²_m (2.9)

Cm= cm

√

2 (2.10)

ϕ_m= arctan(a_m bm

) if b_m≥0 (2.11)

ϕ_m=π+ arctan(am

b_m) if b_m≤0 (2.12)

b_m= 2 Tw

Z Tw

0

f(t) sin(m

Nω₁t+φ_m)dt (2.13)

a_m= 2 Tw

Z T_w

0

f(t) cos(m

Nω₁t+φ_m)dt (2.14)

c0= 1 Tw

Z T_w

0

f(t)dt (2.15)

The signal is decomposed into a finite set of sinusoids of fundamental frequencyω₁and mean valuec_m. T_w is the width of the time window in which the Fourier transform is performed, while N denotes the number of fundamental periods within said window. The harmonic order, or the ordinal number related to the frequency basis(f = _T¹

w), is given bym. Note that as the measurements obtained in the power system are of discrete-time, the FT is not implemented exactly as in the equations above, but rather as an adapted version to fit the FFT previously mentioned. Digital measurement instruments based on 50 Hz frequency systems are recommended by IEEE Standard 519, when applying FT based techniques, to use a window length of 10 cycles i.e., Tw = 0.2[4]. This results in new spectral components every 0.2 seconds.As such, the method is not capable of instantaneous detection and will result in unreliable readings when time-varying frequencies or inter-harmonics are present. The effect of a time-varying frequency on the FFT is referred to asleakage. The energy at one point spreads to the adjacent frequencies, resulting in amplitude peaks spread over a larger frequency domain. In the presence of inter-harmonics i.e., harmonics given by non-integer times the fundamental, a distorted amplitude split between the closest integer frequency-bins can be observed. This is referred to as the picket-fence-effect [5]. To further illustrate the FFT capabilities, it is applied to several synthetic signals featuring different characteristics. The respective signals and their resulting transformation can be seen below in Figure 2.2 to 2.6.

(31)

2.2 Fourier Transform

Figure 2.2:FFT applied to a linear and stationary signalx(t) = 2cos(2π1000t).

Figure 2.3:FFT applied to a signal with varying amplitudex(t) = (1+cos(2π50t))cos(2π1000t).

Figure 2.4: FFT applied to a signal with increasing frequencyx(t) = 2cos(2π((z+ 1000)t)). z being the rate of frequency increase.

(32)

Figure 2.5:FFT applied to a signal with increasing frequency and varying amplitudex(t) = (1 + cos(π50t))cos(2π((z+ 1000)t)). z being the rate of frequency increase.

Figure 2.6:FFT applied two superimposed signalsx(t) =cos(2π800t) + 2cos(2π1000t).

(33)

2.3 Hilbert Transform

The uncertainty principle limits the time-frequency relationship in many of the methods used for time series analysis. Though this issue is addressed in versions such as Short Time Fourier Transform (STFT), where the spectral window can be chosen based on the required time and frequency resolution [16], or Wavelet Transform (WT) were the frequency and time resolution is relative to the frequency of the signal [3], they are not capable of accurate real-time tracking. Integral methods such as these are, by definition, not suited to deal with varying frequency, as the frequency is not a function of time within the integral limit.

The differential equations can only see frequency as a constant, not a varying parameter.

Thus, it is necessary to introduce new characteristics; Instantaneous Frequency (IF) and Instantaneous Amplitude (IA). By using the Hilbert Transform (HT) of a given time series, a strong analytical signal is generated, making it possible to calculate the IA and IF.

The HT is an integral transform defined in the time domain as a convolution between the Hilbert Transformer _πt¹ and a functionf(t)[18]:

fˆ(t) = 1 πP

Z ∞

−∞

f(τ)

t−τdτ (2.16)

Where P indicates the Cauchy principal value.

Contrary to other transformations such as the FT, HT does not involve a domain change.

Thus the resulting transformed time series is not equivalent to the original signal. How- ever, as HT can be represented by a phase shift of90^◦ to every Fourier component of a functionfˆ(cos(t)) =sin(t),(i.e., a linear filter shifting the phase by ^π₂ while maintaining the amplitude), four consecutive HT applied to a time series will result in the original representation of said time series. Figure 2.7 shows how the real part of the HT represents the envelope of the transformed signal.

Figure 2.7:A signal with varying amplitude and its envelope, calculated using the magnitude of the Hilbert Transform

(34)

The HTfˆ(t)is related to the original signalf(t)such that together they create astrong analytic signal, and can be expressed in exponential form as shown in Equation (2.17).

z(t) =f(t) +ifˆ(t) =A(t)e^iϕ(t) (2.17)

where

A(t) = q

f(t)²+ ˆf(t)² (2.18)

ϕ(t) =arctan(

fˆ(t)

f(t)) (2.19)

By expressingz(t)with an amplitudeA(t)as in 2.18 and a phaseϕ(t)as in 2.19, the IF can be expressed as the derivative of the phase, or asrate of change of phase, shown in equation 2.20.

ω(t) = 1 2π

dϕ(t)

dt (2.20)

As can be seen from Equation 2.17, a visualization of the HT can be achieved by viewing the resulting analytical signal as a vector rotating in the complex plane with angular speed ω(t)and magnitudeA(t). A perfect circle represents a mono-variable time series with constant amplitude and frequency, as a linear increase in phase at a constant magnitude results in smooth rotation and a constant IF and IA. Using an arbitrary linearly increasing phase (and magnitude) as a baseline, deviations from this line is the cause of change in the calculated IF. The instantly calculated frequency makes it possible to handle both varying and instantly changing frequencies in a mono-variable signal. However, given a multi-variable signal, the IF will represent a mix of all frequencies in the signal with increasing complexity depending on the number of components. It is thus commonly used in conjunction with EMD to separate the modes of oscillation before the HT. This process is referred to as the Hilbert-Huang Transform [19]. The definition of frequency as a time-varying parameter is, however, highly controversial. Time and frequency are gen- erally viewed as inverse quantities, thus making instantaneous localization of frequency ambiguous at best [16]. Nevertheless, HT has been shown to determine the spectral peaks for non-stationary signals quite well. Examples of the HT applied to different signals can be seen in Figures 2.8 to 2.13.

(35)

2.3 Hilbert Transform

Figure 2.8: The complex phasor representation, instantaneous frequency and instantaneous amplitude resulting from the HT applied to a linear and stationary signal x(t) = (1 + cos(2π50t))cos(2π1000t).

Figure 2.9: The complex phasor representation, instantaneous frequency and instantaneous amplitude resulting from the HT applied to a signal with varying amplitude x(t) = (1 + cos(2π50t))cos(2π1000t).

Figure 2.10: The complex phasor representation, instantaneous frequency and instantaneous amplitude resulting from the HT applied to a signal with increasing frequencyx(t) = 2cos(2π((z+ 1000)t)). z being the rate of frequency increase.

(36)

Figure 2.11: The complex phasor representation, instantaneous frequency and instantaneous amplitude resulting from the HT applied to a signal with increasing frequency and varying amplitude x(t) = (1 +cos(π50t))cos(2π((z+ 1000)t)). z being the rate of frequency increase.

Figure 2.12:The complex phasor representation, instantaneous frequency and instantaneous amplitude resulting from the HT applied to a signal with offsetx(t) = 2 + 2cos(2π1000t).

Figure 2.13: The complex phasor representation, instantaneous frequency and instantaneous amplitude resulting from the HT applied to two superimposed signals x(t) = cos(2π800t) + 2cos(2π1000t)

(37)

2.4 Empirical Mode Decomposition

The Empirical Mode Decomposition (EMD) is a data-driven time-series method capable of handling both non-linear and non-stationary signals. By considering oscillations on a local level, the method separates multi-component time-series into a finite set of modes, each representing zero-mean amplitude and frequency-modulated oscillations [20]. These modes are referred to as intrinsic mode functions (IMFs). Requiring no a priori defined basis system while fully capable of operating unsupervised, the method has proven to be highly effective [21, 22, 23]. However, due to its adaptive nature, the method lacks an analytical formulation, making theoretical analysis and performance evaluation challenging [24]. As mentioned in section 2.3, EMD is referred to asHilbert-Huang Transformation (HHT) when used conjunction with the HT.

2.4.1 Intrinsic Mode Functions

An IMF is a function that satisfy two conditions [19]:

1. In the whole data set, the number of extrema and the number of zero crossings must either equal or differ at most by one.

2. At any point, the mean value of the envelope defined by the local maxima and the envelope defined by the local minima is zero.

The conditions limiting the IMFs are essential for defining a meaningful instantaneous frequency and amplitude. The symmetry concerning the local zero mean, enforced by requiring the mean of the envelopes constructed by local extrema to be zero, changes what traditionally used to be a global requirement to a local one [19]. This reduces asymmetries in the resulting waveform, further limiting fluctuations in the instantaneous frequency. To remove all unwanted fluctuations and asymmetries, one could enforce the ideal requirement of ‘the local mean of the data being zero.’ However, it is immediately evident that non-stationary signals would require a local time scale for each local mean, essentially being impossible to define [19]. However, using envelopes as an approximation to the ideal, in combination with the first requirement, which is the traditional narrow-band requirement for a stationary Gaussian process, an IMF is not restricted to a narrow band signal.

Thus, an IMF can contain both amplitude and frequency modulation and can represent non-stationary signals [19].

2.4.2 The Sifting process

The process in which EMD decomposes a signal into a set of IMFs is referred to assifting.

The original signalx(t)is thus represented by IMFsxn(t)and a residuer(t), as can be seen in Equation (2.21)

x(t) =X

xn(t) +r(t) (2.21)

The decomposition satisfies the perfect reconstruction property as superimposing all IMFs and residue reconstructs the original signal without any loss of information [20].

(38)

Sifting is illustrated in Figure 2.14, and Figure 2.15 shows the envelope and mean calcu- lation for a given signal. The sifting process is summarized by the algorithm given below [25]

• Step 0: Initializen:= 1, r0=x(t)

• Step 1: Extract then-th IMF accordingly:

(a) Seth0(t) =rn−1andk:= 1

(b) Identify all local maxima and minima ofh_k−1(t)

(c) Construct forh_k−1(t)the upper envelopeU_k−1(t)defined by the maxima, and lower envelopeL_k−1(t)defined by the minima, by using a fitting interpolation method (traditionally cubic spline interpolation).

(d) Determine the meanm_k−1(t) = ¹₂(U_k−1(t)−L_k−1(t))of both envelops of h_k−1(t)

(e) Sethk(t) :=hk−1(t)−mk−1(t)

(i) ifh_k(t)does not satisfy IMF criteria, setk→k+1and repeat the process starting at point (b)

(ii) ifh_k(t)satisfy the IMF criteria setx_n(t) :=h_k(t)andr_n(t) :=r_n−1(t)−

x_n(t)

• Step 2: Ifrn(t)represents a residue, stop the process; if not setn→n+ 1and go to step 1.

Figure 2.14:Flow diagram of the EMD algorithm [25].

(39)

Figure 2.15: Example of the mean of a signal calculated by use of local envelopes based on local extrema.

EMD applied to a signal composite of three frequency components is illustrated in Figure 2.16. As expected the sifting process decomposes the signal into three different IMFs.

Figure 2.16:EMD applied to a signal composite of three different frequencies.

(40)

2.4.3 Stoppage Criteria

The sifting process serves two important functions; elimination of riding waves and smooth- ing of uneven amplitudes [19]. Riding waves is defined by [26]

In a singal, if there exists a local minimum greater than zero between two successive local maxima, or if there exists a local maximum less than zero between two successive local minima, the segment between these two successive local maxima (or local minima) is called a riding wave.

The first function translates into ensuing meaningful instantaneous frequency, while the second function handles cases of large amplitude disparity between waves. While these functions are necessary for the decomposition and extraction of IMFs, they can when overused rob IMFs of their physically meaningful characteristics [19]. As such, a signal may be fragmented into physically meaningless IMFs that do not hold any relevance to the original signal. Therefore, it is necessary to determine a stoppage criterion for the sifting process that ensues the IMFs retain as much physical sense as possible while keeping the mean close to zero. Some commonly used stoppage criteria are given below [16]:

1. Specifying a fixed number of siftings

2. Stopping whenxn, or residuern, is smaller than a pre-set value.

3. Standard deviation (SD) between two consecutive siftings is smaller than a pre-set value.

4. Using an evaluating function*

The SD in the 3rd stopping criterion is often referred to as the Rilling stopping criterion and is given by:

SD=

T

X

t=0

|(h_1(k−1)(t)−h1k(t))|²

h²_1(k−1)(t) (2.22)

A typical value given for the SD is set between 0.2 and 0.3 [19].

*An evaluating function proposed by [24] uses thresholdsθ1 andθ2. Defining a mode amplitude as in Equation (2.23) and an evaluating function as in Equation (2.24), sifting is performed whileσ(t)< θ1for a fraction of the signal duration(1−α), and changes to σ(t)< θ2for the remaining fraction. Thus the method handles globally small fluctuations while accounting for locally large excursions. Typical values for this method isα= 0.05, θ1= 0.05andθ2= 10θ1[24].

a(t) := U(t)−L(t)

2 (2.23)

σ(t) :=|m(t)

a(t)| (2.24)

(41)

2.4.4 End effect and boundary conditions

The estimation of envelopsU(t)andL(t)is a critical aspect of the EMD. As the IMFs are dependent on the mean calculated from these values, the technique used to connect extrema can have a tremendous impact on the decomposition quality of the EMD. The most common method used for connecting extrema is spline interpolation. Being piece- wise composed function of polynomial order n, splines are n-1 times differentiable at their knots acting as continuous functions [25]. The spline interpolation recommended by [19]

is thecubic spline, as though computationally costly, it usually yields the best approximation of the signal envelopes [25]. In some situations, other types of interpolation may be better as the vital aspect is how well it is matched to the original signal.

Using spline interpolation to construct envelopes can cause mismatch at the boundaries of the signal, causing increasingly larger fluctuations as the error propagates through the IMFs [25]. The effect is caused by splines treating the first and last data points as knots for both the upper and lower envelope, resulting in a mean with no physical meaning when close to the ends of the signal [25]. As this mean is used for further extraction of IMFs, improper boundary condition can impact all subsequent IMFs extracted. Over-sifting will cause these errors to propagate further into their respective IMFs and can, if not handled, render the resulting IMFs useless. An example of a boundary condition error propagating through IMFs can be seen in Figure 2.17.

(a)EMD of the non-stationary time series x(t) =sin(7t) +sin(4t) + 0.1t−1with improper boundary conditions. Here the first and last data point of the time series have been treated as knots of x(t).

(b)EMD of time series given in(a) using proper boundary conditions.

Figure 2.17:Example of improper and proper boundary conditions [25].

2.4.5 Mode mixing

The inability to distinguish between closely spaced spectral components is a common issue in spectral analysis. This is also the case for EMD, as mode-mode mixing can cause problems when decomposing a signal into a set of IMFs. However, an important aspect of EMD is the IMFs representation of physically important characteristics in the original signal [27]. Therefore, it is necessary to be aware of situations in which the fragmentation

(42)

of a signal by EMD can cause the results to lose their physical meaning. An example of this is the human perception of tones. When presented with two tones sufficiently close in frequency, the ear will perceive these tones not as two separate entities but as a single tone modulated in amplitude. This effect is commonly referred to as thebeat-effectand is shown in Figure 2.18. This can be expressed mathematically, as there is no difference in expressingx(t) =cos(2πf1t) +cos(2πf2t)asx(t) = 2cosπ(f1−f2)tcosπ(f1+f2)t.

As such, one can wonder if a decomposition into two modes yields a more accurate answer if the goal is a representation matched to perception rather than mathematics.

(a) Upper left: cos(2π30t). Upper right: cos(2π5t). Bottom: cos(2π30t) + cos(2π5t)

(b)Upper left: cos(2π30t). Upper right:

cos(2π35t). Bottom: cos(2π30t) + cos(2π35t)

Figure 2.18:Beat effect — When the frequencies of the two superimposed tones are sufficiently far apart(left) the two-tones interpretation is meaningful. When they get closer(right), an interpretation in terms of a single tone modulated in amplitude is favored [27].

To maximize the potential of the EMD, it is necessary to know how and when it manages to retrieve individual tones, and when it considers input as a single tone. Even if a full decomposition is not accurate in a physical sense, a thorough understanding of the process can be useful for understanding complex signals.

Expressing a two-tone composite signal by Equation (2.25)

x[n] =a₁cos(2πf₁+ϕ₁) +a₂cos(2πf₂+ϕ₂), f₁> f₂ (2.25) It has been discovered in [27] that the behavior of the EMD will depend on the relative parameters a ≡ ^a_a²

1,f ≡ ^f_f²

1 and ϕ ≡ ^ϕ_ϕ²

1. Thus the composite signal can instead be expressed as seen in Equation (2.26).

x(t;a, f) =cos2πt+acos(2πt+ϕ), t∈R (2.26) As the EMD extracts IMFs using a mean calculated by envelopes based on empirically located extrema, it seems rather intuitive that the EMD will notice only components im- pacting the final extrema in the superimposed signal. It is found in [27] that the contribu- tion of extrema is dependent on parametersaandf. As there are three possible outcomes i.e., clean separation of modes, no separation of modes, and something in-between, the

(43)

relationship is divided into three outcomes. This has been extensively investigated in [27].

Below is a brief visualized recount of the findings reported.

Case 1:af <1. Ifaf <1the extrema rate in the superimposed signal is the same as the high frequency (HF) component. This is illustrated in 2.19(a)(a). By plotting the derivative of the HF component versus the opposite derivative of the LF component, each crossing represents an extremum in the superimposed signal. Thus, for each pair of extrema in the high amplitude HF component, there is exactly one crossing i.e., one extremum in the composite signal. As the HF component has a consistent impact on the superimposed signal, the EMD is fully able to detect and extract the HF component [27].

Case 2: af >1∪af² <1. Ifaf >1∪af² <1there is no regular pattern to the rate of extrema in the superimposed signal. Looking at 2.19(a)(b) there is no guarantee for a regular distribution of crossings. In this case, the EMD issometimesable to detect the HF component resulting in IMF with partly separated modes.

Case 3:af²>1Ifaf²>1the extrema rate is the same as the low-frequency component (LF). Thus, the EMD cannot detect the HF component, extracting it as a part of the LF component. This is illustrated in Figure 2.19(a)(c), as there is one extremum in the superimposed signal for each pair of extremum in the high amplitude LF component. Thus, the superimposed signal’s extrema are only impacted by the LF component, hiding the HF component from the EMD.

(a)Each graph plots the derivative of the HF component and the opposite derivative of the LF component, each crossing representing an extremum in the superimposed signal.

(b) A 3D representation of when mode- mixing appears in IMFs.

Figure 2.19:The mode-mixing criteria visualized [27].

A visualization of the complete relationship betweenaandf is shown in Figure 2.19(b).

The ground level of the graph represents the complete separation of modes, while the top-level represents no decomposition. As the relationship betweenaandf changes, the degree of mode-mixing changes accordingly. Note that there is a drastic change in the EMDs’ ability to distinguish closely spaced spectral components as the amplitude reaches a threshold value, while a change in frequency results in a more gradual decline in decomposition quality.

Mode mixing due to closely spaced spectral tones can be a problem when applying EMD

(44)

Figure 2.20:Boundary condition map for mode-mixing in EMD. Red area represent mode-mixing, blue area represent clean separation of two signals [28].

to complex signals. A method for mode mixing separation is presented in [28]. By utilizing a 2D projection of the amplitude and frequency ratios given in aBoundary Map[28]

depicted in Figure 2.20 it is possible to construct amasking singal fm(t)that will mix withoneof the tones involved in a mode mix. Adding the masking signal to the original signal will result in a controlled artificial mode mixing, causing the extraction of a mode mixed IMF containing the masking signal and the intended mode. This removes the intended mode from the previous mode-mixed signal, leaving a pure IMF. Furthermore, as the masking signal characteristics are known, it can be easily removed from the new mode-mixed IMF by averaging the results of the EMD with a separate EMD using the opposite masking signal. This process is depicted in Equation (2.27) to (2.4.5).

x₊(t) =x(t) +f_m(t) (2.27)

x₋(t) =x(t)−fm(t) (2.28)

x(t) =Xxn,+(t) +xn,−(t)

2 +r(t) (2.29)

(45)

2.5 Real-time Empirical Mode Decomposition

Though the EMD is a sophisticated method with clear advantages over other empirical decomposition techniques, which are often limited due to conditions that only apply ap- proximately in certain situations [20], there are two apparent shortcomings. First, due to the empirical nature of the method, time series with a high sampling rate over a long duration, i.e., large amounts of data, can be a challenge. Current EMD techniques and computational capacities limit the complete analysis of large amounts of data requiring the segmentation of large data streams. Second, the traditional analysis requires the complete data set, thus the need to wait until the data collection is finished. This is indeed quite limiting if applying EMD for use in the power system as surveillance and diagnostic of current and voltage often requires a high sampling rate and real-time tracking. Thus, there is a need for a real-time version of the EMD capable of handling large amounts of data. In this section, two such methods are presented. These being the weighted sliding Empirical Mode Decomposition (wSEMD) [20] and its extension, the Online EMD [29].

2.5.1 Weighted Sliding Empirical Mode Decomposition

To handle the problem of computer memory and load when dealing with time series of considerable length, one could propose clean segmentation of the data stream. By sepa- rating the data into several distinct segments, applying EMD to each segment would significantly reduce the computational load. However, due to the boundary condition effect on the EMD, the resulting IMFs would have discontinuous and boundary condition errors propagating in the IMFs, as shown in Figure 2.21.

Figure 2.21:EMD decomposition of synthetic signalx(t) =cos(2π10t) +cos(2π5t). The series has been decomposed into four equal segments each with a separate EMD indicated by the color scheme. After adding the resulting IMFs, there are clear boundary errors at t = 1.25, 2.5 and 3.75.

SEMD preforms the traditional EMD or EEMD within a fixed window of sizem, shifting the window forward with step sizekapplying EMD for each shifting until the whole time series has been covered [20]. Bothmandkare chosen apriori, and should satisfy Equation (2.30) and (2.31) to avoid discontinues as in Figure 2.21 and to keep the window a multiple

(46)

of the step size.

k < m (2.30)

E= m

k ∈N (2.31)

This ensures that, if omitting the firstmpoints, there will beEestimations for each point in the original time seriesx(t)for each windowm_i. These estimated values may vary due to the boundary effect on each separate EMD, thus the mean value of each estimation is taken as the result of the SEMD. GivenjIMFs for each segmentm_ithe original time series is decomposed intojdifferent IMFsx_m_i(t)and a residuer_m_i(t)as can be seen in Equation (2.32).

xm_i(t) =X

xm_ij(t) +rm_i(t), (2.32)

Each IMF is then collected in a matrix of E columns such that each value in a given column corresponds to the same point in the original time series. As averaging the amplitudes for each respective point in the IMFs requires an equal amount of data in each column, the columns related to the beginning and end of the times series are not included in further calculations. This algorithm is illustrated in Figure 2.22. For each pointta mean value is calculated such that fort > m

xj(t) = 1 E

i+E−1

X

i

xm_ij(t) (2.33)

r(t) = 1 E

i+E−1

X

i

rm_i(t) (2.34)

i=t−m

k + 2 (2.35)

The resulting functionsx_j(t)are referred to as SIMFs andr(t)as the sliding residuum [20]. As averaging is used to calculate the resulting SIMFs, the number of IMFs uncovered in each window needs to be constant. i.e., chosen apriori and the and stoppage criterium is consequently limited to a fixed amount of siftings.

(47)

2.5 Real-time Empirical Mode Decomposition

Figure 2.22:Schema of the SEMD algorithm [20]

Though averaging the resulting IMFs from each step does hinder discontinues and reduce error, the resulting SIMFs are still affected by the boundary conditions. As illustrated in Figure 2.23, the ends of each IMF is significantly more defective when compared to the rest of the signal.

Figure 2.23:Error in resulting IMFs given signal as in Figure 2.21

This is usually not a severe problem using traditional EMD. However, when a time series is segmented, each segment has individual ends contributing to the boundary condition

Instantaneous Frequency Identification in Microgrids Through Adaptive Data Analysis

Master ’s thesis

Erlend Westad

Instantaneous Frequency Identification in Microgrids Through Adaptive Data Analysis

Erlend Westad

Instantaneous Frequency

Identification in Microgrids Through Adaptive Data Analysis

Master’s thesis in Energy and Environmental Engineering Supervisor: Olav Bjarte Fosso and Marta Molinas

June 2020

Norwegian University of Science and Technology

Faculty of Information Technology and Electrical Engineering

Department of Electric Power Engineering

Abstract

Sammendrag

Preface

Table of Contents

List of Figures

Abbreviations

Chapter 1

Introduction

1.1 Background and Motivation

1.2 Objective and Scope of Work

1.3 Structure of Thesis

Chapter 2

State-of-the-art in Methods for Frequency Identification in

Microgrids

2.1 Concept of Microgrids and Harmonics

2.1.1 Microgrids

2.1.2 Harmonics

2.2 Fourier Transform

2.2.1 Fourier Transform in the Power System

2.3 Hilbert Transform

2.4 Empirical Mode Decomposition

2.4.1 Intrinsic Mode Functions

2.4.2 The Sifting process

2.4.3 Stoppage Criteria

2.4.4 End effect and boundary conditions

2.4.5 Mode mixing

2.5 Real-time Empirical Mode Decomposition

2.5.1 Weighted Sliding Empirical Mode Decomposition