From 1599e58150de5af3e40272360e4ff6155878e196 Mon Sep 17 00:00:00 2001 From: "dugg@lab-desk" Date: Sun, 25 May 2014 12:13:10 +0800 Subject: [PATCH] =?UTF-8?q?=E5=85=88=E7=AE=97=E8=BF=9E=E7=BB=AD=E9=87=8F?= =?UTF-8?q?=EF=BC=8C=E5=86=8D=E7=AE=97=E7=A6=BB=E6=95=A3=E9=87=8F=EF=BC=8C?= =?UTF-8?q?=E6=94=B6=E6=95=9B=E4=BA=86=EF=BC=8C=E8=80=8C=E4=B8=94=E7=AE=97?= =?UTF-8?q?=E7=9A=84=E7=BB=93=E6=9E=9C=E6=98=AF=E5=AF=B9=E7=9A=84=E3=80=82?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Signed-off-by: dugg@lab-desk --- OPF.m | 55 ++++++++++++++++++++++++++++++++++++++++++++++++++++--- 1 file changed, 52 insertions(+), 3 deletions(-) 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);