Etiquetas
If you need to generate synthetic wind speed time series, you may find useful the procedure described in «A Markov method for simulating non-gaussian wind speed time series» by G.M. McNerney and P.S. Veers (Sandia Laboratories, 1985), and «Estimation of extreme wind speeds with very long return periods» by M.D.G Dukes and J.P. Palutikof (Journal of applied meteorology, 1994). This procedure is internally implemented in the software Hybrid2 and Homer. I have implemented this procedure with R. The script is available here and is described below:
First, we have to define the parameters of the wind speed distribution:
MeanSpeed<-6 MaxSpeed<-30 nStates<-30; nRows<-nStates; nColumns<-nStates; ##Position of each state LCateg<-MaxSpeed/nStates; ##A Rayleigh fdp is a weibull with shape factor 2 Shape=2; ##and scale factor sigma*sqrt(2), Scale=2*MeanSpeed/sqrt(pi);
Now, we build a vector of wind speeds centered around the mean value of each state and obtain the density function defined by Shape and Scale:
##Vector of wind speeds WindSpeed=seq(LCateg/2, MaxSpeed-LCateg/2, by=LCateg); ##FDP of wind speeds (P matrix) fdpWind<-dweibull(WindSpeed,shape=Shape, scale=Scale); fdpWind<-fdpWind/sum(fdpWind);
The correlation between the different categories of wind speed is defined by a matrix constructed with a decreasing function:
##Decreasing function g<-function(x){2^(-abs(x))} G<-matrix(nrow=nRows,ncol=nColumns) G <- row(G)-col(G) G <- g(G)
Then, an iterative procedure is needed for the calculation of the matrix of initial probability (P matrix):
##Initial value of the P matrix P0<-diag(fdpWind); ##Initial value of the p vector p0<-fdpWind; #The iterative procedure should stop when reaching a predefined error. However, for simplicity I have only constructed a for loop. To be improved! steps=10; P=P0; p=p0; for (i in 1:steps){ r<-P%*%G%*%p; r<-as.vector(r/sum(r)) p=p+0.5*(p0-r)#vector P=diag(p)}#matrix
Let’s combine the p vector and P and G matrices for obtaining a Markov Transition Matrix:
N=diag(1/as.vector(p%*%G)); MTM=N%*%G%*%P MTMcum<-t(apply(MTM,1,cumsum));
With this cumulated MTM we are able to simulate a wind speed time series:
##One value per minute for 1 day LSerie<-24*60; #Random number for choosing the state RandNum1<-runif(LSerie); State<-InitialState<-1; StatesSeries=InitialState; ##The next iterative procedure chooses the next state whose random number is greater than the cumulated probability defined by the MTM for (i in 2:LSerie) { State=min(which(RandNum1[i]<=MTMcum[State,])); StatesSeries=c(StatesSeries,State)} ##Another random number to choose wind speeds inside each category RandNum2<-runif(LSerie); ##And the wind speed series is SpeedSeries=WindSpeed[StatesSeries]-0.5+RandNum2*LCateg; ##where the 0.5 correction is needed since the the WindSpeed vector is centered around the mean value of each category.
A result of this procedure is displayed in the next figure:
We can check that the distribution of this wind speed series is correct:
library(MASS) print(fitdistr(SpeedSeries, 'weibull'))
graeme dijo:
Hello, I found this code very useful, so I made a version in Python:
http://www.lutralutra.co.uk/2012/07/02/simulating-a-wind-speed-time-series-in-python/
Thanks for the information!
Finnix dijo:
Hi.
You mention that this procedure is implemented in the software Homer, but this one uses three extra parameters, apart from the Mean Speed and the Weibull shape, both included in your R script:
– autocorrelation: «a measure of how strongly the wind speed in one hour depends on the wind speeds in previous hours»
– diurnal pattern strength: «measure of how strongly the wind speed tends to depend on the time of day»
– hour of peak windspeed: «the hour of the day that tends to be the windiest, on average»
I’m very interested in modifying your script to include these three parameters, and wonder if there is a modified version of it.
Oscar Perpiñán Lamigueiro dijo:
Hi,
I wrote this provisional script for a data analysis. Finally I didn’t get involved in that analysis and couldn’t find time to improve this script.
Please, let me know if you modify it to include more parameters.
Thanks!
Oscar
Darío dijo:
Very interesting!! Do you think its posible to obtain the wind rose? and the distribution depending the elevation? I am thinking to include your script in our project to study the feasibility for an offshore wind power plant in different locations.
Oscar Perpiñán Lamigueiro dijo:
The wind speed distribution changes with the altitude. Therefore you can use this procedure tuning the values of MeanSpeed and MaxSpeed (and perhaps the distribution function…).
There are several functions in R to *plot* the wind rose. For example:
– windrose at oce
– rosavent at climatol
– windrose at circular
But first you have to generate the wind directions sequence. You can use the procedure I described provided you know some information about the probabilistic distribution of the wind directions. For wind speeds Rayleigh and Weibull functions are particularly recommendable, but I am not so sure if it’s correct for directions.
Mirco dijo:
Hi!
Did you implemented the wind direction so far?
Cheers,
Mirco
Oscar Perpiñán Lamigueiro dijo:
No, Mirco. I have not worked on this code anymore.
tagesgeldkonto vergleich dijo:
Hey very nice blog!! Man .. Beautiful .. Amazing .. I will bookmark your blog and take the feeds
also…