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-21 14:50:06 +08:00
parent 6b78384380
commit fddb3316bd
3 changed files with 59 additions and 39 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

10
IPMLoop.m
View File

@ -30,7 +30,7 @@ kmax=15;
isSetBound=0; isSetBound=0;
largerBound=0; largerBound=0;
realBound=1; realBound=1;
tic
Ir=real(guessIf1(ismember(Loadi,boundedLoadi)) ); Ir=real(guessIf1(ismember(Loadi,boundedLoadi)) );
pIr=find(Ir>0); pIr=find(Ir>0);
@ -62,7 +62,7 @@ if realBound==1
upper(pIi+length(Ir))=(1.30)*Ii(pIi); upper(pIi+length(Ir))=(1.30)*Ii(pIi);
upper(nIi+length(Ir))=(0.7)*Ii(nIi); upper(nIi+length(Ir))=(0.7)*Ii(nIi);
end end
tic
while(abs(Gap)>1e-5) while(abs(Gap)>1e-5)
if KK>=kmax if KK>=kmax
break; break;
@ -108,11 +108,11 @@ deltG=func_deltGstate1(busNum,Loadi,boundedLoadi,I1r,I1i);
[Init_Z,Init_L,Init_W,Init_U,Init_Y,V1r,V1i,I1r,I1i]=Modification(Init_Z,Init_L,Init_W,Init_U,Init_Y,deltZ,deltL,deltW,deltU,deltX,deltY,V1r,V1i,I1r,I1i,ContrlCount,Balance,busNum,Loadi,Vref); [Init_Z,Init_L,Init_W,Init_U,Init_Y,V1r,V1i,I1r,I1i]=Modification(Init_Z,Init_L,Init_W,Init_U,Init_Y,deltZ,deltL,deltW,deltU,deltX,deltY,V1r,V1i,I1r,I1i,ContrlCount,Balance,busNum,Loadi,Vref);
Gap=(Init_L'*Init_Z-Init_U'*Init_W);+max(abs(deltX)); Gap=(Init_L'*Init_Z-Init_U'*Init_W);+max(abs(deltX));
% Gap=max([max(abs(Mat_H)),max(Lz),max(Lw) ]); % Gap=max([max(abs(Mat_H)),max(Lz),max(Lw) ]);
% fprintf('Gap %f :%d\n',full(Gap),KK+1); fprintf('Gap %f :%d\n',full(Gap),KK+1);
KK=KK+1; KK=KK+1;
end end
toc % toc
calcuTime=toc/KK; calcuTime=toc;
% f=sum(([real(I1measurement);imag(I1measurement)]-[I1r;I1i]).^2)+sum((real(V1measurement)-V1r(Loadi)).^2)+sum((imag(V1measurement)-V1i(Loadi)).^2); % 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-5 && KK<kmax if abs(Gap)<1e-5 && KK<kmax
isConverged=1; isConverged=1;

73
run.m
View File

@ -1,4 +1,4 @@
function [JMeasurement,AME_Volt,AME_VAngle,AME_PD,AME_QD,AME_mVolt,AME_mPD,AME_mQD,isConverged,totalTime1,totalTime2,totalTime0]=run() function [JMeasurement,AME_Volt,AME_VAngle,AME_PD,AME_QD,AME_mVolt,AME_mPD,AME_mQD,isConverged,totalTime1,totalTime2,totalTime0,eV2,eV0,eI2,eI0]=run()
%% 利用先把负荷转换为电流的方法。这个方法要求知道电压量。 %% 利用先把负荷转换为电流的方法。这个方法要求知道电压量。
% %
close all close all
@ -9,36 +9,36 @@ lineZ=readLineZ('modified-feeder69\lineParameter.txt');
phaseASpotLoadQ,phaseBSpotLoadQ,phaseCSpotLoadQ,setIJ,nodeNum,Balance,phaseABCY ... phaseASpotLoadQ,phaseBSpotLoadQ,phaseCSpotLoadQ,setIJ,nodeNum,Balance,phaseABCY ...
cap]=dataRead(lineZ,'modified-feeder69\data.txt'); cap]=dataRead(lineZ,'modified-feeder69\data.txt');
%负荷不平衡 %负荷不平衡
% for I=1:length(phaseASpotLoadP) for I=1:length(phaseASpotLoadP)
% roll=mod(round(rand()*10),3); roll=mod(round(rand()*10),3);
% if roll==0 if roll==0
% phaseBSpotLoadP(I)=phaseBSpotLoadP(I)+phaseASpotLoadP(I)*.45; phaseBSpotLoadP(I)=phaseBSpotLoadP(I)+phaseASpotLoadP(I)*.45;
% phaseCSpotLoadP(I)=phaseCSpotLoadP(I)+phaseASpotLoadP(I)*.45; phaseCSpotLoadP(I)=phaseCSpotLoadP(I)+phaseASpotLoadP(I)*.45;
% phaseASpotLoadP(I)=phaseASpotLoadP(I)*.10; phaseASpotLoadP(I)=phaseASpotLoadP(I)*.10;
%
% phaseBSpotLoadQ(I)=phaseBSpotLoadQ(I)+phaseASpotLoadQ(I)*.45; phaseBSpotLoadQ(I)=phaseBSpotLoadQ(I)+phaseASpotLoadQ(I)*.45;
% phaseCSpotLoadQ(I)=phaseCSpotLoadQ(I)+phaseASpotLoadQ(I)*.45; phaseCSpotLoadQ(I)=phaseCSpotLoadQ(I)+phaseASpotLoadQ(I)*.45;
% phaseASpotLoadQ(I)=phaseASpotLoadQ(I)*.10; phaseASpotLoadQ(I)=phaseASpotLoadQ(I)*.10;
% end end
% if roll==1 if roll==1
% phaseASpotLoadP(I)=phaseASpotLoadP(I)+phaseBSpotLoadP(I)*.45; phaseASpotLoadP(I)=phaseASpotLoadP(I)+phaseBSpotLoadP(I)*.45;
% phaseCSpotLoadP(I)=phaseCSpotLoadP(I)+phaseBSpotLoadP(I)*.45; phaseCSpotLoadP(I)=phaseCSpotLoadP(I)+phaseBSpotLoadP(I)*.45;
% phaseBSpotLoadP(I)=phaseBSpotLoadP(I)*.10; phaseBSpotLoadP(I)=phaseBSpotLoadP(I)*.10;
%
% phaseASpotLoadQ(I)=phaseASpotLoadQ(I)+phaseBSpotLoadQ(I)*.45; phaseASpotLoadQ(I)=phaseASpotLoadQ(I)+phaseBSpotLoadQ(I)*.45;
% phaseCSpotLoadQ(I)=phaseCSpotLoadQ(I)+phaseBSpotLoadQ(I)*.45; phaseCSpotLoadQ(I)=phaseCSpotLoadQ(I)+phaseBSpotLoadQ(I)*.45;
% phaseBSpotLoadQ(I)=phaseBSpotLoadQ(I)*.10; phaseBSpotLoadQ(I)=phaseBSpotLoadQ(I)*.10;
% end end
% if roll==2 if roll==2
% phaseASpotLoadP(I)=phaseASpotLoadP(I)+phaseCSpotLoadP(I)*.45; phaseASpotLoadP(I)=phaseASpotLoadP(I)+phaseCSpotLoadP(I)*.45;
% phaseBSpotLoadP(I)=phaseBSpotLoadP(I)+phaseCSpotLoadP(I)*.45; phaseBSpotLoadP(I)=phaseBSpotLoadP(I)+phaseCSpotLoadP(I)*.45;
% phaseCSpotLoadP(I)=phaseCSpotLoadP(I)*.10; phaseCSpotLoadP(I)=phaseCSpotLoadP(I)*.10;
%
% phaseASpotLoadQ(I)=phaseASpotLoadQ(I)+phaseCSpotLoadQ(I)*.45; phaseASpotLoadQ(I)=phaseASpotLoadQ(I)+phaseCSpotLoadQ(I)*.45;
% phaseBSpotLoadQ(I)=phaseBSpotLoadQ(I)+phaseCSpotLoadQ(I)*.45; phaseBSpotLoadQ(I)=phaseBSpotLoadQ(I)+phaseCSpotLoadQ(I)*.45;
% phaseCSpotLoadQ(I)=phaseCSpotLoadQ(I)*.10; phaseCSpotLoadQ(I)=phaseCSpotLoadQ(I)*.10;
% end end
% end end
a=exp(1j*2*pi/3); a=exp(1j*2*pi/3);
Tp2f=1/3*[1 1 1; Tp2f=1/3*[1 1 1;
@ -150,7 +150,14 @@ end
FortiscueToc=toc; FortiscueToc=toc;
fprintf('Fortiscue法计算时间 %f\n',FortiscueToc); fprintf('Fortiscue法计算时间 %f\n',FortiscueToc);
Vf1=Vmf1.*exp(1j*Vaf1); Vf1=Vmf1.*exp(1j*Vaf1);
%% %% 计算三相不平衡度
%电压不平衡度
eV2=mean(abs(Vf2)./abs(Vf1))/length(Vf1)/3;
eV0=mean(abs(Vf0)./abs(Vf1))/length(Vf1)/3;
%电流不平衡度
eI2=mean(abs(If2(find(If2)))./abs(If1(find(If1))))/length(find(If1))/3;
eI0=mean(abs(If0(find(If0)))./abs(If1(find(If1))))/length(find(If1))/3;
%%
(Vf0.*conj(fsY00*Vf0)+Vf1.*conj(fsY11*Vf1)+Vf2.*conj(fsY22*Vf2))*3;%包含补偿电容的功率 (Vf0.*conj(fsY00*Vf0)+Vf1.*conj(fsY11*Vf1)+Vf2.*conj(fsY22*Vf2))*3;%包含补偿电容的功率
conj(Tf2p*[If0(2);If1(2);If2(2)]).*(Tf2p*[Vf0(2);Vf1(2);Vf2(2)]); conj(Tf2p*[If0(2);If1(2);If2(2)]).*(Tf2p*[Vf0(2);Vf1(2);Vf2(2)]);
IpABC=Tf2p*conj([If0';If1';If2']); IpABC=Tf2p*conj([If0';If1';If2']);
@ -356,7 +363,7 @@ noLoadi=Loadi(1:6:end);
% noLoadi=[1,19]; % noLoadi=[1,19];
% noLoadi=[4,3,2,5,6,7,11,14,22]; % noLoadi=[4,3,2,5,6,7,11,14,22];
noLoadi=Loadi; noLoadi=Loadi;
noLoadi=[5]; % noLoadi=[5];
wI0r(ismember(Loadi,noLoadi))=1e100000; wI0r(ismember(Loadi,noLoadi))=1e100000;
wI0i(ismember(Loadi,noLoadi))=1e100000; wI0i(ismember(Loadi,noLoadi))=1e100000;
wI1r(ismember(Loadi,noLoadi))=1e100000; wI1r(ismember(Loadi,noLoadi))=1e100000;

15
test.m
View File

@ -14,10 +14,14 @@ totalTime2Sum=10000;
totalTime0Sum=10000; totalTime0Sum=10000;
timeOfPrepare=0.00272333333;%×¼±¸Êý¾ÝµÄʱ¼ä timeOfPrepare=0.00272333333;%×¼±¸Êý¾ÝµÄʱ¼ä
timeOfPrepare=0; timeOfPrepare=0;
eV2Sum=0;
eV0Sum=0;
eI2Sum=0;
eI0Sum=0;
N=200; N=200;
loopN=1; loopN=1;
while 1 while 1
[JMeasurement,AME_Volt,AME_VAngle,AME_PD,AME_QD,AME_mVolt,AME_mPD,AME_mQD,isConverged,totalTime1,totalTime2,totalTime0]=run(); [JMeasurement,AME_Volt,AME_VAngle,AME_PD,AME_QD,AME_mVolt,AME_mPD,AME_mQD,isConverged,totalTime1,totalTime2,totalTime0,eV2,eV0,eI2,eI0]=run();
if isConverged==0 if isConverged==0
continue continue
end end
@ -33,6 +37,11 @@ while 1
AME_mVoltSum=AME_mVoltSum+AME_mVolt; AME_mVoltSum=AME_mVoltSum+AME_mVolt;
AME_mPDSum=AME_mPDSum+AME_mPD; AME_mPDSum=AME_mPDSum+AME_mPD;
AME_mQDSum=AME_mQDSum+AME_mQD; AME_mQDSum=AME_mQDSum+AME_mQD;
eV2Sum=eV2Sum+eV2;
eV0Sum=eV0Sum+eV0;
eI2Sum=eI2Sum+eI2;
eI0Sum=eI0Sum+eI0;
%ʱ¼ä %ʱ¼ä
% totalTime1Sum=totalTime1Sum+totalTime1; % totalTime1Sum=totalTime1Sum+totalTime1;
% totalTime2Sum=totalTime2Sum+totalTime2; % totalTime2Sum=totalTime2Sum+totalTime2;
@ -49,6 +58,10 @@ AME_QDSum=AME_QDSum/N;
AME_mVoltSum=AME_mVoltSum/N; AME_mVoltSum=AME_mVoltSum/N;
AME_mPDSum=AME_mPDSum/N; AME_mPDSum=AME_mPDSum/N;
AME_mQDSum=AME_mQDSum/N; AME_mQDSum=AME_mQDSum/N;
eV2Sum=eV2Sum/N;
eV0Sum=eV0Sum/N;
eI2Sum=eI2Sum/N;
eI0Sum=eI0Sum/N;
totalTime1Sum=totalTime1Sum+timeOfPrepare; totalTime1Sum=totalTime1Sum+timeOfPrepare;
totalTime2Sum=totalTime2Sum+timeOfPrepare; totalTime2Sum=totalTime2Sum+timeOfPrepare;