将电压加入目标函数中。

@ForThesis/MaxBranchDeviation.m

@@ -29,8 +29,8 @@ end
 Volt0=U;
 UAngel0=Uangle;
 
-[dispLineloss0 dispTransloss0]=Lineloss(Linei,Linej,Liner,Linex,Lineb,Transfori,Transforj,Transforr,Transforx,Branchi,Branchg,Branchb,Transfork0,Volt0,UAngel0);%测量值
+[dispLineloss0 dispTransloss0]=Lineloss(Linei,Linej,Liner,Linex,Lineb,Transfori,Transforj,Transforr,Transforx,Branchi,Branchg,Branchb,Transfork0,Volt0,UAngel0);
-[dispLineloss dispTransloss]=Lineloss(Linei,Linej,Liner,Linex,Lineb,Transfori,Transforj,Transforr,Transforx,Branchi,Branchg,Branchb,Transfork0,Volt,UAngel);%估计值
+[dispLineloss dispTransloss]=Lineloss(Linei,Linej,Liner,Linex,Lineb,Transfori,Transforj,Transforr,Transforx,Branchi,Branchg,Branchb,Transfork0,Volt,UAngel);
 t1=(dispLineloss0(:,3) - dispLineloss(:,3))./dispLineloss0(:,3);
 t2=(dispTransloss0(:,3) - dispTransloss(:,3))./dispTransloss0(:,3);
 t11=dispLineloss0(:,3)>1e-5;% 太小的值不计算

FormG.m

@@ -1,10 +1,8 @@
-function Mat_G=FormG(Volt,PVi,PGi,PG,QG,PD,QD,Loadi,Vbi)
+function Mat_G=FormG(Volt,PVi,PGi,PG,QG,PD,QD,Loadi)
 
 Mat_G=[
     sparse(PD(Loadi));
     sparse(QD(Loadi));
-    Volt'-Vbi;
+    Volt';
-    Volt'+Vbi;
-    Vbi;
     ];
 end

FormH.m

@@ -1,4 +1,10 @@
-function Mat_H=FormH(Busnum,Volt,PG,PD,QG,QD,Y,UAngel,r,c,Angle,Loadi)
+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.asv

@@ -0,0 +1,25 @@
+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);
+%QU=PVQU(:,1);%发电机无功上界
+QU=5*ones(length(PVQU(:,1)),1);
+VoltU=(1.1)*ones(1,Busnum);
+%VoltU=10*ones(1,Busnum);
+% 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;
+%PDU=10*ones(length(Loadi),1);
+% QDU=QD0(Loadi);
+% QDU(QDU>0)=1.200*QDU(QDU>0);
+% QDU(QDU<0)=0.800*QDU(QDU<0);
+% QDU(QDU==0)=0.200;
+% 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;
+
+end

FormLw.m

@@ -9,22 +9,14 @@ PDU=PD0(Loadi);
 PDU(PDU>0)=1.200*PDU(PDU>0);
 PDU(PDU<0)=0.800*PDU(PDU<0);
 PDU(PDU==0)=0.400;
-realPD=PD0(Loadi);
-indPD=find(realPD>0);
-PDU(indPD(3:12:end))=1.55*realPD(indPD(3:12:end));
-PDU(indPD(9:12:end))=1.05*realPD(indPD(9:12:end));
 %PDU=10*ones(length(Loadi),1);
 QDU=QD0(Loadi);
 QDU(QDU>0)=1.200*QDU(QDU>0);
 QDU(QDU<0)=0.800*QDU(QDU<0);
 QDU(QDU==0)=0.400;
-realQD=QD0(Loadi);
-indQD=find(realQD>0);
-QDU(indQD(3:12:end))=1.55*realQD(indQD(3:12:end));
-QDU(indQD(9:12:end))=1.05*realQD(indQD(9:12:end));
 % PF=0.85;
 % QDU=1.0*PD(Loadi).*sqrt(1 -PF.^2)./PF;
-t1=([PDU',QDU',1.1*VoltU,1.1*VoltU,1*ones(1,Busnum)])';
+t1=([PDU',QDU',VoltU])';
 t2=Mat_G+Init_U'-t1;
 Lw=t2;
 

FormLz.asv

@@ -0,0 +1,20 @@
+function Lz=FormLz(Mat_G,Init_L,GenL,Busnum,PVQL,PD,PD0,QD0,Loadi,KK,PF)
+KK=999;
+
+VoltL=(0.9)*ones(1,Busnum);
+%VoltL=-10*ones(1,Busnum);
+PDL=PD0(Loadi);
+PDL(PDL>0)=0.700*PDL(PDL>0);
+PDL(PDL<0)=1.300*PDL(PDL<0);
+PDL(PDL==0)=-0.400;
+%PDL=-10*ones(length(Loadi),1);
+QDL=QD0(Loadi);
+QDL(QDL>0)=0.700*QDL(QDL>0);
+QDL(QDL<0)=1.300*QDL(QDL<0);
+QDL(QDL==0)=-0.400;
+% QDL=0*PD(Loadi).*sqrt((1-PF.^2))./PF;
+t1=([PDL',QDL',VoltL])';
+t2=Mat_G-Init_L'-t1;
+Lz=t2;
+
+end

FormLz.m

@@ -7,22 +7,13 @@ PDL=PD0(Loadi);
 PDL(PDL>0)=0.800*PDL(PDL>0);
 PDL(PDL<0)=1.200*PDL(PDL<0);
 PDL(PDL==0)=-0.400;
-realPD=PD0(Loadi);
-indPD=find(realPD>0);
-PDL(indPD(3:12:end))=0.95*realPD(indPD(3:12:end));
-PDL(indPD(9:12:end))=0.45*realPD(indPD(9:12:end));
 %PDL=-10*ones(length(Loadi),1);
 QDL=QD0(Loadi);
 QDL(QDL>0)=0.800*QDL(QDL>0);
 QDL(QDL<0)=1.200*QDL(QDL<0);
 QDL(QDL==0)=-0.400;
-realQD=QD0(Loadi);
-indQD=find(realQD>0);
-QDL(indQD(3:12:end))=0.95*realQD(indQD(3:12:end));
-QDL(indQD(9:12:end))=0.95*realQD(indQD(9:12:end));
 % QDL=0*PD(Loadi).*sqrt((1-PF.^2))./PF;
-%t1=([PDL',QDL',VoltL,VoltL,0*ones(1,Busnum)])';
+t1=([PDL',QDL',VoltL])';
-t1=([PDL',QDL',-20000*VoltL,-20000*VoltL,0*ones(1,Busnum)])';
 t2=Mat_G-Init_L'-t1;
 Lz=t2;
 

JSMJZM.asv

@@ -0,0 +1,12 @@
+JSM=[0.0575 0.051 0.0029];
+JZM=[0.003 0.02629 0.0582];
+plot(1:3,JSM,'k');
+axis([0 4 0 0.07])
+text(1,JSM(1)+0.001,'0.0575')
+text(2,JSM(2)+0.001,'0.051')
+text(3,JSM(3)+0.001,'0.0029')
+hold on 
+plot(1:3,JZM,'k--');
+text(1,JZM(1)+0.001,'0.003')
+text(2,JZM(2)+0.001,'0.02629')
+text(3,JZM(3)+0.001,'0.0029')

JSMJZM.m

@@ -1,18 +1,18 @@
-JSM=[0.070244 0.0533 0.004379];
+JSM=[0.0575 0.051 0.0029];
-JZM=[0.002275 0.047311 0.070318];
+JZM=[0.0003 0.02629 0.0582];
-for I=0.01:0.01:0.08
+for I=0.01:0.01:0.06
     line([0.5 3.5],[I I],'Color',[220 220 220]/255)
 end
 hold on
 plot(1:3,JSM,'k');
 ylabel('ͳ¼ÆÎó²î')
-xlabel('ÇéÐÎ')
+xlabel('Çé¿ö')
-axis([0.5 3.5 0 0.075])
+axis([0.5 3.5 0 0.065])
-text(1,JSM(1)+0.001,'0.0702')
+text(1,JSM(1)+0.001,'0.0575')
-text(2,JSM(2)+0.001,'0.0533')
+text(2,JSM(2)+0.001,'0.0510')
-text(3,JSM(3)+0.001,'0.0043')
+text(3,JSM(3)+0.001,'0.0029')
 hold on 
 plot(1:3,JZM,'k--');
-text(1-.3,JZM(1)+0.001,'0.0023')
+text(1-.3,JZM(1)+0.001,'0.0003')
-text(2-0.3,JZM(2)-0.001,'0.0473')
+text(2,JZM(2)-0.001,'0.0263')
-text(3,JZM(3)-0.001,'0.0703')
+text(3,JZM(3)-0.001,'0.0029')

JSNZN.asv

@@ -0,0 +1,14 @@
+JSN=[0.1213 0.1092 0.0587];
+JZN=[0.0971 0.1314 0.1408];
+plot(1:3,JSM,'k');
+ylabel('Öľ')
+xlabel('Çéżö')
+axis([0.5 3.5 0 0.065])
+text(1,JSM(1)+0.001,'0.1213')
+text(2,JSM(2)+0.001,'0.0510')
+text(3,JSM(3)+0.001,'0.0029')
+hold on 
+plot(1:3,JZM,'k--');
+text(1-.3,JZM(1)+0.001,'0.0003')
+text(2,JZM(2)-0.001,'0.0263')
+text(3,JZM(3)-0.001,'0.0029')

Lineloss.asv

@@ -0,0 +1,38 @@
+%% 计算线损
+function [dispLineloss dispTransloss]=Lineloss(Linei,Linej,Liner,Linex,Lineb2,Transi,Transj,Transr,Transx,Branchi,Branchg,Branchb,k0,Volt,Angle)
+%format long
+% fprintf('功率为有名值\n');
+% fprintf('节点号\t节点号\t有功损耗 MW\t无功损耗 MVar')
+cmpVolt=Volt'.*cos(Angle')+1i*Volt'.*sin(Angle');
+% cmpVolt=cmpVolt';
+y0=1i*Lineb2;
+%yj0=1i*standardinput(:,7);
+yij=1./(Liner+1i*Linex);
+%% 线路损耗
+Sij=cmpVolt(Linei).*conj(  cmpVolt(Linei) .* y0 + ( cmpVolt(Linei)- cmpVolt(Linej) ).*yij );
+Sji=cmpVolt(Linej).*conj(  cmpVolt(Linej) .*y0 + ( cmpVolt(Linej)- cmpVolt(Linei) ).*yij );
+%Sij1==Sij2
+deltLineS=Sij+Sji;
+
+%% 另一种计算方式begin
+ss=1*(Volt(Linei)'.^2.*abs(yij).*cos( angle(yij) ) -Volt(Linei)'.*Volt(Linej)'.*cos( Angle(Linei)' - Angle(Linej)' - angle(yij)).*abs(yij));
+ss=(Volt(Linei)'.^2+Volt(Linej)'.^2).*abs(yij).*cos(angle(yij))-2*Volt(Linei)'.*Volt(Linej)'.*cos( Angle(Linei)' - Angle(Linej)').*cos( - angle(yij)).*abs(yij);
+ss=Volt(Linei)'.*Volt(Linej)'.*abs(yij).*cos()
+%% 另一种计算方式end
+dispLineloss=[Linei Linej real(deltLineS)*100 imag(deltLineS)*100];
+%full(dispLineloss)
+% dispLineloss=sortrows(dispLineloss,-3);
+full(dispLineloss);
+%% 以下是变压器损耗
+yij=1./(Transr+1i*Transx);
+Sij=cmpVolt(Transi)./k0.*conj(  ( cmpVolt(Transi)./k0- cmpVolt(Transj) ).*yij );
+Sji=cmpVolt(Transj).*conj( ( cmpVolt(Transj)- cmpVolt(Transi)./k0 ).*yij );
+deltTransS=Sij+Sji;
+%% 接地支路损耗
+% 没有考虑变压器变比
+deltTransS =deltTransS+sum(cmpVolt(Branchi).*conj((cmpVolt(Branchi).*(Branchg+1j*Branchb))));
+%% 
+dispTransloss=[Transi Transj real(deltTransS)*100 imag(deltTransS)*100];
+% dispTransloss=sortrows(dispTransloss,-3);
+full(dispTransloss);
+end

Modification.m

@@ -1,4 +1,4 @@
-function [Init_Z,Init_L,Init_W,Init_U,Init_Y,PG,QG,Volt,UAngel,PD,QD,Vbi]=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,Vbi)
+function [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)
 AlphaP=FormAlphaP(Init_L,deltL,Init_U,deltU);
 %fprintf('AlphaP %f\n',full(AlphaP));
 AlphaD=FormAlphaD(Init_Z,deltZ,Init_W,deltW);
@@ -24,5 +24,4 @@ balVolt=Volt(Balance);
 Volt=Volt+AlphaP*t(1:Busnum);
 Volt(Balance)=balVolt;
 UAngel=UAngel+AlphaP*t(Busnum+1:2*Busnum);
-Vbi=Vbi+AlphaP*t(2*Busnum+1:end)';
 end

OPF.asv

@@ -1,80 +1,69 @@
 tic
 clc
 clear
-%% 存在问题
+thesis=ForThesis(4,8);
-% 变压器变比的位置没有考虑，由于现在用的变比都是1，所以没有影响。 20130123
+for II=1:4
-%%
+    [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]= ...
-thesis=ForThesis(1,62);
+        pf('c:/newFIle.txt');
-[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,Branchi,Branchg,Branchb,Transfork0]= ...
+    %pf('D:\Project\青秀降损项目\最小化潮流\最小潮流算例\原始\津头站津视922（3-1）_0.5_120%.txt');
-    pf('E:\算例\东际911_2751267_2012-09-05\newFIle20-使用22.txt');
+    %pf('D:\Project\最小化潮流\最小潮流算例\仙海919.txt');
-%         pf('E:\算例\柳金Ⅰ926_21671693_2012-09-06\newFIle16.txt');
+    %pf('c:/file31.txt');
-%pf('D:\Project\青秀降损项目\最小化潮流\最小潮流算例\原始\津头站津视922（3-1）_0.5_120%.txt');
+    
-%pf('D:\Project\最小化潮流\最小潮流算例\仙海919.txt');
+    
-%pf('c:/file31.txt');
+    %% 计算功率因数
-%% 计算功率因数
+    Loadi=QD~=0 | PD~=0;
-Loadi=QD~=0 | PD~=0;
+    PF=sqrt(PD(Loadi).^2./(QD(Loadi).^2+PD(Loadi).^2))
-PF=sqrt(PD(Loadi).^2./(QD(Loadi).^2+PD(Loadi).^2));
+    %%
-%%
+    Volt;
-Volt;
+    UAngel*180/3.1415926;
-UAngel*180/3.1415926;
+    %% 通过潮流计算PG
-%% 通过潮流计算PG
+    AngleIJ=sparse(r,c,UAngel(r)-UAngel(c)-Angle',Busnum,Busnum);
-AngleIJ=sparse(r,c,UAngel(r)-UAngel(c)-Angle',Busnum,Busnum);
+    PGBal=PD+diag(Volt)*Y.*cos(AngleIJ)*Volt';
-PGBal=PD+diag(Volt)*Y.*cos(AngleIJ)*Volt';
+    QGBal=QD+diag(Volt)*Y.*sin(AngleIJ)*Volt';
-QGBal=QD+diag(Volt)*Y.*sin(AngleIJ)*Volt';
+    %% 初值-即测量值
-%% 初值-即测量值
+    PG0=PG;
-PG0=PG;
+    QG0=QG;
-QG0=QG;
+    PD0=PD;
-PD0=PD;
+    QD0=QD;
-QD0=QD;
+    %%
-Volt0=Volt;
+    PG0(Balance)=PGBal(Balance);
-UAngel0=UAngel;
+    PG(Balance)=PGBal(Balance);
-%%
+    QG0(Balance)=QGBal(Balance);
-PG0(Balance)=PGBal(Balance);
+    QG0(PVi)=QGBal(PVi);
-PG(Balance)=PGBal(Balance);
+    QG(PVi)=QGBal(PVi);
-QG0(Balance)=QGBal(Balance);
+    %% 真实值
-QG0(PVi)=QGBal(PVi);
+    RealPG=PG0;
-QG(PVi)=QGBal(PVi);
+    RealQG=QG0;
-%% 真实值
+    RealPD=PD0;
-RealPG=PG0;
+    RealQD=QD0;
-RealQG=QG0;
+    %%
-RealPD=PD0;
+    %PGi=zeros(1,1);
-RealQD=QD0;
+    [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,RealPD,RealQD,QD,PD);
-%%
+    Gap=(Init_L*Init_Z'-Init_U*Init_W');
-[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,RealPD,RealQD,QD,PD);
+    KK=0;
-Gap=(Init_L*Init_Z'-Init_U*Init_W');
+    plotGap=zeros(1,60);
-KK=0;
+    ContrlCount=size(Loadi,1)*2+Busnum*2;
-plotGap=zeros(1,60);
+    kmax=60;
-ContrlCount=size(Loadi,1)*2+Busnum*2+Busnum;
+    %% 20120523 临时
-kmax=60;
+    QD_NON_ZERO=QD(PD==0 & QD~=0);
-Precision=Precision/1;
+    QD_NON_ZERO_IND=find(PD==0 & QD~=0);
-%% 加误差
+    %%
-PD0(Loadi)=PD0(Loadi).*(1+normrnd(0,0.05,length(Loadi),1));
+    Precision=Precision/10;
-QD0(Loadi)=QD0(Loadi).*(1+normrnd(0,0.05,length(Loadi),1));
+    
-Vbi=sparse(1*ones(Busnum,1));
+    %% 加误差
-kInit_Z=zeros(length(Init_Z),4);
+    PD0(Loadi)=PD0(Loadi).*(1+normrnd(0,0.05,length(Loadi),1));
-kInit_L=zeros(length(Init_L),4);
+    QD0(Loadi)=QD0(Loadi).*(1+normrnd(0,0.05,length(Loadi),1));
-kInit_W=zeros(length(Init_W),4);
+%     PG0(PGi)=PG0(PGi).*(1+normrnd(0,0.01,length(PGi),1));
-kInit_U=zeros(length(Init_U),4);
+%     QG0(PVi)=QG0(PVi).*(1+normrnd(0,0.01,length(PVi),1));
-kInit_Y=zeros(length(Init_Y),4);
+    %% 读变压器容量
-kPG=zeros(length(PG),4);
+    %[noDataTransNum noDataTransCapacity noDataTransPowerFactor]=ReadNoDataTrans('C:/b/东际911_2751267_2012-09-05/iPso_东际911_2751267_2012-09-05_变压器无负载.txt');
-kQG=zeros(length(QG),4);
+    noDataTransCapacity=0;
-kVolt=zeros(length(Volt),4);
+    
-kUAngel=zeros(length(UAngel),4);
+    while(abs(Gap)>Precision)
-kPD=zeros(length(PD),4);
-kQD=zeros(length(QD),4);
-kVbi=zeros(length(Vbi),4);
-kdeltZ=zeros(length(Init_Z),4);
-kdeltL=zeros(length(Init_L),4);
-kdeltW=zeros(length(Init_W),4);
-kdeltU=zeros(length(Init_U),4);
-kdeltX=zeros(ContrlCount,4);
-kdeltY=zeros(length(Init_Y),4);
-while(abs(Gap)>Precision)
         if KK>kmax
             break;
         end
         plotGap(KK+1)=Gap;
-    for I=1:4
         Init_u=Gap/2/RestraintCount*CenterA;
         AngleIJMat=0;
         %% 开始计算OPF
@@ -93,59 +82,56 @@ while(abs(Gap)>Precision)
         ddg=func_ddg(PGi,PVi,Busnum,RestraintCount,Loadi,PD,QD);
         %%  开始构建deltF
         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,Vbi);
+        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,Loadi);
+        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,RealPD,RealQD,Loadi,KK,PF);
-        Lw=FormLw(Mat_G,Init_U,GenU,Busnum,PVQU,PD,RealPD,RealQD,Loadi,KK,PF);
+        Lw=FormLw(Mat_G,Init_U,GenU,Busnum,PVQU,PD,RealPD,RealQD,Loadi,KK,PF,noDataTransCapacity);
         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);
         %%取各分量
-        [kdeltZ(:,I),kdeltL(:,I),kdeltW(:,I),kdeltU(:,I),kdeltX(:,I),kdeltY(:,I)]=AssignXX(XX,ContrlCount,RestraintCount,Busnum);
+        [deltZ,deltL,deltW,deltU,deltX,deltY]=AssignXX(XX,ContrlCount,RestraintCount,Busnum);
-        if I==1
+        [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);
-            stepT=0.5;
-        end
-        if I==2
-            stepT=0.5;
-        end
-        if I==3
-            stepT=1;
-        end
-        if I<=3
-            [kInit_Z(:,I),kInit_L(:,I),kInit_W(:,I),kInit_U(:,I),kInit_Y(:,I), ...
-                kPG(:,I),kQG(:,I),kVolt(:,I),kUAngel(:,I),kPD(:,I),kQD(:,I),kVbi(:,I)]=Modification(Init_Z,Init_L,Init_W,Init_U,Init_Y,stepT*kdeltZ(:,I),stepT*kdeltL(:,I),stepT*kdeltW(:,I),stepT*kdeltU(:,I),stepT*kdeltX(:,I),stepT*kdeltY(:,I),PG,QG,Volt,UAngel,PVi,ContrlCount,Balance,Busnum,PGi,PD,QD,Loadi,Vbi);
-            Init_Z=kInit_Z(:,I)';
-            Init_L=kInit_L(:,I)';
-            Init_W=kInit_W(:,I)';
-            Init_U=kInit_U(:,I)';
-            Init_Y=kInit_Y(:,I)';
-            PG=kPG(:,I);
-            QG=kQG(:,I);
-            Volt=kVolt(:,I)';
-            UAngel=kUAngel(:,I)';
-            PD=kPD(:,I);
-            QD=kQD(:,I);
-            Vbi=kVbi(:,I);
-        end
-        if I==4
-            deltZ=(kdeltZ(:,1)+kdeltZ(:,2)*2+kdeltZ(:,3)*2+kdeltZ(:,4))/6;
-            deltL=(kdeltL(:,1)+kdeltL(:,2)*2+kdeltL(:,3)*2+kdeltL(:,4))/6;
-            deltW=(kdeltW(:,1)+kdeltW(:,2)*2+kdeltW(:,3)*2+kdeltW(:,4))/6;
-            deltU=kdeltU(:,1)+kdeltU(:,2)*2+kdeltU(:,3)*2+kdeltU(:,4));
-            deltX=kdeltX(:,1)+kdeltX(:,2)*2+kdeltX(:,3)*2+kdeltX(:,4);
-            deltY=kdeltY(:,1)+kdeltY(:,2)*2+kdeltY(:,3)*2+kdeltY(:,4);
-            [Init_Z,Init_L,Init_W,Init_U,Init_Y,PG,QG,Volt,UAngel,PD,QD,Vbi]=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,Vbi);
-        end
-    end
         Gap=(Init_L*Init_Z'-Init_U*Init_W');
         KK=KK+1;
+    end
+    fprintf('迭代次数%d\n',KK);
+    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');
+    %PD=full(PD);
+    %% 统计PD误差
+    % absPDLoad=abs(  (PD(Loadi)-PD0(Loadi))./PD0(Loadi)  );
+    absPDLoad=abs(  (PD(Loadi)-RealPD(Loadi))./RealPD(Loadi)  );
+    maxPDError=max(absPDLoad(absPDLoad<10))
+    Loadi(maxPDError==absPDLoad)
+    absQDLoad=abs(  (QD(Loadi)-RealQD(Loadi))./RealQD(Loadi)  );
+    maxQDError=max(absQDLoad(absQDLoad<10))
+    Loadi(maxQDError==absQDLoad)
+    disp('index');
+    %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);
+    thesis=thesis.AddPDQDPGQG(PD(Loadi),QD(Loadi),PG(Balance),QG(PVi));
+    
 end
-fprintf('%f\n',sum(full(Vbi)));
+PD(Loadi)=thesis.MeanPD();
+QD(Loadi)=thesis.MeanQD();
+PG(Balance)=thesis.MeanPG();
+QG(PVi)=thesis.MeanQG();
+thesis.MaxDeviation(RealPG(Balance),RealQG(PVi),RealPD(Loadi),RealQD(Loadi))
+thesis.StatDeviation(RealPG(Balance),RealQG(PVi),RealPD(Loadi),RealQD(Loadi));
 toc
 

OPF.m

@@ -4,77 +4,71 @@ clear
 %% 存在问题
 % 变压器变比的位置没有考虑，由于现在用的变比都是1，所以没有影响。 20130123
 %% 
-thesis=ForThesis(1,62);
+thesis=ForThesis(1,45);
-[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,Branchi,Branchg,Branchb,Transfork0]= ...
+for II=1:1
-    pf('E:\算例\东际911_2751267_2012-09-05\newFIle20-使用22.txt');
+    [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,Branchi,Branchg,Branchb,Transfork0]= ...
+    pf('E:/算例/柳金Ⅰ926_21671693_2012-09-06/newFIle20 - 计算用2 - (最终用这个做20点的计算).txt');
 %         pf('E:\算例\柳金Ⅰ926_21671693_2012-09-06\newFIle16.txt');
-%pf('D:\Project\青秀降损项目\最小化潮流\最小潮流算例\原始\津头站津视922（3-1）_0.5_120%.txt');
+    %pf('D:\Project\青秀降损项目\最小化潮流\最小潮流算例\原始\津头站津视922（3-1）_0.5_120%.txt');
-%pf('D:\Project\最小化潮流\最小潮流算例\仙海919.txt');
+    %pf('D:\Project\最小化潮流\最小潮流算例\仙海919.txt');
-%pf('c:/file31.txt');
+    %pf('c:/file31.txt'); 
-%% 计算功率因数
+    %% 计算功率因数
-Loadi=QD~=0 | PD~=0;
+    Loadi=QD~=0 | PD~=0;
-PF=sqrt(PD(Loadi).^2./(QD(Loadi).^2+PD(Loadi).^2));
+    PF=sqrt(PD(Loadi).^2./(QD(Loadi).^2+PD(Loadi).^2))
-%%
+    %%
-Volt;
+    Volt;
-UAngel*180/3.1415926;
+    UAngel*180/3.1415926;
-%% 通过潮流计算PG
+    %% 通过潮流计算PG
-AngleIJ=sparse(r,c,UAngel(r)-UAngel(c)-Angle',Busnum,Busnum);
+    AngleIJ=sparse(r,c,UAngel(r)-UAngel(c)-Angle',Busnum,Busnum);
-PGBal=PD+diag(Volt)*Y.*cos(AngleIJ)*Volt';
+    PGBal=PD+diag(Volt)*Y.*cos(AngleIJ)*Volt';
-QGBal=QD+diag(Volt)*Y.*sin(AngleIJ)*Volt';
+    QGBal=QD+diag(Volt)*Y.*sin(AngleIJ)*Volt';
-%% 初值-即测量值
+    %% 初值-即测量值
-PG0=PG;
+    PG0=PG;
-QG0=QG;
+    QG0=QG;
-PD0=PD;
+    PD0=PD;
-QD0=QD;
+    QD0=QD;
-Volt0=Volt;
+    Volt0=Volt;
-UAngel0=UAngel;
+    UAngel0=UAngel;
-%%
+    %%
-PG0(Balance)=PGBal(Balance);
+    PG0(Balance)=PGBal(Balance);
-PG(Balance)=PGBal(Balance);
+    PG(Balance)=PGBal(Balance);
-QG0(Balance)=QGBal(Balance);
+    QG0(Balance)=QGBal(Balance);
-QG0(PVi)=QGBal(PVi);
+    QG0(PVi)=QGBal(PVi);
-QG(PVi)=QGBal(PVi);
+    QG(PVi)=QGBal(PVi);
-%% 真实值
+    %% 真实值
-RealPG=PG0;
+    RealPG=PG0;
-RealQG=QG0;
+    RealQG=QG0;
-RealPD=PD0;
+    RealPD=PD0;
-RealQD=QD0;
+    RealQD=QD0;
-%%
+    %%
-[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,RealPD,RealQD,QD,PD);
+    %PGi=zeros(1,1);
-Gap=(Init_L*Init_Z'-Init_U*Init_W');
+    [Volt,UAngel,Init_Z,Init_W,Init_L,Init_U,Init_Y,PG,QG,RestraintCount,wPG,wQG,wPD,wQD,wVolt,PD,PD0,QD,randPDind,Loadi]=OPF_Init(Busnum,Balance,PG,QG,Volt,GenU,GenL,PVi,PGi,PVQU,PVQL,RealPD,RealQD,QD,PD);
-KK=0;
+    Gap=(Init_L*Init_Z'-Init_U*Init_W');
-plotGap=zeros(1,60);
+    KK=0;
-ContrlCount=size(Loadi,1)*2+Busnum*2+Busnum;
+    plotGap=zeros(1,60);
-kmax=60;
+    ContrlCount=size(Loadi,1)*2+Busnum*2;
-Precision=Precision/1;
+    kmax=60;
-%% 加误差
+    %% 20120523 临时
-PD0(Loadi)=PD0(Loadi).*(1+normrnd(0,0.05,length(Loadi),1));
+    QD_NON_ZERO=QD(PD==0 & QD~=0);
-QD0(Loadi)=QD0(Loadi).*(1+normrnd(0,0.05,length(Loadi),1));
+    QD_NON_ZERO_IND=find(PD==0 & QD~=0);
-Vbi=sparse(1*ones(Busnum,1));
+    %%
-kInit_Z=zeros(length(Init_Z),4);
+    Precision=Precision/1;
-kInit_L=zeros(length(Init_L),4);
+    %% 加误差
-kInit_W=zeros(length(Init_W),4);
+    PD0(Loadi)=PD0(Loadi).*(1+normrnd(0,0.05,length(Loadi),1));
-kInit_U=zeros(length(Init_U),4);
+    QD0(Loadi)=QD0(Loadi).*(1+normrnd(0,0.05,length(Loadi),1));
-kInit_Y=zeros(length(Init_Y),4);
+%     save('20PD0.mat','PD0');
-kPG=zeros(length(PG),4);
+%     save('20QD0.mat','QD0');
-kQG=zeros(length(QG),4);
+    load('20PD0.mat');
-kVolt=zeros(length(Volt),4);
+    load('20QD0.mat');
-kUAngel=zeros(length(UAngel),4);
+%     PG0(PGi)=PG0(PGi).*(1+normrnd(0,0.01,length(PGi),1));
-kPD=zeros(length(PD),4);
+%     QG0(PVi)=QG0(PVi).*(1+normrnd(0,0.01,length(PVi),1));
-kQD=zeros(length(QD),4);
+    %% 读变压器容量
-kVbi=zeros(length(Vbi),4);
+    noDataTransCapacity=0;
-kdeltZ=zeros(length(Init_Z),4);
+    while(abs(Gap)>Precision)
-kdeltL=zeros(length(Init_L),4);
-kdeltW=zeros(length(Init_W),4);
-kdeltU=zeros(length(Init_U),4);
-kdeltX=zeros(ContrlCount,4);
-kdeltY=zeros(length(Init_Y),4);
-while(abs(Gap)>Precision)
         if KK>kmax
             break;
         end
         plotGap(KK+1)=Gap;
-    for I=1:4
         Init_u=Gap/2/RestraintCount*CenterA;
         AngleIJMat=0;
         %% 开始计算OPF
@@ -86,66 +80,69 @@ while(abs(Gap)>Precision)
         L_1Z=diag(Init_Z./Init_L);
         U_1W=diag(Init_W./Init_U);
         %% 形成海森阵
-        deltdeltF=func_deltdeltF(PVi,wPG,wQG,wPD,wQD,ContrlCount);
+        deltdeltF=func_deltdeltF(PVi,wPG,wQG,wVolt,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,PD,QD);
         %%  开始构建deltF
-        deltF=func_deltF(PG,QG,PVi,PGi,wPG,wQG,wPD,wQD,PG0,QG0,PD0,PD,QD,QD0,Busnum,Loadi);
+        deltF=func_deltF(Volt,PG,QG,PVi,PGi,wPG,wQG,wPD,wQD,wVolt,PG0,QG0,PD0,PD,QD,QD0,Volt0,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,Vbi);
+        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,Loadi);
+        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,RealPD,RealQD,Loadi,KK,PF);
-        Lw=FormLw(Mat_G,Init_U,GenU,Busnum,PVQU,PD,RealPD,RealQD,Loadi,KK,PF);
+        Lw=FormLw(Mat_G,Init_U,GenU,Busnum,PVQU,PD,RealPD,RealQD,Loadi,KK,PF,noDataTransCapacity);
         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);
         %%取各分量
-        [kdeltZ(:,I),kdeltL(:,I),kdeltW(:,I),kdeltU(:,I),kdeltX(:,I),kdeltY(:,I)]=AssignXX(XX,ContrlCount,RestraintCount,Busnum);
+        [deltZ,deltL,deltW,deltU,deltX,deltY]=AssignXX(XX,ContrlCount,RestraintCount,Busnum);
-        if I==1
+        [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);
-            stepT=0.5;
-        end
-        if I==2
-            stepT=0.5;
-        end
-        if I==3
-            stepT=1;
-        end
-        if I<=3
-            [kInit_Z(:,I),kInit_L(:,I),kInit_W(:,I),kInit_U(:,I),kInit_Y(:,I), ...
-                kPG(:,I),kQG(:,I),kVolt(:,I),kUAngel(:,I),kPD(:,I),kQD(:,I),kVbi(:,I)]=Modification(Init_Z,Init_L,Init_W,Init_U,Init_Y,stepT*kdeltZ(:,I),stepT*kdeltL(:,I),stepT*kdeltW(:,I),stepT*kdeltU(:,I),stepT*kdeltX(:,I),stepT*kdeltY(:,I),PG,QG,Volt,UAngel,PVi,ContrlCount,Balance,Busnum,PGi,PD,QD,Loadi,Vbi);
-            Init_Z=kInit_Z(:,I)';
-            Init_L=kInit_L(:,I)';
-            Init_W=kInit_W(:,I)';
-            Init_U=kInit_U(:,I)';
-            Init_Y=kInit_Y(:,I)';
-            PG=kPG(:,I);
-            QG=kQG(:,I);
-            Volt=kVolt(:,I)';
-            UAngel=kUAngel(:,I)';
-            PD=kPD(:,I);
-            QD=kQD(:,I);
-            Vbi=kVbi(:,I);
-        end
-        if I==4
-            deltZ=(kdeltZ(:,1)+kdeltZ(:,2)*2+kdeltZ(:,3)*2+kdeltZ(:,4))/6;
-            deltL=(kdeltL(:,1)+kdeltL(:,2)*2+kdeltL(:,3)*2+kdeltL(:,4))/6;
-            deltW=(kdeltW(:,1)+kdeltW(:,2)*2+kdeltW(:,3)*2+kdeltW(:,4))/6;
-            deltU=(kdeltU(:,1)+kdeltU(:,2)*2+kdeltU(:,3)*2+kdeltU(:,4))/6;
-            deltX=(kdeltX(:,1)+kdeltX(:,2)*2+kdeltX(:,3)*2+kdeltX(:,4))/6;
-            deltY=(kdeltY(:,1)+kdeltY(:,2)*2+kdeltY(:,3)*2+kdeltY(:,4))/6;
-            [Init_Z,Init_L,Init_W,Init_U,Init_Y,PG,QG,Volt,UAngel,PD,QD,Vbi]=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,Vbi);
-        end
-    end
         Gap=(Init_L*Init_Z'-Init_U*Init_W');
         KK=KK+1;
+    end
+    fprintf('迭代次数%d\n',KK);
+    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');
+    %PD=full(PD);
+    %% 统计PD误差
+    % absPDLoad=abs(  (PD(Loadi)-PD0(Loadi))./PD0(Loadi)  );
+    absPDLoad=abs(  (PD(Loadi)-RealPD(Loadi))./RealPD(Loadi)  );
+    maxPDError=max(absPDLoad(absPDLoad<10))
+    Loadi(maxPDError==absPDLoad)
+    absQDLoad=abs(  (QD(Loadi)-RealQD(Loadi))./RealQD(Loadi)  );
+    maxQDError=max(absQDLoad(absQDLoad<10))
+    Loadi(maxQDError==absQDLoad)
+    disp('index');
+    %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);
+    thesis=thesis.AddPDQDPGQG(PD(Loadi),QD(Loadi),PG(Balance),QG(PVi));
+    
 end
-fprintf('%f\n',sum(full(Vbi)));
+PD(Loadi)=thesis.MeanPD();
+QD(Loadi)=thesis.MeanQD();
+PG(Balance)=thesis.MeanPG();
+QG(PVi)=thesis.MeanQG();
+thesis.MaxDeviation(RealPG(Balance),RealQG(PVi),RealPD(Loadi),RealQD(Loadi))
+thesis.StatDeviation(RealPG(Balance),RealQG(PVi),RealPD(Loadi),RealQD(Loadi))
+thesis.MaxDeviation(PG0(Balance),QG0(PVi),PD0(Loadi),QD0(Loadi))
+thesis.StatDeviation(PG0(Balance),QG0(PVi),PD0(Loadi),QD0(Loadi))
+% thesis.SqareDeviation(RealPG(Balance(1)),RealQG(PVi(1)),RealPD(Loadi),RealQD(Loadi))
+thesis.SqareDeviation(PG0(Balance(1)),QG0(PVi(1)),PD0(Loadi),QD0(Loadi))
+thesis.PercentOfPass(PG0(Balance(1)),QG0(PVi(1)),PD0(Loadi),QD0(Loadi))
+thesis.MaxBranchDeviation(Linei,Linej,Liner,Linex,Lineb,Transfori,Transforj,Transforr,Transforx,Branchi,Branchg,Branchb,Transfork0,Volt0,UAngel0,Volt,UAngel,'E:\算例\柳金Ⅰ926_21671693_2012-09-06\newFIle20.txt',PD0,QD0)
+thesis.StatBranchDeviation(Linei,Linej,Liner,Linex,Lineb,Transfori,Transforj,Transforr,Transforx,Branchi,Branchg,Branchb,Transfork0,Volt0,UAngel0,Volt,UAngel)
 toc
 

OPFWrongDataCheck.asv

@@ -0,0 +1,106 @@
+clc
+    
+    clear
+ for II=1:16   
+     tic
+    [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('D:\Project\最小化潮流\最小潮流算例\金湖924（2-1）_0.5_85%.txt');
+    %pf('c:/file31.txt');
+    %pf('ieee118PG.dat');
+    
+    %% 计算功率因数
+    atan(PD(QD~=0 | PD~=0)./QD(QD~=0 | PD~=0));
+    Volt;
+    UAngel*180/3.1415926;
+    %% 通过潮流计算PG
+    AngleIJ=sparse(r,c,UAngel(r)-UAngel(c)-Angle',Busnum,Busnum);
+    PGBal=diag(Volt)*Y.*cos(AngleIJ)*Volt';
+    
+    %% 初值-即测量值
+    PG0=PG;
+    PD0=PD;
+    
+    PDReal=PD;%真值
+    
+    %%
+    PG0(Balance)=PGBal(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);
+    PD0(Loadi(II))=PD0(Loadi(II))*(1+0.018);
+    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;
+    kmax=60;
+    %% 20120523 临时
+    QD_NON_ZERO=QD(PD==0 & QD~=0);
+    QD_NON_ZERO_IND=find(PD==0 & QD~=0);
+    %%
+    while(abs(Gap)>Precision)
+        if KK>kmax
+            break;
+        end
+        plotGap(KK+1)=Gap;
+        Init_u=Gap/2/RestraintCount*CenterA;
+        AngleIJMat=0;
+        %% 开始计算OPF
+        %% 形成等式约束的雅克比
+        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(PVi,wG,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);
+        
+        %% 形成方程矩阵
+        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,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);
+        Lw=FormLw(Mat_G,Init_U,GenU,Busnum,PVQU,PD0,Loadi,KK);
+        Lx=FormLx(deltF,deltH,Init_Y,deltG,Init_Z,Init_W);
+        YY=FormYY(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,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);
+        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);
+    %DrawGap(plotGap);
+    %%
+    %Volt=full(Volt');
+    %PD=full(PD);
+    %% 统计PD误差
+    % absPDLoad=abs(  (PD(Loadi)-PD0(Loadi))./PD0(Loadi)  );
+    absPDLoad=abs(  (PD(Loadi)-PDReal(Loadi))./PDReal(Loadi)  );
+    maxPDError=max(absPDLoad);
+    %disp('index')
+    LoadiArray=Loadi(absPDLoad==maxPDError);
+    if length(LoadiArray)>1
+        disp('没找出')
+        Loadi(II)
+    end
+    if length(LoadiArray)==1
+        if LoadiArray~=Loadi(II)
+            disp('没找出')
+            Loadi(II)
+        end
+    end
+    toc;
+    
+end
+

OPF_Init.asv

@@ -0,0 +1,34 @@
+function [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)
+Loadi=find(QD~=0 | PD~=0);
+%Loadi=[1:Busnum]';
+RestraintCount=size(Loadi,1)*2+Busnum*1; %约束条件数,放开所有QD
+t_Bal_volt=Volt(Balance);
+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=1*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));% 无功上限
+tPL=sparse(GenL(:,2));% 发电机有功下限
+tQL=sparse(PVQL(:,1));% 无功下限
+wPG=0;
+wQG=0;
+%randInt=randperm(size(Loadi,1));
+%randPDind=randInt(1:10);
+randPDind=0;
+wPD=1/.05^2*ones(size(Loadi,1),1);
+% wPD(1:2:end)=0;
+wQD=1/.05^2*zeros(size(Loadi,1),1);
+% wQD(1:2:end)=0;
+%wD(randPDind)=0;%一些负荷不约束
+%wD(7)=0;
+% wD(11)=0;
+PD=0.5*PD0;
+%powerFacter=0.98;
+%QD=PD*sqrt((1-powerFacter^2)/powerFacter^2);
+QD=.5*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,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)
+function [Volt,UAngel,Init_Z,Init_W,Init_L,Init_U,Init_Y,PG,QG,RestraintCount,wPG,wQG,wPD,wQD,wVolt,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(Loadi,1)*2+Busnum*2+Busnum; %约束条件数,放开所有QD
+RestraintCount=size(Loadi,1)*2+Busnum*1; %约束条件数,放开所有QD
 t_Bal_volt=Volt(Balance);
 Volt=sparse(1*ones(1,Busnum));
 Volt(Balance)=t_Bal_volt;
@@ -20,15 +20,16 @@ wQG=0;
 %randInt=randperm(size(Loadi,1));
 %randPDind=randInt(1:10);
 randPDind=0;
-wPD=zeros(size(Loadi,1),1);
+wPD=1/.05^2*ones(size(Loadi,1),1);
-wPD(1:2:end)=0;
+% wPD(1:2:end)=0;
-wQD=zeros(size(Loadi,1),1);
+wQD=1/.05^2*ones(size(Loadi,1),1);
-wQD(1:2:end)=0;
+% wQD(1:2:end)=0;
 %wD(randPDind)=0;%一些负荷不约束
 %wD(7)=0;
 % wD(11)=0;
-PD=0.8*PD0;
+PD=0.5*PD0;
 %powerFacter=0.98;
 %QD=PD*sqrt((1-powerFacter^2)/powerFacter^2);
-QD=.8*QD0;
+QD=.5*QD0;
+wVolt=1/.05^2*ones(size(Busnum,1),1);
 end

SolveIt.asv

@@ -0,0 +1,42 @@
+function 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)
+LxComa=FormLxComa(deltF,deltG,deltH,Init_L,Luu,Lul,Init_Z,Init_Y,Lz,Init_U,Init_W,Lw,Lx);
+H=-deltdeltF+ddh;%+ddg*(Init_Z'+Init_W');
+t1=diag(Init_L.\Init_Z-Init_U.\Init_W);
+t2=-deltG*( t1 )*deltG';
+aa=[
+    (H+t2),deltH;
+    deltH',zeros(size(Init_Y,2));
+    ];
+yy=[LxComa;-Ly];
+%% ƽºâ½Úµãµçѹ²»±ä 
+t=size(PVi,1)+size(PGi,1)+size(Loadi,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)),ContrlCount+2*Busnum,ContrlCount+2*Busnum);
+deltG(t+Balance,:)=0;
+%%
+t=size(PVi,1)+size(PGi,1)+size(Loadi,1)+Busnum*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)),ContrlCount,ContrlCount);
+deltG(t+Balance,:)=0;
+%%
+dxdy=aa\yy;
+dX=dxdy(1:ContrlCount);
+dY=dxdy(ContrlCount+1:ContrlCount+2*Busnum);
+dL=Lz+deltG'*dX;
+dU=-Lw-deltG'*dX;
+dZ=-diag(Init_L)\Lul-diag(Init_L)\diag(Init_Z)*dL;
+dW=-diag(Init_U)\Luu-diag(Init_U)\diag(Init_W)*dU;
+XX=[
+    dX;
+    dY;
+    dZ;
+    dW;
+    dL;
+    dU;
+    
+];
+end

func_ddg.asv

@@ -0,0 +1,9 @@
+function ddg=func_ddg(PGi,PVi,Busnum,RestraintCount,Loadi,PD,QD)
+
+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_ddh.m

@@ -51,7 +51,6 @@ t=[ddPdVdV+ddQdVdV,ddPdVdT+ddQdVdT;
 sizeLoadi=size(Loadi,1)*2;
 ddh=[
     sparse(sizeLoadi,ContrlCount);
-    sparse(2*Busnum,sizeLoadi),-t,sparse(2*Busnum,Busnum);
+    sparse(2*Busnum,sizeLoadi),-t;
-    sparse(Busnum,ContrlCount);
     ];
 end

func_deltF.asv

@@ -0,0 +1,19 @@
+function deltF=func_deltF(PG,PVi,GenC,ContrlCount,PGi,wG,wD,PG0,PD0,PD,Busnum)
+%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
+t2=2*wD.*(PD-PD0);
+deltF=[
+    zeros(size(PGi));
+    zeros(size(PVi));
+    t2;
+    zeros(2*Busnum,1);
+];
+
+end

func_deltF.m

@@ -1,10 +1,11 @@
-function deltF=func_deltF(PG,QG,PVi,PGi,wPG,wQG,wPD,wQD,PG0,QG0,PD0,PD,QD,QD0,Busnum,Loadi)
+function deltF=func_deltF(Volt,PG,QG,PVi,PGi,wPG,wQG,wPD,wQD,wVolt,PG0,QG0,PD0,PD,QD,QD0,Volt0,Busnum,Loadi)
 t3=2*wPD.*(PD(Loadi)-PD0(Loadi));
 t4=2*wQD.*(QD(Loadi)-QD0(Loadi));
-deltF=[sparse(length(Loadi),1);
+t5=2*wVolt.*(Volt-Volt0);
-    sparse(length(Loadi),1);
+deltF=[sparse(t3);
-    sparse(2*Busnum,1);
+    sparse(t4);
-    sparse(ones(Busnum,1));
+    sparse(t5');
+    sparse(1*Busnum,1);
 ];
 
 end

func_deltG.asv

@@ -1,34 +1,46 @@
 function deltG=func_deltG(Busnum,PVi,PGi,Loadi,PD,QD)
+%%
+sizePGi=size(PGi,1);
+sizePVi=size(PVi,1);
 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,length(Loadi));
+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,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),length(Loadi));
+dg4_dPD=sparse(size(Loadi,1),1);
 dg5_dPD=sparse(size(Loadi,1),Busnum);
-dg6_dPD=dg5_dPD;
-dg7_dPD=sparse(sizeLoadi,Busnum);
 %% 
-dg3_dQD=sparse(length(Loadi),length(Loadi));
+dg1_dQD=sparse(size(Loadi(1),1),size(PGi,1));
-dg4_dQD=sparse(1:size(Loadi,1),1:size(Loadi,1),ones(size(Loadi,1),1),size(Loadi,1),size(Loadi,1));
+dg2_dQD=sparse(size(Loadi(1),1),size(PVi,1));
-dg5_dQD=sparse(size(Loadi,1),Busnum);
+dg3_dQD=sparse(length(Loadi(1)),length(Loadi));
-dg6_dQD=dg5_dQD;
+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));
-dg7_dQD=sparse(sizeLoadi,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,length(Loadi));
+dg4_dx=sparse(2*Busnum,1);
 dg5_dx=[sparse(1:Busnum,1:Busnum,ones(Busnum,1),Busnum,Busnum);
     sparse(Busnum,Busnum);
     ];
-dg6_dx=dg5_dx;
-dg7_dx=sparse(2*Busnum,Busnum);
 %% 
-dg3_dvbi=sparse(Busnum,sizeLoadi);
+deltG=[dg1_dPg,dg2_dPg,dg3_dPg,dg4_dPg,dg5_dPg;
-dg4_dvbi=sparse(Busnum,length(Loadi));
+       dg1_dQr,dg2_dQr,dg3_dQr,dg4_dQr,dg5_dQr;
-dg5_dvbi=sparse(Busnum,Busnum);
+       dg1_dPD,dg2_dPD,dg3_dPD,dg4_dPD,dg5_dPD;
-dg6_dvbi=sparse(Busnum,Busnum);
+       dg1_dQD,dg2_dQD,dg3_dQD,dg4_dQD,dg5_dQD;
-dg6_dvbi=sparse(Busnum,Busnum);
+       dg1_dx,dg2_dx,dg3_dx,dg4_dx,dg5_dx;
-%%
-deltG=[dg3_dPD,dg4_dPD,dg5_dPD,dg6_dPD;
-       dg3_dQD,dg4_dQD,dg5_dQD,dg6_dQD;
-       dg3_dx,dg4_dx,dg5_dx,dg6_dx;
-       dg3_dvbi,dg4_dvbi,dg5_dvbi,dg6_dvbi;
 ];
 end

func_deltG.m

@@ -1,34 +1,29 @@
 function deltG=func_deltG(Busnum,PVi,PGi,Loadi,PD,QD)
+%%
+
 sizeLoadi=size(Loadi,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),length(Loadi));
 dg5_dPD=sparse(size(Loadi,1),Busnum);
-dg6_dPD=dg5_dPD;
-dg7_dPD=sparse(sizeLoadi,Busnum);
 %% 
 dg3_dQD=sparse(length(Loadi),length(Loadi));
 dg4_dQD=sparse(1:size(Loadi,1),1:size(Loadi,1),ones(size(Loadi,1),1),size(Loadi,1),size(Loadi,1));
 dg5_dQD=sparse(size(Loadi,1),Busnum);
-dg6_dQD=dg5_dQD;
+
-dg7_dQD=sparse(sizeLoadi,Busnum);
 %% 
 dg3_dx=sparse(2*Busnum,sizeLoadi);
 dg4_dx=sparse(2*Busnum,length(Loadi));
 dg5_dx=[sparse(1:Busnum,1:Busnum,ones(Busnum,1),Busnum,Busnum);
     sparse(Busnum,Busnum);
     ];
-dg6_dx=dg5_dx;
-dg7_dx=sparse(2*Busnum,Busnum);
 %% 
-dg3_dvbi=sparse(Busnum,sizeLoadi);
+
-dg4_dvbi=sparse(Busnum,length(Loadi));
+deltG=[dg3_dPD,dg4_dPD,dg5_dPD;
-dg5_dvbi=sparse(-eye(Busnum,Busnum));
+       dg3_dQD,dg4_dQD,dg5_dQD;
-dg6_dvbi=sparse(eye(Busnum,Busnum));
+       dg3_dx,dg4_dx,dg5_dx;
-dg7_dvbi=sparse(eye(Busnum,Busnum));
-%%
-deltG=[dg3_dPD,dg4_dPD,dg5_dPD,dg6_dPD,dg7_dPD;
-       dg3_dQD,dg4_dQD,dg5_dQD,dg6_dQD,dg7_dQD;
-       dg3_dx,dg4_dx,dg5_dx,dg6_dx,dg7_dx;
-       dg3_dvbi,dg4_dvbi,dg5_dvbi,dg6_dvbi,dg7_dvbi;
 ];
 end

func_deltH.asv

@@ -0,0 +1,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_dx = jacobian_M3(Busnum,Volt,Y,Angle,UAngel,r,c);  %形成雅克比矩阵
+deltH=[dH_dPg;dH_dQr;dH_dPD;dH_dx'];
+end

func_deltH.m

@@ -2,7 +2,6 @@ function deltH=func_deltH(Busnum,Volt,PVi,Y,PGi,UAngel,r,c,Angle,Loadi)
 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_dx = jacobian_M(Busnum,Volt,Y,Angle,UAngel,r,c);  %形成雅克比矩阵
-dH_dvbi=sparse(Busnum,2*Busnum);
+deltH=[dH_dPD;dH_dQD;dH_dx'];
-deltH=[dH_dPD;dH_dQD;dH_dx';dH_dvbi];
 
 end

func_deltdeltF.asv

@@ -0,0 +1,11 @@
+function deltdeltF=func_deltdeltF(Busnum,PVi,PGi,wG,wD,ContrlCount)
+
+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);
+deltdeltF=[
+    diagC*2,sparse(sizeC,ContrlCount-sizeC);
+    sparse(ContrlCount-sizeC,ContrlCount);
+    ];
+end

func_deltdeltF.m

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

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

openfile2.asv

@@ -0,0 +1,84 @@
+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)
+%**************************************************************************
+%    程序简介 : 子函数——读取潮流计算所需数据
+%    编    者：
+%    编制时间 ：2010.12
+%**************************************************************************
+data = dlmread(FileName);                                           % 一次读入全部数据
+zeroRow = find(data(:,1)==0); 
+Busnum= data(1,1);                                                       % 节点数
+PQstandard = data(1,3);                                               % 基准容量
+kmax = data(1,4);                                                           %最大迭代次数
+Precision = data(1,4);                                                    % 精度
+%Balance = data(3,2);     
+Balance=data(3:zeroRow(2)-1,2);% 生成1到节点号的列向量
+%CenterA=data(1,5);  %中心参数
+%LineNum=data(1,2); %支路数
+Base=data(1,3);
+%% 各参数矩阵分块
+
+line = data(zeroRow(2)+1:zeroRow(3)-1,:);               % 形成线路参数矩阵
+ground = data(zeroRow(5)+1:zeroRow(6)-1,:);         % 形成对地支路参数矩阵
+tran = data(zeroRow(3)+1:zeroRow(4)-1,:);              % 形成变压器参数矩阵
+buspq = data(zeroRow(8)+1:zeroRow(9)-1,:);          % 形成节点功率参数矩阵
+PV = data(zeroRow(11)+1:zeroRow(12)-1,:);                % 形成pv节点功率参数矩阵
+GenP=data(zeroRow(9)+1:zeroRow(10)-1,:);
+GenQ=data(zeroRow(11)+1:zeroRow(12)-1,:);
+%% 线路参数矩阵分块
+Linei = line(:,2);                                                             % 节点i
+Linej= line(:,3);                                                              % 节点j
+Liner = line(:,4);                                                             % 线路电阻
+Linex = line(:,5);                                                            % 线路电抗
+Lineb = line(:,6);                                                            % b/2
+%% 对地支路参数矩阵
+Branchi = ground(:,2);                                                  % 对地支路节点号
+Branchb = ground(:,4);                                                 % 对地支路的导纳
+%% 变压器参数矩阵
+Transfori = tran(:,3);                                                     % 节点i
+Transforj= tran(:,4);                                                      % 节点j
+Transforr = tran(:,5);                                                     % 变压器电阻
+Transforx= tran(:,6);                                                     % 变压器电抗
+Transfork0 = tran(:,7);                                                  % 变压器变比
+%% 节点功率参数矩阵
+Pointpoweri = buspq(:,3);
+PG=buspq(:,5);                                                             % 发电机有功
+QG=buspq(:,6);                                                             % 发电机无功
+PD=buspq(:,7);                                                              % 负荷有功
+QD=buspq(:,8);                                                             % 负荷无功
+%%除以基值
+PG=PG/Base;
+QG=QG/Base;
+PD=PD/Base;
+QD=QD/Base;
+%%
+PD=sparse(PD);
+QD=sparse(QD);
+PG=sparse(PG);
+QG=sparse(QG);
+%% pv节点功率参数矩阵
+PVi = PV(:,3); % PV节点的节点号
+PVu = PV(:,5); % PV节点电压
+PVQL=PV(:,6);%PV节点无功下限
+PVQL=PVQL/Base;
+PVQU=PV(:,7); %PV节点无功上限
+PVQU=PVQU/Base;
+%% 发电机参数
+%GenU=Gen(:,[1 5 6]);
+%GenL=Gen(:,[1 7 8]);
+GenC=GenP(:,[3 7:9]);
+t=GenC(:,2);
+GenC(:,2)=GenC(:,4);
+GenC(:,4)=t;
+%%%%%%%%%%%%%%%%%%%%
+%GenC(:,2:4)=100*GenC(:,2:4);
+t=GenP(:,[3 5]);
+%GenL=[t,PVQL(PVi)];
+GenL=t;%有功下界
+GenL(:,2)=GenL(:,2)/Base;
+t=GenP(:,[3 6]);
+%GenU=[t,PVQU(PVi)];
+GenU=t;%有功上届
+GenU(:,2)=GenU(:,2)/Base;
+PGi=GenP(:,3);%发电机节点号
+end

pf.asv

@@ -0,0 +1,34 @@
+function [kmax,Precision,Uangle,U,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(FileName)
+%**************************************************************************
+%    程序名称：电力系统潮流计算程序
+%    程序算法：极坐标下的牛顿-拉夫逊法
+%    程序功能：主函数
+%    程序编者：
+%    编制时间：2010.12
+%**************************************************************************
+clc;
+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);
+%% 形成节点导纳矩阵
+[GB,Y,r,c,Angle] = admmatrix(Busnum,Linei,Linej,Liner,Linex,Lineb,Transfori,Transforj,Transforr,...
+    Transforx,Transfork0,Branchi,Branchb);
+[P0,Q0,U,Uangle] = Initial(PG,PD,PQstandard,Pointpoweri,QG,QD,Busnum);                         %求功率不平衡量
+disp('迭代次数i  最大不平衡量');
+%% 循环体计算
+for i = 0:kmax
+    [Jacob,PQ,U,Uangle] = jacobian(Busnum,Balance,PVi,PVu,U,Uangle,Y,Angle,P0,Q0,r,c);  %形成雅克比矩阵
+   % disp('第一次雅克比');
+    m = max(abs(PQ));
+    m=full(m);
+    fprintf('       %u         %.8f \n',i,m);
+    if m > Precision                                                                                                                                 %判断不平衡量是否满足精度要求
+        [Uangle,U] = solvefun(Busnum,Jacob,PQ,Uangle,U);                                                           %求解修正方程，更新电压变量
+    else
+        disp(['收敛，迭代次数为',num2str(i),'次']);
+        break                                                                                                                                             %若满足精度要求，则计算收敛
+    end
+end
+toc;
+end

plotVolt.asv

@@ -0,0 +1,26 @@
+%% 画电压，为了写论文用
+load('case1V.mat');
+CV1=Volt;
+load('case2V.mat');
+CV2=Volt;
+load('case3V.mat');
+CV3=Volt;
+load('caseR.mat');
+CVR=Volt;
+load('caseM.mat');
+CVM=Volt;
+subplot(3,2,1)
+hist(CVR)
+title('真实值')
+
+
+subplot(2,2,3)
+hist(CV1)
+title('Case 1')
+subplot(2,2,4)
+hist(CV2)
+title('Case 2')
+subplot(2,2,5)
+hist(CV3)
+title('Case 3')
+

plotVolt.m

@@ -5,8 +5,6 @@ load('case2V.mat');
 CV2=Volt;
 load('case3V.mat');
 CV3=Volt;
-load('case4V.mat');
-CV4=Volt;%负荷曲线偏差大
 load('caseR.mat');
 CVR=Volt;
 load('caseM.mat');
@@ -25,22 +23,18 @@ subplot(3,2,3)
 hist(CV1)
 xlabel('电压/p.u');
 ylabel('数量/个');
-title('100%可知情形电压分布')
+title('情形1电压分布')
 subplot(3,2,4)
 hist(CV2)
 xlabel('电压/p.u');
 ylabel('数量/个');
-title('50%可知情形电压分布')
+title('情形2电压分布')
 subplot(3,2,5)
 hist(CV3)
 xlabel('电压/p.u');
 ylabel('数量/个');
-title('始端覆盖情形电压分布')
+title('情形3电压分布')
-subplot(3,2,6)
+
-hist(CV4)
-xlabel('电压/p.u');
-ylabel('数量/个');
-title('典型负荷曲线偏差大情形电压分布')
 
 % figure()
 % hold on