From ec8a4e8189cc3ea79a6eaf16a7a6368043aa54c0 Mon Sep 17 00:00:00 2001 From: "dugg@lab-desk" Date: Fri, 27 Mar 2015 10:52:23 +0800 Subject: [PATCH] =?UTF-8?q?1.=E5=8A=A0=E5=85=A5=E7=94=9F=E6=88=90Laplace?= =?UTF-8?q?=E5=88=86=E5=B8=83=E9=9A=8F=E6=9C=BA=E6=95=B0=E7=9A=84=E5=87=BD?= =?UTF-8?q?=E6=95=B0=202.=E6=8A=8A=E4=BB=A3=E7=A0=81=E8=B0=83=E5=9B=9E?= =?UTF-8?q?=E6=9C=80=E5=88=9D=E5=BD=A2=E5=BC=8F=E3=80=82?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Signed-off-by: dugg@lab-desk --- OPF.m | 677 ++++++++++++++++++++++++++++--------------------------- laprnd.m | 30 +++ subOPF.m | 4 +- 3 files changed, 373 insertions(+), 338 deletions(-) create mode 100644 laprnd.m diff --git a/OPF.m b/OPF.m index 9ad9af7..7505c94 100644 --- a/OPF.m +++ b/OPF.m @@ -29,12 +29,12 @@ for badDataNode=1:1 mVolt=rVolt.*(1+normrnd(0,sigma,length(rVolt),1))'; %加载保存的变量 - PD0=load('PD0'); - PD0=PD0.PD0; - QD0=load('QD0'); - QD0=QD0.QD0; - mVolt=load('mVolt'); - mVolt=mVolt.mVolt; +% PD0=load('PD0'); +% PD0=PD0.PD0; +% QD0=load('QD0'); +% QD0=QD0.QD0; +% mVolt=load('mVolt'); +% mVolt=mVolt.mVolt; % mVolt(3)=rVolt(3)*(1-sigma*6); %% 画Case A的图 @@ -43,7 +43,7 @@ for badDataNode=1:1 % badDataResult(I,badDataNode)=sum(Vbi); % badDataLocation(1:33,I)=Vbi; % badDataLocation(34,I)=sum(abs((rVolt-Volt)./rVolt./length(rVolt)))+sum(abs( (UAngel(2:33)-rUAngel(2:33))./rUAngel(2:33)./length(rUAngel(2:33)))); - break; +% break; if isConverge==0 continue; end @@ -62,7 +62,7 @@ for badDataNode=1:1 nodeMaxDVAngle_t([1:18,23:33])=0; nodeMaxDVolt(nodeMaxDVolt=500 + if loopN>=1 % nodeMaxDVolt(badDataNode)=maxDVolt; % nodeMaxDVAngle(badDataNode)=maxDVAngle; break; @@ -72,335 +72,340 @@ for badDataNode=1:1 end end -% end -subplot(4,1,1,'XTick',[1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33]); -CaseAREV=(Volt-rVolt);%Relative Error of Voltage in Case A -% CaseAREV=CaseAREV(2:end)*100; -%真实值的 -plot(1:length(CaseAREV),(CaseAREV),'k.:','Marker','diamond'); -%测量值的 -% plot(1:length(CaseAREV),abs((mVolt-rVolt)*100),'c.:','Marker','diamond'); -box off; -set(gca,'XTick',[1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33]) -xlabel('节点号'); -ylabel('误差%'); -subplot(4,1,2); -% CaseAREA=(UAngel-rUAngel)*100;%Relative Error of Angle in Case A -CaseAREA=(UAngel-rUAngel);%Relative Error of Angle in Case A -CaseAREA(1)=0; -%真实值的 -plot(1:length(CaseAREA),(CaseAREA),'k:','Marker','diamond'); -box off; -set(gca,'XTick',[1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33]) -xlabel('节点号'); -ylabel('误差%'); -subplot(4,1,3); -% CaseAREP=(PD-RealPD)./(RealPD+0.00001)*100;%Relative Error of PD in Case A -CaseAREP=(PD-RealPD)./RealPD*100;%Relative Error of PD in Case A -CaseAREP(1)=0; -%真实值的 -plot(1:length(CaseAREP),(CaseAREP),'k:','Marker','diamond'); -%测量值的 -% plot(1:length(CaseAREV),abs((PD0-RealPD)./(RealPD+0.00001)*100),'c.:','Marker','diamond'); -box off; -set(gca,'XTick',[1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33]) -xlabel('节点号'); -ylabel('误差%'); -subplot(4,1,4); -% CaseAREQ=(QD-RealQD)./(RealQD+0.00001)*100;%Relative Error of QD in Case A -CaseAREQ=(QD-RealQD)./RealQD*100;%Relative Error of QD in Case A -CaseAREQ(1)=0; -%真实值的 -plot(1:length(CaseAREQ),(CaseAREQ),'k:','Marker','diamond'); -%测量值的 -% plot(1:length(CaseAREV),abs((QD0-RealQD)./(RealQD+0.00001)*100),'c.:','Marker','diamond'); -box off; -set(gca,'XTick',[1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33]) -xlabel('节点号'); -ylabel('误差%'); -%计算Case A的误差 -CaseAE=sqrt((sum(CaseAREV.^2)+sum(CaseAREA.^2)+sum(CaseAREP.^2)+sum(CaseAREQ.^2))/132); -objA=full(sum(Vbi)+sum(PDbi)+sum(QDbi)); -notZeros=find(PD0~=0); -CaseA_SE=sum(((RealPD(notZeros)-PD(notZeros))./(RealPD(notZeros)*1)).^2)+sum(((RealQD(notZeros)-QD(notZeros))./(RealQD(notZeros)*1)).^2)+sum(((Volt-rVolt)./(rVolt*1)).^2); -CaseA_SE=(CaseA_SE/(length(notZeros)*2+length(Volt)))^.5; -sumCaseA_SE=sumCaseA_SE+CaseA_SE; -% arrayA(1:19,I)=Vbi; -% arrayA(21,I)=CaseAE*1000; -%% Case B -% [Busnum,Loadi,Volt,PD,QD,rVolt,UAngel,RealPD,RealQD,rUAngel,Vbi,PDbi,QDbi,plotGapB]=subOPF(setdiff(1:Busnum,[18,21,22,29]),PD0,QD0,mVolt,sigma);% -[Busnum,Loadi,Volt,PD,QD,rVolt,UAngel,RealPD,RealQD,rUAngel,Vbi,PDbi,QDbi,plotGapB]=subOPF(setdiff(1:Busnum,[2,3,5,20,24,27,28,10,11,12,13]),PD0,QD0,mVolt,sigma);% -subplot(4,1,1); -hold on; -% CaseBREV=(Volt-rVolt)*100;%Relative Error of Voltage in Case B -CaseBREV=(Volt-rVolt);%Relative Error of Voltage in Case B -plot(1:length(CaseBREV),(CaseBREV),'b.:','Marker','square'); -subplot(4,1,2); -hold on; -% CaseBREA=(UAngel-rUAngel)*100;%Relative Error of Angle in Case B -CaseBREA=(UAngel-rUAngel);%Relative Error of Angle in Case B -plot(1:length(CaseBREA),(CaseBREA),'b:','Marker','square'); -subplot(4,1,3); -hold on; -% CaseBREP=(PD-RealPD)./(RealPD+0.00001)*100;%Relative Error of PD in Case B -CaseBREP=(PD-RealPD)./RealPD*100;%Relative Error of PD in Case B -RealPD(1)=0; -plot(1:length(CaseBREP),(CaseBREP),'b:','Marker','square'); -subplot(4,1,4); -hold on; -% CaseBREQ=(QD-RealQD)./(RealQD+0.00001)*100;%Relative Error of QD in Case B -CaseBREQ=(QD-RealQD)./RealQD*100;%Relative Error of QD in Case B -CaseBREQ(1)=0; -plot(1:length(CaseBREQ),(CaseBREQ),'b:','Marker','square'); -CaseBE=sqrt((sum(CaseBREV.^2)+sum(CaseBREA.^2)+sum(CaseBREP.^2)+sum(CaseBREQ.^2))/132); -objB=full(sum(Vbi)+sum(PDbi)+sum(QDbi)); -noMeasurei=[2,3,5,20,24,27,28,10,11,12,13]; -Measurei=setdiff(2:33,[2,3,5,20,24,27,28,10,11,12,13]); -CaseB_SE=sum(((RealPD(Measurei)-PD(Measurei))./(RealPD(Measurei)*1)).^2)+sum(((RealQD(Measurei)-QD(Measurei))./(RealQD(Measurei)*1)).^2)+sum(((Volt(Measurei)-rVolt(Measurei))./(rVolt(Measurei)*1)).^2); -CaseB_SE=CaseB_SE+sum(((RealPD(noMeasurei)-PD(noMeasurei))./(RealPD(noMeasurei)*1)).^2)+sum(((RealQD(noMeasurei)-QD(noMeasurei))./(RealQD(noMeasurei)*1)).^2)+sum(((Volt(noMeasurei)-rVolt(noMeasurei))./(rVolt(noMeasurei)*1)).^2); -CaseB_SE=(CaseB_SE/(length(notZeros)+length(noMeasurei) +length(Volt)))^.5; -sumCaseB_SE=sumCaseB_SE+CaseB_SE; -%% Case C -[Busnum,Loadi,Volt,PD,QD,rVolt,UAngel,RealPD,RealQD,rUAngel,Vbi,PDbi,QDbi,plotGapC]=subOPF([1:33],PD0,QD0,mVolt,sigma);% -subplot(4,1,1); -hold on; -% CaseCREV=(Volt-rVolt)*100;%Relative Error of Voltage in Case C -CaseCREV=(Volt-rVolt); -plot(1:length(CaseCREV),(CaseCREV),'r.:','Marker','o'); -subplot(4,1,2); -hold on; -% CaseCREA=(UAngel-rUAngel)*100;%Relative Error of Angle in Case C -CaseCREA=(UAngel-rUAngel); -plot(1:length(CaseCREA),(CaseCREA),'r:','Marker','o'); -subplot(4,1,3); -hold on; -% CaseCREP=(PD-RealPD)./(RealPD+0.00001)*100;%Relative Error of PD in Case C -CaseCREP=(PD-RealPD)./RealPD*100;%Relative Error of PD in Case C -CaseCREP(1); -plot(1:length(CaseCREP),(CaseCREP),'r:','Marker','o'); -subplot(4,1,4); -hold on; -% CaseCREQ=(QD-RealQD)./(RealQD+0.00001)*100;%Relative Error of QD in Case C -CaseCREQ=(QD-RealQD)./RealQD*100;%Relative Error of QD in Case C -CaseCREQ(1)=0; -plot(1:length(CaseCREQ),(CaseCREQ),'r:','Marker','o'); -% 画legend -subplot(4,1,1); -% title('Voltage'); -ld=legend('算例A','算例B','算例C'); -set(ld,'Position',[0.847865087908145 0.786094477711244 0.0543595263724435 0.0605455755156354]); -subplot(4,1,2); -% title('Voltage Angle'); -ld=legend('算例A','算例B','算例C'); -set(ld,'Position',[0.847865087908145 0.586094477711244 0.0543595263724435 0.0605455755156354]); -subplot(4,1,3); -% title('Active load power'); -ld=legend('算例A','算例B','算例C'); -set(ld,'Position',[0.847865087908145 0.386094477711244 0.0543595263724435 0.0605455755156354]); -subplot(4,1,4); -% title('Reactive load power'); -ld=legend('算例A','算例B','算例C'); -set(ld,'Position',[0.847865087908145 0.186094477711244 0.0543595263724435 0.0605455755156354]); -CaseCE=sqrt((sum(CaseCREV.^2)+sum(CaseCREA.^2)+sum(CaseCREP.^2)+sum(CaseCREQ.^2))/132); -objC=full(sum(Vbi)+sum(PDbi)+sum(QDbi)); -CaseC_SE=sum(((RealPD(notZeros)-PD(notZeros))./(RealPD(notZeros)*1)).^2)+sum(((RealQD(notZeros)-QD(notZeros))./(RealQD(notZeros)*1)).^2)+sum(((Volt-rVolt)./(rVolt*1)).^2); -CaseC_SE=(CaseC_SE/(length(notZeros)*2+length(Volt)))^.5; -sumCaseC_SE=sumCaseC_SE+CaseC_SE; +%% PLOTING -% fprintf('目标函数值 %.2f\n',full(obj)); -fprintf('Case A Case B Case C 的误差\n') -fprintf('%f %f %f \n',CaseAE,CaseBE,CaseCE); -fprintf('三个Case目标值\n') -fprintf('%f\t%f\t%f \n',objA,objB,objC) -%% 画测量值 -% subplot(4,1,1); -% plot(1:Busnum,mVolt-rVolt,'k.:','Marker','pentagram') +% % end +% subplot(4,1,1,'XTick',[1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33]); +% CaseAREV=(Volt-rVolt);%Relative Error of Voltage in Case A +% % CaseAREV=CaseAREV(2:end)*100; +% %真实值的 +% plot(1:length(CaseAREV),(CaseAREV),'k.:','Marker','diamond'); +% %测量值的 +% % plot(1:length(CaseAREV),abs((mVolt-rVolt)*100),'c.:','Marker','diamond'); +% box off; +% set(gca,'XTick',[1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33]) +% xlabel('节点号'); +% ylabel('误差%'); +% subplot(4,1,2); +% % CaseAREA=(UAngel-rUAngel)*100;%Relative Error of Angle in Case A +% CaseAREA=(UAngel-rUAngel);%Relative Error of Angle in Case A +% CaseAREA(1)=0; +% %真实值的 +% plot(1:length(CaseAREA),(CaseAREA),'k:','Marker','diamond'); +% box off; +% set(gca,'XTick',[1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33]) +% xlabel('节点号'); +% ylabel('误差%'); % subplot(4,1,3); -% plot(1:Busnum,(PD0-RealPD)./(RealPD+0.00001),'k:','Marker','pentagram'); +% % CaseAREP=(PD-RealPD)./(RealPD+0.00001)*100;%Relative Error of PD in Case A +% CaseAREP=(PD-RealPD)./RealPD*100;%Relative Error of PD in Case A +% CaseAREP(1)=0; +% %真实值的 +% plot(1:length(CaseAREP),(CaseAREP),'k:','Marker','diamond'); +% %测量值的 +% % plot(1:length(CaseAREV),abs((PD0-RealPD)./(RealPD+0.00001)*100),'c.:','Marker','diamond'); +% box off; +% set(gca,'XTick',[1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33]) +% xlabel('节点号'); +% ylabel('误差%'); % subplot(4,1,4); -% plot(1:Busnum,(QD0-RealQD)./(RealQD+0.00001),'k:','Marker','pentagram'); -%% 直方图 -% 电压 -figure('Name','电压直方图') -split_number=20; -%Case A -subplot(1,3,1) -y=CaseAREV; -ymin=min(y); -ymax=max(y); -x=linspace(ymin,ymax,split_number); %将最大最小区间分成split_number个等分点(19等分),然后分别计算各个区间的个数 -yy=hist(y,x); %计算各个区间的个数 -% yy=yy/(sum(yy)*(ymax-ymin)/split_number); %计算各个区间的个数,除以总面积，总面积计算的方式为：所有小分割的面积的和即: -bar(x,yy) %画出概率密度分布图 -voltBarCaseAX=x; -voltBarCaseAY=yy; -xlabel('Error'); -ylabel('Number of buses'); -title('算例A'); -% ylim([0 4]) -%Case B -subplot(1,3,2) -y=CaseBREV; -ymin=min(y); -ymax=max(y); -x=linspace(ymin,ymax,split_number); %将最大最小区间分成split_number个等分点(19等分),然后分别计算各个区间的个数 -yy=hist(y,x); %计算各个区间的个数 -% yy=yy/(sum(yy)*(ymax-ymin)/split_number); %计算各个区间的个数,除以总面积，总面积计算的方式为：所有小分割的面积的和即: -bar(x,yy) %画出概率密度分布图 -voltBarCaseBX=x; -voltBarCaseBY=yy; -xlabel('Error'); -ylabel('Number of buses'); -title('算例B'); -% ylim([0 4]) -%Case C -subplot(1,3,3) -y=CaseCREV; -ymin=min(y); -ymax=max(y); -x=linspace(ymin,ymax,split_number); %将最大最小区间分成split_number个等分点(19等分),然后分别计算各个区间的个数 -yy=hist(y,x); %计算各个区间的个数 -% yy=yy/(sum(yy)*(ymax-ymin)/split_number); %计算各个区间的个数,除以总面积，总面积计算的方式为：所有小分割的面积的和即: -bar(x,yy) %画出概率密度分布图 -voltBarCaseCX=x; -voltBarCaseCY=yy; -xlabel('Error'); -ylabel('Number of buses'); -title('算例C'); -% ylim([0 4]) -figure('Name','相角直方图') -% 相角 -split_number=20; -%Case A -subplot(2,2,1) -y=CaseAREA; -ymin=min(y); -ymax=max(y); -x=linspace(ymin,ymax,split_number); %将最大最小区间分成split_number个等分点(19等分),然后分别计算各个区间的个数 -yy=hist(y,x); %计算各个区间的个数 -% yy=yy/(sum(yy)*(ymax-ymin)/split_number); %计算各个区间的个数,除以总面积，总面积计算的方式为：所有小分割的面积的和即: -bar(x,yy) %画出概率密度分布图 -angelBarCaseAX=x; -angelBarCaseAY=yy; -% ylim([0 4]) -%Case B -subplot(2,2,2) -y=CaseBREA; -ymin=min(y); -ymax=max(y); -x=linspace(ymin,ymax,split_number); %将最大最小区间分成split_number个等分点(19等分),然后分别计算各个区间的个数 -yy=hist(y,x); %计算各个区间的个数 -% yy=yy/(sum(yy)*(ymax-ymin)/split_number); %计算各个区间的个数,除以总面积，总面积计算的方式为：所有小分割的面积的和即: -bar(x,yy) %画出概率密度分布图 -angelBarCaseBX=x; -angelBarCaseBY=yy; -% ylim([0 4]) -%Case C -subplot(2,2,3) -y=CaseCREA; -ymin=min(y); -ymax=max(y); -x=linspace(ymin,ymax,split_number); %将最大最小区间分成split_number个等分点(19等分),然后分别计算各个区间的个数 -yy=hist(y,x); %计算各个区间的个数 -angelBarCaseCX=x; -angelBarCaseCY=yy; -% yy=yy/(sum(yy)*(ymax-ymin)/split_number); %计算各个区间的个数,除以总面积，总面积计算的方式为：所有小分割的面积的和即: -bar(x,yy) %画出概率密度分布图 -%PD -figure('Name','有功直方图') -split_number=20; -%Case A -subplot(2,2,1) -y=CaseAREP; -ymin=min(y); -ymax=max(y); -x=linspace(ymin,ymax,split_number); %将最大最小区间分成split_number个等分点(19等分),然后分别计算各个区间的个数 -yy=hist(y,x); %计算各个区间的个数 -% yy=yy/(sum(yy)*(ymax-ymin)/split_number); %计算各个区间的个数,除以总面积，总面积计算的方式为：所有小分割的面积的和即: -bar(x,yy) %画出概率密度分布图 -PDBarCaseAX=x; -PDBarCaseAY=yy; -% ylim([0 4]) -%Case B -subplot(2,2,2) -y=CaseBREP; -ymin=min(y); -ymax=max(y); -x=linspace(ymin,ymax,split_number); %将最大最小区间分成split_number个等分点(19等分),然后分别计算各个区间的个数 -yy=hist(y,x); %计算各个区间的个数 -% yy=yy/(sum(yy)*(ymax-ymin)/split_number); %计算各个区间的个数,除以总面积，总面积计算的方式为：所有小分割的面积的和即: -bar(x,yy) %画出概率密度分布图 -PDBarCaseBX=x; -PDBarCaseBY=yy; -% ylim([0 4]) -%Case C -subplot(2,2,3) -y=CaseCREP; -ymin=min(y); -ymax=max(y); -x=linspace(ymin,ymax,split_number); %将最大最小区间分成split_number个等分点(19等分),然后分别计算各个区间的个数 -yy=hist(y,x); %计算各个区间的个数 -% yy=yy/(sum(yy)*(ymax-ymin)/split_number); %计算各个区间的个数,除以总面积，总面积计算的方式为：所有小分割的面积的和即: -bar(x,yy) %画出概率密度分布图 -PDBarCaseCX=x; -PDBarCaseCY=yy; -%QD -figure('Name','无功直方图') -split_number=20; -%Case A -subplot(2,2,1) -y=CaseAREQ; -ymin=min(y); -ymax=max(y); -x=linspace(ymin,ymax,split_number); %将最大最小区间分成split_number个等分点(19等分),然后分别计算各个区间的个数 -yy=hist(y,x); %计算各个区间的个数 -% yy=yy/(sum(yy)*(ymax-ymin)/split_number); %计算各个区间的个数,除以总面积，总面积计算的方式为：所有小分割的面积的和即: -bar(x,yy) %画出概率密度分布图 -QDBarCaseAX=x; -QDBarCaseAY=yy; -% ylim([0 4]) -%Case B -subplot(2,2,2) -y=CaseBREQ; -ymin=min(y); -ymax=max(y); -x=linspace(ymin,ymax,split_number); %将最大最小区间分成split_number个等分点(19等分),然后分别计算各个区间的个数 -yy=hist(y,x); %计算各个区间的个数 -% yy=yy/(sum(yy)*(ymax-ymin)/split_number); %计算各个区间的个数,除以总面积，总面积计算的方式为：所有小分割的面积的和即: -bar(x,yy) %画出概率密度分布图 -QDBarCaseBX=x; -QDBarCaseBY=yy; -% ylim([0 4]) -%Case C -subplot(2,2,3) -y=CaseCREQ; -ymin=min(y); -ymax=max(y); -x=linspace(ymin,ymax,split_number); %将最大最小区间分成split_number个等分点(19等分),然后分别计算各个区间的个数 -yy=hist(y,x); %计算各个区间的个数 -% yy=yy/(sum(yy)*(ymax-ymin)/split_number); %计算各个区间的个数,除以总面积，总面积计算的方式为：所有小分割的面积的和即: -bar(x,yy) %画出概率密度分布图 -QDBarCaseCX=x; -QDBarCaseCY=yy; -%画收敛曲线 -fz=find(abs(plotGapA)==0); -% fz=fz(1); -figure('Name','互补曲线') -plot(1:fz-1,plotGapA(1:fz-1)); -figure('Name','最大不平衡量'); -%% 最大不平衡量 -% maxDismatchPQ = [0.3123e-10 0.1497e-10 0.7351e-10; 0.6854e-10 0.1973e-10 0.5824e-10]; -% bar(maxDismatchPQ); -%% 计算时间 -calTime=[70.16 29.68; 68.48 31.661; 65.156 30.08;]; -bar(calTime,'stacked'); -% figure(); -% DeviationFigure(2:33,[CaseAREV(2:end);CaseBREV(2:end);CaseCREV(2:end)],[CaseAREA(2:end);CaseBREA(2:end);CaseCREA(2:end)],[CaseAREP(2:end),CaseBREP(2:end),CaseCREP(2:end)],[CaseAREQ(2:end),CaseBREQ(2:end),CaseCREQ(2:end)]);%,[CaseAREA;CaseBREA;CaseCREA],[CaseAREV;CaseBREV;CaseCREV],[CaseAREV;CaseBREV;CaseCREV]); -% VoltBar(voltBarCaseAX,voltBarCaseAY,voltBarCaseBX,voltBarCaseBY,voltBarCaseCX,voltBarCaseCY); -% AngelBar(angelBarCaseAX,angelBarCaseBY,angelBarCaseCX,angelBarCaseAY,angelBarCaseBX,angelBarCaseCY); -% PDBar(PDBarCaseAX,PDBarCaseAY,PDBarCaseBX,PDBarCaseBY,PDBarCaseCX,PDBarCaseCY); -% QDBar(QDBarCaseAX,QDBarCaseAY,QDBarCaseBX,QDBarCaseBY,QDBarCaseCX,QDBarCaseCY); -% MaxErrorFigure() -% MaxBoundErrorFigure(); -% DrawLoadProfile(); \ No newline at end of file +% % CaseAREQ=(QD-RealQD)./(RealQD+0.00001)*100;%Relative Error of QD in Case A +% CaseAREQ=(QD-RealQD)./RealQD*100;%Relative Error of QD in Case A +% CaseAREQ(1)=0; +% %真实值的 +% plot(1:length(CaseAREQ),(CaseAREQ),'k:','Marker','diamond'); +% %测量值的 +% % plot(1:length(CaseAREV),abs((QD0-RealQD)./(RealQD+0.00001)*100),'c.:','Marker','diamond'); +% box off; +% set(gca,'XTick',[1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33]) +% xlabel('节点号'); +% ylabel('误差%'); +% %计算Case A的误差 +% CaseAE=sqrt((sum(CaseAREV.^2)+sum(CaseAREA.^2)+sum(CaseAREP.^2)+sum(CaseAREQ.^2))/132); +% objA=full(sum(Vbi)+sum(PDbi)+sum(QDbi)); +% notZeros=find(PD0~=0); +% CaseA_SE=sum(((RealPD(notZeros)-PD(notZeros))./(RealPD(notZeros)*1)).^2)+sum(((RealQD(notZeros)-QD(notZeros))./(RealQD(notZeros)*1)).^2)+sum(((Volt-rVolt)./(rVolt*1)).^2); +% CaseA_SE=(CaseA_SE/(length(notZeros)*2+length(Volt)))^.5; +% sumCaseA_SE=sumCaseA_SE+CaseA_SE; +% % arrayA(1:19,I)=Vbi; +% % arrayA(21,I)=CaseAE*1000; +% %% Case B +% % [Busnum,Loadi,Volt,PD,QD,rVolt,UAngel,RealPD,RealQD,rUAngel,Vbi,PDbi,QDbi,plotGapB]=subOPF(setdiff(1:Busnum,[18,21,22,29]),PD0,QD0,mVolt,sigma);% +% [Busnum,Loadi,Volt,PD,QD,rVolt,UAngel,RealPD,RealQD,rUAngel,Vbi,PDbi,QDbi,plotGapB]=subOPF(setdiff(1:Busnum,[2,3,5,20,24,27,28,10,11,12,13]),PD0,QD0,mVolt,sigma);% +% subplot(4,1,1); +% hold on; +% % CaseBREV=(Volt-rVolt)*100;%Relative Error of Voltage in Case B +% CaseBREV=(Volt-rVolt);%Relative Error of Voltage in Case B +% plot(1:length(CaseBREV),(CaseBREV),'b.:','Marker','square'); +% subplot(4,1,2); +% hold on; +% % CaseBREA=(UAngel-rUAngel)*100;%Relative Error of Angle in Case B +% CaseBREA=(UAngel-rUAngel);%Relative Error of Angle in Case B +% plot(1:length(CaseBREA),(CaseBREA),'b:','Marker','square'); +% subplot(4,1,3); +% hold on; +% % CaseBREP=(PD-RealPD)./(RealPD+0.00001)*100;%Relative Error of PD in Case B +% CaseBREP=(PD-RealPD)./RealPD*100;%Relative Error of PD in Case B +% RealPD(1)=0; +% plot(1:length(CaseBREP),(CaseBREP),'b:','Marker','square'); +% subplot(4,1,4); +% hold on; +% % CaseBREQ=(QD-RealQD)./(RealQD+0.00001)*100;%Relative Error of QD in Case B +% CaseBREQ=(QD-RealQD)./RealQD*100;%Relative Error of QD in Case B +% CaseBREQ(1)=0; +% plot(1:length(CaseBREQ),(CaseBREQ),'b:','Marker','square'); +% CaseBE=sqrt((sum(CaseBREV.^2)+sum(CaseBREA.^2)+sum(CaseBREP.^2)+sum(CaseBREQ.^2))/132); +% objB=full(sum(Vbi)+sum(PDbi)+sum(QDbi)); +% noMeasurei=[2,3,5,20,24,27,28,10,11,12,13]; +% Measurei=setdiff(2:33,[2,3,5,20,24,27,28,10,11,12,13]); +% CaseB_SE=sum(((RealPD(Measurei)-PD(Measurei))./(RealPD(Measurei)*1)).^2)+sum(((RealQD(Measurei)-QD(Measurei))./(RealQD(Measurei)*1)).^2)+sum(((Volt(Measurei)-rVolt(Measurei))./(rVolt(Measurei)*1)).^2); +% CaseB_SE=CaseB_SE+sum(((RealPD(noMeasurei)-PD(noMeasurei))./(RealPD(noMeasurei)*1)).^2)+sum(((RealQD(noMeasurei)-QD(noMeasurei))./(RealQD(noMeasurei)*1)).^2)+sum(((Volt(noMeasurei)-rVolt(noMeasurei))./(rVolt(noMeasurei)*1)).^2); +% CaseB_SE=(CaseB_SE/(length(notZeros)+length(noMeasurei) +length(Volt)))^.5; +% sumCaseB_SE=sumCaseB_SE+CaseB_SE; +% %% Case C +% [Busnum,Loadi,Volt,PD,QD,rVolt,UAngel,RealPD,RealQD,rUAngel,Vbi,PDbi,QDbi,plotGapC]=subOPF([1:33],PD0,QD0,mVolt,sigma);% +% subplot(4,1,1); +% hold on; +% % CaseCREV=(Volt-rVolt)*100;%Relative Error of Voltage in Case C +% CaseCREV=(Volt-rVolt); +% plot(1:length(CaseCREV),(CaseCREV),'r.:','Marker','o'); +% subplot(4,1,2); +% hold on; +% % CaseCREA=(UAngel-rUAngel)*100;%Relative Error of Angle in Case C +% CaseCREA=(UAngel-rUAngel); +% plot(1:length(CaseCREA),(CaseCREA),'r:','Marker','o'); +% subplot(4,1,3); +% hold on; +% % CaseCREP=(PD-RealPD)./(RealPD+0.00001)*100;%Relative Error of PD in Case C +% CaseCREP=(PD-RealPD)./RealPD*100;%Relative Error of PD in Case C +% CaseCREP(1); +% plot(1:length(CaseCREP),(CaseCREP),'r:','Marker','o'); +% subplot(4,1,4); +% hold on; +% % CaseCREQ=(QD-RealQD)./(RealQD+0.00001)*100;%Relative Error of QD in Case C +% CaseCREQ=(QD-RealQD)./RealQD*100;%Relative Error of QD in Case C +% CaseCREQ(1)=0; +% plot(1:length(CaseCREQ),(CaseCREQ),'r:','Marker','o'); +% % 画legend +% subplot(4,1,1); +% % title('Voltage'); +% ld=legend('算例A','算例B','算例C'); +% set(ld,'Position',[0.847865087908145 0.786094477711244 0.0543595263724435 0.0605455755156354]); +% subplot(4,1,2); +% % title('Voltage Angle'); +% ld=legend('算例A','算例B','算例C'); +% set(ld,'Position',[0.847865087908145 0.586094477711244 0.0543595263724435 0.0605455755156354]); +% subplot(4,1,3); +% % title('Active load power'); +% ld=legend('算例A','算例B','算例C'); +% set(ld,'Position',[0.847865087908145 0.386094477711244 0.0543595263724435 0.0605455755156354]); +% subplot(4,1,4); +% % title('Reactive load power'); +% ld=legend('算例A','算例B','算例C'); +% set(ld,'Position',[0.847865087908145 0.186094477711244 0.0543595263724435 0.0605455755156354]); +% CaseCE=sqrt((sum(CaseCREV.^2)+sum(CaseCREA.^2)+sum(CaseCREP.^2)+sum(CaseCREQ.^2))/132); +% objC=full(sum(Vbi)+sum(PDbi)+sum(QDbi)); +% +% CaseC_SE=sum(((RealPD(notZeros)-PD(notZeros))./(RealPD(notZeros)*1)).^2)+sum(((RealQD(notZeros)-QD(notZeros))./(RealQD(notZeros)*1)).^2)+sum(((Volt-rVolt)./(rVolt*1)).^2); +% CaseC_SE=(CaseC_SE/(length(notZeros)*2+length(Volt)))^.5; +% sumCaseC_SE=sumCaseC_SE+CaseC_SE; +% +% % fprintf('目标函数值 %.2f\n',full(obj)); +% fprintf('Case A Case B Case C 的误差\n') +% fprintf('%f %f %f \n',CaseAE,CaseBE,CaseCE); +% fprintf('三个Case目标值\n') +% fprintf('%f\t%f\t%f \n',objA,objB,objC) +% %% 画测量值 +% % subplot(4,1,1); +% % plot(1:Busnum,mVolt-rVolt,'k.:','Marker','pentagram') +% % subplot(4,1,3); +% % plot(1:Busnum,(PD0-RealPD)./(RealPD+0.00001),'k:','Marker','pentagram'); +% % subplot(4,1,4); +% % plot(1:Busnum,(QD0-RealQD)./(RealQD+0.00001),'k:','Marker','pentagram'); +% %% 直方图 +% % 电压 +% figure('Name','电压直方图') +% split_number=20; +% %Case A +% subplot(1,3,1) +% y=CaseAREV; +% ymin=min(y); +% ymax=max(y); +% x=linspace(ymin,ymax,split_number); %将最大最小区间分成split_number个等分点(19等分),然后分别计算各个区间的个数 +% yy=hist(y,x); %计算各个区间的个数 +% % yy=yy/(sum(yy)*(ymax-ymin)/split_number); %计算各个区间的个数,除以总面积，总面积计算的方式为：所有小分割的面积的和即: +% bar(x,yy) %画出概率密度分布图 +% voltBarCaseAX=x; +% voltBarCaseAY=yy; +% xlabel('Error'); +% ylabel('Number of buses'); +% title('算例A'); +% % ylim([0 4]) +% %Case B +% subplot(1,3,2) +% y=CaseBREV; +% ymin=min(y); +% ymax=max(y); +% x=linspace(ymin,ymax,split_number); %将最大最小区间分成split_number个等分点(19等分),然后分别计算各个区间的个数 +% yy=hist(y,x); %计算各个区间的个数 +% % yy=yy/(sum(yy)*(ymax-ymin)/split_number); %计算各个区间的个数,除以总面积，总面积计算的方式为：所有小分割的面积的和即: +% bar(x,yy) %画出概率密度分布图 +% voltBarCaseBX=x; +% voltBarCaseBY=yy; +% xlabel('Error'); +% ylabel('Number of buses'); +% title('算例B'); +% % ylim([0 4]) +% %Case C +% subplot(1,3,3) +% y=CaseCREV; +% ymin=min(y); +% ymax=max(y); +% x=linspace(ymin,ymax,split_number); %将最大最小区间分成split_number个等分点(19等分),然后分别计算各个区间的个数 +% yy=hist(y,x); %计算各个区间的个数 +% % yy=yy/(sum(yy)*(ymax-ymin)/split_number); %计算各个区间的个数,除以总面积，总面积计算的方式为：所有小分割的面积的和即: +% bar(x,yy) %画出概率密度分布图 +% voltBarCaseCX=x; +% voltBarCaseCY=yy; +% xlabel('Error'); +% ylabel('Number of buses'); +% title('算例C'); +% % ylim([0 4]) +% figure('Name','相角直方图') +% % 相角 +% split_number=20; +% %Case A +% subplot(2,2,1) +% y=CaseAREA; +% ymin=min(y); +% ymax=max(y); +% x=linspace(ymin,ymax,split_number); %将最大最小区间分成split_number个等分点(19等分),然后分别计算各个区间的个数 +% yy=hist(y,x); %计算各个区间的个数 +% % yy=yy/(sum(yy)*(ymax-ymin)/split_number); %计算各个区间的个数,除以总面积，总面积计算的方式为：所有小分割的面积的和即: +% bar(x,yy) %画出概率密度分布图 +% angelBarCaseAX=x; +% angelBarCaseAY=yy; +% % ylim([0 4]) +% %Case B +% subplot(2,2,2) +% y=CaseBREA; +% ymin=min(y); +% ymax=max(y); +% x=linspace(ymin,ymax,split_number); %将最大最小区间分成split_number个等分点(19等分),然后分别计算各个区间的个数 +% yy=hist(y,x); %计算各个区间的个数 +% % yy=yy/(sum(yy)*(ymax-ymin)/split_number); %计算各个区间的个数,除以总面积，总面积计算的方式为：所有小分割的面积的和即: +% bar(x,yy) %画出概率密度分布图 +% angelBarCaseBX=x; +% angelBarCaseBY=yy; +% % ylim([0 4]) +% %Case C +% subplot(2,2,3) +% y=CaseCREA; +% ymin=min(y); +% ymax=max(y); +% x=linspace(ymin,ymax,split_number); %将最大最小区间分成split_number个等分点(19等分),然后分别计算各个区间的个数 +% yy=hist(y,x); %计算各个区间的个数 +% angelBarCaseCX=x; +% angelBarCaseCY=yy; +% % yy=yy/(sum(yy)*(ymax-ymin)/split_number); %计算各个区间的个数,除以总面积，总面积计算的方式为：所有小分割的面积的和即: +% bar(x,yy) %画出概率密度分布图 +% %PD +% figure('Name','有功直方图') +% split_number=20; +% %Case A +% subplot(2,2,1) +% y=CaseAREP; +% ymin=min(y); +% ymax=max(y); +% x=linspace(ymin,ymax,split_number); %将最大最小区间分成split_number个等分点(19等分),然后分别计算各个区间的个数 +% yy=hist(y,x); %计算各个区间的个数 +% % yy=yy/(sum(yy)*(ymax-ymin)/split_number); %计算各个区间的个数,除以总面积，总面积计算的方式为：所有小分割的面积的和即: +% bar(x,yy) %画出概率密度分布图 +% PDBarCaseAX=x; +% PDBarCaseAY=yy; +% % ylim([0 4]) +% %Case B +% subplot(2,2,2) +% y=CaseBREP; +% ymin=min(y); +% ymax=max(y); +% x=linspace(ymin,ymax,split_number); %将最大最小区间分成split_number个等分点(19等分),然后分别计算各个区间的个数 +% yy=hist(y,x); %计算各个区间的个数 +% % yy=yy/(sum(yy)*(ymax-ymin)/split_number); %计算各个区间的个数,除以总面积，总面积计算的方式为：所有小分割的面积的和即: +% bar(x,yy) %画出概率密度分布图 +% PDBarCaseBX=x; +% PDBarCaseBY=yy; +% % ylim([0 4]) +% %Case C +% subplot(2,2,3) +% y=CaseCREP; +% ymin=min(y); +% ymax=max(y); +% x=linspace(ymin,ymax,split_number); %将最大最小区间分成split_number个等分点(19等分),然后分别计算各个区间的个数 +% yy=hist(y,x); %计算各个区间的个数 +% % yy=yy/(sum(yy)*(ymax-ymin)/split_number); %计算各个区间的个数,除以总面积，总面积计算的方式为：所有小分割的面积的和即: +% bar(x,yy) %画出概率密度分布图 +% PDBarCaseCX=x; +% PDBarCaseCY=yy; +% %QD +% figure('Name','无功直方图') +% split_number=20; +% %Case A +% subplot(2,2,1) +% y=CaseAREQ; +% ymin=min(y); +% ymax=max(y); +% x=linspace(ymin,ymax,split_number); %将最大最小区间分成split_number个等分点(19等分),然后分别计算各个区间的个数 +% yy=hist(y,x); %计算各个区间的个数 +% % yy=yy/(sum(yy)*(ymax-ymin)/split_number); %计算各个区间的个数,除以总面积，总面积计算的方式为：所有小分割的面积的和即: +% bar(x,yy) %画出概率密度分布图 +% QDBarCaseAX=x; +% QDBarCaseAY=yy; +% % ylim([0 4]) +% %Case B +% subplot(2,2,2) +% y=CaseBREQ; +% ymin=min(y); +% ymax=max(y); +% x=linspace(ymin,ymax,split_number); %将最大最小区间分成split_number个等分点(19等分),然后分别计算各个区间的个数 +% yy=hist(y,x); %计算各个区间的个数 +% % yy=yy/(sum(yy)*(ymax-ymin)/split_number); %计算各个区间的个数,除以总面积，总面积计算的方式为：所有小分割的面积的和即: +% bar(x,yy) %画出概率密度分布图 +% QDBarCaseBX=x; +% QDBarCaseBY=yy; +% % ylim([0 4]) +% %Case C +% subplot(2,2,3) +% y=CaseCREQ; +% ymin=min(y); +% ymax=max(y); +% x=linspace(ymin,ymax,split_number); %将最大最小区间分成split_number个等分点(19等分),然后分别计算各个区间的个数 +% yy=hist(y,x); %计算各个区间的个数 +% % yy=yy/(sum(yy)*(ymax-ymin)/split_number); %计算各个区间的个数,除以总面积，总面积计算的方式为：所有小分割的面积的和即: +% bar(x,yy) %画出概率密度分布图 +% QDBarCaseCX=x; +% QDBarCaseCY=yy; +% %画收敛曲线 +% fz=find(abs(plotGapA)==0); +% % fz=fz(1); +% figure('Name','互补曲线') +% plot(1:fz-1,plotGapA(1:fz-1)); +% figure('Name','最大不平衡量'); +% %% 最大不平衡量 +% % maxDismatchPQ = [0.3123e-10 0.1497e-10 0.7351e-10; 0.6854e-10 0.1973e-10 0.5824e-10]; +% % bar(maxDismatchPQ); +% %% 计算时间 +% calTime=[70.16 29.68; 68.48 31.661; 65.156 30.08;]; +% bar(calTime,'stacked'); +% % figure(); +% % DeviationFigure(2:33,[CaseAREV(2:end);CaseBREV(2:end);CaseCREV(2:end)],[CaseAREA(2:end);CaseBREA(2:end);CaseCREA(2:end)],[CaseAREP(2:end),CaseBREP(2:end),CaseCREP(2:end)],[CaseAREQ(2:end),CaseBREQ(2:end),CaseCREQ(2:end)]);%,[CaseAREA;CaseBREA;CaseCREA],[CaseAREV;CaseBREV;CaseCREV],[CaseAREV;CaseBREV;CaseCREV]); +% % VoltBar(voltBarCaseAX,voltBarCaseAY,voltBarCaseBX,voltBarCaseBY,voltBarCaseCX,voltBarCaseCY); +% % AngelBar(angelBarCaseAX,angelBarCaseBY,angelBarCaseCX,angelBarCaseAY,angelBarCaseBX,angelBarCaseCY); +% % PDBar(PDBarCaseAX,PDBarCaseAY,PDBarCaseBX,PDBarCaseBY,PDBarCaseCX,PDBarCaseCY); +% % QDBar(QDBarCaseAX,QDBarCaseAY,QDBarCaseBX,QDBarCaseBY,QDBarCaseCX,QDBarCaseCY); +% % MaxErrorFigure() +% % MaxBoundErrorFigure(); +% % DrawLoadProfile(); + +%% PLOTING \ No newline at end of file diff --git a/laprnd.m b/laprnd.m new file mode 100644 index 0000000..ac7a10f --- /dev/null +++ b/laprnd.m @@ -0,0 +1,30 @@ +function y = laprnd(m, n, mu, sigma) +%LAPRND generate i.i.d. laplacian random number drawn from laplacian distribution +% with mean mu and standard deviation sigma. +% mu : mean +% sigma : standard deviation +% [m, n] : the dimension of y. +% Default mu = 0, sigma = 1. +% For more information, refer to +% http://en.wikipedia.org./wiki/Laplace_distribution + +% Author : Elvis Chen (bee33@sjtu.edu.cn) +% Date : 01/19/07 + +%Check inputs +if nargin < 2 + error('At least two inputs are required'); +end + +if nargin == 2 + mu = 0; sigma = 1; +end + +if nargin == 3 + sigma = 1; +end + +% Generate Laplacian noise +u = rand(m, n)-0.5; +b = sigma / sqrt(2); +y = mu - b * sign(u).* log(1- 2* abs(u)); diff --git a/subOPF.m b/subOPF.m index f0aef7c..f13f423 100644 --- a/subOPF.m +++ b/subOPF.m @@ -49,7 +49,7 @@ Gap=(Init_L*Init_Z'-Init_U*Init_W'); KK=0; plotGap=zeros(1,60); ContrlCount=size(Loadi,1)*2+Busnum*2+Busnum+length(Loadi)*2; -kmax=170; +kmax=30; Precision=Precision/1; %% 加误差 %找DG @@ -240,7 +240,7 @@ dQ(Balance)=0; mdP=max(dP) mdQ=max(dQ) %判断是否收敛 -if min(abs(Vbi-.1))>5e-005 +if min(abs(Vbi-0))>5e-005 isConverge=0; end if KK>=kmax