TIME SERIES MODELLING FOR WIND POWER PREDICTION AND CONTROL

(1)

TIME SERIES MODELLING FOR WIND POWER

PREDICTION AND CONTROL

CLUSTERING AND ASSOCIATION RULES OF DATA MINING FOR CFD AND TIME SERIES DATA OF POWER RAMPS

Dissertation in partial fulfillment of the requirements for

the degree of

MASTER OF SCIENCE WITH A MAJOR IN

ENERGY TECHNOLOGY WITH FOCUS ON WIND

POWER

Uppsala University

Dept. of Earth Sciences, Campus Gotland

Dept. of Engineering Science, Division of Electricity, Campus

Gotland

(2)

TIME SERIES MODELLING FOR WIND POWER PREDICTION AND CONTROL CLUSTERING AND ASSOCIATION RULES OF DATA MINING FOR CFD AND TIME SERIES DATA OF POWER RAMPS

Dissertation in partial fulfillment of the requirements for the degree of

MASTER OF SCIENCE WITH A MAJOR IN ENERGY TECHNOLOGY WITH FOCUS ON WIND POWER

Uppsala University

Department of Earth Sciences, Campus Gotland

Approved by:

Supervisor, Associate Professor Bahri UZUNO ˘GLU

Examiner, Professor Jens Nørkær SØRENSEN

(3)

ABSTRACT

Intermittent supply of wind power generation can cause power ramps which are sudden changes of power production in time. This is an important prob-lem in a power system that aims to keep the load and generation balance. Unbalance in the system can lead to price volatility and grid security issues that will create power stability problems and financial losses.

The objectives of this thesis are the improvement of the labeling procedure for extreme power ramp occurrences and the creation of a methodology to forecast power ramp occurrences in future to increase grid security and to help for power balancing procedure between different power production units (such as wind farms and hydroelectric power plants, in this study our case will be the investigation of power ramp balancing procedure between wind turbines) and investigate a test methodology with the help of multivariate meteorological data records to control for estimations’ accuracy of power jump problems that are referred as ramps. Final goal is to investigate and prove the relationship between the power ramp and kinetical energy of node in three dimensional space. Power ramp is defined through gradient of the power production. Power series for power ramp problem in this thesis is modelled with two approaches: First approach is time dependent, single point space cluster analysis, the second method is multi point spatial cluster analysis.

Once the power ramp data is generated through time series analysis by Autoregressive Integrated Moving Average model from SCADA power pro-duction data, near future ramp predictions are predicted. Based on this prediction, meteorological relationships which affect ramp occurrences are examined by association rule discovery technique. This technique is improved for better operational decision making.

To provide control room operators, validation methodology for the accu-racy of power ramp estimations, clustering analysis of data mining is em-ployed to see interrelations between meteorological conditions, space loca-tions and power jumps. Spatio-Temporal clustering analysis for ramp issues is introduced. Data mining methodology is employed with the K-Means clustering application in both space and time dimension.

(4)

better operational decisions.

In space scales, K-means algorithm is used to create spatial-temporal (PRR) based clusters to define most related clusters in space so that the impact of each spatial cluster can be introduced in positive and negative impact of each cluster.

Spatio-temporal clustering power ramp rate forecasting to the best of our knowledge is introduced for the first time. This technique can be used to estimate results of re-powering process or new wind turbine installation for total ramp rate characteristic of a wind farm.

AYYILDIZ,Bandirma site is used for this analysis. To validate the data with the model that will be used in clustering analysis, derivative of annual energy production result of WindSim simulations are employed to test sug-gested methodology accuracy. As a second control mechanism SCADA data is employed to create artificial test cases. WindSim results also generate the input data for space clustering.

(5)

ACKNOWLEDGEMENTS

I would like to thank my supervisor Bahri Uzuno˘glu for his guidance and his support. I would also like to thank TUBITAK,YEGM and Computational Renewables for the financial support and data that was provided for this study. I would like to thank my examiner for the constructive feedback on the thesis manuscript. I am indebted to Hasan Engin Duran who created an opportunity for me to study on wind energy and time series. I am indebted to my former lecturer, Turhan C¸ oban who taught me open source programming. Thanks also go to my friends and colleagues and the department faculty and staff for making my time at Uppsala University Campus Gotland a great experience.

I take this opportunity to record my thanks to Sabahattin Tolun and ¨

(6)

NOMENCLATURE

met mast: Meteorological measurement tower

SCADA: Supervisory control and data acquisition system KE: Kinetic Energy

P RR: Power Ramp Rate

P RRT : Total Power Ramp Rate τji: Stress tensor

Sij: Rate of stress tensor

u: x direction velocity component of wind speed v: y direction velocity component of wind speed w: z direction velocity component of wind speed ∂: Differentiation Operator

µ: Dynamic viscosity δ: Unit stress tensor ρ: Density

UTM: The Universal Transverse Mercator System ED: European Datum

AIC: Akaike Information Criterion SIC: Schwaz Information Criterion Σ: Summation operator

R

(7)

ACF: Auto correlation function

(8)

List of Figures

3.1 a) Wind Farm Layout with objects; triangles are turbines, circle is the met mast b) Power and Thrust Coefficient curves 7

3.2 Turbine details and locations. . . 9

3.3 a)Wind Rose b) Fitted Weibull Distribution of Wind Speed . 10 3.4 Terrain maps for AYYILDIZ a) Vector Map b) Elevation Grid c) Roughness Grid d) Slope Grid e) Slope Angle Grid f) Aspect Grid . . . 11

3.5 Main family of clustering algorithms(Gan, Guojun, Chaoqun Ma, and Jianhong Wu., 2007) . . . 13

4.1 state of art . . . 31

4.2 Power and Ramp (Power gradient) graph for January 2013 . . 33

4.3 a)Cluster number 2 b) Cluster number 3 c)Gradient of power that is ramp as observation number . . . 36

4.4 PACF graph for optimization of q . . . 38

4.5 Fitted value and SCADA data . . . 38

4.6 ARIMA fitted power versus measured power . . . 40

4.7 Future Prediction of Ramp . . . 40

4.8 a)Wind Speed b)Humidity c)Temperature d)Pressure . . . 41

4.9 Convergence results for spot values and residuals from sector 0 to sector 90 . . . 44

4.12 3-D illustration of space clusters . . . 52

4.13 2-D illustration of space clusters, height versus northing coor-dinate . . . 52

4.14 2-D illustration of space clusters, height versus easting coor-dinate . . . 53

4.15 2-D illustration of space clusters, easting coordinate versus northing coordinate . . . 54

(11)

List of Tables

3.1 Power and Thrust Coefficient curves . . . 8

3.2 2.5 M mesh Grid Resolution . . . 12

3.3 4.5 M mesh Grid Resolution . . . 12

3.4 Database . . . 26

3.5 Apriori Process . . . 27

4.1 Turbine 1 . . . 32

4.2 Suggested Cluster Numbers . . . 34

4.3 Cluster Number:2 . . . 35

4.4 Cluster Number:3 . . . 35

4.5 Forecasted Values . . . 39

4.6 Forecasted RAMP Cluster Analysis . . . 41

4.7 Association rules for 3 cluster case . . . 42

4.8 Setting for Simulations . . . 43

4.9 Predicted power comparison to Scada measured power . . . . 48

4.10 Predicted power comparison to Scada measured power for Tur-bines . . . 49

4.11 Index values for different cluster numbers(2) . . . 51

4.12 Ramp-Space Clusters . . . 55

5.1 Model Performance Comparison with Cluster Number=3 . . . 58

5.2 Ramp Estimation Accuracy Test Case 1 . . . 59

(12)

Chapter 1. INTRODUCTION

Four main objectives are followed through this study. First issue is to label extreme power ramp occurrences in the system and second one is the creation of a methodology to forecast power ramp occurrences to increase grid security and to help for power balancing procedure between different power production units (such as wind farms and hydroelectric power plants, in this study our case will be the investigation of power ramp balancing procedure between wind turbines) and investigate a test methodology with the help of multivariate meteorological data records to control for estimations of power jump problems that are referred as ramps. Final goal is to investigate and to prove the relationship between the power ramp and kinetical energy of node in three dimensional space. Power ramp is defined through gradient of the power production. Power series for power ramp problem in this thesis is modelled with two approaches: First approach is time dependent, single point space cluster analysis, the second method is multi point spatial cluster analysis.

Instant change in the power level is defined as ramp event, power level moves up and down in a stochastic fashion due to intermittency of wind speed. The main challenges with ramp are to detect the beginning and the end of the ramp event and to classify similar occurrences of the same phenomena. Ramp up refers to positive ramp issues and ramp down refers to negative ramp issues. Both can be defined through anomalies of the first derivation of power curve. In case of a negative ramp event, this can cause security and contingency issues on power system with financial impacts in reserve electricity market and the operators must have enough backup power to keep the load balanced. Inaccurate prediction of positive ramps is more important to negative ramp wrong predictions. Regardless of wind farm production, power is always positive or zero but PRR can be negative and positive due to change direction of power and larger absolute value of PRR refers to a the faster power surge (drop).

Power production time series for power ramp problem in this thesis is modelled with two approaches: First approach is from SCADA power data of each individual turbine with a clustering approach. The second method is from steady CFD solutions that uses time series data from meteorologi-cal mast, wind speed and direction for space relations of each turbine power output. This technique is improved for better accuracy on wind energy cal-culations for ramp modelling.

(13)

rules are investigated to define the relations behind meteorological events and ramps via association rules, Bandirma region of Turkey is chosen. The power ramp rate prediction estimation is tested for AYYILDIZ wind farm during February, 2013 for approximately three days. Power ramp rate char-acteristics of layout are defined and associated with space clustering analysis. Chapter 2 summarizes the literature review of power ramp estimations. Chapter 3 introduces input data sources (details of site data and CFD data generated by WindSim)and methodologies, formulations and algorithms of ARIMA, Autoregressive Integrated Moving Average and K-Means clustering study, Apriori Rule algorithm are summarized to provide the background for this study. Also the background of computational fluid dynamics, theoretical definitions and root studies are reviewed in this section.

Chapter 4 starts with time series data mining application and contin-ues on power forecasting for ramp prediction. The reader is introduced to Apriori rule set for the different meteorological variables as a forerunner to ramp. Preparation of input data, raw measured data is documented to cre-ate model input (cleaning, organizing, power ramp rcre-ate calculations, output of wind fields, database groups). Model and algorithm selections, model parameter optimization are conducted according to suggestions of literature and statistical tests (K-Means and Autoregressive Integrated Moving Aver-age). Herein also spatio-temporal clustering methodology will be introduced for wind power ramp forecasting to the best of our knowledge for the first time.

Chapter 5 introduces the control mechanisms and validation analysis of methodologies.

(14)

Chapter 2. LITERATURE REVIEWS

We can find direct applications of data mining of ramp events for load balancing and power generation issues of wind farms that will be reviewed below. The methods that have been employed from data mining can be summarized as feature selection, clustering and support vector machine ap-proaches. A literature review for the data mining tools in wind energy can be found in (Colak, Ilhami, et al., 2013).

Using feature selection techniques from data mining, Chandrika Kamath demonstrated that some atmospheric physical parameters are more impor-tant than others as can be expected. This offered the potential of employing data-driven predictive models for days with ramp events. This study also presents an importance level for different weather variables effects which can be caught by data mining selection structure. At the end of the study, au-thor announced a helper key set for control room operators. Control room operators have two basic problems related with ramp issue. Positive ramp when the wind speed increases and causes jump on produced power in a short period that can lead to financial losses and imbalance in the system (Uzunoglu, B. and Dervis., 2013), (Michel S. Menin, 2014), (Ulker M., 2011). The operators must reduce other generation (for example from hydro power plants, thermal power stations) to keep the load balanced. The challenging issue occurs here for the case of intermittent, unpredicted wind power ramp, the other generation sources cannot be reduced as a quick solution. Another challenge was identified as the positive ramp being inaccurately predicted. This will cause financial losses. In different regions, findings of the study has demonstrated different variables can be important for power ramp issues. It was also demonstrated, selected filtering technique has a considerable ef-fect on results. Three different filtering techniques; Distance, Chi-squared, Stump simulated different lists for same locations. In that study, data are taken from the Tehachapi Pass in Southern California (15 minutes intervals power generation during 2007-2008) and the Columbia Basin region (5 min-utes intervals power generation during 2007-2009) on the Oregon-Washington border. General outcomes of the study for TEHACHAPI PASS were that average humidity and wind speed are the main factors for ramp issues. For Columbia Basin wind speed and wind direction are important indicators for ramp issue.

(15)

and power - a clustering approach” (Erdem, Ergin, Jing Shi, and Yidong Peng., 2014). The wind power with wind direction and wind speed variables are used as input variables to model. Wind speed and wind direction are modelled to be merged to reach the power production. The authors have employed traditional ARMA(auto regressive moving average) (Uzunoglu, B. , 2012) and linked ARMA with clustering approach(k-means clustering) for wind direction prediction and they created a mixed ARMA model which has outperformed wind power estimation compared to traditional ARMA model. Another of the recent studies also from 2014 is focused on minimization the impact of wind power fluctuations on electricity grid (Guo,Zhenhai., 2014). As an application case, four different site data and four parameters from each of site (wind speed, humidity, pressure, temperature) are analyzed as input to K-means clustering so the authors created clusters and found clus-ter centroids in their second step. Authors used Apriori algorithm to figure out hidden rules between pressure clusters, temperature clusters, wind speed clusters and humidity clusters for each site. The chaotic time series forecast-ing model gave them 31 daily averaged wind speed estimation. Association rules are used for the correction of the wind speed estimation with mean wind speeds from clusters. Hexi Corridor of China is the case study location for several other studies conducted by the first author (Guo,Zhenhai., 2011) and other two research groups (Liu, Hui, Hong-qi Tian, and Yan-fei Li., 2012), (Zhang, Wenyu, et al., 2012).

There has been research in literature addressing modelling of time series using support vector machine method of data mining by (Kusiak, Andrew, and Haiyang Zheng., 2008). Multivariate time series models can be created for prediction of ramp events by this approach. Five different algorithms were used in this study. To reduce the dimensionality of the input and to enhance prediction accuracy, the boosting tree algorithm was employed as a preliminary process for this study. In that study, time series data from Supervisory Control and Data Acquisition (SCADA) is analyzed and power ramp rate (PRR) is defined for ten minutes intervals.

(16)

different initial conditions or different physical models to generate the un-certainty in the forecast to improve predictions. This approach will not be pursued in this study; however, instead spatio-temporal will be introduced.

(17)

Chapter 3. METHODOLOGY AND DATA

Test procedure and applications are validated with the AYYILDIZ case study. Herein we have three main inputs from the AYYILDIZ wind farm. First one will be introduced as Supervisory control and data acquisition sys-tem, SCADA data which refers as power output and second one will be me-teorological mast data which refers as humidity, wind speed, wind direction, temperature, pressure data and final input can be summarized as terrain maps, elevation data, roughness data.

3.1 Data Input 1: SCADA data

The AYYILDIZ project is a windfarm with 5 VESTAS V90-3.0 MW tur-bines that are located in the town of Bandirma which is closed to sea of Marmara as seen in Figure 3.1a with open sea in West and East directions. From Figure 3.3a the main wind direction is from North to North East in this wind farm.

(18)

(a)

(b)

(19)

Table 3.1: Power and Thrust Coefficient curves Wind Speed (m/s) Power (kW) Thrust coefficient

(20)

Figure 3.2: Turbine details and locations. WECS number: 1

(21)

3.2 Data Input 2: Met Mast data

Meteorological Mast data was provided by Tubitak. The site had a mast that had measurements on three heights at 50m, 65m and 80m. For the wind speed measurements, anemometer which is located at 80m and for the wind direction wind vane which is located at 50m are employed. Cleaning proce-dure of data is done due to NREL (Scott, George N., Dennis L. Elliott, and Marc N. Schwartz., 2010). Shadow effect due to the tower has been taken into consideration within anemometers tolerance. Shadow effects due to tur-bines are ignored since the main wind direction is coming from undisturbed sector. Contrary to the common case for many projects 12 months data with almost 100% availability was provided for measured 10 minutes av-eraged wind speeds with high maintenance practices. The data was not in raw data format but was in ascii format. The data provided exhibited wind behaviour from multi direction and sector. It can be seen from Figure 3.3(a). The given data is reorganized as time series data in the format of .tws file. 5 turbines are located in the given coordinates. Measured wind speed, wind direction data of wind farm is used for simulations wind speeds on turbine locations given UTM ED 50 35T. All the data in the project was given in datum UTM ED 50 35T.

(a) (b)

(22)

3.3 Data Input 3: Map

(a) _(b)

(c) (d)

(e) (f)

(23)

The details of terrain layout of the wind farm can be observed in Figure 3.4 which gives the terrain maps of the wind farm where the objects that are wind turbines and met mast are presented in Figure 3.1. The locations of individual turbines are recorded in Figure 3.2. The characteristics of power and thrust curves are recorded in Table 3.1.

The following details are the grid settings used for 2.5M mesh refined grid in Table 3.2 and 4.5 M mesh refined grid simulations in Table 3.3 while 3.5M mesh refined grid was also simulated.

Terrain Settings of 2.5M Mesh

Table 3.2: 2.5 M mesh Grid Resolution

Easting Northing z Total Grid spacing (m) 59.4-522.2 59.1-645.8 Variable -Number of cells 269 232 40 2496320

Terrain Settings of 4.5M Mesh

Table 3.3: 4.5 M mesh Grid Resolution

Easting Northing z Total Grid spacing (m) 43.7-446.1 43.7-538.5 Variable

Number of cells 362 310 40 4488800

3.4 Methodology Tool 1:K-Means

(24)

Hard clustering algorithms are further divided into hierarchical algo-rithms and partitional algoalgo-rithms. A partitional algorithm divides a data set into a single partition. In contrast hierarchical algorithm divides a data set into a sequence of nested partitions. Hierarchical algorithms can be further divided into agglomerative hierarchical and divisive hierarchical algorithms (see Figure 3.5).

Figure 3.5: Main family of clustering algorithms(Gan, Guojun, Chaoqun Ma, and Jianhong Wu., 2007)

Agglomerative hierarchical clustering starts with single cluster while merg-ing the closest pair of clusters accordmerg-ing to some similarity criteria till one cluster is achieved. Agglomerative hierarchical clustering are inflexible to data points incorrectly grouped at an early stage and different similarity measures for measuring the similarity between clusters can generate differ-ent clusters. Agglomerative hierarchical clustering is a bottom-up clustering method.

(25)

further sections. In the terminology of clustering, distances and similarities are the milestones for the process.

Similarity and dissimilarity measures

These crucial terms are passing in literature as similarity or dissimilarity measures defining quantitatively how similar two data points are or how similar two clusters are so here magnitude of similarity coefficient linearly proportional to two different point similarity degree. Some of the common distance and similarity measures:

(26)

Data Standardization and Transformation

To create dimensionless data from multi-dimensional data set like 3-D coordinates and time, standardization procedure is required. Knowledge of the location, scale units of raw data may be lost during standardization procedure; however, standardization of variables is obligation in cases where the dissimilarity measure is highly dependent to the differences in the scales of the input variables for example the Euclidean distance.(Gan, Guojun, Chaoqun Ma, and Jianhong Wu., 2007). There are various ways to perform this task in literature (Gan, Guojun, Chaoqun Ma, and Jianhong Wu., 2007). In this study z-score standardization method is used to scale input data.

scaled(x) = x− mean(x)

sd(x) (3.6)

sd() refers to standard deviation, mean() is the averaging function (Gan, Guojun, Chaoqun Ma, and Jianhong Wu., 2007).

3.4.1 K-Means clustering method

Clustering algorithms optimize an objective function F . Function value changes with partition of the number of clusters k as defined by C1, C2, C3, ..., Ck

F : Qk(Ω) → R (3.7)

Qk(Ω) is the set of all the partitions of data Ω = w1, w2, w3, ..., wm in K non empty cluster. The K-Means algorithm creates a solution to figure out locally optimal values via clustering criterion F which depends on the sum of the distance between each element and its nearest cluster center (centroid) (Pena, Jos´e Manuel, Jose Antonio Lozano, and Pedro Larranaga. , 1999). We can formulate it as follows where K is the number of clusters and Ki is the number of objects of the cluster i, wij is the jth object of the ith cluster and ¯wi is the centroid of the each cluster ;

(27)

The conventional k-means algorithm pseudo code (Gan, Guojun, Chao-qun Ma, and Jianhong Wu., 2007) can be summarized as below;

1) One of the main problems of clustering methodology is the decision of optimal cluster number which stands for meaningful partition size to investigate each observation without losing information in large clusters or creating too many small clusters. In the literature there is various ways to solve this issue, one of them is the entropy calculation which is not employed in this study and the second one is calculation and optimization of indexes (NbClust., 2014). In this study NbClust pack-age in R program will be employed to perform this optimization task. Majority rule, which is the most frequent result of 30 different index type, based cluster number suggestion will be used in our methodology. 2) Initialization partition (C1, C2,. . . ,Ck) of database is executed via Forgy

Algorithm (Gan, Guojun, Chaoqun Ma, and Jianhong Wu., 2007) 3) Compute centroids of each cluster.

4) Reassign wi to closest cluster centroid. 5) Recalculate centroids for each cluster ¯wi.

6) Reiterate until no further changes of cluster membership occur in a complete iteration and stop.

3.5 Methodology Tool 2:ARIMA

Methodologies, formulations and algorithms of times series analysis of Au-toregressive integrated moving average (ARIMA) and associate clustering methods of data mining will be summarized in this chapter to provide the theoretical background for the power ramp studies. These concepts will be later used in spatio-temporal dimensions.

3.5.1 Forecasting and time series analysis through

Au-toregressive methods

(28)

is the historical record of some certain activity, with measurements taken at defined regular time domains in the activity and without changing the method of measurement. Uncertainties of measurement technique and de-vice capabilities like power generation meter, anemometer sensitivity (cup anemometer,Lidar,Sodar etc.) effects method of measurement and time se-ries characteristics. Autoregressive models give us opportunity to predict power generation, future wind speed and direction so we can model power and we can make accurate prediction on wind resources on a site.

Short-term wind power forecasting can be made in the order of several days and also from minutes to hours. Usually, hourly forecasts of expected power are helpful in decision-making on unit dispatching, daily forecasts of hourly power are useful for the load scheduling strategy. Mainly four cate-gories are accepted in literature for time scales of power forecasting (Soman., 2010). For short term forecasting sub categorization should be useful.

i) Very short-term forecasting: From few seconds to 30 minutes ahead. ii) Short-term forecasting: From 30 minutes to 6 hours ahead.

iii) Medium-term forecasting: From 6 hours to 1 day ahead. iv) Long-term forecasting: From 1 day to 1 week ahead.

For the time scales that are used in the above list, AR (autoregressive), MA (moving average) and a combination of previous ones ARMA , ARIMA models will be introduced below for power forecasting needs.

AR(p) Model Form

AR model (Margherita Gerolimetto, 2010) is an univariate linear predic-tion model. A linear, homogeneous Gaussian process yt, has time reversible characteristic. Other accepted common models with non-Gaussian shocks are linear yet not time reversible. Gaussian Stochastic Process is a process that is completely specified through its first moments (mean, variance and co-variance). So, if mean and variance, covariance matrix is known, the process is known. This implies Gaussian processes have strict stationarity.

(29)

lagged value and moving average parameter. ACF and PACF both are de-fined through Variance and Covariance definitions.

ACF(Autocorrelation Function):In time series analysis, autocorrelation plot; is an outcome of the sample autocorrelations , versus h (the time lags).In other terms the ACF represents the randomness degree in a data set. For random data, autocorrelation function should give a result near zero for whole time domain. If non-random, then at least one of the autocorrelations will give a non-zero outcome. For a stationary process Yt, the autocovariance between Yt and Yt−j can be defined as,

γk = Cov(Yt− Yt−j) = E[(Yt− µ)(Yt−j− µ)] (3.9) and the autocorrelation function is

ρk = Corr(Yt− Yt−j) = γγk0 (3.10)

γ0 = V ar(Yt)

PACF(Partial Autocorrelation Function):PACF is the correlation between Yt and Yt−j after their mutual linear dependency on the intervening variables Yt−1, Yt−2 ..., Yt−j+1 has been removed. The conditional correlation will be,

Corr(Yt− Yt−j)|Yt−1, Yt−2, ...Yt−j+1) = φkk (3.11) which is usually referred as the partial autocorrelation in time series. For example,

φ11 = Corr(Yt− Yt−1) = ρ1 (3.12) φ22= Corr(Yt− Yt−2|Yt−1). (3.13) After defining the basic statistical derivations and definitions we will now introduce the autoregressive models. Autoregressive models give us oppor-tunity to have an idea for future power production.

If y is defined as random univariate time dependent variable and y(t) = y(t1), y(t2)...y(tn), y(t) series can be defined by an AR (p) model where p th is the order past data:

(30)

is error term, model deviation from actual value. t ∼ N(0,σ2) meaning that the t are identically and independently distributed each with a normal distribution having mean equals 0 and the same variance. Properties of error that is are independent of y.

Yt= φ1Yt−1+ φ2Yt−2+ ... + φpYt−p+ t (3.15) φj = φ11, φ22, ...φkk; φkk = φk for t=1,2,... (and conceptually t=0,-1,...) and with s ⊥ t for t 6=s yt is regressed values from p th order past data.The following back-shift operators and characteristic polynomial operators are introduced to define autoregressive processes:

• Operator B such that Bys = ys−1 and so Bk ys= ys−k for all k. • yt−p =Ppj=1ΦjBjyt+ t

• Φ(B)yt = t

• Characteristic Polynomial Φ(u) = 1 − Φ1u− Φ2u2− ... − Φpup. This is a polynomial of order p defined on _{|u| < 1.}

In Autoregressive AR(p) processes: the ACF function is infinite extended and PACF is infinite extended (to order p) for identifying an AR(p) process the PACF plays a very relevant role. Inversion and Moving Averages(MAs) are defined as follows. For MA(1):

yt = µ + t+ α1t−1 (3.16) µ is the mean of yt series.

yt =P∞j=1αt−j+ t (3.17) As in the AR(1) case, iterative substitution of yt−1 then yt−2 and so on in the right hand side of the model equation represents yt as a linear function of more distant past ys and t,t−1.... . Formally, if the inversion can be performed, it will be defined as

(31)

where α: is a weight and this can be further defined through an operator

yt= Φ(B)−1t (3.19) where operator B and characteristic polynomial are defined (Margherita Gerolimetto, 2010). This summarizes the explicit representation of yt as a linear combination of independent innovations, a linear (Gaussian) pro-cess. This is also referred to as a Moving Average (MA) representation. If the weights αj cut-off to zero after some finite lag which is defined as q, we have a finite M A(q) representation. In M A(q) processes: the ACF function is infinite extended (to order q) and PACF is infinite extended for identifying an M A(q) process the ACF plays a very relevant role.

Stationarity

A stationary AR(p) process is generally such a model in which inversion exists(Margherita Gerolimetto, 2010). The moving average weights must decay to zero eventually, otherwise the linear combination of past innovations will not converge.

ARMA Model

ARM A(p, q) models are defined by

Xt= Φ1Xt−1+ Φ2Xt−2+ ...ΦpXt−p

| {z }

autoregressive part

−θ1αt−1− ... − θqαt−q

| {z }

moving average part

(3.20)

where p = q = 1, then we have ARM A(1, 1) defined as

Xt= Φ1Xt−1− θ1αt−1+ αt. (3.21) ARIMA Model

(32)

smooth. If the sequence is not stationary, the smoothing of time series be-come essential for procedure. In reality, many atmospheric variable (wind speed, wind direction, temperature, pressure, humidity, thermal radiation etc.) are nonstationary in first moment (mean). For modelling mentioned atmospheric parameters, the source of variation must be removed. Gener-ally this treatment is done by differencing (I) the series (Friendly, Michael., 2002). As a de-trending treatment, regular differencing is applied together with AR and MA, they are referred to as ARIMA, I means ”integrated” and the differencing procedure (Friendly, Michael., 2002).

Suppose Xt is nonstationary in mean, to set up an ARMA model on the series wt,de-trending process is defined here as the result of the operation of differencing the series d times (in general d = 1): wt = δdXt. (Margherita Gerolimetto, 2010).

So, ARIMA models (where I stays for integrated) are the ARMA models (I=0 so data nature is stationary) defined on the d-th difference of the original process,main steps. Introduce orderly differential operator δy = 1_{− By for} the original non-stationary time series first order difference transform,

δyt− (1 − B)yt= yt− yt−1 (3.22) The original time sequence through d-order difference,

δdyt = (1− B)dyt (3.23) Difference series and final ARIMA equation becomes;

Φ(B)∆dyt = θBαt (3.24) where Φ(B)∆d _{is called generalized autoregressive operator and ∆}d_y

t is a

result after the differentiation operator cleaned the non-stationarity and now this time series can be modelled with an autoregressive- moving average methodology (Hui, Zhao, Li Bin, and Zhao Zhuo-qun., 2011).

Examples:

(33)

3.5.2 Statistical Processes That Impacts ARIMA

The following characteristic of time series will impact predictive perfor-mance of the AR, ARMA, ARIMA (autoregressive time series analysis) namely stability, stationarity and seasonality which is defined below.

Stability: A situation is all of the response parameters that we use to measure the process (first moments of data) that have both constant means and constant variances over time, and also have a constant distribution.

Stationary: A situation is a stochastic process which is based on probabil-ity distribution that does not change when shifted in time. The main impor-tant characteristic of stationary, parameters such as the mean and variance also do not change over time and do not follow any trends.

Seasonality: Cyclic variation.

A time series is considered to be stationary if the roots of the characteristic polynomial (Backshift Operators and Characteristic Polynomials) and are all greater than unity in absolute value. For example a simple AR(1); the characteristic equation is 1− θB = 0 that gives a root of B = 1/θ. If the root is greater than unity in absolute value, then yt is stationary. Thus, the AR(1) model is stationary when _{|θ| < 1. If the root equal to one, we call} such root as a unit root.(Halim, Siana, and Indriati N. Bisono., 2008). Unit root test includes DF (Dickey-Fuller) test, ADF(Augmented Dickey Fuller) test, PP (Philips Perron) test and so on (Gao, Shan, Yu He, and Hao Chen., 2009). For this project ADF test is selected as an analysing methodology. In ADF test, the same equation as DF test is defined (FinMath User Guide, 2014);

∆yt = α + βt + γyt−1+ δ1∆yt−1+· · · + δp−1∆yt−p+1+ εt, (3.25) where

- t is the time index,

- α is an intercept constant called a drift, - β is the coefficient on a time trend,

(34)

- t is an independent identically distributes residual term.

The focus of testing is whether the coefficient γ equals to zero, what means that the original ADF Time Series process has a unit root; hence, the null hypothesis of γ = 0 (random walk process) is tested against the alternative hypothesis γ < 0 of stationarity.

3.5.3 Model parameters p,d,q selection

Model parameter selection is based on Akaike Information Criterion (AIC) (Akaike, H. , 1974)

AIC = ln(ˆσ_ε2) + 2k

N (3.26)

and SIC Schwarz Information Criterion (Schwarz, G. , 1978) SIC = ln(bσ2

ε) +

k

N ln(T ) (3.27)

where σ_bε2 is the penalty term of process that is

b σ2 ε = 1 N N X i=1 (xi− ˆxi)2. (3.28) where xi is the observed data and ˆxi is the estimated data, N refers as total observation number. These values represent maximum meaningful lag value in AR variant applications. Herein k is the representative of total model parameters, in our case the model parameters

(35)

3.6 Methodology Tool 3:APRIORI Rules

It was already mentioned that association rules are the typical application of local-pattern discovery in unsupervised learning systems. Actually tricky name of data ”mining” is coming from this application. ”Rule” or ”rela-tionship” refers as gold. Here actually there is a mining act to find hidden ”knowledge”.

Data mining and discovery rule finding methodologies are highly depen-dent on input data types. For the association rule process, input data is typically transaction data which has time dimension, different objects and numerical equivalent. Transaction is data form of act which you can not di-vide to smaller parts. As an example in a certain time t pressure P, humidity H, direction D of wind are 3-items of a transaction and also as an individual P, pressure is a transaction that can be counted repeatedly but P/2 values can not counted this way.

For example; Positive power ramp rate values for each turbine during a year, negative ramp rate values for each turbine during a year is a transaction there is no way to say 0.5 positive ramp.

Basically for discovery of association rules in literature four common al-gorithms can be found.

(1) Apriori Algorithm (2) Eclat Algorithm

(3) Market Basket Analysis (4) WWW pathtraversal patterns

(36)

Artificial data set for explanation of computational concept will be intro-duced here. (P1,P2,P3,T1,T2,T3,H1,H2,D1) are the elements of the biggest cluster. In first step Apriori algorithm will count them, and according to pre-defined support value, it will filter and select most frequent itemsets.

Support in this context of association rules is the occurrence or the size of this sub-cluster in combination cluster of X and Y. For an ideal strong relationship between X and Y, support must be large and confidence must be high. According to occurrence of counted candidates, their support value is calculated by ratio between candidate repetition number to observation set, in our example it is 4 according to Table 3.4. In final step each candidate will be compared with their own support values to pre-defined threshold support value 50%. The candidates that are exceeding this threshold or catch this threshold will be selected. Here in this example of Table 3.4, we will set our support threshold to 50% for selection procedure. In the Table 3.5 reader can find item-set groups from 1 item set to 3 item set. After the generation of item sets, counting phase will be start to identify frequency of item groups in each item set class.

Confidence level in this context of association rules requires the knowl-edge of X combination Y clusters means, standard deviations for desired confidence level. The rule above X => Y which is defined in database DB with confidence c belongs to both X sub-cluster and Y sub-cluster. For a rule that implies H2, D1 => P 1, it is necessary that both itemset P 1, H2, D1 and H2, D1 are frequent so confidence of the rule is the ratio between two frequent item set support values. In our case it will be 0.50/0.75 = 0.66666. Desired association rules are defined with a confidence value c above a pre-defined threshold. Indeed in our process when we define a confidence level of example 0.1 only the rules which have greater confidence value form 0.1 are investigated. Recalling the previous example our confidence value was 0.66666 for the rule of H2, D1 => P 1. P1 is the right hand side of the rule that has a support value of 0.50, (H2, D1) is the left hand side of the rule has a support value, 0.75. If P 1 => H2, D1 is examined, this rule has a confidence level of 0.50/0.50.

(37)

conf (X− > Y ) = supp(X− > Y )/supp(X) P (XandY )/P (X) = P (Y|X)

lif t(X− > Y ) = lift(Y − > X)

lif t(Y_{− > X) = conf(X− > Y )/supp(Y )} conf (X_{− > Y )/supp(Y ) = conf(Y − > X)/supp(X)} conf (Y_{− > X)/supp(X) = P (XandY )/(P (X)P (Y ))} confidence is ratio between rule’s support and transaction’s support.

Herein performance measurement tool lift (C. Silverstein, S. Brin and R. Motwani, 1998) rules are also employed. Lift rule is a ratio which is defined by relationship of confidence and support values. If we calculate lift values of X = H2, D1 and Y = P 1. The H2, D1 => P 1 = (0.50/0.75)/0.50 = 1/0.75 = 1.333 and if we change the left hand side to right hand side P 1 => H2, D1 = (0.50/0.50)/0.75 = 1/0.75 = 1.333. Lift value of H2, D1 => P 1 equals the lift value of P 1 => H2, D1. But their confidence level is different. So when we set our desired lift level to a certain value only the rules which have greater lift value from our limit value will be selected and presented. In another words if we run our process with pre-setted c, confidence level=0.7 and l, lift level= 1.1 rule P 1 => H2 and D1 will be mined but H2, D1 => P 1 will not be presented or mined due to its low confidence level although it has a proper lift value. Greater magnitude of lift values are the proof for stronger associations between items (stronger rules).

Table 3.4: Database

Pressure Temperature Humidity Direction

P1 T2 H2 D1

P2 T3 H2 D1

P1 T1 H2 D1

(38)

Table 3.5: Apriori Process 1-item set Count s[%] selection

P1 2 50 P1 P2 1 25 P3 1 25 T1 2 50 T1 T2 1 25 T3 1 25 H1 1 25 H2 3 75 H2 D1 4 100 D1

2-item set Count s[%] selection P1,T2 1 25 P1,H2 2 50 P1,H2 P1,D1 2 50 P1,D1 T2,H2 1 25 T2,D1 1 25 H2,D1 3 75 H2,D1 P2,T3 1 25 P2,H2 1 25 P2,D1 1 25 T3,H2 1 25 P1,T1 1 25 T1,H2 1 25 T1,D1 2 50 T1,D1 P3,T1 1 25 P3,H1 1 25 P3,D1 1 25 T1,H1 1 25 H1,D1 1 25

(39)

If we summarize Apriori Algorithm, support threshold is s, DB is the transactions database, Ckis the rule candidate between X and Y sub clusters and k represents level, algorithm starts with support level entry and database, in the second step creation of item-sets is the main focus. From DB database transaction T with the level of k search enters a while loop.(Reminder to previous example:2-item setP1, T2 is a transaction with the level of 2). In while loop until k level reduced zero, Candidate Ck is investigated between X and Y sub clusters which are the members of DB database. If Ck rule is the member of combination cluster of a and b and also member of each one rule will enter a for loop to control support level threshold that is provided by Ck candidate. Here within a second for loop thanks to count[n] support value of candidate will be calculated.

(40)

3.7 Methodology Tool 4:Computational Fluid

Dynamics via WindSim

RANS models, Reynolds Averaged Navier Stokes equations are the mile-stones of computational fluid dynamics approach. For incompressible flows, the general form of the Navier-Stokes (NS) is given by

∂ ∂xj [ρuj] = 0 (3.29) ∂ ∂t(ρui) + ∂ ∂xj [ρuiuj+ pδij − τij] = 0, i = 1, 2, 3 (3.30) where ρ defines density, p defines hydrostatic pressure, uj defines velocity, u1 = u, u2 = v, u3 = w and xj spatial coordinates and the stress tensor τji depends linearly on the rate-of-strain tensor Sij,µ dynamic viscosity.

τij = −pδij + 2 µ Sij, (3.31) and Sij = 1 2 ∂ui ∂xj + ∂uj ∂xi . (3.32) Time Averaging

In an appropriate time T interval if the velocity u is time averaged:

u(x) = 1 T

Z t0+T

t0

u(x, t)dt (3.33) The averaged term and the fluctuation term of the Reynolds decomposition of velocity u = u+u0 that satisfies the properties u0 = 0, and u = u can simplify convective term in the momentum equation 3.35 that can be simplified as uiuj = uiuj + u0iu

0

j = uiuj − τijR where the fluctuation term is Reynolds stress tensor τR

ij = −u

0

iu

0

(41)

(42)

Chapter 4. APPLICATION OF THE

METHOD-OLOGY AND RESULTS

This section describes different sets of results, methodology, set-up and scope for a spatio-temporal power ramp data mining applications. In figure4.1, the flow chart of process can be seen.

Figure 4.1: state of art

4.1 TEMPORAL DATA MINING OF POWER

RAMPS VIA POWER FORECASTING:

Ramp Labels

(43)

P RR = |P ower(t + 10) − P ower(t)|

10 (4.1)

2013, January power production values of AYYILDIZ wind farm Tur-bine 1 presented in Chapter 3 and meteorological variables of this site was employed for PPR analysis while data was scaled for unit conversion as nec-essary(Gan, Guojun, Chaoqun Ma, and Jianhong Wu., 2007).

The PRR values will be directly fed as input to K-Means Algorithm. The details of Turbine 1 is presented in the below Table.

Table 4.1: Turbine 1

WECS number: 1 WECS name: Turbine1 Manufacturer: VESTAS Type: V90/3000 Nominal effect: 3000 Hub height: 80 x (local): 8876 y (local): 9482 x (global): 574923 y (global): 4467686

(44)

Figure 4.2: Power and Ramp (Power gradient) graph for January 2013

4.1.1 Optimization of Cluster Size

(45)

Table 4.2: Suggested Cluster Numbers

Index Optimization Method opt.cl. ”dunn” Max. value of the index 5 (Dunn 1974)

”sdindex” Min. value of the index 14 (Halkidi et al. 2000)

”dindex” Graphical method 5 (Lebart et al. 2000)

”sdbw” Min. value of the index 15 (Halkidi and Vazirgiannis 2001)

”ccc” Max. value of the index 3 (Sarle 1983)

”kl” Max. value of the index 10 (Krzanowski and Lai 1988)

”cindex” Min. value of the index 9 (Hubert and Levin 1976)

”db” Min. value of the index 14 (Davies and Bouldin 1979)

”silhouette” Max. value of the index 3 (Rousseeuw 1987)

”mcclain” Min. value of the index 2 (McClain and Rao 1975)

”ratkowsky” Max. value of the index 4 (Ratkowsky and Lance 1978)

”ball” Max. difference between hierarchy 3 (Ball and Hall 1965) levels of the index

”ptbiserial” Max. value of the index 4 (Milligan 1980)

(46)

4.1.2 Creation of Clusters with K-Means Algorithm

In the previous step, according to index values some critical cluster num-bers were determined. Here according those models, PRR clusters can be generated.

It is obvious that only positive and negative ramp definitions can be seen from cluster number 2 for PRR clusters as in Figure 4.3. In Figure 4.3, cluster number 3 is the final optimum case for our study, it is selected according to majority rule presented above. Here we should focus on characteristics of partitions. If there is an obvious jump between observations in each cluster it is a sign for ramp event as in Figure 4.2. For example if we inspect Figure4.3b, in cluster 3 there is considerable gap between observations of PRR near 500 MW/10 min node with the rest of the values in that cluster which means a power jump. Since ramp is the gradient of Power, a stable ramp value for a period will mean a jump and ramp as can be observed in Figure 4.2. The center of clusters are summarized in in Table 4.3 and 4.4.

Table 4.3: Cluster Number:2

Cluster Partition Size Centroid

1 3892 -8.845

2 564 60.791

Table 4.4: Cluster Number:3

1 3651 0.871

2 351 78.555

(47)

(a) Cluster size 2; x-axis PRR range, y-axis cluster number

(b) Cluster size 3; x-axis PRR range, y-axis cluster number

(c) Gradient of power that is ramp as observation number

(48)

4.2 RAMP CHARACTERISTIC

IDENTIFI-CATION VIA RAMP FORECASTING

(49)

Figure 4.4: PACF graph for optimization of q

3. Visual inspection is conducted for trend issue, optimization of d value, see Figure4.2.

4. The auto.arima forecast function in R is executed with settings p=5,d=1,q=3. Fitted value using these parameters and SCADA data for January 2013 is as follows.

(50)

Forecasted power as in Table 4.5 and in Figure 4.6 for 2.7 day with probability high P90, with low P80 intervals creates different scenarios for forecasted power.

Table 4.5: Forecasted Values

(51)

Fitted values through ARIMA are presented as below.

Figure 4.6: ARIMA fitted power versus measured power

A reminder is necessary here for positive ramp. Positive ramp is more desirable. For 3 cluster scenario we analyze to which clusters forecasted PRR values fall due to Low.P95 as presented in Figure 4.5 and Table 4.6. These clusters are predicted from ARIMA model.

(52)

(a) (b)

(c) (d)

Figure 4.8: a)Wind Speed b)Humidity c)Temperature d)Pressure Table 4.6: Forecasted RAMP Cluster Analysis

1 6 -8.80

2 1 -37.31

3 25 -2.35

4.3 ASSOCIATION RULES

Now we will use our meteorological mast data and cluster labels with the power ramp transaction database to figure out hidden rules.

(53)

this site.

Table 4.7: Association rules for 3 cluster case

Right-hand side Left-hand side Support Confidence Lift {temperature=2} => {clusters=1} 0.1052513 0.9474747 1.156381 {power=500} => {PRR=0} 0.1133303 0.8278689 3.115696 {power=500} => {clusters=1} 0.1342011 0.9803279 1.196478 {power=0} => {PRR=0} 0.1254488 0.9016129 3.393232 {power=0} => {clusters=1} 0.1389138 0.9983871 1.218519 {temperature=6} => {clusters=1} 0.1389138 0.8456284 1.032079 {PRR=0} => {clusters=1} 0.2657092 1 1.220488 {humidity=100} => {clusters=1} 0.2407989 0.8462145 1.032794 {power=500, PRR=0}=> {clusters=1} 0.1133303 1 1.220488 {power=500, clusters=1} => {PRR=0} 0.1133303 0.8444816 3.178218 {power=0, PRR=0} => {clusters=1} 0.1254488 1 1.220488 {power=0,clusters=1} => {PRR=0} 0.1254488 0.9030695 3.398714

4.4 SPATIAL DATA MINING OF POWER

RAMPS

(54)

4.4.1 Numerical convergence analysis

WindSim gives us an opportunity to select different solver types like: seg-regated solver and The General Collocated Velocity method, GCV multigrid solver which was used to generate the simulations due to reach secure faster convergence than Segregated solver. Settings for 15 different annual energy production cases including different vertical profile case are based on the be-low settings while turbulence was also modified to test different models. Our models are based on the following initial and boundary conditions logarith-mic law, at top fixed pressure, disregard temperature so there will not be a stability analysis. In the bottom our roughness and elevation data is the our boundary condition. Herein our purpose to figure out velocity components for each node.

Table 4.8: Setting for Simulations

Height of boundary layer (m) 500 Speed above boundary layer (m/s) 10 Boundary condition at the at top fix pres. Potential temperature No Turbulence model Standard

Solver GCV

Maximum iterations 500

(55)

(56)

(57)

(58)

4.4.2 Comparison of Scada data AEP to numerical AEP

prediction

The following measured value= 45477.76 MWh/year AEP was derived through SCADA supervisory control and data acquisition system.

The accuracy of the wake models are beyond the scope of this research. Default values for the wake models were used. Actuator disc computations, even though likely to improve results were not employed. For the calculations, different wake models are employed for the investigation of an effect of wake assumptions used in Table 4.9. The following list of acronyms will be used in the Table 4.9:

(i) Disregard Wake Model is the assumption which claims there is no effect between wind turbines. Each of them can take undisturbed clean wind without wake. So the highest prediction of annual energy production is occurred in this case.

(ii) Wake Model 1-Jensen Model Wake Loss Calculation (iii) Wake Model 2- Larsen Model Wake Loss Calculation

(iv) Wake Model 3- A turbulent depending rate of wake expansion

(59)

Optimized Wind Resource Map Grid Parameters

Table 4.9: Predicted power comparison to Scada measured power

(60)

These results are collected from WindSim energy module, frequency table annual energy production results. In table4.9 best estimation can be seen for the case of 4.5M mesh so we will take our wind resource map according to this settings.

Table 4.10: Predicted power comparison to Scada measured power for Tur-bines

Measured Power MWh/y 4.5 M Disregard Wake MWh/y

T1 9284.08 8827.60 T2 8643.90 8590.50 T3 9077.21 8919.70 T4 9244.19 9081.00 T5 9228.37 9400.10 %error T1 4.92 %error T2 0.62 %error T3 1.74 %error T4 1.77 %error T5 1.86

We must be careful in the analysis. In Table 4.9 although the best result seems to be given by 4.5M mesh resolution with air density correction to 1.1813 kg/m3 _{with wind field prediction error 1.449% compared to SCADA} data, we need to control every turbine with their original power estimation errors. In table 4.10 Turbine 1 is the challenge for WindSim estimations.

Now we have 4.5 M mesh std. K-Epsilon turbulence model, wake dis-regarded wind resource map at hub height. This means we have easting, northing, height and wind speed values.

(61)

To give control room operators; qualitative result from quantitative re-sults of different size of ramps, based on clustering analysis of data mining, interrelations between meteorological conditions, space locations and power jumps, Spatio-Temporal clustering analysis for ramp issues will be introduced employing the K-Means clustering.

An Euclidean Spatio-Temporal similarity or dissimilarity index to account for the distances in both space and time is defined below; however, since we are using steady state CFD simulations to diffuse the information from single point time series, the ratio for the time series between two points are fixed. This leads to simplification from Spatio-Temporal distance to spatial distance for approaches that use steady state CFD data as is the case.

4.4.3 Data Standardization and Transformation and

Similarity and dissimilarity measures

Herein xi _{= (x}i

1xi2xi3) represents three dimensional space for a point that is easting, northing height and yi

t represents time series of interest herein PRR for point i, t is the total time, C(xi_{, x}i−n_{) is the ratio between two} velocities that have n row differences at time t in space, one can observe that ut _{= (u}t

1ut2ut3) is not time dependent but dependent on only space locations. If we define the following relation for only two time steps

(62)

2, with easting, northing and height coordinates of five wind turbines and meteorological tower. Turbine hub height velocity is mean measured velocity at 80 meters during 2013. Different variables of standardization of data must be performed. Here z0 standardization or scale formula is employed to do the scaling task as discussed in Chapter 3.

If we think about power ramp rate definition we can observe power ramp rate is the derivation of power and power depends on cubic power velocity of wind speed values. So PRR depends on square power of velocity in other words Kinetic Energy. Because of this as a representative of PRR, square power of velocity values in our input data set will be used.

deuc = q

scaled((xi_{− x}i−n₎2_{) + scaled((C(x}i_{, x}i−n_)(µi₎2_{− (µ}i−n₎2₎2_{) (4.2)}

4.4.4 Optimization of cluster size

Following the same procedure to determine cluster size, the results of the optimization set on cluster test via NbClust function of R based on K-Means algorithm is presented below.

Table 4.11: Index values for different cluster numbers(2)

KL CH Hartigan CCC Scott Marriot

Number clusters 9 5 5 20 4 4.00E+00

Value Index 31.7281 3168.937 821.2199000000001 64.3785 7075.753 2.49E+14

TrCovW TraceW Friedman Rubin Cindex DB

Number clusters 3 3 20 18 4 4

Value Index 11407328 1715.817 4.3621 -0.4578 0.2036 1.0783

Silhouette Duda PseudoT2 Beale Ratkowsky Ball

Value Index 0.3523 1.3209 -595.96 1.0151 0.3968 4453.005

PtBiserial Frey McClain Dunn Hubert SDindex

Value Index 0.5846 NA 0.697 0.0071 0 1.6511

SDbw Dindex

Number clusters 20 0

Value Index 0.1525 0

(63)

Figure 4.12: 3-D illustration of space clusters

(64)

Figure 4.14: 2-D illustration of space clusters, height versus easting coordi-nate

Figure 4.12 will create the opportunity to locate useful ramp character-istics of a large region, thanks to it, a ramp optimization between different wind farms can be done successfully.

For a case study we will investigate AYYILDIZ wind farm layout, herein optimum cluster number is defined with the gap index(available in NbClust package) due to avoid problems which occur because of the sample size. It suggested us to select 3 as an optimum number of cluster.

(65)

(66)

Table 4.12: Ramp-Space Clusters

Item Easting Northing Hub Height-Height V2 _{Cluster Label}

(67)

Chapter 5. DISCUSSION AND ANALYSIS

In this chapter we will introduce our control mechanisms. We will start with decision rules which will help to re-consider positive and negative power ramp estimations. Later, affect of wind direction data and clustering methodology selection will be summarized. Finally the reader can find root mean square error calculations of suggested space-ramp cluster methodology application in AYYILDIZ wind farm.

5.1 Apriori Rules:Ramp Estimation Controller

In control room there will be different sets of meteorological variables if the operator considers association rules to act according to future ramp es-timations or to eliminate the suspicious eses-timations like a positive ramp estimation with low temperature forecast with rain expectation.

If we try to summarize rules in table4.7 for an operator that would like to use these results for decision making on ramp ranges, one can conclude cluster number 1 partition has the biggest size among other clusters so more observations are clustered there. Also another characteristics low tempera-ture value is seen in cluster 1. It is logical because we used January data so it can be interpreted as considerable numbers of observations have low tem-perature characteristics. When the power reaches 500 MW, PRR comes as 0 for this site when the power reduced to 0 so the case of storm or maintenance directs us to also PRR=0 value.

For the times relative humidity reached 100% value one can conclude, ramp ranges will be collected in cluster 1. So there should be rain or snow is-sues. So if one reassess again Figure4.3b, information for ramps are clustered in 1.

(68)

5.2 Residuals for forecasted power for future

ramp estimation

As we introduced earlier for future power ramp estimations we start with the forecast of power. Model optimization algorithm is based on a loop and here reader can see the residuals of power forecast is quite close to white noise. Mean of residuals are calculated as -.004 which is close to zero.

Figure 5.1: residuals of power forecast with ARIMA(5,1,3)

5.3 Model sensitivity to input data

In the beginning of study, site manager gave uni-directional wind vane data and our time series climatology input files are employed this informa-tion. Without changing wind speed data just a treatment on direction data caused a % 10 to 15 overestimation results compared to measured value. It is an important result which shows us selection of grid in space-ramp clusters depends on also meteorological mast data with a significant importance level.

5.4 Selection of clustering methodology

(69)

H(C) =− K X k=1

P (k) log P (k) (5.1) where C: Cluster Number, k: variable for our case PRR clusters (Kim H and Park H , 2007). The comparison is presented for a fixed cluster size in the below Table in terms of entropy value and computational time. While creating more information through agglomerate clustering, computational time is increased four folds.

Table 5.1: Model Performance Comparison with Cluster Number=3

Entropy CPU Time K-Means 0.589564 + Agglomerative 0.636863 ++++

After a few tests, it was observed that K-means require much lower com-putational time requirements for this problem.

5.5 Root Means Square Error investigation

for space-ramp clusters

In figure4.12 Turbine 1 and Turbine 2 are in same cluster. If we assume non-existence of Turbine 2, before installation and powering of Turbine 2 we can estimate PRR effect of Turbine 2 to total wind farm PRR values with Turbine 1. RM SE = s_P n i=1(P RRT eˆ i− P RRTi)2 n . (5.2)

(70)

difference (P RRT eˆ i− P RRTi) is defined by the distance between estimated total power ramp from cluster analysis suggestions and calculated power ramp rate from measured power production SCADA data. The index e is the estimation representative. TheP RRT eˆ i estimates total power ramp rate in a certain i th step. The P RRTi is calculated total PRR values from measured power production from SCADA data, total power ramp rate in a certain i th step.

Table notation below uses P RRT − P RR2 + P RR1 = P RRT e1 where test case is PPR from first turbine which is (PRR1) and suggestion of cluster ramp analysis is PRRTe which is (predicted total ramp rate of wind farm).

• P RRT − P RR2 + P RR3 = P RRT e2 • P RRT − P RR2 + P RR4 = P RRT e3 • P RRT − P RR2 + P RR5 = P RRT e4

Table 5.2: Ramp Estimation Accuracy Test Case 1

PRRT-PRR2+PRR1 PRRT-PRR2+PRR3 PRRT-PRR2+PRR4 PRRT-PRR2+PRR5 RMSE 170.2391 174.6876 177.0265 196.7343

Turbine 3 and Turbine 4 are in same cluster. If we assume for the case of non-existence of the Turbine 4 , before installation and powering of Turbine 4 we can estimate PRR effect of Turbine 4 to total wind farm PRR values with Turbine 3.

PRRT-PRR4+PRR1 PRRT-PRR4+PRR3 PRRT-PRR4+PRR2 PRRT-PRR4+PRR5 RMSE 188.8657 157.0205 177.0265 157.77 98

(71)

members address neutralization of ramp issue. So two different turbines can be used as balancing tool.

In another way to examine it we can use direct RMSE values between PRR values of each turbine.

PRR1versusPRR2 PRR1versusPRR3 PRR1versusPRR4 PRR1versusPRR5 RMSE 170.239 174.97 188.96 202.84

PRR4versusPRR1 PRR4versusPRR2 PRR4versusPRR3 PRR4versusPRR5 RMSE 188.965 177.0264 157.020 157.779

The following list summarizes our approach. 1. With climatology data run windsim. 2. Export rotor profiles.

3. Create easting, northing, hub height, square wind speed input data for layout (only items turbines and met tower).

4. Scale multivariable data. 5. Create clusters.

6. Locate each turbine clusters

7. Define PRR values of each wind turbine

(72)

Chapter 6. CONCLUSION

A new methodology is introduced to label power ramp occurrences to point extreme power ramp rates and clustering them according to their significance level. Power ramp rates of measured power data is calculated and then with the help of clustering analysis and K-Means algorithm power ramp rates are classified. For the test case extreme positive ramp occurrences are labeled as 2, extreme negative ramp occurrences are labeled as 3 and remaining ramp occurrences are clustered in cluster 1.

For temporal data mining of power ramps, ARIMA is employed to fore-cast short term power production and from this knowledge, PRR values are calculated and clustered. From previous PRR cluster characteristics, a methodology is created for possible future PRR occurrence ranges and exe-cution plans for operators based on data mining association rules discovery technique. Future power ramp estimations for approximately 2.5 days are done for test case.Operational guidelines based on atmospheric variables for power control room operators are introduced, as a control mechanism of power ramp rate estimations, association discovery rules are employed.

AYYILDIZ BANDIRMA, TURKEY test case was investigated for the year 2013. Annual energy production calculation by WindSim has underes-timated real total production. Different turbulence models like, YAP cor-rection, RNG are employed to see difference between standard K-epsilon turbulence models. 4.5 M mesh K-epsilon disregard wake model gave the closest annual energy production estimation to our measured annual energy production value.Best estimator grid features are selected as settings of Wind resource map at hub height (80 meter) is exported from WindSim to be later used in clustering analysis.

(73)

(74)

Bibliography

A. K. Jain, M. N. Murty, and P. J. Flynn. Data clustering: a review. ACM Comput. Surv., 31(3):264-323, 1999.

Akaike, H. A new look at the statistical model identification, IEEE Transac-tions on Automatic Control AC-19: 716-723. 1978.

Ball GH, Hall DJ . ”ISODATA, A novel method of data analysis and pattern classification”. Menlo Park: Stanford Research Institute. (1965).(NTIS No. AD 699616).

Blonbou, Ruddy, St´ephanie Monjoly, and Jean-Fran¸cois Dorville. ”An adap-tive short-term prediction scheme for wind energy storage management.” Energy Conversion and Management 52.6 (2011): 2412-2416.

C. Silverstein, S. Brin and R. Motwani, Data Mining and Knowledge Discov-ery 2(1), 39 (1998), DOI: 10.1023/A:1009713703947

Colak, Ilhami, et al. ”A data mining approach: Analyzing wind speed and insolation period data in Turkey for installations of wind and solar power plants.” Energy Conversion and Management 65 (2013): 185-197.

Davies DL, Bouldin DW . ”A cluster separation measure”. IEEE Transactions on Pattern Analysis and Machine Intelligence, 1(2), 224-227.(1979) Dunn J . ”Well separated clusters and optimal fuzzy partitions”. Journal

Cybern, pp. 95-104. (1974)

Erdem, Ergin, and Jing Shi. ”ARMA based approaches for forecasting the tuple of wind speed and direction.” Applied Energy 88.4 (2011): 1405-1414. Erdem, Ergin, Jing Shi, and Yidong Peng. ”Short-Term Forecasting of Wind Speed and Power-A Clustering Approach.”Proceedings of the 2014 Indus-trial and Systems Engineering Research Conference.(2014)

FinMath User Guide,”Augmented Dickey-Fuller (ADF)

Test”.,https://rtmath.net/help/html/93a7b7b9-e3c3-4f19-8a57-49c3938d607d.htm.(accessed at January 2014)

TIME SERIES MODELLING FOR WIND POWER PREDICTION AND CONTROL