这是一个不把功率因数作为不等式约束，还是在QD的不等式约束中期作用的版本。

Signed-off-by: unknown <dmy@dmy-PC.(none)>

FormLw.m

@@ -1,8 +1,10 @@
 function Lw=FormLw(Mat_G,Init_U,GenU,Busnum,PVQU,PD0,QD0,Loadi,KK)
 KK=999;
-PU=GenU(:,2);%发电机有功上界
+%PU=GenU(:,2);%发电机有功上界
-QU=PVQU(:,1);%发电机无功上界
+PU=5*ones(length(GenU(:,2)),1);
-VoltU=(1.2)*ones(1,Busnum);
+%QU=PVQU(:,1);%发电机无功上界
+QU=5*ones(length(PVQU(:,1)),1);
+VoltU=(1.1)*ones(1,Busnum);
 %VoltU=10*ones(1,Busnum);
 PDU=PD0(Loadi);
 PDU(PDU>0)=1.200*PDU(PDU>0);

FormLz.m

@@ -1,8 +1,10 @@
 function Lz=FormLz(Mat_G,Init_L,GenL,Busnum,PVQL,PD0,QD0,Loadi,KK)
 KK=999;
-PL=GenL(:,2);%发电机有功下界
+%PL=GenL(:,2);%发电机有功下界
-QL=PVQL(:,1);%发电机无功下界
+PL=-5*ones(length(GenL(:,2)),1);
-VoltL=(0.8)*ones(1,Busnum);
+%QL=PVQL(:,1);%发电机无功下界
+QL=-5*ones(length(PVQL(:,1)),1);
+VoltL=(0.9)*ones(1,Busnum);
 %VoltL=-10*ones(1,Busnum);
 PDL=PD0(Loadi);
 PDL(PDL>0)=0.800*PDL(PDL>0);

OPF.m

@@ -2,7 +2,7 @@ tic
 clc
 clear
 [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,Liner,Linex,Lineb,Transforr,Transforx,Transfork0]= ...
-pf('ieee10471.dat');
+pf('ieee118.dat');
 %pf('D:\Project\青秀降损项目\最小化潮流\最小潮流算例\原始\津头站津视922（3-1）_0.5_120%.txt'); 
 %pf('D:\Project\最小化潮流\最小潮流算例\仙海919.txt'); 
 %pf('c:/file31.txt'); 
@@ -14,8 +14,8 @@ Volt;
 UAngel*180/3.1415926;
 %% 通过潮流计算PG
 AngleIJ=sparse(r,c,UAngel(r)-UAngel(c)-Angle',Busnum,Busnum);
-PGBal=diag(Volt)*Y.*cos(AngleIJ)*Volt';
+PGBal=PD+diag(Volt)*Y.*cos(AngleIJ)*Volt';
-QGBal=diag(Volt)*Y.*sin(AngleIJ)*Volt';
+QGBal=QD+diag(Volt)*Y.*sin(AngleIJ)*Volt';
 %% 初值-即测量值
 PG0=PG;
 QG0=QG;
@@ -25,9 +25,12 @@ PDReal=PD;%
 %PD0(12)=PD0(12)+0.001;
 %%
 PG0(Balance)=PGBal(Balance);
+PG(Balance)=PGBal(Balance);
 QG0(Balance)=QGBal(Balance);
+QG0(PVi)=QGBal(PVi);
+QG(PVi)=QGBal(PVi);
 %%
-[Volt,UAngel,Init_Z,Init_W,Init_L,Init_U,Init_Y,PG,QG,RestraintCount,wPG,wQG,wD,PD,PD0,QD,randPDind,Loadi]=OPF_Init(Busnum,Balance,PG,QG,Volt,GenU,GenL,PVi,PGi,PVQU,PVQL,PD0,QD0,QD,PD);
+[Volt,UAngel,Init_Z,Init_W,Init_L,Init_U,Init_Y,PG,QG,RestraintCount,wPG,wQG,wPD,wQD,PD,PD0,QD,randPDind,Loadi]=OPF_Init(Busnum,Balance,PG,QG,Volt,GenU,GenL,PVi,PGi,PVQU,PVQL,PD0,QD0,QD,PD);
 Gap=(Init_L*Init_Z'-Init_U*Init_W');
 KK=0;
 plotGap=zeros(1,60);
@@ -53,13 +56,13 @@ while(abs(Gap)>Precision)
     L_1Z=diag(Init_Z./Init_L);
     U_1W=diag(Init_W./Init_U);
     %% 形成海森阵
-    deltdeltF=func_deltdeltF(PVi,wPG,wQG,wD,ContrlCount);
+    deltdeltF=func_deltdeltF(PVi,wPG,wQG,wPD,ContrlCount);
     %% 形成ddHy
     ddh=func_ddh(Volt,Init_Y,Busnum,PVi,PGi,Y,UAngel,r,c,Angle,Loadi,ContrlCount);
     %%  开始构建ddg
     ddg=func_ddg(PGi,PVi,Busnum,RestraintCount,Loadi,PD,QD);
     %%  开始构建deltF
-    deltF=func_deltF(PG,QG,PVi,PGi,wPG,wQG,wD,PG0,QG0,PD0,PD,Busnum,Loadi);
+    deltF=func_deltF(PG,QG,PVi,PGi,wPG,wQG,wPD,PG0,QG0,PD0,PD,Busnum,Loadi);
     
     %% 形成方程矩阵
     Luu=Init_U'.*Init_W'+Init_u*ones(RestraintCount,1);
@@ -81,7 +84,7 @@ while(abs(Gap)>Precision)
     KK=KK+1;
 end
 fprintf('迭代次数%d\n',KK);
-%fprintf('Ä¿±êÖµ%f\n',full(ObjectiveFun(PG,PG0,PGi,PD,PD0,wG,wD,Loadi)));
+fprintf('Ä¿±êÖµ%f\n',full(ObjectiveFun(PG,PG0,PGi,QG,QG0,PVi,PD,PD0,QD,QD0,wPG,wQG,wPD,wQD,Loadi)));
 DrawGap(plotGap);
 %% 
 %Volt=full(Volt');
@@ -89,13 +92,13 @@ DrawGap(plotGap);
 %% 统计PD误差
 % absPDLoad=abs(  (PD(Loadi)-PD0(Loadi))./PD0(Loadi)  );
 absPDLoad=abs(  (PD(Loadi)-PDReal(Loadi))./PDReal(Loadi)  );
-maxPDError=max(absPDLoad);
+maxPDError=max(absPDLoad)
 disp('index');
 Loadi(absPDLoad==maxPDError);
 %% 计算总线损
 totalLoss=(sum(PG)-sum(PD(Loadi)))*100;
 fprintf('总的损耗为%f(MW 有名值)\n',full(totalLoss));
 %% 计算各线损
-Lineloss(Linei,Linej,Liner,Linex,Lineb,Transfori,Transforj,Transforr,Transforx,Transfork0,Volt,UAngel);
+%Lineloss(Linei,Linej,Liner,Linex,Lineb,Transfori,Transforj,Transforr,Transforx,Transfork0,Volt,UAngel);
 toc  
 

OPF_Init.m

@@ -1,4 +1,4 @@
-function [Volt,UAngel,Init_Z,Init_W,Init_L,Init_U,Init_Y,PG,QG,RestraintCount,wPG,wQG,wD,PD,PD0,QD,randPDind,Loadi]=OPF_Init(Busnum,Balance,PG,QG,Volt,GenU,GenL,PVi,PGi,PVQU,PVQL,PD0,QD0,QD,PD)
+function [Volt,UAngel,Init_Z,Init_W,Init_L,Init_U,Init_Y,PG,QG,RestraintCount,wPG,wQG,wPD,wQD,PD,PD0,QD,randPDind,Loadi]=OPF_Init(Busnum,Balance,PG,QG,Volt,GenU,GenL,PVi,PGi,PVQU,PVQL,PD0,QD0,QD,PD)
 Loadi=find(QD~=0 | PD~=0);
 %Loadi=[1:Busnum]';
 RestraintCount=size(PVi,1)+size(PGi,1)+size(Loadi,1)*2+Busnum*1; %约束条件数,放开所有QD
@@ -22,7 +22,8 @@ wQG=1*ones(size(PVi,1),1);
 %randInt=randperm(size(Loadi,1));
 %randPDind=randInt(1:10);
 randPDind=0;
-wD=1*ones(size(Loadi,1),1);
+wPD=1*ones(size(Loadi,1),1);
+wQD=1*ones(size(Loadi,1),1);
 %wD(randPDind)=0;%一些负荷不约束
 %wD(7)=0;
 % wD(11)=0;

ObjectiveFun.m

@@ -1,6 +1,9 @@
-function [out_arg]=ObjectiveFun(PG,PG0,PGi,PD,PD0,wG,wD,Loadi)
+function [out_arg]=ObjectiveFun(PG,PG0,PGi,QG,QG0,PVi,PD,PD0,QD,QD0,wPG,wQG,wPD,wQD,Loadi)
 t1=PG(PGi)-PG0(PGi);
-t2=wG.*t1.*t1;
+t2=wPG.*t1.*t1;
-t3=wD.*((PD(Loadi)-PD0(Loadi)).*(PD(Loadi)-PD0(Loadi)));
+t1=QG(PVi)-QG0(PVi);
-out_arg= sum(t2)+sum(t3);
+t3=wQG.*t1.*t1;
+t4=wPD.*((PD(Loadi)-PD0(Loadi)).*(PD(Loadi)-PD0(Loadi)));
+t5=wQD.*((QD(Loadi)-QD0(Loadi)).*(QD(Loadi)-QD0(Loadi)));
+out_arg= sum(t2)+sum(t3)+sum(t4)+sum(t5);
 end

func_deltF.m

@@ -1,8 +1,8 @@
-function deltF=func_deltF(PG,QG,PVi,PGi,wPG,wQG,wD,PG0,QG0,PD0,PD,Busnum,Loadi)
+function deltF=func_deltF(PG,QG,PVi,PGi,wPG,wQG,wPD,PG0,QG0,PD0,PD,Busnum,Loadi)
 
 t1=2*wPG.*(PG(PGi) - PG0(PGi) );
 t2=2*wQG.*(QG(PVi) - QG0(PVi) );
-t3=2*wD.*(PD(Loadi)-PD0(Loadi));
+t3=2*wPD.*(PD(Loadi)-PD0(Loadi));
 deltF=[
     sparse(t1);
     sparse(t2);

func_deltdeltF.m

@@ -1,7 +1,7 @@
-function deltdeltF=func_deltdeltF(PVi,wPG,wQG,wD,ContrlCount)
+function deltdeltF=func_deltdeltF(PVi,wPG,wQG,wPD,ContrlCount)
 
 %ContrlCount=size(PVi,1)+size(PGi,1)+size(wD,1)+Busnum*2; %P,Q,Volt theta这些控制变量数
-C=[wPG' wQG' wD'];
+C=[wPG' wQG' wPD'];
 sizeC=size(C,2);
 diagC=sparse(1:sizeC,1:sizeC,C,sizeC,sizeC);
 deltdeltF=[

pf.m

@@ -31,4 +31,8 @@ for i = 0:kmax
     end
 end
 toc;
-end
+PG=PG/PQstandard;
+QG=QG/PQstandard;
+PD=PD/PQstandard;
+QD=QD/PQstandard;
+end