基于特征选取和相关相量机的负荷预测代码

function [ output_args ] = RVM_Load_Predicting( input_args ) %RVM_LOAD_PREDICTING Summary of this function goes here % Detailed explanation goes here clear clc tic %% 加载数据 load WeatherData; load powerLoadData; load timeOfLoad; weatherData = data; powerLoadData = powerLoad; timeOfload = textdata; %每月的天数 daysPerMonth = [31,28,31,30,31,30,31,31,30,31,30,31]; temperatureSruct = struct('minTemp',0,'maxTemp',0); totalDays = length(powerLoadData)/24; loadDataMatrix = zeros(totalDays,24); weatherMaxtrix = zeros(totalDays,2); %% 特征选取 % 相空间重构 SN=1200; CN=168; Set1=zeros(SN,CN); for i=1:SN t=1; for j=i:CN+i-1 Set1(i,t)=powerLoadData(j); t=t+1; end Target(i)=powerLoadData(CN+i); end % 归一化 到[0,1]期间 n_Set1 = mapminmax(Set1',1,2); n_Set = n_Set1'; [n_Target, PS]=mapminmax(Target,1,2); %归一化完毕 % 找各样本的前 k 位nearest Hit 和 nearest Miss CI1=find(n_Target<=median(n_Target)); %小于中位数以下的索引下标 CI2=find(n_Target>median(n_Target)); %大于中位数以上的索引下标 C_Set1=n_Set(CI1,:); % 根据 CI 的值将样本分成两类 C_Set2=n_Set(CI2,:); k=round(log2(SN)); % k为 要找的同类或非同类邻近点个数(nearest Hit or nearest Miss) % 在两类样本中与所选样本有最小欧氏距离的样本为所选样本的 NH 或NM NH=zeros(SN,CN,k); % 用来存储前K位邻近点的矩阵 NM=zeros(SN,CN,k); for i=1:SN if ~isempty(find(CI1==i)) % 如果 i属于CI1类，则在C_Set1类找其nearest Hit, % 在C_Set2类找其nearest Miss. HitDist=[]; for j=1:length(CI1) a=n_Set(i,:); b=C_Set1(j,:); HitDist(j)=(sum(a.^2)+sum(b.^2)-2*sum(a.*b))^0.5; end MissDist=[]; for j=1:length(CI2) a=n_Set(i,:); b=C_Set2(j,:); MissDist(j)=(sum(a.^2)+sum(b.^2)-2*sum(a.*b))^0.5; end [newHitDist,t1]=sort(HitDist); [newMissDist,t2]=sort(MissDist); indHit=CI1(t1(2:k+1)); indMiss=CI2(t2(1:k)); else % 否则在C_Set2类找其nearest Hit, % 在C_Set1类找其nearest Miss. HitDist=[]; for j=1:length(CI2) a=n_Set(i,:); b=C_Set2(j,:); HitDist(j)=(sum(a.^2)+sum(b.^2)-2*sum(a.*b))^0.5; end MissDist=[]; for j=1:length(CI1) a=n_Set(i,:); b=C_Set1(j,:); MissDist(j)=(sum(a.^2)+sum(b.^2)-2*sum(a.*b))^0.5; end [newHitDist,t1]=sort(HitDist); [newMissDist,t2]=sort(MissDist); indHit=CI2(t1(2:k+1)); indMiss=CI1(t2(1:k)); end % 找出了每个样本的前10位 nearest Hits和nearest Misses for ii=1:k NH(i,:,ii)=n_Set(indHit(ii),:); NM(i,:,ii)=n_Set(indMiss(ii),:); end if(mod(i,50)==0) % 隔20次输出一次i i end end % 计算各个特征输出的权重 % 形成差矩阵 DNH=zeros(SN,CN,k); % 差矩阵 DNM=zeros(SN,CN,k); % for i=1:k DNH(:,:,i)=n_Set-NH(:,:,i); DNM(:,:,i)=n_Set-NM(:,:,i); end % 计算各特征输出在前K个邻近点的权重影响 WNH=zeros(k,CN); WNM=zeros(k,CN); for i=1:k for j=1:CN WNH(i,j)=sum(abs(DNH(:,j,i)))/SN; WNM(i,j)=sum(abs(DNM(:,j,i)))/SN; end end % 计算各特征的综合权重 W=zeros(1,CN); for i=1:CN W(i)=sum(WNM(:,i))/sum(WNH(:,i)); end W=mapminmax(W,0,1); [newW,index]=sort(-W); featureID=index(1:30)'; fprintf('选取算法得到的输入特征量为\n'); for i=1:length(featureID) if featureID(i)<=CN temp=abs(featureID(i)-(CN+1)); fprintf('Load(t-%d)\n',temp); end end % 特征选取结束 fprintf('特征选取结束\n'); %% 踢除冗余特征 % figure(1) % hist(W); % 计算各特征与目标输出的相关系数 X=40; TH=0.92; for i=1:X for j=1:X rho(:,:,i)=corrcoef(n_Set(:,index(i)),n_Set(:,index(j))) correlations(i,j)=abs(rho(1,2,i)); end end correlations; ID=1; Rid=[]; t1=1; t2=2; for i=2:X if ~isempty(find(correlations(i,1:i-1)>TH)) Rid(t1)=i; t1=t1+1; else ID(t2)=i; t2=t2+1; end end Fid=index(ID); % Fid = [168 145 121 97 25 1]; %% 负荷预测 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%生成数据训练集%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% fprintf('电价预测开始\n'); IN=length(Fid); %输入变量的维数 % Fid=newID(1:IN); % feature id 该数组存储特征输入的下标 fprintf('最终输入特征量为\n'); for i=1:IN if Fid(i)<=CN temp=abs(Fid(i)-(CN+1)); fprintf('Load(t-%d)\n',temp); end end NumOfSample = 480; %为训练样本 point = 24; % 预测点数 for i=1:NumOfSample input_train(i,:)=n_Set(i,Fid); end output_train=n_Target(1:NumOfSample)'; test_input = n_Set(NumOfSample + 1 : NumOfSample + 24, Fid); actual_value = Target(NumOfSample+1:NumOfSample+point); %% %2015年整年每天的最高气温和最低气温 t = 0; for i = 1:length(daysPerMonth) for j = 1:daysPerMonth(i) t = t+1; max = weatherData(j,2*i); min = weatherData(j,2*i+1); temperatureSruct(t).minTemp = min; temperatureSruct(t).maxTemp = max; weatherMaxtrix(t,1) = max; weatherMaxtrix(t,2) = min; end end % 对天气数据归一化 [nor_weather, PS_weather] = mapminmax(weatherMaxtrix',1,2); % 加入天气参数 % 最高温度 Tmax, 最低温度Tmin % for i = 1:NumOfSample % integer = ceil(i/point); % input_train(i,length(Fid) + 1) = nor_weather(1,integer+7); % input_train(i,length(Fid) + 2) = nor_weather(2,integer+7); % end % % % 测试集样本加入天气参数 % for i = 1:point % integer = ceil(i/point); % test_input(i,length(Fid) + 1) = nor_weather(1,integer + 7 + NumOfSample/point); % test_input(i,length(Fid) + 2) = nor_weather(2,integer + 7 + NumOfSample/point); % end x = input_train; t = output_train; x_test = test_input; %% 训练样本并返回测试结果 testOut = SparseBayesTrainAndPredict(x,t,x_test); %% 结果处理与显示 predict_value = mapminmax('reverse',testOut',PS); predict_value = predict_value'; [actual_value',predict_value] RMS = sqrt(sum((actual_value' - predict_value).^2)) plot(actual_value,'-b') hold on plot(predict_value','-r') hold off rvm2_predict = predict_value; save RVM actual_value rvm2_predict toc end