收敛了，计算比原来用海森矩阵更慢

Signed-off-by: dmy@lab <dmy@lab.lab>

FormH.m

@@ -10,5 +10,5 @@ dP=PG-sparse(Loadi,1,PD,Busnum*3,1)-diag(Volt)*Y.*cos(AngleIJ)*Volt;
 dQ=QG-sparse(Loadi,1,QD,Busnum*3,1)-diag(Volt)*Y.*sin(AngleIJ)*Volt;
 
 Mat_H=[dP;dQ;];
-
+Mat_H=sparse(Mat_H);
 end

OPF.m

@@ -252,7 +252,6 @@ Precision=1e-5;
 CenterA=0.1;
 %% 加误差
 Volt=1*Vp3m;
-% Volt=mVoltABCV;
 UAngel=Vp3a;
 maxD=100;
 tic
@@ -264,41 +263,18 @@ while(maxD>Precision)
     %% 形成等式约束的雅克比
     deltH=func_deltH(Busnum,Volt,Y,UAngel,r,c,Yangle,Loadi);
     %% 形成不等式约束的雅克比
-%     deltG=func_deltG(Busnum,Loadi);
-    %%
-    %% Ðγɺ£É­Õó
-    deltdeltF=func_deltdeltF(wVolt,wPD,wQD,ContrlCount,Loadi);
-%     deltdeltF=zeros(size(deltdeltF));
-    %% ÐγÉddHy
-    ddh=func_ddh(Volt,Init_Y,Busnum,Y,UAngel,r,c,Yangle,Loadi,ContrlCount);
-%     ddh=zeros(size(ddh));
     %%  开始构建ddg
     %%  开始构建deltF
     deltF=func_deltF(wVolt,wPD,wQD,mPD3P,PD3P,QD3P,mQD3P,Volt,mVoltABCV,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,PD3P,QD3P,Loadi);
     Mat_H=FormH(Busnum,Volt,PG3P,PD3P,QG3P,QD3P,Y,UAngel,r,c,Yangle,Loadi);
-    Ly=Mat_H;
-%     if isSetBound==0
-%         PDL=0.70*rPD3P.*(1+unifrnd(-0.15,0.15,length(rPD3P),1));
-%         QDL=0.70*rQD3P.*(1+unifrnd(-0.15,0.15,length(rQD3P),1));
-%         PDU=1.3*rPD3P.*(1+unifrnd(-0.15,0.15,length(rPD3P),1));
-%         QDU=1.3*rQD3P.*(1+unifrnd(-0.15,0.15,length(rPD3P),1));
-%         isSetBound=1;
-%     end
-%     Lz=FormLz(Mat_G,Init_L,Loadi,rPD3P,rQD3P,rVoltABCV,PDL,QDL);
-%     Lw=FormLw(Mat_G,Init_U,Loadi,rPD3P,rQD3P,rVoltABCV,PDU,QDU);
-    Lx=FormLx(deltF,deltH,Init_Y);
-    %     YY=FormYY(Lul,Lz,Ly,Luu,Lw,Lx);
     %% 开始解方程
     %     fprintf('迭代次数 %d Gap %f\n',KK+1,plotGap(KK+1));
-    XX=SolveIt(deltdeltF,ddh,deltH,Init_Y,Ly,ContrlCount,Lx,Balance,Busnum,Loadi);
+    XX=SolveIt(deltH,ContrlCount,Balance,Busnum,Loadi,deltF,mPD3P,PD3P,mQD3P,QD3P,mVoltABCV,Volt,Mat_H,wVolt,wPD,wQD);
     %%取各分量
     [deltX,deltY]=AssignXX(XX,ContrlCount,Loadi,Balance,Busnum);
     [Init_Y,PG,QG,Volt,UAngel,PD3P,QD3P]=Modification(Init_Y,deltX,deltY,PG,QG,Volt,UAngel,ContrlCount,Balance,Busnum,PD3P,QD3P,Loadi);
-    maxD=max(abs([deltX;deltY]));
+    maxD=max(abs([deltX]));
     fprintf('%f\n',full(maxD));
     KK=KK+1;
 end

README.md

@@ -1 +1,2 @@
+2015.5.19 在之前的程序上改成用Gauss-Newton法，而不是含有海森矩阵的Newton法
 把原来内点法的程序给改了，改成只有等式约束的普通WLS方法。

SolveIt.m

@@ -1,12 +1,28 @@
-function XX=SolveIt(deltdeltF,ddh,deltH,Init_Y,Ly,ContrlCount,Lx,Balance,Busnum,Loadi)
-LxComa=FormLxComa(Lx);%Lx=deltF-deltH*Init_Y';
-Lx=LxComa;
-H=-deltdeltF+ddh;
+function XX=SolveIt(deltH,ContrlCount,Balance,Busnum,Loadi,deltF,PD0,PD,QD0,QD,mVolt,Volt,Mat_H,wVolt,wPD,wQD)
+% LxComa=FormLxComa(Lx);%Lx=deltF-deltH*Init_Y';
+% Lx=LxComa;
+% H=-deltdeltF+ddh;
+% aa=[
+%     H,deltH;
+%     deltH',zeros(size(Init_Y,2));
+%     ];
 aa=[
-    H,deltH;
-    deltH',zeros(size(Init_Y,2));
+    deltF'*diag([wPD ;wQD ;wVolt])*deltF deltH'
+    deltH ones(size(deltH,1),size(deltH',2))
+    
+];
+% yy=[Lx;-Ly];
+% t3=2*wPD.*(PD-PD0);
+% t4=2*wQD.*(QD-QD0);
+% t5=2*wVolt.*(Volt-mVolt);
+deltZ=[(PD0-PD);
+    (QD0-QD);
+    (mVolt-Volt);
     ];
-yy=[Lx;-Ly];
+yy=[
+    deltF'*diag([wPD; wQD; wVolt])*deltZ;
+-Mat_H    ;
+];
 %% 平衡节点电压不变 
 t=size(Loadi,1)*2;
 aa(t+(Balance-1)*3+1,:)=0;

func_deltF.m

@@ -1,11 +1,18 @@
 function deltF=func_deltF(wVolt,wPD,wQD,PD0,PD,QD,QD0,Volt,mVolt,Busnum,Loadi)
-t3=2*wPD.*(PD-PD0);
-t4=2*wQD.*(QD-QD0);
-t5=2*wVolt.*(Volt-mVolt);
-deltF=[sparse(t3);
-    sparse(t4);
-    sparse(t5);
-    sparse(1*Busnum*3,1);
-];
-
+% t3=2*wPD.*(PD-PD0);
+% t4=2*wQD.*(QD-QD0);
+% t5=2*wVolt.*(Volt-mVolt);
+% deltF=[sparse(t3);
+%     sparse(t4);
+%     sparse(t5);
+%     sparse(1*Busnum*3,1);
+% ];
+totalM=length(PD)+length(QD)+length(Volt);
+deltF=[
+    eye(length(PD)),sparse(length(PD),totalM-length(PD));
+    sparse(length(QD),length(PD)),eye(length(PD)),sparse(length(QD),totalM-length(PD)-length(QD));
+    sparse(length(Volt),length(PD)+length(QD)),eye(length(Volt));
+    sparse(length(Volt),totalM);
+]';
+deltF=sparse(deltF);
 end

func_deltH.m

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