diff --git a/OPF.m b/OPF.m
index 479b3b9..8319786 100644
--- a/OPF.m
+++ b/OPF.m
@@ -7,14 +7,17 @@ sumCaseB_SE=0;
 sumCaseC_SE=0;
 badDataResult=zeros(10,33);
 badDataLocation=zeros(34,10);
-maxDVolt=0;
-maxDVAngle=0;
-% for badDataNode=1:33
-    for I=1:1
+nodeMaxDVolt=zeros(33,1);
+nodeMaxDVAngle=zeros(33,1);
+for badDataNode=1:1
+    loopN=1;
+    maxDVolt=0;
+    maxDVAngle=0;
+    while 1
         close
         [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,Liner,Linex,Lineb,Transforr,Transforx,Branchi,Branchg,Branchb,Transfork0]= ...
             pf('E:\算例\feeder33\feeder33.txt');
-        sigma=0.01;
+        sigma=0.001;
         RealPD=PD;
         RealQD=QD;
         rVolt=Volt;
@@ -33,23 +36,35 @@ maxDVAngle=0;
         % mVolt=load('mVolt');
         % mVolt=mVolt.mVolt;
         
-        mVolt(31)=rVolt(31)*(1-sigma*6);
+%         mVolt(badDataNode)=rVolt(badDataNode)*(1-sigma*6);
         %% 画Case A的图
         % figure('Color',[1 1 1]);
-        [Busnum,Loadi,Volt,PD,QD,rVolt,UAngel,RealPD,RealQD,rUAngel,Vbi,PDbi,QDbi,plotGapA]=subOPF([],PD0,QD0,mVolt,sigma);%全部有
-%         badDataResult(I,badDataNode)=sum(Vbi);
-%         badDataLocation(1:33,I)=Vbi;
-%         badDataLocation(34,I)=sum(abs((rVolt-Volt)./rVolt./length(rVolt)))+sum(abs( (UAngel(2:33)-rUAngel(2:33))./rUAngel(2:33)./length(rUAngel(2:33))));
+        [Busnum,Loadi,Volt,PD,QD,rVolt,UAngel,RealPD,RealQD,rUAngel,Vbi,PDbi,QDbi,plotGapA,isConverge]=subOPF([1:33],PD0,QD0,mVolt,sigma);%全部有
+        %         badDataResult(I,badDataNode)=sum(Vbi);
+        %         badDataLocation(1:33,I)=Vbi;
+        %         badDataLocation(34,I)=sum(abs((rVolt-Volt)./rVolt./length(rVolt)))+sum(abs( (UAngel(2:33)-rUAngel(2:33))./rUAngel(2:33)./length(rUAngel(2:33))));
+                break;
+        if isConverge==0
+            continue;
+        end
+        % 加最大偏差
+        maxDVolt_=max(abs((rVolt-Volt)));
+        if maxDVolt_>maxDVolt
+            maxDVolt=maxDVolt_;
+        end
+        maxDVAngle_=max(abs((UAngel(2:33)-rUAngel(2:33))));
+        if maxDVAngle_>maxDVAngle
+            maxDVAngle=maxDVAngle_;
+        end
+        
+        if loopN>=1
+            nodeMaxDVolt(badDataNode)=maxDVolt;
+            nodeMaxDVAngle(badDataNode)=maxDVAngle;
+            break;
+        end
 
-% 加最大偏差
-maxDVolt_=max(abs((rVolt-Volt)));
-if maxDVolt_>maxDVolt
-    maxDVolt=maxDVolt_;
-end
-maxDVAngle_=max(abs((UAngel(2:33)-rUAngel(2:33))));
-if maxDVAngle_>maxDVAngle
-    maxDVAngle=maxDVAngle_;
-end
+        loopN=loopN+1;
+    end
     
 end
 % end
diff --git a/subOPF.m b/subOPF.m
index 41656f0..a6fee55 100644
--- a/subOPF.m
+++ b/subOPF.m
@@ -1,9 +1,10 @@
-function [Busnum,Loadi,Volt,PD,QD,rVolt,UAngel,RealPD,RealQD,rUAngel,Vbi,PDbi,QDbi,plotGap]=subOPF(noMeasurei,PD0,QD0,mVolt,sigma)
+function [Busnum,Loadi,Volt,PD,QD,rVolt,UAngel,RealPD,RealQD,rUAngel,Vbi,PDbi,QDbi,plotGap,isConverge]=subOPF(noMeasurei,PD0,QD0,mVolt,sigma)
 tic
 %% 存在问题
 % 变压器变比的位置没有考虑，由于现在用的变比都是1，所以没有影响。 20130123
 %%
 % thesis=ForThesis(1,62);
+isConverge=1;%判断是否收敛
 [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,Liner,Linex,Lineb,Transforr,Transforx,Branchi,Branchg,Branchb,Transfork0]= ...
     pf('E:\算例\feeder33\feeder33.txt');
 %         pf('E:\算例\柳金Ⅰ926_21671693_2012-09-06\newFIle16.txt');
@@ -45,7 +46,7 @@ Gap=(Init_L*Init_Z'-Init_U*Init_W');
 KK=0;
 plotGap=zeros(1,60);
 ContrlCount=size(Loadi,1)*2+Busnum*2+Busnum+length(Loadi)*2;
-kmax=100;
+kmax=150;
 Precision=Precision/1;
 %% 加误差
 %找DG
@@ -53,6 +54,11 @@ DGi=find(PD0<0);
 %
 mPD=PD0;
 mQD=QD0;
+
+%负荷曲线有误
+% mVolt(3)=0.9*mVolt(3);%15%
+% mQD(13)=0.3*mQD(13);
+
 %% 加估计上下界,用真实值
 % lPD=abs(RealPD*3*sigma);
 % uPD=abs(RealPD*3*sigma);
@@ -86,8 +92,16 @@ lPD(noMeasurei)=0.15*mPD(noMeasurei);%15%
 uPD(noMeasurei)=0.15*mPD(noMeasurei);
 lQD(noMeasurei)=0.15*mQD(noMeasurei);
 uQD(noMeasurei)=0.15*mQD(noMeasurei);
+
+
+%电压伪量测错误
+mVolt(noMeasurei)=0.95;
+
 lVolt(noMeasurei)=0.7*mVolt(noMeasurei);%电压在0.93~1.07
 uVolt(noMeasurei)=0.7*mVolt(noMeasurei);
+
+
+
 %错误数据
 %mVolt(2)=5;
 % bigM=0.000003;
@@ -124,7 +138,7 @@ while(abs(Gap)>Precision*10)
     %% 形成方程矩阵
     Luu=Init_U'.*Init_W'+Init_u*ones(RestraintCount,1);
     Lul=Init_L'.*Init_Z'-Init_u*ones(RestraintCount,1);
-    bigM=0.8;
+    bigM=1;
     Mat_G=FormG(Volt,PD,QD,Loadi,bigM,mVolt,rVolt,sigma,Vbi,PDbi,QDbi,mPD,mQD,uPD,lPD,uQD,lQD,uVolt,lVolt);
     Mat_H=FormH(Busnum,Volt,PG,PD,QG,QD,Y,UAngel,r,c,Angle,Loadi);
     Ly=Mat_H;
@@ -177,7 +191,7 @@ while(abs(Gap)>Precision*10)
     %% 形成方程矩阵
     Luu=Init_U'.*Init_W'+Init_u*ones(RestraintCount,1);
     Lul=Init_L'.*Init_Z'-Init_u*ones(RestraintCount,1);
-    bigM=0.4;
+    bigM=1;
     eps=Gap*0.1;
     Mat_G=FormG(Volt,PD,QD,Loadi,bigM,mVolt,rVolt,sigma,Vbi,PDbi,QDbi,mPD,mQD,uPD,lPD,uQD,lQD,uVolt,lVolt);
     Mat_H=FormH(Busnum,Volt,PG,PD,QG,QD,Y,UAngel,r,c,Angle,Loadi);
@@ -204,5 +218,12 @@ dQ=QD+diag(Volt)*Y.*sin(AngleIJ)*Volt';%
 dQ(Balance)=0;
 mdP=max(dP);
 mdQ=max(dQ);
+%判断是否收敛
+if min(abs(Vbi-.1))>5e-005
+    isConverge=0;
+end
+if KK>=kmax 
+    isConverge=0;
+end
 toc
 end
\ No newline at end of file