小修补

Signed-off-by: dugg@lab-desk <dugg@lab-desk>

.gitignore

@@ -0,0 +1,2 @@
+*.asv
+*.7z

@ForThesis/MaxBranchDeviation.asv

@@ -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
-

@ForThesis/MaxDeviation.asv

@@ -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

@ForThesis/StatDeviation.asv

@@ -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
-

AngelBar.m

@@ -0,0 +1,83 @@
+function createfigure(xvector1, yvector1, xvector2, yvector2, xvector3, yvector3)
+%CREATEFIGURE(XVECTOR1,YVECTOR1,XVECTOR2,YVECTOR2,XVECTOR3,YVECTOR3)
+%  XVECTOR1:  bar xvector
+%  YVECTOR1:  bar yvector
+%  XVECTOR2:  bar xvector
+%  YVECTOR2:  bar yvector
+%  XVECTOR3:  bar xvector
+%  YVECTOR3:  bar yvector
+
+%  Auto-generated by MATLAB on 12-Feb-2015 13:48:54
+
+% Create figure
+figure1 = figure('Name','相角直方图','Color',[1 1 1]);
+
+% Create axes
+axes1 = axes('Parent',figure1,'XTick',[-0.0002 -0.0001 0 0.0001 0.0002],...
+    'Position',[0.0920853080568723 0.470619469026549 0.20649289099526 0.372300884955752],...
+    'FontSize',13);
+% Uncomment the following line to preserve the X-limits of the axes
+xlim(axes1,[-0.0003 0.0003]);
+% Uncomment the following line to preserve the Y-limits of the axes
+ylim(axes1,[0 11]);
+box(axes1,'on');
+hold(axes1,'all');
+
+% Create bar
+bar(xvector1,yvector1,'FaceColor',[0 0 0],'Parent',axes1);
+
+% Create xlabel
+xlabel('偏差/rad','FontSize',13);
+
+% Create ylabel
+ylabel('数量','FontSize',13);
+
+% Create title
+title('情形A','FontSize',13);
+
+% Create axes
+axes2 = axes('Parent',figure1,'XTick',[-0.002 -0.001 0 0.001 0.002],...
+    'Position',[0.350370561164917 0.466814159292035 0.20649289099526 0.372300884955752],...
+    'FontSize',13);
+% Uncomment the following line to preserve the X-limits of the axes
+xlim(axes2,[-0.002 0.002]);
+% Uncomment the following line to preserve the Y-limits of the axes
+ylim(axes2,[0 7]);
+box(axes2,'on');
+hold(axes2,'all');
+
+% Create bar
+bar(xvector2,yvector2,'FaceColor',[0 0 0],'Parent',axes2);
+
+% Create xlabel
+xlabel('偏差/rad','FontSize',13);
+
+% Create title
+title('情形B','FontSize',13);
+
+% Create ylabel
+ylabel('数量','FontSize',13);
+
+% Create axes
+axes3 = axes('Parent',figure1,'XTick',[-0.002 -0.001 0 0.001 0.002],...
+    'Position',[0.606286146026517 0.463982300884956 0.20649289099526 0.372300884955752],...
+    'FontSize',13);
+% Uncomment the following line to preserve the X-limits of the axes
+xlim(axes3,[-0.002 0.002]);
+% Uncomment the following line to preserve the Y-limits of the axes
+ylim(axes3,[0 10]);
+box(axes3,'on');
+hold(axes3,'all');
+
+% Create bar
+bar(xvector3,yvector3,'FaceColor',[0 0 0],'Parent',axes3);
+
+% Create xlabel
+xlabel('偏差/rad','FontSize',13);
+
+% Create ylabel
+ylabel('数量','FontSize',13);
+
+% Create title
+title('情形C','FontSize',13);
+

CalTimeBar.m

@@ -0,0 +1,73 @@
+function createfigure(ymatrix1)
+%CREATEFIGURE(YMATRIX1)
+%  YMATRIX1:  bar matrix data
+
+%  Auto-generated by MATLAB on 11-Feb-2015 20:23:32
+
+% Create figure
+figure1 = figure('Name','最大不平衡量','Color',[1 1 1]);
+
+% Create axes
+axes1 = axes('Parent',figure1,'XTickLabel',{'情形A','情形B','情形C'},...
+    'XTick',[1 2 3],...
+    'Position',[0.13 0.11 0.29654028436019 0.391061571125265],...
+    'FontName','Times New Roman');
+hold(axes1,'all');
+
+% Create multiple lines using matrix input to bar
+bar1 = bar(ymatrix1,'BarLayout','stacked');
+set(bar1(1),...
+    'FaceColor',[0.831372559070587 0.815686285495758 0.7843137383461],...
+    'DisplayName','阶段1');
+set(bar1(2),...
+    'FaceColor',[0.941176474094391 0.941176474094391 0.941176474094391],...
+    'DisplayName','fdsfsd');
+
+% Create xlabel
+xlabel('情形','FontSize',12,'FontName','黑体');
+
+% Create ylabel
+ylabel('时间/ms','FontSize',12,'FontName','黑体');
+
+% Create textbox
+annotation(figure1,'textbox',...
+    [0.166919431279624 0.371781316348195 0.0888625592417062 0.059447983014862],...
+    'String',{'29.68ms'},...
+    'FontName','Times New Roman',...
+    'LineStyle','none');
+
+% Create textbox
+annotation(figure1,'textbox',...
+    [0.16706161137441 0.262269639065817 0.0888625592417062 0.059447983014862],...
+    'String',{'70.16ms'},...
+    'FontName','Times New Roman',...
+    'LineStyle','none');
+
+% Create textbox
+annotation(figure1,'textbox',...
+    [0.243791469194317 0.353904458598726 0.0888625592417062 0.059447983014862],...
+    'String',{'31.66ms'},...
+    'FontName','Times New Roman',...
+    'LineStyle','none');
+
+% Create textbox
+annotation(figure1,'textbox',...
+    [0.243933649289103 0.244392781316348 0.0888625592417062 0.059447983014862],...
+    'String',{'68.48ms'},...
+    'FontName','Times New Roman',...
+    'LineStyle','none');
+
+% Create textbox
+annotation(figure1,'textbox',...
+    [0.314881516587683 0.343288747346072 0.0888625592417062 0.059447983014862],...
+    'String',{'30.08ms'},...
+    'FontName','Times New Roman',...
+    'LineStyle','none');
+
+% Create textbox
+annotation(figure1,'textbox',...
+    [0.315023696682469 0.233777070063694 0.0888625592417062 0.059447983014862],...
+    'String',{'65.16ms'},...
+    'FontName','Times New Roman',...
+    'LineStyle','none');
+

CreateMaxErrorFigure.m

@@ -0,0 +1,40 @@
+function createfigure(ymatrix1)
+%CREATEFIGURE(YMATRIX1)
+%  YMATRIX1:  bar matrix data
+
+%  Auto-generated by MATLAB on 11-Feb-2015 21:08:10
+
+% Create figure
+figure1 = figure('Color',[1 1 1]);
+
+% Create axes
+axes1 = axes('Parent',figure1,...
+    'XTick',[0 2 4 6 8 10 12 14 16 18 20 22 24 26 28 30 32 34 35],...
+    'Position',[0.127153700189753 0.112232142857143 0.430721062618596 0.385535714285714],...
+    'FontName','Times New Roman');
+% Uncomment the following line to preserve the X-limits of the axes
+xlim(axes1,[1 35]);
+% Uncomment the following line to preserve the Y-limits of the axes
+ylim(axes1,[0 0.0021]);
+hold(axes1,'all');
+
+% Create multiple lines using matrix input to bar
+bar1 = bar(ymatrix1,'BarWidth',1,'Parent',axes1);
+set(bar1(1),'FaceColor',[0 0 0],'DisplayName','电压幅值');
+set(bar1(2),...
+    'FaceColor',[0.941176474094391 0.941176474094391 0.941176474094391],...
+    'DisplayName','电压相角');
+
+% Create xlabel
+xlabel('节点号','FontSize',12,'FontName','黑体');
+
+% Create ylabel
+ylabel('最大绝对偏差/(p.u，弧度)','FontSize',12,'FontName','黑体');
+
+% Create legend
+legend1 = legend(axes1,'show');
+set(legend1,...
+    'Position',[0.367488931056294 0.487351190476188 0.115749525616698 0.09375],...
+    'FontSize',12,...
+    'FontName','黑体');
+

DeviationFigure.m

@@ -0,0 +1,115 @@
+function createfigure(X1, YMatrix1, YMatrix2, YMatrix3, YMatrix4)
+%CREATEFIGURE(X1,YMATRIX1,YMATRIX2,YMATRIX3,YMATRIX4)
+%  X1:  vector of x data
+%  YMATRIX1:  matrix of y data
+%  YMATRIX2:  matrix of y data
+%  YMATRIX3:  matrix of y data
+%  YMATRIX4:  matrix of y data
+
+%  Auto-generated by MATLAB on 11-Feb-2015 16:31:59
+
+% Create figure
+figure1 = figure('Color',[1 1 1]);
+
+% Create axes
+axes1 = axes('Parent',figure1,'XTick',[1 4 7 10 13 16 19 22 25 28 31 33],...
+    'Position',[0.140663507109006 0.572508000217951 0.537061611374407 0.365998852384935],...
+    'FontName','Times New Roman');
+hold(axes1,'all');
+
+% Create multiple lines using matrix input to plot
+plot1 = plot(X1,YMatrix1,'Parent',axes1,'LineStyle',':','Color',[0 0 0]);
+set(plot1(1),'Marker','diamond','DisplayName','情形A');
+set(plot1(2),'Marker','square','DisplayName','情形B');
+set(plot1(3),'Marker','o','DisplayName','情形C');
+
+% Create xlabel
+xlabel('节点号','FontName','黑体');
+
+% Create ylabel
+ylabel('偏差/p.u','FontName','Times New Roman');
+
+% Create axes
+axes2 = axes('Parent',figure1,'XTick',[2 5 8 11 14 17 20 23 26 29 32],...
+    'Position',[0.143033175355449 0.120301205914361 0.537914691943128 0.365998852384935],...
+    'FontName','Times New Roman');
+% Uncomment the following line to preserve the X-limits of the axes
+% xlim(axes2,[2 33]);
+hold(axes2,'all');
+
+% Create multiple lines using matrix input to plot
+plot2 = plot(X1,YMatrix3,'Parent',axes2,'LineStyle',':','Color',[0 0 0]);
+set(plot2(1),'Marker','diamond','DisplayName','算例A');
+set(plot2(2),'Marker','square','DisplayName','算例B');
+set(plot2(3),'Marker','o','DisplayName','算例C');
+
+% Create xlabel
+xlabel('节点号','FontName','Times New Roman');
+
+% Create ylabel
+ylabel('相对偏差%','FontName','Times New Roman');
+
+% Create title
+title('有功');
+
+% Create axes
+axes3 = axes('Parent',figure1,...
+    'XTickLabel',{'1','4','7','10','13','16','19','22','25','28','31','33'},...
+    'XTick',[1 4 7 10 13 16 19 22 25 28 31 33],...
+    'Position',[0.440758293838865 0.574607379571982 0.537914691943128 0.365998852384935]);
+hold(axes3,'all');
+
+% Create multiple lines using matrix input to plot
+plot3 = plot(X1,YMatrix2,'Parent',axes3,'LineStyle',':','Color',[0 0 0]);
+set(plot3(1),'Marker','diamond','DisplayName','情形A');
+set(plot3(2),'Marker','square','DisplayName','情形B');
+set(plot3(3),'Marker','o','DisplayName','情形C');
+
+% Create xlabel
+xlabel('节点号');
+
+% Create ylabel
+ylabel('偏差/弧度');
+
+% Create axes
+axes4 = axes('Parent',figure1,'XTick',[2 5 8 11 14 17 20 23 26 29 32],...
+    'Position',[-0.490853080568712 -0.0598281717337482 0.537914691943128 0.365998852384935]);
+% Uncomment the following line to preserve the X-limits of the axes
+% xlim(axes4,[2 33]);
+hold(axes4,'all');
+
+% Create multiple lines using matrix input to plot
+plot4 = plot(X1,YMatrix4,'Parent',axes4,'LineStyle',':','Color',[0 0 0]);
+set(plot4(1),'Marker','diamond','DisplayName','算例A');
+set(plot4(2),'Marker','square','DisplayName','算例B');
+set(plot4(3),'Marker','o','DisplayName','算例C');
+
+% Create xlabel
+xlabel('节点号');
+
+% Create ylabel
+ylabel('相对偏差%');
+
+% Create legend
+legend1 = legend(axes1,'show');
+set(legend1,...
+    'Position',[0.155301294649576 0.633902089695378 0.112559241706161 0.12206572769953],...
+    'FontName','黑体');
+
+% Create legend
+legend2 = legend(axes2,'show');
+set(legend2,...
+    'Position',[0.533220206544619 0.178635859349444 0.112559241706161 0.12206572769953],...
+    'FontName','黑体');
+
+% Create legend
+legend3 = legend(axes3,'show');
+set(legend3,...
+    'Position',[0.450281817919025 0.807479462339755 0.117298578199052 0.12206572769953],...
+    'FontName','黑体');
+
+% Create legend
+legend4 = legend(axes4,'show');
+set(legend4,...
+    'Position',[-0.098296381133102 -0.0539564029130473 0.117298578199052 0.122065727699531]);
+

DrawGap.asv

@@ -1,6 +0,0 @@
-function DrawGap(plotGap)
-x=find(plotGap);
-ts=size(x,2);
-
-
-end

DrawLoadProfile.m

@@ -0,0 +1,20 @@
+function [  ] = DrawLoadProfile( )
+figure();
+data=dlmread('LoadProfile.txt');
+hold on;
+y1=data(:,3);
+y2=data(:,4);
+x=1:96;
+k=1;
+for i=0:96-1
+    x1(k)=x(length(x)-i);
+    Z2(k)=y2(length(x)-i);
+    k=k+1;
+end
+h=patch([x x1],[y1' Z2],[230 230 230]./255)
+set(h,'EdgeColor','None');
+plot(1:96,data(:,2),'k');
+plot(1:96,data(:,3),'--k');
+plot(1:96,data(:,4),'--k');
+end
+

FormAA.asv

@@ -1,12 +0,0 @@
-function AA=FormAA(L_1Z,deltG,U_1W,Hcoma,deltH)
-tOnes=eye(14);
-tZeros=zeros(14);
-AA=[
-    tOnes,L_1Z,tZeros,tZeros,tZeros,zeros(14,10);
-    tZeros,tOnes,tZeros,tZeros,-deltG',zeros(14,10);
-    tZeros,tZeros,tOnes,U_1W,tZeros,zeros(14,10);
-    tZeros,tZeros,tZeros,tOnes,deltG,zeros(14,10);
-    tZeros,tZeros,tZeros,tZeros,Hcoma,deltH;
-    zeros(10,14),zeros(10,14),zeros(10,14),zeros(10,14),deltH',zeros(10,10);
-    ];
-end

FormAlphaP.asv

@@ -1,7 +0,0 @@
-function AlphaP=FormAlphaP(Init_L,deltL,Init_U,deltU)
-ti=deltL(delt)
-t1=-Init_L./deltL';
-t2=-Init_U./deltU';
-t3=[t1,t2];
-t4=t3()
-end

FormG.asv

@@ -1,11 +0,0 @@
-function Mat_G=FormG(Volt,PVi,PGi,PG,QG,PD,QD,Loadi)
-
-Mat_G=[
-    PG(PGi);
-    QG(PVi);
-    sparse(PD(Loadi));
-    sparse(QD(Loadi));
-    Volt';
-    
-    ];
-end

FormG.m

@@ -1,10 +1,17 @@
-function Mat_G=FormG(Volt,PVi,PGi,PG,QG,PD,QD,Loadi,Vbi)
+function Mat_G=FormG(Volt,PD,QD,Loadi,bigM,mVolt,rVolt,sigma,Vbi,PDbi,QDbi,mPD,mQD,uPD,lPD,uQD,lQD,uVolt,lVolt)
 
 Mat_G=[
-    sparse(PD(Loadi));
+    sparse(PD(Loadi))-mPD(Loadi)-bigM*PDbi-uPD(Loadi);
-    sparse(QD(Loadi));
+    sparse(PD(Loadi))-mPD(Loadi)+bigM*PDbi+lPD(Loadi);
-    Volt';
+    sparse(QD(Loadi))-mQD(Loadi)-bigM*QDbi-uQD(Loadi);
-    Volt';
+    sparse(QD(Loadi))-mQD(Loadi)+bigM*QDbi+lQD(Loadi);
+    Volt'-mVolt'-bigM*Vbi-uVolt;
+    Volt'-mVolt'+bigM*Vbi+lVolt;
     Vbi;
+    PDbi;
+    QDbi;
+    (0.1-Vbi).*Vbi;
+    (0.1-PDbi).*PDbi;
+    (0.1-QDbi).*QDbi;
     ];
 end

FormH.asv

@@ -1,31 +0,0 @@
-function Mat_H=FormH(Busnum,GB,AngleIJMat,Volt,PG,PD,QG,QD,Y)
-t1=real(GB).*cos(AngleIJMat)+imag(GB).*sin(AngleIJMat);
-t2=Volt'*Volt;
-t3=t1.*t2;
-t4=sum(-t3,2);%P
-t5=real(GB).*sin(AngleIJMat)-imag(GB).*cos(AngleIJMat);
-t6=t2.*t5;
-t7=sum(-t6,2);%Q
-t8=PG-PD;
-t9=QG-QD;
-%Mat_H=([(PG-PD)',(QG-QD)'])'+([t4',t7'])';
-Mat_H(1:2:2*Busnum)=t8(1:Busnum)+t4(1:Busnum);
-Mat_H(2:2:2*Busnum)=t9(1:Busnum)+t7(1:Busnum);
-Mat_H=Mat_H';
-%% 
-QDcos=textread('300glys.txt');
-t=QD(PD==0 &&);
-aa=QD;
-QD(QD~=0)=PD(QD~=0)./tan(QDcos);
-QD(PD==0)=t;
-%% 
-%%%%一下是学姐给的公式
-AngleIJ=AngleIJMat-angle(GB);
-%dP=PG-PD-diag(Volt)*Y*cos(AngleIJ)*Volt';
-dP=PG-PD-diag(Volt)*Y.*cos(AngleIJ)*Volt';
-dQ=QG-QD-diag(Volt)*Y.*sin(AngleIJ)*Volt';
-%Mat_H(1:2:2*Busnum)=dP(1:Busnum);暂时改一下 20111227
-%Mat_H(2:2:2*Busnum)=dQ(1:Busnum);暂时改一下 20111227
-Mat_H=[dP;dQ;];
-
-end

FormH.m

@@ -2,7 +2,5 @@ function Mat_H=FormH(Busnum,Volt,PG,PD,QG,QD,Y,UAngel,r,c,Angle,Loadi)
 AngleIJ=sparse(r,c,UAngel(r)-UAngel(c)-Angle',Busnum,Busnum);
 dP=PG-PD-diag(Volt)*Y.*cos(AngleIJ)*Volt';
 dQ=QG-QD-diag(Volt)*Y.*sin(AngleIJ)*Volt';
-
 Mat_H=[dP;dQ;];
-
 end

FormLw.m

@@ -1,11 +1,7 @@
-function Lw=FormLw(Mat_G,Init_U,GenU,Busnum,PVQU,PD,PD0,QD0,Loadi,KK,PF,noDataTransCapacity)
+function Lw=FormLw(Mat_G,Init_U,GenU,Busnum,PVQU,PD,PD0,QD0,Loadi,KK,PF,eps)
 KK=999;
-%PU=GenU(:,2);%发电机有功上界
-
 VoltU=(1.1)*ones(1,Busnum);
-%VoltU=10*ones(1,Busnum);
 PDU=PD0(Loadi);
-% PDU=noDataTransCapacity;
 PDU(PDU>0)=1.200*PDU(PDU>0);
 PDU(PDU<0)=0.800*PDU(PDU<0);
 PDU(PDU==0)=0.400;
@@ -13,7 +9,6 @@ realPD=PD0(Loadi);
 indPD=find(realPD>0);
 PDU(indPD(3:12:end))=1.55*realPD(indPD(3:12:end));
 PDU(indPD(9:12:end))=1.05*realPD(indPD(9:12:end));
-%PDU=10*ones(length(Loadi),1);
 QDU=QD0(Loadi);
 QDU(QDU>0)=1.200*QDU(QDU>0);
 QDU(QDU<0)=0.800*QDU(QDU<0);
@@ -22,9 +17,12 @@ realQD=QD0(Loadi);
 indQD=find(realQD>0);
 QDU(indQD(3:12:end))=1.55*realQD(indQD(3:12:end));
 QDU(indQD(9:12:end))=1.05*realQD(indQD(9:12:end));
-% PF=0.85;
+t1=([0*PDU',1*ones(1,length(PDU)),0*QDU',1*ones(1,length(QDU)),0*VoltU,1*ones(1,length(VoltU)),0.1*ones(1,Busnum+length(Loadi)*2)])';
-% QDU=1.0*PD(Loadi).*sqrt(1 -PF.^2)./PF;
+t1=[t1;
-t1=([PDU',QDU',VoltU,VoltU,0*ones(1,Busnum)])';
+    sparse(eps*ones(Busnum,1));
+   sparse(eps*ones(length(Loadi),1));
+    sparse(eps*ones(length(Loadi),1));
+    ];
 t2=Mat_G+Init_U'-t1;
 Lw=t2;
 

FormLxComa.m

@@ -4,8 +4,9 @@ t2=Lul+diag(Init_Z)*Lz;
 t3=inv(diag(Init_L));
 t4=t3*t2;%
 t5=Luu-diag(Init_W)*Lw;
-t6=inv(diag(Init_U));
+% t6=inv(diag(Init_U));
-t7=t6*t5;%
+% t7=t6*t5;%
+t7=diag(Init_U)\t5;
 t8=deltG*(t4+t7);%%
 LxComa=Lx+t8;
 end

FormLz.m

@@ -1,4 +1,4 @@
-function Lz=FormLz(Mat_G,Init_L,GenL,Busnum,PVQL,PD,PD0,QD0,Loadi,KK,PF)
+function Lz=FormLz(Mat_G,Init_L,GenL,Busnum,PVQL,PD,PD0,QD0,Loadi,KK,PF,eps)
 KK=999;
 
 VoltL=(0.9)*ones(1,Busnum);
@@ -21,7 +21,12 @@ indQD=find(realQD>0);
 QDL(indQD(3:12:end))=0.95*realQD(indQD(3:12:end));
 QDL(indQD(9:12:end))=0.95*realQD(indQD(9:12:end));
 % QDL=0*PD(Loadi).*sqrt((1-PF.^2))./PF;
-t1=([PDL',QDL',VoltL,VoltL,0*ones(1,Busnum)])';
+t1=([-1*ones(1,length(PDL)),0*PDL',-1*ones(1,length(QDL)),0*QDL',-1*ones(1,length(VoltL)),0*VoltL,0*ones(1,Busnum+length(Loadi)*2)])';
+t1=[t1;
+    sparse(-eps*ones(Busnum,1));
+   sparse(-eps*ones(length(Loadi),1));
+   sparse(-eps*ones(length(Loadi),1));
+];
 t2=Mat_G-Init_L'-t1;
 Lz=t2;
 

FormYY.asv

@@ -1,5 +0,0 @@
-function FormYY(Init_L,Lul.Lz,Init_U,Luu,Lw,Lz,LxComa)
-t=[
-   -    
-]
-end

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)

MaxBoundErrorFigure.m

@@ -0,0 +1,9 @@
+function MaxBoundErrorFigure()
+figure();
+data=dlmread('maxBoundError.txt');
+% data(1,:)=0;
+maxV=data(:,1);
+maxA=data(:,2);
+% bar([maxV';maxA']');
+CreateMaxErrorFigure([maxV';maxA']');
+end

MaxErrorFigure.m

@@ -0,0 +1,9 @@
+function MaxErrorFigure()
+figure();
+data=dlmread('maxError.txt');
+data(1,:)=0;
+maxV=data(:,1);
+maxA=data(:,2);
+% bar([maxV';maxA']');
+CreateMaxErrorFigure([maxV';maxA']');
+end

Modification.asv

@@ -1,18 +0,0 @@
-function [Init_Z,Init_L,Init_W,Init_U,Init_Y,deltZ,deltL,deltW,deltU,deltX,deltY,PG,QG,Volt,UAngel]=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)
-AlphaP=FormAlphaP(Init_L,deltL,Init_U,deltU);
-fprintf('AlphaP %f\n',AlphaP);
-AlphaD=FormAlphaD(Init_Z,deltZ,Init_W,deltW);
-fprintf('AlphaD %f\n',AlphaD);
-
-Init_Z=Init_Z+AlphaD*deltZ';
-Init_L=Init_L+AlphaP*deltL';
-Init_W=Init_W+AlphaD*deltW';
-Init_U=Init_U+AlphaP*deltU';
-Init_Y=Init_Y+AlphaD*deltY';
-PG(PVi)=PG(PGi)+deltX(size(PGi,1));
-QG(PVi)=QG(PVi)+deltX(size(PGi,1)+1,1:(size(PVi,1)+size(PGi,1)));
-t=deltX(size(PVi,1)+size(PGi,1)+1,ContrlCount)';
-t(2*Balance-1)=0;
-Volt=Volt+AlphaP*t(2:2:2*Busnum);
-UAngel=UAngel+AlphaP*t(1:2:2*Busnum);
-end

Modification.m

@@ -1,8 +1,16 @@
-function [Init_Z,Init_L,Init_W,Init_U,Init_Y,PG,QG,Volt,UAngel,PD,QD,Vbi]=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)
+function [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)
 AlphaP=FormAlphaP(Init_L,deltL,Init_U,deltU);
 %fprintf('AlphaP %f\n',full(AlphaP));
 AlphaD=FormAlphaD(Init_Z,deltZ,Init_W,deltW);
 %fprintf('AlphaD %f\n',full(AlphaD));
+% if AlphaP<1e-4
+%     AlphaP=0.1;
+% end
+% 
+% if AlphaD<1e-4
+%     AlphaD=0.1;
+% end
+
 
 Init_Z=Init_Z+AlphaD*deltZ';
 Init_L=Init_L+AlphaP*deltL';
@@ -24,5 +32,7 @@ balVolt=Volt(Balance);
 Volt=Volt+AlphaP*t(1:Busnum);
 Volt(Balance)=balVolt;
 UAngel=UAngel+AlphaP*t(Busnum+1:2*Busnum);
-Vbi=Vbi+AlphaP*t(2*Busnum+1:end)';
+Vbi=Vbi+AlphaP*t(2*Busnum+1:2*Busnum+Busnum)';
+PDbi=PDbi+AlphaP*t(3*Busnum+1:3*Busnum+length(Loadi))';
+QDbi=QDbi+AlphaP*t(3*Busnum+length(Loadi)+1:end)';
 end

OPF.m

@@ -1,97 +1,406 @@
-tic
 clc
 clear
-%% 存在问题
+close all
-% 变压器变比的位置没有考虑，由于现在用的变比都是1，所以没有影响。 20130123
+arrayA=zeros(21,10);
-%%
+sumCaseA_SE=0;
-thesis=ForThesis(1,62);
+sumCaseB_SE=0;
-[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]= ...
+sumCaseC_SE=0;
-    pf('C:\lz\结题后重新计算\重构结果\78\木材厂堤Ⅰ线_2229880_2013-10-29\木材厂堤Ⅰ线_2229880_2013-10-29_iPso_newFile.txt');
+badDataResult=zeros(10,33);
-%         pf('E:\算例\柳金Ⅰ926_21671693_2012-09-06\newFIle16.txt');
+badDataLocation=zeros(34,10);
-%pf('D:\Project\青秀降损项目\最小化潮流\最小潮流算例\原始\津头站津视922（3-1）_0.5_120%.txt');
+nodeMaxDVolt=zeros(33,1);
-%pf('D:\Project\最小化潮流\最小潮流算例\仙海919.txt');
+nodeMaxDVAngle=zeros(33,1);
-%pf('c:/file31.txt');
+for badDataNode=1:1
-%% 计算功率因数
+    loopN=1;
-Loadi=QD~=0 | PD~=0;
+    maxDVolt=0;
-PF=sqrt(PD(Loadi).^2./(QD(Loadi).^2+PD(Loadi).^2));
+    maxDVAngle=0;
-%%
+    while 1
-Volt;
+        close
-UAngel*180/3.1415926;
+        [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]= ...
-%% 通过潮流计算PG
+            pf('E:\算例\feeder33\feeder33.txt');
-AngleIJ=sparse(r,c,UAngel(r)-UAngel(c)-Angle',Busnum,Busnum);
+        sigma=0.01;
-PGBal=PD+diag(Volt)*Y.*cos(AngleIJ)*Volt';
+        RealPD=PD;
-QGBal=QD+diag(Volt)*Y.*sin(AngleIJ)*Volt';
+        RealQD=QD;
-%% 初值-即测量值
+        rVolt=Volt;
-PG0=PG;
+        Loadi=find(PD~=0);
-QG0=QG;
+        PD0=sparse(Busnum,1);
-PD0=PD;
+        QD0=sparse(Busnum,1);
-QD0=QD;
+        PD0(Loadi)=RealPD(Loadi).*(1+normrnd(0,sigma,length(Loadi),1));
-Volt0=Volt;
+        QD0(Loadi)=RealQD(Loadi).*(1+normrnd(0,sigma,length(Loadi),1));
-UAngel0=UAngel;
+        mVolt=rVolt.*(1+normrnd(0,sigma,length(rVolt),1))';
-%%
+        
-PG0(Balance)=PGBal(Balance);
+        %加载保存的变量
-PG(Balance)=PGBal(Balance);
+        PD0=load('PD0');
-QG0(Balance)=QGBal(Balance);
+        PD0=PD0.PD0;
-QG0(PVi)=QGBal(PVi);
+        QD0=load('QD0');
-QG(PVi)=QGBal(PVi);
+        QD0=QD0.QD0;
-%% 真实值
+        mVolt=load('mVolt');
-RealPG=PG0;
+        mVolt=mVolt.mVolt;
-RealQG=QG0;
+        
-RealPD=PD0;
+%         mVolt(3)=rVolt(3)*(1-sigma*6);
-RealQD=QD0;
+        %% 画Case A的图
-%%
+        % figure('Color',[1 1 1]);
-[Volt,UAngel,Init_Z,Init_W,Init_L,Init_U,Init_Y,PG,QG,RestraintCount,wPG,wQG,wPD,wQD,PD,PD0,QD,randPDind,Loadi]=OPF_Init(Busnum,Balance,PG,QG,Volt,GenU,GenL,PVi,PGi,PVQU,PVQL,RealPD,RealQD,QD,PD);
+        [Busnum,Loadi,Volt,PD,QD,rVolt,UAngel,RealPD,RealQD,rUAngel,Vbi,PDbi,QDbi,plotGapA,isConverge]=subOPF([],PD0,QD0,mVolt,sigma);%全部有
-Gap=(Init_L*Init_Z'-Init_U*Init_W');
+        %         badDataResult(I,badDataNode)=sum(Vbi);
-KK=0;
+        %         badDataLocation(1:33,I)=Vbi;
-plotGap=zeros(1,60);
+        %         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))));
-ContrlCount=size(Loadi,1)*2+Busnum*2+Busnum;
+                break;
-kmax=60;
+        if isConverge==0
-Precision=Precision/1;
+            continue;
-%% 加误差
+        end
-PD0(Loadi)=PD0(Loadi).*(1+normrnd(0,0.05,length(Loadi),1));
+        % 加最大偏差
-QD0(Loadi)=QD0(Loadi).*(1+normrnd(0,0.05,length(Loadi),1));
+        maxDVolt_=max(abs((rVolt-Volt)));
-Vbi=sparse(ones(Busnum,1));
+        if maxDVolt_>maxDVolt
-while(abs(Gap)>Precision)
+            maxDVolt=maxDVolt_;
-    if KK>kmax
+        end
+        maxDVAngle_=max(abs((UAngel(2:33)-rUAngel(2:33))));
+        if maxDVAngle_>maxDVAngle
+            maxDVAngle=maxDVAngle_;
+        end
+        nodeMaxDVolt_t=abs((rVolt-Volt))';
+        nodeMaxDVAngle_t=abs((UAngel-rUAngel))';
+        nodeMaxDVolt_t([1:18,23:33])=0;
+        nodeMaxDVAngle_t([1:18,23:33])=0;
+        nodeMaxDVolt(nodeMaxDVolt<nodeMaxDVolt_t)=nodeMaxDVolt_t(nodeMaxDVolt<nodeMaxDVolt_t);
+        nodeMaxDVAngle(nodeMaxDVAngle<nodeMaxDVAngle_t)=nodeMaxDVAngle_t(nodeMaxDVAngle<nodeMaxDVAngle_t);
+        if loopN>=500
+%             nodeMaxDVolt(badDataNode)=maxDVolt;
+%             nodeMaxDVAngle(badDataNode)=maxDVAngle;
             break;
         end
-    plotGap(KK+1)=Gap;
-    Init_u=Gap/2/RestraintCount*CenterA;
-    AngleIJMat=0;
-    %% 开始计算OPF
-    %% 形成等式约束的雅克比
-    deltH=func_deltH(Busnum,Volt,PVi,Y,PGi,UAngel,r,c,Angle,Loadi);
-    %% 形成不等式约束的雅克比
-    deltG=func_deltG(Busnum,PVi,PGi,Loadi,PD,QD);
-    %%
-    L_1Z=diag(Init_Z./Init_L);
-    U_1W=diag(Init_W./Init_U);
-    %% 形成海森阵
-    deltdeltF=func_deltdeltF(PVi,wPG,wQG,wPD,wQD,ContrlCount);
-    %% 形成ddHy
-    ddh=func_ddh(Volt,Init_Y,Busnum,PVi,PGi,Y,UAngel,r,c,Angle,Loadi,ContrlCount);
-    %%  开始构建ddg
-    ddg=func_ddg(PGi,PVi,Busnum,RestraintCount,Loadi,PD,QD);
-    %%  开始构建deltF
-    deltF=func_deltF(PG,QG,PVi,PGi,wPG,wQG,wPD,wQD,PG0,QG0,PD0,PD,QD,QD0,Busnum,Loadi);
-    %% 形成方程矩阵
-    Luu=Init_U'.*Init_W'+Init_u*ones(RestraintCount,1);
-    Lul=Init_L'.*Init_Z'-Init_u*ones(RestraintCount,1);
-    Mat_G=FormG(Volt,PVi,PGi,PG,QG,PD,QD,Loadi,Vbi);
-    Mat_H=FormH(Busnum,Volt,PG,PD,QG,QD,Y,UAngel,r,c,Angle,Loadi);
-    Ly=Mat_H;
-    Lz=FormLz(Mat_G,Init_L,GenL,Busnum,PVQL,PD,RealPD,RealQD,Loadi,KK,PF);
-    Lw=FormLw(Mat_G,Init_U,GenU,Busnum,PVQU,PD,RealPD,RealQD,Loadi,KK,PF);
-    Lx=FormLx(deltF,deltH,Init_Y,deltG,Init_Z,Init_W);
-    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,ddg,deltH,Init_Y,Ly,Lz,ContrlCount,Lw,Lul,Luu,RestraintCount,Lx,Balance,PVi,PGi,Busnum,Loadi);
-    %%取各分量
-    [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]=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);
-    Gap=(Init_L*Init_Z'-Init_U*Init_W');
-    KK=KK+1;
-end
-toc
 
+        loopN=loopN+1;
+    end
+    
+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);%Relative Error of Voltage in Case A
+% CaseAREV=CaseAREV(2:end)*100;
+%真实值的
+plot(1:length(CaseAREV),(CaseAREV),'k.:','Marker','diamond');
+%测量值的
+% plot(1:length(CaseAREV),abs((mVolt-rVolt)*100),'c.:','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);
+% 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),(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=(PD-RealPD)./RealPD*100;%Relative Error of PD in Case A
+CaseAREP(1)=0;
+%真实值的
+plot(1:length(CaseAREP),(CaseAREP),'k:','Marker','diamond');
+%测量值的
+% plot(1:length(CaseAREV),abs((PD0-RealPD)./(RealPD+0.00001)*100),'c.:','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);
+% 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),(CaseAREQ),'k:','Marker','diamond');
+%测量值的
+% plot(1:length(CaseAREV),abs((QD0-RealQD)./(RealQD+0.00001)*100),'c.:','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));
+notZeros=find(PD0~=0);
+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;
+% arrayA(21,I)=CaseAE*1000;
+%% Case B
+% [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);%
+[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=(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=(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=(PD-RealPD)./RealPD*100;%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=(QD-RealQD)./RealQD*100;%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(((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=(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=(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=(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=(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');
+ld=legend('算例A','算例B','算例C');
+set(ld,'Position',[0.847865087908145 0.786094477711244 0.0543595263724435 0.0605455755156354]);
+subplot(4,1,2);
+% 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');
+ld=legend('算例A','算例B','算例C');
+set(ld,'Position',[0.847865087908145 0.386094477711244 0.0543595263724435 0.0605455755156354]);
+subplot(4,1,4);
+% 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));
+
+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;
+
+% 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) %画出概率密度分布图
+voltBarCaseAX=x;
+voltBarCaseAY=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) %画出概率密度分布图
+voltBarCaseBX=x;
+voltBarCaseBY=yy;
+xlabel('Error');
+ylabel('Number of buses');
+title('算例B');
+% 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) %画出概率密度分布图
+voltBarCaseCX=x;
+voltBarCaseCY=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) %画出概率密度分布图
+angelBarCaseAX=x;
+angelBarCaseAY=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) %画出概率密度分布图
+angelBarCaseBX=x;
+angelBarCaseBY=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); %计算各个区间的个数
+angelBarCaseCX=x;
+angelBarCaseCY=yy;
+% 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) %画出概率密度分布图
+PDBarCaseAX=x;
+PDBarCaseAY=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) %画出概率密度分布图
+PDBarCaseBX=x;
+PDBarCaseBY=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) %画出概率密度分布图
+PDBarCaseCX=x;
+PDBarCaseCY=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) %画出概率密度分布图
+QDBarCaseAX=x;
+QDBarCaseAY=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) %画出概率密度分布图
+QDBarCaseBX=x;
+QDBarCaseBY=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) %画出概率密度分布图
+QDBarCaseCX=x;
+QDBarCaseCY=yy;
+%画收敛曲线
+fz=find(abs(plotGapA)==0);
+% fz=fz(1);
+figure('Name','互补曲线')
+plot(1:fz-1,plotGapA(1:fz-1));
+figure('Name','最大不平衡量');
+%% 最大不平衡量
+% maxDismatchPQ = [0.3123e-10  0.1497e-10 0.7351e-10; 0.6854e-10 0.1973e-10 0.5824e-10];
+% bar(maxDismatchPQ);
+%% 计算时间
+calTime=[70.16 29.68; 68.48 31.661; 65.156 30.08;];
+bar(calTime,'stacked');
+% figure();
+% DeviationFigure(1:33,[CaseAREV;CaseBREV;CaseCREV],[CaseAREA;CaseBREA;CaseCREA],[CaseAREP,CaseBREP,CaseCREP],[CaseAREQ,CaseBREQ,CaseCREQ]);%,[CaseAREA;CaseBREA;CaseCREA],[CaseAREV;CaseBREV;CaseCREV],[CaseAREV;CaseBREV;CaseCREV]);
+% VoltBar(voltBarCaseAX,voltBarCaseAY,voltBarCaseBX,voltBarCaseBY,voltBarCaseCX,voltBarCaseCY);
+% AngelBar(angelBarCaseAX,angelBarCaseBY,angelBarCaseCX,angelBarCaseAY,angelBarCaseBX,angelBarCaseCY);
+% PDBar(PDBarCaseAX,PDBarCaseAY,PDBarCaseBX,PDBarCaseBY,PDBarCaseCX,PDBarCaseCY);
+% QDBar(QDBarCaseAX,QDBarCaseAY,QDBarCaseBX,QDBarCaseBY,QDBarCaseCX,QDBarCaseCY);
+% MaxErrorFigure()
+% MaxBoundErrorFigure();
+% DrawLoadProfile();

OPF_Init.m

@@ -1,7 +1,9 @@
 function [Volt,UAngel,Init_Z,Init_W,Init_L,Init_U,Init_Y,PG,QG,RestraintCount,wPG,wQG,wPD,wQD,PD,PD0,QD,randPDind,Loadi]=OPF_Init(Busnum,Balance,PG,QG,Volt,GenU,GenL,PVi,PGi,PVQU,PVQL,PD0,QD0,QD,PD)
 Loadi=find(QD~=0 | PD~=0);
 %Loadi=[1:Busnum]';
-RestraintCount=size(Loadi,1)*2+Busnum*2+Busnum; %约束条件数,放开所有QD
+RestraintCount=size(Loadi,1)*4+Busnum*2+Busnum+length(Loadi)*2; %约束条件数,放开所有QD
+%互补约束数量
+RestraintCount=RestraintCount+Busnum+length(Loadi)*2;
 t_Bal_volt=Volt(Balance);
 Volt=sparse(1*ones(1,Busnum));
 Volt(Balance)=t_Bal_volt;
@@ -27,8 +29,8 @@ wQD(1:2:end)=0;
 %wD(randPDind)=0;%一些负荷不约束
 %wD(7)=0;
 % wD(11)=0;
-PD=0.8*PD0;
+PD=0.0*PD0;
 %powerFacter=0.98;
 %QD=PD*sqrt((1-powerFacter^2)/powerFacter^2);
-QD=.8*QD0;
+QD=.0*QD0;
 end

PDBar.m

@@ -0,0 +1,83 @@
+function createfigure(xvector1, yvector1, xvector2, yvector2, xvector3, yvector3)
+%CREATEFIGURE(XVECTOR1,YVECTOR1,XVECTOR2,YVECTOR2,XVECTOR3,YVECTOR3)
+%  XVECTOR1:  bar xvector
+%  YVECTOR1:  bar yvector
+%  XVECTOR2:  bar xvector
+%  YVECTOR2:  bar yvector
+%  XVECTOR3:  bar xvector
+%  YVECTOR3:  bar yvector
+
+%  Auto-generated by MATLAB on 11-Feb-2015 17:32:07
+
+% Create figure
+figure1 = figure('Name','相角直方图','Color',[1 1 1]);
+
+% Create axes
+axes1 = axes('Parent',figure1,...
+    'Position',[0.0920853080568723 0.470619469026549 0.20649289099526 0.372300884955752],...
+    'FontSize',13);
+% Uncomment the following line to preserve the X-limits of the axes
+% xlim(axes1,[-4 4]);
+% Uncomment the following line to preserve the Y-limits of the axes
+% ylim(axes1,[0 6]);
+box(axes1,'on');
+hold(axes1,'all');
+
+% Create bar
+bar(xvector1,yvector1,'FaceColor',[0 0 0],'Parent',axes1);
+
+% Create xlabel
+xlabel('相对误差%','FontSize',13);
+
+% Create ylabel
+ylabel('数量','FontSize',13);
+
+% Create title
+title('情形A','FontSize',13);
+
+% Create axes
+axes2 = axes('Parent',figure1,'XTick',[-12 -6 0 6 12],...
+    'Position',[0.350370561164917 0.466814159292035 0.20649289099526 0.372300884955752],...
+    'FontSize',13);
+% Uncomment the following line to preserve the X-limits of the axes
+% xlim(axes2,[-12 12]);
+% Uncomment the following line to preserve the Y-limits of the axes
+% ylim(axes2,[0 8]);
+box(axes2,'on');
+hold(axes2,'all');
+
+% Create bar
+bar(xvector2,yvector2,'FaceColor',[0 0 0],'Parent',axes2);
+
+% Create xlabel
+xlabel('相对误差%','FontSize',13);
+
+% Create title
+title('情形B','FontSize',13);
+
+% Create ylabel
+ylabel('数量','FontSize',13);
+
+% Create axes
+axes3 = axes('Parent',figure1,'XTick',[-12 -6 0 6 12],...
+    'Position',[0.606286146026517 0.463982300884956 0.20649289099526 0.372300884955752],...
+    'FontSize',13);
+% Uncomment the following line to preserve the X-limits of the axes
+% xlim(axes3,[-12 12]);
+% Uncomment the following line to preserve the Y-limits of the axes
+% ylim(axes3,[0 5]);
+box(axes3,'on');
+hold(axes3,'all');
+
+% Create bar
+bar(xvector3,yvector3,'FaceColor',[0 0 0],'Parent',axes3);
+
+% Create xlabel
+xlabel('相对误差%','FontSize',13);
+
+% Create ylabel
+ylabel('数量','FontSize',13);
+
+% Create title
+title('情形C','FontSize',13);
+

QDBar.m

@@ -0,0 +1,77 @@
+function createfigure(xvector1, yvector1, xvector2, yvector2, xvector3, yvector3)
+%CREATEFIGURE(XVECTOR1,YVECTOR1,XVECTOR2,YVECTOR2,XVECTOR3,YVECTOR3)
+%  XVECTOR1:  bar xvector
+%  YVECTOR1:  bar yvector
+%  XVECTOR2:  bar xvector
+%  YVECTOR2:  bar yvector
+%  XVECTOR3:  bar xvector
+%  YVECTOR3:  bar yvector
+
+%  Auto-generated by MATLAB on 11-Feb-2015 17:36:00
+
+% Create figure
+figure1 = figure('Name','相角直方图','Color',[1 1 1]);
+
+% Create axes
+axes1 = axes('Parent',figure1,...
+    'Position',[0.0920853080568723 0.470619469026549 0.20649289099526 0.372300884955752],...
+    'FontSize',13);
+% Uncomment the following line to preserve the X-limits of the axes
+xlim(axes1,[-3 3]);
+box(axes1,'on');
+hold(axes1,'all');
+
+% Create bar
+bar(xvector1,yvector1,'FaceColor',[0 0 0],'Parent',axes1);
+
+% Create xlabel
+xlabel('相对误差%','FontSize',13);
+
+% Create ylabel
+ylabel('数量','FontSize',13);
+
+% Create title
+title('情形A','FontSize',13);
+
+% Create axes
+axes2 = axes('Parent',figure1,'XTick',[-12 -6 0 6 12],...
+    'Position',[0.350370561164917 0.466814159292035 0.20649289099526 0.372300884955752],...
+    'FontSize',13);
+% Uncomment the following line to preserve the X-limits of the axes
+xlim(axes2,[-8 8]);
+box(axes2,'on');
+hold(axes2,'all');
+
+% Create bar
+bar(xvector2,yvector2,'FaceColor',[0 0 0],'Parent',axes2);
+
+% Create xlabel
+xlabel('相对误差%','FontSize',13);
+
+% Create title
+title('情形B','FontSize',13);
+
+% Create ylabel
+ylabel('数量','FontSize',13);
+
+% Create axes
+axes3 = axes('Parent',figure1,'XTick',[-12 -6 0 6 12],...
+    'Position',[0.606286146026517 0.463982300884956 0.20649289099526 0.372300884955752],...
+    'FontSize',13);
+% Uncomment the following line to preserve the X-limits of the axes
+xlim(axes3,[-9 9]);
+box(axes3,'on');
+hold(axes3,'all');
+
+% Create bar
+bar(xvector3,yvector3,'FaceColor',[0 0 0],'Parent',axes3);
+
+% Create xlabel
+xlabel('相对误差%','FontSize',13);
+
+% Create ylabel
+ylabel('数量','FontSize',13);
+
+% Create title
+title('情形C','FontSize',13);
+

SolveIt.m

@@ -1,6 +1,6 @@
-function XX=SolveIt(deltF,deltG,Init_L,Init_Z,Init_U,Init_W,deltdeltF,ddh,ddg,deltH,Init_Y,Ly,Lz,ContrlCount,Lw,Lul,Luu,RestraintCount,Lx,Balance,PVi,PGi,Busnum,Loadi)
+function 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)
 LxComa=FormLxComa(deltF,deltG,deltH,Init_L,Luu,Lul,Init_Z,Init_Y,Lz,Init_U,Init_W,Lw,Lx);
-H=-deltdeltF+ddh;
+H=-(deltdeltF-ddh-ddgzw);
 %t1=diag(Init_L.\Init_Z-Init_U.\Init_W);
 t1=diag(Init_Z./Init_L-Init_W./Init_U);
 t2=-deltG*( t1 )*deltG';

VoltBar.m

@@ -0,0 +1,79 @@
+function createfigure(xvector1, yvector1, xvector2, yvector2, xvector3, yvector3)
+%CREATEFIGURE(XVECTOR1,YVECTOR1,XVECTOR2,YVECTOR2,XVECTOR3,YVECTOR3)
+%  XVECTOR1:  bar xvector
+%  YVECTOR1:  bar yvector
+%  XVECTOR2:  bar xvector
+%  YVECTOR2:  bar yvector
+%  XVECTOR3:  bar xvector
+%  YVECTOR3:  bar yvector
+
+%  Auto-generated by MATLAB on 11-Feb-2015 17:21:34
+
+% Create figure
+figure1 = figure('Name','电压直方图','Color',[1 1 1]);
+
+% Create axes
+axes1 = axes('Parent',figure1,'XTick',[-0.0002 -0.0001 0 0.0001 0.0002],...
+    'Position',[0.0920853080568723 0.470619469026549 0.20649289099526 0.372300884955752],...
+    'FontSize',13);
+% Uncomment the following line to preserve the Y-limits of the axes
+ylim(axes1,[0 7]);
+box(axes1,'on');
+hold(axes1,'all');
+
+% Create bar
+bar(xvector1,yvector1,'FaceColor',[0 0 0],'Parent',axes1);
+
+% Create xlabel
+xlabel('误差/p.u','FontSize',13);
+
+% Create ylabel
+ylabel('数量','FontSize',13);
+
+% Create title
+title('情形A','FontSize',13);
+
+% Create axes
+axes2 = axes('Parent',figure1,...
+    'Position',[0.350370561164917 0.466814159292035 0.20649289099526 0.372300884955752],...
+    'FontSize',13);
+% Uncomment the following line to preserve the X-limits of the axes
+xlim(axes2,[-0.004 0.004]);
+box(axes2,'on');
+hold(axes2,'all');
+
+% Create bar
+bar(xvector2,yvector2,'FaceColor',[0 0 0],'Parent',axes2);
+
+% Create xlabel
+xlabel('误差/p.u','FontSize',13);
+
+% Create title
+title('情形B','FontSize',13);
+
+% Create ylabel
+ylabel('数量','FontSize',13);
+
+% Create axes
+axes3 = axes('Parent',figure1,...
+    'Position',[0.606286146026517 0.463982300884956 0.20649289099526 0.372300884955752],...
+    'FontSize',13);
+% Uncomment the following line to preserve the X-limits of the axes
+xlim(axes3,[-0.004 0.004]);
+% Uncomment the following line to preserve the Y-limits of the axes
+ylim(axes3,[0 7]);
+box(axes3,'on');
+hold(axes3,'all');
+
+% Create bar
+bar(xvector3,yvector3,'FaceColor',[0 0 0],'Parent',axes3);
+
+% Create xlabel
+xlabel('误差/p.u','FontSize',13);
+
+% Create ylabel
+ylabel('数量','FontSize',13);
+
+% Create title
+title('情形C','FontSize',13);
+

admmatrix.asv

@@ -1,31 +0,0 @@
-function [G,B,GB,Y,r,c,Angle]  = admmatrix(Busnum,Linei,Linej,Liner,Linex,Lineb,Transfori...
-    ,Transforj,Transforr,Transforx,Transfork0,Branchi,Branchb)
-%**************************************************************************
-%  程序功能 : 子函数——形成节点导纳矩阵Y
-%  编    者： 
-%  编制时间：2010.12
-%**************************************************************************
-%% 支路导纳计算
-G = -sparse(Linei,Linej,Liner./(Liner.^2+Linex.^2),Busnum,Busnum) - sparse(Linej,Linei,Liner./(Liner.^2+Linex.^2),Busnum,Busnum);
-G = G - sparse(1:Busnum,1:Busnum,sum(G,2)');                                                                                  % 计算各线路支路电导
-B = sparse(Linei,Linej,Linex./(Liner.^2+Linex.^2),Busnum,Busnum)+sparse(Linej,Linei,Linex./(Liner.^2+Linex.^2),Busnum,Busnum);
-B = B - sparse(1:Busnum,1:Busnum,sum(B,2)')+sparse(Linei,Linei,Lineb,Busnum,Busnum)+sparse(Linej,Linej,Lineb,Busnum,Busnum);
-%% 变压器支路计算
-if Transfori>0 
-    mr = Transforr./(Transforr.^2+Transforx.^2);                                                                                  % 计算变压器支路电导
-    mx = -Transforx./(Transforr.^2+Transforx.^2);                                                                                % 计算变压器支路电纳
-    G = G-sparse(Transfori,Transforj,mr./Transfork0,Busnum,Busnum)-sparse(Transforj,Transfori,mr./Transfork0,Busnum,Busnum)...
-        +sparse(Transfori,Transfori,mr./Transfork0./Transfork0,Busnum,Busnum)+sparse(Transforj,Transforj,mr,Busnum,Busnum);
-    B = B-sparse(Transfori,Transforj,mx./Transfork0,Busnum,Busnum)-sparse(Transforj,Transfori,mx./Transfork0,Busnum,Busnum)...
-        +sparse(Transfori,Transfori,mx./Transfork0./Transfork0,Busnum,Busnum)+sparse(Transforj,Transforj,mx,Busnum,Busnum);
-end
-%% 接地支路计算
-if Branchi>0                                                                 % 判断有无接地支路
-    B = B+sparse(Branchi,Branchi,Branchb,Busnum,Busnum);
-end
-%% 化作极坐标形式
-GB = G+B.*1i;                                                              %将电导，电纳合并，写成复数形式
-Y = abs(GB);                                                                %求节点导纳幅值
-[r,c] = find(Y);
-Angle = angle(GB(GB~=0));                                      %求节点导纳角度
-%Angle=angle(GB);

createLoadProfileFigure.m

@@ -0,0 +1,49 @@
+function createfigure(VertexNormals1, YData1, XData1, Vertices1, Faces1, X1, YMatrix1)
+%CREATEFIGURE(VERTEXNORMALS1,YDATA1,XDATA1,VERTICES1,FACES1,X1,YMATRIX1)
+%  VERTEXNORMALS1:  patch vertexnormals
+%  YDATA1:  patch ydata
+%  XDATA1:  patch xdata
+%  VERTICES1:  patch vertices
+%  FACES1:  patch faces
+%  X1:  vector of x data
+%  YMATRIX1:  matrix of y data
+
+%  Auto-generated by MATLAB on 11-Feb-2015 22:03:00
+
+% Create figure
+figure1 = figure('Color',[1 1 1]);
+
+% Create axes
+axes1 = axes('Parent',figure1,...
+    'XTickLabel',{'6:00','12:00','16:00','24:00','0:00'},...
+    'XTick',[24 48 72 96 1],...
+    'Position',[0.13 0.11 0.373795066413662 0.585744680851064],...
+    'FontName','Times New Roman');
+% Uncomment the following line to preserve the X-limits of the axes
+xlim(axes1,[1 96]);
+hold(axes1,'all');
+
+% Create patch
+patch('Parent',axes1,'VertexNormals',VertexNormals1,'YData',YData1,...
+    'XData',XData1,...
+    'Vertices',Vertices1,...
+    'Faces',Faces1,...
+    'FaceColor',[0.901960784313726 0.901960784313726 0.901960784313726],...
+    'EdgeColor','none');
+
+% Create multiple lines using matrix input to plot
+plot1 = plot(X1,YMatrix1,'Parent',axes1,'LineStyle','--','Color',[0 0 0]);
+set(plot1(2),'DisplayName','边界');
+set(plot1(3),'DisplayName','负荷曲线','LineStyle','-');
+
+% Create xlabel
+xlabel('时刻/时','FontSize',12,'FontName','黑体');
+
+% Create ylabel
+ylabel('负荷/kW','FontSize',12,'FontName','黑体');
+
+% Create legend
+legend1 = legend(axes1,'show');
+set(legend1,...
+    'Position',[0.155913978494625 0.517021276595745 0.108159392789374 0.168794326241135]);
+

func_ddg.m

@@ -1,4 +1,17 @@
-function ddg=func_ddg(PGi,PVi,Busnum,RestraintCount,Loadi,PD,QD)
+function ddgzw=func_ddg(Busnum,ContrlCount,Loadi,Init_Z,Init_W)
-
+CCVbiZW=Init_Z(length(Loadi)*4+Busnum*2+Busnum+length(Loadi)*2+1:length(Loadi)*4+Busnum*2+Busnum+length(Loadi)*2+Busnum)+Init_W(length(Loadi)*4+Busnum*2+Busnum+length(Loadi)*2+1:length(Loadi)*4+Busnum*2+Busnum+length(Loadi)*2+Busnum);
-ddg=0;
+CCPDbiZW=Init_Z(length(Loadi)*4+Busnum*2+Busnum+length(Loadi)*2+Busnum+1:length(Loadi)*4+Busnum*2+Busnum+length(Loadi)*2+Busnum+length(Loadi))+Init_W(length(Loadi)*4+Busnum*2+Busnum+length(Loadi)*2+Busnum+1:length(Loadi)*4+Busnum*2+Busnum+length(Loadi)*2+Busnum+length(Loadi));
+CCQDbiZW=Init_Z(length(Loadi)*4+Busnum*2+Busnum+length(Loadi)*2+Busnum+length(Loadi)+1:end)+Init_W(length(Loadi)*4+Busnum*2+Busnum+length(Loadi)*2+Busnum+length(Loadi)+1:end);
+% ddgzw=[
+% sparse(length(Loadi)*2+Busnum*2,ContrlCount);
+% sparse(Busnum,length(Loadi)*2+Busnum*2),sparse(-2*eye(Busnum))*diag(CCVbiZW),sparse(Busnum,length(Loadi)*2);
+% sparse(length(Loadi),length(Loadi)*2+Busnum*2+Busnum),sparse(-2*eye(length(Loadi)))*diag(CCPDbiZW),sparse(length(Loadi),length(Loadi));
+% sparse(length(Loadi),length(Loadi)*2+Busnum*2+Busnum+length(Loadi)),sparse(-2*eye(length(Loadi)))*diag(CCQDbiZW);
+% ];
+ddgzw=[
+sparse(length(Loadi)*2+Busnum*2,ContrlCount);
+sparse(Busnum,length(Loadi)*2+Busnum*2),sparse(-2*eye(Busnum))*diag(CCVbiZW),sparse(Busnum,length(Loadi)*2);
+sparse(length(Loadi),length(Loadi)*2+Busnum*2+Busnum),sparse(-2*eye(length(Loadi)))*diag(CCPDbiZW),sparse(length(Loadi),length(Loadi));
+sparse(length(Loadi),length(Loadi)*2+Busnum*2+Busnum+length(Loadi)),sparse(-2*eye(length(Loadi)))*diag(CCQDbiZW);
+];
 end

func_ddh.m

@@ -51,7 +51,7 @@ t=[ddPdVdV+ddQdVdV,ddPdVdT+ddQdVdT;
 sizeLoadi=size(Loadi,1)*2;
 ddh=[
     sparse(sizeLoadi,ContrlCount);
-    sparse(2*Busnum,sizeLoadi),-t,sparse(2*Busnum,Busnum);
+    sparse(2*Busnum,sizeLoadi),-t,sparse(2*Busnum,Busnum+length(Loadi)*2);
-    sparse(Busnum,ContrlCount);
+    sparse(Busnum+length(Loadi)*2,ContrlCount);
     ];
 end

func_deltF.m

@@ -1,10 +1,10 @@
 function deltF=func_deltF(PG,QG,PVi,PGi,wPG,wQG,wPD,wQD,PG0,QG0,PD0,PD,QD,QD0,Busnum,Loadi)
-t3=2*wPD.*(PD(Loadi)-PD0(Loadi));
+deltF=[sparse(length(Loadi),1);
-t4=2*wQD.*(QD(Loadi)-QD0(Loadi));
+    sparse(length(Loadi),1);
-deltF=[sparse(t3);
-    sparse(t4);
     sparse(2*Busnum,1);
-    sparse(Busnum,1);
+    sparse(ones(Busnum,1));
+    sparse(ones(length(Loadi),1));
+    sparse(ones(length(Loadi),1))
 ];
 
 end

func_deltG.asv

@@ -1,34 +0,0 @@
-function deltG=func_deltG(Busnum,PVi,PGi,Loadi,PD,QD)
-sizeLoadi=size(Loadi,1);
-dg3_dPD=sparse(1:size(Loadi,1),1:size(Loadi,1),ones(size(Loadi,1),1),size(Loadi,1),size(Loadi,1));
-dg4_dPD=sparse(size(Loadi,1),length(Loadi));
-dg5_dPD=sparse(size(Loadi,1),Busnum);
-dg6_dPD=dg5_dPD;
-dg7_dPD=sparse(sizeLoadi,Busnum);
-%% 
-dg3_dQD=sparse(length(Loadi),length(Loadi));
-dg4_dQD=sparse(1:size(Loadi,1),1:size(Loadi,1),ones(size(Loadi,1),1),size(Loadi,1),size(Loadi,1));
-dg5_dQD=sparse(size(Loadi,1),Busnum);
-dg6_dQD=dg5_dQD;
-dg7_dQD=sparse(sizeLoadi,Busnum);
-%% 
-dg3_dx=sparse(2*Busnum,sizeLoadi);
-dg4_dx=sparse(2*Busnum,length(Loadi));
-dg5_dx=[sparse(1:Busnum,1:Busnum,ones(Busnum,1),Busnum,Busnum);
-    sparse(Busnum,Busnum);
-    ];
-dg6_dx=dg5_dx;
-dg7_dx=sparse(2*Busnum,Busnum);
-%%
-dg3_dvbi=sparse(Busnum,sizeLoadi);
-dg4_dvbi=sparse(Busnum,length(Loadi));
-dg5_dvbi=sparse(Busnum,Busnum);
-dg6_dvbi=sparse(Busnum,Busnum);
-dg6_dvbi=sparse(Busnum,Busnum);
-%%
-deltG=[dg3_dPD,dg4_dPD,dg5_dPD,dg6_dPD;
-       dg3_dQD,dg4_dQD,dg5_dQD,dg6_dQD;
-       dg3_dx,dg4_dx,dg5_dx,dg6_dx;
-       dg3_dvbi,dg4_dvbi,dg5_dvbi,dg6_dvbi;
-];
-end

func_deltG.m

@@ -1,34 +1,107 @@
-function deltG=func_deltG(Busnum,PVi,PGi,Loadi,PD,QD)
+function deltG=func_deltG(Busnum,PVi,PGi,Loadi,PD,QD,Vbi,PDbi,QDbi)
 sizeLoadi=size(Loadi,1);
+%dg3 PD-M*b-t-PD0
 dg3_dPD=sparse(1:size(Loadi,1),1:size(Loadi,1),ones(size(Loadi,1),1),size(Loadi,1),size(Loadi,1));
+%dg32 PD+M*B+t+PD0
+dg32_dPD=sparse(1:size(Loadi,1),1:size(Loadi,1),ones(size(Loadi,1),1),size(Loadi,1),size(Loadi,1));
 dg4_dPD=sparse(size(Loadi,1),length(Loadi));
+dg42_dPD=sparse(size(Loadi,1),length(Loadi));
 dg5_dPD=sparse(size(Loadi,1),Busnum);
-dg6_dPD=dg5_dPD;
+dg6_dPD=sparse(size(Loadi,1),Busnum);
 dg7_dPD=sparse(sizeLoadi,Busnum);
+dgPDbi_dPD=sparse(sizeLoadi,sizeLoadi);
+dgQDbi_dPD=sparse(sizeLoadi,sizeLoadi);
+%互补约束
+dgCCVbi_dPD=sparse(length(Loadi),Busnum);
+dgCCPDbi_dPD=sparse(length(Loadi),length(Loadi));
+dgCCQDbi_dPD=sparse(length(Loadi),length(Loadi));
 %% 
 dg3_dQD=sparse(length(Loadi),length(Loadi));
+dg32_dQD=sparse(length(Loadi),length(Loadi));
 dg4_dQD=sparse(1:size(Loadi,1),1:size(Loadi,1),ones(size(Loadi,1),1),size(Loadi,1),size(Loadi,1));
+dg42_dQD=sparse(1:size(Loadi,1),1:size(Loadi,1),ones(size(Loadi,1),1),size(Loadi,1),size(Loadi,1));
 dg5_dQD=sparse(size(Loadi,1),Busnum);
-dg6_dQD=dg5_dQD;
+dg6_dQD=sparse(size(Loadi,1),Busnum);
 dg7_dQD=sparse(sizeLoadi,Busnum);
+dgPDbi_dQD=sparse(sizeLoadi,sizeLoadi);
+dgQDbi_dQD=sparse(sizeLoadi,sizeLoadi);
+%互补约束
+dgCCVbi_dQD=sparse(length(Loadi),Busnum);
+dgCCPDbi_dQD=sparse(length(Loadi),length(Loadi));
+dgCCQDbi_dQD=sparse(length(Loadi),length(Loadi));
 %% 
 dg3_dx=sparse(2*Busnum,sizeLoadi);
+dg32_dx=sparse(2*Busnum,sizeLoadi);
 dg4_dx=sparse(2*Busnum,length(Loadi));
+dg42_dx=sparse(2*Busnum,length(Loadi));
 dg5_dx=[sparse(1:Busnum,1:Busnum,ones(Busnum,1),Busnum,Busnum);
     sparse(Busnum,Busnum);
     ];
-dg6_dx=dg5_dx;
+dg6_dx=[sparse(1:Busnum,1:Busnum,ones(Busnum,1),Busnum,Busnum);
+    sparse(Busnum,Busnum);
+    ];
 dg7_dx=sparse(2*Busnum,Busnum);
+dgPDbi_dx=sparse(2*Busnum,sizeLoadi);
+dgQDbi_dx=sparse(2*Busnum,sizeLoadi);
+%互补约束
+dgCCVbi_dx=sparse(Busnum*2,Busnum);
+dgCCPDbi_dx=sparse(Busnum*2,length(Loadi));
+dgCCQDbi_dx=sparse(Busnum*2,length(Loadi));
 %%
 dg3_dvbi=sparse(Busnum,sizeLoadi);
+dg32_dvbi=sparse(Busnum,sizeLoadi);
 dg4_dvbi=sparse(Busnum,length(Loadi));
-dg5_dvbi=sparse(Busnum,Busnum);
+dg42_dvbi=sparse(Busnum,length(Loadi));
-dg6_dvbi=sparse(Busnum,Busnum);
+dg5_dvbi=sparse(-eye(Busnum,Busnum));
+dg6_dvbi=sparse(eye(Busnum,Busnum));
 dg7_dvbi=sparse(eye(Busnum,Busnum));
+dgPDbi_dvbi=sparse(Busnum,sizeLoadi);
+dgQDbi_dvbi=sparse(Busnum,sizeLoadi);
+%互补约束
+dgCCVbi_dvbi=sparse(diag(0.1-2*Vbi));%互补约束为(1-Vbi).*Vbi
+dgCCPDbi_dvbi=sparse(Busnum,length(Loadi));
+dgCCQDbi_dvbi=sparse(Busnum,length(Loadi));
 %%
-deltG=[dg3_dPD,dg4_dPD,dg5_dPD,dg6_dPD,dg7_dPD;
+dg3_dPDbi=sparse(-eye(sizeLoadi));
-       dg3_dQD,dg4_dQD,dg5_dQD,dg6_dQD,dg7_dQD;
+dg32_dPDbi=sparse(eye(sizeLoadi));
-       dg3_dx,dg4_dx,dg5_dx,dg6_dx,dg7_dx;
+dg4_dPDbi=sparse(length(Loadi),length(Loadi));
-       dg3_dvbi,dg4_dvbi,dg5_dvbi,dg6_dvbi,dg7_dvbi;
+dg42_dPDbi=sparse(length(Loadi),length(Loadi));
+dg5_dPDbi=sparse(length(Loadi),Busnum);
+dg6_dPDbi=sparse(length(Loadi),Busnum);
+dg7_dPDbi=sparse(length(Loadi),Busnum);
+dgPDbi_dPDbi=sparse(eye(length(Loadi)));
+dgQDbi_dPDbi=sparse(zeros(length(Loadi)));
+%互补约束
+dgCCVbi_dPDbi=sparse(length(Loadi),Busnum);
+dgCCPDbi_dPDbi=sparse(diag(0.1-2*PDbi));
+dgCCQDbi_dPDbi=sparse(length(Loadi),length(Loadi));
+%%
+dg3_dQDbi=sparse(length(Loadi),length(Loadi));
+dg32_dQDbi=sparse(length(Loadi),length(Loadi));
+dg4_dQDbi=sparse(-eye(sizeLoadi));
+dg42_dQDbi=sparse(eye(sizeLoadi));
+dg5_dQDbi=sparse(length(Loadi),Busnum);
+dg6_dQDbi=sparse(length(Loadi),Busnum);
+dg7_dQDbi=sparse(length(Loadi),Busnum);
+dgPDbi_dQDbi=sparse(zeros(length(Loadi)));
+dgQDbi_dQDbi=sparse(eye(length(Loadi)));
+%互补约束
+dgCCVbi_dQDbi=sparse(length(Loadi),Busnum);
+dgCCPDbi_dQDbi=sparse(length(Loadi),length(Loadi));
+dgCCQDbi_dQDbi=sparse(diag(0.1-2*QDbi));
+%%
+% deltG=[dg3_dPD,dg32_dPD,dg4_dPD,dg42_dPD,dg5_dPD,dg6_dPD,dg7_dPD,dgPDbi_dPD,dgQDbi_dPD,dgCCVbi_dPD,dgCCPDbi_dPD,dgCCQDbi_dPD;
+%        dg3_dQD,dg32_dQD,dg4_dQD,dg42_dQD,dg5_dQD,dg6_dQD,dg7_dQD,dgPDbi_dQD,dgQDbi_dQD,dgCCVbi_dQD,dgCCPDbi_dQD,dgCCQDbi_dQD;
+%        dg3_dx,dg32_dx,dg4_dx,dg42_dx,dg5_dx,dg6_dx,dg7_dx,dgPDbi_dx,dgQDbi_dx,dgCCVbi_dx,dgCCPDbi_dx,dgCCQDbi_dx;
+%        dg3_dvbi,dg32_dvbi,dg4_dvbi,dg42_dvbi,dg5_dvbi,dg6_dvbi,dg7_dvbi,dgPDbi_dvbi,dgQDbi_dvbi,dgCCVbi_dvbi,dgCCPDbi_dvbi,dgCCQDbi_dvbi;
+%        dg3_dPDbi,dg32_dPDbi,dg4_dPDbi,dg42_dPDbi,dg5_dPDbi,dg6_dPDbi,dg7_dPDbi,dgPDbi_dPDbi,dgQDbi_dPDbi,dgCCVbi_dPDbi,dgCCPDbi_dPDbi,dgCCQDbi_dPDbi;
+%        dg3_dQDbi,dg32_dQDbi,dg4_dQDbi,dg42_dQDbi,dg5_dQDbi,dg6_dQDbi,dg7_dQDbi,dgPDbi_dQDbi,dgQDbi_dQDbi,dgCCVbi_dQDbi,dgCCPDbi_dQDbi,dgCCQDbi_dQDbi;
+% ];
+deltG=[dg3_dPD,dg32_dPD,dg4_dPD,dg42_dPD,dg5_dPD,dg6_dPD,dg7_dPD,dgPDbi_dPD,dgQDbi_dPD,dgCCVbi_dPD,dgCCPDbi_dPD,dgCCQDbi_dPD;
+       dg3_dQD,dg32_dQD,dg4_dQD,dg42_dQD,dg5_dQD,dg6_dQD,dg7_dQD,dgPDbi_dQD,dgQDbi_dQD,dgCCVbi_dQD,dgCCPDbi_dQD,dgCCQDbi_dQD;
+       dg3_dx,dg32_dx,dg4_dx,dg42_dx,dg5_dx,dg6_dx,dg7_dx,dgPDbi_dx,dgQDbi_dx,dgCCVbi_dx,dgCCPDbi_dx,dgCCQDbi_dx;
+       dg3_dvbi,dg32_dvbi,dg4_dvbi,dg42_dvbi,dg5_dvbi,dg6_dvbi,dg7_dvbi,dgPDbi_dvbi,dgQDbi_dvbi,dgCCVbi_dvbi,dgCCPDbi_dvbi,dgCCQDbi_dvbi;
+       dg3_dPDbi,dg32_dPDbi,dg4_dPDbi,dg42_dPDbi,dg5_dPDbi,dg6_dPDbi,dg7_dPDbi,dgPDbi_dPDbi,dgQDbi_dPDbi,dgCCVbi_dPDbi,dgCCPDbi_dPDbi,dgCCQDbi_dPDbi;
+       dg3_dQDbi,dg32_dQDbi,dg4_dQDbi,dg42_dQDbi,dg5_dQDbi,dg6_dQDbi,dg7_dQDbi,dgPDbi_dQDbi,dgQDbi_dQDbi,dgCCVbi_dQDbi,dgCCPDbi_dQDbi,dgCCQDbi_dQDbi;
 ];
 end

func_deltH.m

@@ -3,6 +3,8 @@ dH_dPD=[sparse(1:size(Loadi,1),Loadi,-ones(size(Loadi,1),1),size(Loadi,1),Busnum
 dH_dQD=[sparse(size(Loadi,1),Busnum) sparse(1:size(Loadi,1),Loadi,-ones(size(Loadi,1),1),size(Loadi,1),Busnum)];
 dH_dx = jacobian_M(Busnum,Volt,Y,Angle,UAngel,r,c);  %形成雅克比矩阵
 dH_dvbi=sparse(Busnum,2*Busnum);
-deltH=[dH_dPD;dH_dQD;dH_dx';dH_dvbi];
+dH_dPDbi=sparse(length(Loadi),2*Busnum);
+dH_dQDbi=sparse(length(Loadi),2*Busnum);
+deltH=[dH_dPD;dH_dQD;dH_dx';dH_dvbi;dH_dPDbi;dH_dQDbi];
 
 end

func_deltdeltF.m

@@ -1,11 +1,3 @@
 function deltdeltF=func_deltdeltF(PVi,wPG,wQG,wPD,wQD,ContrlCount)
-
+deltdeltF=0;
-%ContrlCount=size(PVi,1)+size(PGi,1)+size(wD,1)+Busnum*2; %P,Q,Volt theta这些控制变量数
-C=[wPD' wQD'];
-sizeC=size(C,2);
-diagC=sparse(1:sizeC,1:sizeC,C,sizeC,sizeC);
-deltdeltF=[
-    diagC*2,sparse(sizeC,ContrlCount-sizeC);
-    sparse(ContrlCount-sizeC,ContrlCount);
-    ];
 end

graph/badDataGraph.m

@@ -0,0 +1,16 @@
+x=1:33;
+y=1:33;
+x=[x,15,33];
+y=[y,17,23];
+scatter(x,y);
+xlabel('存在坏数据的节点');
+ylabel('辨识出坏数据的点');
+hold on;
+for I=1:33
+% plot([1,1,1,1]-1+I,[[1:3]-2+I,I]);
+plot([I,I],[0,I]);
+end
+% plot([1,1,1,1]-1+15,[1:3]-2+17);
+% plot([1,1,1,1]-1+33,[1,2,12]-2+23);
+plot([15,15],[0,17]);
+plot([33,33],[0,23]);

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
+

loadprofile.txt

@@ -0,0 +1,96 @@
+57697.32	576.9732	663.51918	490.42722
+55270.921	552.70921	635.6155915	469.8028285
+53275.549	532.75549	612.6688135	452.8421665
+51147.685	511.47685	588.1983775	434.7553225
+49132.866	491.32866	565.027959	417.629361
+47594.538	475.94538	547.337187	404.553573
+45822.028	458.22028	526.953322	389.487238
+44106.141	441.06141	507.2206215	374.9021985
+42717.381	427.17381	491.2498815	363.0977385
+42185.436	421.85436	485.132514	358.576206
+40721.922	407.21922	468.302103	346.136337
+40139.269	401.39269	461.6015935	341.1837865
+38901.576	389.01576	447.368124	330.663396
+37535.238	375.35238	431.655237	319.049523
+37234.754	372.34754	428.199671	316.495409
+36277.407	362.77407	417.1901805	308.3579595
+35705.188	357.05188	410.609662	303.494098
+35364.886	353.64886	406.696189	300.601531
+34644.62	346.4462	398.41313	294.47927
+34811.66	348.1166	400.33409	295.89911
+34473.493	344.73493	396.4451695	293.0246905
+35018.258	350.18258	402.709967	297.655193
+35782.508	357.82508	411.498842	304.151318
+35520.4	355.204	408.4846	301.9234
+37126.94	371.2694	426.95981	315.57899
+37852.994	378.52994	435.309431	321.750449
+39096.988	390.96988	449.615362	332.324398
+40615.023	406.15023	467.0727645	345.2276955
+42599.93	425.9993	489.899195	362.099405
+44509.964	445.09964	511.864586	378.334694
+49173.813	491.73813	565.4988495	417.9774105
+53864.107	538.64107	619.4372305	457.8449095
+60260.101	602.60101	692.9911615	512.2108585
+65195.266	651.95266	749.745559	554.159761
+71841.983	718.41983	826.1828045	610.6568555
+77802.345	778.02345	894.7269675	661.3199325
+88879.392	888.79392	1022.113008	755.474832
+100781.676	1007.81676	1158.989274	856.644246
+111517.636	1115.17636	1282.452814	947.899906
+119678.504	1196.78504	1376.302796	1017.267284
+126802.705	1268.02705	1458.231108	1077.822993
+131988.024	1319.88024	1517.862276	1121.898204
+134124.319	1341.24319	1542.429669	1140.056712
+137238.631	1372.38631	1578.244257	1166.528364
+138993.205	1389.93205	1598.421858	1181.442243
+140492.797	1404.92797	1615.667166	1194.188775
+141981.796	1419.81796	1632.790654	1206.845266
+142679.116	1426.79116	1640.809834	1212.772486
+143682.364	1436.82364	1652.347186	1221.300094
+144526.194	1445.26194	1662.051231	1228.472649
+146089.742	1460.89742	1680.032033	1241.762807
+146519.133	1465.19133	1684.97003	1245.412631
+144834.953	1448.34953	1665.60196	1231.097101
+145784.897	1457.84897	1676.526316	1239.171625
+144503.807	1445.03807	1661.793781	1228.28236
+145438.77	1454.3877	1672.545855	1236.229545
+145203.708	1452.03708	1669.842642	1234.231518
+143780.956	1437.80956	1653.480994	1222.138126
+142645.898	1426.45898	1640.427827	1212.490133
+143466.557	1434.66557	1649.865406	1219.465735
+141922.25	1419.2225	1632.105875	1206.339125
+145841.261	1458.41261	1677.174502	1239.650719
+144923.155	1449.23155	1666.616283	1231.846818
+145430.575	1454.30575	1672.451613	1236.159888
+145334.611	1453.34611	1671.348027	1235.344194
+146149.143	1461.49143	1680.715145	1242.267716
+146482.33	1464.8233	1684.546795	1245.099805
+147876.05	1478.7605	1700.574575	1256.946425
+149817.698	1498.17698	1722.903527	1273.450433
+150415.094	1504.15094	1729.773581	1278.528299
+154442.119	1544.42119	1776.084369	1312.758012
+158010.804	1580.10804	1817.124246	1343.091834
+160774.574	1607.74574	1848.907601	1366.583879
+162786.258	1627.86258	1872.041967	1383.683193
+161608.648	1616.08648	1858.499452	1373.673508
+159454.941	1594.54941	1833.731822	1355.366999
+158027.837	1580.27837	1817.320126	1343.236615
+157666.495	1576.66495	1813.164693	1340.165208
+156460.303	1564.60303	1799.293485	1329.912576
+155215.52	1552.1552	1784.97848	1319.33192
+152007.853	1520.07853	1748.09031	1292.066751
+149345.228	1493.45228	1717.470122	1269.434438
+147650.189	1476.50189	1697.977174	1255.026607
+145332.662	1453.32662	1671.325613	1235.327627
+142143.342	1421.43342	1634.648433	1208.218407
+136555.706	1365.55706	1570.390619	1160.723501
+134474.784	1344.74784	1546.460016	1143.035664
+131788.941	1317.88941	1515.572822	1120.205999
+125144.813	1251.44813	1439.16535	1063.730911
+116612.327	1166.12327	1341.041761	991.2047795
+105557.928	1055.57928	1213.916172	897.242388
+90041.36	900.4136	1035.47564	765.35156
+79091.837	790.91837	909.5561255	672.2806145
+68164.237	681.64237	783.8887255	579.3960145
+61305.02	613.0502	705.00773	521.09267
+55158.221	551.58221	634.3195415	468.8448785

maxBoundError.txt

@@ -0,0 +1,4 @@
+0.010149012	0.000145532
+0.010613753	0.001269963
+0.010686824	0.001667111
+0.01098127	0.002552894

maxError.txt

@@ -0,0 +1,33 @@
+0.001940128	0.001207832
+0.002036094	0.00080888
+0.002099042	0.000927452
+0.002424885	0.000935623
+0.001893664	0.001181456
+0.002018829	0.001257517
+0.00210676	0.000743724
+0.002013494	0.000936731
+0.002013003	0.001095999
+0.001842267	0.000931757
+0.002001811	0.000859772
+0.001866919	0.000876312
+0.001834921	0.00075514
+0.001937562	0.000923904
+0.00189722	0.000822962
+0.001933966	0.000901072
+0.001920252	0.000994934
+0.001747249	0.000857745
+0.002008556	0.000823106
+0.002093043	0.000753621
+0.001971762	0.000883769
+0.002987941	0.00129286
+0.002436167	0.000997165
+0.002037594	0.001015502
+0.002066403	0.000954605
+0.00276576	0.001295141
+0.001881581	0.001290146
+0.00193955	0.001351513
+0.002069665	0.001408808
+0.001936768	0.001467117
+0.001780607	0.001485602
+0.001861682	0.00143783
+0.001955221	0.001473342

modifyadmmatrix.asv

@@ -1,18 +0,0 @@
-function [new_G,new_B,GB,Y,r,c,Angle]  = modifyadmmatrix(ii,jj,G,GG,B,BB)
-%**************************************************************************
-%  程序功能 : 子函数——形成节点导纳矩阵Y
-%  编    者： 
-%  编制时间：2010.12
-%**************************************************************************
-%% 支路导纳计算
-new_G=G;
-new_G(ii,jj)=new_G(ii,jj)-G(ii,jj);
-new_G(jj,ii)=new_G(jj,ii)-G(jj,ii);
-new_G(ii,ii)=new_G(ii,ii)+G(ii,jj);
-new_G(jj,jj)=new_G(jj,jj)+G(ii,jj);
-
-%% 化作极坐标形式
-GB = G+B.*1i;                                                              %将电导，电纳合并，写成复数形式
-Y = abs(GB);                                                                %求节点导纳幅值
-[r,c] = find(Y);
-Angle = angle(GB(GB~=0));                                      %求节点导纳角度

openfile.asv

@@ -1,64 +0,0 @@
-function [Busnum,Balance,PQstandard,Precision,Linei,Linej,Liner,Linex,Lineb,kmax,Transfori ,...
-    Transforj,Transforr,Transforx,Transfork0,Branchi,Branchb,Pointpoweri,PG,QG,PD,QD,PVi,PVu,Gen,GenU,GenL,GenC,CenterA,PGi,PVQU,PVQL] = openfile(FileName)
-%**************************************************************************
-%    程序简介 : 子函数——读取潮流计算所需数据
-%    编    者：
-%    编制时间 ：2010.12
-%**************************************************************************
-data = dlmread(FileName);                                           % 一次读入全部数据
-Busnum= data(1,1);                                                       % 节点数
-PQstandard = data(1,3);                                               % 基准容量
-kmax = data(1,4);                                                           %最大迭代次数
-Precision = data(2,1);                                                    % 精度
-Balance = data(3,2);                                                      % 生成1到节点号的列向量
-CenterA=data(1,5);  %中心参数
-LineNum=data(1,2); %支路数
-%% 各参数矩阵分块
-zeroRow = find(data(:,1)==0);                                     %查找第一列元素为零的行号
-line = data(zeroRow(1)+1:zeroRow(2)-1,:);               % 形成线路参数矩阵
-ground = data(zeroRow(2)+1:zeroRow(3)-1,:);         % 形成对地支路参数矩阵
-tran = data(zeroRow(3)+1:zeroRow(4)-1,:);              % 形成变压器参数矩阵
-buspq = data(zeroRow(4)+1:zeroRow(5)-1,:);          % 形成节点功率参数矩阵
-PV = data(zeroRow(5)+1:zeroRow(6)-1,:);                % 形成pv节点功率参数矩阵
-Gen=data(zeroRow(6)+1:zeroRow(7)-1,:);
-%% 线路参数矩阵分块
-Linei = line(:,2);                                                             % 节点i
-Linej= line(:,3);                                                              % 节点j
-Liner = line(:,4);                                                             % 线路电阻
-Linex = line(:,5);                                                            % 线路电抗
-Lineb = line(:,6);                                                            % b/2
-%% 对地支路参数矩阵
-Branchi = ground(:,1);                                                  % 对地支路节点号
-Branchb = ground(:,2);                                                 % 对地支路的导纳
-%% 变压器参数矩阵
-Transfori = tran(:,2);                                                     % 节点i
-Transforj= tran(:,3);                                                      % 节点j
-Transforr = tran(:,4);                                                     % 变压器电阻
-Transforx= tran(:,5);                                                     % 变压器电抗
-Transfork0 = tran(:,6);                                                  % 变压器变比
-%% 节点功率参数矩阵
-Pointpoweri = buspq(:,1);
-PG=buspq(:,2);                                                             % 发电机有功
-QG=buspq(:,3);                                                             % 发电机无功
-PD=buspq(:,4);                                                              % 负荷有功
-QD=buspq(:,5);                                                             % 负荷无功
-%% pv节点功率参数矩阵
-PVi = PV(:,1); % PV节点的节点号
-PVu = PV(:,2); % PV节点电压
-PVQL=PV(:,3);%PV节点无功下限
-PVQU=PV(:,4); %PV节点无功上限
-%% 发电机参数
-%GenU=Gen(:,[1 5 6]);
-%GenL=Gen(:,[1 7 8]);
-GenC=Gen(:,[1 2:4]);
-t=GenC(:,2);
-GenC(:,2)=GenC(:,4);
-GenC(:,4)=t;
-t=Gen(:,[1 5]);
-%GenL=[t,PVQL(PVi)];
-GenL=t;%有功下界
-t=Gen(:,[1 6]);
-%GenU=[t,PVQU(PVi)];
-GenU=t;%有功上届
-PGi=Gen(:,1);%发电机节点号
-end

openfile2.m

@@ -56,6 +56,12 @@ QD=QD/Base;
 PD=sparse(PD);
 QD=sparse(QD);
 %QD=PD*sqrt(1-.85^2)/.85;
+
+%%
+%增加大负荷
+% PD(22)=PD(22)*65;
+%%
+
 PG=sparse(PG);
 QG=sparse(QG);
 %% pv节点功率参数矩阵

subOPF.m

@@ -0,0 +1,250 @@
+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');
+%pf('D:\Project\青秀降损项目\最小化潮流\最小潮流算例\原始\津头站津视922（3-1）_0.5_120%.txt');
+%pf('D:\Project\最小化潮流\最小潮流算例\仙海919.txt');
+%pf('c:/file31.txt');
+%% 计算功率因数
+Loadi=QD~=0 | PD~=0;
+PF=sqrt(PD(Loadi).^2./(QD(Loadi).^2+PD(Loadi).^2));
+%%
+% Volt;
+% UAngel*180/3.1415926;
+%% 通过潮流计算PG
+AngleIJ=sparse(r,c,UAngel(r)-UAngel(c)-Angle',Busnum,Busnum);
+PGBal=PD+diag(Volt)*Y.*cos(AngleIJ)*Volt';
+QGBal=QD+diag(Volt)*Y.*sin(AngleIJ)*Volt';
+%% 真实值
+% PG0=PG;
+% QG0=QG;
+% PD0=PD;
+% QD0=QD;
+% Volt0=Volt;
+% UAngel0=UAngel;
+RealPD=PD;
+RealQD=QD;
+rUAngel=UAngel;
+rVolt=Volt;
+%%
+PG0(Balance)=PGBal(Balance);
+PG(Balance)=PGBal(Balance);
+QG0(Balance)=QGBal(Balance);
+QG0(PVi)=QGBal(PVi);
+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;
+
+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=170;
+Precision=Precision/1;
+%% 加误差
+%找DG
+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);
+% 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(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);
+%% 不测量的数据
+% mPD(22)=mPD(22)/65;
+
+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);
+
+
+% lPD(22)=0.17*mPD(22);%15%
+% uPD(22)=0.17*mPD(22);
+% lQD(22)=0.17*mQD(22);
+% uQD(22)=0.17*mQD(22);
+
+
+%电压伪量测错误
+mVolt(noMeasurei)=1;
+
+lVolt(noMeasurei)=0.12*mVolt(noMeasurei);%电压在0.93~1.07
+uVolt(noMeasurei)=0.02*mVolt(noMeasurei);
+
+
+
+%错误数据
+%mVolt(2)=5;
+% bigM=0.000003;
+Vbi=sparse(0.5*ones(Busnum,1));
+% Vbi(2)=1;
+PDbi=sparse(0.5*ones(length(Loadi),1));
+QDbi=sparse(0.5*ones(length(Loadi),1));
+eps=1;
+% Gap=0;
+% 第一遍，算连续的值
+fprintf('第1次迭代，算连续量。\n');
+tic;
+while(abs(Gap)>Precision*10)
+    if KK>kmax
+        break;
+    end
+    plotGap(KK+1)=Gap;
+    Init_u=Gap/2/RestraintCount*CenterA;
+    AngleIJMat=0;
+    %% 开始计算OPF
+    %% 形成等式约束的雅克比
+    deltH=func_deltH(Busnum,Volt,PVi,Y,PGi,UAngel,r,c,Angle,Loadi);
+    %% 形成不等式约束的雅克比
+    deltG=func_deltG(Busnum,PVi,PGi,Loadi,PD,QD,Vbi,PDbi,QDbi);
+    %%
+    L_1Z=diag(Init_Z./Init_L);
+    U_1W=diag(Init_W./Init_U);
+    %% 形成海森阵
+    deltdeltF=func_deltdeltF(PVi,wPG,wQG,wPD,wQD,ContrlCount);
+    %% 形成ddHy
+    ddh=func_ddh(Volt,Init_Y,Busnum,PVi,PGi,Y,UAngel,r,c,Angle,Loadi,ContrlCount);
+    %%  开始构建ddg
+    ddgzw=func_ddg(Busnum,ContrlCount,Loadi,Init_Z,Init_W);
+    %%  开始构建deltF
+    deltF=func_deltF(PG,QG,PVi,PGi,wPG,wQG,wPD,wQD,PG0,QG0,PD0,PD,QD,QD0,Busnum,Loadi);
+    %% 形成方程矩阵
+    Luu=Init_U'.*Init_W'+Init_u*ones(RestraintCount,1);
+    Lul=Init_L'.*Init_Z'-Init_u*ones(RestraintCount,1);
+    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;
+    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);
+    %% 开始解方程
+    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);
+    %%取各分量
+    [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
+fprintf('第一遍计算时间:%f\n',toc)
+% 第二遍，算离散
+% [~,~,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=1;
+% KK=0;
+% 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);%与学姐一致
+tic;
+fprintf('第2次迭代，算离散量。\n');
+while(abs(Gap)>Precision*10)
+    if KK>kmax
+                    break;
+    end
+%             Vbi(23)=0.1;
+    plotGap(KK+1)=Gap;
+    Init_u=Gap/2/RestraintCount*CenterA;
+    AngleIJMat=0;
+    %% 开始计算OPF
+    %% 形成等式约束的雅克比
+    deltH=func_deltH(Busnum,Volt,PVi,Y,PGi,UAngel,r,c,Angle,Loadi);
+    %% 形成不等式约束的雅克比
+    deltG=func_deltG(Busnum,PVi,PGi,Loadi,PD,QD,Vbi,PDbi,QDbi);
+    %%
+    L_1Z=diag(Init_Z./Init_L);
+    U_1W=diag(Init_W./Init_U);
+    %% 形成海森阵
+    deltdeltF=func_deltdeltF(PVi,wPG,wQG,wPD,wQD,ContrlCount);
+    %% 形成ddHy
+    ddh=func_ddh(Volt,Init_Y,Busnum,PVi,PGi,Y,UAngel,r,c,Angle,Loadi,ContrlCount);
+    %%  开始构建ddg
+    ddgzw=func_ddg(Busnum,ContrlCount,Loadi,Init_Z,Init_W);
+    %%  开始构建deltF
+    deltF=func_deltF(PG,QG,PVi,PGi,wPG,wQG,wPD,wQD,PG0,QG0,PD0,PD,QD,QD0,Busnum,Loadi);
+    %% 形成方程矩阵
+    Luu=Init_U'.*Init_W'+Init_u*ones(RestraintCount,1);
+    Lul=Init_L'.*Init_Z'-Init_u*ones(RestraintCount,1);
+    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);
+    Ly=Mat_H;
+    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);
+    %% 开始解方程
+%         Vbi(Vbi>0.003)=.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);
+    [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
+fprintf('第二遍计算时间:%f\n',toc)
+%% 计算最大不平衡量
+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)
+%判断是否收敛
+if min(abs(Vbi-.1))>5e-005
+    isConverge=0;
+end
+if KK>=kmax 
+    isConverge=0;
+end
+toc
+end