diff --git a/IPMLoop.m b/IPMLoop.m
index eaaa770..421fb72 100644
--- a/IPMLoop.m
+++ b/IPMLoop.m
@@ -30,7 +30,7 @@ kmax=15;
 isSetBound=0;
 largerBound=0;
 realBound=1;
-tic
+
 
 Ir=real(guessIf1(ismember(Loadi,boundedLoadi)) );
 pIr=find(Ir>0);
@@ -62,7 +62,7 @@ if realBound==1
     upper(pIi+length(Ir))=(1.30)*Ii(pIi);
     upper(nIi+length(Ir))=(0.7)*Ii(nIi);
 end
-
+tic
 while(abs(Gap)>1e-5)
     if KK>=kmax
         break;
@@ -108,11 +108,11 @@ deltG=func_deltGstate1(busNum,Loadi,boundedLoadi,I1r,I1i);
     [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);+max(abs(deltX));
 %     Gap=max([max(abs(Mat_H)),max(Lz),max(Lw) ]);
-%     fprintf('Gap %f :%d\n',full(Gap),KK+1);
+    fprintf('Gap %f :%d\n',full(Gap),KK+1);
     KK=KK+1;
 end
-toc
-calcuTime=toc/KK;
+% toc
+calcuTime=toc;
 % 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-5 && KK<kmax
     isConverged=1;
diff --git a/run.m b/run.m
index 1a055f3..5540935 100644
--- a/run.m
+++ b/run.m
@@ -1,4 +1,4 @@
-function [JMeasurement,AME_Volt,AME_VAngle,AME_PD,AME_QD,AME_mVolt,AME_mPD,AME_mQD,isConverged,totalTime1,totalTime2,totalTime0]=run()
+function [JMeasurement,AME_Volt,AME_VAngle,AME_PD,AME_QD,AME_mVolt,AME_mPD,AME_mQD,isConverged,totalTime1,totalTime2,totalTime0,eV2,eV0,eI2,eI0]=run()
 %% 利用先把负荷转换为电流的方法。这个方法要求知道电压量。
 %
 close all
@@ -9,36 +9,36 @@ lineZ=readLineZ('modified-feeder69\lineParameter.txt');
     phaseASpotLoadQ,phaseBSpotLoadQ,phaseCSpotLoadQ,setIJ,nodeNum,Balance,phaseABCY ...
     cap]=dataRead(lineZ,'modified-feeder69\data.txt');
 %负荷不平衡
-% for I=1:length(phaseASpotLoadP)
-%     roll=mod(round(rand()*10),3);
-%     if roll==0
-%         phaseBSpotLoadP(I)=phaseBSpotLoadP(I)+phaseASpotLoadP(I)*.45;
-%         phaseCSpotLoadP(I)=phaseCSpotLoadP(I)+phaseASpotLoadP(I)*.45;
-%         phaseASpotLoadP(I)=phaseASpotLoadP(I)*.10;
-%         
-%         phaseBSpotLoadQ(I)=phaseBSpotLoadQ(I)+phaseASpotLoadQ(I)*.45;
-%         phaseCSpotLoadQ(I)=phaseCSpotLoadQ(I)+phaseASpotLoadQ(I)*.45;
-%         phaseASpotLoadQ(I)=phaseASpotLoadQ(I)*.10;
-%     end
-%     if roll==1
-%         phaseASpotLoadP(I)=phaseASpotLoadP(I)+phaseBSpotLoadP(I)*.45;
-%         phaseCSpotLoadP(I)=phaseCSpotLoadP(I)+phaseBSpotLoadP(I)*.45;
-%         phaseBSpotLoadP(I)=phaseBSpotLoadP(I)*.10;
-%         
-%         phaseASpotLoadQ(I)=phaseASpotLoadQ(I)+phaseBSpotLoadQ(I)*.45;
-%         phaseCSpotLoadQ(I)=phaseCSpotLoadQ(I)+phaseBSpotLoadQ(I)*.45;
-%         phaseBSpotLoadQ(I)=phaseBSpotLoadQ(I)*.10;
-%     end
-%     if roll==2
-%         phaseASpotLoadP(I)=phaseASpotLoadP(I)+phaseCSpotLoadP(I)*.45;
-%         phaseBSpotLoadP(I)=phaseBSpotLoadP(I)+phaseCSpotLoadP(I)*.45;
-%         phaseCSpotLoadP(I)=phaseCSpotLoadP(I)*.10;
-%         
-%         phaseASpotLoadQ(I)=phaseASpotLoadQ(I)+phaseCSpotLoadQ(I)*.45;
-%         phaseBSpotLoadQ(I)=phaseBSpotLoadQ(I)+phaseCSpotLoadQ(I)*.45;
-%         phaseCSpotLoadQ(I)=phaseCSpotLoadQ(I)*.10;
-%     end
-% end
+for I=1:length(phaseASpotLoadP)
+    roll=mod(round(rand()*10),3);
+    if roll==0
+        phaseBSpotLoadP(I)=phaseBSpotLoadP(I)+phaseASpotLoadP(I)*.45;
+        phaseCSpotLoadP(I)=phaseCSpotLoadP(I)+phaseASpotLoadP(I)*.45;
+        phaseASpotLoadP(I)=phaseASpotLoadP(I)*.10;
+        
+        phaseBSpotLoadQ(I)=phaseBSpotLoadQ(I)+phaseASpotLoadQ(I)*.45;
+        phaseCSpotLoadQ(I)=phaseCSpotLoadQ(I)+phaseASpotLoadQ(I)*.45;
+        phaseASpotLoadQ(I)=phaseASpotLoadQ(I)*.10;
+    end
+    if roll==1
+        phaseASpotLoadP(I)=phaseASpotLoadP(I)+phaseBSpotLoadP(I)*.45;
+        phaseCSpotLoadP(I)=phaseCSpotLoadP(I)+phaseBSpotLoadP(I)*.45;
+        phaseBSpotLoadP(I)=phaseBSpotLoadP(I)*.10;
+        
+        phaseASpotLoadQ(I)=phaseASpotLoadQ(I)+phaseBSpotLoadQ(I)*.45;
+        phaseCSpotLoadQ(I)=phaseCSpotLoadQ(I)+phaseBSpotLoadQ(I)*.45;
+        phaseBSpotLoadQ(I)=phaseBSpotLoadQ(I)*.10;
+    end
+    if roll==2
+        phaseASpotLoadP(I)=phaseASpotLoadP(I)+phaseCSpotLoadP(I)*.45;
+        phaseBSpotLoadP(I)=phaseBSpotLoadP(I)+phaseCSpotLoadP(I)*.45;
+        phaseCSpotLoadP(I)=phaseCSpotLoadP(I)*.10;
+        
+        phaseASpotLoadQ(I)=phaseASpotLoadQ(I)+phaseCSpotLoadQ(I)*.45;
+        phaseBSpotLoadQ(I)=phaseBSpotLoadQ(I)+phaseCSpotLoadQ(I)*.45;
+        phaseCSpotLoadQ(I)=phaseCSpotLoadQ(I)*.10;
+    end
+end
 
 a=exp(1j*2*pi/3);
 Tp2f=1/3*[1 1 1;
@@ -150,7 +150,14 @@ end
 FortiscueToc=toc;
 fprintf('Fortiscue法计算时间 %f\n',FortiscueToc);
 Vf1=Vmf1.*exp(1j*Vaf1);
-%%
+%% 计算三相不平衡度
+%电压不平衡度
+eV2=mean(abs(Vf2)./abs(Vf1))/length(Vf1)/3;
+eV0=mean(abs(Vf0)./abs(Vf1))/length(Vf1)/3;
+%电流不平衡度
+eI2=mean(abs(If2(find(If2)))./abs(If1(find(If1))))/length(find(If1))/3;
+eI0=mean(abs(If0(find(If0)))./abs(If1(find(If1))))/length(find(If1))/3;
+%% 
 (Vf0.*conj(fsY00*Vf0)+Vf1.*conj(fsY11*Vf1)+Vf2.*conj(fsY22*Vf2))*3;%包含补偿电容的功率
 conj(Tf2p*[If0(2);If1(2);If2(2)]).*(Tf2p*[Vf0(2);Vf1(2);Vf2(2)]);
 IpABC=Tf2p*conj([If0';If1';If2']);
@@ -356,7 +363,7 @@ noLoadi=Loadi(1:6:end);
 % noLoadi=[1,19];
 % noLoadi=[4,3,2,5,6,7,11,14,22];
 noLoadi=Loadi;
-noLoadi=[5];
+% noLoadi=[5];
 wI0r(ismember(Loadi,noLoadi))=1e100000;
 wI0i(ismember(Loadi,noLoadi))=1e100000;
 wI1r(ismember(Loadi,noLoadi))=1e100000;
diff --git a/test.m b/test.m
index 3f3eb2b..4156aff 100644
--- a/test.m
+++ b/test.m
@@ -14,10 +14,14 @@ totalTime2Sum=10000;
 totalTime0Sum=10000;
 timeOfPrepare=0.00272333333;%准备数据的时间
 timeOfPrepare=0;
+eV2Sum=0;
+eV0Sum=0;
+eI2Sum=0;
+eI0Sum=0;
 N=200;
 loopN=1;
 while 1
-    [JMeasurement,AME_Volt,AME_VAngle,AME_PD,AME_QD,AME_mVolt,AME_mPD,AME_mQD,isConverged,totalTime1,totalTime2,totalTime0]=run();
+    [JMeasurement,AME_Volt,AME_VAngle,AME_PD,AME_QD,AME_mVolt,AME_mPD,AME_mQD,isConverged,totalTime1,totalTime2,totalTime0,eV2,eV0,eI2,eI0]=run();
     if isConverged==0
         continue
     end
@@ -33,6 +37,11 @@ while 1
     AME_mVoltSum=AME_mVoltSum+AME_mVolt;
     AME_mPDSum=AME_mPDSum+AME_mPD;
     AME_mQDSum=AME_mQDSum+AME_mQD;
+    eV2Sum=eV2Sum+eV2;
+    eV0Sum=eV0Sum+eV0;
+    eI2Sum=eI2Sum+eI2;
+    eI0Sum=eI0Sum+eI0;
+    
     %时间 
 %     totalTime1Sum=totalTime1Sum+totalTime1;
 % totalTime2Sum=totalTime2Sum+totalTime2;
@@ -49,6 +58,10 @@ AME_QDSum=AME_QDSum/N;
 AME_mVoltSum=AME_mVoltSum/N;
 AME_mPDSum=AME_mPDSum/N;
 AME_mQDSum=AME_mQDSum/N;
+eV2Sum=eV2Sum/N;
+eV0Sum=eV0Sum/N;
+eI2Sum=eI2Sum/N;
+eI0Sum=eI0Sum/N;
 
         totalTime1Sum=totalTime1Sum+timeOfPrepare;
 totalTime2Sum=totalTime2Sum+timeOfPrepare;