先算连续量,再算离散量,收敛了,而且算的结果是对的。

Signed-off-by: dugg@lab-desk <dugg@lab-desk>
This commit is contained in:
dugg@lab-desk 2014-05-25 12:13:10 +08:00
parent c1727c877f
commit 1599e58150
1 changed files with 52 additions and 3 deletions

55
OPF.m
View File

@ -61,6 +61,8 @@ Vbi=sparse(0.5*ones(Busnum,1));
PDbi=sparse(0.5*ones(length(Loadi),1)); PDbi=sparse(0.5*ones(length(Loadi),1));
QDbi=sparse(0.5*ones(length(Loadi),1)); QDbi=sparse(0.5*ones(length(Loadi),1));
eps=10; eps=10;
%
fprintf('1\n');
while(abs(Gap)>Precision*1) while(abs(Gap)>Precision*1)
if KK>kmax if KK>kmax
break; break;
@ -90,9 +92,56 @@ while(abs(Gap)>Precision*1)
Mat_G=FormG(Volt,PD,QD,Loadi,bigM,mVolt,rVolt,sigma,Vbi,PDbi,QDbi,mPD,mQD); Mat_G=FormG(Volt,PD,QD,Loadi,bigM,mVolt,rVolt,sigma,Vbi,PDbi,QDbi,mPD,mQD);
Mat_H=FormH(Busnum,Volt,PG,PD,QG,QD,Y,UAngel,r,c,Angle,Loadi); Mat_H=FormH(Busnum,Volt,PG,PD,QG,QD,Y,UAngel,r,c,Angle,Loadi);
Ly=Mat_H; Ly=Mat_H;
% if KK>5 Lz=FormLz(Mat_G,Init_L,GenL,Busnum,PVQL,PD,RealPD,RealQD,Loadi,KK,PF,eps);
% eps=eps*1.9; Lw=FormLw(Mat_G,Init_U,GenU,Busnum,PVQU,PD,PD0,QD0,Loadi,KK,PF,eps);
% end 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,ddgzw,deltH,Init_Y,Ly,Lz,ContrlCount,Lw,Lul,Luu,RestraintCount,Lx,Balance,PVi,PGi,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,Vbi,PDbi,QDbi]=Modification(Init_Z,Init_L,Init_W,Init_U,Init_Y,deltZ,deltL,deltW,deltU,deltX,deltY,PG,QG,Volt,UAngel,PVi,ContrlCount,Balance,Busnum,PGi,PD,QD,Loadi,Vbi,PDbi,QDbi);
Gap=(Init_L*Init_Z'-Init_U*Init_W');
KK=KK+1;
end
%
Gap=1000;
KK=0;
fprintf('\n');
fprintf('2\n');
while(abs(Gap)>Precision*1)
if KK>kmax
break;
end
plotGap(KK+1)=Gap;
Init_u=Gap/2/RestraintCount*CenterA;
AngleIJMat=0;
%% OPF
%%
deltH=func_deltH(Busnum,Volt,PVi,Y,PGi,UAngel,r,c,Angle,Loadi);
%%
deltG=func_deltG(Busnum,PVi,PGi,Loadi,PD,QD,Vbi,PDbi,QDbi);
%%
L_1Z=diag(Init_Z./Init_L);
U_1W=diag(Init_W./Init_U);
%%
deltdeltF=func_deltdeltF(PVi,wPG,wQG,wPD,wQD,ContrlCount);
%% ddHy
ddh=func_ddh(Volt,Init_Y,Busnum,PVi,PGi,Y,UAngel,r,c,Angle,Loadi,ContrlCount);
%% ddg
ddgzw=func_ddg(Busnum,ContrlCount,Loadi,Init_Z,Init_W);
%% deltF
deltF=func_deltF(PG,QG,PVi,PGi,wPG,wQG,wPD,wQD,PG0,QG0,PD0,PD,QD,QD0,Busnum,Loadi);
%%
Luu=Init_U'.*Init_W'+Init_u*ones(RestraintCount,1);
Lul=Init_L'.*Init_Z'-Init_u*ones(RestraintCount,1);
Mat_G=FormG(Volt,PD,QD,Loadi,bigM,mVolt,rVolt,sigma,Vbi,PDbi,QDbi,mPD,mQD);
Mat_H=FormH(Busnum,Volt,PG,PD,QG,QD,Y,UAngel,r,c,Angle,Loadi);
Ly=Mat_H;
if KK>5
eps=eps*0.1;
end
Lz=FormLz(Mat_G,Init_L,GenL,Busnum,PVQL,PD,RealPD,RealQD,Loadi,KK,PF,eps); Lz=FormLz(Mat_G,Init_L,GenL,Busnum,PVQL,PD,RealPD,RealQD,Loadi,KK,PF,eps);
Lw=FormLw(Mat_G,Init_U,GenU,Busnum,PVQU,PD,PD0,QD0,Loadi,KK,PF,eps); Lw=FormLw(Mat_G,Init_U,GenU,Busnum,PVQU,PD,PD0,QD0,Loadi,KK,PF,eps);
Lx=FormLx(deltF,deltH,Init_Y,deltG,Init_Z,Init_W); Lx=FormLx(deltF,deltH,Init_Y,deltG,Init_Z,Init_W);