为了论文用的数据，加入了读变压器容量，用变压器容量做无功上限。

FormLw.m

@@ -1,4 +1,4 @@
-function Lw=FormLw(Mat_G,Init_U,GenU,Busnum,PVQU,PD,PD0,QD0,Loadi,KK,PF)
+function Lw=FormLw(Mat_G,Init_U,GenU,Busnum,PVQU,PD,PD0,QD0,Loadi,KK,PF,noDataTransCapacity)
 KK=999;
 %PU=GenU(:,2);%发电机有功上界
 PU=5*ones(length(GenU(:,2)),1);
@@ -6,7 +6,8 @@ PU=5*ones(length(GenU(:,2)),1);
 QU=5*ones(length(PVQU(:,1)),1);
 VoltU=(1.1)*ones(1,Busnum);
 %VoltU=10*ones(1,Busnum);
-PDU=PD0(Loadi);
+% PDU=PD0(Loadi);
+PDU=noDataTransCapacity;
 PDU(PDU>0)=1.200*PDU(PDU>0);
 PDU(PDU<0)=0.800*PDU(PDU<0);
 PDU(PDU==0)=0.400;
@@ -15,7 +16,8 @@ PDU(PDU==0)=0.400;
 % QDU(QDU>0)=1.200*QDU(QDU>0);
 % QDU(QDU<0)=0.800*QDU(QDU<0);
 % QDU(QDU==0)=0.200;
-QDU=1.1*PD(Loadi).*sqrt(1 -PF.^2)./PF;
+PF=0.85;
+QDU=1.0*PD(Loadi).*sqrt(1 -PF.^2)./PF;
 t1=([PU',QU',PDU',QDU',VoltU])';
 t2=Mat_G+Init_U'-t1;
 Lw=t2;

FormLz.m

@@ -15,7 +15,7 @@ PDL(PDL==0)=-0.400;
 % QDL(QDL>0)=0.800*QDL(QDL>0);
 % QDL(QDL<0)=1.200*QDL(QDL<0);
 % QDL(QDL==0)=-0.200;
-QDL=0.9*PD(Loadi).*sqrt((1-PF.^2))./PF;
+QDL=0*PD(Loadi).*sqrt((1-PF.^2))./PF;
 t1=([PL',QL',PDL',QDL',VoltL])';
 t2=Mat_G-Init_L'-t1;
 Lz=t2;

OPF.asv

@@ -2,41 +2,56 @@ 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('ieee4.dat');
+pf('c:/newFIle.txt');
 %pf('D:\Project\青秀降损项目\最小化潮流\最小潮流算例\原始\津头站津视922（3-1）_0.5_120%.txt'); 
 %pf('D:\Project\最小化潮流\最小潮流算例\仙海919.txt'); 
 %pf('c:/file31.txt'); 
 
 
 %% 计算功率因数
-%atan(PD(QD~=0 | PD~=0)./QD(QD~=0 | PD~=0));
+Loadi=QD~=0 | PD~=0;
+PF=sqrt(PD(Loadi).^2./(QD(Loadi).^2+PD(Loadi).^2))
+%% 
 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;
 PD0=PD;
 QD0=QD;
 PDReal=PD;%真值
+QDReal=QD;%真值
 %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,50);
+plotGap=zeros(1,60);
 ContrlCount=size(PVi,1)+size(PGi,1)+size(Loadi,1)*2+Busnum*2;
 kmax=60;
 %% 20120523 临时
 QD_NON_ZERO=QD(PD==0 & QD~=0);
 QD_NON_ZERO_IND=find(PD==0 & QD~=0);
 %%
+Precision=Precision/1;
+
+%% 加误差
+PD0(Loadi)=PD0(Loadi).*(1+normrnd(0,0.05,length(Loadi),1));
+QD0(Loadi)=QD0(Loadi).*(1+normrnd(0,0.05,length(Loadi),1));
+PG0(PGi)=PG0(PGi).*(1+normrnd(0,0.01,length(PGi),1));
+QG0(PVi)=QG0(PVi).*(1+normrnd(0,0.01,length(PVi),1));
+%% 读变压器容量
+[noDataTransNum noDataTransCapacity noDataTransPowerFactor]=ReadNoDataTrans(fileName);
 while(abs(Gap)>Precision)
     if KK>kmax
         break;
@@ -48,31 +63,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,wPG,wQG,wD,ContrlCount);
+    deltdeltF=func_deltdeltF(PVi,wPG,wQG,wPD,wQD,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,QG,PVi,PGi,wPG,wQG,wD,PG0,QG0,PD0,PD,Busnum,Loadi);
+    deltF=func_deltF(PG,QG,PVi,PGi,wPG,wQG,wPD,wQD,PG0,QG0,PD0,PD,QD,QD0,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);
+    Lz=FormLz(Mat_G,Init_L,GenL,Busnum,PVQL,PD,PD0,QD0,Loadi,KK,PF);
-    Lw=FormLw(Mat_G,Init_U,GenU,Busnum,PVQU,PD0,QD0,Loadi,KK);
+    Lw=FormLw(Mat_G,Init_U,GenU,Busnum,PVQU,PD,PD0,QD0,Loadi,KK,PF);
     Lx=FormLx(deltF,deltH,Init_Y,deltG,Init_Z,Init_W);
     YY=FormYY(Lul,Lz,Ly,Luu,Lw,Lx);
     %% 开始解方程
-    fprintf(')
+    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);
@@ -81,7 +96,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 +104,16 @@ 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(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);
+%Lineloss(Linei,Linej,Liner,Linex,Lineb,Transfori,Transforj,Transforr,Transforx,Transfork0,Volt,UAngel);
 toc  
 

OPF.m

@@ -44,6 +44,14 @@ QD_NON_ZERO=QD(PD==0 & QD~=0);
 QD_NON_ZERO_IND=find(PD==0 & QD~=0);
 %%
 Precision=Precision/1;
+
+%% 加误差
+PD0(Loadi)=PD0(Loadi).*(1+normrnd(0,0.05,length(Loadi),1));
+QD0(Loadi)=QD0(Loadi).*(1+normrnd(0,0.05,length(Loadi),1));
+PG0(PGi)=PG0(PGi).*(1+normrnd(0,0.01,length(PGi),1));
+QG0(PVi)=QG0(PVi).*(1+normrnd(0,0.01,length(PVi),1));
+%% 读变压器容量
+[noDataTransNum noDataTransCapacity noDataTransPowerFactor]=ReadNoDataTrans('C:/b/东际911_2751267_2012-09-05/iPso_东际911_2751267_2012-09-05_变压器无负载.txt');
 while(abs(Gap)>Precision)
     if KK>kmax
         break;
@@ -75,7 +83,7 @@ while(abs(Gap)>Precision)
     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,PD,PD0,QD0,Loadi,KK,PF);
-    Lw=FormLw(Mat_G,Init_U,GenU,Busnum,PVQU,PD,PD0,QD0,Loadi,KK,PF);
+    Lw=FormLw(Mat_G,Init_U,GenU,Busnum,PVQU,PD,PD0,QD0,Loadi,KK,PF,noDataTransCapacity);
     Lx=FormLx(deltF,deltH,Init_Y,deltG,Init_Z,Init_W);
     YY=FormYY(Lul,Lz,Ly,Luu,Lw,Lx);
     %% ¿ªÊ¼½â·½³Ì

OPF_Init.m

@@ -22,8 +22,8 @@ wQG=100*ones(size(PVi,1),1);
 %randInt=randperm(size(Loadi,1));
 %randPDind=randInt(1:10);
 randPDind=0;
-wPD=1*ones(size(Loadi,1),1);
+wPD=20*ones(size(Loadi,1),1);
-wQD=1*ones(size(Loadi,1),1);
+wQD=20*ones(size(Loadi,1),1);
 %wD(randPDind)=0;%一些负荷不约束
 %wD(7)=0;
 % wD(11)=0;

ReadNoDataTrans.m

@@ -0,0 +1,9 @@
+function [noDataTransNum noDataTransCapacity noDataTransPowerFactor]=ReadNoDataTrans(fileName)
+data = dlmread(fileName); 
+tmpMat=sortrows(data,1);
+noDataTransNum=tmpMat(:,1);
+%noDataTransCapacity=tmpMat(:,2)/100000;
+noDataTransCapacity=tmpMat(:,2)/1000;
+noDataTransPowerFactor=tmpMat(:,3);
+
+end