diff --git a/@ForThesis/MaxBranchDeviation.asv b/@ForThesis/MaxBranchDeviation.asv
deleted file mode 100644
index ac15e4c..0000000
--- a/@ForThesis/MaxBranchDeviation.asv
+++ /dev/null
@@ -1,41 +0,0 @@
-function [ output_args ] = MaxBranchDeviation( ~, Linei,Linej,Liner,Linex,Lineb,Transfori,Transforj,Transforr,Transforx,Branchi,Branchg,Branchb,Transfork0,Volt0,UAngel0,Volt,UAngel,FileName,PD0,QD0)
-%% 最大支路功率偏差
-% 支路功率包括线路和变压器
-%% 先用加了误差的负荷功率计算潮流值
-[Busnum,Balance,PQstandard,Precision,~,~,~,~,~,kmax,~ ,...
-    ~,~,~,~,~,~,~,Pointpoweri,PG,QG,PD,QD,PVi,PVu,~,~,~,~,~,~,~]= openfile2(FileName);
-PD=PD0;
-QD=QD0;
-%% 形成节点导纳矩阵
-[~,Y,r,c,Angle] = admmatrix(Busnum,Linei,Linej,Liner,Linex,Lineb,Transfori,Transforj,Transforr,...
-    Transforx,Transfork0,Branchi,Branchg,Branchb);
-[P0,Q0,U,Uangle] = Initial(PG,PD,PQstandard,Pointpoweri,QG,QD,Busnum);                         %求功率不平衡量
-%disp('迭代次数i  最大不平衡量');
-%% 循环体计算
-for i = 0:kmax
-    [Jacob,PQ,U,Uangle] = jacobian(Busnum,Balance,PVi,PVu,U,Uangle,Y,Angle,P0,Q0,r,c);  %形成雅克比矩阵
-   % disp('第一次雅克比');
-    m = max(abs(PQ));
-    m=full(m);
-    %fprintf('       %u         %.8f \n',i,m);
-    if m > Precision                                                                                                                                 %判断不平衡量是否满足精度要求
-        [Uangle,U] = solvefun(Busnum,Jacob,PQ,Uangle,U);                                                           %求解修正方程，更新电压变量
-    else
-        disp(['收敛，迭代次数为',num2str(i),'次']);
-        break                                                                                                                                             %若满足精度要求，则计算收敛
-    end
-end
-
-Volt0=U;
-UAngel0=Uangle;
-
-[dispLineloss0 dispTransloss0]=Lineloss(Linei,Linej,Liner,Linex,Lineb,Transfori,Transforj,Transforr,Transforx,Branchi,Branchg,Branchb,Transfork0,Volt0,UAngel0);
-[dispLineloss dispTransloss]=Lineloss(Linei,Linej,Liner,Linex,Lineb,Transfori,Transforj,Transforr,Transforx,Branchi,Branchg,Branchb,Transfork0,Volt,UAngel);
-t1=(dispLineloss0(:,3) - dispLineloss(:,3))./dispLineloss0(:,3);
-t2=(dispTransloss0(:,3) - dispTransloss(:,3))./dispTransloss0(:,3);
-t11=dispLineloss0(:,3)>1e-5;% 太小的值不计算
-t22=dispTransloss0(:,3)>1e-5;% 太小的值不计算
-t3=abs([t1(t11);t2(t22)]);
-output_args=max(t3(t3~=Inf));
-end
-
diff --git a/@ForThesis/MaxDeviation.asv b/@ForThesis/MaxDeviation.asv
deleted file mode 100644
index bcad661..0000000
--- a/@ForThesis/MaxDeviation.asv
+++ /dev/null
@@ -1,7 +0,0 @@
-function MaxDeviation(this,PG0,QG0,PD0,QD0)
-PD0Array=repmat(PD0,1,this.sampleNum);
-QD0Array=repmat(QD0,1,this.sampleNum);
-PDMaxDev=max(abs(this.PDArray-PD0Array),[],2)
-QDMaxDev=max(abs(this.QDArray-QD0Array),[],2)
-PG0Array=repmat()
-end
diff --git a/@ForThesis/StatDeviation.asv b/@ForThesis/StatDeviation.asv
deleted file mode 100644
index 7687bc2..0000000
--- a/@ForThesis/StatDeviation.asv
+++ /dev/null
@@ -1,15 +0,0 @@
-function [ output_args ] = StatDeviation( this,PG0,QG0,PD0,QD0 )%统计误差
-%STATDEVIATION Summary of this function goes here
-%   Detailed explanation goes here
-
-PD0Array=repmat(PD0,1,this.sampleNum);
-QD0Array=repmat(QD0,1,this.sampleNum);
-PDMaxDev=max(abs(this.PDArray-PD0Array),[],2);
-QDMaxDev=max(abs(this.QDArray-QD0Array),[],2);
-PG0Array=repmat(PG0,this.sampleNum,1);
-QG0Array=repmat(QG0,this.sampleNum,1);
-PGMaxDev=max(abs(PG0Array));
-QGMaxDev=max(abs(QG0Array));
-
-end
-
diff --git a/@ForThesis/AddPDQDPGQG.m b/@ForThesis1/AddPDQDPGQG.m
similarity index 100%
rename from @ForThesis/AddPDQDPGQG.m
rename to @ForThesis1/AddPDQDPGQG.m
diff --git a/@ForThesis/ForThesis.m b/@ForThesis1/ForThesis.m
similarity index 100%
rename from @ForThesis/ForThesis.m
rename to @ForThesis1/ForThesis.m
diff --git a/@ForThesis/MaxBranchDeviation.m b/@ForThesis1/MaxBranchDeviation.m
similarity index 100%
rename from @ForThesis/MaxBranchDeviation.m
rename to @ForThesis1/MaxBranchDeviation.m
diff --git a/@ForThesis/MaxDeviation.m b/@ForThesis1/MaxDeviation.m
similarity index 100%
rename from @ForThesis/MaxDeviation.m
rename to @ForThesis1/MaxDeviation.m
diff --git a/@ForThesis/MeanPD.m b/@ForThesis1/MeanPD.m
similarity index 100%
rename from @ForThesis/MeanPD.m
rename to @ForThesis1/MeanPD.m
diff --git a/@ForThesis/MeanPG.m b/@ForThesis1/MeanPG.m
similarity index 100%
rename from @ForThesis/MeanPG.m
rename to @ForThesis1/MeanPG.m
diff --git a/@ForThesis/MeanQD.m b/@ForThesis1/MeanQD.m
similarity index 100%
rename from @ForThesis/MeanQD.m
rename to @ForThesis1/MeanQD.m
diff --git a/@ForThesis/MeanQG.m b/@ForThesis1/MeanQG.m
similarity index 100%
rename from @ForThesis/MeanQG.m
rename to @ForThesis1/MeanQG.m
diff --git a/@ForThesis/PercentOfPass.m b/@ForThesis1/PercentOfPass.m
similarity index 100%
rename from @ForThesis/PercentOfPass.m
rename to @ForThesis1/PercentOfPass.m
diff --git a/@ForThesis/SqareDeviation.m b/@ForThesis1/SqareDeviation.m
similarity index 100%
rename from @ForThesis/SqareDeviation.m
rename to @ForThesis1/SqareDeviation.m
diff --git a/@ForThesis/StatBranchDeviation.m b/@ForThesis1/StatBranchDeviation.m
similarity index 100%
rename from @ForThesis/StatBranchDeviation.m
rename to @ForThesis1/StatBranchDeviation.m
diff --git a/@ForThesis/StatDeviation.m b/@ForThesis1/StatDeviation.m
similarity index 100%
rename from @ForThesis/StatDeviation.m
rename to @ForThesis1/StatDeviation.m
diff --git a/LoadProfile.m b/LoadProfile.m
new file mode 100644
index 0000000..d264e5a
--- /dev/null
+++ b/LoadProfile.m
@@ -0,0 +1,18 @@
+txt=dlmread('E:\书籍\论文\配电网状态估计\负荷.txt');
+txt=txt./10000;
+x=1:length(txt);
+plot(x,txt,'k');
+hold on
+plot(x,txt*0.8,'k--');
+plot(x,txt*1.2,'k--');
+legend('Load profile','Lower/Upper bound')
+
+ylabel('Load /kW')
+xlabel('Time /h')
+k=1;
+for i=0:length(x)-1
+    x1(k)=x(length(x)-i);
+    Z2(k)=txt(length(x)-i)*1.2;
+    k=k+1;
+end
+fill([x x1],[txt'*0.8 Z2],[240, 240, 240]./255,'edgealpha',0)
diff --git a/OPF.m b/OPF.m
index 1f530cc..b1f8dab 100644
--- a/OPF.m
+++ b/OPF.m
@@ -1,8 +1,11 @@
 clc
 clear
+close all
+arrayA=zeros(21,10);
+for I=1: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:\算例\17\17.csv');
+    pf('E:\算例\feeder33\feeder33.txt');
 sigma=0.01;
 RealPD=PD;
 RealQD=QD;
@@ -13,59 +16,277 @@ QD0=sparse(Busnum,1);
 PD0(Loadi)=RealPD(Loadi).*(1+normrnd(0,sigma,length(Loadi),1));
 QD0(Loadi)=RealQD(Loadi).*(1+normrnd(0,sigma,length(Loadi),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;
+
+% mVolt(19)=rVolt(19)*(1+sigma*4);
 %% 画Case A的图
-figV=figure();
-[Busnum,Loadi,Volt,PD,QD,rVolt,UAngel,RealPD,RealQD,rUAngel,Vbi,PDbi,QDbi]=subOPF([],PD0,QD0,mVolt,sigma);%全部有
-subplot(4,1,1);
-plot(1:Busnum,Volt-rVolt,'b.:','Marker','diamond');
+figure('Color',[1 1 1]);
+[Busnum,Loadi,Volt,PD,QD,rVolt,UAngel,RealPD,RealQD,rUAngel,Vbi,PDbi,QDbi,plotGapA]=subOPF([],PD0,QD0,mVolt,sigma);%全部有
+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=CaseAREV(2:end)*100;
+plot(1:length(CaseAREV),CaseAREV,'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,2);
-plot(1:Busnum,UAngel-rUAngel,'b:','Marker','diamond');
+CaseAREA=(UAngel-rUAngel)*100;%Relative Error of Angle in Case A
+% CaseAREA=CaseAREA(2:end)*100;
+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);
-plot(1:Busnum,(PD-RealPD)./(RealPD+0.00001),'b:','Marker','diamond');
+CaseAREP=(PD-RealPD)./(RealPD+0.00001)*100;%Relative Error of PD in Case A
+% CaseAREP=CaseAREP(2:end)*100;
+plot(1:length(CaseAREP),CaseAREP,'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,4);
-plot(1:Busnum,(QD-RealQD)./(RealQD+0.00001),'b:','Marker','diamond');
+CaseAREQ=(QD-RealQD)./(RealQD+0.00001)*100;%Relative Error of QD in Case A
+% CaseAREQ=CaseAREQ(2:end)*100;
+plot(1:length(CaseAREQ),CaseAREQ,'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('误差%');
+%计算Case A的误差
+CaseAE=sqrt((sum(CaseAREV.^2)+sum(CaseAREA.^2)+sum(CaseAREP.^2)+sum(CaseAREQ.^2))/132);
+objA=full(sum(Vbi)+sum(PDbi)+sum(QDbi));
+
+% arrayA(1:19,I)=Vbi;
+% arrayA(21,I)=CaseAE*1000;
+
+end
 %% Case B
-[Busnum,Loadi,Volt,PD,QD,rVolt,UAngel,RealPD,RealQD,rUAngel,Vbi,PDbi,QDbi]=subOPF(cat(2,[1,2,3,4,5,11,12,13,16],[18,19]),PD0,QD0,mVolt,sigma);%
-% [Busnum,Loadi,Volt,PD,QD,rVolt,UAngel,RealPD,RealQD,rUAngel,Vbi,PDbi,QDbi]=subOPF(cat(2,[1,2,3,4,6,7],[]),PD0,QD0,mVolt,sigma);%
+[Busnum,Loadi,Volt,PD,QD,rVolt,UAngel,RealPD,RealQD,rUAngel,Vbi,PDbi,QDbi,plotGapB]=subOPF(setdiff(1:Busnum,[18,21,22,29]),PD0,QD0,mVolt,sigma);%
 subplot(4,1,1);
 hold on;
-plot(1:Busnum,Volt-rVolt,'g.:','Marker','square');
+CaseBREV=(Volt-rVolt)*100;%Relative Error of Voltage in Case B
+% CaseBREV=CaseBREV(2:end)*100;
+plot(1:length(CaseBREV),CaseBREV,'b.:','Marker','square');
 subplot(4,1,2);
 hold on;
-plot(1:Busnum,UAngel-rUAngel,'g:','Marker','square');
+CaseBREA=(UAngel-rUAngel)*100;%Relative Error of Angle in Case B
+% CaseBREA=CaseBREA(2:end)*100;
+plot(1:length(CaseBREA),CaseBREA,'b:','Marker','square');
 subplot(4,1,3);
 hold on;
-plot(1:Busnum,(PD-RealPD)./(RealPD+0.00001),'g:','Marker','square');
+CaseBREP=(PD-RealPD)./(RealPD+0.00001)*100;%Relative Error of PD in Case B
+% CaseBREP=CaseBREP(2:end)*100;
+plot(1:length(CaseBREP),CaseBREP,'b:','Marker','square');
 subplot(4,1,4);
 hold on;
-plot(1:Busnum,(QD-RealQD)./(RealQD+0.00001),'g:','Marker','square');
+CaseBREQ=(QD-RealQD)./(RealQD+0.00001)*100;%Relative Error of QD in Case B
+% CaseBREQ=CaseBREQ(2:end)*100;
+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));
 %% Case C
-[Busnum,Loadi,Volt,PD,QD,rVolt,UAngel,RealPD,RealQD,rUAngel,Vbi,PDbi,QDbi]=subOPF([1:17],PD0,QD0,mVolt,sigma);%
+[Busnum,Loadi,Volt,PD,QD,rVolt,UAngel,RealPD,RealQD,rUAngel,Vbi,PDbi,QDbi,plotGapC]=subOPF([1:17],PD0,QD0,mVolt,sigma);%
 subplot(4,1,1);
 hold on;
-plot(1:Busnum,Volt-rVolt,'r.:','Marker','o');
+CaseCREV=(Volt-rVolt)*100;%Relative Error of Voltage in Case C
+% CaseCREV=CaseCREV(2:end)*100;
+plot(1:length(CaseCREV),CaseCREV,'r.:','Marker','o');
 subplot(4,1,2);
 hold on;
-plot(1:Busnum,UAngel-rUAngel,'r:','Marker','o');
+CaseCREA=(UAngel-rUAngel)*100;%Relative Error of Angle in Case C
+% CaseCREA=CaseCREA(2:end)*100;
+plot(1:length(CaseCREA),CaseCREA,'r:','Marker','o');
 subplot(4,1,3);
 hold on;
-plot(1:Busnum,(PD-RealPD)./(RealPD+0.00001),'r:','Marker','o');
+CaseCREP=(PD-RealPD)./(RealPD+0.00001)*100;%Relative Error of PD in Case C
+% CaseCREP=CaseCREP(2:end)*100;
+plot(1:length(CaseCREP),CaseCREP,'r:','Marker','o');
 subplot(4,1,4);
 hold on;
-plot(1:Busnum,(QD-RealQD)./(RealQD+0.00001),'r:','Marker','o');
+CaseCREQ=(QD-RealQD)./(RealQD+0.00001)*100;%Relative Error of QD in Case C
+% CaseCREQ=CaseCREQ(2:end)*100;
+plot(1:length(CaseCREQ),CaseCREQ,'r:','Marker','o');
 % 画legend
 subplot(4,1,1);
-title('Voltage');
-legend('Case A','Case B','Case C');
+% title('Voltage');
+ld=legend('算例A','算例B','算例C');
+set(ld,'Position',[0.847865087908145 0.786094477711244 0.0543595263724435 0.0605455755156354]);
 subplot(4,1,2);
-title('Voltage Angle');
-legend('Case A','Case B','Case C');
+% title('Voltage Angle');
+ld=legend('算例A','算例B','算例C');
+set(ld,'Position',[0.847865087908145 0.586094477711244 0.0543595263724435 0.0605455755156354]);
 subplot(4,1,3);
-title('Active load power');
-legend('Case A','Case B','Case C');
+% title('Active load power');
+ld=legend('算例A','算例B','算例C');
+set(ld,'Position',[0.847865087908145 0.386094477711244 0.0543595263724435 0.0605455755156354]);
 subplot(4,1,4);
-title('Reactive load  power');
-legend('Case A','Case B','Case C');
-obj=sum(Vbi)+sum(PDbi)+sum(QDbi);
-fprintf('目标函数值 %.2f\n',full(obj));
-
+% title('Reactive load  power');
+ld=legend('算例A','算例B','算例C');
+set(ld,'Position',[0.847865087908145 0.186094477711244 0.0543595263724435 0.0605455755156354]);
+CaseCE=sqrt((sum(CaseCREV.^2)+sum(CaseCREA.^2)+sum(CaseCREP.^2)+sum(CaseCREQ.^2))/132);
+objC=full(sum(Vbi)+sum(PDbi)+sum(QDbi));
+% fprintf('目标函数值 %.2f\n',full(obj));
+fprintf('Case A    Case B    Case C 的误差\n')
+fprintf('%f    %f    %f    \n',CaseAE,CaseBE,CaseCE);
+fprintf('三个Case目标值\n')
+fprintf('%f\t%f\t%f \n',objA,objB,objC)
+%% 画测量值
+% subplot(4,1,1);
+% plot(1:Busnum,mVolt-rVolt,'k.:','Marker','pentagram')
+% subplot(4,1,3);
+% plot(1:Busnum,(PD0-RealPD)./(RealPD+0.00001),'k:','Marker','pentagram');
+% subplot(4,1,4);
+% plot(1:Busnum,(QD0-RealQD)./(RealQD+0.00001),'k:','Marker','pentagram');
+%% 直方图
+% 电压
+figure('Name','电压直方图')
+split_number=20;
+%Case A
+subplot(1,3,1)
+y=CaseAREV;
+ymin=min(y);
+ymax=max(y);
+x=linspace(ymin,ymax,split_number); %将最大最小区间分成split_number个等分点(19等分),然后分别计算各个区间的个数
+yy=hist(y,x); %计算各个区间的个数
+% yy=yy/(sum(yy)*(ymax-ymin)/split_number); %计算各个区间的个数,除以总面积，总面积计算的方式为：所有小分割的面积的和即:
+bar(x,yy) %画出概率密度分布图
+xlabel('Error');
+ylabel('Number of buses');
+title('算例A');
+% ylim([0 4])
+%Case B
+subplot(1,3,2)
+y=CaseBREV;
+ymin=min(y);
+ymax=max(y);
+x=linspace(ymin,ymax,split_number); %将最大最小区间分成split_number个等分点(19等分),然后分别计算各个区间的个数
+yy=hist(y,x); %计算各个区间的个数
+% yy=yy/(sum(yy)*(ymax-ymin)/split_number); %计算各个区间的个数,除以总面积，总面积计算的方式为：所有小分割的面积的和即:
+bar(x,yy) %画出概率密度分布图
+xlabel('Error');
+ylabel('Number of buses');
+title('算例D');
+% ylim([0 4])
+%Case C
+subplot(1,3,3)
+y=CaseCREV;
+ymin=min(y);
+ymax=max(y);
+x=linspace(ymin,ymax,split_number); %将最大最小区间分成split_number个等分点(19等分),然后分别计算各个区间的个数
+yy=hist(y,x); %计算各个区间的个数
+% yy=yy/(sum(yy)*(ymax-ymin)/split_number); %计算各个区间的个数,除以总面积，总面积计算的方式为：所有小分割的面积的和即:
+bar(x,yy) %画出概率密度分布图
+xlabel('Error');
+ylabel('Number of buses');
+title('算例C');
+% ylim([0 4])
+figure('Name','相角直方图')
+% 相角
+split_number=20;
+%Case A
+subplot(2,2,1)
+y=CaseAREA;
+ymin=min(y);
+ymax=max(y);
+x=linspace(ymin,ymax,split_number); %将最大最小区间分成split_number个等分点(19等分),然后分别计算各个区间的个数
+yy=hist(y,x); %计算各个区间的个数
+% yy=yy/(sum(yy)*(ymax-ymin)/split_number); %计算各个区间的个数,除以总面积，总面积计算的方式为：所有小分割的面积的和即:
+bar(x,yy) %画出概率密度分布图
+% ylim([0 4])
+%Case B
+subplot(2,2,2)
+y=CaseBREA;
+ymin=min(y);
+ymax=max(y);
+x=linspace(ymin,ymax,split_number); %将最大最小区间分成split_number个等分点(19等分),然后分别计算各个区间的个数
+yy=hist(y,x); %计算各个区间的个数
+% yy=yy/(sum(yy)*(ymax-ymin)/split_number); %计算各个区间的个数,除以总面积，总面积计算的方式为：所有小分割的面积的和即:
+bar(x,yy) %画出概率密度分布图
+% ylim([0 4])
+%Case C
+subplot(2,2,3)
+y=CaseCREA;
+ymin=min(y);
+ymax=max(y);
+x=linspace(ymin,ymax,split_number); %将最大最小区间分成split_number个等分点(19等分),然后分别计算各个区间的个数
+yy=hist(y,x); %计算各个区间的个数
+% yy=yy/(sum(yy)*(ymax-ymin)/split_number); %计算各个区间的个数,除以总面积，总面积计算的方式为：所有小分割的面积的和即:
+bar(x,yy) %画出概率密度分布图
+%PD
+figure('Name','有功直方图')
+split_number=20;
+%Case A
+subplot(2,2,1)
+y=CaseAREP;
+ymin=min(y);
+ymax=max(y);
+x=linspace(ymin,ymax,split_number); %将最大最小区间分成split_number个等分点(19等分),然后分别计算各个区间的个数
+yy=hist(y,x); %计算各个区间的个数
+% yy=yy/(sum(yy)*(ymax-ymin)/split_number); %计算各个区间的个数,除以总面积，总面积计算的方式为：所有小分割的面积的和即:
+bar(x,yy) %画出概率密度分布图
+% ylim([0 4])
+%Case B
+subplot(2,2,2)
+y=CaseBREP;
+ymin=min(y);
+ymax=max(y);
+x=linspace(ymin,ymax,split_number); %将最大最小区间分成split_number个等分点(19等分),然后分别计算各个区间的个数
+yy=hist(y,x); %计算各个区间的个数
+% yy=yy/(sum(yy)*(ymax-ymin)/split_number); %计算各个区间的个数,除以总面积，总面积计算的方式为：所有小分割的面积的和即:
+bar(x,yy) %画出概率密度分布图
+% ylim([0 4])
+%Case C
+subplot(2,2,3)
+y=CaseCREP;
+ymin=min(y);
+ymax=max(y);
+x=linspace(ymin,ymax,split_number); %将最大最小区间分成split_number个等分点(19等分),然后分别计算各个区间的个数
+yy=hist(y,x); %计算各个区间的个数
+% yy=yy/(sum(yy)*(ymax-ymin)/split_number); %计算各个区间的个数,除以总面积，总面积计算的方式为：所有小分割的面积的和即:
+bar(x,yy) %画出概率密度分布图
+%QD
+figure('Name','无功直方图')
+split_number=20;
+%Case A
+subplot(2,2,1)
+y=CaseAREQ;
+ymin=min(y);
+ymax=max(y);
+x=linspace(ymin,ymax,split_number); %将最大最小区间分成split_number个等分点(19等分),然后分别计算各个区间的个数
+yy=hist(y,x); %计算各个区间的个数
+% yy=yy/(sum(yy)*(ymax-ymin)/split_number); %计算各个区间的个数,除以总面积，总面积计算的方式为：所有小分割的面积的和即:
+bar(x,yy) %画出概率密度分布图
+% ylim([0 4])
+%Case B
+subplot(2,2,2)
+y=CaseBREQ;
+ymin=min(y);
+ymax=max(y);
+x=linspace(ymin,ymax,split_number); %将最大最小区间分成split_number个等分点(19等分),然后分别计算各个区间的个数
+yy=hist(y,x); %计算各个区间的个数
+% yy=yy/(sum(yy)*(ymax-ymin)/split_number); %计算各个区间的个数,除以总面积，总面积计算的方式为：所有小分割的面积的和即:
+bar(x,yy) %画出概率密度分布图
+% ylim([0 4])
+%Case C
+subplot(2,2,3)
+y=CaseCREQ;
+ymin=min(y);
+ymax=max(y);
+x=linspace(ymin,ymax,split_number); %将最大最小区间分成split_number个等分点(19等分),然后分别计算各个区间的个数
+yy=hist(y,x); %计算各个区间的个数
+% yy=yy/(sum(yy)*(ymax-ymin)/split_number); %计算各个区间的个数,除以总面积，总面积计算的方式为：所有小分割的面积的和即:
+bar(x,yy) %画出概率密度分布图
+%画收敛曲线
+fz=find(abs(plotGapA)==0);
+fz=fz(1);
+figure('Name','互补曲线')
+plot(1:fz-1,plotGapA(1:fz-1));
diff --git a/PD0.mat b/PD0.mat
new file mode 100644
index 0000000..9caa533
Binary files /dev/null and b/PD0.mat differ
diff --git a/QD0.mat b/QD0.mat
new file mode 100644
index 0000000..bb89563
Binary files /dev/null and b/QD0.mat differ
diff --git a/figure1.fig b/figure1.fig
new file mode 100644
index 0000000..7b274fb
Binary files /dev/null and b/figure1.fig differ
diff --git a/figure2.fig b/figure2.fig
new file mode 100644
index 0000000..f2a0cb4
Binary files /dev/null and b/figure2.fig differ
diff --git a/histA.fig b/histA.fig
new file mode 100644
index 0000000..75089ff
Binary files /dev/null and b/histA.fig differ
diff --git a/histP.fig b/histP.fig
new file mode 100644
index 0000000..00019b1
Binary files /dev/null and b/histP.fig differ
diff --git a/histQ.fig b/histQ.fig
new file mode 100644
index 0000000..c926b7e
Binary files /dev/null and b/histQ.fig differ
diff --git a/histV.fig b/histV.fig
new file mode 100644
index 0000000..7408880
Binary files /dev/null and b/histV.fig differ
diff --git a/hs_err_pid4476.log b/hs_err_pid4476.log
new file mode 100644
index 0000000..80d257d
--- /dev/null
+++ b/hs_err_pid4476.log
@@ -0,0 +1,258 @@
+#
+# An unexpected error has been detected by Java Runtime Environment:
+#
+#  EXCEPTION_ACCESS_VIOLATION (0xc0000005) at pc=0x77cb80dc, pid=4476, tid=6716
+#
+# Java VM: Java HotSpot(TM) Client VM (11.2-b01 mixed mode windows-x86)
+# Problematic frame:
+# C  [ntdll.dll+0x580dc]
+#
+# If you would like to submit a bug report, please visit:
+#   http://java.sun.com/webapps/bugreport/crash.jsp
+# The crash happened outside the Java Virtual Machine in native code.
+# See problematic frame for where to report the bug.
+#
+
+---------------  T H R E A D  ---------------
+
+Current thread (0x09b9c800):  JavaThread "main" [_thread_in_native, id=6716, stack(0x00430000,0x00c30000)]
+
+siginfo: ExceptionCode=0xc0000005, reading address 0xffffffff
+
+Registers:
+EAX=0x00000000, EBX=0x00000000, ECX=0x00000008, EDX=0xfffffff8
+ESP=0x00c29b0c, EBP=0x00c29b14, ESI=0x09b90000, EDI=0x79cd34e4
+EIP=0x77cb80dc, EFLAGS=0x00010246
+
+Top of Stack: (sp=0x00c29b0c)
+0x00c29b0c:   00c29c14 79d135a4 00c29b60 7bfc20d6
+0x00c29b1c:   09b90000 00000000 00000000 79cd34e4
+0x00c29b2c:   79d135a4 00c29c14 00000000 056b4b88
+0x00c29b3c:   000007fd 00c29bc0 00c29bc4 00c29b28
+0x00c29b4c:   00c2954c 00c29b8c 7bfc240d 7bffa368
+0x00c29b5c:   ffffffff 00c29b9c 7bfcc0b4 00000000
+0x00c29b6c:   79cd34e4 00000003 00c29c14 002b0023
+0x00c29b7c:   0053002b 002b002b 00c29b6c 00c2954c 
+
+Instructions: (pc=0x77cb80dc)
+0x77cb80cc:   02 00 8b 45 10 a8 07 0f 85 41 f7 02 00 8d 50 f8
+0x77cb80dc:   80 7a 07 05 0f 84 1d f7 02 00 f6 42 07 3f 0f 84 
+
+
+Stack: [0x00430000,0x00c30000],  sp=0x00c29b0c,  free space=8166k
+Native frames: (J=compiled Java code, j=interpreted, Vv=VM code, C=native code)
+C  [ntdll.dll+0x580dc]
+C  [MSVCR71.dll+0x20d6]
+C  [MSVCR71.dll+0xc0b4]
+V  [jvm.dll+0x1c8b74]
+
+Java frames: (J=compiled Java code, j=interpreted, Vv=VM code)
+j  com.mathworks.jmi.NativeMatlab.SendMatlabMessage(Ljava/lang/Object;)Ljava/lang/Object;+0
+j  com.mathworks.jmi.NativeMatlab.sendMatlabMessage(Ljava/lang/Object;)Ljava/lang/Object;+22
+j  com.mathworks.jmi.MatlabLooper.sendMatlabMessage(Lcom/mathworks/services/message/MWMessage;)Ljava/lang/Object;+20
+j  com.mathworks.jmi.Matlab.mtFevalConsoleOutput(Ljava/lang/String;[Ljava/lang/Object;I)Ljava/lang/Object;+58
+j  com.mathworks.mde.desk.MLDesktop$9.run()V+14
+j  com.mathworks.jmi.NativeMatlab.dispatchMTRequests(Z)V+50
+v  ~StubRoutines::call_stub
+
+---------------  P R O C E S S  ---------------
+
+Java Threads: ( => current thread )
+  0x22e17400 JavaThread "Inactive RequestProcessor thread [Was:TimedSoftReference/org.openide.util.TimedSoftReference]" daemon [_thread_blocked, id=1756, stack(0x05f20000,0x05fa0000)]
+  0x25229400 JavaThread "Thread-1175" [_thread_in_native, id=6188, stack(0x2fd00000,0x30500000)]
+  0x250c0000 JavaThread "pool-2-thread-1" [_thread_blocked, id=8112, stack(0x05bc0000,0x05c40000)]
+  0x231b5800 JavaThread "Timer-120" [_thread_blocked, id=8128, stack(0x04710000,0x04790000)]
+  0x22ff0800 JavaThread "Active Reference Queue Daemon" daemon [_thread_blocked, id=7892, stack(0x1eee0000,0x1ef60000)]
+  0x22ff3c00 JavaThread "Timer-4" daemon [_thread_blocked, id=5616, stack(0x1ebc0000,0x1ec40000)]
+  0x22ff3800 JavaThread "Timer-3" daemon [_thread_blocked, id=4204, stack(0x1eaf0000,0x1eb70000)]
+  0x22ff2c00 JavaThread "Prefs Updater" daemon [_thread_blocked, id=7656, stack(0x24300000,0x24380000)]
+  0x22ff3000 JavaThread "Timer-1" [_thread_blocked, id=2792, stack(0x23fe0000,0x24060000)]
+  0x22c6c400 JavaThread "TimerQueue" daemon [_thread_blocked, id=6316, stack(0x23df0000,0x23e70000)]
+  0x22c66c00 JavaThread "AWT-EventQueue-0" [_thread_blocked, id=2764, stack(0x23d20000,0x23da0000)]
+  0x22b78000 JavaThread "AWT-Windows" daemon [_thread_in_native, id=7904, stack(0x23300000,0x23380000)]
+  0x22bffc00 JavaThread "AWT-Shutdown" [_thread_blocked, id=7840, stack(0x23280000,0x23300000)]
+  0x22bff800 JavaThread "Java2D Disposer" daemon [_thread_blocked, id=8164, stack(0x23200000,0x23280000)]
+  0x1c7c9800 JavaThread "Timer-0" [_thread_blocked, id=7820, stack(0x1c990000,0x1ca10000)]
+  0x1c793800 JavaThread "JMI Unnamed Thread" [_thread_in_native, id=7988, stack(0x061c0000,0x069c0000)]
+  0x09e2c800 JavaThread "Low Memory Detector" daemon [_thread_blocked, id=7900, stack(0x1c680000,0x1c700000)]
+  0x09e2b800 JavaThread "CompilerThread0" daemon [_thread_blocked, id=8152, stack(0x1f950000,0x20150000)]
+  0x09e0d400 JavaThread "Attach Listener" daemon [_thread_blocked, id=7876, stack(0x1c5c0000,0x1c640000)]
+  0x09df8400 JavaThread "Finalizer" daemon [_thread_blocked, id=8136, stack(0x1e3f0000,0x1e470000)]
+  0x09df6c00 JavaThread "Reference Handler" daemon [_thread_blocked, id=7848, stack(0x0ff80000,0x10000000)]
+=>0x09b9c800 JavaThread "main" [_thread_in_native, id=6716, stack(0x00430000,0x00c30000)]
+
+Other Threads:
+  0x09df1000 VMThread [stack: 0x1db90000,0x1e390000] [id=3524]
+  0x09e32400 WatcherThread [stack: 0x21e00000,0x22600000] [id=5768]
+
+VM state:not at safepoint (normal execution)
+
+VM Mutex/Monitor currently owned by a thread: None
+
+Heap
+ def new generation   total 29504K, used 23735K [0x10390000, 0x12390000, 0x12390000)
+  eden space 26240K,  83% used [0x10390000, 0x118fcfb8, 0x11d30000)
+  from space 3264K,  55% used [0x12060000, 0x12220f70, 0x12390000)
+  to   space 3264K,   0% used [0x11d30000, 0x11d30000, 0x12060000)
+ tenured generation   total 98304K, used 69305K [0x12390000, 0x18390000, 0x18390000)
+   the space 98304K,  70% used [0x12390000, 0x1673e6f8, 0x1673e800, 0x18390000)
+ compacting perm gen  total 39680K, used 39594K [0x18390000, 0x1aa50000, 0x1c390000)
+   the space 39680K,  99% used [0x18390000, 0x1aa3a830, 0x1aa3aa00, 0x1aa50000)
+No shared spaces configured.
+
+Dynamic libraries:
+0x00400000 - 0x00425000 	C:\Program Files (x86)\MATLAB\R2010a\bin\win32\MATLAB.exe
+0x77c60000 - 0x77db7000 	C:\Windows\SYSTEM32\ntdll.dll
+0x75bf0000 - 0x75d20000 	C:\Windows\SYSTEM32\KERNEL32.DLL
+0x76360000 - 0x76406000 	C:\Windows\SYSTEM32\KERNELBASE.dll
+0x74720000 - 0x747c7000 	C:\Windows\system32\apphelp.dll
+0x7b940000 - 0x7b9ec000 	C:\Program Files (x86)\MATLAB\R2010a\bin\win32\libut.dll
+0x78770000 - 0x787ee000 	C:\Program Files (x86)\MATLAB\R2010a\bin\win32\libmwfl.dll
+0x7ba50000 - 0x7bc63000 	C:\Program Files (x86)\MATLAB\R2010a\bin\win32\libmwservices.dll
+0x7a0e0000 - 0x7a144000 	C:\Program Files (x86)\MATLAB\R2010a\bin\win32\mcr.dll
+0x00140000 - 0x001c7000 	C:\Windows\WinSxS\x86_microsoft.vc80.crt_1fc8b3b9a1e18e3b_8.0.50727.6910_none_d089c358442de345\MSVCP80.dll
+0x001d0000 - 0x0026b000 	C:\Windows\WinSxS\x86_microsoft.vc80.crt_1fc8b3b9a1e18e3b_8.0.50727.6910_none_d089c358442de345\MSVCR80.dll
+0x76900000 - 0x76912000 	C:\Windows\SYSTEM32\imagehlp.dll
+0x75d20000 - 0x75d26000 	C:\Windows\SYSTEM32\PSAPI.DLL
+0x7bfa0000 - 0x7bfbf000 	C:\Program Files (x86)\MATLAB\R2010a\bin\win32\LIBEXPAT.dll
+0x7bd50000 - 0x7be38000 	C:\Program Files (x86)\MATLAB\R2010a\bin\win32\icuuc40.dll
+0x7b4c0000 - 0x7b4cc000 	C:\Program Files (x86)\MATLAB\R2010a\bin\win32\icuio40.dll
+0x76490000 - 0x765a6000 	C:\Windows\SYSTEM32\USER32.dll
+0x76940000 - 0x769ee000 	C:\Windows\SYSTEM32\ADVAPI32.dll
+0x7d160000 - 0x7d16e000 	C:\Program Files (x86)\MATLAB\R2010a\bin\win32\boost_date_time-vc80-mt-1_36.dll
+0x7d330000 - 0x7d346000 	C:\Program Files (x86)\MATLAB\R2010a\bin\win32\boost_filesystem-vc80-mt-1_36.dll
+0x7d050000 - 0x7d060000 	C:\Program Files (x86)\MATLAB\R2010a\bin\win32\boost_signals-vc80-mt-1_36.dll
+0x7c500000 - 0x7c507000 	C:\Program Files (x86)\MATLAB\R2010a\bin\win32\boost_system-vc80-mt-1_36.dll
+0x7b4a0000 - 0x7b4ac000 	C:\Program Files (x86)\MATLAB\R2010a\bin\win32\boost_thread-vc80-mt-1_36.dll
+0x7b9f0000 - 0x7ba4d000 	C:\Program Files (x86)\MATLAB\R2010a\bin\win32\libmx.dll
+0x7bc70000 - 0x7bd47000 	C:\Program Files (x86)\MATLAB\R2010a\bin\win32\libmwmathutil.dll
+0x79e30000 - 0x79e6a000 	C:\Program Files (x86)\MATLAB\R2010a\bin\win32\mpath.dll
+0x7c0a0000 - 0x7c166000 	C:\Program Files (x86)\MATLAB\R2010a\bin\win32\mlutil.dll
+0x78720000 - 0x7874e000 	C:\Program Files (x86)\MATLAB\R2010a\bin\win32\tbb.dll
+0x73d00000 - 0x73d87000 	C:\Windows\WinSxS\x86_microsoft.windows.common-controls_6595b64144ccf1df_5.82.9200.16658_none_bf1359a245f1cd12\COMCTL32.dll
+0x76870000 - 0x768f9000 	C:\Windows\SYSTEM32\comdlg32.dll
+0x75030000 - 0x75042000 	C:\Windows\SYSTEM32\NETAPI32.dll
+0x753a0000 - 0x753f0000 	C:\Windows\SYSTEM32\WS2_32.dll
+0x76ad0000 - 0x77b96000 	C:\Windows\SYSTEM32\SHELL32.dll
+0x7e890000 - 0x7e8a9000 	C:\Program Files (x86)\MATLAB\R2010a\bin\win32\libmwi18n.dll
+0x7b4d0000 - 0x7b539000 	C:\Program Files (x86)\MATLAB\R2010a\bin\win32\iqm.dll
+0x7e850000 - 0x7e871000 	C:\Program Files (x86)\MATLAB\R2010a\bin\win32\libmwbridge.dll
+0x7b550000 - 0x7b561000 	C:\Program Files (x86)\MATLAB\R2010a\bin\win32\libmex.dll
+0x79fb0000 - 0x7a017000 	C:\Program Files (x86)\MATLAB\R2010a\bin\win32\m_dispatcher.dll
+0x787f0000 - 0x7894e000 	C:\Program Files (x86)\MATLAB\R2010a\bin\win32\mcos.dll
+0x7be40000 - 0x7bf3f000 	C:\Program Files (x86)\MATLAB\R2010a\bin\win32\libmwgui.dll
+0x78cf0000 - 0x7914f000 	C:\Program Files (x86)\MATLAB\R2010a\bin\win32\hg.dll
+0x785c0000 - 0x78628000 	C:\Program Files (x86)\MATLAB\R2010a\bin\win32\hgutils.dll
+0x7a6c0000 - 0x7ab7b000 	C:\Program Files (x86)\MATLAB\R2010a\bin\win32\m_interpreter.dll
+0x79d40000 - 0x79d96000 	C:\Program Files (x86)\MATLAB\R2010a\bin\win32\jmi.dll
+0x7a150000 - 0x7a1cc000 	C:\Program Files (x86)\MATLAB\R2010a\bin\win32\udd_mi.dll
+0x7b210000 - 0x7b492000 	C:\Program Files (x86)\MATLAB\R2010a\bin\win32\uiw.dll
+0x78c20000 - 0x78c37000 	C:\Program Files (x86)\MATLAB\R2010a\bin\win32\mwoles05.DLL
+0x79e70000 - 0x79eb9000 	C:\Program Files (x86)\MATLAB\R2010a\bin\win32\comcli.dll
+0x7b600000 - 0x7b60b000 	C:\Program Files (x86)\MATLAB\R2010a\bin\win32\mlautoregister.dll
+0x760b0000 - 0x761ad000 	C:\Windows\SYSTEM32\GDI32.dll
+0x75e80000 - 0x75f99000 	C:\Windows\SYSTEM32\ole32.dll
+0x767b0000 - 0x76861000 	C:\Windows\SYSTEM32\msvcrt.dll
+0x7b4b0000 - 0x7b4b4000 	C:\Program Files (x86)\MATLAB\R2010a\bin\win32\icudt40.dll
+0x7b670000 - 0x7b768000 	C:\Program Files (x86)\MATLAB\R2010a\bin\win32\icuin40.dll
+0x75400000 - 0x75434000 	C:\Windows\SYSTEM32\sechost.dll
+0x75d30000 - 0x75ddc000 	C:\Windows\SYSTEM32\RPCRT4.dll
+0x78c90000 - 0x78ca3000 	C:\Program Files (x86)\MATLAB\R2010a\bin\win32\zlib1.dll
+0x7cfb0000 - 0x7d00e000 	C:\Program Files (x86)\MATLAB\R2010a\bin\win32\xmlcore.dll
+0x75800000 - 0x75840000 	C:\Windows\SYSTEM32\SHLWAPI.dll
+0x75020000 - 0x7502b000 	C:\Windows\SYSTEM32\netutils.dll
+0x75000000 - 0x7501c000 	C:\Windows\SYSTEM32\srvcli.dll
+0x74ff0000 - 0x75000000 	C:\Windows\SYSTEM32\wkscli.dll
+0x75de0000 - 0x75de8000 	C:\Windows\SYSTEM32\NSI.dll
+0x756c0000 - 0x757f6000 	C:\Windows\SYSTEM32\combase.dll
+0x7b180000 - 0x7b209000 	C:\Program Files (x86)\MATLAB\R2010a\bin\win32\udd.dll
+0x79720000 - 0x798e2000 	C:\Program Files (x86)\MATLAB\R2010a\bin\win32\xerces-c_2_7.dll
+0x7d280000 - 0x7d31d000 	C:\Program Files (x86)\MATLAB\R2010a\bin\win32\boost_regex-vc80-mt-1_36.dll
+0x78b60000 - 0x78b92000 	C:\Program Files (x86)\MATLAB\R2010a\bin\win32\libmat.dll
+0x79a40000 - 0x79a70000 	C:\Program Files (x86)\MATLAB\R2010a\bin\win32\libmwhardcopy.dll
+0x79da0000 - 0x79de3000 	C:\Program Files (x86)\MATLAB\R2010a\bin\win32\libuij.dll
+0x78630000 - 0x7871f000 	C:\Program Files (x86)\MATLAB\R2010a\bin\win32\hgdatatypes.dll
+0x78ba0000 - 0x78bd4000 	C:\Program Files (x86)\MATLAB\R2010a\bin\win32\libmwlapack.dll
+0x00cb0000 - 0x00dbf000 	C:\Windows\WinSxS\x86_microsoft.vc80.mfc_1fc8b3b9a1e18e3b_8.0.50727.6195_none_cbf5e994470a1a8f\MFC80.DLL
+0x73470000 - 0x734d0000 	C:\Windows\SYSTEM32\WINSPOOL.DRV
+0x76720000 - 0x767ab000 	C:\Windows\SYSTEM32\OLEAUT32.dll
+0x79ec0000 - 0x79f1b000 	C:\Program Files (x86)\MATLAB\R2010a\bin\win32\profiler.dll
+0x7b5d0000 - 0x7b5f2000 	C:\Program Files (x86)\MATLAB\R2010a\bin\win32\libmwmathrng.dll
+0x78c00000 - 0x78c12000 	C:\Program Files (x86)\MATLAB\R2010a\bin\win32\m_pcodeio.dll
+0x79df0000 - 0x79e27000 	C:\Program Files (x86)\MATLAB\R2010a\bin\win32\m_ir.dll
+0x7a1d0000 - 0x7a6bd000 	C:\Program Files (x86)\MATLAB\R2010a\bin\win32\m_parser.dll
+0x78be0000 - 0x78bfa000 	C:\Program Files (x86)\MATLAB\R2010a\bin\win32\m_pcodegen.dll
+0x7e810000 - 0x7e844000 	C:\Program Files (x86)\MATLAB\R2010a\bin\win32\libmwspmatrix.dll
+0x7b660000 - 0x7b669000 	C:\Program Files (x86)\MATLAB\R2010a\bin\win32\uinone.dll
+0x00f20000 - 0x00f3b000 	C:\Windows\WinSxS\x86_microsoft.vc80.atl_1fc8b3b9a1e18e3b_8.0.50727.6195_none_d1cb102c435421de\ATL80.DLL
+0x75250000 - 0x7526c000 	C:\Windows\SYSTEM32\SspiCli.dll
+0x7b770000 - 0x7b939000 	C:\Program Files (x86)\MATLAB\R2010a\bin\win32\libhdf5.dll
+0x7b580000 - 0x7b58e000 	C:\Program Files (x86)\MATLAB\R2010a\bin\win32\libmwbinder.dll
+0x7b540000 - 0x7b54e000 	C:\Program Files (x86)\MATLAB\R2010a\bin\win32\ir_xfmr.dll
+0x7b610000 - 0x7b61a000 	C:\Program Files (x86)\MATLAB\R2010a\bin\win32\mtok.dll
+0x75240000 - 0x75249000 	C:\Windows\SYSTEM32\CRYPTBASE.dll
+0x75140000 - 0x751b3000 	C:\Windows\SYSTEM32\SHCORE.DLL
+0x751e0000 - 0x75231000 	C:\Windows\SYSTEM32\bcryptPrimitives.dll
+0x76920000 - 0x76940000 	C:\Windows\system32\IMM32.DLL
+0x769f0000 - 0x76acc000 	C:\Windows\SYSTEM32\MSCTF.dll
+0x03a20000 - 0x03a2a000 	C:\Windows\WinSxS\x86_microsoft.vc80.mfcloc_1fc8b3b9a1e18e3b_8.0.50727.6195_none_03ce2c72205943d3\MFC80CHS.DLL
+0x76470000 - 0x76481000 	C:\Windows\SYSTEM32\profapi.dll
+0x73ac0000 - 0x73b48000 	C:\Windows\system32\uxtheme.dll
+0x70cb0000 - 0x70e70000 	C:\Program Files (x86)\360\360safe\safemon\safemon.dll
+0x73dc0000 - 0x73fb7000 	C:\Windows\WinSxS\x86_microsoft.windows.common-controls_6595b64144ccf1df_6.0.9200.16384_none_893961408605e985\COMCTL32.dll
+0x751d0000 - 0x751d8000 	C:\Windows\SYSTEM32\VERSION.dll
+0x72b50000 - 0x72b90000 	C:\Program Files (x86)\TeamViewer\Version9\tv_w32.dll
+0x761b0000 - 0x7635f000 	C:\Windows\SYSTEM32\SETUPAPI.dll
+0x75460000 - 0x754a6000 	C:\Windows\SYSTEM32\CFGMGR32.dll
+0x75440000 - 0x7545e000 	C:\Windows\SYSTEM32\DEVOBJ.dll
+0x75fa0000 - 0x76014000 	C:\Windows\SYSTEM32\clbcatq.dll
+0x747e0000 - 0x748f9000 	C:\Windows\system32\propsys.dll
+0x74f40000 - 0x74f62000 	C:\Windows\SYSTEM32\iphlpapi.dll
+0x74f30000 - 0x74f38000 	C:\Windows\SYSTEM32\WINNSI.DLL
+0x74cd0000 - 0x74d45000 	C:\Windows\SYSTEM32\DNSAPI.dll
+0x74b90000 - 0x74ba0000 	C:\Windows\SYSTEM32\dhcpcsvc6.DLL
+0x74df0000 - 0x74e02000 	C:\Windows\SYSTEM32\dhcpcsvc.DLL
+0x042e0000 - 0x042e3000 	C:\Windows\SYSTEM32\icmp.Dll
+0x79ae0000 - 0x79d36000 	C:\Program Files (x86)\MATLAB\R2010a\sys\java\jre\win32\jre\bin\client\jvm.dll
+0x74e40000 - 0x74e61000 	C:\Windows\SYSTEM32\WINMM.dll
+0x7bfc0000 - 0x7c016000 	C:\Program Files (x86)\MATLAB\R2010a\bin\win32\MSVCR71.dll
+0x74e10000 - 0x74e3a000 	C:\Windows\SYSTEM32\WINMMBASE.dll
+0x749f0000 - 0x74b19000 	C:\Windows\SYSTEM32\dbghelp.dll
+0x79f70000 - 0x79f78000 	C:\Program Files (x86)\MATLAB\R2010a\sys\java\jre\win32\jre\bin\hpi.dll
+0x72930000 - 0x729a1000 	C:\Windows\SYSTEM32\SogouTSF.ime
+0x735b0000 - 0x735b6000 	C:\Windows\SYSTEM32\MSIMG32.dll
+0x73780000 - 0x737a6000 	C:\Windows\SYSTEM32\ntmarta.dll
+0x79f90000 - 0x79f9c000 	C:\Program Files (x86)\MATLAB\R2010a\sys\java\jre\win32\jre\bin\verify.dll
+0x78cb0000 - 0x78ccf000 	C:\Program Files (x86)\MATLAB\R2010a\sys\java\jre\win32\jre\bin\java.dll
+0x79fa0000 - 0x79faf000 	C:\Program Files (x86)\MATLAB\R2010a\sys\java\jre\win32\jre\bin\zip.dll
+0x10000000 - 0x1038d000 	C:\Windows\system32\SogouPy.ime
+0x73340000 - 0x73390000 	C:\Windows\SYSTEM32\OLEACC.dll
+0x72900000 - 0x72923000 	C:\Program Files (x86)\SogouInput\Components\SgAppender\1.0.0.207\SgAppender_Dll.dll
+0x0fc70000 - 0x0fd2e000 	C:\Program Files (x86)\SogouInput\7.1.0.2005\Resource.dll
+
+VM Arguments:
+jvm_args: -Xss512k -XX:PermSize=32m -Xms64m -XX:NewRatio=3 -XX:MaxPermSize=64m -Xmx128m -XX:MaxDirectMemorySize=1200000000 -Dsun.java2d.noddraw=true -Dsun.awt.nopixfmt=true -Xshare:off -Xrs -Djava.library.path=C:\Program Files (x86)\MATLAB\R2010a\bin\win32 vfprintf abort 
+java_command: <unknown>
+Launcher Type: generic
+
+Environment Variables:
+PATH=C:\Program Files (x86)\MATLAB\R2010a\sys\java\jre\win32\jre\bin;C:\Program Files (x86)\MATLAB\R2010a\sys\webrenderer\windows\corecomponents;C:\Program Files (x86)\MATLAB\R2010a\sys\webrenderer\windows;C:\Program Files\Java\jre7\bin;C:\Program Files (x86)\AMD APP\bin\x86_64;C:\Program Files (x86)\AMD APP\bin\x86;C:\Program Files (x86)\Git\bin;c:/go/bin;C:\Python34\;C:\Python34\Scripts;C:\Program Files (x86)\Common Files\Intel\Shared Libraries\redist\ia32\mpirt;C:\Program Files (x86)\Common Files\Intel\Shared Libraries\redist\ia32\compiler;C:\Windows\system32;C:\Windows;C:\Windows\System32\Wbem;C:\Windows\System32\WindowsPowerShell\v1.0\;C:\Program Files (x86)\MATLAB\R2010a\runtime\win32;C:\Program Files (x86)\MATLAB\R2010a\bin;C:\Program Files (x86)\Bitvise SSH Client;C:\CTEX\UserData\miktex\bin;C:\CTEX\MiKTeX\miktex\bin;C:\CTEX\CTeX\ctex\bin;C:\CTEX\CTeX\cct\bin;C:\CTEX\CTeX\ty\bin;C:\CTEX\Ghostscript\gs9.05\bin;C:\CTEX\GSview\gsview;C:\CTEX\WinEdt;C:\Program Files\TortoiseSVN\bin;C:\Program Files (x86)\ATI Technologies\ATI.ACE\Core-Static;
+USERNAME=dmy
+OS=Windows_NT
+PROCESSOR_IDENTIFIER=Intel64 Family 6 Model 42 Stepping 7, GenuineIntel
+
+
+
+---------------  S Y S T E M  ---------------
+
+OS: Windows NT 6.2 Build 9200 
+
+CPU:total 4 (8 cores per cpu, 2 threads per core) family 6 model 10 stepping 7, cmov, cx8, fxsr, mmx, sse, sse2, sse3, ssse3, ht
+
+Memory: 4k page, physical 4194303k(4194303k free), swap 4194303k(4194303k free)
+
+vm_info: Java HotSpot(TM) Client VM (11.2-b01) for windows-x86 JRE (1.6.0_12-b04), built on Jan 17 2009 09:57:14 by "java_re" with MS VC++ 7.1
+
+time: Wed Sep 03 22:06:16 2014
+elapsed time: 6267 seconds
+
diff --git a/mVolt.mat b/mVolt.mat
new file mode 100644
index 0000000..34197ec
Binary files /dev/null and b/mVolt.mat differ
diff --git a/subOPF.m b/subOPF.m
index e13702c..f61e831 100644
--- a/subOPF.m
+++ b/subOPF.m
@@ -1,11 +1,11 @@
-function [Busnum,Loadi,Volt,PD,QD,rVolt,UAngel,RealPD,RealQD,rUAngel,Vbi,PDbi,QDbi]=subOPF(noMeasurei,PD0,QD0,mVolt,sigma)
+function [Busnum,Loadi,Volt,PD,QD,rVolt,UAngel,RealPD,RealQD,rUAngel,Vbi,PDbi,QDbi,plotGap]=subOPF(noMeasurei,PD0,QD0,mVolt,sigma)
 tic
 %% 存在问题
 % 变压器变比的位置没有考虑，由于现在用的变比都是1，所以没有影响。 20130123
 %%
-thesis=ForThesis(1,62);
+% thesis=ForThesis(1,62);
 [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:\算例\17\17.csv');
+    pf('E:\算例\feeder33\feeder33.txt');
 %         pf('E:\算例\柳金Ⅰ926_21671693_2012-09-06\newFIle16.txt');
 %pf('D:\Project\青秀降损项目\最小化潮流\最小潮流算例\原始\津头站津视922（3-1）_0.5_120%.txt');
 %pf('D:\Project\最小化潮流\最小潮流算例\仙海919.txt');
@@ -48,42 +48,49 @@ ContrlCount=size(Loadi,1)*2+Busnum*2+Busnum+length(Loadi)*2;
 kmax=100;
 Precision=Precision/1;
 %% 加误差
-% sigma=0.03;
-% PD0(Loadi)=RealPD(Loadi).*(1+normrnd(0,sigma,length(Loadi),1));
-% QD0(Loadi)=RealQD(Loadi).*(1+normrnd(0,sigma,length(Loadi),1));
-% mVolt=rVolt.*(1+normrnd(0,sigma,length(rVolt),1))';
 %找DG
 DGi=find(PD0<0);
-% PD0(DGi)=PD0(DGi).*(1+normrnd(0,sigma*0.1,length(DGi),1));
-% QD0(DGi)=QD0(DGi).*(1+normrnd(0,sigma*0.1,length(DGi),1));
-% mVolt(DGi)=mVolt(DGi).*(1+normrnd(0,sigma*0.1,length(DGi),1))';
 %
 mPD=PD0;
 mQD=QD0;
-%加估计上下界
-lPD=abs(RealPD*3*sigma);
-uPD=abs(RealPD*3*sigma);
-lQD=abs(RealQD*3*sigma);
-uQD=abs(RealQD*3*sigma);
+%% 加估计上下界,用真实值
+% lPD=abs(RealPD*3*sigma);
+% uPD=abs(RealPD*3*sigma);
+% lQD=abs(RealQD*3*sigma);
+% uQD=abs(RealQD*3*sigma);
+% lVolt=abs(mVolt'*3*sigma);
+% uVolt=abs(mVolt'*3*sigma);
+% %DG
+% lPD(DGi)=abs(RealPD(DGi)*3*sigma*0.1);
+% uPD(DGi)=abs(RealPD(DGi)*3*sigma*0.1);
+% lQD(DGi)=abs(RealQD(DGi)*3*sigma*0.1);
+% uQD(DGi)=abs(RealQD(DGi)*3*sigma*0.1);
+% lVolt(DGi)=abs(mVolt(DGi)'*3*sigma*0.1);
+% uVolt(DGi)=abs(mVolt(DGi)'*3*sigma*0.1);
+%% 加估计上下界,用测量值
+lPD=abs(mPD*3*sigma);
+uPD=abs(mPD*3*sigma);
+lQD=abs(mQD*3*sigma);
+uQD=abs(mQD*3*sigma);
 lVolt=abs(mVolt'*3*sigma);
 uVolt=abs(mVolt'*3*sigma);
 %DG
-lPD(DGi)=abs(RealPD(DGi)*3*sigma*0.1);
-uPD(DGi)=abs(RealPD(DGi)*3*sigma*0.1);
-lQD(DGi)=abs(RealQD(DGi)*3*sigma*0.1);
-uQD(DGi)=abs(RealQD(DGi)*3*sigma*0.1);
-lVolt(DGi)=abs(mVolt(DGi)'*3*sigma*0.1);
-uVolt(DGi)=abs(mVolt(DGi)'*3*sigma*0.1);
+lPD(DGi)=abs(mPD(DGi)*3*sigma*0.3);
+uPD(DGi)=abs(mPD(DGi)*3*sigma*0.3);
+lQD(DGi)=abs(mQD(DGi)*3*sigma*0.3);
+uQD(DGi)=abs(mQD(DGi)*3*sigma*0.3);
+lVolt(DGi)=abs(mVolt(DGi)'*3*sigma*0.3);
+uVolt(DGi)=abs(mVolt(DGi)'*3*sigma*0.3);
 %% 不测量的数据
-lPD(noMeasurei)=0.1*mPD(noMeasurei);
-uPD(noMeasurei)=0.1*mPD(noMeasurei);
-lQD(noMeasurei)=0.1*mQD(noMeasurei);
-uQD(noMeasurei)=0.1*mPD(noMeasurei);
-lVolt(noMeasurei)=0.1*mVolt(noMeasurei);
-uVolt(noMeasurei)=0.1*mVolt(noMeasurei);
+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);
+lVolt(noMeasurei)=0.7*mVolt(noMeasurei);%电压在0.93~1.07
+uVolt(noMeasurei)=0.7*mVolt(noMeasurei);
 %错误数据
 %mVolt(2)=5;
-bigM=10;
+bigM=1;
 Vbi=sparse(0.5*ones(Busnum,1));
 PDbi=sparse(0.5*ones(length(Loadi),1));
 QDbi=sparse(0.5*ones(length(Loadi),1));
@@ -122,7 +129,7 @@ while(abs(Gap)>Precision*1)
     Lz=FormLz(Mat_G,Init_L,GenL,Busnum,PVQL,PD,RealPD,RealQD,Loadi,KK,PF,eps);
     Lw=FormLw(Mat_G,Init_U,GenU,Busnum,PVQU,PD,PD0,QD0,Loadi,KK,PF,eps);
     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);
     %% 开始解方程
     fprintf('迭代次数 %d Gap %f\n',KK+1,plotGap(KK+1));
     XX=SolveIt(deltF,deltG,Init_L,Init_Z,Init_U,Init_W,deltdeltF,ddh,ddgzw,deltH,Init_Y,Ly,Lz,ContrlCount,Lw,Lul,Luu,RestraintCount,Lx,Balance,PVi,PGi,Busnum,Loadi);
@@ -133,8 +140,10 @@ while(abs(Gap)>Precision*1)
     KK=KK+1;
 end
 % 第二遍，算离散
+% [~,~,Init_Z,Init_W,Init_L,Init_U,~,~,~,RestraintCount,~,~,~,~,~,~,~,~,~]=OPF_Init(Busnum,Balance,PG,QG,Volt,GenU,GenL,PVi,PGi,PVQU,PVQL,RealPD,RealQD,QD,PD);
+% Gap=(Init_L*Init_Z'-Init_U*Init_W');
 Gap=1000;
-KK=0;
+% KK=0;
 fprintf('\n');
 fprintf('第2次迭代，算离散量。\n');
 while(abs(Gap)>Precision*1)
@@ -166,12 +175,15 @@ while(abs(Gap)>Precision*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;
-    if KK>8
-        eps=eps*0.3;
+    if KK>22
+        eps=eps*0.1;
         if abs(eps)<1e-6
             eps=1e-5;
         end
         eps;
+       if any(Vbi>0.002)
+           Vbi(Vbi>0.002)=1;
+       end
     end
     Lz=FormLz(Mat_G,Init_L,GenL,Busnum,PVQL,PD,RealPD,RealQD,Loadi,KK,PF,eps);
     Lw=FormLw(Mat_G,Init_U,GenU,Busnum,PVQU,PD,PD0,QD0,Loadi,KK,PF,eps);
@@ -181,10 +193,21 @@ while(abs(Gap)>Precision*1)
     fprintf('迭代次数 %d Gap %f\n',KK+1,plotGap(KK+1));
     XX=SolveIt(deltF,deltG,Init_L,Init_Z,Init_U,Init_W,deltdeltF,ddh,ddgzw,deltH,Init_Y,Ly,Lz,ContrlCount,Lw,Lul,Luu,RestraintCount,Lx,Balance,PVi,PGi,Busnum,Loadi);
     %%取各分量
+%     Vbi_=Vbi;
+%     PDbi_=PDbi;
+%     QDbi_=QDbi;
     [deltZ,deltL,deltW,deltU,deltX,deltY]=AssignXX(XX,ContrlCount,RestraintCount,Busnum);
     [Init_Z,Init_L,Init_W,Init_U,Init_Y,PG,QG,Volt,UAngel,PD,QD,Vbi,PDbi,QDbi]=Modification(Init_Z,Init_L,Init_W,Init_U,Init_Y,deltZ,deltL,deltW,deltU,deltX,deltY,PG,QG,Volt,UAngel,PVi,ContrlCount,Balance,Busnum,PGi,PD,QD,Loadi,Vbi,PDbi,QDbi);
     Gap=(Init_L*Init_Z'-Init_U*Init_W');
     KK=KK+1;
 end
+%% 计算最大不平衡量
+AngleIJ=sparse(r,c,UAngel(r)-UAngel(c)-Angle',Busnum,Busnum);
+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);
 toc
 end
\ No newline at end of file