将电压加入目标函数中。

2013-04-04 22:50:20 +08:00
106 changed files with 1202 additions and 1934 deletions
--- a/.gitignore
+++ b/.gitignore
@@ -1,2 +0,0 @@
-*.asv
-*.7z
--- a/0PD0.mat
+++ b/0PD0.mat
--- a/0QD0.mat
+++ b/0QD0.mat
--- a/12PD0.mat
+++ b/12PD0.mat
--- a/12QD0.mat
+++ b/12QD0.mat
--- a/16PD0.mat
+++ b/16PD0.mat
--- a/16QD0.mat
+++ b/16QD0.mat
--- a/20PD0.mat
+++ b/20PD0.mat
--- a/20QD0.mat
+++ b/20QD0.mat
--- a/4PD0.mat
+++ b/4PD0.mat
--- a/4QD0.mat
+++ b/4QD0.mat
--- a/8PD0.mat
+++ b/8PD0.mat
--- a/8QD0.mat
+++ b/8QD0.mat
--- a/@ForThesis1/AddPDQDPGQG.m
+++ b/@ForThesis1/AddPDQDPGQG.m
--- a/@ForThesis1/ForThesis.m
+++ b/@ForThesis1/ForThesis.m
--- a/@ForThesis/MaxBranchDeviation.asv
+++ b/@ForThesis/MaxBranchDeviation.asv
@@ -0,0 +1,41 @@
+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
+
--- a/@ForThesis1/MaxBranchDeviation.m
+++ b/@ForThesis1/MaxBranchDeviation.m
@@ -29,8 +29,8 @@ 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);%¹À¼ÆÖµ
+[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;% 太小的值不计算
--- a/@ForThesis/MaxDeviation.asv
+++ b/@ForThesis/MaxDeviation.asv
@@ -0,0 +1,7 @@
+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
--- a/@ForThesis1/MaxDeviation.m
+++ b/@ForThesis1/MaxDeviation.m
--- a/@ForThesis1/MeanPD.m
+++ b/@ForThesis1/MeanPD.m
--- a/@ForThesis1/MeanPG.m
+++ b/@ForThesis1/MeanPG.m
--- a/@ForThesis1/MeanQD.m
+++ b/@ForThesis1/MeanQD.m
--- a/@ForThesis1/MeanQG.m
+++ b/@ForThesis1/MeanQG.m
--- a/@ForThesis1/PercentOfPass.m
+++ b/@ForThesis1/PercentOfPass.m
--- a/@ForThesis1/SqareDeviation.m
+++ b/@ForThesis1/SqareDeviation.m
--- a/@ForThesis1/StatBranchDeviation.m
+++ b/@ForThesis1/StatBranchDeviation.m
--- a/@ForThesis/StatDeviation.asv
+++ b/@ForThesis/StatDeviation.asv
@@ -0,0 +1,15 @@
+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
+
--- a/@ForThesis1/StatDeviation.m
+++ b/@ForThesis1/StatDeviation.m
--- a/AngelBar.m
+++ b/AngelBar.m
@@ -1,83 +0,0 @@
-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);
-
--- a/CalTimeBar.m
+++ b/CalTimeBar.m
@@ -1,73 +0,0 @@
-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');
-
--- a/CreateMaxErrorFigure.m
+++ b/CreateMaxErrorFigure.m
@@ -1,40 +0,0 @@
-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','黑体');
-
--- a/DeviationFigure.m
+++ b/DeviationFigure.m
@@ -1,115 +0,0 @@
-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]);
-
--- a/DrawGap.asv
+++ b/DrawGap.asv
@@ -0,0 +1,6 @@
+function DrawGap(plotGap)
+x=find(plotGap);
+ts=size(x,2);
+
+
+end
--- a/DrawLoadProfile.m
+++ b/DrawLoadProfile.m
@@ -1,20 +0,0 @@
-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
-
--- a/FormAA.asv
+++ b/FormAA.asv
@@ -0,0 +1,12 @@
+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
--- a/FormAlphaP.asv
+++ b/FormAlphaP.asv
@@ -0,0 +1,7 @@
+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
--- a/FormG.asv
+++ b/FormG.asv
@@ -0,0 +1,11 @@
+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
--- a/FormG.m
+++ b/FormG.m
@@ -1,17 +1,8 @@
-function Mat_G=FormG(Volt,PD,QD,Loadi,bigM,mVolt,rVolt,sigma,Vbi,PDbi,QDbi,mPD,mQD,uPD,lPD,uQD,lQD,uVolt,lVolt)
+function Mat_G=FormG(Volt,PVi,PGi,PG,QG,PD,QD,Loadi)

 Mat_G=[
-    sparse(PD(Loadi))-mPD(Loadi)-bigM*PDbi-uPD(Loadi);
-    sparse(PD(Loadi))-mPD(Loadi)+bigM*PDbi+lPD(Loadi);
-    sparse(QD(Loadi))-mQD(Loadi)-bigM*QDbi-uQD(Loadi);
-    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;
+    sparse(PD(Loadi));
+    sparse(QD(Loadi));
+    Volt';
    ];
 end
--- a/FormH.asv
+++ b/FormH.asv
@@ -0,0 +1,31 @@
+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
--- a/FormH.m
+++ b/FormH.m
@@ -1,6 +1,14 @@
-function Mat_H=FormH(Busnum,Volt,PG,PD,QG,QD,Y,UAngel,r,c,Angle,Loadi)
+function Mat_H=FormH(Busnum,Volt,PG,PD,QG,QD,Y,UAngel,r,c,Angle,QD_NON_ZERO,QD_NON_ZERO_IND,Loadi)
+%% 
+%QDcos=textread('D:\Project\×îÐ¡»¯³±Á÷\×îÐ¡³±Á÷ËãÀý\ÏÉº£919PDQDglys.txt');
+%QD(QD~=0)=PD(QD~=0)./tan(QDcos);
+%QD(QD_NON_ZERO_IND)=QD_NON_ZERO;
+%% 
+%PD(Loadi)=QD(Loadi)./tan(acos(0.98));
 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
--- a/FormLw.asv
+++ b/FormLw.asv
@@ -0,0 +1,25 @@
+function Lw=FormLw(Mat_G,Init_U,GenU,Busnum,PVQU,PD,PD0,QD0,Loadi,KK,PF,noDataTransCapacity)
+KK=999;
+%PU=GenU(:,2);%发电机有功上界
+PU=5*ones(length(GenU(:,2)),1);
+%QU=PVQU(:,1);%发电机无功上界
+QU=5*ones(length(PVQU(:,1)),1);
+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;
+%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);
+% QDU(QDU==0)=0.200;
+% PF=0.85;
+% QDU=1.0*PD(Loadi).*sqrt(1 -PF.^2)./PF;
+t1=([PU',QU',PDU',QDU',VoltU])';
+t2=Mat_G+Init_U'-t1;
+Lw=t2;
+
+end
--- a/FormLw.m
+++ b/FormLw.m
@@ -1,28 +1,22 @@
-function Lw=FormLw(Mat_G,Init_U,GenU,Busnum,PVQU,PD,PD0,QD0,Loadi,KK,PF,eps)
+function Lw=FormLw(Mat_G,Init_U,GenU,Busnum,PVQU,PD,PD0,QD0,Loadi,KK,PF,noDataTransCapacity)
 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;
-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);
 QDU(QDU==0)=0.400;
-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));
-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)])';
-t1=[t1;
-    sparse(eps*ones(Busnum,1));
-   sparse(eps*ones(length(Loadi),1));
-    sparse(eps*ones(length(Loadi),1));
-    ];
+% PF=0.85;
+% QDU=1.0*PD(Loadi).*sqrt(1 -PF.^2)./PF;
+t1=([PDU',QDU',VoltU])';
 t2=Mat_G+Init_U'-t1;
 Lw=t2;

--- a/FormLxComa.m
+++ b/FormLxComa.m
@@ -4,9 +4,8 @@ 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));
-% t7=t6*t5;%
-t7=diag(Init_U)\t5;
+t6=inv(diag(Init_U));
+t7=t6*t5;%
 t8=deltG*(t4+t7);%%
 LxComa=Lx+t8;
 end
--- a/FormLz.asv
+++ b/FormLz.asv
@@ -0,0 +1,20 @@
+function Lz=FormLz(Mat_G,Init_L,GenL,Busnum,PVQL,PD,PD0,QD0,Loadi,KK,PF)
+KK=999;
+
+VoltL=(0.9)*ones(1,Busnum);
+%VoltL=-10*ones(1,Busnum);
+PDL=PD0(Loadi);
+PDL(PDL>0)=0.700*PDL(PDL>0);
+PDL(PDL<0)=1.300*PDL(PDL<0);
+PDL(PDL==0)=-0.400;
+%PDL=-10*ones(length(Loadi),1);
+QDL=QD0(Loadi);
+QDL(QDL>0)=0.700*QDL(QDL>0);
+QDL(QDL<0)=1.300*QDL(QDL<0);
+QDL(QDL==0)=-0.400;
+% QDL=0*PD(Loadi).*sqrt((1-PF.^2))./PF;
+t1=([PDL',QDL',VoltL])';
+t2=Mat_G-Init_L'-t1;
+Lz=t2;
+
+end
--- a/FormLz.m
+++ b/FormLz.m
@@ -1,4 +1,4 @@
-function Lz=FormLz(Mat_G,Init_L,GenL,Busnum,PVQL,PD,PD0,QD0,Loadi,KK,PF,eps)
+function Lz=FormLz(Mat_G,Init_L,GenL,Busnum,PVQL,PD,PD0,QD0,Loadi,KK,PF)
 KK=999;

 VoltL=(0.9)*ones(1,Busnum);
@@ -7,26 +7,13 @@ PDL=PD0(Loadi);
 PDL(PDL>0)=0.800*PDL(PDL>0);
 PDL(PDL<0)=1.200*PDL(PDL<0);
 PDL(PDL==0)=-0.400;
-realPD=PD0(Loadi);
-indPD=find(realPD>0);
-PDL(indPD(3:12:end))=0.95*realPD(indPD(3:12:end));
-PDL(indPD(9:12:end))=0.45*realPD(indPD(9:12:end));
 %PDL=-10*ones(length(Loadi),1);
 QDL=QD0(Loadi);
 QDL(QDL>0)=0.800*QDL(QDL>0);
 QDL(QDL<0)=1.200*QDL(QDL<0);
 QDL(QDL==0)=-0.400;
-realQD=QD0(Loadi);
-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=([-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));
-];
+t1=([PDL',QDL',VoltL])';
 t2=Mat_G-Init_L'-t1;
 Lz=t2;

--- a/FormYY.asv
+++ b/FormYY.asv
@@ -0,0 +1,5 @@
+function FormYY(Init_L,Lul.Lz,Init_U,Luu,Lw,Lz,LxComa)
+t=[
+   -    
+]
+end
--- a/JSMJZM.asv
+++ b/JSMJZM.asv
@@ -0,0 +1,12 @@
+JSM=[0.0575 0.051 0.0029];
+JZM=[0.003 0.02629 0.0582];
+plot(1:3,JSM,'k');
+axis([0 4 0 0.07])
+text(1,JSM(1)+0.001,'0.0575')
+text(2,JSM(2)+0.001,'0.051')
+text(3,JSM(3)+0.001,'0.0029')
+hold on 
+plot(1:3,JZM,'k--');
+text(1,JZM(1)+0.001,'0.003')
+text(2,JZM(2)+0.001,'0.02629')
+text(3,JZM(3)+0.001,'0.0029')
--- a/JSMJZM.m
+++ b/JSMJZM.m
@@ -1,18 +1,18 @@
-JSM=[0.070244 0.0533 0.004379];
-JZM=[0.002275 0.047311 0.070318];
-for I=0.01:0.01:0.08
+JSM=[0.0575 0.051 0.0029];
+JZM=[0.0003 0.02629 0.0582];
+for I=0.01:0.01:0.06
    line([0.5 3.5],[I I],'Color',[220 220 220]/255)
 end
 hold on
 plot(1:3,JSM,'k');
 ylabel('Í³¼ÆÎó²î')
-xlabel('ÇéÐÎ')
-axis([0.5 3.5 0 0.075])
-text(1,JSM(1)+0.001,'0.0702')
-text(2,JSM(2)+0.001,'0.0533')
-text(3,JSM(3)+0.001,'0.0043')
+xlabel('Çé¿ö')
+axis([0.5 3.5 0 0.065])
+text(1,JSM(1)+0.001,'0.0575')
+text(2,JSM(2)+0.001,'0.0510')
+text(3,JSM(3)+0.001,'0.0029')
 hold on 
 plot(1:3,JZM,'k--');
-text(1-.3,JZM(1)+0.001,'0.0023')
-text(2-0.3,JZM(2)-0.001,'0.0473')
-text(3,JZM(3)-0.001,'0.0703')
+text(1-.3,JZM(1)+0.001,'0.0003')
+text(2,JZM(2)-0.001,'0.0263')
+text(3,JZM(3)-0.001,'0.0029')
--- a/JSNZN.asv
+++ b/JSNZN.asv
@@ -0,0 +1,14 @@
+JSN=[0.1213 0.1092 0.0587];
+JZN=[0.0971 0.1314 0.1408];
+plot(1:3,JSM,'k');
+ylabel('Öľ')
+xlabel('Çéżö')
+axis([0.5 3.5 0 0.065])
+text(1,JSM(1)+0.001,'0.1213')
+text(2,JSM(2)+0.001,'0.0510')
+text(3,JSM(3)+0.001,'0.0029')
+hold on 
+plot(1:3,JZM,'k--');
+text(1-.3,JZM(1)+0.001,'0.0003')
+text(2,JZM(2)-0.001,'0.0263')
+text(3,JZM(3)-0.001,'0.0029')
--- a/Lineloss.asv
+++ b/Lineloss.asv
@@ -0,0 +1,38 @@
+%% 计算线损
+function [dispLineloss dispTransloss]=Lineloss(Linei,Linej,Liner,Linex,Lineb2,Transi,Transj,Transr,Transx,Branchi,Branchg,Branchb,k0,Volt,Angle)
+%format long
+% fprintf('功率为有名值\n');
+% fprintf('节点号\t节点号\t有功损耗 MW\t无功损耗 MVar')
+cmpVolt=Volt'.*cos(Angle')+1i*Volt'.*sin(Angle');
+% cmpVolt=cmpVolt';
+y0=1i*Lineb2;
+%yj0=1i*standardinput(:,7);
+yij=1./(Liner+1i*Linex);
+%% 线路损耗
+Sij=cmpVolt(Linei).*conj(  cmpVolt(Linei) .* y0 + ( cmpVolt(Linei)- cmpVolt(Linej) ).*yij );
+Sji=cmpVolt(Linej).*conj(  cmpVolt(Linej) .*y0 + ( cmpVolt(Linej)- cmpVolt(Linei) ).*yij );
+%Sij1==Sij2
+deltLineS=Sij+Sji;
+
+%% 另一种计算方式begin
+ss=1*(Volt(Linei)'.^2.*abs(yij).*cos( angle(yij) ) -Volt(Linei)'.*Volt(Linej)'.*cos( Angle(Linei)' - Angle(Linej)' - angle(yij)).*abs(yij));
+ss=(Volt(Linei)'.^2+Volt(Linej)'.^2).*abs(yij).*cos(angle(yij))-2*Volt(Linei)'.*Volt(Linej)'.*cos( Angle(Linei)' - Angle(Linej)').*cos( - angle(yij)).*abs(yij);
+ss=Volt(Linei)'.*Volt(Linej)'.*abs(yij).*cos()
+%% 另一种计算方式end
+dispLineloss=[Linei Linej real(deltLineS)*100 imag(deltLineS)*100];
+%full(dispLineloss)
+% dispLineloss=sortrows(dispLineloss,-3);
+full(dispLineloss);
+%% 以下是变压器损耗
+yij=1./(Transr+1i*Transx);
+Sij=cmpVolt(Transi)./k0.*conj(  ( cmpVolt(Transi)./k0- cmpVolt(Transj) ).*yij );
+Sji=cmpVolt(Transj).*conj( ( cmpVolt(Transj)- cmpVolt(Transi)./k0 ).*yij );
+deltTransS=Sij+Sji;
+%% 接地支路损耗
+% 没有考虑变压器变比
+deltTransS =deltTransS+sum(cmpVolt(Branchi).*conj((cmpVolt(Branchi).*(Branchg+1j*Branchb))));
+%% 
+dispTransloss=[Transi Transj real(deltTransS)*100 imag(deltTransS)*100];
+% dispTransloss=sortrows(dispTransloss,-3);
+full(dispTransloss);
+end
--- a/LoadProfile.m
+++ b/LoadProfile.m
@@ -1,18 +0,0 @@
-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)
--- a/MaxBoundErrorFigure.m
+++ b/MaxBoundErrorFigure.m
@@ -1,9 +0,0 @@
-function MaxBoundErrorFigure()
-figure();
-data=dlmread('maxBoundError.txt');
-% data(1,:)=0;
-maxV=data(:,1);
-maxA=data(:,2);
-% bar([maxV';maxA']');
-CreateMaxErrorFigure([maxV';maxA']');
-end
--- a/MaxErrorFigure.m
+++ b/MaxErrorFigure.m
@@ -1,9 +0,0 @@
-function MaxErrorFigure()
-figure();
-data=dlmread('maxError.txt');
-data(1,:)=0;
-maxV=data(:,1);
-maxA=data(:,2);
-% bar([maxV';maxA']');
-CreateMaxErrorFigure([maxV';maxA']');
-end
--- a/Modification.asv
+++ b/Modification.asv
@@ -0,0 +1,18 @@
+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
--- a/Modification.m
+++ b/Modification.m
@@ -1,16 +1,8 @@
-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)
+function [Init_Z,Init_L,Init_W,Init_U,Init_Y,PG,QG,Volt,UAngel,PD,QD]=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)
 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';
@@ -32,7 +24,4 @@ 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: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
--- a/OPF.asv
+++ b/OPF.asv
@@ -0,0 +1,137 @@
+tic
+clc
+clear
+thesis=ForThesis(4,8);
+for II=1:4
+    [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,Transfork0]= ...
+        pf('c:/newFIle.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;
+    %%
+    PG0(Balance)=PGBal(Balance);
+    PG(Balance)=PGBal(Balance);
+    QG0(Balance)=QGBal(Balance);
+    QG0(PVi)=QGBal(PVi);
+    QG(PVi)=QGBal(PVi);
+    %% 真实值
+    RealPG=PG0;
+    RealQG=QG0;
+    RealPD=PD0;
+    RealQD=QD0;
+    %%
+    %PGi=zeros(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);
+    Gap=(Init_L*Init_Z'-Init_U*Init_W');
+    KK=0;
+    plotGap=zeros(1,60);
+    ContrlCount=size(Loadi,1)*2+Busnum*2;
+    kmax=60;
+    %% 20120523 临时
+    QD_NON_ZERO=QD(PD==0 & QD~=0);
+    QD_NON_ZERO_IND=find(PD==0 & QD~=0);
+    %%
+    Precision=Precision/10;
+    
+    %% 加误差
+    PD0(Loadi)=PD0(Loadi).*(1+normrnd(0,0.05,length(Loadi),1));
+    QD0(Loadi)=QD0(Loadi).*(1+normrnd(0,0.05,length(Loadi),1));
+%     PG0(PGi)=PG0(PGi).*(1+normrnd(0,0.01,length(PGi),1));
+%     QG0(PVi)=QG0(PVi).*(1+normrnd(0,0.01,length(PVi),1));
+    %% 读变压器容量
+    %[noDataTransNum noDataTransCapacity noDataTransPowerFactor]=ReadNoDataTrans('C:/b/东际911_2751267_2012-09-05/iPso_东际911_2751267_2012-09-05_变压器无负载.txt');
+    noDataTransCapacity=0;
+    
+    while(abs(Gap)>Precision)
+        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);
+        %%
+        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);
+        Mat_H=FormH(Busnum,Volt,PG,PD,QG,QD,Y,UAngel,r,c,Angle,QD_NON_ZERO,QD_NON_ZERO_IND,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,noDataTransCapacity);
+        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]=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);
+        Gap=(Init_L*Init_Z'-Init_U*Init_W');
+        KK=KK+1;
+    end
+    fprintf('迭代次数%d\n',KK);
+    fprintf('目标值%f\n',full(ObjectiveFun(PG,PG0,PGi,QG,QG0,PVi,PD,PD0,QD,QD0,wPG,wQG,wPD,wQD,Loadi)));
+    % DrawGap(plotGap);
+    %%
+    %Volt=full(Volt');
+    %PD=full(PD);
+    %% 统计PD误差
+    % absPDLoad=abs(  (PD(Loadi)-PD0(Loadi))./PD0(Loadi)  );
+    absPDLoad=abs(  (PD(Loadi)-RealPD(Loadi))./RealPD(Loadi)  );
+    maxPDError=max(absPDLoad(absPDLoad<10))
+    Loadi(maxPDError==absPDLoad)
+    absQDLoad=abs(  (QD(Loadi)-RealQD(Loadi))./RealQD(Loadi)  );
+    maxQDError=max(absQDLoad(absQDLoad<10))
+    Loadi(maxQDError==absQDLoad)
+    disp('index');
+    %Loadi(absPDLoad==maxPDError);
+    %% 计算总线损
+    totalLoss=(sum(PG)-sum(PD(Loadi)))*100;
+    fprintf('总的损耗为%f(MW 有名值)\n',full(totalLoss));
+    fprintf('线损率为 %f\n',full(totalLoss/sum(PG)));
+    %% 计算各线损
+    %Lineloss(Linei,Linej,Liner,Linex,Lineb,Transfori,Transforj,Transforr,Transforx,Transfork0,Volt,UAngel);
+    thesis=thesis.AddPDQDPGQG(PD(Loadi),QD(Loadi),PG(Balance),QG(PVi));
+    
+end
+PD(Loadi)=thesis.MeanPD();
+QD(Loadi)=thesis.MeanQD();
+PG(Balance)=thesis.MeanPG();
+QG(PVi)=thesis.MeanQG();
+thesis.MaxDeviation(RealPG(Balance),RealQG(PVi),RealPD(Loadi),RealQD(Loadi))
+thesis.StatDeviation(RealPG(Balance),RealQG(PVi),RealPD(Loadi),RealQD(Loadi));
+toc
+
--- a/OPF.m
+++ b/OPF.m
@@ -1,406 +1,148 @@
+tic
 clc
 clear
-close all
-arrayA=zeros(21,10);
-sumCaseA_SE=0;
-sumCaseB_SE=0;
-sumCaseC_SE=0;
-badDataResult=zeros(10,33);
-badDataLocation=zeros(34,10);
-nodeMaxDVolt=zeros(33,1);
-nodeMaxDVAngle=zeros(33,1);
-for badDataNode=1:1
-    loopN=1;
-    maxDVolt=0;
-    maxDVAngle=0;
-    while 1
-        close
-        [kmax,Precision,UAngel,Volt,Busnum,PVi,PVu,Balance,Y,Angle,P0,Q0,r,c,GB,Linei,Linej,Transfori,Transforj,GenU,GenL,GenC,PG,QG,PD,QD,CenterA,PGi,PVQU,PVQL,Liner,Linex,Lineb,Transforr,Transforx,Branchi,Branchg,Branchb,Transfork0]= ...
-            pf('E:\算例\feeder33\feeder33.txt');
-        sigma=0.01;
-        RealPD=PD;
-        RealQD=QD;
-        rVolt=Volt;
-        Loadi=find(PD~=0);
-        PD0=sparse(Busnum,1);
-        QD0=sparse(Busnum,1);
-        PD0(Loadi)=RealPD(Loadi).*(1+normrnd(0,sigma,length(Loadi),1));
-        QD0(Loadi)=RealQD(Loadi).*(1+normrnd(0,sigma,length(Loadi),1));
-        mVolt=rVolt.*(1+normrnd(0,sigma,length(rVolt),1))';
-        
-        %加载保存的变量
-        PD0=load('PD0');
-        PD0=PD0.PD0;
-        QD0=load('QD0');
-        QD0=QD0.QD0;
-        mVolt=load('mVolt');
-        mVolt=mVolt.mVolt;
-        
-%         mVolt(3)=rVolt(3)*(1-sigma*6);
-        %% 画Case A的图
-        % figure('Color',[1 1 1]);
-        [Busnum,Loadi,Volt,PD,QD,rVolt,UAngel,RealPD,RealQD,rUAngel,Vbi,PDbi,QDbi,plotGapA,isConverge]=subOPF([],PD0,QD0,mVolt,sigma);%全部有
-        %         badDataResult(I,badDataNode)=sum(Vbi);
-        %         badDataLocation(1:33,I)=Vbi;
-        %         badDataLocation(34,I)=sum(abs((rVolt-Volt)./rVolt./length(rVolt)))+sum(abs( (UAngel(2:33)-rUAngel(2:33))./rUAngel(2:33)./length(rUAngel(2:33))));
-                break;
-        if isConverge==0
-            continue;
-        end
-        % 加最大偏差
-        maxDVolt_=max(abs((rVolt-Volt)));
-        if maxDVolt_>maxDVolt
-            maxDVolt=maxDVolt_;
-        end
-        maxDVAngle_=max(abs((UAngel(2:33)-rUAngel(2:33))));
-        if maxDVAngle_>maxDVAngle
-            maxDVAngle=maxDVAngle_;
-        end
-        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;
+%% 存在问题
+% 变压器变比的位置没有考虑，由于现在用的变比都是1，所以没有影响。 20130123
+%% 
+thesis=ForThesis(1,45);
+for II=1: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:/算例/柳金Ⅰ926_21671693_2012-09-06/newFIle20 - 计算用2 - (最终用这个做20点的计算).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;
+    %%
+    PG0(Balance)=PGBal(Balance);
+    PG(Balance)=PGBal(Balance);
+    QG0(Balance)=QGBal(Balance);
+    QG0(PVi)=QGBal(PVi);
+    QG(PVi)=QGBal(PVi);
+    %% 真实值
+    RealPG=PG0;
+    RealQG=QG0;
+    RealPD=PD0;
+    RealQD=QD0;
+    %%
+    %PGi=zeros(1,1);
+    [Volt,UAngel,Init_Z,Init_W,Init_L,Init_U,Init_Y,PG,QG,RestraintCount,wPG,wQG,wPD,wQD,wVolt,PD,PD0,QD,randPDind,Loadi]=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');
+    KK=0;
+    plotGap=zeros(1,60);
+    ContrlCount=size(Loadi,1)*2+Busnum*2;
+    kmax=60;
+    %% 20120523 临时
+    QD_NON_ZERO=QD(PD==0 & QD~=0);
+    QD_NON_ZERO_IND=find(PD==0 & QD~=0);
+    %%
+    Precision=Precision/1;
+    %% 加误差
+    PD0(Loadi)=PD0(Loadi).*(1+normrnd(0,0.05,length(Loadi),1));
+    QD0(Loadi)=QD0(Loadi).*(1+normrnd(0,0.05,length(Loadi),1));
+%     save('20PD0.mat','PD0');
+%     save('20QD0.mat','QD0');
+    load('20PD0.mat');
+    load('20QD0.mat');
+%     PG0(PGi)=PG0(PGi).*(1+normrnd(0,0.01,length(PGi),1));
+%     QG0(PVi)=QG0(PVi).*(1+normrnd(0,0.01,length(PVi),1));
+    %% 读变压器容量
+    noDataTransCapacity=0;
+    while(abs(Gap)>Precision)
+        if KK>kmax
            break;
        end
-
-        loopN=loopN+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);
+        %%
+        L_1Z=diag(Init_Z./Init_L);
+        U_1W=diag(Init_W./Init_U);
+        %% 形成海森阵
+        deltdeltF=func_deltdeltF(PVi,wPG,wQG,wVolt,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(Volt,PG,QG,PVi,PGi,wPG,wQG,wPD,wQD,wVolt,PG0,QG0,PD0,PD,QD,QD0,Volt0,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);
+        Mat_H=FormH(Busnum,Volt,PG,PD,QG,QD,Y,UAngel,r,c,Angle,QD_NON_ZERO,QD_NON_ZERO_IND,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,noDataTransCapacity);
+        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]=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);
+        Gap=(Init_L*Init_Z'-Init_U*Init_W');
+        KK=KK+1;
    end
+    fprintf('迭代次数%d\n',KK);
+    fprintf('目标值%f\n',full(ObjectiveFun(PG,PG0,PGi,QG,QG0,PVi,PD,PD0,QD,QD0,wPG,wQG,wPD,wQD,Loadi)));
+    % DrawGap(plotGap);
+    %%
+    %Volt=full(Volt');
+    %PD=full(PD);
+    %% 统计PD误差
+    % absPDLoad=abs(  (PD(Loadi)-PD0(Loadi))./PD0(Loadi)  );
+    absPDLoad=abs(  (PD(Loadi)-RealPD(Loadi))./RealPD(Loadi)  );
+    maxPDError=max(absPDLoad(absPDLoad<10))
+    Loadi(maxPDError==absPDLoad)
+    absQDLoad=abs(  (QD(Loadi)-RealQD(Loadi))./RealQD(Loadi)  );
+    maxQDError=max(absQDLoad(absQDLoad<10))
+    Loadi(maxQDError==absQDLoad)
+    disp('index');
+    %Loadi(absPDLoad==maxPDError);
+    %% 计算总线损
+    totalLoss=(sum(PG)-sum(PD(Loadi)))*100;
+    fprintf('总的损耗为%f(MW 有名值)\n',full(totalLoss));
+    fprintf('线损率为 %f\n',full(totalLoss/sum(PG)));
+    %% 计算各线损
+    %Lineloss(Linei,Linej,Liner,Linex,Lineb,Transfori,Transforj,Transforr,Transforx,Transfork0,Volt,UAngel);
+    thesis=thesis.AddPDQDPGQG(PD(Loadi),QD(Loadi),PG(Balance),QG(PVi));
    
 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));
+PD(Loadi)=thesis.MeanPD();
+QD(Loadi)=thesis.MeanQD();
+PG(Balance)=thesis.MeanPG();
+QG(PVi)=thesis.MeanQG();
+thesis.MaxDeviation(RealPG(Balance),RealQG(PVi),RealPD(Loadi),RealQD(Loadi))
+thesis.StatDeviation(RealPG(Balance),RealQG(PVi),RealPD(Loadi),RealQD(Loadi))
+thesis.MaxDeviation(PG0(Balance),QG0(PVi),PD0(Loadi),QD0(Loadi))
+thesis.StatDeviation(PG0(Balance),QG0(PVi),PD0(Loadi),QD0(Loadi))
+% thesis.SqareDeviation(RealPG(Balance(1)),RealQG(PVi(1)),RealPD(Loadi),RealQD(Loadi))
+thesis.SqareDeviation(PG0(Balance(1)),QG0(PVi(1)),PD0(Loadi),QD0(Loadi))
+thesis.PercentOfPass(PG0(Balance(1)),QG0(PVi(1)),PD0(Loadi),QD0(Loadi))
+thesis.MaxBranchDeviation(Linei,Linej,Liner,Linex,Lineb,Transfori,Transforj,Transforr,Transforx,Branchi,Branchg,Branchb,Transfork0,Volt0,UAngel0,Volt,UAngel,'E:\算例\柳金Ⅰ926_21671693_2012-09-06\newFIle20.txt',PD0,QD0)
+thesis.StatBranchDeviation(Linei,Linej,Liner,Linex,Lineb,Transfori,Transforj,Transforr,Transforx,Branchi,Branchg,Branchb,Transfork0,Volt0,UAngel0,Volt,UAngel)
+toc

-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();
--- a/OPFWrongDataCheck.asv
+++ b/OPFWrongDataCheck.asv
@@ -0,0 +1,106 @@
+clc
+    
+    clear
+ for II=1:16   
+     tic
+    [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]= ...
+        pf('D:\Project\最小化潮流\最小潮流算例\金湖924（2-1）_0.5_85%.txt');
+    %pf('c:/file31.txt');
+    %pf('ieee118PG.dat');
+    
+    %% 计算功率因数
+    atan(PD(QD~=0 | PD~=0)./QD(QD~=0 | PD~=0));
+    Volt;
+    UAngel*180/3.1415926;
+    %% 通过潮流计算PG
+    AngleIJ=sparse(r,c,UAngel(r)-UAngel(c)-Angle',Busnum,Busnum);
+    PGBal=diag(Volt)*Y.*cos(AngleIJ)*Volt';
+    
+    %% 初值-即测量值
+    PG0=PG;
+    PD0=PD;
+    
+    PDReal=PD;%真值
+    
+    %%
+    PG0(Balance)=PGBal(Balance);
+    %%
+    [Volt,UAngel,Init_Z,Init_W,Init_L,Init_U,Init_Y,PG,QG,RestraintCount,wG,wD,PD,PD0,randPDind,Loadi]=OPF_Init(Busnum,Balance,PG,QG,Volt,GenU,GenL,PVi,PGi,PVQU,PVQL,PD0,QD,PD);
+    PD0(Loadi(II))=PD0(Loadi(II))*(1+0.018);
+    Gap=(Init_L*Init_Z'-Init_U*Init_W');
+    KK=0;
+    plotGap=zeros(1,50);
+    ContrlCount=size(PVi,1)+size(PGi,1)+size(Loadi,1)+Busnum*2;
+    kmax=60;
+    %% 20120523 临时
+    QD_NON_ZERO=QD(PD==0 & QD~=0);
+    QD_NON_ZERO_IND=find(PD==0 & QD~=0);
+    %%
+    while(abs(Gap)>Precision)
+        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);
+        %%
+        L_1Z=diag(Init_Z./Init_L);
+        U_1W=diag(Init_W./Init_U);
+        %% 形成海森阵
+        deltdeltF=func_deltdeltF(PVi,wG,wD,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);
+        %%  开始构建deltF
+        deltF=func_deltF(PG,PVi,PGi,wG,wD,PG0,PD0,PD,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,Loadi);
+        Mat_H=FormH(Busnum,Volt,PG,PD,QG,QD,Y,UAngel,r,c,Angle,QD_NON_ZERO,QD_NON_ZERO_IND);
+        Ly=Mat_H;
+        Lz=FormLz(Mat_G,Init_L,GenL,Busnum,PVQL,PD0,Loadi,KK);
+        Lw=FormLw(Mat_G,Init_U,GenU,Busnum,PVQU,PD0,Loadi,KK);
+        Lx=FormLx(deltF,deltH,Init_Y,deltG,Init_Z,Init_W);
+        YY=FormYY(Lul,Lz,Ly,Luu,Lw,Lx);
+        %% 开始解方程
+        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]=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,Loadi);
+        Gap=(Init_L*Init_Z'-Init_U*Init_W');
+        KK=KK+1;
+    end
+    fprintf('迭代次数%d\n',KK);
+    ObjectiveFun(PG,PG0,PGi,PD,PD0,wG,wD,Loadi);
+    %DrawGap(plotGap);
+    %%
+    %Volt=full(Volt');
+    %PD=full(PD);
+    %% 统计PD误差
+    % absPDLoad=abs(  (PD(Loadi)-PD0(Loadi))./PD0(Loadi)  );
+    absPDLoad=abs(  (PD(Loadi)-PDReal(Loadi))./PDReal(Loadi)  );
+    maxPDError=max(absPDLoad);
+    %disp('index')
+    LoadiArray=Loadi(absPDLoad==maxPDError);
+    if length(LoadiArray)>1
+        disp('没找出')
+        Loadi(II)
+    end
+    if length(LoadiArray)==1
+        if LoadiArray~=Loadi(II)
+            disp('没找出')
+            Loadi(II)
+        end
+    end
+    toc;
+    
+end
+
--- a/OPF_Init.asv
+++ b/OPF_Init.asv
@@ -0,0 +1,34 @@
+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*1; %约束条件数,放开所有QD
+t_Bal_volt=Volt(Balance);
+Volt=sparse(1*ones(1,Busnum));
+Volt(Balance)=t_Bal_volt;
+UAngel=sparse(1,Busnum);
+Init_Z=sparse(ones(1,RestraintCount));
+Init_W=sparse(-1*ones(1,RestraintCount));
+Init_L=1*sparse(ones(1,RestraintCount));
+Init_U=1*sparse(ones(1,RestraintCount));
+Init_Y=sparse(1,2*Busnum);%与学姐一致
+tPU=sparse(GenU(:,2));% 发电机有功上限
+tQU=sparse(PVQU(:,1));% 无功上限
+tPL=sparse(GenL(:,2));% 发电机有功下限
+tQL=sparse(PVQL(:,1));% 无功下限
+wPG=0;
+wQG=0;
+%randInt=randperm(size(Loadi,1));
+%randPDind=randInt(1:10);
+randPDind=0;
+wPD=1/.05^2*ones(size(Loadi,1),1);
+% wPD(1:2:end)=0;
+wQD=1/.05^2*zeros(size(Loadi,1),1);
+% wQD(1:2:end)=0;
+%wD(randPDind)=0;%一些负荷不约束
+%wD(7)=0;
+% wD(11)=0;
+PD=0.5*PD0;
+%powerFacter=0.98;
+%QD=PD*sqrt((1-powerFacter^2)/powerFacter^2);
+QD=.5*QD0;
+end
--- a/OPF_Init.m
+++ b/OPF_Init.m
@@ -1,9 +1,7 @@
-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)
+function [Volt,UAngel,Init_Z,Init_W,Init_L,Init_U,Init_Y,PG,QG,RestraintCount,wPG,wQG,wPD,wQD,wVolt,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)*4+Busnum*2+Busnum+length(Loadi)*2; %约束条件数,放开所有QD
-%互补约束数量
-RestraintCount=RestraintCount+Busnum+length(Loadi)*2;
+RestraintCount=size(Loadi,1)*2+Busnum*1; %约束条件数,放开所有QD
 t_Bal_volt=Volt(Balance);
 Volt=sparse(1*ones(1,Busnum));
 Volt(Balance)=t_Bal_volt;
@@ -22,15 +20,16 @@ wQG=0;
 %randInt=randperm(size(Loadi,1));
 %randPDind=randInt(1:10);
 randPDind=0;
-wPD=zeros(size(Loadi,1),1);
-wPD(1:2:end)=0;
-wQD=zeros(size(Loadi,1),1);
-wQD(1:2:end)=0;
+wPD=1/.05^2*ones(size(Loadi,1),1);
+% wPD(1:2:end)=0;
+wQD=1/.05^2*ones(size(Loadi,1),1);
+% wQD(1:2:end)=0;
 %wD(randPDind)=0;%一些负荷不约束
 %wD(7)=0;
 % wD(11)=0;
-PD=0.0*PD0;
+PD=0.5*PD0;
 %powerFacter=0.98;
 %QD=PD*sqrt((1-powerFacter^2)/powerFacter^2);
-QD=.0*QD0;
+QD=.5*QD0;
+wVolt=1/.05^2*ones(size(Busnum,1),1);
 end
--- a/PD0.mat
+++ b/PD0.mat
--- a/PDBar.m
+++ b/PDBar.m
@@ -1,83 +0,0 @@
-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);
-
--- a/QD0.mat
+++ b/QD0.mat
--- a/QDBar.m
+++ b/QDBar.m
@@ -1,77 +0,0 @@
-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);
-
--- a/SolveIt.asv
+++ b/SolveIt.asv
@@ -0,0 +1,42 @@
+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)
+LxComa=FormLxComa(deltF,deltG,deltH,Init_L,Luu,Lul,Init_Z,Init_Y,Lz,Init_U,Init_W,Lw,Lx);
+H=-deltdeltF+ddh;%+ddg*(Init_Z'+Init_W');
+t1=diag(Init_L.\Init_Z-Init_U.\Init_W);
+t2=-deltG*( t1 )*deltG';
+aa=[
+    (H+t2),deltH;
+    deltH',zeros(size(Init_Y,2));
+    ];
+yy=[LxComa;-Ly];
+%% Æ½ºâ½ÚµãµçÑ¹²»±ä 
+t=size(PVi,1)+size(PGi,1)+size(Loadi,1);
+aa(t+Balance,:)=0;
+aa(:,t+Balance)=0;
+%aa(t+Balance,t+Balance)=1;
+aa=aa+sparse(t+Balance,t+Balance,ones(length(Balance)),ContrlCount+2*Busnum,ContrlCount+2*Busnum);
+deltG(t+Balance,:)=0;
+%%
+t=size(PVi,1)+size(PGi,1)+size(Loadi,1)+Busnum*1;
+aa(t+Balance,:)=0;
+aa(:,t+Balance)=0;
+%aa(t+Balance,t+Balance)=1;
+aa=aa+sparse(t+Balance,t+Balance,ones(length(Balance)),ContrlCount,ContrlCount);
+deltG(t+Balance,:)=0;
+%%
+dxdy=aa\yy;
+dX=dxdy(1:ContrlCount);
+dY=dxdy(ContrlCount+1:ContrlCount+2*Busnum);
+dL=Lz+deltG'*dX;
+dU=-Lw-deltG'*dX;
+dZ=-diag(Init_L)\Lul-diag(Init_L)\diag(Init_Z)*dL;
+dW=-diag(Init_U)\Luu-diag(Init_U)\diag(Init_W)*dU;
+XX=[
+    dX;
+    dY;
+    dZ;
+    dW;
+    dL;
+    dU;
+    
+];
+end
--- a/SolveIt.m
+++ b/SolveIt.m
@@ -1,6 +1,6 @@
-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)
+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)
 LxComa=FormLxComa(deltF,deltG,deltH,Init_L,Luu,Lul,Init_Z,Init_Y,Lz,Init_U,Init_W,Lw,Lx);
-H=-(deltdeltF-ddh-ddgzw);
+H=-deltdeltF+ddh;
 %t1=diag(Init_L.\Init_Z-Init_U.\Init_W);
 t1=diag(Init_Z./Init_L-Init_W./Init_U);
 t2=-deltG*( t1 )*deltG';
--- a/VoltBar.m
+++ b/VoltBar.m
@@ -1,79 +0,0 @@
-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);
-
--- a/admmatrix.asv
+++ b/admmatrix.asv
@@ -0,0 +1,31 @@
+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);
--- a/case1V.mat
+++ b/case1V.mat
--- a/case2V.mat
+++ b/case2V.mat
--- a/case3V.mat
+++ b/case3V.mat
--- a/caseM.mat
+++ b/caseM.mat
--- a/caseR.mat
+++ b/caseR.mat
--- a/createLoadProfileFigure.m
+++ b/createLoadProfileFigure.m
@@ -1,49 +0,0 @@
-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]);
-
--- a/figure1.fig
+++ b/figure1.fig
--- a/figure2.fig
+++ b/figure2.fig
--- a/func_ddg.asv
+++ b/func_ddg.asv
@@ -0,0 +1,9 @@
+function ddg=func_ddg(PGi,PVi,Busnum,RestraintCount,Loadi,PD,QD)
+
+ddg1=sparse(size(PVi,1)+size(PGi,1)+Loadi*2+Busnum,RestraintCount);
+t1=sparse(length(Loadi),size(PVi,1)+size(PGi,1)+Busnum)
+ddg2=sparse(length(Loadi),diag(-2./PD.^2-2*(QD.^2-PD.^2)/(PD.^2+QD.^2)));
+ddg3=diag(2*(PD.^2-QD.^2)/(PD.^2+QD.^2));
+ddg4=sparse(Busnum,RestraintCount);
+ddg=[ddg1;ddg2;ddg3;ddg4];
+end
--- a/func_ddg.m
+++ b/func_ddg.m
@@ -1,17 +1,4 @@
-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);
-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);
-];
+function ddg=func_ddg(PGi,PVi,Busnum,RestraintCount,Loadi,PD,QD)
+
+ddg=0;
 end
--- a/func_ddh.m
+++ b/func_ddh.m
@@ -51,7 +51,6 @@ t=[ddPdVdV+ddQdVdV,ddPdVdT+ddQdVdT;
 sizeLoadi=size(Loadi,1)*2;
 ddh=[
    sparse(sizeLoadi,ContrlCount);
-    sparse(2*Busnum,sizeLoadi),-t,sparse(2*Busnum,Busnum+length(Loadi)*2);
-    sparse(Busnum+length(Loadi)*2,ContrlCount);
+    sparse(2*Busnum,sizeLoadi),-t;
    ];
 end
--- a/func_deltF.asv
+++ b/func_deltF.asv
@@ -0,0 +1,19 @@
+function deltF=func_deltF(PG,PVi,GenC,ContrlCount,PGi,wG,wD,PG0,PD0,PD,Busnum)
+%t1=PG(setdiff(PVi,Balance));
+% t2=Volt'*Volt;
+% t3=real(GB).*cos(AngleIJMat)+imag(GB).*sin(AngleIJMat);
+% t4=t2.*t3;
+% t5=sum(t4,2);
+% PBal=t5(Balance);
+% PPG=([PQ(1),PBal])';%暂时用土办法处理一下
+%% 
+t1=2
+t2=2*wD.*(PD-PD0);
+deltF=[
+    zeros(size(PGi));
+    zeros(size(PVi));
+    t2;
+    zeros(2*Busnum,1);
+];
+
+end
--- a/func_deltF.m
+++ b/func_deltF.m
@@ -1,10 +1,11 @@
-function deltF=func_deltF(PG,QG,PVi,PGi,wPG,wQG,wPD,wQD,PG0,QG0,PD0,PD,QD,QD0,Busnum,Loadi)
-deltF=[sparse(length(Loadi),1);
-    sparse(length(Loadi),1);
-    sparse(2*Busnum,1);
-    sparse(ones(Busnum,1));
-    sparse(ones(length(Loadi),1));
-    sparse(ones(length(Loadi),1))
+function deltF=func_deltF(Volt,PG,QG,PVi,PGi,wPG,wQG,wPD,wQD,wVolt,PG0,QG0,PD0,PD,QD,QD0,Volt0,Busnum,Loadi)
+t3=2*wPD.*(PD(Loadi)-PD0(Loadi));
+t4=2*wQD.*(QD(Loadi)-QD0(Loadi));
+t5=2*wVolt.*(Volt-Volt0);
+deltF=[sparse(t3);
+    sparse(t4);
+    sparse(t5');
+    sparse(1*Busnum,1);
 ];

 end
--- a/func_deltG.asv
+++ b/func_deltG.asv
@@ -0,0 +1,46 @@
+function deltG=func_deltG(Busnum,PVi,PGi,Loadi,PD,QD)
+%%
+sizePGi=size(PGi,1);
+sizePVi=size(PVi,1);
+sizeLoadi=size(Loadi,1);
+%%
+dg1_dPg=sparse(1:sizePGi,1:sizePGi,ones(sizePGi,1),sizePGi,sizePGi);
+dg2_dPg=sparse(sizePGi,sizePVi);
+dg3_dPg=sparse(sizePGi,sizeLoadi);
+dg4_dPg=sparse(sizePGi,length(Loadi));
+dg5_dPg=sparse(sizePGi,Busnum);
+%% 
+dg1_dQr=sparse(sizePVi,sizePGi);
+dg2_dQr=sparse(1:sizePVi,1:sizePVi,ones(sizePVi,1),sizePVi,sizePVi);
+dg3_dQr=sparse(sizePVi,sizeLoadi);
+dg4_dQr=sparse(sizePVi,1);
+dg5_dQr=sparse(sizePVi,Busnum);
+%%
+dg1_dPD=sparse(size(Loadi,1),size(PGi,1));
+dg2_dPD=sparse(size(Loadi,1),size(PVi,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),1);
+dg5_dPD=sparse(size(Loadi,1),Busnum);
+%% 
+dg1_dQD=sparse(size(Loadi(1),1),size(PGi,1));
+dg2_dQD=sparse(size(Loadi(1),1),size(PVi,1));
+dg3_dQD=sparse(length(Loadi(1)),length(Loadi));
+dg4_dQD=sparse(1:size(Loadi(1),1),1:size(Loadi(1),1),ones(size(Loadi(1),1),1),size(Loadi(1),1),size(Loadi(1),1));
+dg5_dQD=sparse(size(Loadi(1),1),Busnum);
+
+%% 
+dg1_dx=sparse(2*Busnum,sizePGi);
+dg2_dx=sparse(2*Busnum,sizePVi);
+dg3_dx=sparse(2*Busnum,sizeLoadi);
+dg4_dx=sparse(2*Busnum,1);
+dg5_dx=[sparse(1:Busnum,1:Busnum,ones(Busnum,1),Busnum,Busnum);
+    sparse(Busnum,Busnum);
+    ];
+%% 
+deltG=[dg1_dPg,dg2_dPg,dg3_dPg,dg4_dPg,dg5_dPg;
+       dg1_dQr,dg2_dQr,dg3_dQr,dg4_dQr,dg5_dQr;
+       dg1_dPD,dg2_dPD,dg3_dPD,dg4_dPD,dg5_dPD;
+       dg1_dQD,dg2_dQD,dg3_dQD,dg4_dQD,dg5_dQD;
+       dg1_dx,dg2_dx,dg3_dx,dg4_dx,dg5_dx;
+];
+end
--- a/func_deltG.m
+++ b/func_deltG.m
@@ -1,107 +1,29 @@
-function deltG=func_deltG(Busnum,PVi,PGi,Loadi,PD,QD,Vbi,PDbi,QDbi)
+function deltG=func_deltG(Busnum,PVi,PGi,Loadi,PD,QD)
+%%
+
 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=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=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=[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));
-dg42_dvbi=sparse(Busnum,length(Loadi));
-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));
-%%
-dg3_dPDbi=sparse(-eye(sizeLoadi));
-dg32_dPDbi=sparse(eye(sizeLoadi));
-dg4_dPDbi=sparse(length(Loadi),length(Loadi));
-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;
+%% 
+
+deltG=[dg3_dPD,dg4_dPD,dg5_dPD;
+       dg3_dQD,dg4_dQD,dg5_dQD;
+       dg3_dx,dg4_dx,dg5_dx;
 ];
 end
--- a/func_deltH.asv
+++ b/func_deltH.asv
@@ -0,0 +1,8 @@
+function deltH=func_deltH(Busnum,Volt,PVi,Y,PGi,UAngel,r,c,Angle,Loadi)
+
+dH_dPg=sparse(1:size(PGi,1),PGi,ones(size(PGi,1),1),size(PGi,1),2*Busnum);
+dH_dQr=sparse(1:size(PVi,1),PVi+Busnum,ones(size(PVi,1),1),size(PVi,1),2*Busnum);
+dH_dPD=[sparse(1:size(Loadi,1),Loadi,-ones(size(Loadi,1),1),size(Loadi,1),Busnum) sparse(size(Loadi,1),Busnum)];
+dH_dx = jacobian_M3(Busnum,Volt,Y,Angle,UAngel,r,c);  %形成雅克比矩阵
+deltH=[dH_dPg;dH_dQr;dH_dPD;dH_dx'];
+end
--- a/func_deltH.m
+++ b/func_deltH.m
@@ -2,9 +2,6 @@ function deltH=func_deltH(Busnum,Volt,PVi,Y,PGi,UAngel,r,c,Angle,Loadi)
 dH_dPD=[sparse(1:size(Loadi,1),Loadi,-ones(size(Loadi,1),1),size(Loadi,1),Busnum) sparse(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);
-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];
+deltH=[dH_dPD;dH_dQD;dH_dx'];

 end
--- a/func_deltdeltF.asv
+++ b/func_deltdeltF.asv
@@ -0,0 +1,11 @@
+function deltdeltF=func_deltdeltF(Busnum,PVi,PGi,wG,wD,ContrlCount)
+
+ContrlCount=size(PVi,1)+size(PGi,1)+size(wD,1)+Busnum*2; %P,Q,Volt theta这些控制变量数
+C=[wG' zeros(size(PVi))' wD'];
+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
--- a/func_deltdeltF.m
+++ b/func_deltdeltF.m
@@ -1,3 +1,11 @@
-function deltdeltF=func_deltdeltF(PVi,wPG,wQG,wPD,wQD,ContrlCount)
-deltdeltF=0;
+function deltdeltF=func_deltdeltF(PVi,wPG,wQG,wVolt,wPD,wQD,ContrlCount)
+
+%ContrlCount=size(PVi,1)+size(PGi,1)+size(wD,1)+Busnum*2; %P,Q,Volt theta这些控制变量数
+C=[wPD' wQD' wVolt];
+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
--- a/graph/badDataGraph.m
+++ b/graph/badDataGraph.m
@@ -1,16 +0,0 @@
-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]);
--- a/histA.fig
+++ b/histA.fig
--- a/histP.fig
+++ b/histP.fig
--- a/histQ.fig
+++ b/histQ.fig
--- a/histV.fig
+++ b/histV.fig
--- a/hs_err_pid4476.log
+++ b/hs_err_pid4476.log
@@ -1,258 +0,0 @@
-#
-# 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
-
--- a/jacobian.asv
+++ b/jacobian.asv
@@ -0,0 +1,40 @@
+function [Jacob,PQ,U,Uangle]=jacobian(Busnum,Balance,PVi,PVu,U,Uangle,Y,Angle,P0,Q0,r,c)   
+%**************************************************************************
+%  程序功能 : 子函数——形成雅可比矩阵Jacobian
+%  编    者： 
+%  编制时间：2010.12
+%**************************************************************************
+%% 分别求雅克比矩阵的子阵H,L,N,J及有功无功分量P,Q
+AngleIJ = Uangle(r) - Uangle(c)- Angle';
+U(PVi) = PVu;
+U(Balance)=;
+temp1= -sparse(1:Busnum,1:Busnum,U,Busnum,Busnum)*Y*sparse(1:Busnum,1:Busnum,U,Busnum,Busnum);                                                          % 计算雅克比矩阵可利用的中间变量
+temp2 = sum(temp1.*sparse(r,c,sin(AngleIJ)),2);                     
+temp3 = sum(temp1.*sparse(r,c,cos(AngleIJ)),2);
+temp4=sparse(1:Busnum,1:Busnum,temp2,Busnum,Busnum);
+temp5=sparse(1:Busnum,1:Busnum,temp3,Busnum,Busnum);
+H = temp1.*sparse(r,c,sin(AngleIJ))-temp4;
+L = temp1.*sparse(r,c,sin(AngleIJ))+temp4;
+N = temp1.*sparse(r,c,cos(AngleIJ))+temp5;
+J = -temp1.*sparse(r,c,cos(AngleIJ))+temp5;
+
+Q = Q0+temp2';                                                                              %求有功分量P
+P = P0+temp3';                                                                               %求无功分量Q
+%% 处理平衡节点和pv节点
+H(:,Balance) = 0;                                                     
+H(Balance,:) = 0;
+%H(Balance,Balance) = 100;                                                           % 平衡节点对应的对角元素置一个有限数
+H=H+sparse(Balance,Balance,ones(1,length(Balance)),Busnum,Busnum);
+L(:,PVi) = 0;
+L(PVi,:) = 0;
+L = L+sparse(PVi,PVi,ones(1,length(PVi)),Busnum,Busnum); % PV节点对应的对角元素置为1
+J(:,Balance) = 0;
+J(PVi,:) = 0;
+N(:,PVi) = 0;
+N(Balance,:) = 0; 
+Q(PVi) = 0;                                                                                       % 将pv节点的无功不平衡分量置零   
+P(Balance) = 0;                                                                               % 平衡节点的有功功率不平衡分量置零
+%% 合成PQ和雅可比矩阵
+PQ = cat(2,P,Q);                                                                             % 形成功率不平衡分量列向量
+Jacob = cat(1,cat(2,H,N),cat(2,J,L));                                            % 形成Jacobian矩阵
+end
--- a/loadprofile.txt
+++ b/loadprofile.txt
@@ -1,96 +0,0 @@
-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
--- a/mVolt.mat
+++ b/mVolt.mat
--- a/maxBoundError.txt
+++ b/maxBoundError.txt
@@ -1,4 +0,0 @@
-0.010149012	0.000145532
-0.010613753	0.001269963
-0.010686824	0.001667111
-0.01098127	0.002552894
--- a/maxError.txt
+++ b/maxError.txt
@@ -1,33 +0,0 @@
-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
--- a/modifyadmmatrix.asv
+++ b/modifyadmmatrix.asv
@@ -0,0 +1,18 @@
+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));                                      %求节点导纳角度
--- a/openfile.asv
+++ b/openfile.asv
@@ -0,0 +1,64 @@
+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
--- a/Show More
+++ b/Show More