From a091d9f502def92d3d9f2f192dfe24ad5cec5006 Mon Sep 17 00:00:00 2001 From: "dugg@lab-desk" Date: Tue, 3 Feb 2015 21:52:51 +0800 Subject: [PATCH] =?UTF-8?q?1.=E5=8A=A0=E4=BA=86=E4=B8=80=E4=B8=AA=E7=94=BB?= =?UTF-8?q?=E9=94=99=E8=AF=AF=E6=95=B0=E6=8D=AE=E7=9A=84=E5=9B=BE=202.?= =?UTF-8?q?=E5=BC=80=E5=A7=8B=E4=B8=80=E7=82=B9=E7=82=B9=E5=81=9A=E6=B5=8B?= =?UTF-8?q?=E8=AF=95?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Signed-off-by: dugg@lab-desk --- OPF.m | 72 ++++++++++++++++++++++++++++---------------- graph/badDataGraph.m | 16 ++++++++++ subOPF.m | 26 +++------------- 3 files changed, 67 insertions(+), 47 deletions(-) create mode 100644 graph/badDataGraph.m diff --git a/OPF.m b/OPF.m index fef5845..28913b0 100644 --- a/OPF.m +++ b/OPF.m @@ -5,34 +5,54 @@ arrayA=zeros(21,10); sumCaseA_SE=0; sumCaseB_SE=0; sumCaseC_SE=0; -for I=1:1 -close -[kmax,Precision,UAngel,Volt,Busnum,PVi,PVu,Balance,Y,Angle,P0,Q0,r,c,GB,Linei,Linej,Transfori,Transforj,GenU,GenL,GenC,PG,QG,PD,QD,CenterA,PGi,PVQU,PVQL,Liner,Linex,Lineb,Transforr,Transforx,Branchi,Branchg,Branchb,Transfork0]= ... - pf('E:\算例\feeder33\feeder33.txt'); -sigma=0.01; -RealPD=PD; -RealQD=QD; -rVolt=Volt; -Loadi=find(PD~=0); -PD0=sparse(Busnum,1); -QD0=sparse(Busnum,1); -PD0(Loadi)=RealPD(Loadi).*(1+normrnd(0,sigma,length(Loadi),1)); -QD0(Loadi)=RealQD(Loadi).*(1+normrnd(0,sigma,length(Loadi),1)); -mVolt=rVolt.*(1+normrnd(0,sigma,length(rVolt),1))'; +badDataResult=zeros(10,33); +badDataLocation=zeros(34,10); +maxDVolt=0; +maxDVAngle=0; +% for badDataNode=1:33 + for I=1:5 + close + [kmax,Precision,UAngel,Volt,Busnum,PVi,PVu,Balance,Y,Angle,P0,Q0,r,c,GB,Linei,Linej,Transfori,Transforj,GenU,GenL,GenC,PG,QG,PD,QD,CenterA,PGi,PVQU,PVQL,Liner,Linex,Lineb,Transforr,Transforx,Branchi,Branchg,Branchb,Transfork0]= ... + pf('E:\算例\feeder33\feeder33.txt'); + sigma=0.01; + RealPD=PD; + RealQD=QD; + rVolt=Volt; + Loadi=find(PD~=0); + PD0=sparse(Busnum,1); + QD0=sparse(Busnum,1); + PD0(Loadi)=RealPD(Loadi).*(1+normrnd(0,sigma,length(Loadi),1)); + QD0(Loadi)=RealQD(Loadi).*(1+normrnd(0,sigma,length(Loadi),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; + + mVolt(18)=rVolt(18)*(1-sigma*6); + %% 画Case A的图 + % figure('Color',[1 1 1]); + [Busnum,Loadi,Volt,PD,QD,rVolt,UAngel,RealPD,RealQD,rUAngel,Vbi,PDbi,QDbi,plotGapA]=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)))); -%加载保存的变量 -% PD0=load('PD0'); -% PD0=PD0.PD0; -% QD0=load('QD0'); -% QD0=QD0.QD0; -% mVolt=load('mVolt'); -% mVolt=mVolt.mVolt; - -mVolt(2)=rVolt(2)*(1+sigma*8); -%% 画Case A的图 -% figure('Color',[1 1 1]); -[Busnum,Loadi,Volt,PD,QD,rVolt,UAngel,RealPD,RealQD,rUAngel,Vbi,PDbi,QDbi,plotGapA]=subOPF([],PD0,QD0,mVolt,sigma);%全部有 +% 加最大偏差 +maxDVolt_=max(abs((rVolt-Volt))); +if maxDVolt_>maxDVolt + maxDVolt=maxDVolt_; end +maxDVAngle_=max(abs((UAngel(2:33)-rUAngel(2:33)))); +if maxDVAngle_>maxDVAngle + maxDVAngle=maxDVAngle_; +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)*100;%Relative Error of Voltage in Case A % CaseAREV=CaseAREV(2:end)*100; diff --git a/graph/badDataGraph.m b/graph/badDataGraph.m new file mode 100644 index 0000000..82f0da9 --- /dev/null +++ b/graph/badDataGraph.m @@ -0,0 +1,16 @@ +x=1:33; +y=1:33; +x=[x,15,33]; +y=[y,17,23]; +scatter(x,y); +xlabel('存在坏数据的节点'); +ylabel('辨识出坏数据的点'); +hold on; +for I=1:33 +% plot([1,1,1,1]-1+I,[[1:3]-2+I,I]); +plot([I,I],[0,I]); +end +% plot([1,1,1,1]-1+15,[1:3]-2+17); +% plot([1,1,1,1]-1+33,[1,2,12]-2+23); +plot([15,15],[0,17]); +plot([33,33],[0,23]); \ No newline at end of file diff --git a/subOPF.m b/subOPF.m index 2ded212..25e78c7 100644 --- a/subOPF.m +++ b/subOPF.m @@ -45,7 +45,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=100; +kmax=500; Precision=Precision/1; %% 加误差 %找DG @@ -124,7 +124,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=1; + bigM=0.6; 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; @@ -150,15 +150,9 @@ Gap=1; fprintf('\n'); fprintf('第2次迭代，算离散量。\n'); -% while eps>0.000001 -% 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);%与学姐一致 while(abs(Gap)>Precision*10) if KK>kmax - % break; + break; end plotGap(KK+1)=Gap; Init_u=Gap/2/RestraintCount*CenterA; @@ -182,7 +176,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.7; + bigM=0.4; 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); @@ -192,24 +186,14 @@ while(abs(Gap)>Precision*10) Lx=FormLx(deltF,deltH,Init_Y,deltG,Init_Z,Init_W); YY=FormYY(Lul,Lz,Ly,Luu,Lw,Lx); %% 开始解方程 +% Vbi(Vbi>0.003)=.1; fprintf('迭代次数 %d Gap %f\n',KK+1,plotGap(KK+1)); XX=SolveIt(deltF,deltG,Init_L,Init_Z,Init_U,Init_W,deltdeltF,ddh,ddgzw,deltH,Init_Y,Ly,Lz,ContrlCount,Lw,Lul,Luu,RestraintCount,Lx,Balance,PVi,PGi,Busnum,Loadi); - %%取各分量 - % Vbi_=Vbi; - % PDbi_=PDbi; - % QDbi_=QDbi; - % Vbi=1./(1+exp( (0.06-Vbi) /0.01))*0.1; [deltZ,deltL,deltW,deltU,deltX,deltY]=AssignXX(XX,ContrlCount,RestraintCount,Busnum); [Init_Z,Init_L,Init_W,Init_U,Init_Y,PG,QG,Volt,UAngel,PD,QD,Vbi,PDbi,QDbi]=Modification(Init_Z,Init_L,Init_W,Init_U,Init_Y,deltZ,deltL,deltW,deltU,deltX,deltY,PG,QG,Volt,UAngel,PVi,ContrlCount,Balance,Busnum,PGi,PD,QD,Loadi,Vbi,PDbi,QDbi); Gap=(Init_L*Init_Z'-Init_U*Init_W'); - - % Vbi(Vbi>0.01)=1; KK=KK+1; end - -% eps=Gap -% Gap=100; -% end %% 计算最大不平衡量 AngleIJ=sparse(r,c,UAngel(r)-UAngel(c)-Angle',Busnum,Busnum); dP=PD+diag(Volt)*Y.*cos(AngleIJ)*Volt';