diff --git a/OPF.m b/OPF.m
index 26ddf60..f20e22d 100644
--- a/OPF.m
+++ b/OPF.m
@@ -29,21 +29,21 @@ for badDataNode=1:1
         mVolt=rVolt.*(1+normrnd(0,sigma,length(rVolt),1))';
         
         %加载保存的变量
-        % PD0=load('PD0');
-        % PD0=PD0.PD0;
-        % QD0=load('QD0');
-        % QD0=QD0.QD0;
-        % mVolt=load('mVolt');
-        % mVolt=mVolt.mVolt;
+        PD0=load('PD0');
+        PD0=PD0.PD0;
+        QD0=load('QD0');
+        QD0=QD0.QD0;
+        mVolt=load('mVolt');
+        mVolt=mVolt.mVolt;
         
 %         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,isConverge]=subOPF([19:22],PD0,QD0,mVolt,sigma);%全部有
+        [Busnum,Loadi,Volt,PD,QD,rVolt,UAngel,RealPD,RealQD,rUAngel,Vbi,PDbi,QDbi,plotGapA,isConverge]=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))));
-%                 break;
+                break;
         if isConverge==0
             continue;
         end
@@ -74,10 +74,10 @@ for badDataNode=1:1
 end
 % end
 subplot(4,1,1,'XTick',[1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33]);
-CaseAREV=(Volt-rVolt)*100;%Relative Error of Voltage in Case A
+CaseAREV=(Volt-rVolt);%Relative Error of Voltage in Case A
 % CaseAREV=CaseAREV(2:end)*100;
 %真实值的
-plot(1:length(CaseAREV),abs(CaseAREV),'k.:','Marker','diamond');
+plot(1:length(CaseAREV),(CaseAREV),'k.:','Marker','diamond');
 %测量值的
 % plot(1:length(CaseAREV),abs((mVolt-rVolt)*100),'c.:','Marker','diamond');
 box off;
@@ -85,19 +85,21 @@ set(gca,'XTick',[1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24
 xlabel('节点号');
 ylabel('误差%');
 subplot(4,1,2);
-CaseAREA=(UAngel-rUAngel)*100;%Relative Error of Angle in Case A
-% CaseAREA=CaseAREA(2:end)*100;
+% CaseAREA=(UAngel-rUAngel)*100;%Relative Error of Angle in Case A
+CaseAREA=(UAngel-rUAngel);%Relative Error of Angle in Case A
+CaseAREA(1)=0;
 %真实值的
-plot(1:length(CaseAREA),abs(CaseAREA),'k:','Marker','diamond');
+plot(1:length(CaseAREA),(CaseAREA),'k:','Marker','diamond');
 box off;
 set(gca,'XTick',[1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33])
 xlabel('节点号');
 ylabel('误差%');
 subplot(4,1,3);
-CaseAREP=(PD-RealPD)./(RealPD+0.00001)*100;%Relative Error of PD in Case A
-% CaseAREP=CaseAREP(2:end)*100;
+% CaseAREP=(PD-RealPD)./(RealPD+0.00001)*100;%Relative Error of PD in Case A
+CaseAREP=(PD-RealPD)./RealPD*100;%Relative Error of PD in Case A
+CaseAREP(1)=0;
 %真实值的
-plot(1:length(CaseAREP),abs(CaseAREP),'k:','Marker','diamond');
+plot(1:length(CaseAREP),(CaseAREP),'k:','Marker','diamond');
 %测量值的
 % plot(1:length(CaseAREV),abs((PD0-RealPD)./(RealPD+0.00001)*100),'c.:','Marker','diamond');
 box off;
@@ -105,10 +107,11 @@ set(gca,'XTick',[1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24
 xlabel('节点号');
 ylabel('误差%');
 subplot(4,1,4);
-CaseAREQ=(QD-RealQD)./(RealQD+0.00001)*100;%Relative Error of QD in Case A
-% CaseAREQ=CaseAREQ(2:end)*100;
+% CaseAREQ=(QD-RealQD)./(RealQD+0.00001)*100;%Relative Error of QD in Case A
+CaseAREQ=(QD-RealQD)./RealQD*100;%Relative Error of QD in Case A
+CaseAREQ(1)=0;
 %真实值的
-plot(1:length(CaseAREQ),abs(CaseAREQ),'k:','Marker','diamond');
+plot(1:length(CaseAREQ),(CaseAREQ),'k:','Marker','diamond');
 %测量值的
 % plot(1:length(CaseAREV),abs((QD0-RealQD)./(RealQD+0.00001)*100),'c.:','Marker','diamond');
 box off;
@@ -119,7 +122,7 @@ ylabel('
 CaseAE=sqrt((sum(CaseAREV.^2)+sum(CaseAREA.^2)+sum(CaseAREP.^2)+sum(CaseAREQ.^2))/132);
 objA=full(sum(Vbi)+sum(PDbi)+sum(QDbi));
 notZeros=find(PD0~=0);
-CaseA_SE=sum(((PD0(notZeros)-PD(notZeros))./(PD0(notZeros)*sigma)).^2)+sum(((QD0(notZeros)-QD(notZeros))./(QD0(notZeros)*sigma)).^2)+sum(((Volt-mVolt)./(mVolt*sigma)).^2);
+CaseA_SE=sum(((RealPD(notZeros)-PD(notZeros))./(RealPD(notZeros)*1)).^2)+sum(((RealQD(notZeros)-QD(notZeros))./(RealQD(notZeros)*1)).^2)+sum(((Volt-rVolt)./(rVolt*1)).^2);
 CaseA_SE=(CaseA_SE/(length(notZeros)*2+length(Volt)))^.5;
 sumCaseA_SE=sumCaseA_SE+CaseA_SE;
 % arrayA(1:19,I)=Vbi;
@@ -129,54 +132,58 @@ sumCaseA_SE=sumCaseA_SE+CaseA_SE;
 [Busnum,Loadi,Volt,PD,QD,rVolt,UAngel,RealPD,RealQD,rUAngel,Vbi,PDbi,QDbi,plotGapB]=subOPF(setdiff(1:Busnum,[2,3,5,20,24,27,28,10,11,12,13]),PD0,QD0,mVolt,sigma);%
 subplot(4,1,1);
 hold on;
-CaseBREV=(Volt-rVolt)*100;%Relative Error of Voltage in Case B
-% CaseBREV=CaseBREV(2:end)*100;
-plot(1:length(CaseBREV),abs(CaseBREV),'b.:','Marker','square');
+% CaseBREV=(Volt-rVolt)*100;%Relative Error of Voltage in Case B
+CaseBREV=(Volt-rVolt);%Relative Error of Voltage in Case B
+plot(1:length(CaseBREV),(CaseBREV),'b.:','Marker','square');
 subplot(4,1,2);
 hold on;
-CaseBREA=(UAngel-rUAngel)*100;%Relative Error of Angle in Case B
-% CaseBREA=CaseBREA(2:end)*100;
-plot(1:length(CaseBREA),abs(CaseBREA),'b:','Marker','square');
+% CaseBREA=(UAngel-rUAngel)*100;%Relative Error of Angle in Case B
+CaseBREA=(UAngel-rUAngel);%Relative Error of Angle in Case B
+plot(1:length(CaseBREA),(CaseBREA),'b:','Marker','square');
 subplot(4,1,3);
 hold on;
-CaseBREP=(PD-RealPD)./(RealPD+0.00001)*100;%Relative Error of PD in Case B
-% CaseBREP=CaseBREP(2:end)*100;
-plot(1:length(CaseBREP),abs(CaseBREP),'b:','Marker','square');
+% CaseBREP=(PD-RealPD)./(RealPD+0.00001)*100;%Relative Error of PD in Case B
+CaseBREP=(PD-RealPD)./RealPD;%Relative Error of PD in Case B
+RealPD(1)=0;
+plot(1:length(CaseBREP),(CaseBREP),'b:','Marker','square');
 subplot(4,1,4);
 hold on;
-CaseBREQ=(QD-RealQD)./(RealQD+0.00001)*100;%Relative Error of QD in Case B
-% CaseBREQ=CaseBREQ(2:end)*100;
-plot(1:length(CaseBREQ),abs(CaseBREQ),'b:','Marker','square');
+% CaseBREQ=(QD-RealQD)./(RealQD+0.00001)*100;%Relative Error of QD in Case B
+CaseBREQ=(QD-RealQD)./RealQD;%Relative Error of QD in Case B
+CaseBREQ(1)=0;
+plot(1:length(CaseBREQ),(CaseBREQ),'b:','Marker','square');
 CaseBE=sqrt((sum(CaseBREV.^2)+sum(CaseBREA.^2)+sum(CaseBREP.^2)+sum(CaseBREQ.^2))/132);
 objB=full(sum(Vbi)+sum(PDbi)+sum(QDbi));
 noMeasurei=[2,3,5,20,24,27,28,10,11,12,13];
 Measurei=setdiff(2:33,[2,3,5,20,24,27,28,10,11,12,13]);
-CaseB_SE=sum(((PD0(Measurei)-PD(Measurei))./(PD0(Measurei)*sigma)).^2)+sum(((QD0(Measurei)-QD(Measurei))./(QD0(Measurei)*sigma)).^2)+sum(((Volt-mVolt)./(mVolt*sigma)).^2);
-CaseB_SE=CaseB_SE+sum(((PD0(noMeasurei)-PD(noMeasurei))./(PD0(noMeasurei)*20*sigma)).^2)+sum(((QD0(noMeasurei)-QD(noMeasurei))./(QD0(noMeasurei)*20*sigma)).^2);
-CaseB_SE=(CaseB_SE/(length(notZeros)*2+length(Volt)))^.5;
+CaseB_SE=sum(((RealPD(Measurei)-PD(Measurei))./(RealPD(Measurei)*1)).^2)+sum(((RealQD(Measurei)-QD(Measurei))./(RealQD(Measurei)*1)).^2)+sum(((Volt(Measurei)-rVolt(Measurei))./(rVolt(Measurei)*1)).^2);
+CaseB_SE=CaseB_SE+sum(((RealPD(noMeasurei)-PD(noMeasurei))./(RealPD(noMeasurei)*1)).^2)+sum(((RealQD(noMeasurei)-QD(noMeasurei))./(RealQD(noMeasurei)*1)).^2)+sum(((Volt(noMeasurei)-rVolt(noMeasurei))./(rVolt(noMeasurei)*1)).^2);
+CaseB_SE=(CaseB_SE/(length(notZeros)+length(noMeasurei) +length(Volt)))^.5;
 sumCaseB_SE=sumCaseB_SE+CaseB_SE;
 %% Case C
 [Busnum,Loadi,Volt,PD,QD,rVolt,UAngel,RealPD,RealQD,rUAngel,Vbi,PDbi,QDbi,plotGapC]=subOPF([1:33],PD0,QD0,mVolt,sigma);%
 subplot(4,1,1);
 hold on;
-CaseCREV=(Volt-rVolt)*100;%Relative Error of Voltage in Case C
-% CaseCREV=CaseCREV(2:end)*100;
-plot(1:length(CaseCREV),abs(CaseCREV),'r.:','Marker','o');
+% CaseCREV=(Volt-rVolt)*100;%Relative Error of Voltage in Case C
+CaseCREV=(Volt-rVolt);
+plot(1:length(CaseCREV),(CaseCREV),'r.:','Marker','o');
 subplot(4,1,2);
 hold on;
-CaseCREA=(UAngel-rUAngel)*100;%Relative Error of Angle in Case C
-% CaseCREA=CaseCREA(2:end)*100;
-plot(1:length(CaseCREA),abs(CaseCREA),'r:','Marker','o');
+% CaseCREA=(UAngel-rUAngel)*100;%Relative Error of Angle in Case C
+CaseCREA=(UAngel-rUAngel);
+plot(1:length(CaseCREA),(CaseCREA),'r:','Marker','o');
 subplot(4,1,3);
 hold on;
-CaseCREP=(PD-RealPD)./(RealPD+0.00001)*100;%Relative Error of PD in Case C
-% CaseCREP=CaseCREP(2:end)*100;
-plot(1:length(CaseCREP),abs(CaseCREP),'r:','Marker','o');
+% CaseCREP=(PD-RealPD)./(RealPD+0.00001)*100;%Relative Error of PD in Case C
+CaseCREP=(PD-RealPD)./RealPD*100;%Relative Error of PD in Case C
+CaseCREP(1);
+plot(1:length(CaseCREP),(CaseCREP),'r:','Marker','o');
 subplot(4,1,4);
 hold on;
-CaseCREQ=(QD-RealQD)./(RealQD+0.00001)*100;%Relative Error of QD in Case C
-% CaseCREQ=CaseCREQ(2:end)*100;
-plot(1:length(CaseCREQ),abs(CaseCREQ),'r:','Marker','o');
+% CaseCREQ=(QD-RealQD)./(RealQD+0.00001)*100;%Relative Error of QD in Case C
+CaseCREQ=(QD-RealQD)./RealQD*100;%Relative Error of QD in Case C
+CaseCREQ(1)=0;
+plot(1:length(CaseCREQ),(CaseCREQ),'r:','Marker','o');
 % 画legend
 subplot(4,1,1);
 % title('Voltage');
@@ -197,7 +204,7 @@ set(ld,'Position',[0.847865087908145 0.186094477711244 0.0543595263724435 0.0605
 CaseCE=sqrt((sum(CaseCREV.^2)+sum(CaseCREA.^2)+sum(CaseCREP.^2)+sum(CaseCREQ.^2))/132);
 objC=full(sum(Vbi)+sum(PDbi)+sum(QDbi));
 
-CaseC_SE=sum(((PD0(notZeros)-PD(notZeros))./(PD0(notZeros)*sigma)).^2)+sum(((QD0(notZeros)-QD(notZeros))./(QD0(notZeros)*sigma)).^2)+sum(((Volt-mVolt)./(mVolt*sigma)).^2);
+CaseC_SE=sum(((RealPD(notZeros)-PD(notZeros))./(RealPD(notZeros)*1)).^2)+sum(((RealQD(notZeros)-QD(notZeros))./(RealQD(notZeros)*1)).^2)+sum(((Volt-rVolt)./(rVolt*1)).^2);
 CaseC_SE=(CaseC_SE/(length(notZeros)*2+length(Volt)))^.5;
 sumCaseC_SE=sumCaseC_SE+CaseC_SE;
 
@@ -354,7 +361,7 @@ yy=hist(y,x); %
 bar(x,yy) %画出概率密度分布图
 %画收敛曲线
 fz=find(abs(plotGapA)==0);
-fz=fz(1);
+% fz=fz(1);
 figure('Name','互补曲线')
 plot(1:fz-1,plotGapA(1:fz-1));
 figure('Name','最大不平衡量');
diff --git a/PD0.mat b/PD0.mat
index db6d1f3..b345fc1 100644
Binary files a/PD0.mat and b/PD0.mat differ
diff --git a/QD0.mat b/QD0.mat
index ed2f708..fc08891 100644
Binary files a/QD0.mat and b/QD0.mat differ
diff --git a/mVolt.mat b/mVolt.mat
index 6969da9..e6c2925 100644
Binary files a/mVolt.mat and b/mVolt.mat differ
diff --git a/openfile2.m b/openfile2.m
index 10ac88e..d0b7693 100644
--- a/openfile2.m
+++ b/openfile2.m
@@ -59,7 +59,7 @@ QD=sparse(QD);
 
 %%
 %增加大负荷
-PD(22)=PD(22)*65;
+% PD(22)=PD(22)*65;
 %%
 
 PG=sparse(PG);
diff --git a/subOPF.m b/subOPF.m
index 9beec7b..5eff289 100644
--- a/subOPF.m
+++ b/subOPF.m
@@ -43,7 +43,7 @@ QG(PVi)=QGBal(PVi);
 %%
 [Volt,UAngel,Init_Z,Init_W,Init_L,Init_U,Init_Y,PG,QG,RestraintCount,wPG,wQG,wPD,wQD,PD,PD00,QD,randPDind,Loadi]=OPF_Init(Busnum,Balance,PG,QG,Volt,GenU,GenL,PVi,PGi,PVQU,PVQL,RealPD,RealQD,QD,PD);
 
-Volt=0.99*Volt;
+% Volt=0.99*Volt;
 
 Gap=(Init_L*Init_Z'-Init_U*Init_W');
 KK=0;
@@ -106,10 +106,10 @@ uQD(noMeasurei)=0.15*mQD(noMeasurei);
 
 
 %电压伪量测错误
-mVolt(noMeasurei)=.95;
+mVolt(noMeasurei)=1;
 
-lVolt(noMeasurei)=0.07*mVolt(noMeasurei);%电压在0.93~1.07
-uVolt(noMeasurei)=0.07*mVolt(noMeasurei);
+lVolt(noMeasurei)=0.12*mVolt(noMeasurei);%电压在0.93~1.07
+uVolt(noMeasurei)=0.02*mVolt(noMeasurei);
 
 
 
@@ -150,7 +150,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.5;
+    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;
@@ -175,11 +175,11 @@ Gap=1;
 % eps=0.00001;
 fprintf('\n');
 
-Init_Z=sparse(ones(1,RestraintCount));
-Init_W=sparse(-1*ones(1,RestraintCount));
-Init_L=1*sparse(ones(1,RestraintCount));
-Init_U=1*sparse(ones(1,RestraintCount));
-Init_Y=sparse(1,2*Busnum);%与学姐一致
+% Init_Z=sparse(ones(1,RestraintCount));
+% Init_W=sparse(-1*ones(1,RestraintCount));
+% Init_L=1*sparse(ones(1,RestraintCount));
+% Init_U=1*sparse(ones(1,RestraintCount));
+% Init_Y=sparse(1,2*Busnum);%与学姐一致
 
 fprintf('第2次迭代，算离散量。\n');
 while(abs(Gap)>Precision*10)
@@ -209,7 +209,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.5;
+    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);
@@ -234,8 +234,8 @@ dP=PD+diag(Volt)*Y.*cos(AngleIJ)*Volt';
 dP(Balance)=0;
 dQ=QD+diag(Volt)*Y.*sin(AngleIJ)*Volt';%不为0 ，因为没有考虑电容的QG
 dQ(Balance)=0;
-mdP=max(dP);
-mdQ=max(dQ);
+mdP=max(dP)
+mdQ=max(dQ)
 %判断是否收敛
 if min(abs(Vbi-.1))>5e-005
     isConverge=0;