clear all
clc;
trainnum=12;
%訓練數據（輸入）
%用的是2008年1月份15，16，17號的負荷值和電價          po用來訓練的負荷值   p1用來預測負荷值  t0用來訓練的電價
p0=[29969 	29183 	28920 	29087 	29946 	32376 	36616 	38698 	38189 	37691 	37384 	37183 	37195 	37286 	37145 	37248 	38457 	41535 	41886 	41363 	40354 	38514 	35758 	35758 ;31376 	30557 	30196 	30170 	30932 	33330 	37361 	39319 	38756 	38173 	37766 	37250 	36674 	36346 	35897 	35816 	36945 	40256 	40836 	40419 	39624 	38023 	35455 	35455 ;31156 	30436 	30185 	30120 	30894 	33059 	37053 	39057 	38778 	38587 	38695 	38666 	38580 	38497 	38399 	38660 	39691 	41982 	42083 	41344 	40081 	37971 	34891 	34891 ];
p10=[31376 	30557 	30196 	30170 	30932 	33330 	37361 	39319 	38756 	38173 	37766 	37250 	36674 	36346 	35897 	35816 	36945 	40256 	40836 	40419 	39624 	38023 	35455 	35455;31156 	30436 	30185 	30120 	30894 	33059 	37053 	39057 	38778 	38587 	38695 	38666 	38580 	38497 	38399 	38660 	39691 	41982 	42083 	41344 	40081 	37971 	34891 	34891;31829 	30852 	30389 	30233 	30561 	31439 	32843 	34032 	35388 	35979 	35687 	34952 	34189 	33544 	33231 	33285 	34818 	37048 	37096 	36498 	35935 	34715 	32926 	32926 ];
t0=[46.57	42.06	40.46	41.11	46.09	66.76	100.91	88.73	81.57	76.09	79.45	71.05	63.88	60.27	56.39	52.78	66.28	92.62	96.5	83.21	75.23	64.56	48.5	44.32];
for tt=1:2;            %分兩次訓練

for i=1:3;             %每次分別取6*3個負荷p0，再取6個電價t0數據
    k=1;
    for j=(1+(tt-1)*trainnum):(trainnum*tt);
     p1(i,k)=p0(i,j);
     p2(i,k)=p10(i,j);
     t1(k)=t0(j);
     k=k+1;
    end;
end;

uu=max(max(p1));
p=p1./uu; %歸一化處理
%訓練數據（輸出）

uu=max(t1);
t=t1./uu;%歸一化處理
%%%----------------------------訓練部分----------------------------------------
s1=5;%隱層的神經元數
[r,q]=size(p);[s2,q]=size(t);
w1=rand(s1,r);w2=rand(s2,s1);%權值矩陣W1的初始化； %權值矩陣W2的初始化
b1=zeros(s1,1);b2=zeros(s2,1);%b1,b2的初始化
max_epoch=10000;lr=0.01;err_goal=0.005;%訓練的最大輪數； %學習因子； %誤差允許
dw1=0;db1=0;dw2=0;db2=0;
for epoch=1:max_epoch
    A1=tansig(w1*p,b1);A2=purelin(w2*A1,b2);%隱層實際輸出； %輸出層實際輸出
    E=t-A2;
    d2=deltalin(A2,E);%分別計算實際誤差
    d1=deltatan(A1,d2,w2);
    [dw2,db2]=learnbpm(A1,d2,lr,0.5,dw2,db2);%得到權值增量
    [dw1,db1]=learnbpm(p,d1,lr,0.5,dw1,db1);
   w2=w2+dw2;b2=b2+db2;%修正權值和閾值
   w1=w1+dw1;b1=b1+db1;
    c=purelin(w2*tansig(w1*p,b1),b2);%輸出層實際輸出
    sse=sumsqr(t-c);%檢驗誤差是否滿足要求
    if sse<err_goal
        break;end
end
%-------------------------------訓練結束-------------------------------------
%output 數組存放最后的預測結果
p3=p2./max(max(p2));%歸一化處理
A1=tansig(w1*p3,b1);
A3=purelin(w2*A1,b2);%計算輸出（歸一化的值）
A3=A3*uu;
for i=1:trainnum;
    output(i+(tt-1)*trainnum)=A3(i);
end;
end;      %對應的是第9行的for

%----------------------------------------一以下是畫圖--------------------------
zs=[42.41	38.24	37.21	37.33	38.72	48.12	82.43	88.96	81.86	78.01	69.64	62.71	57.36	54.92	48.57	44.72	54.52	87.36	95.41	81.35	73.07	59.23	44.84	44.95];
x=1:1:24;
plot(x,output,'g:',x,zs,'b');
title('比較,綠色的是預測值,藍色的是歷史值');
xlabel('time');
ylabel('price');