diff --git a/OPF.m b/OPF.m index aa29bcc..cfdf2a9 100644 --- a/OPF.m +++ b/OPF.m @@ -1,7 +1,7 @@ tic clear %[kmax,Precision,UAngel,Volt,Busnum,PVi,PVu,Balance,Y,Angle,P0,Q0,r,c,GB,Linei,Linej,Transfori,Transforj,GenU,GenL,GenC,LineLimti,LineLimtj,LinePLimt,PG,QG,PD,QD,CenterA,LineCount,PGi,PVQU,PVQL]=pf('5sj.txt'); -[kmax,Precision,UAngel,Volt,Busnum,PVi,PVu,Balance,Y,Angle,P0,Q0,r,c,GB,Linei,Linej,Transfori,Transforj,GenU,GenL,GenC,PG,QG,PD,QD,CenterA,PGi,PVQU,PVQL]=pf('ieee118.dat'); +[kmax,Precision,UAngel,Volt,Busnum,PVi,PVu,Balance,Y,Angle,P0,Q0,r,c,GB,Linei,Linej,Transfori,Transforj,GenU,GenL,GenC,PG,QG,PD,QD,CenterA,PGi,PVQU,PVQL]=pf('ieee10471PG.dat'); GB=full(GB); %PVi电压节点序号 %PVu电压节点电压标幺值 @@ -10,8 +10,9 @@ UAngel*180/3.1415926; %sprintf('%f\n',Volt); %sprintf('%f\n',Angel); %% 初值 +PG0=PG; PD0=PD; -[Volt,UAngel,Init_Z,Init_W,Init_L,Init_U,Init_Y,PG,QG,RestraintCount,wD,PD]=OPF_Init(Busnum,Balance,PG,QG,Volt,GenU,GenL,PVi,PGi,PVQU,PVQL,PD0); +[Volt,UAngel,Init_Z,Init_W,Init_L,Init_U,Init_Y,PG,QG,RestraintCount,wG,wD,PD]=OPF_Init(Busnum,Balance,PG,QG,Volt,GenU,GenL,PVi,PGi,PVQU,PVQL,PD0); Gap=(Init_L*Init_Z'-Init_U*Init_W'); KK=0; plotGap=zeros(1,50); @@ -33,13 +34,13 @@ while(abs(Gap)>Precision) L_1Z=diag(Init_Z./Init_L); U_1W=diag(Init_W./Init_U); %% 形成海森阵 - deltdeltF=func_deltdeltF(Busnum,GenC,PVi,PGi,wD,PD0); + deltdeltF=func_deltdeltF(Busnum,GenC,PVi,PGi,wG,wD,PD0); %% 形成ddHy ddh=func_ddh3(AngleIJMat,GB,Volt,Init_Y,Busnum,PVi,PGi,Y); %% 开始构建ddg ddg=func_ddg(PGi,PVi,Busnum,RestraintCount); %% 开始构建deltF - deltF=func_deltF(PG,PVi,GenC,ContrlCount,PGi,wD,PD0,PD,Busnum); + deltF=func_deltF(PG,PVi,GenC,ContrlCount,PGi,wG,wD,PG0,PD0,PD,Busnum); %% 形成方程矩阵 Luu=Init_U'.*Init_W'+Init_u*ones(RestraintCount,1); @@ -61,7 +62,7 @@ while(abs(Gap)>Precision) end fprintf('迭代次数%d\n',KK); %CalCost(GenC,PG,PGi); -ObjectiveFun(PG,PGi,PD,PD0,wD) +ObjectiveFun(PG,PG0,PGi,PD,PD0,wG,wD) DrawGap(plotGap); Volt=Volt'; toc diff --git a/OPF_Init.m b/OPF_Init.m index 2451b89..97142d6 100644 --- a/OPF_Init.m +++ b/OPF_Init.m @@ -1,4 +1,4 @@ -function [Volt,UAngel,Init_Z,Init_W,Init_L,Init_U,Init_Y,PG,QG,RestraintCount,wD,PD]=OPF_Init(Busnum,Balance,PG,QG,Volt,GenU,GenL,PVi,PGi,PVQU,PVQL,PD0) +function [Volt,UAngel,Init_Z,Init_W,Init_L,Init_U,Init_Y,PG,QG,RestraintCount,wG,wD,PD]=OPF_Init(Busnum,Balance,PG,QG,Volt,GenU,GenL,PVi,PGi,PVQU,PVQL,PD0) RestraintCount=size(PVi,1)+size(PGi,1)+Busnum*2; %约束条件数 t_Bal_volt=Volt(Balance); Volt=ones(1,Busnum); @@ -22,8 +22,9 @@ tQL=PVQL(:,1);% PG(PGi)=(tPU+tPL)/2; %QG(4:5)=[0 1.45]; QG(PVi)=(tQU+tQL)/2; +wG=ones(size(PGi,1),1); wD=ones(Busnum,1); -wD(Balance)=0; +%wD(Balance)=0; PD=1*PD0; %PD(PD==0)=.2; end \ No newline at end of file diff --git a/ObjectiveFun.m b/ObjectiveFun.m index 2eb1fcd..d1a21d6 100644 --- a/ObjectiveFun.m +++ b/ObjectiveFun.m @@ -1,5 +1,6 @@ -function [out_arg]=ObjectiveFun(PG,PGi,PD,PD0,wD) - +function [out_arg]=ObjectiveFun(PG,PG0,PGi,PD,PD0,wG,wD) +t1=PG(PGi)-PG0(PGi); +t2=wG.*t1.*t1; t3=wD.*((PD-PD0).*(PD-PD0)); -out_arg= sum(t3); +out_arg= sum(t2)+sum(t3); end \ No newline at end of file diff --git a/func_deltF.asv b/func_deltF.asv index 74784d3..8fb2e42 100644 --- a/func_deltF.asv +++ b/func_deltF.asv @@ -1,4 +1,4 @@ -function deltF=func_deltF(PG,PVi,GenC,ContrlCount,PGi,wD,PD0,PD) +function deltF=func_deltF(PG,PVi,GenC,ContrlCount,PGi,wG,wD,PG0,PD0,PD,Busnum) %t1=PG(setdiff(PVi,Balance)); % t2=Volt'*Volt; % t3=real(GB).*cos(AngleIJMat)+imag(GB).*sin(AngleIJMat); @@ -7,12 +7,13 @@ function deltF=func_deltF(PG,PVi,GenC,ContrlCount,PGi,wD,PD0,PD) % PBal=t5(Balance); % PPG=([PQ(1),PBal])';%暂时用土办法处理一下 %% -t1=2*wD.*(PD-PD0); +t1=2 +t2=2*wD.*(PD-PD0); deltF=[ zeros(size(PGi)); zeros(size(PVi)); - t1; - zeros(ContrlCount-size(PGi,1)-,1); + t2; + zeros(2*Busnum,1); ]; end \ No newline at end of file diff --git a/func_deltF.m b/func_deltF.m index b00a169..ca22864 100644 --- a/func_deltF.m +++ b/func_deltF.m @@ -1,4 +1,4 @@ -function deltF=func_deltF(PG,PVi,GenC,ContrlCount,PGi,wD,PD0,PD,Busnum) +function deltF=func_deltF(PG,PVi,GenC,ContrlCount,PGi,wG,wD,PG0,PD0,PD,Busnum) %t1=PG(setdiff(PVi,Balance)); % t2=Volt'*Volt; % t3=real(GB).*cos(AngleIJMat)+imag(GB).*sin(AngleIJMat); @@ -7,11 +7,12 @@ function deltF=func_deltF(PG,PVi,GenC,ContrlCount,PGi,wD,PD0,PD,Busnum) % PBal=t5(Balance); % PPG=([PQ(1),PBal])';%暂时用土办法处理一下 %% -t1=2*wD.*(PD-PD0); +t1=2*wG.*(PG(PGi) - PG0(PGi) ); +t2=2*wD.*(PD-PD0); deltF=[ - zeros(size(PGi)); - zeros(size(PVi)); t1; + zeros(size(PVi)); + t2; zeros(2*Busnum,1); ]; diff --git a/func_deltdeltF.m b/func_deltdeltF.m index 55773d3..4661b57 100644 --- a/func_deltdeltF.m +++ b/func_deltdeltF.m @@ -1,7 +1,7 @@ -function deltdeltF=func_deltdeltF(Busnum,GenC,PVi,PGi,wD,PD0) +function deltdeltF=func_deltdeltF(Busnum,GenC,PVi,PGi,wG,wD,PD0) ContrlCount=size(PVi,1)+size(PGi,1)+Busnum*3; %P,Q,Volt theta这些控制变量数 -C=[zeros(size(PGi))' zeros(size(PVi))' wD']; +C=[wG' zeros(size(PVi))' wD']; deltdeltF=[diag(C)*2,zeros(size(C,2),ContrlCount-size(C,2)); zeros(ContrlCount-size(C,2),ContrlCount); ];