From 562b209ec89ee010b9eecd55cfc2ef9693867437 Mon Sep 17 00:00:00 2001 From: "dmy@lab" Date: Mon, 6 Apr 2015 11:51:26 +0800 Subject: [PATCH] =?UTF-8?q?1.=E6=94=B9=E7=94=A8=E6=96=B0=E7=AE=97=E4=BE=8B?= =?UTF-8?q?=E3=80=82=202.=E8=80=83=E8=99=91=E8=B4=9F=E8=8D=B7=E6=9B=B2?= =?UTF-8?q?=E7=BA=BF=E4=B8=8D=E5=87=86=E7=A1=AE=E6=83=85=E5=86=B5=E3=80=82?= =?UTF-8?q?=203.=E4=BF=AE=E6=94=B9=E7=BB=9F=E8=AE=A1=E4=B8=AD=E7=9A=84?= =?UTF-8?q?=E5=A4=A7=E5=BE=AA=E7=8E=AF=EF=BC=8C=E4=BF=AE=E5=A4=8D=E7=BB=9F?= =?UTF-8?q?=E8=AE=A1=E9=94=99=E8=AF=AF=E7=9A=84=E9=97=AE=E9=A2=98=E3=80=82?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Signed-off-by: dmy@lab --- FormLwstate1.m | 18 ++++++---- FormLzstate1.m | 25 +++++++++---- IPMLoop.m | 8 ++--- feeder13/data1.txt | 47 ++++++++++++++----------- run.m | 88 ++++++++++++++++++++-------------------------- test.m | 21 ++++++++--- 6 files changed, 116 insertions(+), 91 deletions(-) diff --git a/FormLwstate1.m b/FormLwstate1.m index ae3459d..2788fd7 100644 --- a/FormLwstate1.m +++ b/FormLwstate1.m @@ -1,9 +1,9 @@ -function Lw=FormLwstate1(Loadi,Mat_G,Init_U,I1measurement,dI_F,flag) +function Lw=FormLwstate1(Loadi,Mat_G,Init_U,I1measurement,dI_F,flag,guessIf1) -Ir=real(I1measurement); +Ir=real(guessIf1); pIr=find(Ir>0); nIr=find(Ir<0); -Ii=imag(I1measurement); +Ii=imag(guessIf1); pIi=find(Ii>0); nIi=find(Ii<0); @@ -16,10 +16,14 @@ nIi=find(Ii<0); % upper(pIi+length(Ir))=1.4*Ii(pIi); % upper(nIi+length(Ir))=0.6*Ii(nIi); upper=ones(length(Loadi)*2,1); -upper(pIr)=abs(real(dI_F(pIr,flag)))+Ir(pIr); -upper(nIr)=abs(real(dI_F(nIr,flag)))+Ir(nIr); -upper(pIi+length(Ir))=abs(imag(dI_F(pIi,flag)))+Ii(pIi); -upper(nIi+length(Ir))=abs(imag(dI_F(nIi,flag)))+Ii(nIi); +% upper(pIr)=abs(real(dI_F(pIr,flag)))+Ir(pIr); +% upper(nIr)=abs(real(dI_F(nIr,flag)))+Ir(nIr); +% upper(pIi+length(Ir))=abs(imag(dI_F(pIi,flag)))+Ii(pIi); +% upper(nIi+length(Ir))=abs(imag(dI_F(nIi,flag)))+Ii(nIi); +upper(pIr)=1.3*Ir(pIr); +upper(nIr)=0.7*Ir(nIr); +upper(pIi+length(Ir))=1.3*Ii(pIi); +upper(nIi+length(Ir))=0.7*Ii(nIi); %太小的数都要放宽一些 diff --git a/FormLzstate1.m b/FormLzstate1.m index de28d2a..1886842 100644 --- a/FormLzstate1.m +++ b/FormLzstate1.m @@ -1,8 +1,8 @@ -function Lz=FormLzstate1(Loadi,Mat_G,Init_L,I1measurement,dI_F,flag) -Ir=real(I1measurement); +function Lz=FormLzstate1(Loadi,Mat_G,Init_L,I1measurement,dI_F,flag,guessIf1) +Ir=real(guessIf1); pIr=find(Ir>0); nIr=find(Ir<0); -Ii=imag(I1measurement); +Ii=imag(guessIf1); pIi=find(Ii>0); nIi=find(Ii<0); @@ -14,10 +14,21 @@ lower=ones(length(Loadi)*2,1); % lower(pIi+length(Ir))=0.6*Ii(pIi); % lower(nIi+length(Ir))=1.4*Ii(nIi); -lower(pIr)=-abs(real(dI_F(pIr,flag)))+Ir(pIr); -lower(nIr)=-abs(real(dI_F(nIr,flag)))+Ir(nIr); -lower(pIi+length(Ir))=-abs(imag(dI_F(pIi,flag)))+Ii(pIi); -lower(nIi+length(Ir))=-abs(imag(dI_F(nIi,flag)))+Ii(nIi); +% lower(pIr)=-abs(real(dI_F(pIr,flag)))+Ir(pIr); +% lower(nIr)=-abs(real(dI_F(nIr,flag)))+Ir(nIr); +% lower(pIi+length(Ir))=-abs(imag(dI_F(pIi,flag)))+Ii(pIi); +% lower(nIi+length(Ir))=-abs(imag(dI_F(nIi,flag)))+Ii(nIi); + +% lower(pIr)=-abs(real(dI_F(pIr,flag)))+Ir(pIr); +% lower(nIr)=-abs(real(dI_F(nIr,flag)))+Ir(nIr); +% lower(pIi+length(Ir))=-abs(imag(dI_F(pIi,flag)))+Ii(pIi); +% lower(nIi+length(Ir))=-abs(imag(dI_F(nIi,flag)))+Ii(nIi); + +lower(pIr)=0.7*Ir(pIr); +lower(nIr)=1.3*Ir(nIr); +lower(pIi+length(Ir))=0.7*Ii(pIi); +lower(nIi+length(Ir))=1.3*Ii(nIi); + %太小的数都要放宽一些 % tooSmall=find(abs(Ir)<0.0005); diff --git a/IPMLoop.m b/IPMLoop.m index 8138c96..3bd8053 100644 --- a/IPMLoop.m +++ b/IPMLoop.m @@ -1,4 +1,4 @@ -function [ V1r,V1i,I1r,I1i,isConverged ] = IPMLoop(V1measurement,wV1r,wV1i,I1measurement,wI1r,wI1i,BalI1r,BalI1i,busNum,Loadi,fsY1,Balance,Vref,dI_F,flag ) +function [ V1r,V1i,I1r,I1i,isConverged ] = IPMLoop(V1measurement,wV1r,wV1i,I1measurement,wI1r,wI1i,BalI1r,BalI1i,busNum,Loadi,fsY1,Balance,Vref,dI_F,flag,guessIf1 ) %把每个序的循环写在这个函数中。其实也就是内点法循环。 V1r=Vref*ones(busNum,1); V1i=0*ones(busNum,1); @@ -17,7 +17,7 @@ Init_L=1*sparse(ones(RestraintCount,1)); Init_U=1*sparse(ones(RestraintCount,1)); Init_Y=sparse(2*busNum,1);%等式约束乘子 Gap=(Init_L'*Init_Z-Init_U'*Init_W); -kmax=10; +kmax=30; while(abs(Gap)>0.00001) if KK>=kmax break; @@ -49,8 +49,8 @@ while(abs(Gap)>0.00001) Mat_G=FormGstate1(I1r,I1i); Mat_H=FormH(fsY1,Loadi,V1r,V1i,I1r,I1i,BalI1r,BalI1i,Balance); Ly=Mat_H; - Lz=FormLzstate1(Loadi,Mat_G,Init_L,I1measurement,dI_F,flag); - Lw=FormLwstate1(Loadi,Mat_G,Init_U,I1measurement,dI_F,flag); + Lz=FormLzstate1(Loadi,Mat_G,Init_L,I1measurement,dI_F,flag,guessIf1); + Lw=FormLwstate1(Loadi,Mat_G,Init_U,I1measurement,dI_F,flag,guessIf1); Lx=FormLx(deltF,deltH,Init_Y,deltG,Init_Z,Init_W); YY=FormYY(Lul,Lz,Ly,Luu,Lw,Lx); %% 开始解方程 diff --git a/feeder13/data1.txt b/feeder13/data1.txt index e7c35c7..27e0b6e 100644 --- a/feeder13/data1.txt +++ b/feeder13/data1.txt @@ -1,26 +1,27 @@ 650 4.16 0 -600 632 645 500 -602 632 633 500 -601 633 634 600 -600 645 646 300 -601 650 632 2000 -602 684 652 800 -601 632 671 2000 -600 671 684 300 -602 671 680 1000 -600 671 692 300 -602 684 611 300 -601 692 675 500 +600 632 645 50 +602 632 633 50 +601 633 634 60 +600 645 646 30 +601 650 632 200 +602 684 652 80 +601 632 671 200 +600 671 684 30 +600 671 692 30 +602 684 611 30 +601 692 675 50 0 -632 0 0 0 +692 0 0 0 0 -634 160 110 120 90 120 90 -645 123 102 170 125 133 110 -646 172 150 230 132 211 160 -652 128 86 30 10 80 70 -675 85 10 68 60 20 12 -611 52 42 85 60 70 40 +634 120 90 120 90 160 110 +645 40 30 50 40 60 50 +646 100 50 70 40 60 40 +652 20 10 60 30 40 40 +671 385 220 385 220 385 220 +675 68 60 290 212 285 190 +692 260 200 40 40 70 40 +611 80 20 60 40 30 20 0 600 632 645 500 602 632 633 500 @@ -46,5 +47,11 @@ 692 142 140 153 121 170 151 611 156 142 185 60 170 80 0 - +634 160 110 120 90 120 90 +645 123 102 170 125 133 110 +646 17 15 23 13 21 16 +652 128 86 30 10 80 70 +675 85 10 68 60 20 12 +611 52 42 85 60 70 40 +680 12 4 10 6 7 4 diff --git a/run.m b/run.m index 726a2d9..c09945a 100644 --- a/run.m +++ b/run.m @@ -1,4 +1,4 @@ -function [JMeasurement,AME_Volt,AME_VAngle,AME_PD,AME_QD,isConverged]=run() +function [JMeasurement,AME_Volt,AME_VAngle,AME_PD,AME_QD,AME_mVolt,AME_mPD,AME_mQD,isConverged]=run() %% 利用先把负荷转换为电流的方法。这个方法要求知道电压量。 % close all @@ -32,7 +32,7 @@ PD=Pabc/3; QD=Qabc/3; Loadi=find(PD~=0); maxD=100000;% 最大不平衡量 -EPS=1e-5; +EPS=1e-6; k=0; kmax=20; fsY11=fsY1; @@ -142,13 +142,13 @@ fprintf(' rIf0=If0; rIf1=If1; rIf2=If2; -sigma=0.01; -iterPhaseASpotLoadP=phaseASpotLoadP; -iterPhaseBSpotLoadP=phaseBSpotLoadP; -iterPhaseCSpotLoadP=phaseCSpotLoadP; -iterPhaseASpotLoadQ=phaseASpotLoadQ; -iterPhaseBSpotLoadQ=phaseBSpotLoadQ; -iterPhaseCSpotLoadQ=phaseCSpotLoadQ; +sigma=0.03; +% iterPhaseASpotLoadP=phaseASpotLoadP; +% iterPhaseBSpotLoadP=phaseBSpotLoadP; +% iterPhaseCSpotLoadP=phaseCSpotLoadP; +% iterPhaseASpotLoadQ=phaseASpotLoadQ; +% iterPhaseBSpotLoadQ=phaseBSpotLoadQ; +% iterPhaseCSpotLoadQ=phaseCSpotLoadQ; %加误差 mphaseASpotLoadP=phaseASpotLoadP.*(1+normrnd(0,sigma,length(phaseASpotLoadP),1)); mphaseBSpotLoadP=phaseBSpotLoadP.*(1+normrnd(0,sigma,length(phaseBSpotLoadP),1)); @@ -157,13 +157,13 @@ mphaseASpotLoadQ=phaseASpotLoadQ.*(1+normrnd(0,sigma,length(phaseASpotLoadQ),1)) mphaseBSpotLoadQ=phaseBSpotLoadQ.*(1+normrnd(0,sigma,length(phaseBSpotLoadQ),1)); mphaseCSpotLoadQ=phaseCSpotLoadQ.*(1+normrnd(0,sigma,length(phaseCSpotLoadQ),1)); %电压量测量 -mVoltpA=VoltpA.*(1+normrnd(0,sigma,length(VoltpA),1)); -mVoltpB=VoltpB.*(1+normrnd(0,sigma,length(VoltpB),1)); -mVoltpC=VoltpC.*(1+normrnd(0,sigma,length(VoltpC),1)); +% mVoltpA=VoltpA.*(1+normrnd(0,sigma,length(VoltpA),1)); +% mVoltpB=VoltpB.*(1+normrnd(0,sigma,length(VoltpB),1)); +% mVoltpC=VoltpC.*(1+normrnd(0,sigma,length(VoltpC),1)); % -mVoltpA=abs(VoltpA).*(1+normrnd(0,sigma,length(VoltpA),1)); -mVoltpB=abs(VoltpB).*(1+normrnd(0,sigma,length(VoltpB),1)).*exp(1j*-120/180*pi); -mVoltpC=abs(VoltpC).*(1+normrnd(0,sigma,length(VoltpC),1)).*exp(1j*+120/180*pi); +mVoltpA=abs(rVoltpA).*(1+normrnd(0,sigma,length(VoltpA),1)); +mVoltpB=abs(rVoltpB).*(1+normrnd(0,sigma,length(VoltpB),1)).*exp(1j*-120/180*pi); +mVoltpC=abs(rVoltpC).*(1+normrnd(0,sigma,length(VoltpC),1)).*exp(1j*+120/180*pi); % % mVoltpA=sparse(ones(busNum,1)); % mVoltpB=sparse(ones(busNum,1)).*exp(1j*-120/180*pi); @@ -328,42 +328,24 @@ wV0i=wV0i(Loadi); % dI_F=conj(dI_F)'; % dI_F=dI_F(Loadi,:); % 换一个方式来计算 -deltmCurpA=conj((mphaseASpotLoadP+1j*mphaseASpotLoadQ)./mVoltpA)*0.2; -deltmCurpB=conj((mphaseBSpotLoadP+1j*mphaseBSpotLoadQ)./mVoltpB)*0.2; -deltmCurpC=conj((mphaseCSpotLoadP+1j*mphaseCSpotLoadQ)./mVoltpC)*0.2; +deltmCurpA=conj((mphaseASpotLoadP+1j*mphaseASpotLoadQ)./mVoltpA)*0.3; +deltmCurpB=conj((mphaseBSpotLoadP+1j*mphaseBSpotLoadQ)./mVoltpB)*0.3; +deltmCurpC=conj((mphaseCSpotLoadP+1j*mphaseCSpotLoadQ)./mVoltpC)*0.3; %转换为序电流 deltF012=Tp2f*conj([deltmCurpA';deltmCurpB';deltmCurpC']); dI_F=conj(deltF012)'; dI_F=dI_F(Loadi,:); - -% deltmCurpA1=conj((mphaseASpotLoadP+1j*mphaseASpotLoadQ)./mVoltpA)*0.8; -% deltmCurpB1=conj((mphaseBSpotLoadP+1j*mphaseBSpotLoadQ)./mVoltpB)*0.8; -% deltmCurpC1=conj((mphaseCSpotLoadP+1j*mphaseCSpotLoadQ)./mVoltpC)*0.8; -% -% deltF012=Tp2f*conj([-deltmCurpA1'+mCurpA';-deltmCurpB1'+mCurpB';-deltmCurpC1'+mCurpC']); -% dI_F=conj(deltF012)'; -% dI_F=dI_F(Loadi,:); -% -% -% deltmCurpA1=conj((mphaseASpotLoadP+1j*mphaseASpotLoadQ)./mVoltpA)*1.2; -% deltmCurpB1=conj((mphaseBSpotLoadP+1j*mphaseBSpotLoadQ)./mVoltpB)*1.2; -% deltmCurpC1=conj((mphaseCSpotLoadP+1j*mphaseCSpotLoadQ)./mVoltpC)*1.2; -% -% deltF012=Tp2f*conj([deltmCurpA1'-mCurpA';deltmCurpB1'-mCurpB';deltmCurpC1'-mCurpC']); -% dI_F=conj(deltF012)'; -% dI_F=dI_F(Loadi,:); - %考虑负荷曲线不准确的情况 - -% curveIf1=rIf1.*(1+unifrnd(-0.2,0.2,length(rIf1),1)); -% curveIf2=rIf2.*(1+unifrnd(-0.2,0.2,length(rIf2),1)); -% curveIf0=rIf0.*(1+unifrnd(-0.2,0.2,length(rIf0),1)); -% -% curveCurpA=conj((phaseASpotLoadP.*(1+unifrnd(-0.2,0.2,length(phaseASpotLoadP),1))+1j*phaseASpotLoadQ.*(1+unifrnd(-0.2,0.2,length(phaseASpotLoadQ),1)))./rVoltpA)*0.2; -% curveCurpB=conj((phaseBSpotLoadP.*(1+unifrnd(-0.2,0.2,length(phaseBSpotLoadP),1))+1j*phaseBSpotLoadQ.*(1+unifrnd(-0.2,0.2,length(phaseBSpotLoadQ),1)))./rVoltpB)*0.2; -% curveCurpC=conj((phaseCSpotLoadP.*(1+unifrnd(-0.2,0.2,length(phaseCSpotLoadP),1))+1j*phaseCSpotLoadQ.*(1+unifrnd(-0.2,0.2,length(phaseCSpotLoadQ),1)))./rVoltpC)*0.2; -% deltF012=Tp2f*conj([curveCurpA';curveCurpB';curveCurpC']); - +guessCurpA=conj((phaseASpotLoadP+1j*phaseASpotLoadQ)./mVoltpA).*(1+unifrnd(-0.15,0.15,length(phaseASpotLoadQ),1)); +guessCurpB=conj((phaseBSpotLoadP+1j*phaseBSpotLoadQ)./mVoltpB).*(1+unifrnd(-0.15,0.15,length(phaseASpotLoadQ),1)); +guessCurpC=conj((phaseCSpotLoadP+1j*phaseCSpotLoadQ)./mVoltpC).*(1+unifrnd(-0.15,0.15,length(phaseASpotLoadQ),1)); +guessIf012=Tp2f*conj([-guessCurpA';-guessCurpB';-guessCurpC']); +guessIf0=conj(guessIf012(1,:)'); +guessIf1=conj(guessIf012(2,:)'); +guessIf2=conj(guessIf012(3,:)'); +guessIf0=guessIf0(Loadi); +guessIf1=guessIf1(Loadi); +guessIf2=guessIf2(Loadi); tic for II=1:3 @@ -371,19 +353,19 @@ for II=1:3 if II==1 fprintf('正序\n'); tic - [ V1r,V1i,I1r,I1i,isConverged1 ]=IPMLoop(V1measurement,wV1r,wV1i,I1measurement,wI1r,wI1i,BalI1r,BalI1i,busNum,Loadi,fsY11,Balance,1,dI_F,2 ); + [ V1r,V1i,I1r,I1i,isConverged1 ]=IPMLoop(V1measurement,wV1r,wV1i,I1measurement,wI1r,wI1i,BalI1r,BalI1i,busNum,Loadi,fsY11,Balance,1,dI_F,2,guessIf1 ); toc end if II==2 fprintf('负序\n'); tic - [ V2r,V2i,I2r,I2i,isConverged2 ]=IPMLoop(V2measurement,wV2r,wV2i,I2measurement,wI2r,wI2i,BalI2r,BalI2i,busNum,Loadi,fsY22,Balance,0,dI_F,3 ); + [ V2r,V2i,I2r,I2i,isConverged2 ]=IPMLoop(V2measurement,wV2r,wV2i,I2measurement,wI2r,wI2i,BalI2r,BalI2i,busNum,Loadi,fsY22,Balance,0,dI_F,3,guessIf2 ); toc end if II==3 fprintf('零序\n'); tic - [ V0r,V0i,I0r,I0i,isConverged0 ]=IPMLoop(V0measurement,wV0r,wV0i,I0measurement,wI0r,wI0i,BalI0r,BalI0i,busNum,Loadi,fsY00,Balance,0,dI_F,1 ); + [ V0r,V0i,I0r,I0i,isConverged0 ]=IPMLoop(V0measurement,wV0r,wV0i,I0measurement,wI0r,wI0i,BalI0r,BalI0i,busNum,Loadi,fsY00,Balance,0,dI_F,1,guessIf0 ); toc end @@ -414,6 +396,10 @@ rThreeLoad=[ phaseASpotLoadP'+1j*phaseASpotLoadQ'; phaseBSpotLoadP'+1j*phaseBSpotLoadQ'; phaseCSpotLoadP'+1j*phaseCSpotLoadQ'; ]; +mThreeLoad=[ mphaseASpotLoadP'+1j*mphaseASpotLoadQ'; + mphaseBSpotLoadP'+1j*mphaseBSpotLoadQ'; + mphaseCSpotLoadP'+1j*mphaseCSpotLoadQ'; + ]; % rThreeLoad=rThreeLoad(:,setxor(1:size(SEVoltpABC,2),Balance)); SEThreeLoad=SEVoltpABC(:,Loadi).*conj(-SEIpABC); % phaseLoadPError=real(rThreeLoad-SEThreeLoad)./real(rThreeLoad)*100; @@ -460,6 +446,10 @@ AME_Volt=sum(sum(abs( abs(rVoltABCV)-abs(SEVoltpABC)))); AME_VAngle=sum(sum(abs( angle(rVoltABCV)-angle(SEVoltpABC)))); AME_PD=sum(sum(abs(real(SEThreeLoad-rThreeLoad(:,Loadi))))); AME_QD=sum(sum(abs(imag(SEThreeLoad-rThreeLoad(:,Loadi))))); +%计算与量测值的 +AME_mVolt=sum(sum(abs( abs(mVoltABCV)-abs(rVoltABCV)))); +AME_mPD=sum(sum(abs(real(rThreeLoad(:,Loadi)-mThreeLoad(:,Loadi))))); +AME_mQD=sum(sum(abs(imag(rThreeLoad(:,Loadi)-mThreeLoad(:,Loadi))))); %返回收敛信息 isConverged=isConverged1*isConverged2*isConverged0; end \ No newline at end of file diff --git a/test.m b/test.m index b3c5027..26b41c1 100644 --- a/test.m +++ b/test.m @@ -6,21 +6,34 @@ AME_VAngleSum=0; JMeasurementSum=0; AME_PDSum=0; AME_QDSum=0; -N=1000; -for I=1:N - [JMeasurement,AME_Volt,AME_VAngle,AME_PD,AME_QD,isConverged]=run(); +AME_mVoltSum=0; +AME_mPDSum=0; +AME_mQDSum=0; +N=200; +loopN=1; +while 1 + [JMeasurement,AME_Volt,AME_VAngle,AME_PD,AME_QD,AME_mVolt,AME_mPD,AME_mQD,isConverged]=run(); if isConverged==0 - I=I-1; continue end + loopN=loopN+1; + if loopN>N + break; + end % [JMeasurement]=run(); AME_VoltSum=AME_VoltSum+AME_Volt; AME_VAngleSum=AME_VAngleSum+AME_VAngle; AME_PDSum=AME_PDSum+AME_PD; AME_QDSum=AME_QDSum+AME_QD; + AME_mVoltSum=AME_mVoltSum+AME_mVolt; + AME_mPDSum=AME_mPDSum+AME_mPD; + AME_mQDSum=AME_mQDSum+AME_mQD; end JMeasurementSum=JMeasurementSum/N; AME_VoltSum=AME_VoltSum/N; AME_VAngleSum=AME_VAngleSum/N; AME_PDSum=AME_PDSum/N; AME_QDSum=AME_QDSum/N; +AME_mVoltSum=AME_mVoltSum/N; +AME_mPDSum=AME_mPDSum/N; +AME_mQDSum=AME_mQDSum/N; \ No newline at end of file