风速威布尔分布和ARMA预测模型matlab程序

clc clear

%% 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.5

prob1=y/8760/0.5;%计算原始数据概率密度，频数除以数据种数，除以组距 prob2=(k/c)*(x/c).^(k-1).*exp(-(x/c).^k);%威布尔分布

figure(1)

title('Weibull分布拟合图');

bar(x,prob1,1) hold on

plot(x,prob2,'r')

legend('历史数据','Weibull拟合结果') % legend('Weibull拟合结果') hold off

save('result_weibull.mat')

%% 2.ARMA模型预测风速 clc clear load data

y=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 ~H

SaveDiffData=[SaveDiffData,TempData(1,1)];

TempData=diff(TempData);%差分，平稳化时间序列 difftime=difftime+1;%差分次数

H=adftest(TempData);-f检验，判断时间序列是否平稳化 end

%---------模型定阶或识别-------------- u = iddata(TempData); test = [];

for p = 1:5 %自回归对应PACF,给定滞后长度上限p和q，一般取为T/10、ln(T)或T^(1/2),这里取T/10=12

for q = 1:5 %移动平均对应ACF m = armax(u,[p q]);

AIC = aic(m); %armax(p,q),计算AIC test = [test;p q AIC]; end end

for k = 1:size(test,1)

if test(k,3) == min(test(:,3)) %选择AIC值最小的模型 p_test = test(k,1); q_test = test(k,2); break; end end

%------1阶预测-----------------

TempData=[TempData;zeros(step,1)]; n=iddata(TempData);

%m = armax(u(1:ls),[p_testq_test]); %armax(p,q),[p_testq_test]对应AIC值最小，自动回归滑动平均模型

m = armax(u,[p_testq_test]);

% ------------------------------------------- P1=predict(m,n,1); PreR=P1.OutputData; PreR=PreR';

Noise.std=sqrt(m.NoiseVariance); e=normrnd(0,Noise.std,1,300);

fori=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)~=0

for index=size(SaveDiffData,2):-1:1

PreR=cumsum([SaveDiffData(index),PreR]); end end

%-------------------预测趋势并返回结果---------------- mp1=polyfit([1:size(TrendData',2)],TrendData',1); xt=[];

for j=1:step

xt=[xt,size(TrendData',2)+j]; end

TrendResult=polyval(mp1,xt);

PreData=TrendResult+PreR(size(SourceData',2)+1:size(PreR,2)); tempx=[TrendData',TrendResult]+PreR; % tempx为预测结果 plot(tempx,'r-.'); hold on

plot(Data,'b');

legend('ARMA拟合时序曲线','实际时序风速'); save('resultarma.mat');