diff --git a/OPF.m b/OPF.m index 7ac72af..a7f6b14 100644 --- a/OPF.m +++ b/OPF.m @@ -61,6 +61,8 @@ Vbi=sparse(0.5*ones(Busnum,1)); PDbi=sparse(0.5*ones(length(Loadi),1)); QDbi=sparse(0.5*ones(length(Loadi),1)); eps=10; +% 第一遍,算连续的值 +fprintf('第1次迭代,算连续量。\n'); while(abs(Gap)>Precision*1) if KK>kmax 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_H=FormH(Busnum,Volt,PG,PD,QG,QD,Y,UAngel,r,c,Angle,Loadi); Ly=Mat_H; -% if KK>5 -% eps=eps*1.9; -% end + 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); + 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); 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);