删掉多余的函数参数，统一了ContrlCount

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

FormG.m

@@ -1,21 +1,9 @@
 function Mat_G=FormG(Volt,PVi,PGi,PG,QG,PD,Loadi)
-%t1=PG(PVi);
+
-%GP=t1;%发电机P
-%GP=[4.5 4.5]';
-%%线路
-%发电机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);
     QG(PVi);
     sparse(PD(Loadi));
-    %GQ(PVi);
-    %[0 1.45]';
     Volt';
     ];
 end

FormH.m

@@ -1,4 +1,4 @@
-function Mat_H=FormH(Busnum,GB,AngleIJMat,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)
 %% 
 QDcos=textread('1047glys.txt');
 QD(QD~=0)=PD(QD~=0)./tan(QDcos);

FormLw.m

@@ -1,4 +1,4 @@
-function Lw=FormLw(Mat_G,Init_U,GenU,Busnum,PVQU,PD0,randPDind,Loadi)
+function Lw=FormLw(Mat_G,Init_U,GenU,Busnum,PVQU,PD0,Loadi)
 
 PU=1*GenU(:,2);%发电机有功上界
 QU=PVQU(:,1);%发电机无功上界

FormLz.m

@@ -1,4 +1,4 @@
-function Lz=FormLz(Mat_G,Init_L,GenL,Busnum,PVQL,PD0,randPDind,Loadi)
+function Lz=FormLz(Mat_G,Init_L,GenL,Busnum,PVQL,PD0,Loadi)
 
 PL=1*GenL(:,2);%发电机有功下界
 QL=PVQL(:,1);%发电机无功下界

OPF.asv

@@ -1,6 +1,6 @@
 tic
 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]=pf('ieee3001PG.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');
 %PVi电压节点序号
 %PVu电压节点电压标幺值
 Volt;
@@ -8,11 +8,11 @@ UAngel*180/3.1415926;
 %% 初值
 PG0=PG;
 PD0=PD;
-[Volt,UAngel,Init_Z,Init_W,Init_L,Init_U,Init_Y,PG,QG,RestraintCount,wG,wD,PD,randPDind]=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,PD0,randPDind,Loadi]=OPF_Init(Busnum,Balance,PG,QG,Volt,GenU,GenL,PVi,PGi,PVQU,PVQL,PD0,QD);
 Gap=(Init_L*Init_Z'-Init_U*Init_W');
 KK=0;
 plotGap=zeros(1,50);
-ContrlCount=size(PVi,1)+size(PGi,1)+Busnum*3;
+ContrlCount=size(PVi,1)+size(PGi,1)+size(Loadi,1)+Busnum*2;
 kmax=600;
 %% 20120523 临时
 QD_NON_ZERO=QD(PD==0 & QD~=0);
@@ -27,42 +27,44 @@ while(abs(Gap)>Precision)
     AngleIJMat=0;
     %% 开始计算OPF
     %% 形成等式约束的雅克比
-    deltH=func_deltH(Busnum,Volt,PVi,AngleIJMat,Y,GB,PGi,UAngel,r,c,Angle);
+    deltH=func_deltH(Busnum,Volt,PVi,Y,PGi,UAngel,r,c,Angle,Loadi);
     %% 形成不等式约束的雅克比
-    deltG=func_deltG(Busnum,PVi,PGi);
+    deltG=func_deltG(Busnum,PVi,PGi,Loadi);
     %%
     L_1Z=diag(Init_Z./Init_L);
     U_1W=diag(Init_W./Init_U);
     %% 形成海森阵
-    deltdeltF=func_deltdeltF(Busnum,GenC,PVi,PGi,wG,wD,PD0);
+    deltdeltF=func_deltdeltF(Busnum,PVi,PGi,wG,wD,ContrlCount);
     %% 形成ddHy
-    ddh=func_ddh3(AngleIJMat,GB,Volt,Init_Y,Busnum,PVi,PGi,Y,UAngel,r,c,Angle);
+    ddh=func_ddh3(AngleIJMat,GB,Volt,Init_Y,Busnum,PVi,PGi,Y,UAngel,r,c,Angle,Loadi);
     %%  开始构建ddg
-    ddg=func_ddg(PGi,PVi,Busnum,RestraintCount);
+    ddg=func_ddg(PGi,PVi,Busnum,RestraintCount,Loadi);
     %%  开始构建deltF
-    deltF=func_deltF(PG,PVi,GenC,ContrlCount,PGi,wG,wD,PG0,PD0,PD,Busnum);
+    deltF=func_deltF(PG,PVi,GenC,ContrlCount,PGi,wG,wD,PG0,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);
+    Mat_G=FormG(Volt,PVi,PGi,PG,QG,PD,Loadi);
     Mat_H=FormH(Busnum,GB,AngleIJMat,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);
+    Lz=FormLz(Mat_G,Init_L,GenL,Busnum,PVQL,PD0,randPDind,Loadi);
-    Lw=FormLw(Mat_G,Init_U,GenU,Busnum,PVQU,PD0);
+    Lw=FormLw(Mat_G,Init_U,GenU,Busnum,PVQU,PD0,randPDind,Loadi);
     Lx=FormLx(deltF,deltH,Init_Y,deltG,Init_Z,Init_W);
     YY=FormYY1(Lul,Lz,Ly,Luu,Lw,Lx);
     %% 开始解方程
-    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);
+    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]=AssignXX1(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);
+    [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);
     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)
+ObjectiveFun(PG,PG0,PGi,PD,PD0,wG,wD,Loadi)
 DrawGap(plotGap);
-Volt=Volt';
+%% 
+Volt=full(Volt');
+PD=full(PD);
 toc  
 

OPF.m

@@ -27,29 +27,29 @@ while(abs(Gap)>Precision)
     AngleIJMat=0;
     %% 开始计算OPF
     %% 形成等式约束的雅克比
-    deltH=func_deltH(Busnum,Volt,PVi,AngleIJMat,Y,GB,PGi,UAngel,r,c,Angle,Loadi);
+    deltH=func_deltH(Busnum,Volt,PVi,Y,PGi,UAngel,r,c,Angle,Loadi);
     %% 形成不等式约束的雅克比
     deltG=func_deltG(Busnum,PVi,PGi,Loadi);
     %%
     L_1Z=diag(Init_Z./Init_L);
     U_1W=diag(Init_W./Init_U);
     %% 形成海森阵
-    deltdeltF=func_deltdeltF(Busnum,GenC,PVi,PGi,wG,wD,PD0);
+    deltdeltF=func_deltdeltF(PVi,wG,wD,ContrlCount);
     %% 形成ddHy
-    ddh=func_ddh3(AngleIJMat,GB,Volt,Init_Y,Busnum,PVi,PGi,Y,UAngel,r,c,Angle,Loadi);
+    ddh=func_ddh3(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,GenC,ContrlCount,PGi,wG,wD,PG0,PD0,PD,Busnum,Loadi);
+    deltF=func_deltF(PG,PVi,PGi,wG,wD,PG0,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_H=FormH(Busnum,GB,AngleIJMat,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);
     Ly=Mat_H;
-    Lz=FormLz(Mat_G,Init_L,GenL,Busnum,PVQL,PD0,randPDind,Loadi);
+    Lz=FormLz(Mat_G,Init_L,GenL,Busnum,PVQL,PD0,Loadi);
-    Lw=FormLw(Mat_G,Init_U,GenU,Busnum,PVQU,PD0,randPDind,Loadi);
+    Lw=FormLw(Mat_G,Init_U,GenU,Busnum,PVQU,PD0,Loadi);
     Lx=FormLx(deltF,deltH,Init_Y,deltG,Init_Z,Init_W);
     YY=FormYY1(Lul,Lz,Ly,Luu,Lw,Lx);
     %% 开始解方程

OPF_Init.m

@@ -19,16 +19,10 @@ PG(PGi)=(tPU+tPL)/2;
 QG(PVi)=(tQU+tQL)/2;
 wG=ones(size(PGi,1),1);
 randInt=randperm(size(Loadi,1));
-%asc_randInt=sort(randInt,2,'ascend' );
+randPDind=randInt(1:262);
-%randPDind=randInt(1:200);
-
-randPDind=randInt(1:0);
 
 wD=ones(size(Loadi,1),1);
 wD(randPDind)=0;%一些负荷不约束
-%wD(Balance)=0;
-%PD0(randPDind)=0;
 PD=1*PD0;
-%PD(PD==0)=.2;
 
 end

func_ddh3.m

@@ -1,6 +1,6 @@
-function ddh=func_ddh3(AngleIJMat,GB,Volt,Init_Y,Busnum,PVi,PGi,Y,UAngel,r,c,Angle,Loadi)
+function ddh=func_ddh3(Volt,Init_Y,Busnum,PVi,PGi,Y,UAngel,r,c,Angle,Loadi,ContrlCount)
 %决定用循环重写
-ContrlCount=size(PVi,1)+size(PGi,1)+size(Loadi,1)+Busnum*2;
+%ContrlCount=size(PVi,1)+size(PGi,1)+size(Loadi,1)+Busnum*2;
 
 %AngleIJ=AngleIJMat-angle(GB);
 mat_AngleIJ=sparse(r,c,UAngel(r)-UAngel(c)-Angle',Busnum,Busnum);

func_deltF.m

@@ -1,12 +1,5 @@
-function deltF=func_deltF(PG,PVi,GenC,ContrlCount,PGi,wG,wD,PG0,PD0,PD,Busnum,Loadi)
+function deltF=func_deltF(PG,PVi,PGi,wG,wD,PG0,PD0,PD,Busnum,Loadi)
-%t1=PG(setdiff(PVi,Balance));
+
-% t2=Volt'*Volt;
-% t3=real(GB).*cos(AngleIJMat)+imag(GB).*sin(AngleIJMat);
-% t4=t2.*t3;
-% t5=sum(t4,2);
-% PBal=t5(Balance);
-% PPG=([PQ(1),PBal])';%暂时用土办法处理一下
-%% 
 t1=2*wG.*(PG(PGi) - PG0(PGi) );
 t2=2*wD.*(PD(Loadi)-PD0(Loadi));
 deltF=[

func_deltH.asv

@@ -1,20 +1,8 @@
-function deltH=func_deltH(Busnum,Volt,PVi,AngleIJMat,Y,GB,PGi)
+function deltH=func_deltH(Busnum,Volt,PVi,Y,PGi,UAngel,r,c,Angle,Loadi)
-dH_dPg=zeros(size(PGi,1)+size(PVi,1),2*Busnum);
 
-for I=1:size(PGi,1)
+dH_dPg=sparse(1:size(PGi,1),PGi,ones(size(PGi,1),1),size(PGi,1),2*Busnum);
-    %dH_dPg(I,2*PVi(I)-1)=-1;王锡凡书上的公式
+dH_dQr=sparse(1:size(PVi,1),PVi+Busnum,ones(size(PVi,1),1),size(PVi,1),2*Busnum);
-    %dH_dPg(I,2*PGi(I)-1)=1;暂时改一下20111227
+dH_dPD=[sparse(1:size(Loadi,1),Loadi,-ones(size(Loadi,1),1),size(Loadi,1),Busnum) sparse(size(Loadi,1),Busnum)];
-    dH_dPg(I,PGi(I))=1;
+dH_dx = jacobian_M3(Busnum,Volt,Y,Angle,UAngel,r,c);  %形成雅克比矩阵
-end
+deltH=[dH_dPg;dH_dQr;dH_dPD;dH_dx'];
-
-dH_dQr=zeros(size(PVi,1)+size(PGi,1),2*Busnum);
-for I=1:size(PVi,1)
-    %dH_dQr(I,2*PVi(I))=-1;王锡凡书上的公式
-    %dH_dQr(I,2*PVi(I))=1;暂时改一下20111227
-    dH_dQr(I,PVi(I))=1;
-end
-Angle=angle(GB);
-dH_dx = jacobian_M3(Busnum,Volt,Y,Angle,AngleIJMat);  %形成雅克比矩阵
-%deltH=[dH_dPg;dH_dQr;dH_dx'];%dH_dx 需要使用一下转置 暂时改一下
-deltH=[dH_dPg;dH_dQr;dH_dx'];
 end

func_deltH.m

@@ -1,4 +1,4 @@
-function deltH=func_deltH(Busnum,Volt,PVi,AngleIJMat,Y,GB,PGi,UAngel,r,c,Angle,Loadi)
+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);

func_deltdeltF.asv

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

func_deltdeltF.m

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