facat
/
distributionnetwork-power2current
1
0
Fork 0
Code Issues Pull Requests Releases Wiki Activity

改用误差传递的方差确定方式 

Signed-off-by: dmy@lab <dmy@lab.lab>
This commit is contained in:
dmy@lab 2015-05-08 17:36:59 +08:00
parent 3dd10280c7
commit fa1a29c86b
4 changed files with 80 additions and 103 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

128
IPMLoop.m
View File

@ -30,7 +30,39 @@ isSetBound=0;
largerBound=0;
realBound=1;
tic
while(abs(Gap)>1e-4)
Ir=real(guessIf1(ismember(Loadi,boundedLoadi)) );
pIr=find(Ir>0);
nIr=find(Ir<0);
Ii=imag(guessIf1(ismember(Loadi,boundedLoadi)));
pIi=find(Ii>0);
nIi=find(Ii<0);
lower=ones(length(boundedLoadi)*2,1);
if realBound==1
isSetBound=1;
lower(pIr)=(0.7)*Ir(pIr);
lower(nIr)=(1.300)*Ir(nIr);
lower(pIi+length(Ir))=(0.7)*Ii(pIi);
lower(nIi+length(Ir))=(1.30)*Ii(nIi);
end
Ir=real(guessIf1(ismember(Loadi,boundedLoadi)));
pIr=find(Ir>0);
nIr=find(Ir<0);
Ii=imag(guessIf1(ismember(Loadi,boundedLoadi)));
pIi=find(Ii>0);
nIi=find(Ii<0);
upper=ones(length(boundedLoadi)*2,1);
% if abs(Gap)<1 || isSetBound==1
if realBound==1
isSetBound=1;
upper(pIr)=(1.30)*Ir(pIr);
upper(nIr)=(0.7)*Ir(nIr);
upper(pIi+length(Ir))=(1.30)*Ii(pIi);
upper(nIi+length(Ir))=(0.7)*Ii(nIi);
end
while(abs(Gap)>1e-5)
if KK>=kmax
break;
end
@ -63,99 +95,7 @@ deltG=func_deltGstate1(busNum,Loadi,boundedLoadi,I1r,I1i);
Mat_G=FormGstate1(I1r,I1i,ismember(Loadi,boundedLoadi) );
Mat_H=FormH(fsY1,Loadi,V1r,V1i,I1r,I1i,BalI1r,BalI1i,Balance);
Ly=Mat_H;
% if isSetBound==0
Ir=real(guessIf1(ismember(Loadi,boundedLoadi)) );
pIr=find(Ir>0);
nIr=find(Ir<0);
Ii=imag(guessIf1(ismember(Loadi,boundedLoadi)));
pIi=find(Ii>0);
nIi=find(Ii<0);
% lower=ones(length(Loadi)*2,1);
lower=ones(length(boundedLoadi)*2,1);
% if abs(Gap)<1 || isSetBound==1
% if abs(Gap)<0.001
% largerBound=0;
% realBound=1;
% end
if realBound==1
% lower(pIr)=(0.9-3/(KK+1))*Ir(pIr);
% lower(nIr)=(1.100+3/(KK+1))*Ir(nIr);
% lower(pIi+length(Ir))=(0.9-3/(KK+1))*Ii(pIi);
% lower(nIi+length(Ir))=(1.10+3/(KK+1))*Ii(nIi);
isSetBound=1;
lower(pIr)=(0.7)*Ir(pIr);
lower(nIr)=(1.300)*Ir(nIr);
lower(pIi+length(Ir))=(0.7)*Ii(pIi);
lower(nIi+length(Ir))=(1.30)*Ii(nIi);
% lower( ismember(Loadi,noBoundedLoadi))=-1.9;
% lower(ismember(Loadi,noBoundedLoadi)+length(Ir))=-1.9;
end
if largerBound==1
lower(pIr)=-2.8;
lower(nIr)=-2.8;
lower(pIi+length(Ir))=-2.8;
lower(nIi+length(Ir))=-2.8;
% lower(pIr)=(0.7)*Ir(pIr)-0.1;
% lower(nIr)=(1.300)*Ir(nIr)-0.1;
% lower(pIi+length(Ir))=(0.7)*Ii(pIi)-0.1;
% lower(nIi+length(Ir))=(1.30)*Ii(nIi)-0.1;
end
% lowerR=lower(1:length(lower)/2);
% lowerR(I1r./lowerR>0.98)=lowerR(I1r./lowerR>0.98)-0.001;
% lowerI=lower(length(lower)/2+1:end);
% lowerI(I1i./lowerI>0.98)=lowerI(I1i./lowerI>0.98)-0.001;
% lower=[lowerR;lowerI];
% I1r(I1r./lowerR>0.9998)=lowerR(I1r./lowerR>0.9998);
% I1i(I1i./lowerI>0.9998)=lowerI(I1i./lowerI>0.9998);
Ir=real(guessIf1(ismember(Loadi,boundedLoadi)));
pIr=find(Ir>0);
nIr=find(Ir<0);
Ii=imag(guessIf1(ismember(Loadi,boundedLoadi)));
pIi=find(Ii>0);
nIi=find(Ii<0);
upper=ones(length(boundedLoadi)*2,1);
% if abs(Gap)<1 || isSetBound==1
if realBound==1
% upper(pIr)=(1.10+3/(KK+1))*Ir(pIr);
% upper(nIr)=(0.9-3/(KK+1))*Ir(nIr);
% upper(pIi+length(Ir))=(1.10+3/(KK+1))*Ii(pIi);
% upper(nIi+length(Ir))=(0.9-3/(KK+1))*Ii(nIi);
isSetBound=1;
upper(pIr)=(1.30)*Ir(pIr);
upper(nIr)=(0.7)*Ir(nIr);
upper(pIi+length(Ir))=(1.30)*Ii(pIi);
upper(nIi+length(Ir))=(0.7)*Ii(nIi);
% upper(ismember(Loadi,noBoundedLoadi))=1.9;
% upper(ismember(Loadi,noBoundedLoadi)+length(Ir))=1.9;
end
if largerBound==1
upper(pIr)=2.8;
upper(nIr)=2.8;
upper(pIi+length(Ir))=2.8;
upper(nIi+length(Ir))=2.8;
% upper(pIr)=(1.30)*Ir(pIr)+0.1;
% upper(nIr)=(0.7)*Ir(nIr)+0.1;
% upper(pIi+length(Ir))=(1.30)*Ii(pIi)+0.1;
% upper(nIi+length(Ir))=(0.7)*Ii(nIi)+0.1;
end
% end
Lz=FormLzstate1(Loadi,Mat_G,Init_L,I1measurement,dI_F,flag,guessIf1,rIf1,lower);
Lw=FormLwstate1(Loadi,Mat_G,Init_U,I1measurement,dI_F,flag,guessIf1,rIf1,upper);
Lx=FormLx(deltF,deltH,Init_Y,deltG,Init_Z,Init_W);
@ -173,7 +113,7 @@ end
toc
calcuTime=toc/KK;
% f=sum(([real(I1measurement);imag(I1measurement)]-[I1r;I1i]).^2)+sum((real(V1measurement)-V1r(Loadi)).^2)+sum((imag(V1measurement)-V1i(Loadi)).^2);
if abs(Gap)<1e-4 && KK<kmax
if abs(Gap)<1e-5 && KK<kmax
isConverged=1;
else
isConverged=0;

8
calIWi.m Normal file
View File

@ -0,0 +1,8 @@
function [ output_args ] = calIWi( Vr,Vi,sigmaP,sigmaQ )
dIrdP=(Vi./(Vr.^2+Vi.^2)).^2.*(sigmaP.^2);
dIrdQ=(Vr./(Vr.^2+Vi.^2)).^2.*(sigmaQ.^2);
%dIrdVr=(P.*(Vr.^2+Vi.^2)-2*(P.*Vr+Q.*Vi).*Vr)./(Vr.^2+Vi.^2);
%dIrdVr=(Q.*(Vr.^2+Vi.^2)-2*(P.*Vr+Q.*Vi).*Vi)./(Vr.^2+Vi.^2);
output_args=max(abs([dIrdP,dIrdQ]),[],2);
end

9
calIWr.m Normal file
View File

@ -0,0 +1,9 @@
function [ output_args ] = calIWr(Vr,Vi,sigmaP,sigmaQ )
%计算电流实部的方差
dIrdP=(Vr./(Vr.^2+Vi.^2)).^2.*(sigmaP.^2);
dIrdQ=(Vi./(Vr.^2+Vi.^2)).^2.*(sigmaQ.^2);
%dIrdVr=(P.*(Vr.^2+Vi.^2)-2*(P.*Vr+Q.*Vi).*Vr)./(Vr.^2+Vi.^2);
%dIrdVr=(Q.*(Vr.^2+Vi.^2)-2*(P.*Vr+Q.*Vi).*Vi)./(Vr.^2+Vi.^2);
output_args=max(abs([dIrdP,dIrdQ]),[],2);
end

38
run.m
View File

@ -214,13 +214,21 @@ I1measurement=mIf1(Loadi);
% wI1i=abs( imag(mIf1(Loadi).*sigmaI1) );
%正序测量误差的第二种形式
wI1r=abs(real(I1measurement)).*sigma;
% wI1r=abs(real(I1measurement)).*sigma;
% wI1i=abs(imag(I1measurement)).*sigma;
% dIrdPr=real(V1measurement)./(abs(V1measurement).^2);
% dIrdQr=imag(V1measurement)./(abs(V1measurement).^2);
% dIrdV=
IArW=calIWr(real(mVoltpA),imag(mVoltpA),mphaseASpotLoadP*sigma,mphaseASpotLoadQ*sigma );
IAiW=calIWi(real(mVoltpA),imag(mVoltpA),mphaseASpotLoadP*sigma,mphaseASpotLoadQ*sigma );
IBrW=calIWr(real(mVoltpB),imag(mVoltpB),mphaseBSpotLoadP*sigma,mphaseBSpotLoadQ*sigma );
IBiW=calIWi(real(mVoltpB),imag(mVoltpB),mphaseBSpotLoadP*sigma,mphaseBSpotLoadQ*sigma );
ICrW=calIWr(real(mVoltpC),imag(mVoltpC),mphaseCSpotLoadP*sigma,mphaseCSpotLoadQ*sigma );
ICiW=calIWi(real(mVoltpC),imag(mVoltpC),mphaseCSpotLoadP*sigma,mphaseCSpotLoadQ*sigma );
wI1r=sqrt((IArW+ICrW+IBrW)/9 );
wI1r=wI1r(Loadi);
wI1i=sqrt((IAiW+ICiW+IBiW)/9 );
wI1i=wI1i(Loadi);
wI1i=abs(imag(I1measurement)).*sigma;
wV1r=sqrt(2*( (abs(mVoltpA).*sigma).^2+(abs(mVoltpB).*sigma).^2+(abs(mVoltpC).*sigma).^2 )/9 );
wV1i=sqrt(2*( (abs(mVoltpA).*sigma).^2+(abs(mVoltpB).*sigma).^2+(abs(mVoltpC).*sigma).^2 )/9 );
wV1r=wV1r(Loadi);
@ -249,8 +257,14 @@ I2measurement=mIf2(Loadi);
% wI2r=abs( real(mIf2(Loadi).*sigmaI2) );
% wI2i=abs( imag(mIf2(Loadi).*sigmaI2) );
%负序误差的第二种形式
wI2r=abs(real(I2measurement)).*sigma;
wI2i=abs(imag(I2measurement)).*sigma;
% wI2r=abs(real(I2measurement)).*sigma;
% wI2i=abs(imag(I2measurement)).*sigma;
wI2r=sqrt((IArW+ICrW+IBrW)/9 );
wI2r=wI2r(Loadi);
wI2i=sqrt((IAiW+ICiW+IBiW)/9 );
wI2i=wI2i(Loadi);
wV2r=sqrt(2*( (abs(mVoltpA).*sigma).^2+(abs(mVoltpB).*sigma).^2+(abs(mVoltpC).*sigma).^2 )/9 );
wV2i=sqrt(2*( (abs(mVoltpA).*sigma).^2+(abs(mVoltpB).*sigma).^2+(abs(mVoltpC).*sigma).^2 )/9 );
wV2r=wV2r(Loadi);
@ -280,9 +294,15 @@ I0measurement=mIf0(Loadi);
% wI0r=abs( real(mIf0(Loadi).*sigmaI0) );
% wI0i=abs( imag(mIf0(Loadi).*sigmaI0) );
%零序误差的第二种形式
%
% wI0r=abs(real(I0measurement)).*sigma;
% wI0i=abs(imag(I0measurement)).*sigma;
wI0r=sqrt((IArW+ICrW+IBrW)/9 );
wI0r=wI0r(Loadi);
wI0i=sqrt((IAiW+ICiW+IBiW)/9 );
wI0i=wI0i(Loadi);
wI0r=abs(real(I0measurement)).*sigma;
wI0i=abs(imag(I0measurement)).*sigma;
wV0r=sqrt(2*( (abs(mVoltpA).*sigma).^2+(abs(mVoltpB).*sigma).^2+(abs(mVoltpC).*sigma).^2 )/9 );
wV0i=sqrt(2*( (abs(mVoltpA).*sigma).^2+(abs(mVoltpB).*sigma).^2+(abs(mVoltpC).*sigma).^2 )/9 );
wV0r=wV0r(Loadi);