diff --git a/FormG.m b/FormG.m
index d8e17e8..e336b9a 100644
--- a/FormG.m
+++ b/FormG.m
@@ -1,4 +1,4 @@
-function Mat_G=FormG(Volt,PVi,PGi,PG,QG,PD)
+function Mat_G=FormG(Volt,PVi,PGi,PG,QG,PD,Loadi)
 %t1=PG(PVi);
 %GP=t1;%发电机P
 %GP=[4.5 4.5]';
@@ -13,7 +13,7 @@ Mat_G=[
     %GP;
     PG(PGi);
     QG(PVi);
-    sparse(PD);
+    sparse(PD(Loadi));
     %GQ(PVi);
     %[0 1.45]';
     Volt';
diff --git a/FormH.m b/FormH.m
index f38d2e8..e0f68a6 100644
--- a/FormH.m
+++ b/FormH.m
@@ -5,7 +5,6 @@ QD(QD~=0)=PD(QD~=0)./tan(QDcos);
 QD(QD_NON_ZERO_IND)=QD_NON_ZERO;
 %% 
 AngleIJ=sparse(r,c,UAngel(r)-UAngel(c)-Angle',Busnum,Busnum);
-%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';
 
diff --git a/FormLw.m b/FormLw.m
index 5071ee3..c2312dd 100644
--- a/FormLw.m
+++ b/FormLw.m
@@ -1,14 +1,14 @@
-function Lw=FormLw(Mat_G,Init_U,GenU,Busnum,PVQU,PD0)
+function Lw=FormLw(Mat_G,Init_U,GenU,Busnum,PVQU,PD0,randPDind,Loadi)
 
 PU=1*GenU(:,2);%发电机有功上界
 QU=PVQU(:,1);%发电机无功上界
 VoltU=1.1*ones(1,Busnum);
-PDU=PD0;
-PDU(PD0>0)=1.2*PD0(PD0>0);
+PDU=PD0(Loadi);
+PDU(PDU>0)=1.2*PDU(PDU>0);
+PDU(PDU<0)=0.8*PDU(PDU<0);
 PDU(PDU==0)=.2;
-PDU(PDU<0)=0.8*PD0(PD0<0);
-%PDU=25*ones(Busnum,1)';
 t1=([PU',QU',PDU',VoltU])';
 t2=Mat_G+Init_U'-t1;
 Lw=t2;
+
 end
\ No newline at end of file
diff --git a/FormLz.asv b/FormLz.asv
index 5e7ddd4..b4e5e7c 100644
--- a/FormLz.asv
+++ b/FormLz.asv
@@ -1,10 +1,28 @@
-function Lz=FormLz(Mat_G,Init_L,GenL,LinePLimt,LineCount,Busnum,PVQL)
+function Lz=FormLz(Mat_G,Init_L,GenL,Busnum,PVQL,PD0,randPDind,Loadi)
 
-PL=GenL(:,2);%发电机有功下界
+PL=1*GenL(:,2);%发电机有功下界
 QL=PVQL(:,1);%发电机无功下界
 VoltL=0.9*ones(1,Busnum);
-t1=([PL',QL',VoltL])';
+PDL=PD0(Loadi);
+PDL(PD0>0)=0.8*PD0(PDL>0);
+PDL(PD0==0)=-.2;
+PDL(PDL<0)=1.2*PD0(PDL<0);
+t1=([PL',QL',PDL',VoltL])';
 
 t2=Mat_G-Init_L'-t1;
 Lz=t2;
+
+% PL=1*GenL(:,2);%发电机有功下界
+% QL=PVQL(:,1);%发电机无功下界
+% VoltL=0.9*ones(1,Busnum);
+% PDL=PD0;
+% PDL(PD0>0)=0.8*PD0(PD0>0);
+% PDL(PD0==0)=-.2;
+% PDL(PD0<0)=1.2*PD0(PD0<0);
+% 
+% t1=([PL',QL',PDL',VoltL])';
+% 
+% t2=Mat_G-Init_L'-t1;
+% Lz=t2;
+
 end
\ No newline at end of file
diff --git a/FormLz.m b/FormLz.m
index 4da82f8..b85ce45 100644
--- a/FormLz.m
+++ b/FormLz.m
@@ -1,15 +1,14 @@
-function Lz=FormLz(Mat_G,Init_L,GenL,Busnum,PVQL,PD0)
+function Lz=FormLz(Mat_G,Init_L,GenL,Busnum,PVQL,PD0,randPDind,Loadi)
 
 PL=1*GenL(:,2);%发电机有功下界
 QL=PVQL(:,1);%发电机无功下界
 VoltL=0.9*ones(1,Busnum);
-PDL=PD0;
-PDL(PD0>0)=0.8*PD0(PD0>0);
-PDL(PD0==0)=-.2;
-PDL(PD0<0)=1.2*PD0(PD0<0);
-
+PDL=PD0(Loadi);
+PDL(PDL>0)=0.8*PDL(PDL>0);
+PDL(PDL<0)=1.2*PDL(PDL<0);
+PDL(PDL==0)=-.2;
 t1=([PL',QL',PDL',VoltL])';
-
 t2=Mat_G-Init_L'-t1;
 Lz=t2;
+
 end
\ No newline at end of file
diff --git a/Modification.m b/Modification.m
index 9dfbfa6..65f89ef 100644
--- 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]=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)
+function [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)
 AlphaP=FormAlphaP(Init_L,deltL,Init_U,deltU);
 fprintf('AlphaP %f\n',full(AlphaP));
 AlphaD=FormAlphaD(Init_Z,deltZ,Init_W,deltW);
@@ -13,9 +13,9 @@ Init_Y=Init_Y+AlphaD*deltY';
 PG(PGi)=PG(PGi)+AlphaP*deltX(1:size(PGi,1));
 %QG(PVi)=QG(PVi)+deltX(size(PGi,1)+1:size(PVi,1)+size(PGi,1) );
 QG(PVi)=QG(PVi)+AlphaP*deltX(size(PGi,1)+1:size(PVi,1)+size(PGi,1) );
-t=deltX(size(PVi,1)+size(PGi,1)+1:size(PVi,1)+size(PGi,1)+Busnum);
-PD=PD+AlphaP*t;
-t=deltX(size(PVi,1)+size(PGi,1)+Busnum+1:ContrlCount)';
+t=deltX(size(PVi,1)+size(PGi,1)+1:size(PVi,1)+size(PGi,1)+size(Loadi,1));
+PD(Loadi)=PD(Loadi)+AlphaP*t;
+t=deltX(size(PVi,1)+size(PGi,1)+size(Loadi,1)+1:ContrlCount)';
 t(Busnum+Balance)=0;
 %Volt=Volt+AlphaP*t(2:2:2*Busnum);暂时改一下20111227
 %UAngel=UAngel+AlphaP*t(1:2:2*Busnum);暂时改一下20111227
diff --git a/OPF.asv b/OPF.asv
index 59e3631..73b28d1 100644
--- a/OPF.asv
+++ b/OPF.asv
@@ -1,18 +1,14 @@
 tic
 clear
-%[kmax,Precision,UAngel,Volt,Busnum,PVi,PVu,Balance,Y,Angle,P0,Q0,r,c,GB,Linei,Linej,Transfori,Transforj,GenU,GenL,GenC,LineLimti,LineLimtj,LinePLimt,PG,QG,PD,QD,CenterA,LineCount,PGi,PVQU,PVQL]=pf('5sj.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]=pf('ieee3001PG.dat');
-GB=full(GB);
 %PVi电压节点序号
 %PVu电压节点电压标幺值
 Volt;
 UAngel*180/3.1415926;
-%sprintf('%f\n',Volt);
-%sprintf('%f\n',Angel);
 %% 初值
 PG0=PG;
 PD0=PD;
-[Volt,UAngel,Init_Z,Init_W,Init_L,Init_U,Init_Y,PG,QG,RestraintCount,wG,wD,PD]=OPF_Init(Busnum,Balance,PG,QG,Volt,GenU,GenL,PVi,PGi,PVQU,PVQL,PD0);
+[Volt,UAngel,Init_Z,Init_W,Init_L,Init_U,Init_Y,PG,QG,RestraintCount,wG,wD,PD,randPDind]=OPF_Init(Busnum,Balance,PG,QG,Volt,GenU,GenL,PVi,PGi,PVQU,PVQL,PD0);
 Gap=(Init_L*Init_Z'-Init_U*Init_W');
 KK=0;
 plotGap=zeros(1,50);
@@ -28,10 +24,10 @@ while(abs(Gap)>Precision)
     end
     plotGap(KK+1)=Gap;
     Init_u=Gap/2/RestraintCount*CenterA;
-    AngleIJMat=repmat(UAngel',1,Busnum)-repmat(UAngel,Busnum,1);
+    AngleIJMat=0;
     %% 开始计算OPF
     %% 形成等式约束的雅克比
-    deltH=func_deltH(Busnum,Volt,PVi,AngleIJMat,Y,GB,PGi);
+    deltH=func_deltH(Busnum,Volt,PVi,AngleIJMat,Y,GB,PGi,UAngel,r,c,Angle);
     %% 形成不等式约束的雅克比
     deltG=func_deltG(Busnum,PVi,PGi);
     %%
@@ -40,7 +36,7 @@ while(abs(Gap)>Precision)
     %% 形成海森阵
     deltdeltF=func_deltdeltF(Busnum,GenC,PVi,PGi,wG,wD,PD0);
     %% 形成ddHy
-    ddh=func_ddh3(AngleIJMat,GB,Volt,Init_Y,Busnum,PVi,PGi,Y);
+    ddh=func_ddh3(AngleIJMat,GB,Volt,Init_Y,Busnum,PVi,PGi,Y,UAngel,r,c,Angle);
     %%  开始构建ddg
     ddg=func_ddg(PGi,PVi,Busnum,RestraintCount);
     %%  开始构建deltF
@@ -50,7 +46,7 @@ while(abs(Gap)>Precision)
     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);
-    Mat_H=FormH(Busnum,GB,AngleIJMat,Volt,PG,PD,QG,QD,Y,QD_NON_ZERO,QD_NON_ZERO_IND);
+    Mat_H=FormH(Busnum,GB,AngleIJMat,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);
     Lw=FormLw(Mat_G,Init_U,GenU,Busnum,PVQU,PD0);
@@ -65,7 +61,6 @@ while(abs(Gap)>Precision)
     KK=KK+1;
 end
 fprintf('迭代次数%d\n',KK);
-%CalCost(GenC,PG,PGi);
 ObjectiveFun(PG,PG0,PGi,PD,PD0,wG,wD)
 DrawGap(plotGap);
 Volt=Volt';
diff --git a/OPF.m b/OPF.m
index 44668c6..dc3c6f8 100644
--- a/OPF.m
+++ b/OPF.m
@@ -8,11 +8,11 @@ UAngel*180/3.1415926;
 %% 初值
 PG0=PG;
 PD0=PD;
-[Volt,UAngel,Init_Z,Init_W,Init_L,Init_U,Init_Y,PG,QG,RestraintCount,wG,wD,PD]=OPF_Init(Busnum,Balance,PG,QG,Volt,GenU,GenL,PVi,PGi,PVQU,PVQL,PD0);
+[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);
 Gap=(Init_L*Init_Z'-Init_U*Init_W');
 KK=0;
 plotGap=zeros(1,50);
-ContrlCount=size(PVi,1)+size(PGi,1)+Busnum*3;
+ContrlCount=size(PVi,1)+size(PGi,1)+size(Loadi,1)+Busnum*2;
 kmax=600;
 %% 20120523 临时
 QD_NON_ZERO=QD(PD==0 & QD~=0);
@@ -27,42 +27,44 @@ while(abs(Gap)>Precision)
     AngleIJMat=0;
     %% 开始计算OPF
     %% 形成等式约束的雅克比
-    deltH=func_deltH(Busnum,Volt,PVi,AngleIJMat,Y,GB,PGi,UAngel,r,c,Angle);
+    deltH=func_deltH(Busnum,Volt,PVi,AngleIJMat,Y,GB,PGi,UAngel,r,c,Angle,Loadi);
     %% 形成不等式约束的雅克比
-    deltG=func_deltG(Busnum,PVi,PGi);
+    deltG=func_deltG(Busnum,PVi,PGi,Loadi);
     %%
     L_1Z=diag(Init_Z./Init_L);
     U_1W=diag(Init_W./Init_U);
     %% 形成海森阵
     deltdeltF=func_deltdeltF(Busnum,GenC,PVi,PGi,wG,wD,PD0);
     %% 形成ddHy
-    ddh=func_ddh3(AngleIJMat,GB,Volt,Init_Y,Busnum,PVi,PGi,Y,UAngel,r,c,Angle);
+    ddh=func_ddh3(AngleIJMat,GB,Volt,Init_Y,Busnum,PVi,PGi,Y,UAngel,r,c,Angle,Loadi);
     %%  开始构建ddg
-    ddg=func_ddg(PGi,PVi,Busnum,RestraintCount);
+    ddg=func_ddg(PGi,PVi,Busnum,RestraintCount,Loadi);
     %%  开始构建deltF
-    deltF=func_deltF(PG,PVi,GenC,ContrlCount,PGi,wG,wD,PG0,PD0,PD,Busnum);
+    deltF=func_deltF(PG,PVi,GenC,ContrlCount,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);
+    Mat_G=FormG(Volt,PVi,PGi,PG,QG,PD,Loadi);
     Mat_H=FormH(Busnum,GB,AngleIJMat,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);
-    Lw=FormLw(Mat_G,Init_U,GenU,Busnum,PVQU,PD0);
+    Lz=FormLz(Mat_G,Init_L,GenL,Busnum,PVQL,PD0,randPDind,Loadi);
+    Lw=FormLw(Mat_G,Init_U,GenU,Busnum,PVQU,PD0,randPDind,Loadi);
     Lx=FormLx(deltF,deltH,Init_Y,deltG,Init_Z,Init_W);
     YY=FormYY1(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);
+    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]=AssignXX1(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);
+    [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)
+ObjectiveFun(PG,PG0,PGi,PD,PD0,wG,wD,Loadi)
 DrawGap(plotGap);
-Volt=Volt';
+%% 
+Volt=full(Volt');
+PD=full(PD);
 toc  
 
diff --git a/OPF_Init.asv b/OPF_Init.asv
index 9461ddf..6cff26a 100644
--- a/OPF_Init.asv
+++ b/OPF_Init.asv
@@ -1,27 +1,29 @@
-function [Volt,UAngel,Init_Z,Init_W,Init_L,Init_U,Init_Y,PG,QG,RestraintCount,wD,PD]=OPF_Init(Busnum,Balance,PG,QG,Volt,GenU,GenL,PVi,PGi,PVQU,PVQL,PD0)
+function [Volt,UAngel,Init_Z,Init_W,Init_L,Init_U,Init_Y,PG,QG,RestraintCount,wG,wD,PD,PD0,randPDind]=OPF_Init(Busnum,Balance,PG,QG,Volt,GenU,GenL,PVi,PGi,PVQU,PVQL,PD0)
 RestraintCount=size(PVi,1)+size(PGi,1)+Busnum*2; %约束条件数
 t_Bal_volt=Volt(Balance);
-Volt=ones(1,Busnum);
+Volt=sparse(ones(1,Busnum));
 Volt(Balance)=t_Bal_volt;
-%Volt(Balance)=1;
-UAngel=zeros(1,Busnum);
-Init_Z=ones(1,RestraintCount);
-Init_W=-1*ones(1,RestraintCount);
-Init_L=ones(1,RestraintCount);
-Init_U=ones(1,RestraintCount);
-%Init_Y=zeros(1,2*Busnum);
-%Init_Y=ones(1,2*Busnum);
-Init_Y=zeros(1,2*Busnum);%与学姐一致
-%Init_Y(1:2:2*Busnum)=1e-10;
-%Init_Y(2:2:2*Busnum)=-1e-10;
-tPU=GenU(:,2);% 发电机有功上限
-tQU=PVQU(:,1);% 无功上限
-tPL=GenL(:,2);% 发电机有功下限
-tQL=PVQL(:,1);% 无功下限
-%PG(4:5)=[4.5 4.5];
+UAngel=sparse(1,Busnum);
+Init_Z=sparse(ones(1,RestraintCount));
+Init_W=sparse(-1*ones(1,RestraintCount));
+Init_L=sparse(ones(1,RestraintCount));
+Init_U=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));% 无功下限
 PG(PGi)=(tPU+tPL)/2;
-%QG(4:5)=[0 1.45];
 QG(PVi)=(tQU+tQL)/2;
+wG=ones(size(PGi,1),1);
+randInt=randperm(300);
+asc_randInt=sort(randInt,2,'ascend' );
+randPDind=randInt(1:200);
 wD=ones(Busnum,1);
-PD=.5*PD0
+wD(randPDind)=0;%一些负荷不约束
+%wD(Balance)=0;
+PD0(randPDind)=0;
+PD=1*PD0;
+%PD(PD==0)=.2;
+
 end
\ No newline at end of file
diff --git a/OPF_Init.m b/OPF_Init.m
index 28eaa75..7f18f70 100644
--- a/OPF_Init.m
+++ b/OPF_Init.m
@@ -1,28 +1,33 @@
-function [Volt,UAngel,Init_Z,Init_W,Init_L,Init_U,Init_Y,PG,QG,RestraintCount,wG,wD,PD]=OPF_Init(Busnum,Balance,PG,QG,Volt,GenU,GenL,PVi,PGi,PVQU,PVQL,PD0)
-RestraintCount=size(PVi,1)+size(PGi,1)+Busnum*2; %约束条件数
+function [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)
+Loadi=find(QD~=0);
+%Loadi=[1:Busnum]';
+RestraintCount=size(PVi,1)+size(PGi,1)+size(Loadi,1)+Busnum*1; %约束条件数
 t_Bal_volt=Volt(Balance);
 Volt=sparse(ones(1,Busnum));
 Volt(Balance)=t_Bal_volt;
-%Volt(Balance)=1;
 UAngel=sparse(1,Busnum);
 Init_Z=sparse(ones(1,RestraintCount));
 Init_W=sparse(-1*ones(1,RestraintCount));
 Init_L=sparse(ones(1,RestraintCount));
 Init_U=sparse(ones(1,RestraintCount));
 Init_Y=sparse(1,2*Busnum);%与学姐一致
-%Init_Y(1:2:2*Busnum)=1e-10;
-%Init_Y(2:2:2*Busnum)=-1e-10;
 tPU=sparse(GenU(:,2));% 发电机有功上限
 tQU=sparse(PVQU(:,1));% 无功上限
 tPL=sparse(GenL(:,2));% 发电机有功下限
 tQL=sparse(PVQL(:,1));% 无功下限
-%PG(4:5)=[4.5 4.5];
 PG(PGi)=(tPU+tPL)/2;
-%QG(4:5)=[0 1.45];
 QG(PVi)=(tQU+tQL)/2;
 wG=ones(size(PGi,1),1);
-wD=ones(Busnum,1);
+randInt=randperm(size(Loadi,1));
+%asc_randInt=sort(randInt,2,'ascend' );
+%randPDind=randInt(1:200);
+
+randPDind=randInt(1:0);
+
+wD=ones(size(Loadi,1),1);
+wD(randPDind)=0;%一些负荷不约束
 %wD(Balance)=0;
+%PD0(randPDind)=0;
 PD=1*PD0;
 %PD(PD==0)=.2;
 
diff --git a/ObjectiveFun.m b/ObjectiveFun.m
index d1a21d6..967e682 100644
--- a/ObjectiveFun.m
+++ b/ObjectiveFun.m
@@ -1,6 +1,6 @@
-function [out_arg]=ObjectiveFun(PG,PG0,PGi,PD,PD0,wG,wD)
+function [out_arg]=ObjectiveFun(PG,PG0,PGi,PD,PD0,wG,wD,Loadi)
 t1=PG(PGi)-PG0(PGi);
 t2=wG.*t1.*t1;
-t3=wD.*((PD-PD0).*(PD-PD0));
+t3=wD.*((PD(Loadi)-PD0(Loadi)).*(PD(Loadi)-PD0(Loadi)));
 out_arg= sum(t2)+sum(t3);
 end
\ No newline at end of file
diff --git a/SolveIt.m b/SolveIt.m
index 0c61fd1..b1d72e7 100644
--- a/SolveIt.m
+++ b/SolveIt.m
@@ -1,4 +1,4 @@
-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)
+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);
@@ -9,13 +9,13 @@ aa=[
     ];
 yy=[LxComa;-Ly];
 %% 平衡节点电压不变 
-t=size(PVi,1)+size(PGi,1)+Busnum;
+t=size(PVi,1)+size(PGi,1)+size(Loadi,1);
 aa(t+Balance,:)=0;
 aa(:,t+Balance)=0;
 aa(t+Balance,t+Balance)=1;
 deltG(t+Balance,:)=0;
 %%
-t=size(PVi,1)+size(PGi,1)+Busnum*2;
+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;
diff --git a/func_ddg.m b/func_ddg.m
index bc47135..1dbc89d 100644
--- a/func_ddg.m
+++ b/func_ddg.m
@@ -1,5 +1,5 @@
-function ddg=func_ddg(PGi,PVi,Busnum,RestraintCount)
+function ddg=func_ddg(PGi,PVi,Busnum,RestraintCount,Loadi)
 
-ddg=sparse(size(PVi,1)+size(PGi,1)+2*Busnum,RestraintCount);
+ddg=sparse(size(PVi,1)+size(PGi,1)+size(Loadi,1)+1*Busnum,RestraintCount);
 
 end
\ No newline at end of file
diff --git a/func_ddh3.m b/func_ddh3.m
index a9b01ec..0afb739 100644
--- a/func_ddh3.m
+++ b/func_ddh3.m
@@ -1,6 +1,6 @@
-function ddh=func_ddh3(AngleIJMat,GB,Volt,Init_Y,Busnum,PVi,PGi,Y,UAngel,r,c,Angle)
+function ddh=func_ddh3(AngleIJMat,GB,Volt,Init_Y,Busnum,PVi,PGi,Y,UAngel,r,c,Angle,Loadi)
 %决定用循环重写
-ContrlCount=size(PVi,1)+size(PGi,1)+Busnum*3;
+ContrlCount=size(PVi,1)+size(PGi,1)+size(Loadi,1)+Busnum*2;
 
 %AngleIJ=AngleIJMat-angle(GB);
 mat_AngleIJ=sparse(r,c,UAngel(r)-UAngel(c)-Angle',Busnum,Busnum);
@@ -47,8 +47,9 @@ t=[ddPdVdV+ddQdVdV,ddPdTdV+ddQdTdV ;
 ];
 sizePGi=size(PGi,1);
 sizePVi=size(PVi,1);
+sizeLoadi=size(Loadi,1);
 ddh=[
-    sparse(sizePGi+sizePVi+Busnum,ContrlCount);
-    sparse(2*Busnum,sizePVi+sizePGi+Busnum),-t;
+    sparse(sizePGi+sizePVi+sizeLoadi,ContrlCount);
+    sparse(2*Busnum,sizePVi+sizePGi+sizeLoadi),-t;
     ];
 end
\ No newline at end of file
diff --git a/func_deltF.m b/func_deltF.m
index 1de0334..db6a5a2 100644
--- a/func_deltF.m
+++ b/func_deltF.m
@@ -1,4 +1,4 @@
-function deltF=func_deltF(PG,PVi,GenC,ContrlCount,PGi,wG,wD,PG0,PD0,PD,Busnum)
+function deltF=func_deltF(PG,PVi,GenC,ContrlCount,PGi,wG,wD,PG0,PD0,PD,Busnum,Loadi)
 %t1=PG(setdiff(PVi,Balance));
 % t2=Volt'*Volt;
 % t3=real(GB).*cos(AngleIJMat)+imag(GB).*sin(AngleIJMat);
@@ -8,7 +8,7 @@ function deltF=func_deltF(PG,PVi,GenC,ContrlCount,PGi,wG,wD,PG0,PD0,PD,Busnum)
 % PPG=([PQ(1),PBal])';%暂时用土办法处理一下
 %% 
 t1=2*wG.*(PG(PGi) - PG0(PGi) );
-t2=2*wD.*(PD-PD0);
+t2=2*wD.*(PD(Loadi)-PD0(Loadi));
 deltF=[
     sparse(t1);
     sparse(size(PVi,1),size(PVi,2));
diff --git a/func_deltG.m b/func_deltG.m
index 56ce9e4..26d7ca0 100644
--- a/func_deltG.m
+++ b/func_deltG.m
@@ -1,46 +1,27 @@
-function deltG=func_deltG(Busnum,PVi,PGi)
+function deltG=func_deltG(Busnum,PVi,PGi,Loadi)
 %%
 sizePGi=size(PGi,1);
 sizePVi=size(PVi,1);
+sizeLoadi=size(Loadi,1);
 %%
-% dg1_dPg=eye(size(PGi,1));
-% dg2_dPg=zeros(size(PGi,1),size(PVi,1));
-% dg3_dPg=zeros(size(PGi,1),Busnum);
-% dg4_dPg=zeros(size(PGi,1),Busnum);
 dg1_dPg=sparse(1:sizePGi,1:sizePGi,ones(sizePGi,1),sizePGi,sizePGi);
 dg2_dPg=sparse(sizePGi,sizePVi);
-dg3_dPg=sparse(sizePGi,Busnum);
+dg3_dPg=sparse(sizePGi,sizeLoadi);
 dg4_dPg=sparse(sizePGi,Busnum);
 %% 
-% dg1_dQr=zeros(size(PVi,1),size(PGi,1));
-% dg2_dQr=eye(size(PVi,1));
-% dg3_dQr=zeros(size(PVi,1),Busnum);
-% dg4_dQr=zeros(size(PVi,1),Busnum);
 dg1_dQr=sparse(sizePVi,sizePGi);
 dg2_dQr=sparse(1:sizePVi,1:sizePVi,ones(sizePVi,1),sizePVi,sizePVi);
-dg3_dQr=sparse(sizePVi,Busnum);
+dg3_dQr=sparse(sizePVi,sizeLoadi);
 dg4_dQr=sparse(sizePVi,Busnum);
 %%
-% dg1_dPD=zeros(Busnum,size(PGi,1));
-% dg2_dPD=zeros(Busnum,size(PVi,1));
-% dg3_dPD=eye(Busnum,Busnum);
-% dg4_dPD=zeros(Busnum,Busnum);
-dg1_dPD=sparse(Busnum,size(PGi,1));
-dg2_dPD=sparse(Busnum,size(PVi,1));
-dg3_dPD=sparse(1:Busnum,1:Busnum,ones(Busnum,1),Busnum,Busnum);
-dg4_dPD=sparse(Busnum,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),Busnum);
 %% 
-% dg1_dx=zeros(2*Busnum,size(PGi,1));
-% dg2_dx=zeros(2*Busnum,size(PVi,1));
-% dg3_dx=zeros(2*Busnum,Busnum);
-% dg4_dx=zeros(2*Busnum,Busnum);
-% for I=1:Busnum
-%     %dg3_dx(2*I,I)=1;暂时改一下
-%     dg4_dx(I,I)=1;
-% end
 dg1_dx=sparse(2*Busnum,sizePGi);
 dg2_dx=sparse(2*Busnum,sizePVi);
-dg3_dx=sparse(2*Busnum,Busnum);
+dg3_dx=sparse(2*Busnum,sizeLoadi);
 dg4_dx=[sparse(1:Busnum,1:Busnum,ones(Busnum,1),Busnum,Busnum);
     sparse(Busnum,Busnum);
     ];
diff --git a/func_deltH.m b/func_deltH.m
index 1c1ea2a..871346d 100644
--- a/func_deltH.m
+++ b/func_deltH.m
@@ -1,24 +1,8 @@
-function deltH=func_deltH(Busnum,Volt,PVi,AngleIJMat,Y,GB,PGi,UAngel,r,c,Angle)
-% dH_dPg=zeros(size(PGi,1),2*Busnum);
-% 
-% for I=1:size(PGi,1)
-%     %dH_dPg(I,2*PVi(I)-1)=-1;王锡凡书上的公式
-%     %dH_dPg(I,2*PGi(I)-1)=1;暂时改一下20111227
-%     dH_dPg(I,PGi(I))=1;
-% end
+function deltH=func_deltH(Busnum,Volt,PVi,AngleIJMat,Y,GB,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=zeros(size(PVi,1),2*Busnum);
-% for I=1:size(PVi,1)
-%     %dH_dQr(I,2*PVi(I))=-1;王锡凡书上的公式
-%     %dH_dQr(I,2*PVi(I))=1;暂时改一下20111227
-%     dH_dQr(I,PVi(I)+Busnum)=1;
-% end
 dH_dQr=sparse(1:size(PVi,1),PVi+Busnum,ones(size(PVi,1),1),size(PVi,1),2*Busnum);
-dH_dPD=[sparse(1:Busnum,1:Busnum,-ones(Busnum,1),Busnum,Busnum) sparse(Busnum,Busnum)];
-%Angle=angle(GB);
-dH_dx = jacobian_M3(Busnum,Volt,Y,Angle,AngleIJMat,UAngel,r,c);  %形成雅克比矩阵
-%deltH=[dH_dPg;dH_dQr;dH_dx'];%dH_dx 需要使用一下转置 暂时改一下
+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
\ No newline at end of file
diff --git a/func_deltdeltF.m b/func_deltdeltF.m
index 935d370..db1db61 100644
--- a/func_deltdeltF.m
+++ b/func_deltdeltF.m
@@ -1,6 +1,6 @@
 function deltdeltF=func_deltdeltF(Busnum,GenC,PVi,PGi,wG,wD,PD0)
 
-ContrlCount=size(PVi,1)+size(PGi,1)+Busnum*3; %P,Q,Volt theta这些控制变量数
+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);
diff --git a/jacobian_M3.m b/jacobian_M3.m
index 00053a9..c224cf2 100644
--- a/jacobian_M3.m
+++ b/jacobian_M3.m
@@ -1,4 +1,4 @@
-function [Jacob]=jacobian_M3(Busnum,Volt,Y,Angle,AngleIJMat,UAngel,r,c)   
+function [Jacob]=jacobian_M3(Busnum,Volt,Y,Angle,UAngel,r,c)   
 %**************************************************************************
 %  程序功能 : 子函数——形成雅可比矩阵Jacobian
 %  编    者：