diff --git a/FormLw.m b/FormLw.m
index f403971..ecd112c 100644
--- a/FormLw.m
+++ b/FormLw.m
@@ -1,5 +1,5 @@
 function Lw=FormLw(Loadi,Mat_G,Init_U)
 
-upper=0.015*sparse(ones(length(Loadi)*1,1));
+upper=5*sparse(ones(length(Loadi)*1,1));
 Lw=Mat_G+Init_U-upper;
 end
\ No newline at end of file
diff --git a/IPMLoop.m b/IPMLoop.m
index aed3020..c85ab45 100644
--- a/IPMLoop.m
+++ b/IPMLoop.m
@@ -2,8 +2,8 @@ function [ V1r,V1i,I1r,I1i ] = IPMLoop(measurement,BalI1r,BalI1i,busNum,Loadi,fs
 %把每个序的循环写在这个函数中。其实也就是内点法循环。
 V1r=1*ones(busNum,1);
 V1i=0*ones(busNum,1);
-I1r=0.1*ones(length(Loadi),1);%注入电流，相当于放电机电流，与负荷电流负荷相反
-I1i=0.1*ones(length(Loadi),1);%注入电流，相当于放电机电流，与负荷电流负荷相反
+I1r=0.0*ones(length(Loadi),1);%注入电流，相当于放电机电流，与负荷电流负荷相反
+I1i=0.0*ones(length(Loadi),1);%注入电流，相当于放电机电流，与负荷电流负荷相反
 KK=0;
 plotGap=zeros(1,60);
 %初始化
@@ -34,6 +34,7 @@ while(abs(Gap)>0.000001)
     U_1W=diag(Init_W./Init_U);
     %% 形成海森阵
     deltdeltF=func_deltdeltF(busNum,fsY1,Loadi);
+%     deltdeltF=0;
     %% 形成ddHy
 %     ddh=func_ddh(busNum,Loadi,Init_Z,Init_W);
     ddh=0;
@@ -41,6 +42,7 @@ while(abs(Gap)>0.000001)
     ddg=func_ddg(busNum,Loadi,Init_Z,Init_W);
     %%  开始构建deltF
     deltF=func_deltF(measurement,busNum,fsY1,Loadi,V1r,V1i,I1r,I1i);
+%     deltF=0;
     %%
     Luu=Init_U.*Init_W+Init_u*ones(RestraintCount,1);
     Lul=Init_L.*Init_Z-Init_u*ones(RestraintCount,1);
@@ -55,7 +57,7 @@ while(abs(Gap)>0.000001)
     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,busNum);
     %%取各分量
     [deltZ,deltL,deltW,deltU,deltX,deltY]=AssignXX(XX,ContrlCount,RestraintCount,busNum);
-    [Init_Z,Init_L,Init_W,Init_U,Init_Y,V1r,V1i,I1r,I1i]=Modification(Init_Z,Init_L,Init_W,Init_U,Init_Y,deltZ,deltL,deltW,deltU,deltX,deltY,V1r,V1i,I1r,I1i,ContrlCount,Balance,busNum,Loadi);
+    [Init_Z,Init_L,Init_W,Init_U,Init_Y,V1r,V1i,I1r,I1i]=Modification(Init_Z,Init_L,Init_W,Init_U,Init_Y,deltZ,deltL,deltW,deltU,deltX,deltY,V1r,V1i,I1r,I1i,ContrlCount,Balance,busNum,Loadi,Vref);
     Gap=(Init_L'*Init_Z-Init_U'*Init_W);
     fprintf('Gap %f\n',full(Gap));
     KK=KK+1;
diff --git a/Modification.m b/Modification.m
index a1a94e5..6c0afcf 100644
--- a/Modification.m
+++ b/Modification.m
@@ -1,4 +1,4 @@
-function [Init_Z,Init_L,Init_W,Init_U,Init_Y,V1r,V1i,I1r,I1i]=Modification(Init_Z,Init_L,Init_W,Init_U,Init_Y,deltZ,deltL,deltW,deltU,deltX,deltY,V1r,V1i,I1r,I1i,ContrlCount,Balance,Busnum,Loadi)
+function [Init_Z,Init_L,Init_W,Init_U,Init_Y,V1r,V1i,I1r,I1i]=Modification(Init_Z,Init_L,Init_W,Init_U,Init_Y,deltZ,deltL,deltW,deltU,deltX,deltY,V1r,V1i,I1r,I1i,ContrlCount,Balance,Busnum,Loadi,Vref)
 AlphaP=FormAlphaP(Init_L,deltL,Init_U,deltU);
 % fprintf('AlphaP %f\n',full(AlphaP));
 AlphaD=FormAlphaD(Init_Z,deltZ,Init_W,deltW);
@@ -12,7 +12,7 @@ Init_Y=Init_Y+AlphaD*deltY;
 %PG(PGi)=PG(PGi)+deltX(1:size(PGi,1));
 % PG(PGi)=PG(PGi)+AlphaP*deltX(1:size(PGi,1));
 V1r=V1r+AlphaP*deltX(1:Busnum);
-V1r(Balance)=1;
+V1r(Balance)=Vref;
 V1i=V1i+AlphaP*deltX(Busnum+1:2*Busnum);
 V1i(Balance)=0;
 %QG(PVi)=QG(PVi)+deltX(size(PGi,1)+1:size(PVi,1)+size(PGi,1) );
diff --git a/run.m b/run.m
index 9bff95f..165ce92 100644
--- a/run.m
+++ b/run.m
@@ -232,19 +232,30 @@ mIf2=-If2;
 
 %% 先算正序的
 %平衡节点电流
+fprintf('正序\n');
 BalI1r=real(-sum(mIf1));
 BalI1i=imag(-sum(mIf1));
 % inv(fsY11)*(mIf1);
 measurement=-mIf1(Loadi);
-[ V1r,V1i,I1r,I1i ]=IPMLoop(measurement,BalI1r,BalI1i,busNum,Loadi,fsY1,Balance,1 );
+[ V1r,V1i,I1r,I1i ]=IPMLoop(measurement,BalI1r,BalI1i,busNum,Loadi,fsY11,Balance,1 );
 f=sum(([real(measurement);imag(measurement)]-[-I1r;-I1i]).^2);
-
+fprintf('目标值 %f\n',full(f));
 %% 算负序的
+fprintf('负序\n');
 BalI2r=real(-sum(mIf2));
 BalI2i=imag(-sum(mIf2));
 measurement=-mIf2(Loadi);
-[ V2r,V2i,I2r,I2i ]=IPMLoop(measurement,BalI2r,BalI2i,busNum,Loadi,fsY2,Balance,1 );
+[ V2r,V2i,I2r,I2i ]=IPMLoop(measurement,BalI2r,BalI2i,busNum,Loadi,fsY22,Balance,0 );
 f=sum([real(measurement);imag(measurement)]-[-I2r;-I2i]);
+fprintf('目标值 %f\n',full(f));
+%% 算零序
+fprintf('零序\n');
+BalI0r=real(-sum(mIf0));
+BalI0i=imag(-sum(mIf0));
+measurement=-mIf0(Loadi);
+[ V0r,V0i,I0r,I0i ]=IPMLoop(measurement,BalI0r,BalI0i,busNum,Loadi,fsY00,Balance,0 );
+f=sum([real(measurement);imag(measurement)]-[-I0r;-I0i]);
+fprintf('目标值 %f\n',full(f));
 %状态量
 % SEVoltpA=sparse(ones(busNum,1));
 % SEVoltpB=sparse(ones(busNum,1)).*exp(1j*-120/180*pi);