From beaa6855d4282adebb77c98050e8cf74a26792f3 Mon Sep 17 00:00:00 2001 From: "dugg@lab-desk" Date: Mon, 2 Feb 2015 21:28:56 +0800 Subject: [PATCH] =?UTF-8?q?=E5=8F=88=E8=B0=83=E4=BA=86=E4=B8=80=E4=B8=8B?= =?UTF-8?q?=EF=BC=8C=E6=9B=B4=E7=90=86=E6=83=B3=E4=BA=86=E3=80=82?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Signed-off-by: dugg@lab-desk --- subOPF.m | 111 ++++++++++++++++++++++++++++--------------------------- 1 file changed, 56 insertions(+), 55 deletions(-) diff --git a/subOPF.m b/subOPF.m index 2ed2d85..2ded212 100644 --- a/subOPF.m +++ b/subOPF.m @@ -90,7 +90,7 @@ lVolt(noMeasurei)=0.7*mVolt(noMeasurei);% uVolt(noMeasurei)=0.7*mVolt(noMeasurei); %错误数据 %mVolt(2)=5; -bigM=0.000003; +% bigM=0.000003; Vbi=sparse(0.5*ones(Busnum,1)); % Vbi(2)=1; PDbi=sparse(0.5*ones(length(Loadi),1)); @@ -124,7 +124,7 @@ while(abs(Gap)>Precision*10) %% 形成方程矩阵 Luu=Init_U'.*Init_W'+Init_u*ones(RestraintCount,1); Lul=Init_L'.*Init_Z'-Init_u*ones(RestraintCount,1); - bigM=0.7; + bigM=1; Mat_G=FormG(Volt,PD,QD,Loadi,bigM,mVolt,rVolt,sigma,Vbi,PDbi,QDbi,mPD,mQD,uPD,lPD,uQD,lQD,uVolt,lVolt); Mat_H=FormH(Busnum,Volt,PG,PD,QG,QD,Y,UAngel,r,c,Angle,Loadi); Ly=Mat_H; @@ -146,69 +146,70 @@ end % Gap=(Init_L*Init_Z'-Init_U*Init_W'); Gap=1; % KK=0; -eps=0.00001; +% eps=0.00001; fprintf('\n'); fprintf('第2次迭代，算离散量。\n'); -while eps>0.000001 +% while eps>0.000001 % Init_Z=sparse(ones(1,RestraintCount)); % Init_W=sparse(-1*ones(1,RestraintCount)); % Init_L=1*sparse(ones(1,RestraintCount)); % Init_U=1*sparse(ones(1,RestraintCount)); % Init_Y=sparse(1,2*Busnum);%与学姐一致 - while(abs(Gap)>Precision*100) - 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); - bigM=0.7; - Mat_G=FormG(Volt,PD,QD,Loadi,bigM,mVolt,rVolt,sigma,Vbi,PDbi,QDbi,mPD,mQD,uPD,lPD,uQD,lQD,uVolt,lVolt); - Mat_H=FormH(Busnum,Volt,PG,PD,QG,QD,Y,UAngel,r,c,Angle,Loadi); - Ly=Mat_H; - 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); - %%取各分量 - % Vbi_=Vbi; - % PDbi_=PDbi; - % QDbi_=QDbi; -% Vbi=1./(1+exp( (0.06-Vbi) /0.01))*0.1; - [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'); - -% Vbi(Vbi>0.01)=1; - KK=KK+1; +while(abs(Gap)>Precision*10) + if KK>kmax + % break; end - eps=eps*0.4; -% eps=Gap - Gap=100; + 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); + bigM=0.7; + eps=Gap*0.1; + Mat_G=FormG(Volt,PD,QD,Loadi,bigM,mVolt,rVolt,sigma,Vbi,PDbi,QDbi,mPD,mQD,uPD,lPD,uQD,lQD,uVolt,lVolt); + Mat_H=FormH(Busnum,Volt,PG,PD,QG,QD,Y,UAngel,r,c,Angle,Loadi); + Ly=Mat_H; + 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); + %%取各分量 + % Vbi_=Vbi; + % PDbi_=PDbi; + % QDbi_=QDbi; + % Vbi=1./(1+exp( (0.06-Vbi) /0.01))*0.1; + [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'); + + % Vbi(Vbi>0.01)=1; + KK=KK+1; end + +% eps=Gap +% Gap=100; +% end %% 计算最大不平衡量 AngleIJ=sparse(r,c,UAngel(r)-UAngel(c)-Angle',Busnum,Busnum); dP=PD+diag(Volt)*Y.*cos(AngleIJ)*Volt';