基于SVM电力系统短期负荷预测的其中一个例程

clc; clear C = 30; theta = 2;%C为最小二乘支持向量机的正则化参数，theta为高斯径向基的核函数参数，两个需要进行优化选择调试 NumOfPre = 1;%预测天数，在此预测本季度最后七天 %以负荷的加权为参考1 al = 2;%温度的加权 be = 1;%湿度的加权 th = 1;%星期的加权 Time = 48; Data = xlsread('input.xls');%此为从excel表格读数据的命令，表示将表格的数据读到Data数组中，省略表格中的第一行第一列文字部分 可输入你要预测的表格名称 Data = Data(2:end,:); [M,N] = size(Data);%计算读入数据的行和列 M行N列 Data1 = Data; for i = 1:2 maxData = max(Data(:,i)); minData = min(Data(:,i)); Data1(:,i) = (Data(:,i) - minData)/(maxData-minData);%对温度进行归一化处理 end for i = 4:N Data1(:,i) = log10(Data(:,i)) ;%对负荷进行对数处理 温度和负荷的预处理 可采用不同的方法 可不必拘泥 end Dim = M - 2 - NumOfPre;%训练样本数%拥有完整特征描述的数据共有M-2组！！ Input = zeros(M,8,Time);%预先分配处理后的输入向量空间%8是特征空间的维度！也是x的元素数！！ y = zeros(Dim,Time); for i = 3:M for j = 1:Time %%选取前一天温度、同一时刻的负荷，前两天的负荷，当天的温度作为输入特征 x = [Data1(i-1,1:3), Data1(i-1,j+3), Data1(i-2,j+3),Data1(i,1:3)];%x就是每一个负荷对应的特征向量！！！ x(1)=al*x(1); x(6)=al*x(6); x(2)=be*x(2); x(7)=be*x(7); x(3)=th*x(3); x(8)=th*x(8); Input(i-2,:,j) = x;%Input(i,：,k)的含义为：第i天第j个时刻的负荷所对应的特征向量 y(i-2,j) = Data1(i,j+3); end end Dist = zeros(Dim,Dim,Time);%预先分配距离空间 for i=1:Time for j=1:Dim for k=1:Dim Dist(j,k,i) = (Input(j,:,i) - Input(k,:,i))*(Input(j,:,i) - Input(k,:,i))'; end end end %Dist(i,j,k)含义为：第i天第k个时刻的负荷与第j天同一时刻的负荷之间的特征距离 Dist1 = exp(-Dist/(2*theta));%RBF K = zeros(1,Dim); % Pre = zeros(M-2,Time); for i=1:Time H = Dist1(:,:,i) + eye(Dim)/C;%最小二乘支持向量的H矩阵 f = -y(1:Dim,i); Aeq = ones(Dim,1)'; beq = [0]; option.MaxIter=1000; [a,fval]=quadprog(H,f,[],[],Aeq,beq);%,[],[],[],option); b = zeros(1,Dim); for j = 1:Dim b(1,j) = y(j,i) - a(j)/C - a'* Dist1(:,j,i);%求每个输入特征对应的b end b = sum(b)/Dim;%求平均b，消除误差 for j = Dim + 1:M-2 for k = 1:Dim K(1,k) = exp(-(Input(j,:,i) - Input(k,:,i))*(Input(j,:,i) - Input(k,:,i))'/(2*theta));%预测输入特征与训练特征的RBF距离 end Pre(j-Dim,i) = sum(a'*K') + b; %求解预测值 end end % Len = M - (Dim + 3) + 1;%预测的天数 取本季度最后Len天 Len = NumOfPre; Pre = 10.^Pre; for i = 1:Len figure plot(1:Time,Data(i+Dim+2,4:N),'r',1:Time,Pre(i,:),'k');%画出每一天的预测值和真实值 hold on scatter(1:Time,Data(i+Dim+2,4:N),'o') scatter(1:Time,Pre(i,:),'^') % axis([0 25 0 100])%坐标范围 hold off end Acu = (Pre - Data(Dim+3:M,4:N))./Data(Dim+3:M,4:N);%相对误差 save Acu.mat Acu s=0; for i=1:Time s=abs(Acu(1,i))+s; end acu=s/Time; save acu.mat acu; Result=[C,theta,acu]; disp(Result);