bp1.m

短期负荷预测通常是指24小时的日负荷预测和168小时的周负荷预测

%BP_NET
clear
tic, % 开始计时
a=load('D:\课设资料\数据\样本数据\星期一\x.txt')
y=load('D:\课设资料\数据\样本数据\星期一\t.txt')
a=a/max(a);
y=y/max(y);
 y=diag(y);
[nn,n]=size(a); 
p=15;%隐群12个
[qq,q]=size(y);
w=2*rand(n,p)-ones(n,p); 
xw=zeros(n,p);
v=2*rand(p,q)-ones(p,q);
xv=zeros(p,q);
sita=rand(1,p);
xsita=zeros(1,p);
gama=rand(1,q);
xgama=zeros(1,q);
afa=0.7;  bita=0.7;
fai=0.8;
c=zeros(qq,q);
nx=0;
for N=1:1000 
    kn=N       %迭代次数
   allsee(N)=0;
   sallsee=0;
   for ii=1:nn
      a1=a(ii,:);
      see=0;  
            for j=1:p
               s(j)=0;b(j)=0;
            end
            for j=1:p
               for i=1:n
                s1=0;
                s1=w(i,j)*a1(i);
                 s(j)=s(j)+s1;
               end
               s(j)=s(j)-sita(1,j);
            x=exp(-s(j));
            b(j)=1/(1+x);
            end
   %the second step
            for i=1:q
               st(i)=0;
               c(ii,i)=0;
            end
            for t=1:q
               for j=1:p
                  st1=0;
                  st1=v(j,t)*b(j);
                  st(t)=st(t)+st1;
               end
               st(t)=st(t)-gama(1,t);
               x=exp(-st(t));
               c(ii,t)=1/(1+x);
            end
            
            for t=1:q
               d(t)=(y(ii,t)-c(ii,t))*c(ii,t)*(1-c(ii,t));
               for j=1:p
                  yitav=afa*d(t)*b(j)+fai*xv(j,t);
                  v(j,t)=v(j,t)+yitav;
                  xv(j,t)=yitav;
               end
               yitagama=-afa*d(t)+fai*xgama(1,t);
               gama(1,t)=gama(1,t)+yitagama;
               xgama(1,t)=yitagama;
            end
            for j=1:p
               e(j)=0;
               for t=1:q
                  x=d(t)*v(j,t);
                  e(j)=e(j)+x;
               end
               e(j)=e(j)*b(j)*(1-b(j));
               for i=1:n
                  yitaw=bita*e(j)*a1(i)+fai*xw(i,j);
                  w(i,j)=w(i,j)+yitaw;
                  xw(i,j)=yitaw;
               end
               yitasita=-bita*e(j)+fai*xsita(1,j);
               sita(1,j)=sita(1,j)+yitasita;
               xsita(1,j)=yitasita;
            end
           
            for t=1:q
               x=(y(ii,t)-c(ii,t))^2;
               see=see+x;
            end
            osee=see;
            sallsee=sallsee+osee;    
         end
         allsee(N)=sallsee/2;
      end
      c,y
      plot(allsee)
      for i = 1:size(c,1)
      cc(i)=c(i,i);
  end
      for i = 1:size(y,1)
      yy(i)=y(i,i);
  end
      E=sum((cc-yy).^2)./10;
      E=sqrt(E)
      e1=(yy-cc)./yy
      toc  % 结束计时
%##########################################
  xk=input('     please input .....xk==')
   xk
   %jj=input('     please input .....jj=')
   %ij=jj
  % a1=a(3,:)
   %a1(3,ij)=a1(3,ij)+xk
         for j=1:p
            s(j)=0;
            b(j)=0;
         end
            for j=1:p
               for i=1:n
                s1=0;
                s1=w(i,j)*a1(i);
                s(j)=s(j)+s1;
               end
               s(j)=s(j)-sita(1,j);
            x=exp(-s(j));
            b(j)=1/(1+x);
            end
            %the second step
            for i=1:q
               st(i)=0;
               cjd(i)=0;
            end
            for t=1:q
               for j=1:p
                  st1=0;
                  st1=v(j,t)*b(j);
                  st(t)=st(t)+st1;
               end
               st(t)=st(t)-gama(1,t);
               x=exp(-st(t));
               cjd(t)=1/(1+x);
            end
        FD=cjd