clcclear% 1.计算风速weibull分布% 数据处理load data;mu=mean(speed);%原始数据的统计参数sigma=sqrt(var(speed);% 计算威布尔分布参数parmhat=wblfit(speed);k=parmhat(2);c=parmhat(1);% k=(sigma/mu)-1.086;% c=mu/gamma(1+1/k);% 威布尔分布拟合y,x=hist(speed,ceil(max(speed)/0.5);%x是区间中心数,组距-1.5prob1=y/8760/0.5;%计算原始数据概率密度 ，频数除以数据种数，除以组距prob2=(k/c)*(x/c).(k-1).*exp(-(x/c).k);%威布尔分布figure(1)title(Weibull分布拟合图);bar(x,prob1,1) hold onplot(x,prob2,r)legend(历史数据,Weibull拟合结果)% legend(Weibull拟合结果)hold offsave(result_weibull.mat)% 2.ARMA模型预测风速clcclearload datay=speed(1:300);Data=y; %共300个数据SourceData=Data(1:250,1); %前250个训练集step=50; %后50个测试TempData=SourceData;TempData=detrend(TempData);%去趋势线TrendData=SourceData-TempData;%趋势函数%-差分，平稳化时间序列-H=adftest(TempData);difftime=0;SaveDiffData=;while HSaveDiffData=SaveDiffData,TempData(1,1);TempData=diff(TempData);%差分，平稳化时间序列difftime=difftime+1;%差分次数H=adftest(TempData);%adf检验，判断时间序列是否平稳化end%-模型定阶或识别-u = iddata(TempData);test = ;for p = 1:5 %自回归对应PACF,给定滞后长度上限p和q，一般取为T/10、ln(T)或T(1/2),这里取T/10=12for q = 1:5 %移动平均对应ACFm = armax(u,p q); AIC = aic(m); %armax(p,q),计算AICtest = test;p q AIC;endendfor k = 1:size(test,1)if test(k,3) = min(test(:,3) %选择AIC值最小的模型p_test = test(k,1);q_test = test(k,2);break;endend%-1阶预测-TempData=TempData;zeros(step,1);n=iddata(TempData); %m = armax(u(1:ls),p_test q_test); %armax(p,q),p_test q_test对应AIC值最小，自动回归滑动平均模型m = armax(u,p_test q_test); % -P1=predict(m,n,1);PreR=P1.OutputData;PreR=PreR;Noise.std=sqrt(m.NoiseVariance);e=normrnd(0,Noise.std,1,300);for i=251:300 PreR(i)=-m.A(2:p_test+1)*PreR(i-1:-1:i-p_test)+m.C(1:q_test+1)*e(i:-1:i-q_test);end% -%-还原差分-if size(SaveDiffData,2)=0for index=size(SaveDiffData,2):-1:1PreR=cumsum(SaveDiffData(index),PreR);endend %-预测趋势并返回结果-mp1=polyfit(1:size(TrendData,2),TrendData,1);xt=;for j=1:stepxt=xt,size(TrendData,2)+j;endTrendResult=polyval(mp1,xt);PreData=TrendResult+PreR(size(SourceData,2)+1:size(PreR,2);tempx=TrendData,TrendResult+PreR; % tempx为预测结果plot(tempx,r-.);hold onplot(Data,b);legend(ARMA拟合时序曲线,实际时序风速);save(resultarma.mat);