clc

clear

year=50 %模拟的年限

for l=1:year;

%%%%%%%%%%%%%%%%%%%%%%%%% 第一步：数据导入与预处理 %%%%%%%%%%%%%%%%%%%

%SW0=load('windspeed.txt'); %载入原始风速数据

SW0=xlsread('windspeed.xls');

SW0=SW0';

SW0=SW0/10*3.6; %原始数据的风速单位为0.1m/s，这里转化为km/h

N=size(SW0,2);

mu=mean(SW0);

sigma=var(SW0);

sigma=sigma^0.5; %求样本的平均值和标准差

y=(SW0-mu)./sigma; %数据预处理

% figure(1);

% subplot(211);

% autocorr(y); %画出自相关图

% title('自相关图');

% subplot(212);

% parcorr(y); %画出偏自相关图

% title('偏相关图');

%%%%%%%%%%%%%% 第二步：根据AIC准则确定ARMA模型的阶数 %%%%%%%%%%%%%%%%%%

for n=2:7;

m=armax(y',[n,n-1]);

fai=-m.a;

theta=m.c; %把armax函数得到的参数，取出来

for i=1:n;

y1(i)=y(i);

end

Noise=m.NoiseVariance^0.5;

e=normrnd(0,Noise,1,N);

for t=n+1:1:N;

y1(t)=0;

for j=2:n+1;

y1(t)=y1(t)+fai(j)*y1(t-(j-1));

end

for k=1:n;

y1(t)=y1(t)+theta(k)*e(t-(k-1));

end

end %y1(t)为预测值

s(n)=0;

for i1=1:N;

residual=y1(i1)-y(i1);

s(n)=s(n)+residual^2; %求取残差平方和

end

AIC(n)=N*log(s(n))+2*n-1; %求AIC

end

arma=AIC(1,2:7); %n=2-7时，各AIC的值

[AIC,n1]=min(arma);

n1=n1+1; %n1为得到的ARMA模型的阶数

%%%%%%%%%%%%%%%%% 第三步：用ARMA模型预测风速并确定风速分布 %%%%%%%%%%%%%%%%%

m=armax(y',[n1,n1-1]);

fai=-m.a;

theta=m.c;

for i=1:n1;

y2(i)=y(i);

end

Noise=m.NoiseVariance^0.5;

e=normrnd(0,Noise,1,8736);

for t=n1+1:1:8736;

y2(t)=0;

for j=2:n1+1;

y2(t)=y2(t)+fai(j)*y2(t-(j-1));

end

for k=1:n1;

y2(t)=y2(t)+theta(k)*e(t-(k-1));

end

end %y1(t)为预测值

y2;

SW=y2.*sigma+mu;

% figure(2)

% subplot(121);

% hist(SW0,100);

% xlabel('风速');ylabel('频数');title('原始风速的分布');

% subplot(122);

% hist(SW,100);

% xlabel('风速');ylabel('频数');title('预测风速的分布');

%得到8736个小时的预测风速：SW 1*8736

%%%%%%% 第四步：风电场的转移过程，确定其一年中三种状态分别存在的时长%%%%%%%

%风电机组的三状态模型

lambdaRD=5.84;lambdaRF=7.96;lambdaDR=48.3;lambdaFR=58.4;lambdaDF=0;lambdaFD=0; %风力发电机的3个状态的转移率

T=[1 1 1;lambdaRD -lambdaDR-lambdaDF lambdaFD;