facat
/
pes2014-debug2-laplace
1
0
Fork 0
Code Issues Pull Requests Releases Wiki Activity

1.加入生成Laplace分布随机数的函数 

2.把代码调回最初形式。

Signed-off-by: dugg@lab-desk <dugg@lab-desk>
This commit is contained in:
dugg@lab-desk 2015-03-27 10:52:23 +08:00
parent 6b24827aee
commit ec8a4e8189
3 changed files with 373 additions and 338 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

677
OPF.m
View File

@ -29,12 +29,12 @@ for badDataNode=1:1
mVolt=rVolt.*(1+normrnd(0,sigma,length(rVolt),1))'; mVolt=rVolt.*(1+normrnd(0,sigma,length(rVolt),1))';
%¼ÓÔر£´æµÄ±äÁ¿ %¼ÓÔر£´æµÄ±äÁ¿
PD0=load('PD0'); % PD0=load('PD0');
PD0=PD0.PD0; % PD0=PD0.PD0;
QD0=load('QD0'); % QD0=load('QD0');
QD0=QD0.QD0; % QD0=QD0.QD0;
mVolt=load('mVolt'); % mVolt=load('mVolt');
mVolt=mVolt.mVolt; % mVolt=mVolt.mVolt;
% mVolt(3)=rVolt(3)*(1-sigma*6); % mVolt(3)=rVolt(3)*(1-sigma*6);
%% »­Case AµÄͼ %% »­Case AµÄͼ
@ -43,7 +43,7 @@ for badDataNode=1:1
% badDataResult(I,badDataNode)=sum(Vbi); % badDataResult(I,badDataNode)=sum(Vbi);
% badDataLocation(1:33,I)=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)))); % 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 if isConverge==0
continue; continue;
end end
@ -62,7 +62,7 @@ for badDataNode=1:1
nodeMaxDVAngle_t([1:18,23:33])=0; nodeMaxDVAngle_t([1:18,23:33])=0;
nodeMaxDVolt(nodeMaxDVolt<nodeMaxDVolt_t)=nodeMaxDVolt_t(nodeMaxDVolt<nodeMaxDVolt_t); nodeMaxDVolt(nodeMaxDVolt<nodeMaxDVolt_t)=nodeMaxDVolt_t(nodeMaxDVolt<nodeMaxDVolt_t);
nodeMaxDVAngle(nodeMaxDVAngle<nodeMaxDVAngle_t)=nodeMaxDVAngle_t(nodeMaxDVAngle<nodeMaxDVAngle_t); nodeMaxDVAngle(nodeMaxDVAngle<nodeMaxDVAngle_t)=nodeMaxDVAngle_t(nodeMaxDVAngle<nodeMaxDVAngle_t);
if loopN>=500 if loopN>=1
% nodeMaxDVolt(badDataNode)=maxDVolt; % nodeMaxDVolt(badDataNode)=maxDVolt;
% nodeMaxDVAngle(badDataNode)=maxDVAngle; % nodeMaxDVAngle(badDataNode)=maxDVAngle;
break; break;
@ -72,335 +72,340 @@ for badDataNode=1:1
end end
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); %% PLOTING
CaseC_SE=(CaseC_SE/(length(notZeros)*2+length(Volt)))^.5;
sumCaseC_SE=sumCaseC_SE+CaseC_SE;
% fprintf('目标函数值 %.2f\n',full(obj)); % % end
fprintf('Case A Case B Case C 的误差\n') % 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]);
fprintf('%f %f %f \n',CaseAE,CaseBE,CaseCE); % CaseAREV=(Volt-rVolt);%Relative Error of Voltage in Case A
fprintf('三个Case目标值\n') % % CaseAREV=CaseAREV(2:end)*100;
fprintf('%f\t%f\t%f \n',objA,objB,objC) % %真实值的
%% 画测量值 % plot(1:length(CaseAREV),(CaseAREV),'k.:','Marker','diamond');
% subplot(4,1,1); % %测量值的
% plot(1:Busnum,mVolt-rVolt,'k.:','Marker','pentagram') % % 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); % 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); % subplot(4,1,4);
% plot(1:Busnum,(QD0-RealQD)./(RealQD+0.00001),'k:','Marker','pentagram'); % % 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;
figure('Name','电压直方图') % %真实值的
split_number=20; % plot(1:length(CaseAREQ),(CaseAREQ),'k:','Marker','diamond');
%Case A % %测量值的
subplot(1,3,1) % % plot(1:length(CaseAREV),abs((QD0-RealQD)./(RealQD+0.00001)*100),'c.:','Marker','diamond');
y=CaseAREV; % box off;
ymin=min(y); % 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])
ymax=max(y); % xlabel('节点号');
x=linspace(ymin,ymax,split_number); %将最大最小区间分成split_number个等分点(19等分),然后分别计算各个区间的个数 % ylabel('误差%');
yy=hist(y,x); %计算各个区间的个数 % %计算Case A的误差
% yy=yy/(sum(yy)*(ymax-ymin)/split_number); %计算各个区间的个数,除以总面积，总面积计算的方式为：所有小分割的面积的和即: % CaseAE=sqrt((sum(CaseAREV.^2)+sum(CaseAREA.^2)+sum(CaseAREP.^2)+sum(CaseAREQ.^2))/132);
bar(x,yy) %画出概率密度分布图 % objA=full(sum(Vbi)+sum(PDbi)+sum(QDbi));
voltBarCaseAX=x; % notZeros=find(PD0~=0);
voltBarCaseAY=yy; % 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);
xlabel('Error'); % CaseA_SE=(CaseA_SE/(length(notZeros)*2+length(Volt)))^.5;
ylabel('Number of buses'); % sumCaseA_SE=sumCaseA_SE+CaseA_SE;
title('算例A'); % % arrayA(1:19,I)=Vbi;
% ylim([0 4]) % % arrayA(21,I)=CaseAE*1000;
%Case B % %% Case B
subplot(1,3,2) % % [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);%
y=CaseBREV; % [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);%
ymin=min(y); % subplot(4,1,1);
ymax=max(y); % hold on;
x=linspace(ymin,ymax,split_number); %将最大最小区间分成split_number个等分点(19等分),然后分别计算各个区间的个数 % % CaseBREV=(Volt-rVolt)*100;%Relative Error of Voltage in Case B
yy=hist(y,x); %计算各个区间的个数 % CaseBREV=(Volt-rVolt);%Relative Error of Voltage in Case B
% yy=yy/(sum(yy)*(ymax-ymin)/split_number); %计算各个区间的个数,除以总面积，总面积计算的方式为：所有小分割的面积的和即: % plot(1:length(CaseBREV),(CaseBREV),'b.:','Marker','square');
bar(x,yy) %画出概率密度分布图 % subplot(4,1,2);
voltBarCaseBX=x; % hold on;
voltBarCaseBY=yy; % % CaseBREA=(UAngel-rUAngel)*100;%Relative Error of Angle in Case B
xlabel('Error'); % CaseBREA=(UAngel-rUAngel);%Relative Error of Angle in Case B
ylabel('Number of buses'); % plot(1:length(CaseBREA),(CaseBREA),'b:','Marker','square');
title('算例B'); % subplot(4,1,3);
% ylim([0 4]) % hold on;
%Case C % % CaseBREP=(PD-RealPD)./(RealPD+0.00001)*100;%Relative Error of PD in Case B
subplot(1,3,3) % CaseBREP=(PD-RealPD)./RealPD*100;%Relative Error of PD in Case B
y=CaseCREV; % RealPD(1)=0;
ymin=min(y); % plot(1:length(CaseBREP),(CaseBREP),'b:','Marker','square');
ymax=max(y); % subplot(4,1,4);
x=linspace(ymin,ymax,split_number); %将最大最小区间分成split_number个等分点(19等分),然后分别计算各个区间的个数 % hold on;
yy=hist(y,x); %计算各个区间的个数 % % CaseBREQ=(QD-RealQD)./(RealQD+0.00001)*100;%Relative Error of QD in Case B
% yy=yy/(sum(yy)*(ymax-ymin)/split_number); %计算各个区间的个数,除以总面积，总面积计算的方式为：所有小分割的面积的和即: % CaseBREQ=(QD-RealQD)./RealQD*100;%Relative Error of QD in Case B
bar(x,yy) %画出概率密度分布图 % CaseBREQ(1)=0;
voltBarCaseCX=x; % plot(1:length(CaseBREQ),(CaseBREQ),'b:','Marker','square');
voltBarCaseCY=yy; % CaseBE=sqrt((sum(CaseBREV.^2)+sum(CaseBREA.^2)+sum(CaseBREP.^2)+sum(CaseBREQ.^2))/132);
xlabel('Error'); % objB=full(sum(Vbi)+sum(PDbi)+sum(QDbi));
ylabel('Number of buses'); % noMeasurei=[2,3,5,20,24,27,28,10,11,12,13];
title('算例C'); % Measurei=setdiff(2:33,[2,3,5,20,24,27,28,10,11,12,13]);
% ylim([0 4]) % 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);
figure('Name','相角直方图') % 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;
split_number=20; % sumCaseB_SE=sumCaseB_SE+CaseB_SE;
%Case A % %% Case C
subplot(2,2,1) % [Busnum,Loadi,Volt,PD,QD,rVolt,UAngel,RealPD,RealQD,rUAngel,Vbi,PDbi,QDbi,plotGapC]=subOPF([1:33],PD0,QD0,mVolt,sigma);%
y=CaseAREA; % subplot(4,1,1);
ymin=min(y); % hold on;
ymax=max(y); % % CaseCREV=(Volt-rVolt)*100;%Relative Error of Voltage in Case C
x=linspace(ymin,ymax,split_number); %将最大最小区间分成split_number个等分点(19等分),然后分别计算各个区间的个数 % CaseCREV=(Volt-rVolt);
yy=hist(y,x); %计算各个区间的个数 % plot(1:length(CaseCREV),(CaseCREV),'r.:','Marker','o');
% yy=yy/(sum(yy)*(ymax-ymin)/split_number); %计算各个区间的个数,除以总面积，总面积计算的方式为：所有小分割的面积的和即: % subplot(4,1,2);
bar(x,yy) %画出概率密度分布图 % hold on;
angelBarCaseAX=x; % % CaseCREA=(UAngel-rUAngel)*100;%Relative Error of Angle in Case C
angelBarCaseAY=yy; % CaseCREA=(UAngel-rUAngel);
% ylim([0 4]) % plot(1:length(CaseCREA),(CaseCREA),'r:','Marker','o');
%Case B % subplot(4,1,3);
subplot(2,2,2) % hold on;
y=CaseBREA; % % CaseCREP=(PD-RealPD)./(RealPD+0.00001)*100;%Relative Error of PD in Case C
ymin=min(y); % CaseCREP=(PD-RealPD)./RealPD*100;%Relative Error of PD in Case C
ymax=max(y); % CaseCREP(1);
x=linspace(ymin,ymax,split_number); %将最大最小区间分成split_number个等分点(19等分),然后分别计算各个区间的个数 % plot(1:length(CaseCREP),(CaseCREP),'r:','Marker','o');
yy=hist(y,x); %计算各个区间的个数 % subplot(4,1,4);
% yy=yy/(sum(yy)*(ymax-ymin)/split_number); %计算各个区间的个数,除以总面积，总面积计算的方式为：所有小分割的面积的和即: % hold on;
bar(x,yy) %画出概率密度分布图 % % CaseCREQ=(QD-RealQD)./(RealQD+0.00001)*100;%Relative Error of QD in Case C
angelBarCaseBX=x; % CaseCREQ=(QD-RealQD)./RealQD*100;%Relative Error of QD in Case C
angelBarCaseBY=yy; % CaseCREQ(1)=0;
% ylim([0 4]) % plot(1:length(CaseCREQ),(CaseCREQ),'r:','Marker','o');
%Case C % % 画legend
subplot(2,2,3) % subplot(4,1,1);
y=CaseCREA; % % title('Voltage');
ymin=min(y); % ld=legend('算例A','算例B','算例C');
ymax=max(y); % set(ld,'Position',[0.847865087908145 0.786094477711244 0.0543595263724435 0.0605455755156354]);
x=linspace(ymin,ymax,split_number); %将最大最小区间分成split_number个等分点(19等分),然后分别计算各个区间的个数 % subplot(4,1,2);
yy=hist(y,x); %计算各个区间的个数 % % title('Voltage Angle');
angelBarCaseCX=x; % ld=legend('算例A','算例B','算例C');
angelBarCaseCY=yy; % set(ld,'Position',[0.847865087908145 0.586094477711244 0.0543595263724435 0.0605455755156354]);
% yy=yy/(sum(yy)*(ymax-ymin)/split_number); %计算各个区间的个数,除以总面积，总面积计算的方式为：所有小分割的面积的和即: % subplot(4,1,3);
bar(x,yy) %画出概率密度分布图 % % title('Active load power');
%PD % ld=legend('算例A','算例B','算例C');
figure('Name','有功直方图') % set(ld,'Position',[0.847865087908145 0.386094477711244 0.0543595263724435 0.0605455755156354]);
split_number=20; % subplot(4,1,4);
%Case A % % title('Reactive load power');
subplot(2,2,1) % ld=legend('算例A','算例B','算例C');
y=CaseAREP; % set(ld,'Position',[0.847865087908145 0.186094477711244 0.0543595263724435 0.0605455755156354]);
ymin=min(y); % CaseCE=sqrt((sum(CaseCREV.^2)+sum(CaseCREA.^2)+sum(CaseCREP.^2)+sum(CaseCREQ.^2))/132);
ymax=max(y); % objC=full(sum(Vbi)+sum(PDbi)+sum(QDbi));
x=linspace(ymin,ymax,split_number); %将最大最小区间分成split_number个等分点(19等分),然后分别计算各个区间的个数 %
yy=hist(y,x); %计算各个区间的个数 % 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);
% yy=yy/(sum(yy)*(ymax-ymin)/split_number); %计算各个区间的个数,除以总面积，总面积计算的方式为：所有小分割的面积的和即: % CaseC_SE=(CaseC_SE/(length(notZeros)*2+length(Volt)))^.5;
bar(x,yy) %画出概率密度分布图 % sumCaseC_SE=sumCaseC_SE+CaseC_SE;
PDBarCaseAX=x; %
PDBarCaseAY=yy; % % fprintf('目标函数值 %.2f\n',full(obj));
% ylim([0 4]) % fprintf('Case A Case B Case C 的误差\n')
%Case B % fprintf('%f %f %f \n',CaseAE,CaseBE,CaseCE);
subplot(2,2,2) % fprintf('三个Case目标值\n')
y=CaseBREP; % fprintf('%f\t%f\t%f \n',objA,objB,objC)
ymin=min(y); % %% 画测量值
ymax=max(y); % % subplot(4,1,1);
x=linspace(ymin,ymax,split_number); %将最大最小区间分成split_number个等分点(19等分),然后分别计算各个区间的个数 % % plot(1:Busnum,mVolt-rVolt,'k.:','Marker','pentagram')
yy=hist(y,x); %计算各个区间的个数 % % subplot(4,1,3);
% yy=yy/(sum(yy)*(ymax-ymin)/split_number); %计算各个区间的个数,除以总面积，总面积计算的方式为：所有小分割的面积的和即: % % plot(1:Busnum,(PD0-RealPD)./(RealPD+0.00001),'k:','Marker','pentagram');
bar(x,yy) %画出概率密度分布图 % % subplot(4,1,4);
PDBarCaseBX=x; % % plot(1:Busnum,(QD0-RealQD)./(RealQD+0.00001),'k:','Marker','pentagram');
PDBarCaseBY=yy; % %% 直方图
% ylim([0 4]) % % 电压
%Case C % figure('Name','电压直方图')
subplot(2,2,3) % split_number=20;
y=CaseCREP; % %Case A
ymin=min(y); % subplot(1,3,1)
ymax=max(y); % y=CaseAREV;
x=linspace(ymin,ymax,split_number); %将最大最小区间分成split_number个等分点(19等分),然后分别计算各个区间的个数 % ymin=min(y);
yy=hist(y,x); %计算各个区间的个数 % ymax=max(y);
% yy=yy/(sum(yy)*(ymax-ymin)/split_number); %计算各个区间的个数,除以总面积，总面积计算的方式为：所有小分割的面积的和即: % x=linspace(ymin,ymax,split_number); %将最大最小区间分成split_number个等分点(19等分),然后分别计算各个区间的个数
bar(x,yy) %画出概率密度分布图 % yy=hist(y,x); %计算各个区间的个数
PDBarCaseCX=x; % % yy=yy/(sum(yy)*(ymax-ymin)/split_number); %计算各个区间的个数,除以总面积，总面积计算的方式为：所有小分割的面积的和即:
PDBarCaseCY=yy; % bar(x,yy) %画出概率密度分布图
%QD % voltBarCaseAX=x;
figure('Name','无功直方图') % voltBarCaseAY=yy;
split_number=20; % xlabel('Error');
%Case A % ylabel('Number of buses');
subplot(2,2,1) % title('算例A');
y=CaseAREQ; % % ylim([0 4])
ymin=min(y); % %Case B
ymax=max(y); % subplot(1,3,2)
x=linspace(ymin,ymax,split_number); %将最大最小区间分成split_number个等分点(19等分),然后分别计算各个区间的个数 % y=CaseBREV;
yy=hist(y,x); %计算各个区间的个数 % ymin=min(y);
% yy=yy/(sum(yy)*(ymax-ymin)/split_number); %计算各个区间的个数,除以总面积，总面积计算的方式为：所有小分割的面积的和即: % ymax=max(y);
bar(x,yy) %画出概率密度分布图 % x=linspace(ymin,ymax,split_number); %将最大最小区间分成split_number个等分点(19等分),然后分别计算各个区间的个数
QDBarCaseAX=x; % yy=hist(y,x); %计算各个区间的个数
QDBarCaseAY=yy; % % yy=yy/(sum(yy)*(ymax-ymin)/split_number); %计算各个区间的个数,除以总面积，总面积计算的方式为：所有小分割的面积的和即:
% ylim([0 4]) % bar(x,yy) %画出概率密度分布图
%Case B % voltBarCaseBX=x;
subplot(2,2,2) % voltBarCaseBY=yy;
y=CaseBREQ; % xlabel('Error');
ymin=min(y); % ylabel('Number of buses');
ymax=max(y); % title('算例B');
x=linspace(ymin,ymax,split_number); %将最大最小区间分成split_number个等分点(19等分),然后分别计算各个区间的个数 % % ylim([0 4])
yy=hist(y,x); %计算各个区间的个数 % %Case C
% yy=yy/(sum(yy)*(ymax-ymin)/split_number); %计算各个区间的个数,除以总面积，总面积计算的方式为：所有小分割的面积的和即: % subplot(1,3,3)
bar(x,yy) %画出概率密度分布图 % y=CaseCREV;
QDBarCaseBX=x; % ymin=min(y);
QDBarCaseBY=yy; % ymax=max(y);
% ylim([0 4]) % x=linspace(ymin,ymax,split_number); %将最大最小区间分成split_number个等分点(19等分),然后分别计算各个区间的个数
%Case C % yy=hist(y,x); %计算各个区间的个数
subplot(2,2,3) % % yy=yy/(sum(yy)*(ymax-ymin)/split_number); %计算各个区间的个数,除以总面积，总面积计算的方式为：所有小分割的面积的和即:
y=CaseCREQ; % bar(x,yy) %画出概率密度分布图
ymin=min(y); % voltBarCaseCX=x;
ymax=max(y); % voltBarCaseCY=yy;
x=linspace(ymin,ymax,split_number); %将最大最小区间分成split_number个等分点(19等分),然后分别计算各个区间的个数 % xlabel('Error');
yy=hist(y,x); %计算各个区间的个数 % ylabel('Number of buses');
% yy=yy/(sum(yy)*(ymax-ymin)/split_number); %计算各个区间的个数,除以总面积，总面积计算的方式为：所有小分割的面积的和即: % title('算例C');
bar(x,yy) %画出概率密度分布图 % % ylim([0 4])
QDBarCaseCX=x; % figure('Name','相角直方图')
QDBarCaseCY=yy; % % 相角
%画收敛曲线 % split_number=20;
fz=find(abs(plotGapA)==0); % %Case A
% fz=fz(1); % subplot(2,2,1)
figure('Name','互补曲线') % y=CaseAREA;
plot(1:fz-1,plotGapA(1:fz-1)); % ymin=min(y);
figure('Name','最大不平衡量'); % ymax=max(y);
%% 最大不平衡量 % x=linspace(ymin,ymax,split_number); %将最大最小区间分成split_number个等分点(19等分),然后分别计算各个区间的个数
% maxDismatchPQ = [0.3123e-10 0.1497e-10 0.7351e-10; 0.6854e-10 0.1973e-10 0.5824e-10]; % yy=hist(y,x); %计算各个区间的个数
% bar(maxDismatchPQ); % % yy=yy/(sum(yy)*(ymax-ymin)/split_number); %计算各个区间的个数,除以总面积，总面积计算的方式为：所有小分割的面积的和即:
%% 计算时间 % bar(x,yy) %画出概率密度分布图
calTime=[70.16 29.68; 68.48 31.661; 65.156 30.08;]; % angelBarCaseAX=x;
bar(calTime,'stacked'); % angelBarCaseAY=yy;
% figure(); % % ylim([0 4])
% 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]); % %Case B
% VoltBar(voltBarCaseAX,voltBarCaseAY,voltBarCaseBX,voltBarCaseBY,voltBarCaseCX,voltBarCaseCY); % subplot(2,2,2)
% AngelBar(angelBarCaseAX,angelBarCaseBY,angelBarCaseCX,angelBarCaseAY,angelBarCaseBX,angelBarCaseCY); % y=CaseBREA;
% PDBar(PDBarCaseAX,PDBarCaseAY,PDBarCaseBX,PDBarCaseBY,PDBarCaseCX,PDBarCaseCY); % ymin=min(y);
% QDBar(QDBarCaseAX,QDBarCaseAY,QDBarCaseBX,QDBarCaseBY,QDBarCaseCX,QDBarCaseCY); % ymax=max(y);
% MaxErrorFigure() % x=linspace(ymin,ymax,split_number); %将最大最小区间分成split_number个等分点(19等分),然后分别计算各个区间的个数
% MaxBoundErrorFigure(); % yy=hist(y,x); %计算各个区间的个数
% DrawLoadProfile(); % % 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

30
laprnd.m Normal file
View File

@ -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));

4
subOPF.m
View File

@ -49,7 +49,7 @@ Gap=(Init_L*Init_Z'-Init_U*Init_W');
KK=0; KK=0;
plotGap=zeros(1,60); plotGap=zeros(1,60);
ContrlCount=size(Loadi,1)*2+Busnum*2+Busnum+length(Loadi)*2; ContrlCount=size(Loadi,1)*2+Busnum*2+Busnum+length(Loadi)*2;
kmax=170; kmax=30;
Precision=Precision/1; Precision=Precision/1;
%% 加误差 %% 加误差
%找DG %找DG
@ -240,7 +240,7 @@ dQ(Balance)=0;
mdP=max(dP) mdP=max(dP)
mdQ=max(dQ) mdQ=max(dQ)
%判断是否收敛 %判断是否收敛
if min(abs(Vbi-.1))>5e-005 if min(abs(Vbi-0))>5e-005
isConverge=0; isConverge=0;
end end
if KK>=kmax if KK>=kmax