修复一个bug。潮流形成导纳矩阵的时候忘记加节点电导（G）了，所以老是不对。现在修复了。

Signed-off-by: facat <dugg@21cn.com>

Initial.m

@@ -8,6 +8,6 @@ function [P0,Q0,U,Uangle] = Initial(PG,PD,PQstandard,Pointpoweri,QG,QD,Busnum)
 P0 = sparse(1, Pointpoweri,(PG-PD)/PQstandard);                         % 求取节点注入有功功率的标幺值
 Q0 = sparse(1, Pointpoweri,(QG-QD)/PQstandard);                       % 求取节点注入无功功率的标幺值
 %% 平启动赋电压初值
-U = ones(1,Busnum);                                  % 按照平启动给电压幅值赋值
+U = 1*ones(1,Busnum);                                  % 按照平启动给电压幅值赋值
 Uangle = zeros(1,Busnum);                       % 按照平启动给电压相角赋值
 end

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('ieee30.dat');
+pf('c:/newFIle.txt');
 %pf('D:\Project\青秀降损项目\最小化潮流\最小潮流算例\原始\津头站津视922（3-1）_0.5_120%.txt'); 
 %pf('D:\Project\最小化潮流\最小潮流算例\仙海919.txt'); 
 %pf('c:/file31.txt'); 
@@ -96,14 +96,15 @@ 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(absPDLoad<10))
 absQDLoad=abs(  (QD(Loadi)-QDReal(Loadi))./QDReal(Loadi)  );
-maxQDError=max(absQDLoad)
+maxQDError=max(absQDLoad(absQDLoad<10))
 disp('index');
-Loadi(absPDLoad==maxPDError);
+%Loadi(absPDLoad==maxPDError);
 %% 计算总线损
 totalLoss=(sum(PG)-sum(PD(Loadi)))*100;
 fprintf('总的损耗为%f(MW 有名值)\n',full(totalLoss));
+fprintf('线损率为 %f\n',full(totalLoss/sum(PG)));
 %% 计算各线损
 %Lineloss(Linei,Linej,Liner,Linex,Lineb,Transfori,Transforj,Transforr,Transforx,Transfork0,Volt,UAngel);
 toc  

admmatrix.m

@@ -1,5 +1,5 @@
 function [GB,Y,r,c,Angle]  = admmatrix(Busnum,Linei,Linej,Liner,Linex,Lineb,Transfori...
-    ,Transforj,Transforr,Transforx,Transfork0,Branchi,Branchb)
+    ,Transforj,Transforr,Transforx,Transfork0,Branchi,Branchg,Branchb)
 %**************************************************************************
 %  程序功能 : 子函数——形成节点导纳矩阵Y
 %  编    者： 
@@ -22,6 +22,7 @@ end
 %% 接地支路计算
 if Branchi>0                                                                 % 判断有无接地支路
     B = B+sparse(Branchi,Branchi,Branchb,Busnum,Busnum);
+    G = G+sparse(Branchi,Branchi,Branchg,Busnum,Busnum);
 end
 %% 化作极坐标形式
 GB = G+B.*1i;                                                              %将电导，电纳合并，写成复数形式

jacobian.asv

@@ -0,0 +1,40 @@
+function [Jacob,PQ,U,Uangle]=jacobian(Busnum,Balance,PVi,PVu,U,Uangle,Y,Angle,P0,Q0,r,c)   
+%**************************************************************************
+%  程序功能 : 子函数——形成雅可比矩阵Jacobian
+%  编    者： 
+%  编制时间：2010.12
+%**************************************************************************
+%% 分别求雅克比矩阵的子阵H,L,N,J及有功无功分量P,Q
+AngleIJ = Uangle(r) - Uangle(c)- Angle';
+U(PVi) = PVu;
+U(Balance)=;
+temp1= -sparse(1:Busnum,1:Busnum,U,Busnum,Busnum)*Y*sparse(1:Busnum,1:Busnum,U,Busnum,Busnum);                                                          % 计算雅克比矩阵可利用的中间变量
+temp2 = sum(temp1.*sparse(r,c,sin(AngleIJ)),2);                     
+temp3 = sum(temp1.*sparse(r,c,cos(AngleIJ)),2);
+temp4=sparse(1:Busnum,1:Busnum,temp2,Busnum,Busnum);
+temp5=sparse(1:Busnum,1:Busnum,temp3,Busnum,Busnum);
+H = temp1.*sparse(r,c,sin(AngleIJ))-temp4;
+L = temp1.*sparse(r,c,sin(AngleIJ))+temp4;
+N = temp1.*sparse(r,c,cos(AngleIJ))+temp5;
+J = -temp1.*sparse(r,c,cos(AngleIJ))+temp5;
+
+Q = Q0+temp2';                                                                              %求有功分量P
+P = P0+temp3';                                                                               %求无功分量Q
+%% 处理平衡节点和pv节点
+H(:,Balance) = 0;                                                     
+H(Balance,:) = 0;
+%H(Balance,Balance) = 100;                                                           % 平衡节点对应的对角元素置一个有限数
+H=H+sparse(Balance,Balance,ones(1,length(Balance)),Busnum,Busnum);
+L(:,PVi) = 0;
+L(PVi,:) = 0;
+L = L+sparse(PVi,PVi,ones(1,length(PVi)),Busnum,Busnum); % PV节点对应的对角元素置为1
+J(:,Balance) = 0;
+J(PVi,:) = 0;
+N(:,PVi) = 0;
+N(Balance,:) = 0; 
+Q(PVi) = 0;                                                                                       % 将pv节点的无功不平衡分量置零   
+P(Balance) = 0;                                                                               % 平衡节点的有功功率不平衡分量置零
+%% 合成PQ和雅可比矩阵
+PQ = cat(2,P,Q);                                                                             % 形成功率不平衡分量列向量
+Jacob = cat(1,cat(2,H,N),cat(2,J,L));                                            % 形成Jacobian矩阵
+end

jacobian.m

@@ -7,6 +7,7 @@ function [Jacob,PQ,U,Uangle]=jacobian(Busnum,Balance,PVi,PVu,U,Uangle,Y,Angle,P0
 %% 分别求雅克比矩阵的子阵H,L,N,J及有功无功分量P,Q
 AngleIJ = Uangle(r) - Uangle(c)- Angle';
 U(PVi) = PVu;
+U(Balance)=1;
 temp1= -sparse(1:Busnum,1:Busnum,U,Busnum,Busnum)*Y*sparse(1:Busnum,1:Busnum,U,Busnum,Busnum);                                                          % 计算雅克比矩阵可利用的中间变量
 temp2 = sum(temp1.*sparse(r,c,sin(AngleIJ)),2);                     
 temp3 = sum(temp1.*sparse(r,c,cos(AngleIJ)),2);

openfile.m

@@ -51,10 +51,10 @@ QG=QG/Base;
 PD=PD/Base;
 QD=QD/Base;
 %%
-PD=sparse(PD);
+PD=sparse(PD)/Base;
-QD=sparse(QD);
+QD=sparse(QD)/Base;
-PG=sparse(PG);
+PG=sparse(PG)/Base;
-QG=sparse(QG);
+QG=sparse(QG)/Base;
 %% pv节点功率参数矩阵
 PVi = PV(:,1); % PV节点的节点号
 PVu = PV(:,2); % PV节点电压

openfile2.m

@@ -1,5 +1,5 @@
 function [Busnum,Balance,PQstandard,Precision,Linei,Linej,Liner,Linex,Lineb,kmax,Transfori ,...
-    Transforj,Transforr,Transforx,Transfork0,Branchi,Branchb,Pointpoweri,PG,QG,PD,QD,PVi,PVu,GenU,GenL,GenC,CenterA,PGi,PVQU,PVQL] = openfile2(FileName)
+    Transforj,Transforr,Transforx,Transfork0,Branchi,Branchg,Branchb,Pointpoweri,PG,QG,PD,QD,PVi,PVu,GenU,GenL,GenC,CenterA,PGi,PVQU,PVQL] = openfile2(FileName)
 %**************************************************************************
 %    程序简介 : 子函数——读取潮流计算所需数据
 %    编    者：
@@ -34,6 +34,7 @@ Lineb = line(:,6);                                                            %
 %% 对地支路参数矩阵
 Branchi = ground(:,2);                                                  % 对地支路节点号
 Branchb = ground(:,4);                                                 % 对地支路的导纳
+Branchg = ground(:,3);                                                 % 对地支路的导纳
 %% 变压器参数矩阵
 Transfori = tran(:,3);                                                     % 节点i
 Transforj= tran(:,4);                                                      % 节点j
@@ -54,6 +55,7 @@ QD=QD/Base;
 %%
 PD=sparse(PD);
 QD=sparse(QD);
+%QD=PD*sqrt(1-.85^2)/.85;
 PG=sparse(PG);
 QG=sparse(QG);
 %% pv节点功率参数矩阵

pf.m

@@ -10,10 +10,10 @@ function [kmax,Precision,Uangle,U,Busnum,PVi,PVu,Balance,Y,Angle,P0,Q0,r,c,GB,Li
 tic;
 %% 读取数据文件
 [Busnum,Balance,PQstandard,Precision,Linei,Linej,Liner,Linex,Lineb,kmax,Transfori ,...
-    Transforj,Transforr,Transforx,Transfork0,Branchi,Branchb,Pointpoweri,PG,QG,PD,QD,PVi,PVu,GenU,GenL,GenC,CenterA,PGi,PVQU,PVQL]= openfile(FileName);
+    Transforj,Transforr,Transforx,Transfork0,Branchi,Branchg,Branchb,Pointpoweri,PG,QG,PD,QD,PVi,PVu,GenU,GenL,GenC,CenterA,PGi,PVQU,PVQL]= openfile2(FileName);
 %% 形成节点导纳矩阵
 [GB,Y,r,c,Angle] = admmatrix(Busnum,Linei,Linej,Liner,Linex,Lineb,Transfori,Transforj,Transforr,...
-    Transforx,Transfork0,Branchi,Branchb);
+    Transforx,Transfork0,Branchi,Branchg,Branchb);
 [P0,Q0,U,Uangle] = Initial(PG,PD,PQstandard,Pointpoweri,QG,QD,Busnum);                         %求功率不平衡量
 %disp('迭代次数i  最大不平衡量');
 %% 循环体计算
@@ -26,13 +26,9 @@ for i = 0:kmax
     if m > Precision                                                                                                                                 %判断不平衡量是否满足精度要求
         [Uangle,U] = solvefun(Busnum,Jacob,PQ,Uangle,U);                                                           %求解修正方程，更新电压变量
     else
-        %disp(['收敛，迭代次数为',num2str(i),'次']);
+        disp(['收敛，迭代次数为',num2str(i),'次']);
         break                                                                                                                                             %若满足精度要求，则计算收敛
     end
 end
 toc;
-PG=PG/PQstandard;
-QG=QG/PQstandard;
-PD=PD/PQstandard;
-QD=QD/PQstandard;
 end