diff --git a/IPMLoop.m b/IPMLoop.m
index c3654e5..ec3d661 100644
--- a/IPMLoop.m
+++ b/IPMLoop.m
@@ -1,4 +1,4 @@
-function [ V1r,V1i,I1r,I1i ] = IPMLoop(measurement,BalI1r,BalI1i,busNum,Loadi,fsY1,Balance )
+function [ V1r,V1i,I1r,I1i ] = IPMLoop(measurement,BalI1r,BalI1i,busNum,Loadi,fsY1,Balance,Vref )
 %把每个序的循环写在这个函数中。其实也就是内点法循环。
 V1r=1*ones(busNum,1);
 V1i=0*ones(busNum,1);
@@ -29,13 +29,13 @@ while Gap>1e-5 && KK<20
     deltH=func_deltH(busNum,fsY1,Loadi,Balance);
     deltG=func_deltG(busNum,Loadi,I1r,I1i);
 %     deltG=0;
-    L_1Z=diag(Init_Z./Init_L);
-    U_1W=diag(Init_W./Init_U);
+%     L_1Z=diag(Init_Z./Init_L);
+%     U_1W=diag(Init_W./Init_U);
     deltdeltF=func_deltdeltF(busNum,fsY1,Loadi);
 %     deltdeltF=0;
 %     ddh=func_ddh(SEVmf1,Init_Y,busNum,fsY1amp,SEVaf1,r,c,fsY1ang,Loadi,ContrlCount);
     ddh=0;
-    ddg=func_ddg();
+    ddg=func_ddg(busNum,Loadi,Init_Z,Init_W);
     deltF=func_deltF(measurement,busNum,fsY1,Loadi,V1r,V1i,I1r,I1i);
 %     deltF=0;
     Luu=Init_U.*Init_W+Init_u*ones(RestraintCount,1);
@@ -46,13 +46,13 @@ while Gap>1e-5 && KK<20
     Lz=FormLz(Loadi,Mat_G,Init_L);
     Lw=FormLw(Loadi,Mat_G,Init_U);
     Lx=FormLx(deltF,deltH,Init_Y,deltG,Init_Z,Init_W);
-    YY=FormYY(Lul,Lz,Ly,Luu,Lw,Lx);
+%     YY=FormYY(Lul,Lz,Ly,Luu,Lw,Lx);
     %% 开始解方程
     plotGap(KK)=Gap;
     fprintf('迭代次数 %d Gap %f\n',KK,plotGap(KK));
     XX=SolveIt(deltF,deltG,Init_L,Init_Z,Init_U,Init_W,deltdeltF,ddh,ddg,deltH,Init_Y,Ly,Lz,ContrlCount,Lw,Lul,Luu,Lx,Balance,busNum,Loadi);
     [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);
 end
 
diff --git a/Modification.m b/Modification.m
index e3d4c40..fd47839 100644
--- a/Modification.m
+++ b/Modification.m
@@ -1,5 +1,5 @@
 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)
+    V1r,V1i,I1r,I1i,ContrlCount,Balance,busNum,Loadi,Vref)
 AlphaP=FormAlphaP(Init_L,deltL,Init_U,deltU);
 AlphaD=FormAlphaD(Init_Z,deltZ,Init_W,deltW);
 Init_Z=Init_Z+AlphaD*deltZ;
@@ -9,7 +9,7 @@ Init_U=Init_U+AlphaP*deltU;
 Init_Y=Init_Y+AlphaD*deltY;
 t=deltX(1:busNum*2);
 V1r=V1r+AlphaP*t(1:busNum);
-V1r(Balance)=1;
+V1r(Balance)=Vref;
 V1i=V1i+AlphaP*t(busNum+1:end);
 V1i(Balance)=0;
 t=deltX(busNum*2+1:ContrlCount);
diff --git a/SolveIt.m b/SolveIt.m
index 4e22f8d..a7d573a 100644
--- a/SolveIt.m
+++ b/SolveIt.m
@@ -1,6 +1,6 @@
 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,Lx,Balance,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;
+H=-deltdeltF+ddh+ddg;
 t1=diag(Init_Z./Init_L-Init_W./Init_U);
 t2=-deltG*( t1 )*deltG';
 aa=[
diff --git a/func_ddg.m b/func_ddg.m
index 7406e6d..f79563d 100644
--- a/func_ddg.m
+++ b/func_ddg.m
@@ -1,4 +1,11 @@
-function ddg=func_ddg()
-
-ddg=0;
+function ddg=func_ddg(busNum,Loadi,Init_Z,Init_W)
+% ddg=0;
+Init_ZW=Init_Z+Init_W;
+t=[2*diag( ones(length(Loadi),1).*Init_ZW),zeros(length(Loadi));
+    zeros(length(Loadi)),2*diag( ones(length(Loadi),1).*Init_ZW);
+];
+ddg=[
+    zeros(busNum*2,busNum*2+length(Loadi)*2);
+    zeros(length(Loadi)*2,busNum*2), t;
+];
 end
\ No newline at end of file
diff --git a/run.m b/run.m
index f85eb23..9bff95f 100644
--- a/run.m
+++ b/run.m
@@ -200,6 +200,7 @@ fprintf('
 % VoltpC=Vp3m(3:3:end).*exp(1j*Vp3a(3:3:end));
 %% 用牛顿法求解end
 %% 开始进入状态估计
+clear PD QD PG QG;
 %准备量测量
 iterPhaseASpotLoadP=phaseASpotLoadP;
 iterPhaseBSpotLoadP=phaseBSpotLoadP;
@@ -220,21 +221,30 @@ If2=conj(f012(3,:)');
 %试着算一下正序电流
 %     fsY11*V1;
 %形成负荷序电流的测量值
-mIf0=If0;
+mIf0=-If0;
 mIf1=-If1;
 % mIf1(3)=-mIf1(2);
-mIf2=If2;
+mIf2=-If2;
 % fsY11(:,Balance)=0;
 % fsY11(Balance,:)=0;
 % fsY11=fsY11+sparse(Balance,Balance,ones(length(Balance),1),busNum,busNum);
 % mIf1(3)=1;
+
+%% 先算正序的
 %平衡节点电流
 BalI1r=real(-sum(mIf1));
 BalI1i=imag(-sum(mIf1));
-inv(fsY11)*(mIf1);
-
+% inv(fsY11)*(mIf1);
 measurement=-mIf1(Loadi);
-clear PD QD PG QG;
+[ V1r,V1i,I1r,I1i ]=IPMLoop(measurement,BalI1r,BalI1i,busNum,Loadi,fsY1,Balance,1 );
+f=sum(([real(measurement);imag(measurement)]-[-I1r;-I1i]).^2);
+
+%% 算负序的
+BalI2r=real(-sum(mIf2));
+BalI2i=imag(-sum(mIf2));
+measurement=-mIf2(Loadi);
+[ V2r,V2i,I2r,I2i ]=IPMLoop(measurement,BalI2r,BalI2i,busNum,Loadi,fsY2,Balance,1 );
+f=sum([real(measurement);imag(measurement)]-[-I2r;-I2i]);
 %状态量
 % SEVoltpA=sparse(ones(busNum,1));
 % SEVoltpB=sparse(ones(busNum,1)).*exp(1j*-120/180*pi);
@@ -250,6 +260,5 @@ clear PD QD PG QG;
 % SEVaf1=sparse(zeros(busNum,1));
 % SEPD=sparse(zeros(busNum,1));
 % SEQD=sparse(zeros(busNum,1));
-[ V1r,V1i,I1r,I1i ]=IPMLoop(measurement,BalI1r,BalI1i,busNum,Loadi,fsY1,Balance );
+
 %检查目标函数
-f=sum([real(measurement);imag(measurement)]-[-I1r;-I1i]);
\ No newline at end of file