Saisonale Variationen der Windgeschwindigkeits-Zeitreihe hinzufügen

Question

Es folgt eine R blog, die interessant und recht nützlich ist, um die Zeitreihe eines unbekannten Gebiets mit seinen Weibull-Parametern zu simulieren.

Obwohl diese Methode eine recht gute Schätzung der Zeitreihe als Ganzes gibt, leidet sie sehr unter saisonalen Änderungen. Um saisonalen Veränderungen Rechnung zu tragen, möchte ich saisonale Höchstwindgeschwindigkeiten einsetzen und die Zeitreihen-Synthese so durchführen, dass die jährliche Verteilung konstant bleibt, d. Form- und Skalierungsparameter (Jahreswerte).

Ich möchte saisonale maximale Windgeschwindigkeiten nach dem folgenden Code verwenden, indem ich 12 verschiedene maximale Windgeschwindigkeiten, eine für jeden Monat, verwende. Dies ermöglicht größere Windgeschwindigkeiten in bestimmten Monaten und niedriger in anderen und sollte die resultierenden Zeitreihen ausgleichen.

Der Code folgt wie folgt aus:

MeanSpeed<-7.29 ## Mean Yearly Wind Speed at the site. 

Shape=2; ## Input Shape parameter (yearly). 
Scale=8 ##Calculated Scale Parameter (yearly). 

MaxSpeed<-17 (##yearly) 
## $$$ 12 values of these wind speed one for each month to be used. The resultant time series should satisfy shape and scale parameters $$ ### 
nStates<-16 

nRows<-nStates; 
nColumns<-nStates; 


LCateg<-MaxSpeed/nStates; 


WindSpeed=seq(LCateg/2,MaxSpeed-LCateg/2,by=LCateg) ## Fine the velocity vector-centered on the average value of each category. 

##Determine Weibull Probability Distribution. 
wpdWind<-dweibull(WindSpeed,shape=Shape, scale=Scale); # Freqency distribution. 

plot(wpdWind,type = "b", ylab= "frequency", xlab = "Wind Speed") ##Plot weibull probability distribution. 

norm_wpdWind<-wpdWind/sum(wpdWind); ## Convert weibull/Gaussian distribution to normal distribution. 

## Correlation between states (Matrix G) 
g<-function(x){2^(-abs(x))} ## decreasing correlation function between states. 
G<-matrix(nrow=nRows,ncol=nColumns) 
G <- row(G)-col(G) 
G <- g(G) 

##-------------------------------------------------------- 


## iterative process to calculate the matrix P (initial probability) 
P0<-diag(norm_wpdWind); ## Initial value of the MATRIX P. 
P1<-norm_wpdWind; ## Initial value of the VECTOR p. 


## This iterative calculation must be done until a certain error is exceeded 
## Now, as something tentative, I set the number of iterations 

steps=1000; 
P=P0; 
p=P1; 

for (i in 1:steps){ 
    r<-P%*%G%*%p; 
    r<-as.vector(r/sum(r)); ## The above result is in matrix form. I change it to vector 
    p=p+0.5*(P1-r) 
    P=diag(p)} 

    ## $$ ----Markov Transition Matrix --- $$ ## 

N=diag(1/as.vector(p%*%G));## normalization matrix 

MTM=N%*%G%*%P ## Markov Transition Matrix 

MTMcum<-t(apply(MTM,1,cumsum));## From the MTM generated the accumulated 

##------------------------------------------- 
## Calculating the series from the MTMcum 

##Insert number of data sets. 
LSerie<-52560; Wind Speed every 10 minutes for a year. 

RandNum1<-runif(LSerie);## Random number to choose between states 
State<-InitialState<-1;## assumes that the initial state is 1 (this must be changed when concatenating days) 
StatesSeries=InitialState; 

## Initallise---- 

## The next state is selected to the one in which the random number exceeds the accumulated probability value 
##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) { 
    ## i has to start on 2 !! 
    State=min(which(RandNum1[i]<=MTMcum[State,])); 

    ## if (is.infinite (State)) {State = 1}; ## when the above condition is not met max -Inf 
    StatesSeries=c(StatesSeries,State)} 

RandNum2<-runif(LSerie); ## Random number to choose between speeds within a state 

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. 


print(fitdistr(SpeedSeries, 'weibull')) ##MLE fitting of SpeedSeries 

Kann jemand empfehlen, wo und welche Änderungen muss ich den Code machen?

Answer 1

Ich weiß nicht viel über Windgeschwindigkeit Zeitreihen zu erzeugen, aber vielleicht können diese Richtlinien helfen Ihnen, Ihre Lesbarkeit des Codes/Wiederverwertbarkeit zu verbessern:

# 1 Sie wahrscheinlich eine Funktion haben wollen, die einen Wind erzeugen Geschwindigkeitszeit Serie eine Reihe von Beobachtungen und eine saisonale maximale Windgeschwindigkeit gegeben. Also versuchen Sie zuerst Ihren Code in einem Block wie diese zu definieren:

wind_time_serie <- function(nobs, max_speed){ 
    #some code here 
} 

# 2 so tun, wenn es scheint, dass einige Teile des Codes sind nützlich, Windgeschwindigkeit Zeitreihen zu erzeugen, aber nicht über Windgeschwindigkeits-Zeitreihe, versuchen Sie, sie in Funktionen zu setzen (zB der Teil, den Sie berechnen norm_wpdWind, der Teil, den Sie berechnen MTMcum, ...).

# 3 Dann der Teil Ihres Codes am Anfang, wenn Ihre globale Variable definieren sollte verschwinden und Standard-Argumente in Funktionen werden.

# 4 Vermeiden Sie die Verwendung von Endzeilenkommentaren, wenn Ihre Zeile bereits lang ist, und löschen Sie die Endsemicolumns.

#This 
    State<-InitialState<-1;## assumes that the initial state is 1 (this must be changed when concatenating days) 

#Would become this: 
    #Assumes that the initial state is 1 (this must be changed when concatenating days) 
    State<-InitialState<-1 

Dann sollte Ihr Code mehr wiederverwendbar/lesbar sein von anderen Menschen. Sie haben ein Beispiel unten für die Richtlinien, die auf den rnorm-Teil angewendet wurden:

norm_distrib<-function(maxSpeed, states = 16, shape = 2, scale = 8){ 

    #Fine the velocity vector-centered on the average value of each category. 
    LCateg<-maxSpeed/states 
    WindSpeed=seq(LCateg/2,maxSpeed-LCateg/2,by=LCateg) 

    #Determine Weibull Probability Distribution. 
    wpdWind<-dweibull(WindSpeed,shape=shape, scale=scale) 

    #Convert weibull/Gaussian distribution to normal distribution. 
    return(wpdWind/sum(wpdWind)) 

    } 

    #Plot normal distribution with the max speed you want (e.g. 17) 
    plot(norm_distrib(17),type = "b", ylab= "frequency", xlab = "Wind Speed")