完成了QD，并修复了完成PD时的一个求导错误。

Signed-off-by: dmy <dugg@21cn.com>
1.加入.gitignore
2014-05-14 16:37:56 +08:00 · 2014-05-14 11:55:15 +08:00 · 2014-05-14 10:12:34 +08:00 · 2014-05-14 10:00:41 +08:00
23 changed files with 133 additions and 514 deletions
--- a/.gitignore
+++ b/.gitignore
@@ -0,0 +1 @@
+*.asv
--- a/DrawGap.asv
+++ b/DrawGap.asv
@@ -1,6 +0,0 @@
-function DrawGap(plotGap)
-x=find(plotGap);
-ts=size(x,2);
-
-
-end
--- a/FormAA.asv
+++ b/FormAA.asv
@@ -1,12 +0,0 @@
-function AA=FormAA(L_1Z,deltG,U_1W,Hcoma,deltH)
-tOnes=eye(14);
-tZeros=zeros(14);
-AA=[
-    tOnes,L_1Z,tZeros,tZeros,tZeros,zeros(14,10);
-    tZeros,tOnes,tZeros,tZeros,-deltG',zeros(14,10);
-    tZeros,tZeros,tOnes,U_1W,tZeros,zeros(14,10);
-    tZeros,tZeros,tZeros,tOnes,deltG,zeros(14,10);
-    tZeros,tZeros,tZeros,tZeros,Hcoma,deltH;
-    zeros(10,14),zeros(10,14),zeros(10,14),zeros(10,14),deltH',zeros(10,10);
-    ];
-end
--- a/FormAlphaP.asv
+++ b/FormAlphaP.asv
@@ -1,7 +0,0 @@
-function AlphaP=FormAlphaP(Init_L,deltL,Init_U,deltU)
-ti=deltL(delt)
-t1=-Init_L./deltL';
-t2=-Init_U./deltU';
-t3=[t1,t2];
-t4=t3()
-end
--- a/FormG.asv
+++ b/FormG.asv
@@ -1,11 +0,0 @@
-function Mat_G=FormG(Volt,PVi,PGi,PG,QG,PD,QD,Loadi)
-
-Mat_G=[
-    PG(PGi);
-    QG(PVi);
-    sparse(PD(Loadi));
-    sparse(QD(Loadi));
-    Volt';
-    
-    ];
-end
--- a/FormG.m
+++ b/FormG.m
@@ -1,10 +1,14 @@
-function Mat_G=FormG(Volt,PVi,PGi,PG,QG,PD,QD,Loadi,Vbi)
+function Mat_G=FormG(Volt,PD,QD,Loadi,Vbi,mVolt,RealVolt,bigM,sigma,PDbi,QDbi,RealPD,RealQD)

 Mat_G=[
-    sparse(PD(Loadi));
-    sparse(QD(Loadi));
-    Volt'-Vbi;
-    Volt'+Vbi;
+    sparse(PD(Loadi))-bigM*PDbi-RealPD(Loadi)*sigma;
+    sparse(PD(Loadi))+bigM*PDbi+RealPD(Loadi)*sigma;
+    sparse(QD(Loadi))-bigM*QDbi-RealQD(Loadi)*sigma;
+    sparse(QD(Loadi))+bigM*QDbi+RealQD(Loadi)*sigma;
+    Volt'-bigM*Vbi-mVolt'-sigma*RealVolt';
+    Volt'+bigM*Vbi+mVolt'+sigma*RealVolt';
    Vbi;
+    PDbi;
+    QDbi;
    ];
 end
--- a/FormH.asv
+++ b/FormH.asv
@@ -1,31 +0,0 @@
-function Mat_H=FormH(Busnum,GB,AngleIJMat,Volt,PG,PD,QG,QD,Y)
-t1=real(GB).*cos(AngleIJMat)+imag(GB).*sin(AngleIJMat);
-t2=Volt'*Volt;
-t3=t1.*t2;
-t4=sum(-t3,2);%P
-t5=real(GB).*sin(AngleIJMat)-imag(GB).*cos(AngleIJMat);
-t6=t2.*t5;
-t7=sum(-t6,2);%Q
-t8=PG-PD;
-t9=QG-QD;
-%Mat_H=([(PG-PD)',(QG-QD)'])'+([t4',t7'])';
-Mat_H(1:2:2*Busnum)=t8(1:Busnum)+t4(1:Busnum);
-Mat_H(2:2:2*Busnum)=t9(1:Busnum)+t7(1:Busnum);
-Mat_H=Mat_H';
-%% 
-QDcos=textread('300glys.txt');
-t=QD(PD==0 &&);
-aa=QD;
-QD(QD~=0)=PD(QD~=0)./tan(QDcos);
-QD(PD==0)=t;
-%% 
-%%%%一下是学姐给的公式
-AngleIJ=AngleIJMat-angle(GB);
-%dP=PG-PD-diag(Volt)*Y*cos(AngleIJ)*Volt';
-dP=PG-PD-diag(Volt)*Y.*cos(AngleIJ)*Volt';
-dQ=QG-QD-diag(Volt)*Y.*sin(AngleIJ)*Volt';
-%Mat_H(1:2:2*Busnum)=dP(1:Busnum);暂时改一下 20111227
-%Mat_H(2:2:2*Busnum)=dQ(1:Busnum);暂时改一下 20111227
-Mat_H=[dP;dQ;];
-
-end
--- a/FormLw.m
+++ b/FormLw.m
@@ -1,11 +1,7 @@
-function Lw=FormLw(Mat_G,Init_U,GenU,Busnum,PVQU,PD,PD0,QD0,Loadi,KK,PF,noDataTransCapacity)
+function [Lw,UpperLimit]=FormLw(Mat_G,Init_U,Busnum,PD0,QD0,Loadi,bigM)
 KK=999;
-%PU=GenU(:,2);%·¢µç»úÓÐ¹¦ÉÏ½ç
-
 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;
@@ -13,7 +9,6 @@ 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);
@@ -22,9 +17,8 @@ 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=([0*PDU',bigM*ones(1,length(Loadi)),0*QDU',bigM*ones(1,length(Loadi)),0*VoltU,bigM*ones(1,Busnum),1*ones(1,Busnum),1*ones(1,length(Loadi)*2)])';
+UpperLimit=t1;
 t2=Mat_G+Init_U'-t1;
 Lw=t2;

--- a/FormLz.m
+++ b/FormLz.m
@@ -1,8 +1,5 @@
-function Lz=FormLz(Mat_G,Init_L,GenL,Busnum,PVQL,PD,PD0,QD0,Loadi,KK,PF)
-KK=999;
-
+function [Lz,LoweLimit]=FormLz(Mat_G,Init_L,Busnum,PD0,QD0,Loadi,bigM)
 VoltL=(0.9)*ones(1,Busnum);
-%VoltL=-10*ones(1,Busnum);
 PDL=PD0(Loadi);
 PDL(PDL>0)=0.800*PDL(PDL>0);
 PDL(PDL<0)=1.200*PDL(PDL<0);
@@ -11,7 +8,6 @@ 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);
@@ -20,10 +16,9 @@ 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',-20000*VoltL,-20000*VoltL,0*ones(1,Busnum)])';
+%t1=([-1*ones(1,length(Loadi)),0*PDL',-1*ones(1,length(Loadi)),0*QDL',-5*ones(1,length(VoltL)),0*VoltL,0*ones(1,Busnum),0*ones(1,length(Loadi)*2)])';
+t1=([-bigM*ones(1,length(Loadi)),0*PDL',-bigM*ones(1,length(Loadi)),0*QDL',-bigM*ones(1,length(VoltL)),0*VoltL,0*ones(1,Busnum),0*ones(1,length(Loadi)*2)])';
+LoweLimit=t1;
 t2=Mat_G-Init_L'-t1;
 Lz=t2;
-
 end
--- a/FormYY.asv
+++ b/FormYY.asv
@@ -1,5 +0,0 @@
-function FormYY(Init_L,Lul.Lz,Init_U,Luu,Lw,Lz,LxComa)
-t=[
-   -    
-]
-end
--- a/Modification.asv
+++ b/Modification.asv
@@ -1,18 +0,0 @@
-function [Init_Z,Init_L,Init_W,Init_U,Init_Y,deltZ,deltL,deltW,deltU,deltX,deltY,PG,QG,Volt,UAngel]=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)
-AlphaP=FormAlphaP(Init_L,deltL,Init_U,deltU);
-fprintf('AlphaP %f\n',AlphaP);
-AlphaD=FormAlphaD(Init_Z,deltZ,Init_W,deltW);
-fprintf('AlphaD %f\n',AlphaD);
-
-Init_Z=Init_Z+AlphaD*deltZ';
-Init_L=Init_L+AlphaP*deltL';
-Init_W=Init_W+AlphaD*deltW';
-Init_U=Init_U+AlphaP*deltU';
-Init_Y=Init_Y+AlphaD*deltY';
-PG(PVi)=PG(PGi)+deltX(size(PGi,1));
-QG(PVi)=QG(PVi)+deltX(size(PGi,1)+1,1:(size(PVi,1)+size(PGi,1)));
-t=deltX(size(PVi,1)+size(PGi,1)+1,ContrlCount)';
-t(2*Balance-1)=0;
-Volt=Volt+AlphaP*t(2:2:2*Busnum);
-UAngel=UAngel+AlphaP*t(1:2:2*Busnum);
-end
--- a/Modification.m
+++ b/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,Vbi,PDbi,QDbi]=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,PDbi,QDbi)
 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,7 @@ 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)';
+Vbi=Vbi+AlphaP*t(2*Busnum+1:2*Busnum+Busnum)';
+PDbi=PDbi+AlphaP*t(2*Busnum+Busnum+1:2*Busnum+Busnum+length(Loadi))';
+QDbi=QDbi+AlphaP*t(2*Busnum+Busnum+length(Loadi)+1:end)';
 end
--- a/OPF.asv
+++ b/OPF.asv
@@ -1,151 +0,0 @@
-tic
-clc
-clear
-%% 存在问题
-% 变压器变比的位置没有考虑，由于现在用的变比都是1，所以没有影响。 20130123
-%%
-thesis=ForThesis(1,62);
-[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:\算例\东际911_2751267_2012-09-05\newFIle20-使用22.txt');
-%         pf('E:\算例\柳金Ⅰ926_21671693_2012-09-06\newFIle16.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;
-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=PD+diag(Volt)*Y.*cos(AngleIJ)*Volt';
-QGBal=QD+diag(Volt)*Y.*sin(AngleIJ)*Volt';
-%% 初值-即测量值
-PG0=PG;
-QG0=QG;
-PD0=PD;
-QD0=QD;
-Volt0=Volt;
-UAngel0=UAngel;
-%%
-PG0(Balance)=PGBal(Balance);
-PG(Balance)=PGBal(Balance);
-QG0(Balance)=QGBal(Balance);
-QG0(PVi)=QGBal(PVi);
-QG(PVi)=QGBal(PVi);
-%% 真实值
-RealPG=PG0;
-RealQG=QG0;
-RealPD=PD0;
-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');
-KK=0;
-plotGap=zeros(1,60);
-ContrlCount=size(Loadi,1)*2+Busnum*2+Busnum;
-kmax=60;
-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));
-Vbi=sparse(1*ones(Busnum,1));
-kInit_Z=zeros(length(Init_Z),4);
-kInit_L=zeros(length(Init_L),4);
-kInit_W=zeros(length(Init_W),4);
-kInit_U=zeros(length(Init_U),4);
-kInit_Y=zeros(length(Init_Y),4);
-kPG=zeros(length(PG),4);
-kQG=zeros(length(QG),4);
-kVolt=zeros(length(Volt),4);
-kUAngel=zeros(length(UAngel),4);
-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
-        %% 形成等式约束的雅克比
-        deltH=func_deltH(Busnum,Volt,PVi,Y,PGi,UAngel,r,c,Angle,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,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);
-        %% 形成方程矩阵
-        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_H=FormH(Busnum,Volt,PG,PD,QG,QD,Y,UAngel,r,c,Angle,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);
-        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);
-        if I==1
-            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('%f\n',sum(full(Vbi)));
-toc
-
--- a/OPF.m
+++ b/OPF.m
@@ -6,7 +6,7 @@ clear
 %%
 thesis=ForThesis(1,62);
 [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:\算例\东际911_2751267_2012-09-05\newFIle20-使用22.txt');
+    pf('E:\算例\东际911_2751267_2012-09-05\newFIle9-2.txt');
 %         pf('E:\算例\柳金Ⅰ926_21671693_2012-09-06\newFIle16.txt');
 %pf('D:\Project\青秀降损项目\最小化潮流\最小潮流算例\原始\津头站津视922（3-1）_0.5_120%.txt');
 %pf('D:\Project\最小化潮流\最小潮流算例\仙海919.txt');
@@ -44,108 +44,67 @@ RealQD=QD0;
 Gap=(Init_L*Init_Z'-Init_U*Init_W');
 KK=0;
 plotGap=zeros(1,60);
-ContrlCount=size(Loadi,1)*2+Busnum*2+Busnum;
-kmax=60;
+ContrlCount=size(Loadi,1)*2+Busnum*2+Busnum+length(Loadi)*2;
+kmax=600;
 Precision=Precision/1;
+sigma=0.03;
 %% 加误差
-PD0(Loadi)=PD0(Loadi).*(1+normrnd(0,0.05,length(Loadi),1));
-QD0(Loadi)=QD0(Loadi).*(1+normrnd(0,0.05,length(Loadi),1));
+PD0(Loadi)=PD0(Loadi).*(1+normrnd(0,sigma,length(Loadi),1));
+QD0(Loadi)=QD0(Loadi).*(1+normrnd(0,sigma,length(Loadi),1));
+mPD=PD0;
+mQD=QD0;
+RealVolt=Volt0;
+mVolt=Volt0.*(1+normrnd(0,sigma,1,length(Volt0)));
 Vbi=sparse(1*ones(Busnum,1));
-kInit_Z=zeros(length(Init_Z),4);
-kInit_L=zeros(length(Init_L),4);
-kInit_W=zeros(length(Init_W),4);
-kInit_U=zeros(length(Init_U),4);
-kInit_Y=zeros(length(Init_Y),4);
-kPG=zeros(length(PG),4);
-kQG=zeros(length(QG),4);
-kVolt=zeros(length(Volt),4);
-kUAngel=zeros(length(UAngel),4);
-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);
+PDbi=sparse(1*ones(length(Loadi),1));
+QDbi=sparse(1*ones(length(Loadi),1));
+bigM=5;
 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
-        %% 形成等式约束的雅克比
-        deltH=func_deltH(Busnum,Volt,PVi,Y,PGi,UAngel,r,c,Angle,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,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);
-        %% 形成方程矩阵
-        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_H=FormH(Busnum,Volt,PG,PD,QG,QD,Y,UAngel,r,c,Angle,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);
-        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);
-        if I==1
-            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
+    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,PD,QD);
+    %%
+    L_1Z=diag(Init_Z./Init_L);
+    U_1W=diag(Init_W./Init_U);
+    %% 形成海森阵
+    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,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,PD,QD,Loadi,Vbi,mVolt,RealVolt,bigM,sigma,PDbi,QDbi,RealPD,RealQD);
+    Mat_H=FormH(Busnum,Volt,PG,PD,QG,QD,Y,UAngel,r,c,Angle,Loadi);
+    Ly=Mat_H;
+    [Lz,LowerLimit]=FormLz(Mat_G,Init_L,Busnum,PD0,QD0,Loadi,bigM);%下界
+    [Lw,UpperLimit]=FormLw(Mat_G,Init_U,Busnum,PD0,QD0,Loadi,bigM);%上界
+    Lx=FormLx(deltF,deltH,Init_Y,deltG,Init_Z,Init_W);
+    if any(LowerLimit>UpperLimit)
+        warning('不等式设置错误');
    end
+    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);
+    %%取各分量
+    [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,QD,Vbi,PDbi,QDbi]=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,PDbi,QDbi);
    Gap=(Init_L*Init_Z'-Init_U*Init_W');
    KK=KK+1;
 end
-fprintf('%f\n',sum(full(Vbi)));
+fprintf('目标函数: %f\n',sum(full(Vbi)));
+fprintf('目标函数: %f\n',sum(full(QDbi)));
 toc

--- a/OPF_Init.m
+++ b/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)
 Loadi=find(QD~=0 | PD~=0);
 %Loadi=[1:Busnum]';
-RestraintCount=size(Loadi,1)*2+Busnum*2+Busnum; %约束条件数,放开所有QD
+RestraintCount=size(Loadi,1)*2+size(Loadi,1)*2+Busnum*2+Busnum+length(Loadi)*2; %约束条件数,放开所有QD
 t_Bal_volt=Volt(Balance);
 Volt=sparse(1*ones(1,Busnum));
 Volt(Balance)=t_Bal_volt;
--- a/admmatrix.asv
+++ b/admmatrix.asv
@@ -1,31 +0,0 @@
-function [G,B,GB,Y,r,c,Angle]  = admmatrix(Busnum,Linei,Linej,Liner,Linex,Lineb,Transfori...
-    ,Transforj,Transforr,Transforx,Transfork0,Branchi,Branchb)
-%**************************************************************************
-%  程序功能 : 子函数——形成节点导纳矩阵Y
-%  编    者： 
-%  编制时间：2010.12
-%**************************************************************************
-%% 支路导纳计算
-G = -sparse(Linei,Linej,Liner./(Liner.^2+Linex.^2),Busnum,Busnum) - sparse(Linej,Linei,Liner./(Liner.^2+Linex.^2),Busnum,Busnum);
-G = G - sparse(1:Busnum,1:Busnum,sum(G,2)');                                                                                  % 计算各线路支路电导
-B = sparse(Linei,Linej,Linex./(Liner.^2+Linex.^2),Busnum,Busnum)+sparse(Linej,Linei,Linex./(Liner.^2+Linex.^2),Busnum,Busnum);
-B = B - sparse(1:Busnum,1:Busnum,sum(B,2)')+sparse(Linei,Linei,Lineb,Busnum,Busnum)+sparse(Linej,Linej,Lineb,Busnum,Busnum);
-%% 变压器支路计算
-if Transfori>0 
-    mr = Transforr./(Transforr.^2+Transforx.^2);                                                                                  % 计算变压器支路电导
-    mx = -Transforx./(Transforr.^2+Transforx.^2);                                                                                % 计算变压器支路电纳
-    G = G-sparse(Transfori,Transforj,mr./Transfork0,Busnum,Busnum)-sparse(Transforj,Transfori,mr./Transfork0,Busnum,Busnum)...
-        +sparse(Transfori,Transfori,mr./Transfork0./Transfork0,Busnum,Busnum)+sparse(Transforj,Transforj,mr,Busnum,Busnum);
-    B = B-sparse(Transfori,Transforj,mx./Transfork0,Busnum,Busnum)-sparse(Transforj,Transfori,mx./Transfork0,Busnum,Busnum)...
-        +sparse(Transfori,Transfori,mx./Transfork0./Transfork0,Busnum,Busnum)+sparse(Transforj,Transforj,mx,Busnum,Busnum);
-end
-%% 接地支路计算
-if Branchi>0                                                                 % 判断有无接地支路
-    B = B+sparse(Branchi,Branchi,Branchb,Busnum,Busnum);
-end
-%% 化作极坐标形式
-GB = G+B.*1i;                                                              %将电导，电纳合并，写成复数形式
-Y = abs(GB);                                                                %求节点导纳幅值
-[r,c] = find(Y);
-Angle = angle(GB(GB~=0));                                      %求节点导纳角度
-%Angle=angle(GB);
--- a/func_ddh.m
+++ b/func_ddh.m
@@ -51,7 +51,7 @@ t=[ddPdVdV+ddQdVdV,ddPdVdT+ddQdVdT;
 sizeLoadi=size(Loadi,1)*2;
 ddh=[
    sparse(sizeLoadi,ContrlCount);
-    sparse(2*Busnum,sizeLoadi),-t,sparse(2*Busnum,Busnum);
-    sparse(Busnum,ContrlCount);
+    sparse(2*Busnum,sizeLoadi),-t,sparse(2*Busnum,Busnum+length(Loadi)*2);
+    sparse(Busnum+length(Loadi)*2,ContrlCount);
    ];
 end
--- a/func_deltF.m
+++ b/func_deltF.m
@@ -5,6 +5,7 @@ deltF=[sparse(length(Loadi),1);
    sparse(length(Loadi),1);
    sparse(2*Busnum,1);
    sparse(ones(Busnum,1));
+    sparse(ones(length(Loadi),1));
+    sparse(ones(length(Loadi),1));
 ];
-
 end
--- a/func_deltG.asv
+++ b/func_deltG.asv
@@ -1,34 +0,0 @@
-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));
-dg5_dvbi=sparse(Busnum,Busnum);
-dg6_dvbi=sparse(Busnum,Busnum);
-dg6_dvbi=sparse(Busnum,Busnum);
-%%
-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
--- a/func_deltG.m
+++ b/func_deltG.m
@@ -1,34 +1,83 @@
 function deltG=func_deltG(Busnum,PVi,PGi,Loadi,PD,QD)
 sizeLoadi=size(Loadi,1);
+%dg3 PD-M*b-t-PD0
 dg3_dPD=sparse(1:size(Loadi,1),1:size(Loadi,1),ones(size(Loadi,1),1),size(Loadi,1),size(Loadi,1));
+%dg32 PD+M*B+t+PD0
+dg32_dPD=dg3_dPD;
+%dg4 QD-M*B-t-QD0
 dg4_dPD=sparse(size(Loadi,1),length(Loadi));
+%dg42 QD+M*B+t+QD0
+dg42_dPD=sparse(size(Loadi,1),length(Loadi));
 dg5_dPD=sparse(size(Loadi,1),Busnum);
 dg6_dPD=dg5_dPD;
+%dg7 ÊÇ VbiµÄÔ¼Êø
 dg7_dPD=sparse(sizeLoadi,Busnum);
+dpdbi_dPD=sparse(sizeLoadi,sizeLoadi);
+dqdbi_dPD=sparse(sizeLoadi,sizeLoadi);
 %% 
 dg3_dQD=sparse(length(Loadi),length(Loadi));
+dg32_dQD=sparse(sizeLoadi,sizeLoadi);
 dg4_dQD=sparse(1:size(Loadi,1),1:size(Loadi,1),ones(size(Loadi,1),1),size(Loadi,1),size(Loadi,1));
+dg42_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);
+dpdbi_dQD=sparse(sizeLoadi,sizeLoadi);
+dqdbi_dQD=sparse(sizeLoadi,sizeLoadi);
 %% 
 dg3_dx=sparse(2*Busnum,sizeLoadi);
+dg32_dx=sparse(2*Busnum,sizeLoadi);
 dg4_dx=sparse(2*Busnum,length(Loadi));
+dg42_dx=sparse(2*Busnum,sizeLoadi);
 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);
+dpdbi_dx=sparse(2*Busnum,sizeLoadi);
+dqdbi_dx=sparse(2*Busnum,sizeLoadi);
 %%
 dg3_dvbi=sparse(Busnum,sizeLoadi);
+dg32_dvbi=sparse(Busnum,sizeLoadi);
 dg4_dvbi=sparse(Busnum,length(Loadi));
+dg42_dvbi=sparse(Busnum,sizeLoadi);
 dg5_dvbi=sparse(-eye(Busnum,Busnum));
 dg6_dvbi=sparse(eye(Busnum,Busnum));
 dg7_dvbi=sparse(eye(Busnum,Busnum));
+dpdbi_dvbi=sparse(Busnum,sizeLoadi);
+dqdbi_dvbi=sparse(Busnum,sizeLoadi);
 %%
-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;
+dg3_dpdbi=sparse(-eye(sizeLoadi,sizeLoadi));
+dg32_dpdbi=sparse(eye(sizeLoadi,sizeLoadi));
+dg4_dpdbi=sparse(sizeLoadi,sizeLoadi);
+dg42_dpdbi=sparse(sizeLoadi,sizeLoadi);
+dg5_dpdbi=sparse(sizeLoadi,Busnum);
+dg6_dpdbi=sparse(sizeLoadi,Busnum);
+dg7_dpdbi=sparse(sizeLoadi,Busnum);
+dpdbi_dpdbi=sparse(eye(sizeLoadi,sizeLoadi));
+dqdbi_dpdbi=sparse(sizeLoadi,sizeLoadi);
+%%
+dg3_dqdbi=sparse(sizeLoadi,sizeLoadi);
+dg32_dqdbi=sparse(sizeLoadi,sizeLoadi);
+dg4_dqdbi=sparse(-eye(sizeLoadi,sizeLoadi));
+dg42_dqdbi=sparse(eye(sizeLoadi,sizeLoadi));
+dg5_dqdbi=sparse(sizeLoadi,Busnum);
+dg6_dqdbi=sparse(sizeLoadi,Busnum);
+dg7_dqdbi=sparse(sizeLoadi,Busnum);
+dpdbi_dqdbi=sparse(sizeLoadi,sizeLoadi);
+dqdbi_dqdbi=sparse(eye(sizeLoadi,sizeLoadi));
+%%
+deltG=[dg3_dPD,dg32_dPD,dg4_dPD,dg42_dPD,dg5_dPD,dg6_dPD,dg7_dPD,dpdbi_dPD,dqdbi_dPD;
+       dg3_dQD,dg32_dQD,dg4_dQD,dg42_dQD,dg5_dQD,dg6_dQD,dg7_dQD,dpdbi_dQD,dqdbi_dQD;
+       dg3_dx,dg32_dx,dg4_dx,dg42_dx,dg5_dx,dg6_dx,dg7_dx,dpdbi_dx,dqdbi_dx;
+       dg3_dvbi,dg32_dvbi,dg4_dvbi,dg42_dvbi,dg5_dvbi,dg6_dvbi,dg7_dvbi,dpdbi_dvbi,dqdbi_dvbi;
+       dg3_dpdbi,dg32_dpdbi,dg4_dpdbi,dg42_dpdbi,dg5_dpdbi,dg6_dpdbi,dg7_dpdbi,dpdbi_dpdbi,dqdbi_dpdbi;
+       dg3_dqdbi,dg32_dqdbi,dg4_dqdbi,dg42_dqdbi,dg5_dqdbi,dg6_dqdbi,dg7_dqdbi,dpdbi_dqdbi,dqdbi_dqdbi;
 ];
+% deltG=[dg3_dPD,dg32_dPD,dg4_dPD,dg5_dPD,dg6_dPD,dg7_dPD;
+%        dg3_dQD,dg32_dQD,dg4_dQD,dg5_dQD,dg6_dQD,dg7_dQD;
+%        dg3_dx,dg32_dx,dg4_dx,dg5_dx,dg6_dx,dg7_dx;
+%        dg3_dvbi,dg32_dvbi,dg4_dvbi,dg5_dvbi,dg6_dvbi,dg7_dvbi;
+%        dg3_dpdbi,dg32_dpdbi,dg4_dpdbi,dg5_dpdbi,dg6_dpdbi,dg7_dpdbi;
+% ];
 end
--- a/func_deltH.m
+++ b/func_deltH.m
@@ -3,6 +3,8 @@ dH_dPD=[sparse(1:size(Loadi,1),Loadi,-ones(size(Loadi,1),1),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';dH_dvbi];
+dH_dpdbi=sparse(length(Loadi),2*Busnum);
+dH_dqdbi=sparse(length(Loadi),2*Busnum);
+deltH=[dH_dPD;dH_dQD;dH_dx';dH_dvbi;dH_dpdbi;dH_dqdbi];

 end
--- a/modifyadmmatrix.asv
+++ b/modifyadmmatrix.asv
@@ -1,18 +0,0 @@
-function [new_G,new_B,GB,Y,r,c,Angle]  = modifyadmmatrix(ii,jj,G,GG,B,BB)
-%**************************************************************************
-%  程序功能 : 子函数——形成节点导纳矩阵Y
-%  编    者： 
-%  编制时间：2010.12
-%**************************************************************************
-%% 支路导纳计算
-new_G=G;
-new_G(ii,jj)=new_G(ii,jj)-G(ii,jj);
-new_G(jj,ii)=new_G(jj,ii)-G(jj,ii);
-new_G(ii,ii)=new_G(ii,ii)+G(ii,jj);
-new_G(jj,jj)=new_G(jj,jj)+G(ii,jj);
-
-%% 化作极坐标形式
-GB = G+B.*1i;                                                              %将电导，电纳合并，写成复数形式
-Y = abs(GB);                                                                %求节点导纳幅值
-[r,c] = find(Y);
-Angle = angle(GB(GB~=0));                                      %求节点导纳角度
--- a/openfile.asv
+++ b/openfile.asv
@@ -1,64 +0,0 @@
-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,Gen,GenU,GenL,GenC,CenterA,PGi,PVQU,PVQL] = openfile(FileName)
-%**************************************************************************
-%    程序简介 : 子函数——读取潮流计算所需数据
-%    编    者：
-%    编制时间 ：2010.12
-%**************************************************************************
-data = dlmread(FileName);                                           % 一次读入全部数据
-Busnum= data(1,1);                                                       % 节点数
-PQstandard = data(1,3);                                               % 基准容量
-kmax = data(1,4);                                                           %最大迭代次数
-Precision = data(2,1);                                                    % 精度
-Balance = data(3,2);                                                      % 生成1到节点号的列向量
-CenterA=data(1,5);  %中心参数
-LineNum=data(1,2); %支路数
-%% 各参数矩阵分块
-zeroRow = find(data(:,1)==0);                                     %查找第一列元素为零的行号
-line = data(zeroRow(1)+1:zeroRow(2)-1,:);               % 形成线路参数矩阵
-ground = data(zeroRow(2)+1:zeroRow(3)-1,:);         % 形成对地支路参数矩阵
-tran = data(zeroRow(3)+1:zeroRow(4)-1,:);              % 形成变压器参数矩阵
-buspq = data(zeroRow(4)+1:zeroRow(5)-1,:);          % 形成节点功率参数矩阵
-PV = data(zeroRow(5)+1:zeroRow(6)-1,:);                % 形成pv节点功率参数矩阵
-Gen=data(zeroRow(6)+1:zeroRow(7)-1,:);
-%% 线路参数矩阵分块
-Linei = line(:,2);                                                             % 节点i
-Linej= line(:,3);                                                              % 节点j
-Liner = line(:,4);                                                             % 线路电阻
-Linex = line(:,5);                                                            % 线路电抗
-Lineb = line(:,6);                                                            % b/2
-%% 对地支路参数矩阵
-Branchi = ground(:,1);                                                  % 对地支路节点号
-Branchb = ground(:,2);                                                 % 对地支路的导纳
-%% 变压器参数矩阵
-Transfori = tran(:,2);                                                     % 节点i
-Transforj= tran(:,3);                                                      % 节点j
-Transforr = tran(:,4);                                                     % 变压器电阻
-Transforx= tran(:,5);                                                     % 变压器电抗
-Transfork0 = tran(:,6);                                                  % 变压器变比
-%% 节点功率参数矩阵
-Pointpoweri = buspq(:,1);
-PG=buspq(:,2);                                                             % 发电机有功
-QG=buspq(:,3);                                                             % 发电机无功
-PD=buspq(:,4);                                                              % 负荷有功
-QD=buspq(:,5);                                                             % 负荷无功
-%% pv节点功率参数矩阵
-PVi = PV(:,1); % PV节点的节点号
-PVu = PV(:,2); % PV节点电压
-PVQL=PV(:,3);%PV节点无功下限
-PVQU=PV(:,4); %PV节点无功上限
-%% 发电机参数
-%GenU=Gen(:,[1 5 6]);
-%GenL=Gen(:,[1 7 8]);
-GenC=Gen(:,[1 2:4]);
-t=GenC(:,2);
-GenC(:,2)=GenC(:,4);
-GenC(:,4)=t;
-t=Gen(:,[1 5]);
-%GenL=[t,PVQL(PVi)];
-GenL=t;%有功下界
-t=Gen(:,[1 6]);
-%GenU=[t,PVQU(PVi)];
-GenU=t;%有功上届
-PGi=Gen(:,1);%发电机节点号
-end
Author	SHA1	Message	Date
dmy	b071a58e33	完成了QD，并修复了完成PD时的一个求导错误。 Signed-off-by: dmy <dugg@21cn.com>	2014-05-14 16:37:56 +08:00
dmy	17f274f818	1.加入.gitignore 2.完成了PD Signed-off-by: dmy <dugg@21cn.com>	2014-05-14 11:55:15 +08:00
dmy	ac0bcf1da8	复制了一遍PD的约束 Signed-off-by: dmy <dugg@21cn.com>	2014-05-14 10:12:34 +08:00
dmy	7bde55ceb8	完成了电压的部分 Signed-off-by: dmy <dugg@21cn.com>	2014-05-14 10:00:41 +08:00