只加一个QD的收敛了，准备加所有QD。

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

FormG.asv

@@ -1,26 +1,11 @@
-function Mat_G=FormG(Volt,PVi,PG,GB,AngleIJMat,indexi,indexj)
+function Mat_G=FormG(Volt,PVi,PGi,PG,QG,PD,QD,Loadi)
-t1=PG(PVi);
+
-%GP=t1;%发电机P
-GP=[4.5 4.5]';
-%%线路
-t1=Volt'*Volt;
-t2=real(GB).*cos(AngleIJMat)+imag(GB).*sin(AngleIJMat);
-t3=t1.*t2;
-t4=zeros(size(indexi,2),1);
-for I=1:size(indexi,2)
-    t4(I)=Volt(indexi(I))*real(GB(indexi(I),indexj(I)))+t3(indexi(I),indexj(I));%线路
-end
-LP=t4;
-%发电机Q
-t1=Volt'*Volt;
-t2=real(GB).*sin(AngleIJMat)-imag(GB).*cos(AngleIJMat);
-t3=t1.*t2;
-t4=sum(t3,2);%发电机Q
-GQ=t4;
 Mat_G=[
-    GP;
+    PG(PGi);
-    GQ(PVi);
+    QG(PVi);
+    sparse(PD(Loadi));
+    sparse(QD(Loadi));
     Volt';
-    LP;
+    
     ];
 end

FormG.m

@@ -4,7 +4,7 @@ Mat_G=[
     PG(PGi);
     QG(PVi);
     sparse(PD(Loadi));
-    sparse(QD(Loadi));
+    sparse(QD(Loadi(1)));
     Volt';
     ];
 end

FormH.m

@@ -1,9 +1,10 @@
-function Mat_H=FormH(Busnum,Volt,PG,PD,QG,QD,Y,UAngel,r,c,Angle,QD_NON_ZERO,QD_NON_ZERO_IND)
+function Mat_H=FormH(Busnum,Volt,PG,PD,QG,QD,Y,UAngel,r,c,Angle,QD_NON_ZERO,QD_NON_ZERO_IND,Loadi)
 %% 
 %QDcos=textread('D:\Project\×îС»¯³±Á÷\×îС³±Á÷ËãÀý\Ïɺ£919PDQDglys.txt');
 %QD(QD~=0)=PD(QD~=0)./tan(QDcos);
 %QD(QD_NON_ZERO_IND)=QD_NON_ZERO;
 %% 
+%PD(Loadi)=QD(Loadi)./tan(acos(0.98));
 AngleIJ=sparse(r,c,UAngel(r)-UAngel(c)-Angle',Busnum,Busnum);
 dP=PG-PD-diag(Volt)*Y.*cos(AngleIJ)*Volt';
 dQ=QG-QD-diag(Volt)*Y.*sin(AngleIJ)*Volt';

FormLw.m

@@ -9,10 +9,10 @@ PDU(PDU>0)=1.200*PDU(PDU>0);
 PDU(PDU<0)=0.800*PDU(PDU<0);
 PDU(PDU==0)=0.400;
 %PDU=10*ones(length(Loadi),1);
-QDU=QD0(Loadi);
+QDU=QD0(Loadi(1));
 QDU(QDU>0)=1.200*QDU(QDU>0);
 QDU(QDU<0)=0.800*QDU(QDU<0);
-QDU(QDU==0)=0.400;
+QDU(QDU==0)=0.200;
 t1=([PU',QU',PDU',QDU',VoltU])';
 t2=Mat_G+Init_U'-t1;
 Lw=t2;

FormLz.m

@@ -9,10 +9,13 @@ PDL(PDL>0)=0.800*PDL(PDL>0);
 PDL(PDL<0)=1.200*PDL(PDL<0);
 PDL(PDL==0)=-0.400;
 %PDL=-10*ones(length(Loadi),1);
-QDL=QD0(Loadi);
+QDL=QD0(Loadi(1));
 QDL(QDL>0)=0.800*QDL(QDL>0);
 QDL(QDL<0)=1.200*QDL(QDL<0);
-QDL(QDL==0)=-0.400;
+QDL(QDL==0)=-0.200;
+% QDL(QDL>0)=1;
+% QDL(QDL<0)=-1;
+% QDL(QDL==0)=-0.200;
 t1=([PL',QL',PDL',QDL',VoltL])';
 t2=Mat_G-Init_L'-t1;
 Lz=t2;

Modification.m

@@ -13,10 +13,10 @@ Init_Y=Init_Y+AlphaD*deltY';
 PG(PGi)=PG(PGi)+AlphaP*deltX(1:size(PGi,1));
 %QG(PVi)=QG(PVi)+deltX(size(PGi,1)+1:size(PVi,1)+size(PGi,1) );
 QG(PVi)=QG(PVi)+AlphaP*deltX(size(PGi,1)+1:size(PVi,1)+size(PGi,1) );
-t=deltX(size(PVi,1)+size(PGi,1)+1:size(PVi,1)+size(PGi,1)+size(Loadi,1)*2);
+t=deltX(size(PVi,1)+size(PGi,1)+1:size(PVi,1)+size(PGi,1)+size(Loadi,1)*1+1);
 PD(Loadi)=PD(Loadi)+AlphaP*t(1:length(Loadi));
-QD(Loadi)=QD(Loadi)+AlphaP*t(length(Loadi)+1:2*length(Loadi));
+QD(Loadi(1))=QD(Loadi(1))+AlphaP*t(length(Loadi)+1:length(Loadi)+1);
-t=deltX(size(PVi,1)+size(PGi,1)+size(Loadi,1)*2+1:ContrlCount)';
+t=deltX(size(PVi,1)+size(PGi,1)+size(Loadi,1)*1+2:ContrlCount)';
 t(Busnum+Balance)=0;
 %Volt=Volt+AlphaP*t(2:2:2*Busnum);暂时改一下20111227
 %UAngel=UAngel+AlphaP*t(1:2:2*Busnum);暂时改一下20111227

OPF.asv

@@ -1,11 +1,12 @@
 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]= ...
+[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('D:\Project\×îС»¯³±Á÷\×îС³±Á÷ËãÀý\ԭʼ\standard.txt'); 
+pf('ieee4.dat');
+%pf('D:\Project\青秀降损项目\最小化潮流\最小潮流算例\原始\津头站津视922（3-1）_0.5_120%.txt'); 
 %pf('D:\Project\最小化潮流\最小潮流算例\仙海919.txt'); 
 %pf('c:/file31.txt'); 
-%pf('ieee118PG.dat');
+
 
 %% 计算功率因数
 %atan(PD(QD~=0 | PD~=0)./QD(QD~=0 | PD~=0));
@@ -14,21 +15,23 @@ UAngel*180/3.1415926;
 %% 通过潮流计算PG
 AngleIJ=sparse(r,c,UAngel(r)-UAngel(c)-Angle',Busnum,Busnum);
 PGBal=diag(Volt)*Y.*cos(AngleIJ)*Volt';
-
+QGBal=diag(Volt)*Y.*sin(AngleIJ)*Volt';
 %% 初值-即测量值
 PG0=PG;
+QG0=QG;
 PD0=PD;
-
+QD0=QD;
 PDReal=PD;%真值
 %PD0(12)=PD0(12)+0.001;
 %%
 PG0(Balance)=PGBal(Balance);
+QG0(Balance)=QGBal(Balance);
 %%
-[Volt,UAngel,Init_Z,Init_W,Init_L,Init_U,Init_Y,PG,QG,RestraintCount,wG,wD,PD,PD0,randPDind,Loadi]=OPF_Init(Busnum,Balance,PG,QG,Volt,GenU,GenL,PVi,PGi,PVQU,PVQL,PD0,QD,PD);
+[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);
 Gap=(Init_L*Init_Z'-Init_U*Init_W');
 KK=0;
 plotGap=zeros(1,50);
-ContrlCount=size(PVi,1)+size(PGi,1)+size(Loadi,1)+Busnum*2;
+ContrlCount=size(PVi,1)+size(PGi,1)+size(Loadi,1)*2+Busnum*2;
 kmax=60;
 %% 20120523 临时
 QD_NON_ZERO=QD(PD==0 & QD~=0);
@@ -50,34 +53,35 @@ while(abs(Gap)>Precision)
     L_1Z=diag(Init_Z./Init_L);
     U_1W=diag(Init_W./Init_U);
     %% 形成海森阵
-    deltdeltF=func_deltdeltF(PVi,wG,wD,ContrlCount);
+    deltdeltF=func_deltdeltF(PVi,wPG,wQG,wD,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);
     %%  开始构建deltF
-    deltF=func_deltF(PG,PVi,PGi,wG,wD,PG0,PD0,PD,Busnum,Loadi);
+    deltF=func_deltF(PG,QG,PVi,PGi,wPG,wQG,wD,PG0,QG0,PD0,PD,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,PVi,PGi,PG,QG,PD,Loadi);
+    Mat_G=FormG(Volt,PVi,PGi,PG,QG,PD,QD,Loadi);
     Mat_H=FormH(Busnum,Volt,PG,PD,QG,QD,Y,UAngel,r,c,Angle,QD_NON_ZERO,QD_NON_ZERO_IND);
     Ly=Mat_H;
-    Lz=FormLz(Mat_G,Init_L,GenL,Busnum,PVQL,PD0,Loadi,KK);
+    Lz=FormLz(Mat_G,Init_L,GenL,Busnum,PVQL,PD0,QD0,Loadi,KK);
-    Lw=FormLw(Mat_G,Init_U,GenU,Busnum,PVQU,PD0,Loadi,KK);
+    Lw=FormLw(Mat_G,Init_U,GenU,Busnum,PVQU,PD0,QD0,Loadi,KK);
     Lx=FormLx(deltF,deltH,Init_Y,deltG,Init_Z,Init_W);
     YY=FormYY(Lul,Lz,Ly,Luu,Lw,Lx);
     %% 开始解方程
+    fprintf(')
     XX=SolveIt(deltF,deltG,Init_L,Init_Z,Init_U,Init_W,deltdeltF,ddh,ddg,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]=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,Loadi);
+    [Init_Z,Init_L,Init_W,Init_U,Init_Y,PG,QG,Volt,UAngel,PD,QD]=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);
     Gap=(Init_L*Init_Z'-Init_U*Init_W');
     KK=KK+1;
 end
 fprintf('迭代次数%d\n',KK);
-ObjectiveFun(PG,PG0,PGi,PD,PD0,wG,wD,Loadi)
+%fprintf('目标值%f\n',full(ObjectiveFun(PG,PG0,PGi,PD,PD0,wG,wD,Loadi)));
 DrawGap(plotGap);
 %% 
 %Volt=full(Volt');
@@ -88,7 +92,10 @@ absPDLoad=abs(  (PD(Loadi)-PDReal(Loadi))./PDReal(Loadi)  );
 maxPDError=max(absPDLoad);
 disp('index');
 Loadi(absPDLoad==maxPDError);
-%% ¼ÆËãÏßËð
+%% 计算总线损
-Lineloss(Linei,Linej,Liner,Linex,Lineb2,Transfori,Transj,Transx,k0,Volt,Angle)
+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);
 toc  
 

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('ieee4PG.dat');
+pf('ieee3001.dat');
 %pf('D:\Project\青秀降损项目\最小化潮流\最小潮流算例\原始\津头站津视922（3-1）_0.5_120%.txt'); 
 %pf('D:\Project\最小化潮流\最小潮流算例\仙海919.txt'); 
 %pf('c:/file31.txt'); 
@@ -15,21 +15,23 @@ UAngel*180/3.1415926;
 %% 通过潮流计算PG
 AngleIJ=sparse(r,c,UAngel(r)-UAngel(c)-Angle',Busnum,Busnum);
 PGBal=diag(Volt)*Y.*cos(AngleIJ)*Volt';
-
+QGBal=diag(Volt)*Y.*sin(AngleIJ)*Volt';
 %% 初值-即测量值
 PG0=PG;
+QG0=QG;
 PD0=PD;
 QD0=QD;
 PDReal=PD;%真值
 %PD0(12)=PD0(12)+0.001;
 %%
 PG0(Balance)=PGBal(Balance);
+QG0(Balance)=QGBal(Balance);
 %%
-[Volt,UAngel,Init_Z,Init_W,Init_L,Init_U,Init_Y,PG,QG,RestraintCount,wG,wD,PD,PD0,randPDind,Loadi]=OPF_Init(Busnum,Balance,PG,QG,Volt,GenU,GenL,PVi,PGi,PVQU,PVQL,PD0,QD,PD);
+[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);
 Gap=(Init_L*Init_Z'-Init_U*Init_W');
 KK=0;
-plotGap=zeros(1,50);
+plotGap=zeros(1,60);
-ContrlCount=size(PVi,1)+size(PGi,1)+size(Loadi,1)*2+Busnum*2;
+ContrlCount=size(PVi,1)+size(PGi,1)+size(Loadi,1)*1+Busnum*2+1;
 kmax=60;
 %% 20120523 临时
 QD_NON_ZERO=QD(PD==0 & QD~=0);
@@ -46,30 +48,31 @@ while(abs(Gap)>Precision)
     %% 形成等式约束的雅克比
     deltH=func_deltH(Busnum,Volt,PVi,Y,PGi,UAngel,r,c,Angle,Loadi);
     %% 形成不等式约束的雅克比
-    deltG=func_deltG(Busnum,PVi,PGi,Loadi);
+    deltG=func_deltG(Busnum,PVi,PGi,Loadi,PD,QD);
     %%
     L_1Z=diag(Init_Z./Init_L);
     U_1W=diag(Init_W./Init_U);
     %% 形成海森阵
-    deltdeltF=func_deltdeltF(PVi,wG,wD,ContrlCount);
+    deltdeltF=func_deltdeltF(PVi,wPG,wQG,wD,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);
+    ddg=func_ddg(PGi,PVi,Busnum,RestraintCount,Loadi,PD,QD);
     %%  开始构建deltF
-    deltF=func_deltF(PG,PVi,PGi,wG,wD,PG0,PD0,PD,Busnum,Loadi);
+    deltF=func_deltF(PG,QG,PVi,PGi,wPG,wQG,wD,PG0,QG0,PD0,PD,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,PVi,PGi,PG,QG,PD,QD,Loadi);
-    Mat_H=FormH(Busnum,Volt,PG,PD,QG,QD,Y,UAngel,r,c,Angle,QD_NON_ZERO,QD_NON_ZERO_IND);
+    Mat_H=FormH(Busnum,Volt,PG,PD,QG,QD,Y,UAngel,r,c,Angle,QD_NON_ZERO,QD_NON_ZERO_IND,Loadi);
     Ly=Mat_H;
     Lz=FormLz(Mat_G,Init_L,GenL,Busnum,PVQL,PD0,QD0,Loadi,KK);
     Lw=FormLw(Mat_G,Init_U,GenU,Busnum,PVQU,PD0,QD0,Loadi,KK);
     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,ddg,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);
@@ -78,7 +81,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,PD,PD0,wG,wD,Loadi)));
 DrawGap(plotGap);
 %% 
 %Volt=full(Volt');

OPF_Init.asv

@@ -1,13 +1,15 @@
-function [Volt,UAngel,Init_Z,Init_W,Init_L,Init_U,Init_Y,PG,QG,RestraintCount,wG,wD,PD,PD0,randPDind]=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,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)
-RestraintCount=size(PVi,1)+size(PGi,1)+Busnum*2; %约束条件数
+Loadi=find(QD~=0 | PD~=0);
+%Loadi=[1:Busnum]';
+RestraintCount=size(PVi,1)+size(PGi,1)+size(Loadi,1)*1+Busnum*1+1; %约束条件数,只放开一个QD
 t_Bal_volt=Volt(Balance);
-Volt=sparse(ones(1,Busnum));
+Volt=sparse(1*ones(1,Busnum));
 Volt(Balance)=t_Bal_volt;
 UAngel=sparse(1,Busnum);
 Init_Z=sparse(ones(1,RestraintCount));
 Init_W=sparse(-1*ones(1,RestraintCount));
-Init_L=sparse(ones(1,RestraintCount));
+Init_L=1*sparse(ones(1,RestraintCount));
-Init_U=sparse(ones(1,RestraintCount));
+Init_U=1*sparse(ones(1,RestraintCount));
 Init_Y=sparse(1,2*Busnum);%与学姐一致
 tPU=sparse(GenU(:,2));% 发电机有功上限
 tQU=sparse(PVQU(:,1));% 无功上限
@@ -15,15 +17,17 @@ tPL=sparse(GenL(:,2));%
 tQL=sparse(PVQL(:,1));% 无功下限
 PG(PGi)=(tPU+tPL)/2;
 QG(PVi)=(tQU+tQL)/2;
-wG=ones(size(PGi,1),1);
+wPG=1*ones(size(PGi,1),1);
-randInt=randperm(300);
+wQG=1*zeros(size(PVi,1),1);
-asc_randInt=sort(randInt,2,'ascend' );
+%randInt=randperm(size(Loadi,1));
-randPDind=randInt(1:200);
+%randPDind=randInt(1:10);
-wD=ones(Busnum,1);
+randPDind=0;
-wD(randPDind)=0;%一些负荷不约束
+wD=1*zeros(size(Loadi,1),1);
-%wD(Balance)=0;
+%wD(randPDind)=0;%一些负荷不约束
-PD0(randPDind)=0;
+%wD(7)=0;
+% wD(11)=0;
 PD=1*PD0;
-%PD(PD==0)=.2;
+%powerFacter=0.98;
-
+%QD=PD*sqrt((1-powerFacter^2)/powerFacter^2);
+QD=QD0;
 end

OPF_Init.m

@@ -1,7 +1,7 @@
-function [Volt,UAngel,Init_Z,Init_W,Init_L,Init_U,Init_Y,PG,QG,RestraintCount,wG,wD,PD,PD0,randPDind,Loadi]=OPF_Init(Busnum,Balance,PG,QG,Volt,GenU,GenL,PVi,PGi,PVQU,PVQL,PD0,QD,PD)
+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)
 Loadi=find(QD~=0 | PD~=0);
 %Loadi=[1:Busnum]';
-RestraintCount=size(PVi,1)+size(PGi,1)+size(Loadi,1)*2+Busnum*1; %约束条件数
+RestraintCount=size(PVi,1)+size(PGi,1)+size(Loadi,1)*1+Busnum*1+1; %约束条件数,只放开一个QD
 t_Bal_volt=Volt(Balance);
 Volt=sparse(1*ones(1,Busnum));
 Volt(Balance)=t_Bal_volt;
@@ -17,7 +17,8 @@ tPL=sparse(GenL(:,2));%
 tQL=sparse(PVQL(:,1));% 无功下限
 PG(PGi)=(tPU+tPL)/2;
 QG(PVi)=(tQU+tQL)/2;
-wG=1*ones(size(PGi,1),1);
+wPG=1*ones(size(PGi,1),1);
+wQG=1*ones(size(PVi,1),1);
 %randInt=randperm(size(Loadi,1));
 %randPDind=randInt(1:10);
 randPDind=0;
@@ -26,5 +27,7 @@ wD=1*ones(size(Loadi,1),1);
 %wD(7)=0;
 % wD(11)=0;
 PD=1*PD0;
-
+%powerFacter=0.98;
+%QD=PD*sqrt((1-powerFacter^2)/powerFacter^2);
+QD=QD0;
 end

SolveIt.m

@@ -10,14 +10,14 @@ aa=[
     ];
 yy=[LxComa;-Ly];
 %% 平衡节点电压不变 
-t=size(PVi,1)+size(PGi,1)+size(Loadi,1);
+t=size(PVi,1)+size(PGi,1)+size(Loadi,1)+1;
 aa(t+Balance,:)=0;
 aa(:,t+Balance)=0;
 %aa(t+Balance,t+Balance)=1;
 aa=aa+sparse(t+Balance,t+Balance,ones(length(Balance),1),ContrlCount+2*Busnum,ContrlCount+2*Busnum);
 deltG(t+Balance,:)=0;
 %%
-t=size(PVi,1)+size(PGi,1)+size(Loadi,1)+Busnum*1;
+t=size(PVi,1)+size(PGi,1)+size(Loadi,1)+1+Busnum*1;
 aa(t+Balance,:)=0;
 aa(:,t+Balance)=0;
 %aa(t+Balance,t+Balance)=1;

func_ddg.asv

@@ -1,93 +1,9 @@
-function ddg=func_ddg(AngleIJMat,GB,Volt,Init_W,Init_Z,Busnum,indexi,indexj,PVi,RestraintCount)
+function ddg=func_ddg(PGi,PVi,Busnum,RestraintCount,Loadi,PD,QD)
-c=Init_W+Init_Z;
-%% dg4_dTdT 对角元素
-%t1=-Volt'*Volt;
-t1=Volt'*Volt;
-t2=real(GB).*cos(AngleIJMat)+imag(GB).*sin(AngleIJMat);
-t3=t1.*t2;
-t4=t3(sub2ind(size(t3),indexi,indexj));
-t5=zeros(Busnum);
-offset=2*size(PVi,1)+Busnum;
-for I=1:size(indexi,2)
-    t5(indexi(I),indexi(I))=t4(I)*c(offset+I);
-    t5(indexj(I),indexj(I))=t4(I)*c(offset+I);
-end
-dPdTidTi=t5; %@@@
-%% dg4_dTdT 非对角元素
-%t1=Volt'*Volt;
-t1=-Volt'*Volt;
-t2=real(GB).*cos(AngleIJMat)+imag(GB).*sin(AngleIJMat);
-t3=t1.*t2;
-t5=zeros(Busnum);
-for I=1:size(indexi,2)
-    t5(indexi(I),indexj(I))=t3(indexi(I),indexj(I))*c(offset+I);
-    t5(indexj(I),indexi(I))=t3(indexi(I),indexj(I))*c(offset+I);
-end
-dPdTidTj=t5;% @@@
-%% dg4_dVdV 对角元素
-t1=real(GB)*2;
-%dPijdVidVi=t1;
-t2=zeros(Busnum);
-for I=1:size(indexi,2)
-    t2(indexi(I),indexi(I))=t1(indexi(I),indexj(I))*c(offset+I);
-end
-dPdVidVi=t2; % @
-%% dg4_dVdV 非对角元素
-t2=real(GB).*cos(AngleIJMat)+imag(GB).*sin(AngleIJMat);
-t4=zeros(Busnum);
-for I=1:size(indexi,2)
-    %t4(indexi(I),indexj(I))=t2(indexi(I),indexj(I))*c(offset+I);
-    t4(indexi(I),indexj(I))=-t2(indexi(I),indexj(I))*c(offset+I);
-    %t4(indexj(I),indexi(I))=t2(indexi(I),indexj(I))*c(offset+I);
-    t4(indexj(I),indexi(I))=-t2(indexi(I),indexj(I))*c(offset+I);
-end
-dPdVidVj=t4; %@@
-%% dg4_dTdV 对角元素
-t2=real(GB).*sin(AngleIJMat)-imag(GB).*cos(AngleIJMat);
-t4=zeros(Busnum);
-for I=1:size(indexi,2)
-    %t4(indexi(I),indexi(I))=-Volt(indexj(I))*t2(indexi(I),indexi(I))*c(offset+I);
-    t4(indexi(I),indexi(I))=Volt(indexj(I))*t2(indexi(I),indexj(I))*c(offset+I);
-    %t4(indexj(I),indexj(I))=Volt(indexi(I))*t2(indexj(I),indexj(I))*c(offset+I);
-    t4(indexj(I),indexj(I))=-Volt(indexi(I))*t2(indexi(I),indexj(I))*c(offset+I);
-end
-dPdTidVi=t4; % @@
-%% dg4_dTdV 非对角元素
-t1=real(GB).*sin(AngleIJMat)-imag(GB).*cos(AngleIJMat);
-t2=zeros(Busnum);
-for I=1:size(indexi,2)
-    %t2(indexi(I),indexj(I))=-Volt(indexi(I))*t1(indexi(I),indexj(I))*c(offset+I);
-    t2(indexi(I),indexj(I))=Volt(indexi(I))*t1(indexi(I),indexj(I))*c(offset+I);
-    %t2(indexj(I),indexi(I))=Volt(indexj(I))*t1(indexi(I),indexj(I))*c(offset+I);
-    t2(indexj(I),indexi(I))=-Volt(indexj(I))*t1(indexi(I),indexj(I))*c(offset+I);
-end
-dPdTidVj=t2; % @@
-%% dg4_dVdT 对角元素
-dPdVidTi=dPdTidVi; 
-%% dg4_dVdT 非对角元素
-t1=real(GB).*sin(AngleIJMat)-imag(GB).*cos(AngleIJMat);
-t2=zeros(Busnum);
-for I=1:size(indexi,2)
-    %t2(indexi(I),indexj(I))=Volt(indexj(I))*t1(indexi(I),indexj(I))*c(offset+I);
-    t2(indexi(I),indexj(I))=-Volt(indexj(I))*t1(indexi(I),indexj(I))*c(offset+I);
-    %t2(indexj(I),indexi(I))=Volt(indexi(I))*t1(indexi(I),indexj(I))*c(offset+I);
-    t2(indexj(I),indexi(I))=Volt(indexi(I))*t1(indexi(I),indexj(I))*c(offset+I);
-end
-dPdVidTj=t2; % @
-%% 生成ddg4ddx
-ddg4ddx=zeros(2*Busnum,2*Busnum);
-ddg4ddx(1:2:2*Busnum,1:2:2*Busnum)=dPdTidTj;%%非对角 TT
-ddg4ddx(1:2:2*Busnum,2:2:2*Busnum)=dPdTidVj;%%非对角 TV
-ddg4ddx(2:2:2*Busnum,1:2:2*Busnum)=dPdVidTj;%%非对角 VT
-ddg4ddx(2:2:2*Busnum,2:2:2*Busnum)=dPdVidVj;%%非对角 VV
-ddg4ddx(1:2:2*Busnum,1:2:2*Busnum)=dPdTidTi;%%对角
-ddg4ddx(1:2:2*Busnum,2:2:2*Busnum)=dPdTidVi;%%对角
-ddg4ddx(2:2:2*Busnum,1:2:2*Busnum)=dPdVidTi;%%对角
-ddg4ddx(2:2:2*Busnum,2:2:2*Busnum)=dPdVidVi;%%对角
-%% 生成ddg
-t=[zeros(2*size(PVi,1),RestraintCount);
-    zeros(2*Busnum,RestraintCount-size(PVi,1)-2*Busnum),ddg4ddx;
-];
 
-ddg=t;
+ddg1=sparse(size(PVi,1)+size(PGi,1)+Loadi*2+Busnum,RestraintCount);
+t1=sparse(length(Loadi),size(PVi,1)+size(PGi,1)+Busnum)
+ddg2=sparse(length(Loadi),diag(-2./PD.^2-2*(QD.^2-PD.^2)/(PD.^2+QD.^2)));
+ddg3=diag(2*(PD.^2-QD.^2)/(PD.^2+QD.^2));
+ddg4=sparse(Busnum,RestraintCount);
+ddg=[ddg1;ddg2;ddg3;ddg4];
 end

func_ddg.m

@@ -1,5 +1,4 @@
-function ddg=func_ddg(PGi,PVi,Busnum,RestraintCount,Loadi)
+function ddg=func_ddg(PGi,PVi,Busnum,RestraintCount,Loadi,PD,QD)
-
-ddg=sparse(size(PVi,1)+size(PGi,1)+size(Loadi,1)*2+1*Busnum,RestraintCount);
 
+ddg=0;
 end

func_ddh.m

@@ -50,7 +50,7 @@ t=[ddPdVdV+ddQdVdV,ddPdVdT+ddQdVdT;
 ];
 sizePGi=size(PGi,1);
 sizePVi=size(PVi,1);
-sizeLoadi=size(Loadi,1)*2;
+sizeLoadi=size(Loadi,1)*1+1;
 ddh=[
     sparse(sizePGi+sizePVi+sizeLoadi,ContrlCount);
     sparse(2*Busnum,sizePVi+sizePGi+sizeLoadi),-t;

func_deltF.m

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

func_deltG.asv

@@ -1,30 +1,46 @@
-function deltG=func_deltG(Busnum,PVi,PGi)
+function deltG=func_deltG(Busnum,PVi,PGi,Loadi,PD,QD)
-dg1_dPg=eye(size(PGi,1));
-dg2_dPg=zeros(size(PGi,1),size(PVi,1));
-dg3_dPg=zeros(size(PGi,1),Busnum);
-dg4_dPg=zeros(size(PGi,1),Busnum);
 %%
-dg1_dQr=zeros(size(PVi,1),size(PGi,1));
+sizePGi=size(PGi,1);
-dg2_dQr=eye(size(PVi,1));
+sizePVi=size(PVi,1);
-dg3_dQr=zeros(size(PVi,1),Busnum);
+sizeLoadi=size(Loadi,1);
-dg4_dQr=zeros(size(PVi,1),Busnum);
 %%
-dg1_dPD=zeros(Busnum,size(PGi,1));
+dg1_dPg=sparse(1:sizePGi,1:sizePGi,ones(sizePGi,1),sizePGi,sizePGi);
-dg2_dPD=zeros(Busnum,size(PVi,1));
+dg2_dPg=sparse(sizePGi,sizePVi);
-dg3_dPD=zeros(Busnum,Busnum);
+dg3_dPg=sparse(sizePGi,sizeLoadi);
-dg4_dPD=zeros(Busnum,2*Busnum);
+dg4_dPg=sparse(sizePGi,sizeLoadi);
+dg5_dPg=sparse(sizePGi,Busnum);
 %% 
-dg1_dx=zeros(2*Busnum,size(PGi,1));
+dg1_dQr=sparse(sizePVi,sizePGi);
-dg2_dx=zeros(2*Busnum,size(PVi,1));
+dg2_dQr=sparse(1:sizePVi,1:sizePVi,ones(sizePVi,1),sizePVi,sizePVi);
-dg3_dx=zeros(2*Busnum,Busnum);
+dg3_dQr=sparse(sizePVi,sizeLoadi);
-dg4_dx=zeros(2*Busnum,Busnum);
+dg4_dQr=sparse(sizePVi,sizeLoadi);
-for I=1:Busnum
+dg5_dQr=sparse(sizePVi,Busnum);
-    %dg3_dx(2*I,I)=1;ÔÝĘą¸ÄŇťĎÂ
-    dg4_dx(I,I)=1;
-end
 %%
-deltG=[dg1_dPg,dg2_dPg,dg3_dPg;
+dg1_dPD=sparse(size(Loadi,1),size(PGi,1));
-       dg1_dQr,dg2_dQr,dg3_dQr;
+dg2_dPD=sparse(size(Loadi,1),size(PVi,1));
-       dg1_dPD,dg1_dPD,dg3_dPD,dg4_dPD;
+dg3_dPD=sparse(1:size(Loadi,1),1:size(Loadi,1),ones(size(Loadi,1),1),size(Loadi,1),size(Loadi,1));
-       dg1_dx,dg2_dx,dg3_dx;
+dg4_dPD=sparse(size(Loadi,1),sizeLoadi);
+dg5_dPD=sparse(size(Loadi,1),Busnum);
+%% 
+dg1_dQD=sparse(size(Loadi(1),1),size(PGi,1));
+dg2_dQD=sparse(size(Loadi(1),1),size(PVi,1));
+dg3_dQD=sparse(length(Loadi(1)),length(Loadi));
+dg4_dQD=sparse(1:size(Loadi(1),1),1:size(Loadi(1),1),ones(size(Loadi(1),1),1),size(Loadi(1),1),size(Loadi(1),1));
+dg5_dQD=sparse(size(Loadi(1),1),Busnum);
+
+%% 
+dg1_dx=sparse(2*Busnum,sizePGi);
+dg2_dx=sparse(2*Busnum,sizePVi);
+dg3_dx=sparse(2*Busnum,sizeLoadi);
+dg4_dx=sparse(2*Busnum,sizeLoadi);
+dg5_dx=[sparse(1:Busnum,1:Busnum,ones(Busnum,1),Busnum,Busnum);
+    sparse(Busnum,Busnum);
     ];
+%% 
+deltG=[dg1_dPg,dg2_dPg,dg3_dPg,dg4_dPg,dg5_dPg;
+       dg1_dQr,dg2_dQr,dg3_dQr,dg4_dQr,dg5_dQr;
+       dg1_dPD,dg2_dPD,dg3_dPD,dg4_dPD,dg5_dPD;
+       dg1_dQD,dg2_dQD,dg3_dQD,dg4_dQD,dg5_dQD,dg6_dQD;
+       dg1_dx,dg2_dx,dg3_dx,dg4_dx,dg5_dx,dg6_dx;
+];
+end

func_deltG.m

@@ -1,4 +1,4 @@
-function deltG=func_deltG(Busnum,PVi,PGi,Loadi)
+function deltG=func_deltG(Busnum,PVi,PGi,Loadi,PD,QD)
 %%
 sizePGi=size(PGi,1);
 sizePVi=size(PVi,1);
@@ -7,31 +7,32 @@ sizeLoadi=size(Loadi,1);
 dg1_dPg=sparse(1:sizePGi,1:sizePGi,ones(sizePGi,1),sizePGi,sizePGi);
 dg2_dPg=sparse(sizePGi,sizePVi);
 dg3_dPg=sparse(sizePGi,sizeLoadi);
-dg4_dPg=sparse(sizePGi,sizeLoadi);
+dg4_dPg=sparse(sizePGi,1);
 dg5_dPg=sparse(sizePGi,Busnum);
 %% 
 dg1_dQr=sparse(sizePVi,sizePGi);
 dg2_dQr=sparse(1:sizePVi,1:sizePVi,ones(sizePVi,1),sizePVi,sizePVi);
 dg3_dQr=sparse(sizePVi,sizeLoadi);
-dg4_dQr=sparse(sizePVi,sizeLoadi);
+dg4_dQr=sparse(sizePVi,1);
 dg5_dQr=sparse(sizePVi,Busnum);
 %%
 dg1_dPD=sparse(size(Loadi,1),size(PGi,1));
 dg2_dPD=sparse(size(Loadi,1),size(PVi,1));
 dg3_dPD=sparse(1:size(Loadi,1),1:size(Loadi,1),ones(size(Loadi,1),1),size(Loadi,1),size(Loadi,1));
-dg4_dPD=sparse(size(Loadi,1),sizeLoadi);
+dg4_dPD=sparse(size(Loadi,1),1);
 dg5_dPD=sparse(size(Loadi,1),Busnum);
 %% 
-dg1_dQD=sparse(size(Loadi,1),size(PGi,1));
+dg1_dQD=sparse(size(Loadi(1),1),size(PGi,1));
-dg2_dQD=sparse(size(Loadi,1),size(PVi,1));
+dg2_dQD=sparse(size(Loadi(1),1),size(PVi,1));
-dg3_dQD=sparse(length(Loadi),length(Loadi));
+dg3_dQD=sparse(length(Loadi(1)),length(Loadi));
-dg4_dQD=sparse(1:size(Loadi,1),1:size(Loadi,1),ones(size(Loadi,1),1),size(Loadi,1),size(Loadi,1));
+dg4_dQD=sparse(1:size(Loadi(1),1),1:size(Loadi(1),1),ones(size(Loadi(1),1),1),size(Loadi(1),1),size(Loadi(1),1));
-dg5_dQD=sparse(size(Loadi,1),Busnum);
+dg5_dQD=sparse(size(Loadi(1),1),Busnum);
+
 %% 
 dg1_dx=sparse(2*Busnum,sizePGi);
 dg2_dx=sparse(2*Busnum,sizePVi);
 dg3_dx=sparse(2*Busnum,sizeLoadi);
-dg4_dx=sparse(2*Busnum,sizeLoadi);
+dg4_dx=sparse(2*Busnum,1);
 dg5_dx=[sparse(1:Busnum,1:Busnum,ones(Busnum,1),Busnum,Busnum);
     sparse(Busnum,Busnum);
     ];
@@ -42,3 +43,4 @@ deltG=[dg1_dPg,dg2_dPg,dg3_dPg,dg4_dPg,dg5_dPg;
        dg1_dQD,dg2_dQD,dg3_dQD,dg4_dQD,dg5_dQD;
        dg1_dx,dg2_dx,dg3_dx,dg4_dx,dg5_dx;
 ];
+end

func_deltH.m

@@ -3,7 +3,8 @@ function deltH=func_deltH(Busnum,Volt,PVi,Y,PGi,UAngel,r,c,Angle,Loadi)
 dH_dPg=sparse(1:size(PGi,1),PGi,ones(size(PGi,1),1),size(PGi,1),2*Busnum);
 dH_dQr=sparse(1:size(PVi,1),PVi+Busnum,ones(size(PVi,1),1),size(PVi,1),2*Busnum);
 dH_dPD=[sparse(1:size(Loadi,1),Loadi,-ones(size(Loadi,1),1),size(Loadi,1),Busnum) sparse(size(Loadi,1),Busnum)];
-dH_dQD=[sparse(size(Loadi,1),Busnum) sparse(1:size(Loadi,1),Loadi,-ones(size(Loadi,1),1),size(Loadi,1),Busnum)];
+dH_dQD=[sparse(size(Loadi(1),1),Busnum) sparse(1:size(Loadi(1),1),Loadi(1),-ones(size(Loadi(1),1),1),size(Loadi(1),1),Busnum)];
 dH_dx = jacobian_M(Busnum,Volt,Y,Angle,UAngel,r,c);  %形成雅克比矩阵
 deltH=[dH_dPg;dH_dQr;dH_dPD;dH_dQD;dH_dx'];
+
 end

func_deltdeltF.m

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