diff --git a/AssignXX.m b/AssignXX.m index 2ff0f36..870f433 100644 --- a/AssignXX.m +++ b/AssignXX.m @@ -1,10 +1,4 @@ function [deltZ,deltL,deltW,deltU,deltX,deltY]=AssignXX1(XX,ContrlCount,RestraintCount,Busnum) -% deltX=XX(1:14); -% deltY=XX(15:24); -% deltZ=XX(25:38); -% deltW=XX(39:52); -% deltL=XX(53:66); -% deltU=XX(67:80); deltX=XX(1:ContrlCount); k1=ContrlCount+2*Busnum; deltY=XX(ContrlCount+1:k1); diff --git a/FormAlphaD.m b/FormAlphaD.m index bd264d2..f5b2ff1 100644 --- a/FormAlphaD.m +++ b/FormAlphaD.m @@ -3,7 +3,7 @@ tdeltZinx=find(deltZ<0); tdeltWinx=find(deltW>0); t1=-Init_Z(tdeltZinx)./deltZ(tdeltZinx); t2=-Init_W(tdeltWinx)./deltW(tdeltWinx); -t3=[t1,t2]; +t3=[t1;t2]; t4=min(t3); AlphaD=0.9995*min([t4 1]); end \ No newline at end of file diff --git a/SolveIt.m b/SolveIt.m index 039d223..c41122f 100644 --- a/SolveIt.m +++ b/SolveIt.m @@ -12,14 +12,12 @@ yy=[LxComa;-Ly]; t=size(Loadi,1)*2; aa(t+Balance,:)=0; aa(:,t+Balance)=0; -%aa(t+Balance,t+Balance)=1; aa=aa+sparse(t+Balance,t+Balance,ones(length(Balance),1),ContrlCount+2*Busnum,ContrlCount+2*Busnum); deltG(t+Balance,:)=0; %% t=size(Loadi,1)*2+Busnum*1; aa(t+Balance,:)=0; aa(:,t+Balance)=0; -%aa(t+Balance,t+Balance)=1; aa=aa+sparse(t+Balance,t+Balance,ones(length(Balance),1),ContrlCount+2*Busnum,ContrlCount+2*Busnum); deltG(t+Balance,:)=0; %% diff --git a/func_ddh.m b/func_ddh.m index 92f4f86..5905136 100644 --- a/func_ddh.m +++ b/func_ddh.m @@ -1,8 +1,8 @@ function ddh=func_ddh(Volt,Init_Y,Busnum,Y,UAngel,r,c,Angle,Loadi,ContrlCount) mat_AngleIJ=sparse(r,c,UAngel(r)-UAngel(c)-Angle,Busnum,Busnum); mat_INV_AngleIJ=mat_AngleIJ'; -yP=Init_Y(1:length(Init_Y)/2);%暂时改这里 20111227 -yQ=Init_Y(length(Init_Y)/2+1:length(Init_Y));%暂时改这里 20111227 +yP=Init_Y(1:length(Init_Y)/2); +yQ=Init_Y(length(Init_Y)/2+1:end); t1=-sparse(1:Busnum,1:Busnum,Y.*cos(mat_INV_AngleIJ)*sparse(1:Busnum,1:Busnum,Volt,Busnum,Busnum)*yP,Busnum,Busnum); t2=sparse(1:Busnum,1:Busnum,sparse(1:Busnum,1:Busnum,Volt,Busnum,Busnum)*yP,Busnum,Busnum)*Y.*cos(mat_AngleIJ); t3=(t1+t2)*sparse(1:Busnum,1:Busnum,Volt,Busnum,Busnum); @@ -15,7 +15,7 @@ t1=diag( Y.*sin(mat_INV_AngleIJ)*diag(Volt)*yP); t2=diag(yP)*Y.*sin(mat_AngleIJ)*diag(Volt); t3=-diag(yP)*diag(Y.*sin(mat_AngleIJ)*Volt); t4=-Y.*sin(mat_INV_AngleIJ)*diag( diag(Volt)*yP ); -ddPdTdV=t1+t2+t3+t4;%存疑与我的不一样 +ddPdTdV=t1+t2+t3+t4; t1=Y.*cos(mat_INV_AngleIJ)*diag(yP); t2=diag(yP)*Y.*cos(mat_AngleIJ); ddPdVdV=t1+t2; diff --git a/run.m b/run.m index 4c671ea..1305c13 100644 --- a/run.m +++ b/run.m @@ -181,6 +181,7 @@ SEphaseCSpotLoadP=zeros(length(phaseCSpotLoadP),1); SEphaseASpotLoadQ=zeros(length(phaseASpotLoadQ),1); SEphaseBSpotLoadQ=zeros(length(phaseBSpotLoadQ),1); SEphaseCSpotLoadQ=zeros(length(phaseCSpotLoadQ),1); +% SEVmf1=sparse(ones(busNum,1)); SEVaf1=sparse(zeros(busNum,1)); SEPD=sparse(length(Loadi),1); @@ -188,7 +189,7 @@ SEQD=sparse(length(Loadi),1); KK=0; plotGap=zeros(1,60); %初始化 -%状态量为 SEPD SEQD SEVmf1 SEVaf1 +%状态量为 SEPD SEQD SEVmf1 SEVaf1 RestraintCount=length(SEVmf1)+length(Loadi)*2; ContrlCount=length(SEVmf1)*2+length(Loadi)*2; CenterA=0.1; @@ -198,27 +199,35 @@ 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); -Init_u=Gap/2/RestraintCount*CenterA; -deltH=func_deltH(busNum,SEVmf1,fsY1amp,SEVaf1,r,c,fsY1ang,Loadi); -deltG=func_deltG(busNum,Loadi); -L_1Z=diag(Init_Z./Init_L); -U_1W=diag(Init_W./Init_U); -deltdeltF=func_deltdeltF(SEPD,ContrlCount); -ddh=func_ddh(SEVmf1,Init_Y,busNum,fsY1amp,SEVaf1,r,c,fsY1ang,Loadi,ContrlCount); -ddg=func_ddg(); -deltF=func_deltF(SEPD,ContrlCount); -Luu=Init_U.*Init_W+Init_u*ones(RestraintCount,1); -Lul=Init_L.*Init_Z-Init_u*ones(RestraintCount,1); -Mat_G=FormG(SEVaf1,SEPD,SEQD); -Mat_H=FormH(busNum,SEVmf1,PG,SEPD,QG,SEQD,fsY1amp,SEVaf1,r,c,fsY1ang); -Ly=Mat_H; -Lz=FormLz(Mat_G,Init_L,busNum,Loadi); -Lw=FormLw(Mat_G,Init_U,busNum,Loadi); -Lx=FormLx(deltF,deltH,Init_Y,deltG,Init_Z,Init_W); -YY=FormYY(Lul,Lz,Ly,Luu,Lw,Lx); -%% 开始解方程 -fprintf('迭代次数 %d Gap %f\n',KK+1,plotGap(KK+1)); -XX=SolveIt(deltF,deltG,Init_L,Init_Z,Init_U,Init_W,deltdeltF,ddh,ddg,deltH,Init_Y,Ly,Lz,ContrlCount,Lw,Lul,Luu,Lx,Balance,busNum,Loadi); -[deltZ,deltL,deltW,deltU,deltX,deltY]=AssignXX(XX,ContrlCount,RestraintCount,busNum); -[Init_Z,Init_L,Init_W,Init_U,Init_Y,PG,QG,Volt,UAngel,PD,QD]=Modification(Init_Z,Init_L,Init_W,Init_U,Init_Y,deltZ,deltL,deltW,deltU,deltX,deltY,PG,QG,SEVmf1,SEVaf1,ContrlCount,Balance,busNum,PD,QD,Loadi); -Gap=(Init_L'*Init_Z-Init_U'*Init_W); \ No newline at end of file +PG=sparse(busNum,1); +PG(Balance)=0.1105; +QD=sparse(busNum,1); +QG(Balance)=0.0984; +while Gap>1e-5 && KK<60 + KK=KK+1; + Init_u=Gap/2/RestraintCount*CenterA; + deltH=func_deltH(busNum,SEVmf1,fsY1amp,SEVaf1,r,c,fsY1ang,Loadi); + deltG=func_deltG(busNum,Loadi); + L_1Z=diag(Init_Z./Init_L); + U_1W=diag(Init_W./Init_U); + deltdeltF=func_deltdeltF(SEPD,ContrlCount); + ddh=func_ddh(SEVmf1,Init_Y,busNum,fsY1amp,SEVaf1,r,c,fsY1ang,Loadi,ContrlCount); + ddg=func_ddg(); + deltF=func_deltF(SEPD,ContrlCount); + Luu=Init_U.*Init_W+Init_u*ones(RestraintCount,1); + Lul=Init_L.*Init_Z-Init_u*ones(RestraintCount,1); + Mat_G=FormG(SEVaf1,SEPD,SEQD); + Mat_H=FormH(busNum,SEVmf1,PG,SEPD,QG,SEQD,fsY1amp,SEVaf1,r,c,fsY1ang); + Ly=Mat_H; + Lz=FormLz(Mat_G,Init_L,busNum,Loadi); + Lw=FormLw(Mat_G,Init_U,busNum,Loadi); + Lx=FormLx(deltF,deltH,Init_Y,deltG,Init_Z,Init_W); + YY=FormYY(Lul,Lz,Ly,Luu,Lw,Lx); + %% 开始解方程 + plotGap(KK)=Gap; + fprintf('迭代次数 %d Gap %f\n',KK,plotGap(KK)); + XX=SolveIt(deltF,deltG,Init_L,Init_Z,Init_U,Init_W,deltdeltF,ddh,ddg,deltH,Init_Y,Ly,Lz,ContrlCount,Lw,Lul,Luu,Lx,Balance,busNum,Loadi); + [deltZ,deltL,deltW,deltU,deltX,deltY]=AssignXX(XX,ContrlCount,RestraintCount,busNum); + [Init_Z,Init_L,Init_W,Init_U,Init_Y,PG,QG,Volt,UAngel,PD,QD]=Modification(Init_Z,Init_L,Init_W,Init_U,Init_Y,deltZ,deltL,deltW,deltU,deltX,deltY,PG,QG,SEVmf1,SEVaf1,ContrlCount,Balance,busNum,PD,QD,Loadi); + Gap=(Init_L'*Init_Z-Init_U'*Init_W); +end \ No newline at end of file