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

最终选择用这个组数据和这个评价指标 

Signed-off-by: dugg@lab-desk <dugg@lab-desk>
This commit is contained in:
dugg@lab-desk 2015-02-11 14:35:07 +08:00
parent 5a752ca4fb
commit 8f5fe67a6a
6 changed files with 70 additions and 63 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

105
OPF.m
View File

@ -29,21 +29,21 @@ 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(badDataNode)=rVolt(badDataNode)*(1-sigma*6);
%% 画Case A的图
% figure('Color',[1 1 1]);
[Busnum,Loadi,Volt,PD,QD,rVolt,UAngel,RealPD,RealQD,rUAngel,Vbi,PDbi,QDbi,plotGapA,isConverge]=subOPF([19:22],PD0,QD0,mVolt,sigma);%ČŤ˛żÓĐ
[Busnum,Loadi,Volt,PD,QD,rVolt,UAngel,RealPD,RealQD,rUAngel,Vbi,PDbi,QDbi,plotGapA,isConverge]=subOPF([],PD0,QD0,mVolt,sigma);%ČŤ˛żÓĐ
% 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
@ -74,10 +74,10 @@ for badDataNode=1:1
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)*100;%Relative Error of Voltage in Case A
CaseAREV=(Volt-rVolt);%Relative Error of Voltage in Case A
% CaseAREV=CaseAREV(2:end)*100;
%真实值的
plot(1:length(CaseAREV),abs(CaseAREV),'k.:','Marker','diamond');
plot(1:length(CaseAREV),(CaseAREV),'k.:','Marker','diamond');
%测量值的
% plot(1:length(CaseAREV),abs((mVolt-rVolt)*100),'c.:','Marker','diamond');
box off;
@ -85,19 +85,21 @@ 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
xlabel('节点号');
ylabel('误差%');
subplot(4,1,2);
CaseAREA=(UAngel-rUAngel)*100;%Relative Error of Angle in Case A
% CaseAREA=CaseAREA(2:end)*100;
% 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),abs(CaseAREA),'k:','Marker','diamond');
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=CaseAREP(2:end)*100;
% 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),abs(CaseAREP),'k:','Marker','diamond');
plot(1:length(CaseAREP),(CaseAREP),'k:','Marker','diamond');
%测量值的
% plot(1:length(CaseAREV),abs((PD0-RealPD)./(RealPD+0.00001)*100),'c.:','Marker','diamond');
box off;
@ -105,10 +107,11 @@ 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
xlabel('节点号');
ylabel('误差%');
subplot(4,1,4);
CaseAREQ=(QD-RealQD)./(RealQD+0.00001)*100;%Relative Error of QD in Case A
% CaseAREQ=CaseAREQ(2:end)*100;
% 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),abs(CaseAREQ),'k:','Marker','diamond');
plot(1:length(CaseAREQ),(CaseAREQ),'k:','Marker','diamond');
%测量值的
% plot(1:length(CaseAREV),abs((QD0-RealQD)./(RealQD+0.00001)*100),'c.:','Marker','diamond');
box off;
@ -119,7 +122,7 @@ ylabel('
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(((PD0(notZeros)-PD(notZeros))./(PD0(notZeros)*sigma)).^2)+sum(((QD0(notZeros)-QD(notZeros))./(QD0(notZeros)*sigma)).^2)+sum(((Volt-mVolt)./(mVolt*sigma)).^2);
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;
@ -129,54 +132,58 @@ sumCaseA_SE=sumCaseA_SE+CaseA_SE;
[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=CaseBREV(2:end)*100;
plot(1:length(CaseBREV),abs(CaseBREV),'b.:','Marker','square');
% 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=CaseBREA(2:end)*100;
plot(1:length(CaseBREA),abs(CaseBREA),'b:','Marker','square');
% 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=CaseBREP(2:end)*100;
plot(1:length(CaseBREP),abs(CaseBREP),'b:','Marker','square');
% CaseBREP=(PD-RealPD)./(RealPD+0.00001)*100;%Relative Error of PD in Case B
CaseBREP=(PD-RealPD)./RealPD;%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=CaseBREQ(2:end)*100;
plot(1:length(CaseBREQ),abs(CaseBREQ),'b:','Marker','square');
% CaseBREQ=(QD-RealQD)./(RealQD+0.00001)*100;%Relative Error of QD in Case B
CaseBREQ=(QD-RealQD)./RealQD;%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(((PD0(Measurei)-PD(Measurei))./(PD0(Measurei)*sigma)).^2)+sum(((QD0(Measurei)-QD(Measurei))./(QD0(Measurei)*sigma)).^2)+sum(((Volt-mVolt)./(mVolt*sigma)).^2);
CaseB_SE=CaseB_SE+sum(((PD0(noMeasurei)-PD(noMeasurei))./(PD0(noMeasurei)*20*sigma)).^2)+sum(((QD0(noMeasurei)-QD(noMeasurei))./(QD0(noMeasurei)*20*sigma)).^2);
CaseB_SE=(CaseB_SE/(length(notZeros)*2+length(Volt)))^.5;
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=CaseCREV(2:end)*100;
plot(1:length(CaseCREV),abs(CaseCREV),'r.:','Marker','o');
% 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=CaseCREA(2:end)*100;
plot(1:length(CaseCREA),abs(CaseCREA),'r:','Marker','o');
% 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=CaseCREP(2:end)*100;
plot(1:length(CaseCREP),abs(CaseCREP),'r:','Marker','o');
% 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=CaseCREQ(2:end)*100;
plot(1:length(CaseCREQ),abs(CaseCREQ),'r:','Marker','o');
% 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');
@ -197,7 +204,7 @@ set(ld,'Position',[0.847865087908145 0.186094477711244 0.0543595263724435 0.0605
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(((PD0(notZeros)-PD(notZeros))./(PD0(notZeros)*sigma)).^2)+sum(((QD0(notZeros)-QD(notZeros))./(QD0(notZeros)*sigma)).^2)+sum(((Volt-mVolt)./(mVolt*sigma)).^2);
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;
@ -354,7 +361,7 @@ yy=hist(y,x); %
bar(x,yy) %画出概率密度分布图
%画收敛曲线
fz=find(abs(plotGapA)==0);
fz=fz(1);
% fz=fz(1);
figure('Name','互补曲线')
plot(1:fz-1,plotGapA(1:fz-1));
figure('Name','最大不平衡量');

BIN
PD0.mat
View File

Binary file not shown.

BIN
QD0.mat
View File

Binary file not shown.

BIN
mVolt.mat
View File

Binary file not shown.

2
openfile2.m
View File

@ -59,7 +59,7 @@ QD=sparse(QD);
%%
%增加大负荷
PD(22)=PD(22)*65;
% PD(22)=PD(22)*65;
%%
PG=sparse(PG);

26
subOPF.m
View File

@ -43,7 +43,7 @@ QG(PVi)=QGBal(PVi);
%%
[Volt,UAngel,Init_Z,Init_W,Init_L,Init_U,Init_Y,PG,QG,RestraintCount,wPG,wQG,wPD,wQD,PD,PD00,QD,randPDind,Loadi]=OPF_Init(Busnum,Balance,PG,QG,Volt,GenU,GenL,PVi,PGi,PVQU,PVQL,RealPD,RealQD,QD,PD);
Volt=0.99*Volt;
% Volt=0.99*Volt;
Gap=(Init_L*Init_Z'-Init_U*Init_W');
KK=0;
@ -106,10 +106,10 @@ uQD(noMeasurei)=0.15*mQD(noMeasurei);
%电压伪量测错误
mVolt(noMeasurei)=.95;
mVolt(noMeasurei)=1;
lVolt(noMeasurei)=0.07*mVolt(noMeasurei);%电压在0.93~1.07
uVolt(noMeasurei)=0.07*mVolt(noMeasurei);
lVolt(noMeasurei)=0.12*mVolt(noMeasurei);%电压在0.93~1.07
uVolt(noMeasurei)=0.02*mVolt(noMeasurei);
@ -150,7 +150,7 @@ while(abs(Gap)>Precision*10)
%% 形成方程矩阵
Luu=Init_U'.*Init_W'+Init_u*ones(RestraintCount,1);
Lul=Init_L'.*Init_Z'-Init_u*ones(RestraintCount,1);
bigM=0.5;
bigM=1;
Mat_G=FormG(Volt,PD,QD,Loadi,bigM,mVolt,rVolt,sigma,Vbi,PDbi,QDbi,mPD,mQD,uPD,lPD,uQD,lQD,uVolt,lVolt);
Mat_H=FormH(Busnum,Volt,PG,PD,QG,QD,Y,UAngel,r,c,Angle,Loadi);
Ly=Mat_H;
@ -175,11 +175,11 @@ Gap=1;
% eps=0.00001;
fprintf('\n');
Init_Z=sparse(ones(1,RestraintCount));
Init_W=sparse(-1*ones(1,RestraintCount));
Init_L=1*sparse(ones(1,RestraintCount));
Init_U=1*sparse(ones(1,RestraintCount));
Init_Y=sparse(1,2*Busnum);%与学姐一致
% Init_Z=sparse(ones(1,RestraintCount));
% Init_W=sparse(-1*ones(1,RestraintCount));
% Init_L=1*sparse(ones(1,RestraintCount));
% Init_U=1*sparse(ones(1,RestraintCount));
% Init_Y=sparse(1,2*Busnum);%与学姐一致
fprintf('第2次迭代，算离散量。\n');
while(abs(Gap)>Precision*10)
@ -209,7 +209,7 @@ while(abs(Gap)>Precision*10)
%% 形成方程矩阵
Luu=Init_U'.*Init_W'+Init_u*ones(RestraintCount,1);
Lul=Init_L'.*Init_Z'-Init_u*ones(RestraintCount,1);
bigM=0.5;
bigM=1;
eps=Gap*0.1;
Mat_G=FormG(Volt,PD,QD,Loadi,bigM,mVolt,rVolt,sigma,Vbi,PDbi,QDbi,mPD,mQD,uPD,lPD,uQD,lQD,uVolt,lVolt);
Mat_H=FormH(Busnum,Volt,PG,PD,QG,QD,Y,UAngel,r,c,Angle,Loadi);
@ -234,8 +234,8 @@ dP=PD+diag(Volt)*Y.*cos(AngleIJ)*Volt';
dP(Balance)=0;
dQ=QD+diag(Volt)*Y.*sin(AngleIJ)*Volt';%不为0 ，因为没有考虑电容的QG
dQ(Balance)=0;
mdP=max(dP);
mdQ=max(dQ);
mdP=max(dP)
mdQ=max(dQ)
%判断是否收敛
if min(abs(Vbi-.1))>5e-005
isConverge=0;