aryc3.m

开发的指数平滑发用于短期负荷预测

clear
clc
load 2005gxxdata   %原数据构成矩阵a
% 对原始数据进行有效性处理
[i j]=size(a);
for m=1:i
    for n=2:j
        if ((a(m,n)-a(m,n-1))>0)&(a(m,n)/a(m,n-1)>1.1)
    a(m,n)=a(m,n-1)*1.1;
elseif ((a(m,n)-a(m,n-1))<0)&((a(m,n)/a(m,n-1))<0.9)
    a(m,n)=a(m,n-1)*0.9;
end
    end
end


%平稳化数据 
x=zeros(24,j);     %初始化x矩阵
for m=1:24
    for n=9:j
        x(m,n)=a(m,n)-a(m,n-1)+a(m,n-8)-a(m,n-7);   %序列平稳化  xt=(1-B)(1-B^7)yt
    end
end
%对预测日每一点进行自动搜索定阶
%
 date=input('date=');

if date<=23
    'error:please input date>=24 '
else
for dot=1:24
    for p=5:10
        for m=1:p
            A(:,m)=(x(dot,date-m-1:-1:date-m-14))';
        end
        b=x(dot,date-1:-1:date-14)';
      g=A\b;
      for j=1:p
          xx(j)=x(dot,date-j);
      end
      xxt(p)=xx*g;
        bt=A*g;
        da2(p-4)=sum((b-bt).^2);
    end
    e=min(da2);
    pp=find(da2==e)+4;
    if length(pp)>=2
        pp=pp(1);
        
    end
    ppp(dot)=pp;
        
    switch pp
        case 1 
            mt(dot)=xxt(1);
        case 2
            mt(dot)=xxt(2);
        case 3
            mt(dot)=xxt(3);
        case 4
             mt(dot)=xxt(4);
        case 5
             mt(dot)=xxt(5);
        case 6
             mt(dot)=xxt(6);
        case 7
             mt(dot)=xxt(7);
        case 8
             mt(dot)=xxt(8);
        case 9
             mt(dot)=xxt(9);
        case 10
             mt(dot)=xxt(10);
    
    end
end
xt=mt';
  yc=xt+a(:,date-1)+a(:,date-7)-a(:,date-8);
  
  %对每一预测点的值与前一时刻比较，误差超过10%的进行修正
for dot=1:24
if ((yc(dot)-a(dot,date-1))>0)&(yc(dot)/a(dot,date-1)>1.02)
    yc(dot)=a(dot,date-1)*1.02;
elseif ((yc(dot)-a(dot,date-1))<0)&(abs(yc(dot)/a(dot,date-1))<0.98)
    yc(dot)=a(dot,date-1)*0.98;
end
end
yc;
% 计算精度 
error=a(:,date)-yc;
re=error./a(:,date);
ed=sqrt(sum(abs(re).^2)/24);
Ad=1-ed
% 对实际值与预测值进行曲线拟合 
z=[1:24];
plot(z,a(:,date),'-o',z,yc,'-v');
end