Non stationary-ARMA ﬂicker model for squirrel cage induction generators based wind farms

Time varying nature of the wind power is studied previously considering different time intervals from seconds to days. However for power quality problems such as ﬂicker, a model which considers the extremely fast variations is essential. Here by using large number of actual records, a time-varying model is proposed which considers the extremely short-time variations of wind active and reactive powers. The wind farm is modelled as a current source with time varying amplitude and phase which change every 0.01 s. Autoregressive moving-average (ARMA) models are utilized to model the variations and ARMA coefﬁcients are calculated for every record. Same to the actual behaviour, the proposed model is non-stationary as the ARMA order and coefﬁcients are different at every run of the model. The proposed model is conﬁrmed through utilizing several applications which need the power time series with extremely short sampling intervals. The following studies are performed by using the actual data and then the proposed model: power spectral density of active and reactive power variations, instantaneous ﬂicker, short term ﬂicker (Pst), estimating the Pst using the maximum value of the instantaneous ﬂicker, estimation of cumulative Pst for multiple wind turbines, and the impact of SVC on ﬂicker mitigation.


INTRODUCTION
In recent years, the importance of renewable energy sources has been increasing while the supply capability of fossil fuels began to decrease [1,2]. Wind energy is one of the most attractive renewable energy sources which have significantly considered all over the world [3]. There is the most abundant wind energy on the sea and among the various forms of wind energy; offshore wind farms have become more popular than onshore wind farms [4]. Increasing the penetration of wind power in the power system has been associated with many challenges. A main obstacle to the integration of wind power into the grid is its variability [5][6][7]. Assess and analyses as well as provide solutions for problems caused by wind farms, especially connected to the grid, require an appropriate and accurate model of wind generation systems.
Although there are many reports for flicker assessment caused by wind farms, yet there are few studies regarding wind farm models which can be used for flicker studies. In [8], the This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited. © 2021 The Authors. IET Renewable Power Generation published by John Wiley & Sons Ltd on behalf of The Institution of Engineering and Technology mechanical system of a class of wind turbine generators is modelled. The fluctuations of wind speed are modelled by a sinusoidal series with frequencies in range of 0.01-10 Hz. A wind model is based on a power spectral description of the turbulence, which includes the coherence between wind speeds at different wind turbines in a wind farm, is presented in [9]. Both Van der Hoven model and von Karman spectrum are used to model the long, medium and short term variations of wind speed in [10]. The wind speed is modelled using a sinusoidal series with 30 frequencies in the range of 0.001-10 Hz. The wind turbine current is utilized for calculation and assessment of the flicker in [11]. However, the wind turbine speed (or power) fluctuations are not modelled in [11]. An analytical formulation of the generated aerodynamic torque of a three-bladed wind turbine including the effects of wind shear and tower shadow is presented in [12]. The proposed model in [12] was utilized in several applications as the volt/var control of powerelectronic interfaced WTGs [13] and reactive power coordination in DFIG based wind farms [14]. A neural network based model for wind turbine flicker emission calculations is presented in [15]. The neural network is trained by varying all wind and network parameters that might affect the expected flicker levels. Wind speed time-series is modelled based on the Fourier synthesis method without using the actual wind speed data. Therefore, the output power obtained for turbines and as a result the calculated flicker cannot be indicative of actual and accurate quantities. In [16], the variations of the voltage are modelled by a white noise signal to study the flicker produced by wind farm, while voltage variations may not comply with white noise. The flicker caused by a wind farm is minimized in [17] by using an energy capacitor system. However the sampling time of the used wind power time series is not mentioned. The maximum allowable injected apparent power and number of connected wind turbines in a bus is calculated in [18].
The studies [8][9][10][11][12] are some of the well-known models which are suitable for modelling of wind farm flicker. However they are complicated, hard to implement and need a lot of mechanical and electrical details. In the other hand, the caused flicker is due to the stochastic and time-varying nature of the wind and its nonlinear effect on the wind power. The best presentation of a stochastic process with a semi-known PSD is utilizing time series models such as ARMA. Also the process is non-stationary and presents different behaviours in different times which are less focused in the previous studies. However, in the proposed models in [8,10] the wind speed was modelled as a sinusoidal series with some specific frequencies in range of flicker frequencies. Despite the mentioned studies, the fast variations of wind farms output powers in extremely short time periods have not been considered in the most of previous wind power models. The flicker frequencies range is between 0.5 and 25 Hz. So the sampling time of systems quantities such as active and reactive powers should be less than (0.02 s = 1/(2 × 25 Hz). However in this case, instantaneous voltage and current data with high sampling frequency are required [19].
Forecasting and modelling the time variations of wind farm quantities such as wind speed or wind power is always a hot topic and many studies recently were published in this regard [20][21][22][23][24][25][26]. However none of them were about the extremely short time forecasting or modelling this kind of variations. For example in [20] a short-term forecasting with 1 h time intervals is performed by using an artificial neural network. Also in [22] by 1 h time intervals a short-term prediction of wind farm is performed by using a support vector machine. A survey on wind power forecasting methods based on deep learning is performed in [22] where the time intervals are varying from 5 min to several days. A time interval equal to 15 min is used in [23] where the data quality is improved by using a clustering analysis data processing method. In [24] hidden Markov models are applied to model the wind farm power with 20 min time intervals. An empirical mode decomposition is used in [25] for wind power forecasting with 1 h time intervals.
In this study, with a novel approach, actual recorded data are applied to build a simple but accurate model for wind farms which considers the extremely fast variations of active and reactive powers. Wind farm is modelled as a sinusoidal current source with time varying magnitude and angle which are updated every half-cycle. In other words, the sampling time of wind active and reactive powers is equal to 0.01 s. A large number of the actual voltage and current waveform data was collected from a wind farm in Manjil, Iran. The data are gathered in different weather conditions during the winter and summer. In the first step, characteristics of time-varying nature of wind farm output powers are obtained using the recorded data. The most suitable ARMA model is determined and the ARMA coefficients are calculated for all data records. The reason for utilizing ARMA models in this study is to keep the proposed model, as simple as possible while preserve the stochastic nature of the process. In the second step, using the characteristics attained from the first step, a probabilistic-time series model is proposed for wind farms. In this model, variations of wind farm active and reactive powers are attained from ARMA time series which are specified in the first step. ARMA order and coefficients are changed in a probabilistic manner at every run of the model to preserve the non-stationary nature of the model. So at every run of the model, a different ARMA model is chosen randomly which results in a non-stationary model. Even in the case of similar ARMA model for two model runs, the output active and reactive powers will be different due to the existence of the white noise term in the ARMA models. The proposed model electrically is simple which consists of one sinusoidal current source with time varying magnitude and angle. The attained powers are used to determinate the values of the magnitude and angle of the proposed current source every half cycle.
The proposed model can be implemented in several applications which need the power time series with extremely short sampling intervals. The accuracy of the proposed model is proved by implementing its output powers time series on six applications and comparing the results in case of implementing the actual time series. The implemented applications are as follows: -PSD of active and reactive powers variations -Instantaneous flicker sensation -Short term flicker (Pst) -Approximation of Pst using the maximum value of the instantaneous flicker sensation -The relation between the Pst caused by a single wind turbine and the cumulative Pst for multiple wind turbines -Impact of SVC on flicker mitigation of wind farms The main contribution of the present paper can be summarized in proposing a practical model for squirrel cage induction generators wind farms based on a large number of actual data. It is simple for implementation and can be used in power quality studies instead of actual records which always are not available. Same to the actual cases the proposed model gives stochasticnon stationary outputs at every run.

THE ACTUAL RECORDS
In this research, a large amount of the actual instantaneous voltage and current waveform data was collected from Manjil wind farm in north of Iran (coordinates: 36 • 44′18.1″ N 49 • 23′51.5″ E. This plant, in total, includes 105 wind turbines with wound rotor induction generators rated at 660 KW and 46 wind turbines with squirrel cage induction generator rated at 330, 500 and 550 KW. The studies in the present paper are related to records gathered from squirrel cage induction generators and their related substation. The single-line diagram of the studied substation is shown Figure 1(a). Turbines are connected to dedicated transformers with ratings 2000, 1250, and 630 kVA, 20 kV/ 690 V, Y/∆ and then they are connected to medium voltage network via a switching substations. The short circuit power at PCC point is 500 MVA which presents a stiff 20 kV power system. Hence the adjacent loads are not affecting the wind farm time varying active and reactive powers. Each data record includes instantaneous voltage and current with 10 s time period and sampling time equal to 128 µs. The data are recorded through a portable recorder shown in (Figure 1(b)) which has the capability to record the instantaneous voltage and current with a sampling time equal to 128 µs [27].
Two sets of records are gathered. The first set is recorded from the 550 kW squirrel cage induction generators and consists of 50 three phase (150 single phase) records. The details of the 550 kW squirrel cage induction generators are presented in Table 1. Also Table 2 shows the details of the 10 s records from the 550 kW wind turbines.
The second set is gathered from a substation which feeds the squirrel cage induction generators and includes 109 three-phase (327 single-phase) records. The records details are presented in Table 2. Totally 20 wind turbines are connected to the substation, however in the time of recording, eight of turbines were out of service. The substation includes 12 squirrel cage induction generator wind turbine, rated at 330, 500 and 550 KW. The single line diagram of this substation is given in Figure 1(a). The data records are measured at the PCC point presented in this figure.
Active and reactive powers are calculated according to full cycle integration method which is updated every half cycle [28,29]. Equations (1) and (2) are realizing this method for the 128 µs sampling time where 156 samples exist in the 50 Hz   cycle.
where, n = 1,2,…,1000 and k = 1,2,…,78125 Equations (1) and (2) converts the v and i signals with the 128 µs sampling time to active and reactive powers time series with sampling time equal to 0.01 s. Figure 2 shows the calculated active and reactive power for a record from a single wind turbine and another one from the substation.

MODELLING OF WIND FARM'S ACTIVE AND REACTIVE POWERS USING THE ARMA MODELS
ARMA models are based on the idea that the present value of the time series can be explained as a function of the past values. ARMA models can be used for modelling and predicting the stochastic phenomena such as wind power. So far several studies have been done for modelling the wind power and wind speed, based on the ARMA model, considering different  [1,[30][31][32]. In this paper, based on the actual data records, a probabilistic-time series model is proposed to model the extremely short time variations of wind farm active and reactive powers. For this purpose, actual data records are considered in two cases: data records with 10 s length and data records with 1 s length which are attained by splitting the 10 s time series to ten 1 s time series.
An ARMA process of order (p,q) which is denoted as ARMA(p,q) is given as follows The process value at sample t is denoted by z t . a t is a Gaussian noise with zero mean and specified variance at sample t and a t −q is the value of a at sample t-q. Model orders p and q belong to autoregressive (AR) and moving average parts of the ARMA model. Coefficients 1 ,…, p , 1 ,…, q are the model parameters. In case of q = 0 the model is called Auto regressive (AR) and if p = 0 the model is called moving average (MA).
In this section the best ARMA model order is attained for active and reactive power time series. To this end, for every record first the union of sets 1-6 is considered as the set of candidate models to select the best ARMA order from it.
The Considered set: {Set 1, Set 2, Set 3, Set 4, Set 5, Set 6} So, totally for every record, 40 ARMA order candidates exist. Then the best ARMA order was specified by using the Schwarz criteria (SC) [33] as model adequacy checking test as follows: where, 2 a is the noise variance, m is the number of parameters estimated in the ARMA model that is equal to p+q and N is the time series length that is equal to 1000 and 100 for 10 and 1 s data record respectively. From the 40 candidate models, the model which minimizes Equation (1) will be selected as the most appropriate model for the considered time series.
After selecting the best ARMA modelling, in the next step by utilizing the MATLAB system identification toolbox, the ARMA model coefficients besides the residuals standard deviation are attained. Tables 3-6 show the selected ARMA orders and their related coefficients for records with length of 1 and 10 s belong the single turbines and substation measurements. Because of the limited paper size, only ten records are presented in Tables 1-4. However the enclosed MATLAB files include the ARMA coefficients details of all records. The number of records which is used in the proposed model is equal to 150, 1500, 327 and 3270 for 10 and 1 s individual turbine records and for 10 and 1 s substation records.

THE PROPOSED ELECTRICAL WIND FARM MODEL
Wind farm is modelled as a current source with time varying amplitude and phase which change every half cycle. Figure 3(a) shows the modelling procedure.
In this procedure, the most suitable ARMA model for active and reactive output powers is determined and the ARMA coefficients are calculated for all data records. So there will be a different ARMA model for every record and the total number of the attained ARMA models will be equal to the recorded data. At every run of the proposed model, an ARMA model is selected randomly from the attained set of ARMA models and active and reactive powers time series are simulated based on the selected ARMA model. They are used to attain the amplitude and phase (I and δ) of the current source by solving Equation (5) for every half cycle. Equation (5) simply can be attained by applying active and reactive powers balance in Figure 3(b).
where, P (k) and Q(k) are the wind farm active and reactive power consumption time series in the kth half cycle. E, R, and X are parameters of the distribution network Thevenin equivalent circuit. By summing square of the two equations in Equation (5), the resulted equation will be as follows: where: Equation (6) has two positive roots. Based on the actual records, the smaller root is the answer which is calculated by Equation (8): Then I(k) is calculated as the square root of the positive answer in Equation (6). Now (k) simply can be calculated by Equation (9): The modelling procedure is as follows: 1. Generate a number between one and the total number of samples, randomly. 2. Select the ARMA parameters set (coefficients, mean values of records and STD of residuals) in a row of parameters table that is corresponding to the generated random number. 3. Simulation of the P and Q time series using the proposed suitable ARMA model and the selected ARMA coefficients set. 4. Calculation of the amplitude and phase of the current source model of the wind farm by Equation (2).

ASSESSMENT OF THE PROPOSED MODEL
In this section, the accuracy of the proposed model is evaluated through some applications which need the wind farm's active and reactive powers data in extremely short time periods.
In the following applications first the actual recorded data are implemented then the corresponding output data of the proposed model are applied and results are compared.
Active and reactive powers are modelled for two categories of data. The First category includes active and reactive powers of individual squirrel cage turbines with 150 records and the second category consists of 330 records belong to the substation connected to the squirrel cage turbines.

5.1
Power spectral density where, M is the total number of records. Figures 4 and 5 show the resulted MPSD for the actual and simulated records (from the proposed model) for 1 and 10 s data records belong to individual turbines and the substation. These figures demonstrate that the PSD of simulated records is close to the PSD of actual records. This fact represents the high accuracy of the proposed model.

Instantaneous flicker sensation
Voltage flicker is one of the main power quality disturbances caused by the wind farms. One of the most important applications of the proposed model is modelling the flicker caused by wind farms. Here, flicker meter proposed in IEC 868 [35] which shown in Figure 6 is implemented in to calculate the flicker caused by the wind farm by using the actual data and the outputs of the proposed model. Figure 7 shows the calculated instantaneous flicker sensation based on an actual record and its corresponding simulated output from the proposed model. The short-circuit power is considered equal to 200 MVA.
In Sections 5.2-5.6, the active and reactive powers related to the wind farm substation records are multiplied by five and  Table 7 corresponding to each actual record, the model is simulated for 20× and the average values are written in the related tables. The first section of Table 7 shows the mean values of instantaneous flicker sensation Mean(S) for ten actual records and their corresponding simulated records by the   actual and simulated data corresponding to turbine records is 0.00220 and 0.00232 respectively. The same accuracy also can be seen for other turbine records and also the records regarding to the substation. These results show the appropriate accuracy of the proposed model for flicker analysis so, it can be used to evaluate and estimate the flicker emission of wind farms.

Short term flicker (P st )
In this section, the short term flicker is calculated using instantaneous flicker obtained from Section 5.2 based on IEC 868. The values of actual and modelled records are shown in the second section of Table 7. For example for the first substation record, P st is 0.1473 and 0.1504 corresponding to actual and the simulated data from the proposed model. Since the results are similar, it can be concluded that the proposed method for time series modelling has is a good accuracy for estimating P st .

Short term flicker using of maximum instantaneous flicker
In this section, the short term flicker is approximated using the maximum instantaneous flicker obtained from Section 5.2 by the following equation [36]: S max represents maximum instantaneous flicker. The P st approximated with Equation (12) using the actual and simulated data from the proposed model are depicted in third section of Table 7. The results again present the good accuracy of the proposed model. For example for the first turbine record, the P st estimated by Equation (12) is 0.0512 and 0.0529 for actual and modelled data respectively.

Estimation of cumulative P st for multiple wind turbines
When two or more wind turbines or wind farms operate simultaneously the cumulative P st at PCC can be estimated using Equation (13) [37].
where, P st,i is the flicker severity due to one wind turbine and n indicates flicker summation factor. In case of wind turbines with same rating, n can be calculated by Equation (14). n = ln (m) ln (P st cum ) − ln (P st ind ) (14) where, P st ind is flicker severity due to standalone operation of each turbine. In case of known P st of two or more turbines, P st ind can be approximated by the average value of the known P st of the single turbines.
As another application of the proposed model the aim here is attaining the flicker summation factor (n) in case of two wind turbines or two wind farms. The flicker summation factor (n) is calculated by applying the actual and simulated data. To attain the two simultaneously active and reactive powers time series regarding to the two wind turbines or wind farms, the active and reactive powers time series of one record was splatted in two other time series with equal lengths. The time series can be considered as stationary time series in such short time period. So, the attained time series can be considered they are taken from two adjacent wind turbines or wind farms. Now using the attained two time series, P st,1 and P st,2 are calculated. To calculate the cumulative P st , first the simple sum of the attained two powers time series in the previous step is calculated and so, the cumulative power time series are attained. Then the attained time series is used to calculate the P st cum . Now using Equation (13) n is calculated.
The results are presented in the fourth section of Table 7 reveals that the results corresponding to the actual records are close to the modelled records.

5.6
Impact of SVC on flicker mitigation of wind farms SVC is one of the most used reactive power compensators in wind farms [38][39][40][41] as it is cheaper than the STATCOM. Figure 8(a) shows an SVC connected in parallel with a wind farm besides its control system. The control system of the SVC consists of two major loops, that is, a feedback loop and feedforward loop. In the feedback loop which is supplied by reactive power from the upstream grid, a PI-controller adjusts the state error and so on, the average value of reactive power to zero. It is worthwhile to note that this block is dynamically slow and the delay is about 200 ms.
On the other hand, the feedforward loop uses reactive power of the wind farm to determine the reactive power value supplying by SVC. This loop compensates the fast variations of wind farm reactive power and thereby mitigation of voltage flicker. However the time delays relevant to the reactive power calculation and the inherent thyristor ignition limit SVC performance in following the fast variations of the wind farm [42]. Hence due to this time delay which is about 10 ms, SVC is unable to mitigate the flicker completely. To address this problem, a forecasting engine based on the selected ARMA models in the previous sections is added to the feedforward loop [42]. Utilizing the forecasted reactive power can prevent from reduction of the control system bandwidth causing by the time delay hence the inappropriate compensation of the reactive power variations.
As another application of the proposed model, the impact of SVC on flicker mitigation of wind farms is evaluated in this section using the actual data and the proposed model outputs. The instantaneous flicker sensation and the P st are calculated for actual and modelled records for the following three cases: Case 1: Without SVC Case 2: Using SVC without reactive power prediction The results of these three cases for an actual record and its corresponding model output are shown in Figure 8. Figure 8(b) shows the instantaneous flicker sensation signal for the three mentioned cases based on an actual record. It is evident that the SVC could reduce the flicker of the wind farm. However as it is clear from Figure 8(a), further reduction in flicker occurred when the reactive power forecasting engine is used inside SVC's control system. Figure 8(c) shows the instantaneous flicker sensation signal in case of using the modelled records. In comparison between Cases 1 and 2, Figure 8(c) same as Figure 8   shows the similar flicker reduction which proves the performance of SVC in flicker mitigation. Considering Case 3 again the same reduction in flicker is observed in Figure 8(c) similar to Figure 8(b). These similarities between Figure 8(a,b) reveals that the proposed model of wind farm can accurately estimate the performance of SVC in flicker mitigation. Furthermore, the short term flicker (P st ) for ten actual and their modelled records is presented in the fifth section of Table 7. Two type of comparisons can be performed in this section of Table 7. The first comparison can be between the three mentioned cases. Regardless using the results from the actual data or those extracted from the proposed model, results reveals that utilizing the forecasting of reactive power can significantly improve SVC performance in flicker mitigation. For example based on first actual record the P st is equal to 0.9921, 0.4560 and 0.1285 for Cases 1-3 respectively. On the other hand the second comparison can be between the results of the same case but between the actual and the modelled records. The results regarding the actual and modelled data are close to each other for the three cases. For example in case of using SVC with the prediction block regarding to record 1, the P st is 0.1285 and 0.1213 regarding to actual and modelled data.
So in summary the proposed model has the ability to evaluate the SVC performance in flicker mitigation of wind farms. Also, the results show that the SVC performance is modified considerably by implementing the predicted reactive power values compare to the regular method.

Numeric comparison with other methods
Here the results of applying the proposed model are compared with two other modeming approaches; sinusoidal series and grey model.

5.7.1
Sinusoidal series The active and reactive powers variations are modelled by a weighted sinusoidal series with Equation (15).
where, Δt = 0.01s, and k = 1 ∶ 1000. In this equation, there are 20 coefficients related to frequencies between 1 and 10 Hz which are calculated using the least squares method.

Grey Model
Grey models were used in many studies for prediction purposes. In the grey model GM(n,m), n denotes the order of the used differential equation and m is the variables number. Generally, GM(1,1) is used, as it is fast, accurate and requires less computation. The last section of Table 7 presents the resulted P st of the actual records besides the proposed model and the aforementioned modelling approaches. It is clear that compared to the sinusoidal series and grey model, the values regarding the proposed model are closer to the actual records.

CONCLUSION
A simple practical model which consists of a current source with extremely fast time varying amplitude and phase is presented. They change every half cycle (10 ms) based on appropriate ARMA models. The order and coefficients of ARMA models is changed at every run of the model. Hence same to the actual data the proposed model is non-stationary. Even when a same ARMA model is selected in two different runs, the output time series will be different due to the existence of the noise term in the model. The actual data and the outputs of the proposed model are implemented on several applications such as PSD of active and reactive power variations, instantaneous short term flickers, estimation of cumulative P st for multiple wind turbines and SVC impact on flicker mitigation of wind farms. Results confirm that the proposed model gives results similar to the actual records.