风电负荷聚类——出力预测

clc clear all %% 导入数据 X=xlsread('风电'); x0=X; y=xlsread('负荷数据'); figure; subplot(2,1,1); plot(X);xlabel('时间/h');ylabel('标幺值');title('全年风电数据');axis([1 8761 0 1]) subplot(2,1,2);hold on plot(y);xlabel('时间/h');ylabel('功率/KW');title('全年负荷用电数据');axis([1 8761 400 1100]); %% 变量设定 k = 3; Eps = 2; data = [X,y]; %% 准备变量，输出原始结果 [m,n] = size(data); data=[(1:m)',data]; n = n + 1; type = zeros(1,m); cluster_No = 1; visited = zeros(m,1); class = zeros(1,m)-2; figure(2); subplot(2,1,1); plot(data(:,2),data(:,3),'k.'); grid on %daspect([1 1 1]); xlabel('x');ylabel('y'); title('原始输入点'); hold on; %% 运行DBSCAN算法 Kdtree = KDTreeSearcher(data(:,2:3)); for i = 1:m % 抽取一个未访问点 if visited(i)==0 % 标为访问 visited(i) = 1; point_now = data(i,:); Idx_range = rangesearch(Kdtree, point_now(2:3), Eps); index = Idx_range{1}; if length(index) > k class(i) = cluster_No; while index if visited(index(1)) == 0 visited(index(1)) = 1; if class(index(1)) <= 0 class(index(1)) = cluster_No; end point_now = data(index(1),:); Idx_range = rangesearch(Kdtree, point_now(2:3), Eps); index_temp = Idx_range{1}; index(1) = []; if length(index_temp) > k index = [index, index_temp]; end else index(1) = []; end end cluster_No = cluster_No + 1; end end end %% DBSCAN聚类结果 %figure; for i = 1: cluster_No color = [rand(),rand(),rand()]; data_class = data(find(class==i),:); plot(data_class(:,2),data_class(:,3),'.','Color',color,'MarkerFaceColor',color); hold on end data_class = data(find(class<=0),:); subplot(2,1,2); plot(data_class(:,2),data_class(:,3),'r*'); hold on plot(data(:,2),data(:,3),'.'); grid on xlabel('x');ylabel('y'); title('DBscan聚类结果'); %% 数据预处理 去除噪声数据 X(data_class(:,1)')=[]; y(data_class(:,1)')=[]; %% 电负荷有序聚类 M=floor(sqrt(length(y)));%确定分割点个数最大取值 [P1,D1]=Fisher(y,92);%有序聚类分割 O11=sort(P1); P2=zeros(92); P1=P1';C0=0;y=y'; for ll=2:M %得出不同分割点个数所对应的分割位置 P2(ll,:)=[sort(P1(:,1:ll)),zeros(1,92-ll)]; end for l=2:15 %参照论文讨论到分割15组，出于程序运行时间考虑，后续自行尝试（92） P3=fenlie(l,P2,y);%根据分割位置对负荷数据进行分组 CCd(l,:)= xishu(l,P3);%计算不同分割个数的轮廓系数 end figure(3); plot(CCd*2);axis([1 15 0 0.7]);xlabel('聚类个数');ylabel('轮廓系数'); [C0,mm]=max(CCd);%确定最佳分组数 PP3=fenlie(mm,P2,y); for qi=1:mm+1 qiwang(qi,:)=sum(PP3(qi,:))./length(find(PP3(qi,:)~=0));%期望值作为典型负荷场景 gehsu(qi,:)=length(find(PP3(qi,:)~=0)); pyfu(qi,:)=qiwang(qi,:)-50;pyzheng(qi,:)=qiwang(qi,:)+50; end gehsu(2,:)=gehsu(2,:)-5165; py1=[pyfu(1,:)*ones(1,gehsu(1,:))-50 pyfu(2,:)*ones(1,gehsu(2,:)) pyfu(3,:)*ones(1,gehsu(3,:))+60 pyfu(4,:)*ones(1,gehsu(4,:))+65 pyfu(5,:)*ones(1,gehsu(5,:))... pyfu(6,:)*ones(1,gehsu(6,:)) pyfu(7,:)*ones(1,gehsu(7,:))-60]; py2=[pyzheng(1,:)*ones(1,gehsu(1,:))-30 pyzheng(2,:)*ones(1,gehsu(2,:)) pyzheng(3,:)*ones(1,gehsu(3,:))+80 pyzheng(4,:)*ones(1,gehsu(4,:))+110 pyzheng(5,:)*ones(1,gehsu(5,:))... pyzheng(6,:)*ones(1,gehsu(6,:)) pyzheng(7,:)*ones(1,gehsu(7,:))-40]; figure(4); pic01 = fill([tt,fliplr(tt)],[py1,fliplr(py2)],'y');axis([1 8726 400 1100]) hold on plot(y);xlabel('时间/h');ylabel('功率/KW');title('全年负荷用电数据');axis([1 8761 400 1100]); %% 风电K-means聚类 X=X'; PW3=fenlie(mm,P2,X);%依据负荷分组对风电出力进行分组 Pw=PW3(2,:); geh=length(find(PW3(2,:)~=0))-5165; PWb=zeros(58,24);%这个相当于你申请一个空间，即m*n的数组 %%%%选取第二组场景进行聚类%%%%%%%%%%%% for i=1:58%m行 for j=1:24%n列 PWb(i,:)=Pw((1:24)+(i-1)*24);%取每n个数作为一行存入b数组 end end [IDW,CW,sumdw,DW] = kmeans(PWb,6); clustw1=find(IDW==1); clustw1=PWb(clustw1,:); numw1=length(clustw1(:,1)); clustw2=find(IDW==2); clustw2=PWb(clustw2,:); numw2=length(clustw2(:,1)); clustw3=find(IDW==3); clustw3=PWb(clustw3,:); numw3=length(clustw3(:,1)); clustw4=find(IDW==4); clustw4=PWb(clustw4,:); numw4=length(clustw4(:,1)); clustw5=find(IDW==5); clustw5=PWb(clustw5,:); numw5=length(clustw5(:,1)); clustw6=find(IDW==6); clustw6=PWb(clustw6,:); numw6=length(clustw6(:,1)); figure(5); subplot(3,2,1); plot((clustw1')/(max(max(clustw1))),'-','LineWidth',1.5);xlabel('时间/h');ylabel('标幺值');title('场景聚类1');axis([1 24 0 1]) subplot(3,2,2); plot((clustw2')/(max(max(clustw2))),'-','LineWidth',1.5);xlabel('时间/h');ylabel('标幺值');title('场景聚类2');axis([1 24 0 1]) subplot(3,2,3); plot((clustw3')/(max(max(clustw3))),'-','LineWidth',1.5);xlabel('时间/h');ylabel('标幺值');title('场景聚类3');axis([1 24 0 1]) subplot(3,2,4); plot((clustw4')/(max(max(clustw4))),'-','LineWidth',1.5);xlabel('时间/h');ylabel('标幺值');title('场景聚类4');axis([1 24 0 1]) subplot(3,2,5); plot((clustw5')/(max(max(clustw5))),'-','LineWidth',1.5);xlabel('时间/h');ylabel('标幺值');title('场景聚类5');axis([1 24 0 1]) subplot(3,2,6); plot((clustw6')/(max(max(clustw6))),'-','LineWidth',1.6);xlabel('时间/h');ylabel('标幺值');title('场景聚类6');axis([1 24 0 1])