时间序列法预测风速

matlab风速预测程序,基于时间序列的风速预测程序

function main
clear;clf;clc;
load data_wind wind;y=chafeng（wind2）;%x做一次差分得y
[cqcs]=fangcha（y）;%求y的自，偏自相关系数存于cqcs中并画他们的序列图
%选ARMA（4，0）模型；
%[p1p2]=ARMA（62y）;%求已知阶数ARMA模型的参数
[p1p2]=ARMA（60y）;%求已知阶数AR模型的参数
%[p1p2]=ARMA（01y）;%求已知阶数MA模型的参数
result=[];N=size（y2）;%p1=[-0.0172 0.2932 0.2161 0.1603];
n=size（p11）;
for k=1:n                              %%%%验证模型对y的符合度
    result（k）=y（k）; 
end;
for k=n+1:N
    sum1=0;
    for i=1:n
        sum1=sum1+p1（i）*y（k-i）;
    end;
    result（k）=sum1;
end;
for k=1:10                              %%%%对y的预测
    sum1=0;
    for i=1:n
        sum1=sum1+p1（i）*result（k+N-i）;
    end;
    result（k+N）=sum1;
 
end;
 % subplot（223）;plot（y）;hold on;
 
  plot（result‘r‘）;title（‘一次差分‘）;%%%画图
result1（1）=wind（1）;
for k=2:N               %%%%回代得到x的预测，这是一次差分的逆过程
    result1（k）=wind（k-1）+result（k-1）;
end;
%for k=N+1:N+10
 for k=2:N
   result1（k）=result1（k-1）+result（k-1）   
end;
%subplot（224）;
plot（result1‘*-‘）;
hold on;
plot（wind‘r+-‘）;
grid on
xlabel（‘Time（h）‘）
ylabel（‘Vw（m/s）‘）
legend（‘预测‘‘实际‘）;
%title（‘拟合后‘）;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [faitheta]=ARMA（pqwind）
%%函数[faitheta]=ARMA（pqwind）
%%%pq:ARMA模型的阶数，wind:原始序列值，
%%fai:对应p的参数，theta:对应q的参数

N=size（wind2）;
theta=[];
gama=[];G=[];g=[];
for k=0:20%%%%%%%%%%%%%%%%%%%%%%%%%求fai
    sum=0;
    for i=1:N-k
        sum=sum+wind（i）*wind（i+k）;
    end; 
    gama（k+1）=sum/N;
end;
for k=1:p
    g（k）=gama（q+k+1）;%%生成g向量
    for i=1:p
        G（ki）=gama（abs（q+k-i）+1）;%%%生成G矩阵
    end;
end;
fai=inv（G）*g‘;%%%%公式 5-83
%由于求MA模型参数时涉及到非线形方程组的求解，没有统一的公式，本程序只就q=123时给出程序
y=[];                 %%%求theta
k=0;i=0;
for k=p+1:N
    sum=0;
    for i=1:p
        sum=sum+fai（i）*wind（k-i）;
    end;
    y（k-p）=wind（k）-sum;
end;
k=0;i=0;
gama1=[];
for k=0:20
    sum=0;
    for i=1:N-k-p
        sum=sum+y（i）*y（i+k）;
    end;
    gama1（k+1）=sum/N;
end;
syms a1 a2 a3 g1 g2 g3 g4 delta;
switch q%%%%%分别根据q=123求非线性方程组   公式5-77
    case 1
        s=solve（‘g1=delta*（1+a1）‘‘g2=delta*（-a1）‘deltaa1）;
        theta（1）=subs（s.a1{g1gama1（1）}{g2gama1（2）}）;
    case 2
        s=solve（‘g1=delta*（1+a1+a2）‘‘g2=delta*（-a1+a1*a2）‘‘g3=delta*（-a2）‘deltaa1a2）;
        kk=subs（s.a1（1）{g1gama1（1）}{g2gama1（2）}{g3gama1（3）}）;
        kk=subs（s.a2（1）{g1gama1（1）}{g2gama1（2）}{g3gama1（3）}）;
    case 3
        s=solve（‘g1=delta*（1+a1+a2+a3）‘‘g2=delta*（-a1+a1*a2+a2*a3）‘‘g3=delta*（-a2+a2*a3）‘‘g4=delta*（-a3）‘deltaa1a2a3）;
        theta（1）=subs（s.a1{g1gama1（1）}{g2gama1（2）}{g3gama1（3）}{g4gama1（4）}）;
        theta（2）=subs（s.a2{g1gama1（1）}{g2gama1（2）}{g3gama1（3）}{g4gama1（4）}）;
        theta（3）=subs（s.a3{g1gama1（1）}{g2gama1（2）}{g3gama1（3）}{g4gama1（4）}）;
end;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function y=chafeng（windN）%%函数差分
%%%%x:原始序列，N：差分的次数
%%%y:差分后的序列
p=[];
i=size（wind2）;
for k=N+1:i
    p（k-N）=wind（k）-wind（k-N）;
end;
y=p;
%%%%%%%%

 属性            大小     日期    时间   名称
----------- ---------  ---------- -----  ----

     文件       4311  2009-05-05 03:20  time_series.m

----------- ---------  ---------- -----  ----

                 4311                    1