diff --git a/IPMLoop.m b/IPMLoop.m
index d62d9b5..8957b5a 100644
--- a/IPMLoop.m
+++ b/IPMLoop.m
@@ -30,7 +30,39 @@ isSetBound=0;
 largerBound=0;
 realBound=1;
 tic
-while(abs(Gap)>1e-4)
+
+Ir=real(guessIf1(ismember(Loadi,boundedLoadi)) );
+pIr=find(Ir>0);
+nIr=find(Ir<0);
+Ii=imag(guessIf1(ismember(Loadi,boundedLoadi)));
+pIi=find(Ii>0);
+nIi=find(Ii<0);
+lower=ones(length(boundedLoadi)*2,1);
+if realBound==1
+    isSetBound=1;
+    lower(pIr)=(0.7)*Ir(pIr);
+    lower(nIr)=(1.300)*Ir(nIr);
+    lower(pIi+length(Ir))=(0.7)*Ii(pIi);
+    lower(nIi+length(Ir))=(1.30)*Ii(nIi);
+end
+
+Ir=real(guessIf1(ismember(Loadi,boundedLoadi)));
+pIr=find(Ir>0);
+nIr=find(Ir<0);
+Ii=imag(guessIf1(ismember(Loadi,boundedLoadi)));
+pIi=find(Ii>0);
+nIi=find(Ii<0);
+upper=ones(length(boundedLoadi)*2,1);
+%     if abs(Gap)<1 || isSetBound==1
+if realBound==1
+    isSetBound=1;
+    upper(pIr)=(1.30)*Ir(pIr);
+    upper(nIr)=(0.7)*Ir(nIr);
+    upper(pIi+length(Ir))=(1.30)*Ii(pIi);
+    upper(nIi+length(Ir))=(0.7)*Ii(nIi);
+end
+
+while(abs(Gap)>1e-5)
     if KK>=kmax
         break;
     end
@@ -63,99 +95,7 @@ deltG=func_deltGstate1(busNum,Loadi,boundedLoadi,I1r,I1i);
     Mat_G=FormGstate1(I1r,I1i,ismember(Loadi,boundedLoadi) );
     Mat_H=FormH(fsY1,Loadi,V1r,V1i,I1r,I1i,BalI1r,BalI1i,Balance);
     Ly=Mat_H;
-    
-%     if isSetBound==0
-    
-        Ir=real(guessIf1(ismember(Loadi,boundedLoadi)) );
-        pIr=find(Ir>0);
-        nIr=find(Ir<0);
-        Ii=imag(guessIf1(ismember(Loadi,boundedLoadi)));
-        pIi=find(Ii>0);
-        nIi=find(Ii<0);
-%         lower=ones(length(Loadi)*2,1);
-lower=ones(length(boundedLoadi)*2,1);
-%     if abs(Gap)<1 || isSetBound==1
 
-%     if abs(Gap)<0.001
-%         largerBound=0;
-%         realBound=1;
-%     end
-
-     if realBound==1
-%         lower(pIr)=(0.9-3/(KK+1))*Ir(pIr);
-%         lower(nIr)=(1.100+3/(KK+1))*Ir(nIr);
-%         lower(pIi+length(Ir))=(0.9-3/(KK+1))*Ii(pIi);
-%         lower(nIi+length(Ir))=(1.10+3/(KK+1))*Ii(nIi);
-        isSetBound=1;
-        lower(pIr)=(0.7)*Ir(pIr);
-        lower(nIr)=(1.300)*Ir(nIr);
-        lower(pIi+length(Ir))=(0.7)*Ii(pIi);
-        lower(nIi+length(Ir))=(1.30)*Ii(nIi);
-%         lower( ismember(Loadi,noBoundedLoadi))=-1.9;
-%         lower(ismember(Loadi,noBoundedLoadi)+length(Ir))=-1.9;
-     end
-        
-    if largerBound==1
-        lower(pIr)=-2.8;
-        lower(nIr)=-2.8;
-        lower(pIi+length(Ir))=-2.8;
-        lower(nIi+length(Ir))=-2.8;
-        
-%         lower(pIr)=(0.7)*Ir(pIr)-0.1;
-%         lower(nIr)=(1.300)*Ir(nIr)-0.1;
-%         lower(pIi+length(Ir))=(0.7)*Ii(pIi)-0.1;
-%         lower(nIi+length(Ir))=(1.30)*Ii(nIi)-0.1;
-     end
-
-     
-     
-%     lowerR=lower(1:length(lower)/2);
-%     lowerR(I1r./lowerR>0.98)=lowerR(I1r./lowerR>0.98)-0.001;
-%     lowerI=lower(length(lower)/2+1:end);
-%     lowerI(I1i./lowerI>0.98)=lowerI(I1i./lowerI>0.98)-0.001;    
-%     lower=[lowerR;lowerI];
-
-%         I1r(I1r./lowerR>0.9998)=lowerR(I1r./lowerR>0.9998);
-% I1i(I1i./lowerI>0.9998)=lowerI(I1i./lowerI>0.9998);
-    
-        Ir=real(guessIf1(ismember(Loadi,boundedLoadi)));
-        pIr=find(Ir>0);
-        nIr=find(Ir<0);
-        Ii=imag(guessIf1(ismember(Loadi,boundedLoadi)));
-        pIi=find(Ii>0);
-        nIi=find(Ii<0);
-        upper=ones(length(boundedLoadi)*2,1);
-%     if abs(Gap)<1 || isSetBound==1
-   if realBound==1
-%         upper(pIr)=(1.10+3/(KK+1))*Ir(pIr);
-%         upper(nIr)=(0.9-3/(KK+1))*Ir(nIr);
-%         upper(pIi+length(Ir))=(1.10+3/(KK+1))*Ii(pIi);
-%         upper(nIi+length(Ir))=(0.9-3/(KK+1))*Ii(nIi);
-        isSetBound=1;
-        upper(pIr)=(1.30)*Ir(pIr);
-        upper(nIr)=(0.7)*Ir(nIr);
-        upper(pIi+length(Ir))=(1.30)*Ii(pIi);
-        upper(nIi+length(Ir))=(0.7)*Ii(nIi);
-%         upper(ismember(Loadi,noBoundedLoadi))=1.9;
-%         upper(ismember(Loadi,noBoundedLoadi)+length(Ir))=1.9;
-   end
-
-    if largerBound==1
-
-        upper(pIr)=2.8;
-        upper(nIr)=2.8;
-        upper(pIi+length(Ir))=2.8;
-        upper(nIi+length(Ir))=2.8;
-
-%         upper(pIr)=(1.30)*Ir(pIr)+0.1;
-%         upper(nIr)=(0.7)*Ir(nIr)+0.1;
-%         upper(pIi+length(Ir))=(1.30)*Ii(pIi)+0.1;
-%         upper(nIi+length(Ir))=(0.7)*Ii(nIi)+0.1;
-    end
-        
-        
-%     end
-    
     Lz=FormLzstate1(Loadi,Mat_G,Init_L,I1measurement,dI_F,flag,guessIf1,rIf1,lower);
     Lw=FormLwstate1(Loadi,Mat_G,Init_U,I1measurement,dI_F,flag,guessIf1,rIf1,upper);
     Lx=FormLx(deltF,deltH,Init_Y,deltG,Init_Z,Init_W);
@@ -173,7 +113,7 @@ end
 toc
 calcuTime=toc/KK;
 % f=sum(([real(I1measurement);imag(I1measurement)]-[I1r;I1i]).^2)+sum((real(V1measurement)-V1r(Loadi)).^2)+sum((imag(V1measurement)-V1i(Loadi)).^2);
-if abs(Gap)<1e-4 && KK<kmax
+if abs(Gap)<1e-5 && KK<kmax
     isConverged=1;
 else
     isConverged=0;
diff --git a/calIWi.m b/calIWi.m
new file mode 100644
index 0000000..7a6a787
--- /dev/null
+++ b/calIWi.m
@@ -0,0 +1,8 @@
+function [ output_args ] = calIWi( Vr,Vi,sigmaP,sigmaQ )
+dIrdP=(Vi./(Vr.^2+Vi.^2)).^2.*(sigmaP.^2);
+dIrdQ=(Vr./(Vr.^2+Vi.^2)).^2.*(sigmaQ.^2);
+%dIrdVr=(P.*(Vr.^2+Vi.^2)-2*(P.*Vr+Q.*Vi).*Vr)./(Vr.^2+Vi.^2);
+%dIrdVr=(Q.*(Vr.^2+Vi.^2)-2*(P.*Vr+Q.*Vi).*Vi)./(Vr.^2+Vi.^2);
+output_args=max(abs([dIrdP,dIrdQ]),[],2);
+end
+
diff --git a/calIWr.m b/calIWr.m
new file mode 100644
index 0000000..a4c9012
--- /dev/null
+++ b/calIWr.m
@@ -0,0 +1,9 @@
+function [ output_args ] = calIWr(Vr,Vi,sigmaP,sigmaQ )
+%计算电流实部的方差
+dIrdP=(Vr./(Vr.^2+Vi.^2)).^2.*(sigmaP.^2);
+dIrdQ=(Vi./(Vr.^2+Vi.^2)).^2.*(sigmaQ.^2);
+%dIrdVr=(P.*(Vr.^2+Vi.^2)-2*(P.*Vr+Q.*Vi).*Vr)./(Vr.^2+Vi.^2);
+%dIrdVr=(Q.*(Vr.^2+Vi.^2)-2*(P.*Vr+Q.*Vi).*Vi)./(Vr.^2+Vi.^2);
+output_args=max(abs([dIrdP,dIrdQ]),[],2);
+end
+
diff --git a/run.m b/run.m
index 72bbf9e..fd5d8dc 100644
--- a/run.m
+++ b/run.m
@@ -214,13 +214,21 @@ I1measurement=mIf1(Loadi);
 % wI1i=abs( imag(mIf1(Loadi).*sigmaI1) );
 %正序测量误差的第二种形式
 
-wI1r=abs(real(I1measurement)).*sigma;
+% wI1r=abs(real(I1measurement)).*sigma;
+% wI1i=abs(imag(I1measurement)).*sigma;
 
-% dIrdPr=real(V1measurement)./(abs(V1measurement).^2);
-% dIrdQr=imag(V1measurement)./(abs(V1measurement).^2);
-% dIrdV=
+IArW=calIWr(real(mVoltpA),imag(mVoltpA),mphaseASpotLoadP*sigma,mphaseASpotLoadQ*sigma );
+IAiW=calIWi(real(mVoltpA),imag(mVoltpA),mphaseASpotLoadP*sigma,mphaseASpotLoadQ*sigma );
+IBrW=calIWr(real(mVoltpB),imag(mVoltpB),mphaseBSpotLoadP*sigma,mphaseBSpotLoadQ*sigma );
+IBiW=calIWi(real(mVoltpB),imag(mVoltpB),mphaseBSpotLoadP*sigma,mphaseBSpotLoadQ*sigma );
+ICrW=calIWr(real(mVoltpC),imag(mVoltpC),mphaseCSpotLoadP*sigma,mphaseCSpotLoadQ*sigma );
+ICiW=calIWi(real(mVoltpC),imag(mVoltpC),mphaseCSpotLoadP*sigma,mphaseCSpotLoadQ*sigma );
+
+wI1r=sqrt((IArW+ICrW+IBrW)/9 );
+wI1r=wI1r(Loadi);
+wI1i=sqrt((IAiW+ICiW+IBiW)/9 );
+wI1i=wI1i(Loadi);
 
-wI1i=abs(imag(I1measurement)).*sigma;
 wV1r=sqrt(2*( (abs(mVoltpA).*sigma).^2+(abs(mVoltpB).*sigma).^2+(abs(mVoltpC).*sigma).^2 )/9 );
 wV1i=sqrt(2*( (abs(mVoltpA).*sigma).^2+(abs(mVoltpB).*sigma).^2+(abs(mVoltpC).*sigma).^2 )/9 );
 wV1r=wV1r(Loadi);
@@ -249,8 +257,14 @@ I2measurement=mIf2(Loadi);
 % wI2r=abs( real(mIf2(Loadi).*sigmaI2) );
 % wI2i=abs( imag(mIf2(Loadi).*sigmaI2) );
 %负序误差的第二种形式
-wI2r=abs(real(I2measurement)).*sigma;
-wI2i=abs(imag(I2measurement)).*sigma;
+% wI2r=abs(real(I2measurement)).*sigma;
+% wI2i=abs(imag(I2measurement)).*sigma;
+
+wI2r=sqrt((IArW+ICrW+IBrW)/9 );
+wI2r=wI2r(Loadi);
+wI2i=sqrt((IAiW+ICiW+IBiW)/9 );
+wI2i=wI2i(Loadi);
+
 wV2r=sqrt(2*( (abs(mVoltpA).*sigma).^2+(abs(mVoltpB).*sigma).^2+(abs(mVoltpC).*sigma).^2 )/9 );
 wV2i=sqrt(2*( (abs(mVoltpA).*sigma).^2+(abs(mVoltpB).*sigma).^2+(abs(mVoltpC).*sigma).^2 )/9 );
 wV2r=wV2r(Loadi);
@@ -280,9 +294,15 @@ I0measurement=mIf0(Loadi);
 % wI0r=abs( real(mIf0(Loadi).*sigmaI0) );
 % wI0i=abs( imag(mIf0(Loadi).*sigmaI0) );
 %零序误差的第二种形式
+% 
+% wI0r=abs(real(I0measurement)).*sigma;
+% wI0i=abs(imag(I0measurement)).*sigma;
+
+wI0r=sqrt((IArW+ICrW+IBrW)/9 );
+wI0r=wI0r(Loadi);
+wI0i=sqrt((IAiW+ICiW+IBiW)/9 );
+wI0i=wI0i(Loadi);
 
-wI0r=abs(real(I0measurement)).*sigma;
-wI0i=abs(imag(I0measurement)).*sigma;
 wV0r=sqrt(2*( (abs(mVoltpA).*sigma).^2+(abs(mVoltpB).*sigma).^2+(abs(mVoltpC).*sigma).^2 )/9 );
 wV0i=sqrt(2*( (abs(mVoltpA).*sigma).^2+(abs(mVoltpB).*sigma).^2+(abs(mVoltpC).*sigma).^2 )/9 );
 wV0r=wV0r(Loadi);