diff --git a/FormGstate1.m b/FormGstate1.m index 7972cd7..3c9a919 100644 --- a/FormGstate1.m +++ b/FormGstate1.m @@ -1,4 +1,4 @@ -function Mat_G=FormGstate1(I1r,I1i) +function Mat_G=FormGstate1(I1r,I1i,Loadi) %t1=PG(PVi); %GP=t1;%发电机P %GP=[4.5 4.5]'; @@ -10,7 +10,7 @@ function Mat_G=FormGstate1(I1r,I1i) % t4=sum(t3,2);%发电机Q %GQ=t4; Mat_G=[ - I1r; - I1i; + I1r(Loadi); + I1i(Loadi); ]; end \ No newline at end of file diff --git a/IPMLoop.m b/IPMLoop.m index 16cc370..d62d9b5 100644 --- a/IPMLoop.m +++ b/IPMLoop.m @@ -8,12 +8,15 @@ I1r=real(guessIf1); I1i=imag(guessIf1); KK=0; plotGap=zeros(1,60); -%初始化 -%状态量为 SEPD SEQD SEVmf1 SEVaf1 -RestraintCount=length(Loadi)*2; + % ContrlCount=busNum*2+length(Loadi)*2; noBoundedLoadi=setdiff(Loadi,noLoadi); % noBoundedLoadi=[]; +boundedLoadi=setdiff(Loadi,noBoundedLoadi); +%初始化 +%状态量为 SEPD SEQD SEVmf1 SEVaf1 +% RestraintCount=length(Loadi)*2; +RestraintCount=length(boundedLoadi)*2; ContrlCount=busNum*2+length(Loadi)*2; CenterA=0.1; Init_Z=1*sparse(ones(RestraintCount,1)); @@ -22,7 +25,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=100; +kmax=20; isSetBound=0; largerBound=0; realBound=1; @@ -37,7 +40,8 @@ while(abs(Gap)>1e-4) %% 形成等式约束的雅克比 deltH=func_deltH(busNum,fsY1,Loadi,Balance); %% 形成不等式约束的雅克比 - deltG=func_deltGstate1(busNum,Loadi,I1r,I1i); +% deltG=func_deltGstate1(busNum,Loadi,I1r,I1i); +deltG=func_deltGstate1(busNum,Loadi,boundedLoadi,I1r,I1i); %% % L_1Z=diag(Init_Z./Init_L); % U_1W=diag(Init_W./Init_U); @@ -55,22 +59,24 @@ while(abs(Gap)>1e-4) %% Luu=Init_U.*Init_W+Init_u*ones(RestraintCount,1); Lul=Init_L.*Init_Z-Init_u*ones(RestraintCount,1); - Mat_G=FormGstate1(I1r,I1i); +% Mat_G=FormGstate1(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); + Ir=real(guessIf1(ismember(Loadi,boundedLoadi)) ); pIr=find(Ir>0); nIr=find(Ir<0); - Ii=imag(guessIf1); + Ii=imag(guessIf1(ismember(Loadi,boundedLoadi))); pIi=find(Ii>0); nIi=find(Ii<0); - lower=ones(length(Loadi)*2,1); +% lower=ones(length(Loadi)*2,1); +lower=ones(length(boundedLoadi)*2,1); % if abs(Gap)<1 || isSetBound==1 -% if abs(Gap)<0.01 +% if abs(Gap)<0.001 % largerBound=0; % realBound=1; % end @@ -85,8 +91,8 @@ while(abs(Gap)>1e-4) 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))=-0.9; - lower(ismember(Loadi,noBoundedLoadi)+length(Ir))=-0.9; +% lower( ismember(Loadi,noBoundedLoadi))=-1.9; +% lower(ismember(Loadi,noBoundedLoadi)+length(Ir))=-1.9; end if largerBound==1 @@ -112,13 +118,13 @@ while(abs(Gap)>1e-4) % I1r(I1r./lowerR>0.9998)=lowerR(I1r./lowerR>0.9998); % I1i(I1i./lowerI>0.9998)=lowerI(I1i./lowerI>0.9998); - Ir=real(guessIf1); + Ir=real(guessIf1(ismember(Loadi,boundedLoadi))); pIr=find(Ir>0); nIr=find(Ir<0); - Ii=imag(guessIf1); + Ii=imag(guessIf1(ismember(Loadi,boundedLoadi))); pIi=find(Ii>0); nIi=find(Ii<0); - upper=ones(length(Loadi)*2,1); + upper=ones(length(boundedLoadi)*2,1); % if abs(Gap)<1 || isSetBound==1 if realBound==1 % upper(pIr)=(1.10+3/(KK+1))*Ir(pIr); @@ -130,8 +136,8 @@ while(abs(Gap)>1e-4) 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))=0.9; - upper(ismember(Loadi,noBoundedLoadi)+length(Ir))=0.9; +% upper(ismember(Loadi,noBoundedLoadi))=1.9; +% upper(ismember(Loadi,noBoundedLoadi)+length(Ir))=1.9; end if largerBound==1 @@ -160,13 +166,14 @@ while(abs(Gap)>1e-4) [deltZ,deltL,deltW,deltU,deltX,deltY]=AssignXX(XX,ContrlCount,RestraintCount,busNum,Balance); [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=max([max(abs(Mat_H)),max(Lz),max(Lw) ]); fprintf('Gap %f :%d\n',full(Gap),KK+1); KK=KK+1; 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)<0.00001 && KK