diff --git a/FormH.m b/FormH.m
index 87a5378..5b735eb 100644
--- a/FormH.m
+++ b/FormH.m
@@ -1,8 +1,8 @@
 function Mat_H=FormH(Busnum,Volt,PG,PD,QG,QD,Y,UAngel,r,c,Angle,QD_NON_ZERO,QD_NON_ZERO_IND)
 %% 
-QDcos=textread('1047glys.txt');
-QD(QD~=0)=PD(QD~=0)./tan(QDcos);
-QD(QD_NON_ZERO_IND)=QD_NON_ZERO;
+% QDcos=textread('1047glys.txt');
+% 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';
diff --git a/FormLw.asv b/FormLw.asv
index 150c0a5..1c9e4b4 100644
--- a/FormLw.asv
+++ b/FormLw.asv
@@ -1,13 +1,14 @@
-function Lw=FormLw(Mat_G,Init_U,GenU,Busnum,PVQU,PD0)
+function Lw=FormLw(Mat_G,Init_U,GenU,Busnum,PVQU,PD0,Loadi,KK)
 
 PU=1*GenU(:,2);%发电机有功上界
 QU=PVQU(:,1);%发电机无功上界
-VoltU=1.2*ones(1,Busnum);
-PDU=PD0;
-PDU(PD0>0)=1.2*PD0()';
-PDU(PDU==0)=20;
-PDU=25*ones(Busnum,1)';
-t1=([PU',QU',PDU,VoltU])';
+VoltU=(1.2+1/ekk)*ones(1,Busnum);
+PDU=PD0(Loadi);
+PDU(PDU>0)=1.2*PDU(PDU>0);
+PDU(PDU<0)=0.8*PDU(PDU<0);
+PDU(PDU==0)=.2;
+t1=([PU',QU',PDU',VoltU])';
 t2=Mat_G+Init_U'-t1;
 Lw=t2;
+
 end
\ No newline at end of file
diff --git a/FormLw.m b/FormLw.m
index 2029450..37351ca 100644
--- a/FormLw.m
+++ b/FormLw.m
@@ -1,12 +1,12 @@
-function Lw=FormLw(Mat_G,Init_U,GenU,Busnum,PVQU,PD0,Loadi)
-
+function Lw=FormLw(Mat_G,Init_U,GenU,Busnum,PVQU,PD0,Loadi,KK)
+%KK=999;
 PU=1*GenU(:,2);%发电机有功上界
 QU=PVQU(:,1);%发电机无功上界
-VoltU=1.1*ones(1,Busnum);
+VoltU=(1.2+1/exp(KK))*ones(1,Busnum);
 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)=1.20*PDU(PDU>0);
+PDU(PDU<0)=.80*PDU(PDU<0);
+PDU(PDU==0)=.20;
 t1=([PU',QU',PDU',VoltU])';
 t2=Mat_G+Init_U'-t1;
 Lw=t2;
diff --git a/FormLz.asv b/FormLz.asv
index b4e5e7c..d73c882 100644
--- a/FormLz.asv
+++ b/FormLz.asv
@@ -1,28 +1,14 @@
-function Lz=FormLz(Mat_G,Init_L,GenL,Busnum,PVQL,PD0,randPDind,Loadi)
-
+function Lz=FormLz(Mat_G,Init_L,GenL,Busnum,PVQL,PD0,Loadi,KK)
+%KK=999;
 PL=1*GenL(:,2);%发电机有功下界
 QL=PVQL(:,1);%发电机无功下界
-VoltL=0.9*ones(1,Busnum);
+VoltL=(0.8-1/exp(KK))*ones(1,Busnum);
 PDL=PD0(Loadi);
-PDL(PD0>0)=0.8*PD0(PDL>0);
-PDL(PD0==0)=-.2;
-PDL(PDL<0)=1.2*PD0(PDL<0);
+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;
 
-% 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 31263fb..eedd55c 100644
--- a/FormLz.m
+++ b/FormLz.m
@@ -1,12 +1,12 @@
-function Lz=FormLz(Mat_G,Init_L,GenL,Busnum,PVQL,PD0,Loadi)
-
+function Lz=FormLz(Mat_G,Init_L,GenL,Busnum,PVQL,PD0,Loadi,KK)
+%KK=999;
 PL=1*GenL(:,2);%发电机有功下界
 QL=PVQL(:,1);%发电机无功下界
-VoltL=0.9*ones(1,Busnum);
+VoltL=(0.8-1/exp(KK))*ones(1,Busnum);
 PDL=PD0(Loadi);
-PDL(PDL>0)=0.8*PDL(PDL>0);
-PDL(PDL<0)=1.2*PDL(PDL<0);
-PDL(PDL==0)=-.2;
+PDL(PDL>0)=.80*PDL(PDL>0);
+PDL(PDL<0)=1.20*PDL(PDL<0);
+PDL(PDL==0)=-.20;
 t1=([PL',QL',PDL',VoltL])';
 t2=Mat_G-Init_L'-t1;
 Lz=t2;
diff --git a/OPF.asv b/OPF.asv
index 9e397ab..a5d5ebb 100644
--- a/OPF.asv
+++ b/OPF.asv
@@ -1,6 +1,6 @@
 tic
 clear
-[kmax,Precision,UAngel,Volt,Busnum,PVi,PVu,Balance,Y,Angle,P0,Q0,r,c,GB,Linei,Linej,Transfori,Transforj,GenU,GenL,GenC,PG,QG,PD,QD,CenterA,PGi,PVQU,PVQL]=pf('ieee10471PG.dat');
+[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(c:/iee.txt);
 %PVi电压节点序号
 %PVu电压节点电压标幺值
 Volt;
@@ -34,28 +34,28 @@ while(abs(Gap)>Precision)
     L_1Z=diag(Init_Z./Init_L);
     U_1W=diag(Init_W./Init_U);
     %% 形成海森阵
-    deltdeltF=func_deltdeltF(Busnum,PVi,PGi,wG,wD,ContrlCount);
+    deltdeltF=func_deltdeltF(PVi,wG,wD,ContrlCount);
     %% 形成ddHy
-    ddh=func_ddh3(AngleIJMat,GB,Volt,Init_Y,Busnum,PVi,PGi,Y,UAngel,r,c,Angle,Loadi);
+    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,GenC,ContrlCount,PGi,wG,wD,PG0,PD0,PD,Busnum,Loadi);
+    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,GB,AngleIJMat,Volt,PG,PD,QG,QD,Y,UAngel,r,c,Angle,QD_NON_ZERO,QD_NON_ZERO_IND);
+    Mat_H=FormH(Busnum,Volt,PG,PD,QG,QD,Y,UAngel,r,c,Angle,QD_NON_ZERO,QD_NON_ZERO_IND);
     Ly=Mat_H;
-    Lz=FormLz(Mat_G,Init_L,GenL,Busnum,PVQL,PD0,randPDind,Loadi);
-    Lw=FormLw(Mat_G,Init_U,GenU,Busnum,PVQU,PD0,randPDind,Loadi);
+    Lz=FormLz(Mat_G,Init_L,GenL,Busnum,PVQL,PD0,Loadi);
+    Lw=FormLw(Mat_G,Init_U,GenU,Busnum,PVQU,PD0,Loadi);
     Lx=FormLx(deltF,deltH,Init_Y,deltG,Init_Z,Init_W);
-    YY=FormYY1(Lul,Lz,Ly,Luu,Lw,Lx);
+    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]=AssignXX1(XX,ContrlCount,RestraintCount,Busnum);
+    [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;
@@ -64,7 +64,7 @@ fprintf('
 ObjectiveFun(PG,PG0,PGi,PD,PD0,wG,wD,Loadi)
 DrawGap(plotGap);
 %% 
-Volt=full(Volt');
-PD=full(PD);
+% Volt=full(Volt');
+% PD=full(PD);
 toc  
 
diff --git a/OPF.m b/OPF.m
index cb32244..18613dc 100644
--- a/OPF.m
+++ b/OPF.m
@@ -1,6 +1,6 @@
 tic
 clear
-[kmax,Precision,UAngel,Volt,Busnum,PVi,PVu,Balance,Y,Angle,P0,Q0,r,c,GB,Linei,Linej,Transfori,Transforj,GenU,GenL,GenC,PG,QG,PD,QD,CenterA,PGi,PVQU,PVQL]=pf('ieee10471PG.dat');
+[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('e:/算例/standard.txt');
 %PVi电压节点序号
 %PVu电压节点电压标幺值
 Volt;
@@ -13,7 +13,7 @@ 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=600;
+kmax=60;
 %% 20120523 临时
 QD_NON_ZERO=QD(PD==0 & QD~=0);
 QD_NON_ZERO_IND=find(PD==0 & QD~=0);
@@ -48,8 +48,8 @@ while(abs(Gap)>Precision)
     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);
-    Lw=FormLw(Mat_G,Init_U,GenU,Busnum,PVQU,PD0,Loadi);
+    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);
     %% 开始解方程
@@ -64,7 +64,7 @@ fprintf('
 ObjectiveFun(PG,PG0,PGi,PD,PD0,wG,wD,Loadi)
 DrawGap(plotGap);
 %% 
-Volt=full(Volt');
-PD=full(PD);
+% Volt=full(Volt');
+% PD=full(PD);
 toc  
 
diff --git a/OPF_Init.m b/OPF_Init.m
index f0c4745..1f04e02 100644
--- a/OPF_Init.m
+++ b/OPF_Init.m
@@ -3,13 +3,13 @@ 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=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=sparse(ones(1,RestraintCount));
-Init_U=sparse(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));% 无功上限
@@ -19,10 +19,11 @@ PG(PGi)=(tPU+tPL)/2;
 QG(PVi)=(tQU+tQL)/2;
 wG=ones(size(PGi,1),1);
 randInt=randperm(size(Loadi,1));
-randPDind=randInt(1:262);
+randPDind=randInt(1:0);
 
 wD=ones(size(Loadi,1),1);
 wD(randPDind)=0;%一些负荷不约束
+%wD(12)=0;
 PD=1*PD0;
 
 end
\ No newline at end of file
diff --git a/SolveIt.asv b/SolveIt.asv
index 1558133..dceba29 100644
--- a/SolveIt.asv
+++ b/SolveIt.asv
@@ -1,25 +1,31 @@
-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);
 t2=-deltG*( t1 )*deltG';
 aa=[
     (H+t2),deltH;
-    deltH',zeros(size(Init_Y,),ContrlCount);
+    deltH',zeros(size(Init_Y,2));
     ];
 yy=[LxComa;-Ly];
-% t=size(PVi,1)+size(PGi,1);
-% aa(t+2*Balance-1,:)=0;
-% aa(:,t+2*Balance-1)=0;
-% aa(t+2*Balance-1,t+2*Balance-1)=1;
-%暂时改一下
-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;
+%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:RestraintCount+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;
diff --git a/SolveIt.m b/SolveIt.m
index b1d72e7..672f754 100644
--- a/SolveIt.m
+++ b/SolveIt.m
@@ -12,13 +12,15 @@ 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(t+Balance,t+Balance)=1;
+aa=aa+sparse(t+Balance,t+Balance,ones(length(Balance),1),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(t+Balance,t+Balance)=1;
+aa=aa+sparse(t+Balance,t+Balance,ones(length(Balance),1),ContrlCount+2*Busnum,ContrlCount+2*Busnum);
 deltG(t+Balance,:)=0;
 %%
 dxdy=aa\yy;
diff --git a/jacobian.m b/jacobian.m
index 4c8d657..c7ecfb9 100644
--- a/jacobian.m
+++ b/jacobian.m
@@ -22,7 +22,8 @@ P = P0+temp3';
 %% 处理平衡节点和pv节点
 H(:,Balance) = 0;                                                     
 H(Balance,:) = 0;
-H(Balance,Balance) = 100;                                                           % 平衡节点对应的对角元素置一个有限数
+%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
diff --git a/openfile.m b/openfile.m
index 7545b7c..154d8ff 100644
--- a/openfile.m
+++ b/openfile.m
@@ -6,16 +6,18 @@ function [Busnum,Balance,PQstandard,Precision,Linei,Linej,Liner,Linex,Lineb,kmax
 %    编制时间 ：2010.12
 %**************************************************************************
 data = dlmread(FileName);                                           % 一次读入全部数据
+zeroRow = find(data(:,1)==0); 
 Busnum= data(1,1);                                                       % 节点数
 PQstandard = data(1,3);                                               % 基准容量
 kmax = data(1,4);                                                           %最大迭代次数
 Precision = data(2,1);                                                    % 精度
-Balance = data(3,2);                                                      % 生成1到节点号的列向量
+%Balance = data(3,2);     
+Balance=data(3:zeroRow(1)-1,2);% 生成1到节点号的列向量
 CenterA=data(1,5);  %中心参数
 LineNum=data(1,2); %支路数
 Base=data(1,3);
 %% 各参数矩阵分块
-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,:);              % 形成变压器参数矩阵