arma_wind.m

For generating wind speed model in ARMA

t_sim = 300;    
t_sampling=100;   

% ARMA series, original equation, ref: IEEE paper
%  y(t)=1.7901*y(t-1)-0.9087*y(t-2)+0.0948*y(t-3)+alpha(t)-1.0929*alpha(t-1)+0.2892*alpha(t-2) ;
%%%%  
% y(1)=1.7901*y0 + alpha(1); 
% y(2)=1.7901*y(1) - 0.9087*y0  + alpha(2) - 1.0929*alpha(1);
% y(3)=1.7901*y(2) - 0.9087*y(1) + 0.0948*y0 + alpha(3) -1.0929*alpha(2) + 0.2892*alpha(1) ;
% %y(i)= (1.7901*y(i-1)-0.9087*y(i-2)+0.0948*y(i-3)+alpha(i)-1.0929*alpha(i-1)+0.2892*alpha(i-2)) ;

alpha=rand(1,t_sim * t_sampling) ; % random number matrix

%%%%%%%%%%%%% generation of initial values
y0 = 0.005;  %0.01  %    6.1,  let initial wind speed y0 = y(0)=5.75

%%%%%%%%%%%%  weighting values 
a=1.7901;
b=0.9087;
c=0.09848;
d=1.0929;
e=0.2892;

y(1)=a*y0 + alpha(1) ;
y(2)=a*y(1) - a*y0  + alpha(2) - d*alpha(1);
y(3)=a*y(2) - a*y(1) + c*y0 + alpha(3) -d*alpha(2) + e*alpha(1) ;

for i = 4 : t_sim*t_sampling    
y(i)= (a*y(i-1)-b*y(i-2)+c*y(i-3)+alpha(i)-d*alpha(i-1)+e*alpha(i-2)) ;
end

%%% calculation of wind speed from ARMA series
y_mean=mean(y);   % average wind speed
y_std=std(y);     % standard deviation
wind_speed = (y_mean + y_std * y)

plot(wind_speed)