RoCoF Calculation Using Low-Cost Hardware in the Loop: Multi-area Nordic Power System

This research paper presents a rate of change of frequency (RoCoF) implementation using low-cost hardware in the loop (HIL) with application to Nordic Power System (NPS). Two methods to calculate the RoCoF are presented: Incremental difference one step and moving window (MW) or rolling window. HIL) approach is used to reach natural noise levels found in the frequency signal obtained from real power devices. The low-cost HIL implementation is based on two Arduinos® in the loop. An Arduino® working as standalone is used to generate an analogue signal representative of one area of the NPS, a second Arduino® in the loop contains the RoCoF calculation methods and filters. This paper analyses the effect of including lowpass filters (LPFs) because the noise of the frequency affects the RoCoF calculation. Three cases are implemented in the Arduino® and compared with Simulink®/MTALAB™ simulations: Base Case: no filter and noise signal, Case I: LPF applied to the input frequency signal, Case II: LPF at the frequency signal and LFP at the output of the RoCoF (smoothing). Results demonstrate that RoCoF calculation methods worked correctly in the HIL and using two LPFs clean the signal from the noise produced by the Arduino®.


I. INTRODUCTION
Nowadays, power system (PS) frequency and rate of change frequency (RoCoF) are growingly essential concepts for load management, protection and system control [1], [2]. Ideally, it is desired to have error-free and noise data that is available with minimal delay [3]. This is not the case in real PS, because the frequency and RoCoF measurement are sensitive to noise and disturbances in the electrical system. Typically to fix this, it is needed heavy filtering to deliver robust and usable data [4]. The RoCoF is an important indicator to evaluate the stability of electric power networks, and they are also used as input for the control and protection devices that are used in electric power systems (EPS) [5]. It is known that RoCoF in modern electric power networks increases as the inertia of the system decreases; this caused by the integration of wind and solar energy, which are energy sources without inertia [6]. RoCoF is the nearly all widely used anti-islanding protection method [7]. Nevertheless, the safety of the relays based on this scheme is continuously disputed [8], as it is susceptible to network disturbances, causing unwanted trips. The solution to this is to compromise the safety and reliability of the relay, this is not simple labour and is further fiddly when the efficiency of marketable available RoCoF relays varies significantly, even having similar configurations [8], [9].
This research paper presents a rate of change of frequency (RoCoF) implementation using low-cost hardware in the loop (HIL) with application to Nordic Power System (NPS). The implementations focus on two algorithms to calculate the derivative and the noise reduction of the signals (input and output) by the use of a lowpass filters (LPFs); these are implemented in low-cost hardware by combining Simulink ® /MATLAB ® . The impact of these calculations contributes to the safety with which the relay operates, and it is also possible to see how the disturbances affect the local frequency and the RoCoF in an EPS. The main contribution to the development of this paper is the implementation of the RoCoF calculation implemented in a low-cost HIL and applying LPFs to eliminate noise from the signals.

II. RATE OF CHANGE OF FREQUENCY (ROCOF)
Frequency drifts subsequently to the islanding be conditional on active power discordance, inertia among the time instance and island. Therefore, calculating the derivative of frequency, moreover called RoCoF, can be accelerated the islanding detection [10]. The variation of the load and input mechanical power produces a change of frequency. Thanks to this variation posterior to disconnect from the primary grid, a dynamic change in frequency is expected, besides the rate of change of frequency usually can be calculated by the following equation [11], [12]: where PL is the power demand, PG is the total active power generation in the PS, SGR is distributed generator rating, H is the total system inertia, f is the system frequency, and f0 is the rated frequency. Network incidents may also inflict changes in system frequency, culminating in the inappropriate operation of RoCoF relays [8]. The frequency protection relays are founded on RoCoF for their respond against incidents in the EPS, and the foundation of the method is the Fourier transform and the study of the zero-crossings. RoCoF relays measure the rate of change of frequency, and after the rate of change of frequency oversteps the preestablished setting, a trip signal is beginning. In respect of response time and to prevent false activations, the response time of the relay will be within 50 and 500 ms, this makes the device more secure [13]. The number of frequency points atop which the frequency change will be determined, define the measurement window, the interlocks need to be considered due to the low voltage since it could obstruct the activity of the frequency relay [8]. The block diagram, of a typical structure of a frequency relay, including the RoCoF action, is shown in Fig. 1. The initial block in the relay model is called frequency determination; It is where the frequency is calculated. Then an anti-aliasing LPF is applied to the frequency; anti-aliasing filter will commonly either allow some aliasing to take place or else reduced some in-band frequencies close to the Nyquist limit. The signal from the filter feeds to a RoCoF calculation block, and the output of the RoCoF block is eventually smoothed by the last LPF, and that respond is compared against value setting and commands trip just after the calculated values surpass pre-set threshold usually within 0.1 and 1.2 Hz/s, these parameters are mentioned in [14], [15].

A. RoCoF definition
RoCoF function is one of the protective functions used for fast load shedding, to detect loss of grid and to accelerate operation time in under and over frequency conditions. The magnitude of the power imbalance generated from the disturbance defines how fast the system frequency and generator speed changes. The RoCoF is measured in Hertz per second (Hz/s), this might indicate as a measure of the gravity of a disturbance. Disturbances that cause high RoCoF events can have detrimental effects on an EPS. When severe cases exist, the generator protection devices can operate and disconnect the generator from the EPS; synchronous generators can lose synchronism with the EPS, and the RoCoF anti-island relays can operate, disconnecting distributed generation from the PS [16], [17]. Therefore, they could initiate a cascading effect and aggravate the primary disturbance, and this could lead to the complete shutdown of the EPS. Consequently, it is essential that EPS are designed to operate in such a way that the severity of high RoCoF events is reduced or avoided [18], [19]. The most representative way of calculating the RoCoF in electrical power systems after a distributing is by using the centre of inertia (definition is presented in [20]). The RoCoF calculated as [21]: where Δ is the change in power, f0 is the system rated frequency, df/dt is the initial RoCoF, SGR is the rated power of the system and H is the remaining inertia constant of the whole system after the disturbance.

B. Methods to calculate RoCoF
In this scientific paper, two methods are used to calculate the RoCoF and explaining in the next subsections.

1) Incremental difference one step
Within the numerical derivatives, there is the method forward difference derivative, which is a simple approximate calculation; therefore it is only carried out the evaluation of the equation (3) for a small, but finite, ti+1−ti [22].
It is also known as the forward difference derivative. Specified the m(t,f) points, it is possible to evaluate df/dt, at (m-1) points applying the equation (3). The forward difference derivative could be converted into a backward difference derivative by changing to a negative value for ti+1−ti. This approach can approximate a limit by evaluating what is expected to be close to the limit and producing a difference in the numerator of the expression [23].

2) Moving window or rolling window
The effects of actual RoCoF measurements on a moving window can be analysed by implementing the following equation: where f/dt is the RoCoF at sample t, N is the number of samples in the moving average window, T is the duration of the moving average window and f(t) is the frequency at sample t.

C. Filter definition
The function of any filter is to select frequencies within the desired frequency range and to filter out any other frequencies.
Thus, using filter techniques is possible to obtain spectra that are free of any ambiguous signals and avoid misinterpretation in the spectra. The essential idea is that the raw signal, which contains noise and the desired information, will pass through a filter, then the filter output contains the desired information. The analog filters use analogue electronic circuits, and the digital filters simply require a digital processor to perform numerical calculations on the sampled values of an analog signal [24]. The advantages of digital filters are as follows: digital filters are programmable, especially they are extremely stable with respect to both time and temperature, a digital filter is much more versatile because signals can be processed in any variety of ways [24], [25].

1) First-order Filter
LPFs have been used in this research paper, as shown in Fig. 1. These filters take away part of frequency, generating a smooth response, besides providing small variations in the signal obtained and easy to appreciate the trend by increasing proportions in global signal and noise, with little loss at the output [26]- [28]. LPF, mainly moving average filters, also known as Savitzky Golay filters, commonly applied to eliminate unwanted sections from a signal, design decimators and interpolators, perform data averaging, discover patterns, and remove noise. The methodology to create filters infinite impulse response-based, incorporate the Chebyshev, elliptic and Butterworth [24], [29], [30].
As mentioned in RoCoF relay structure, the generated signal from RoCoF calculation is eventually filtered by a function, (T as+1) -1 , create to remove high-frequency transients. Thereby, Ta represents the time constant, used for the selected measurement window and time constant of the filters. The result is eventually compared against the threshold value. Thus, when the previous value oversteps the last, it starts a trigger signal. When a trigger decision is performed, the simulation takes it as the operating time. The LPF blocks are shown in Fig. 1 implement the following transfer function: where, s is the Laplace operator and Tx the time constant of the filter and related to the cut-off frequency. In measurement action, filters are employed the frequency be required to be into a narrow band, as well as the RoCoF have to not overstep the specific values for the system under study. Consequently, it is usually supposed that RoCoF and frequency and can be determined within 100 ms.

III. NORDIC POWER SYSTEM (NPS) MODEL
The model presented in this research paper is a simplification of the Nordic synchronous system. In Fig 2, the arrows indicate the power flow direction relationship and connections with the countries neighbouring the NPS. The NPS was modelled using one load and one generator by area. The DK2 area does not have hydropower, therefore was modelled as a constant value. The Nordic region possesses an interconnected synchronous system formed of eleven areas. These areas are SE1 to SE4 in Sweden, NO1 to NO5 in Norway, DK2 in Denmark and F1 in Finland. Denmark DK1 belongs to the Nord Pool system but not in the Nordic synchronous area [31]- [33]. For this model, the control structure is a procedure established from a merging of the DC power flow and oscillation equation. Thus, it is possible to expand the system to some arbitrary size. In each area of the electrical system, there are a certain number of energy charges and generating units. The system works at 50 Hz frequency, and the per-unit method is used in all calculations. Therefore, when the system frequency is 50 Hz the velocity will be zero, and working in steady-state. The apparent power of the system (Sbase) is 1,000 MVA; there is a global base value. The areas have local bases, associated with the power of the generators. Each modelled area have to the output the generator voltage angles. The turbine governors carry out the primary control in a PS, and the proportional integral derivative (PID) controllers are used to the governing systems. In Fig. 3 shows the block diagram used in this research paper. More details can be seen in [33].
The controller uses the speed change fixed point obtained by the AGC, ΔPref, the difference in reference speed ωref and the system speed, Δω, to the input. The integral time constant, Ti and proportional gain Kp is used to control the signal in the valve position of the turbine Δc [24]. The derivative of the PID controller is KD, in the model used in this paper was set to zero, as it was also used in [33].
During a power deviation produced by a disturbance, the controller is incorporated to make sure that the power production is the same as the power demand after the disturbance occurs. The deviation from planned power is calculated and is used as input to the PID integrator block. The droop rho (ρ) is the steady-state feedback which multiplies the deviation in the mechanical power, ΔPm, through a feedback loop.
The frequency is not an input to the turbine governor, hence does not possible that integrator recover to the reference value. To prevent generators in the system from operating affecting each other is to employ the droop feedback loop. Lastly, the water start-up time of the turbine is Tw. The block diagram of the turbine governor, turbine and the generator in area F1 in Finland, are shown in Fig. 3, the rest generator subsystems are identical, the difference lies in parameter values which depends on each specific area. The full details of the model implementation can be found in [34].

IV. HARDWARE-IN-THE-LOOP IMPLEMENTATION
Matlab TM has the Simulink ® support package for Arduino ® hardware; This can be used to develop and simulate algorithms that run independently on Arduino ® . It also has simultaneous link blocks for configuring and accessing the Arduino ® inputs and outputs and external mode for setting interactive parameters and monitoring the signal while the algorithm runs on the device [35]. A structure for the implementation of RoCoF using a lowcost HIL applied to the NPS model is shown in Fig. 4. In this scheme DUE1 is in standalone and contains the NPS of 11 areas, it outputs an analogue signal with the frequency. On the other hand, DUE2 receives the frequency signal and performs RoCoF calculations, which are displayed on the computer. Arduino ® Due was chosen for this research paper because it has analogue I/O, the processing and storage capacity is more significant compared to other Arduinos ® ; these characteristics are essential for HIL interaction. The communication of the DUE1 with another DUE2 is not direct, so it is necessary to add some operations so that the analog signals entering and leaving the Arduino ® are compatible with each other, as seen in Fig. 5. The changes in the system frequency of area 5 (Finland) is obtained from the NPS; it is necessary to add an offset to pass the signal to the positive side; this will depend on the size of this, 1.4 is used for this model. When the analog block is used, it accepts a value uint16 but considers only the 12 most significant bits for the conversion. This block emits voltage (Vout) on the DAC pin, such that:
substituting have: Moreover, before the Arduino ® sends the signal to the outside, an offset block of 1.4 must be added to have an entirely positive signal, a gain of (2 12 -1)/2.2; also a block of data type conversion uint16, as can be seen in Fig. 5. Fig. 6 shows the output signal from DUE1, has been called analogue frequency, and the fully simulated frequency signal in Simulink TM has been called the synthetic frequency.  The signal received by DUE2 in the analog input block (pin 0) defined as system frequency of area 5 comes from DUE1. If the measured voltage is equal to the analog reference voltage, the output of the block outputs 1,023. A gain of 3.3/1,023 is added to perform the coupling of the analog input signal with the RoCoF calculation blocks. The offset correction is added with a constant of 1.96, as seen in Fig. 7.

A. Base Case: No filtering
This section is the beginning of the cases established to simulate the derivative methods, it was called Base Case, because filters are not included in this model.

1) RoCoF1: Incremental difference one step
Included in this scheme is the frequency signal that comes from the NPS, area 5, the derivative block receives four inputs: frequency, the frequency with delay, time, time with delay, at the output a setting gain and the RoCoF, Fig. 7.

2) RoCoF 2: Moving window or rolling window
In the same way as RoCoF1, this scheme includes the frequency signal that comes from the NPS, area 5, the differentiating block receives two inputs: Frequency, Frequency with delay, at the output a setting gain and the RoCoF, Fig. 8. Fig. 9 shows an example of how the derivatives of the Base Case return the RoCoF signal without using filters, noise effects, and it is not possible to appreciate the expected signal.

B. Case I: LPF-Frequency
This section shows Case I, in this case, in addition to implementing the derivative, an LPF is also implemented for the frequency signal.

1) RoCoF1: Incremental difference one step
In this scheme, an LPF is added for the frequency signal. The derivative receives 4 inputs: Frequency, Frequency with delay, time, time with delay, at the output a setting gain and the RoCoF, Fig. 10.  Fig. 11 shows in black the RoCoF1 signal using a LPF that filters the Frequency signal. Besides, in red shows the RoCoF obtained in Simulink by the differentiating block, it has been called RoCoF Real. 2) RoCoF2: Moving window or rolling window In this scheme also an LPF is added for the frequency signal. The derivative receives two inputs: Frequency, the frequency with delay, at the output a setting gain and the RoCoF, Fig. 12. Fig. 12. LPF-Frequency RoCoF 2. Fig. 13 shows in black the RoCoF 2 signal using an LPF that filters the frequency signal. Besides, in red shows the RoCoF Real. From Fig. 11, and Fig. 13, it is possible to conclude the effect of filtering noise coming from the Arduino ® , unlike Fig.  10 where it is not possible to appreciate a signal, much less the behaviour, with the filter the signal is already possible to appreciate the behaviour, some of the noise has been removed.

C. Case II: LPF-Frequency + LPF-RoCoF
Finally, this section shows Case II, in this case, in addition to implementing the derivative, an LPF is also implemented for the frequency signal and another LPF at the output of the derivative block.

1) RoCoF1: Incremental difference one step
In this scheme, two LPFs have been added, one for the frequency signal and the other at the output of the derivative block Fig. 14.  2) RoCoF2: Moving window or rolling window In this scheme, two LPFs have been added, for the frequency signal and the output of the derived block, Fig. 16.  The effect of attenuating both signals, the one from the frequency and the other from the derivative block, gives a response with a much smoother signal, compared to the signals obtained in the models without filter and with a filter.

V. DISCUSSION
As can be seen in Fig. 5, the cyber implementation 1 is the NPS block diagram implemented in DUE1, working as standalone and is used to generate an analogue frequency signal, in Fig 6. from a representative of area 5 Finland, part of the NPS. As would be expected in Fig 6. the synthetic frequency signal implemented in DUE1 it contains noise, but it is still possible to appreciate the behaviour of the analog frequency signal. The base case filters are not included in this model; Fig. 9 shows the behaviour of the RoCoF1 and RoCoF2 models; the noise does not allow a visual comparison with real RoCoF. In case I, in addition to implementing the derivative, an LPF is also implemented for the analog frequency. As can be seen in Fig. 11, and Fig. 13, the LPF applied to the analog frequency, allows us to appreciate in more detail the behaviour of RoCoF1 and RoCoF2, the result is better compared to the base case, however, the noise continues present.
The Case II, the derivative is applied and an LPF is also implemented for the analog frequency and another LPF at the output of RoCoF1 and RoCoF2. The response to this implementation shows a very favourable performance in noise removal compared to the base case and Case I. Despite using low-cost hardware, the HIL simulation was successfully carried out, the test models, the RoCoF techniques, and the filters were implemented correctly, which could be verified with the results previously presented in Case II.
VI. CONCLUSIONS This paper presented the implementation of a RoCoF calculation algorithm using low-cost HIL. For illustrative purposes, the eleven areas model of the Nordic Power System (NPS) has been used. Two methods to calculate the RoCoF are presented: Incremental difference one step and Moving window (MW) or rolling window. HIL approach is used to obtain natural noise found in the frequency signal obtained from real power devices. The low-cost HIL implementation is based on two Arduinos ® in the loop, creating a cyber-physical system. Using two LPF helps to clean the signal from the noise coming from the DUE1; therefore, the expected signal can be appreciated in Case II. It was concluded that it is possible to implement HIL simulations in low-cost hardware, the correct design of the test models will ensure that the derivative methods work correctly calculating the RoCoF, the filters used in the analog signal and in RoCoF, ensure that the signal obtained is low in noise levels.