1.改用新算例。

2.考虑负荷曲线不准确情况。
3.修改统计中的大循环,修复统计错误的问题。

Signed-off-by: dmy@lab <dmy@lab.lab>
This commit is contained in:
dmy@lab 2015-04-06 11:51:26 +08:00
parent a7141734c0
commit 562b209ec8
6 changed files with 116 additions and 91 deletions

View File

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

View File

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

View File

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

View File

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

88
run.m
View File

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

21
test.m
View File

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