From 8f5fe67a6a5df3f57fe5b785add5e4e13dad1de7 Mon Sep 17 00:00:00 2001 From: "dugg@lab-desk" Date: Wed, 11 Feb 2015 14:35:07 +0800 Subject: [PATCH] =?UTF-8?q?=E6=9C=80=E7=BB=88=E9=80=89=E6=8B=A9=E7=94=A8?= =?UTF-8?q?=E8=BF=99=E4=B8=AA=E7=BB=84=E6=95=B0=E6=8D=AE=E5=92=8C=E8=BF=99?= =?UTF-8?q?=E4=B8=AA=E8=AF=84=E4=BB=B7=E6=8C=87=E6=A0=87?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Signed-off-by: dugg@lab-desk --- OPF.m | 105 ++++++++++++++++++++++++++++------------------------ PD0.mat | Bin 534 -> 533 bytes QD0.mat | Bin 530 -> 533 bytes mVolt.mat | Bin 458 -> 455 bytes openfile2.m | 2 +- subOPF.m | 26 ++++++------- 6 files changed, 70 insertions(+), 63 deletions(-) diff --git a/OPF.m b/OPF.m index 26ddf60..f20e22d 100644 --- a/OPF.m +++ b/OPF.m @@ -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','最大不平衡量'); diff --git a/PD0.mat b/PD0.mat index db6d1f31d988fce42933f5328c4155ae4eab4401..b345fc1d8ffe8e9f30614d2bd2b54ca553ad6445 100644 GIT binary patch delta 323 zcmbQnGL>b5kwkcEih^5el7gY3f}x3(v5A#|fr62Nq3Ogx?THEO6KgmadndNm*8lvw zbH)R~-IwiNU-%@#Gu_YuS@RQZ=!_~%p(m4_7_4P+D>TmyY&FZ&wuO!>P{o8J)=gs?c z&M19%ZRJB3w{HI1Z>EcJ{=fKZ$632S2U|N0BX7SAom(!yOLgzvx}SMl?{3Q3{`PVB zXOB;}?rrnG``PS6$!V+KzGuZ_wy%##*H0`9s`)&>^xNsR#unk>?>7H0R%kw4_|)_F lZ{4juo`R9Ff{~$>k%^U&se+M#q3Ogx?THEO6Kgma`zE&5*3a$S zb^p=fjhFw;sW1MM()jk<>}6?fal5ykkJE~*T6g4H>9>FHw?F^8W7D~RPVLtZ&-*g@ zoBpNGF3;cDX6h^a-uCC)n$mB4GOXR}XRh09O0?~mGv@?!4QT z{f_JYm;ILXnZ2#+@BOs*1*In^ar}Gxt*@tTq3X%oZ}0o=`7>)}%I?2M3_h2Y-o5?J z&42rryLP3g>-;B8kN3T9RlT!8wQ%SB>o@Js-`uXVG-dAJKTn@4EeZN|Q-1HNt+pFp pF8QnV`tNuF diff --git a/QD0.mat b/QD0.mat index ed2f708bb001cfea969fcfdf2822844efa27655b..fc088914251904cf39b166b82c19fa19acfd97ba 100644 GIT binary patch delta 323 zcmbQlGL>b5kwkcEih^5el7gY3f}x3(v9Xn@se+M#q3Ogx?THEO6KgmadndNm)<2cF zH}~q!s=2o3Pu*O#JZSn2_O1Fqe*b+WR{inoAI-g$pWod7w&Z=d*ZidQmruq9+SD%p zd;Nmi`6JENZ=b&25$tc_x3NDa=j4hLrU5oUrqe` ml2x5=jm5r;_fCqwFWYnU(cEwB%e&QL-p3~KGq_HuT?qhVz^xMi delta 320 zcmbQrGKpn^k%U`Orh->uo`R9Ff{~$>k%^VDfr62Nq3Ogx?THEO6KgmayC$~Q*2i8l zoo=tT;htUcri5K_zlDBRvn}U)aXWvbeW}BijlZwk{I+#-7rQsh=IFl(uQrB0S9`z3 z{!?Vj&c%7l|L#5UwMVfq_zm~@7Or#a=OnHFym+_8wxW5pzhh>WuL&2s#ocy$E4yBN zM*j6(zV*(3#D1InG+u9YDC*7J)2%^U()>1_uU|POF(-@GR?zlBelzx7W!`>yXNgd5*pH@kIx!_=ufLrM?XSq!KfC?qwhO1rzj>~Y`x&CKXS?{mcq7%+*VezcM_z4AZ@qVO zxyv-!8au`l^;sJ>f4LuB%{+VOEdQVPZ%W&~S-Yq9*IL`(@yrunUdw;WyzM@hdCi|| z|6XsF-NQZk{w@2o&HIcmzF7NrZqJ%^4EMKxORZln^8RG{{;mBL`!^O;|BwDXxhm^t z_M7{+-bdGN%SgIW{d@AHiGlHvb~ip>-F|+D@z1?`kM4cVHC0Of&D~b9cUE`*-MhEq z=Y^LKCVc<)X>N{Xf5W@lbnD7$gA29)5AB_M_h!SpJ@I>9M81vw9XEaT_P>uBpKYHv t|L=!q(Z8Ek@|o|J;=O;Xe53!|3({Y%{R`OfaF?q3{aefot8Vu`0ssivxh4Pr delta 331 zcmX@ke2RI3k%U`Orh->uo`R9Ff{~$>k*SrTfr62Nq3Ogx?THEO6KgmaohG(w)faqQ z#Q(g;)$YdpuAmP!pP%midwBVs(C6#+Z7=_&z3+JN>%FxJ-K&nJ9Q^)myTtu}bJzU2 z_wUnWaWDU8`!|cv&e^N}=~r~M?6FM`g!|sV;h$Gk)3fuPrjA!nsjaN#``G; zKdpcFQ@Z-;ov>D+rTpL8pN2hKy{0<9ev?{9fqsbPi{%y}@{+{v;w(owQi~Ij1@#yNZ_iw`YO{skF;nUi`Ip2~*4{qE3 vZM*5;dGB|9z7{8&vu*K;E|Ckhi~2>BFWPrecision*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;