diff --git a/@ForThesis/AddPDQDPGQG.m b/@ForThesis/AddPDQDPGQG.m
index ead1720..f5a179c 100644
--- a/@ForThesis/AddPDQDPGQG.m
+++ b/@ForThesis/AddPDQDPGQG.m
@@ -9,6 +9,6 @@ this.currentPos=this.currentPos+1;
 this.PDArray(:,this.currentPos)=PD;
 this.QDArray(:,this.currentPos)=QD;
 this.PGArray(this.currentPos)=PG;
-% this.QGArray(this.currentPos)=QG;
+this.QGArray(this.currentPos)=QG;
 end
 
diff --git a/@ForThesis/MaxBranchDeviation.asv b/@ForThesis/MaxBranchDeviation.asv
new file mode 100644
index 0000000..aeaa522
--- /dev/null
+++ b/@ForThesis/MaxBranchDeviation.asv
@@ -0,0 +1,10 @@
+function [ output_args ] = MaxBranchDeviation( this, Linei,Linej,Liner,Linex,Lineb,Transfori,Transforj,Transforr,Transforx,Branchg,Branchb,Transfork0,Volt0,UAngel0,Volt,UAngel)
+%% 最大支路功率偏差
+% 支路功率包括线路和变压器
+%% 
+[dispLineloss0 dispTransloss0]=Lineloss(Linei,Linej,Liner,Linex,Lineb2,Transi,Transj,Transr,Transx,Branchg,Branchb,k0,Volt0,Angle0);
+[dispLineloss dispTransloss]=Lineloss(Linei,Linej,Liner,Linex,Lineb2,Transi,Transj,Transr,Transx,Branchg,Branchb,k0,Volt,Angle);
+real(dispLineloss0 - dispLineloss)./real(dispLineloss0)
+
+end
+
diff --git a/@ForThesis/MaxBranchDeviation.m b/@ForThesis/MaxBranchDeviation.m
new file mode 100644
index 0000000..dfdfca4
--- /dev/null
+++ b/@ForThesis/MaxBranchDeviation.m
@@ -0,0 +1,14 @@
+function [ output_args ] = MaxBranchDeviation( ~, Linei,Linej,Liner,Linex,Lineb,Transfori,Transforj,Transforr,Transforx,Branchi,Branchg,Branchb,Transfork0,Volt0,UAngel0,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;% 太小的值不计算
+t22=dispTransloss0(:,3)>1e-5;% 太小的值不计算
+t3=abs([t1(t11);t2(t22)]);
+output_args=max(t3(t3~=Inf));
+end
+
diff --git a/@ForThesis/MaxDeviation.m b/@ForThesis/MaxDeviation.m
index 3686bef..006d74e 100644
--- a/@ForThesis/MaxDeviation.m
+++ b/@ForThesis/MaxDeviation.m
@@ -6,6 +6,6 @@ QDMaxDev=max(abs(  (this.QDArray-QD0Array)./QD0Array ),[],2);
 PG0Array=repmat(PG0,this.sampleNum,1);
 QG0Array=repmat(QG0,this.sampleNum,1);
 PGMaxDev=max(abs( (PG0Array-this.PGArray)./PG0Array  ));
-% QGMaxDev=max( abs( (QG0Array-this.QGArray)./QG0Array  ) );
-output_arg=max([PDMaxDev;QDMaxDev;PGMaxDev]);
+QGMaxDev=max( abs( (QG0Array-this.QGArray)./QG0Array  ) );
+output_arg=max([PDMaxDev;QDMaxDev;PGMaxDev;QGMaxDev]);
 end
diff --git a/@ForThesis/PercentOfPass.m b/@ForThesis/PercentOfPass.m
new file mode 100644
index 0000000..86044bb
--- /dev/null
+++ b/@ForThesis/PercentOfPass.m
@@ -0,0 +1,12 @@
+function [ output_args ] = PercentOfPass( this, PG0,QG0,PD0,QD0 )
+%% 估计值合格率
+PD=this.PDArray(:,end);%用最后一列
+QD=this.QDArray(:,end);%用最后一列
+t1=PD-PD0;
+t2=QD-QD0;
+t3=abs(t1)/3*0.05;
+t4=abs(t2)/3*0.05;
+output_args=sum(t3<=1)+sum(t4<=1);
+
+end
+
diff --git a/@ForThesis/SqareDeviation.m b/@ForThesis/SqareDeviation.m
new file mode 100644
index 0000000..be5876f
--- /dev/null
+++ b/@ForThesis/SqareDeviation.m
@@ -0,0 +1,13 @@
+function [ output_args ] = SqareDeviation( this, PG0,QG0,PD0,QD0)
+%% 其实就是负荷系数不为0的目标函数
+
+PD=this.PDArray(:,end);%用最后一列
+QD=this.QDArray(:,end);%用最后一列
+t1=PD-PD0;
+t2=t1./0.05;
+t3=QD-QD0;
+t4=t3./0.05;
+
+output_args=sqrt(sum(t2.^2)+sum(t4.^2));
+end
+
diff --git a/@ForThesis/StatBranchDeviation.m b/@ForThesis/StatBranchDeviation.m
new file mode 100644
index 0000000..4acc382
--- /dev/null
+++ b/@ForThesis/StatBranchDeviation.m
@@ -0,0 +1,19 @@
+function [ output_args ] = StatBranchDeviation(  ~, Linei,Linej,Liner,Linex,Lineb,Transfori,Transforj,Transforr,Transforx,Branchi,Branchg,Branchb,Transfork0,Volt0,UAngel0,Volt,UAngel )
+%STATBRANCHDEVIATION Summary of this function goes here
+%   Detailed explanation goes here
+
+%% 支路功率统计偏差 L2 范数
+% 支路功率包括线路和变压器
+%% 
+[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=[t1(t11);t2(t22)];
+t4=t3.^2;
+t5=sum(t4).^.5;
+output_args=t5
+end
+
diff --git a/@ForThesis/StatDeviation.m b/@ForThesis/StatDeviation.m
index ff43cb2..5e5c442 100644
--- a/@ForThesis/StatDeviation.m
+++ b/@ForThesis/StatDeviation.m
@@ -6,15 +6,15 @@ PD0Array=repmat(PD0,1,this.sampleNum);
 QD0Array=repmat(QD0,1,this.sampleNum);
 PDDev=(this.PDArray-PD0Array)/0.05;
 QDDev=(this.QDArray-QD0Array)/0.05;
-PG0Array=repmat(PG0,this.sampleNum,1);
-QG0Array=repmat(QG0,this.sampleNum,1);
-PGDev=(PG0Array-this.PGArray)/0.01;
+% PG0Array=repmat(PG0,this.sampleNum,1);
+% QG0Array=repmat(QG0,this.sampleNum,1);
+% PGDev=(PG0Array-this.PGArray)/0.01;
 % QGDev=(QG0Array-this.QGArray)/0.01;
-wholeMat=[PDDev;QDDev;PGDev'];
+wholeMat=[PDDev;QDDev;];
 t1=wholeMat.^2;
 t2=sum(t1,1);
 t3=t2/size(t1,1);
-t4=t3.^2;
+t4=t3.^.5;
 output_args=sum(t4)/length(t4);
 end
 
diff --git a/FormLw.m b/FormLw.m
index 1ee1fc2..5088230 100644
--- a/FormLw.m
+++ b/FormLw.m
@@ -1,9 +1,7 @@
 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);
diff --git a/FormLz.m b/FormLz.m
index d527f03..48a736e 100644
--- a/FormLz.m
+++ b/FormLz.m
@@ -1,9 +1,6 @@
 function Lz=FormLz(Mat_G,Init_L,GenL,Busnum,PVQL,PD,PD0,QD0,Loadi,KK,PF)
 KK=999;
-%PL=GenL(:,2);%发电机有功下界
-PL=-5*ones(length(GenL(:,2)),1);
-%QL=PVQL(:,1);%发电机无功下界
-QL=-5*ones(length(PVQL(:,1)),1);
+
 VoltL=(0.9)*ones(1,Busnum);
 %VoltL=-10*ones(1,Busnum);
 PDL=PD0(Loadi);
diff --git a/Lineloss.m b/Lineloss.m
index 4512077..3af8837 100644
--- a/Lineloss.m
+++ b/Lineloss.m
@@ -1,8 +1,8 @@
 %% 计算线损
-function Lineloss(Linei,Linej,Liner,Linex,Lineb2,Transi,Transj,Transr,Transx,k0,Volt,Angle)
+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')
+% fprintf('功率为有名值\n');
+% fprintf('节点号\t节点号\t有功损耗 MW\t无功损耗 MVar')
 cmpVolt=Volt.*cos(Angle)+1i*Volt.*sin(Angle);
 cmpVolt=cmpVolt';
 y0=1i*Lineb2;
@@ -16,13 +16,17 @@ deltLineS=Sij+Sji;
 dispLineloss=[Linei Linej real(deltLineS)*100 imag(deltLineS)*100];
 %full(dispLineloss)
 dispLineloss=sortrows(dispLineloss,-3);
-full(dispLineloss)
+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)
+full(dispTransloss);
 end
\ No newline at end of file
diff --git a/OPF.asv b/OPF.asv
index 551113e..fe989df 100644
--- a/OPF.asv
+++ b/OPF.asv
@@ -1,120 +1,137 @@
 tic
 clc
 clear
-[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;
-PDReal=PD;%真值
-QDReal=QD;%真值
-%PD0(12)=PD0(12)+0.001;
-%%
-PG0(Balance)=PGBal(Balance);
-PG(Balance)=PGBal(Balance);
-QG0(Balance)=QGBal(Balance);
-QG0(PVi)=QGBal(PVi);
-QG(PVi)=QGBal(PVi);
-%%
-[Volt,UAngel,Init_Z,Init_W,Init_L,Init_U,Init_Y,PG,QG,RestraintCount,wPG,wQG,wPD,wQD,PD,PD0,QD,randPDind,Loadi]=OPF_Init(Busnum,Balance,PG,QG,Volt,GenU,GenL,PVi,PGi,PVQU,PVQL,PD0,QD0,QD,PD);
-Gap=(Init_L*Init_Z'-Init_U*Init_W');
-KK=0;
-plotGap=zeros(1,60);
-ContrlCount=size(PVi,1)+size(PGi,1)+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));
-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);
+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');
     
-    %% 形成方程矩阵
-    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,PD0,QD0,Loadi,KK,PF);
-    Lw=FormLw(Mat_G,Init_U,GenU,Busnum,PVQU,PD,PD0,QD0,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);
+    
+    %% 计算功率因数
+    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=KK+1;
+    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
-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)-PDReal(Loadi))./PDReal(Loadi)  );
-maxPDError=max(absPDLoad(absPDLoad<10))
-absQDLoad=abs(  (QD(Loadi)-QDReal(Loadi))./QDReal(Loadi)  );
-maxQDError=max(absQDLoad(absQDLoad<10))
-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);
-toc  
+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
 
diff --git a/OPF.m b/OPF.m
index 5b06934..76f63d0 100644
--- a/OPF.m
+++ b/OPF.m
@@ -1,10 +1,13 @@
 tic
 clc
 clear
+%% 存在问题
+% 变压器变比的位置没有考虑，由于现在用的变比都是1，所以没有影响。 20130123
+%% 
 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');
+    [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:\算例\东际911_2751267_2012-09-05\newFIle20.txt');
     %pf('D:\Project\青秀降损项目\最小化潮流\最小潮流算例\原始\津头站津视922（3-1）_0.5_120%.txt');
     %pf('D:\Project\最小化潮流\最小潮流算例\仙海919.txt');
     %pf('c:/file31.txt');
@@ -25,19 +28,22 @@ for II=1:4
     QG0=QG;
     PD0=PD;
     QD0=QD;
-    PDReal=PD;%真值
-    QDReal=QD;%真值
-    %PD0(12)=PD0(12)+0.001;
+    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;
     %%
-    PVi=zeros(0,0);
-    PGi=zeros(0,0);
-    [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);
+    %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);
@@ -88,8 +94,8 @@ for II=1:4
         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,PD0,QD0,Loadi,KK,PF);
-        Lw=FormLw(Mat_G,Init_U,GenU,Busnum,PVQU,PD,PD0,QD0,Loadi,KK,PF,noDataTransCapacity);
+        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);
         %% 开始解方程
@@ -109,10 +115,12 @@ for II=1:4
     %PD=full(PD);
     %% 统计PD误差
     % absPDLoad=abs(  (PD(Loadi)-PD0(Loadi))./PD0(Loadi)  );
-    absPDLoad=abs(  (PD(Loadi)-PDReal(Loadi))./PDReal(Loadi)  );
+    absPDLoad=abs(  (PD(Loadi)-RealPD(Loadi))./RealPD(Loadi)  );
     maxPDError=max(absPDLoad(absPDLoad<10))
-    absQDLoad=abs(  (QD(Loadi)-QDReal(Loadi))./QDReal(Loadi)  );
+    Loadi(maxPDError==absPDLoad)
+    absQDLoad=abs(  (QD(Loadi)-RealQD(Loadi))./RealQD(Loadi)  );
     maxQDError=max(absQDLoad(absQDLoad<10))
+    Loadi(maxQDError==absQDLoad)
     disp('index');
     %Loadi(absPDLoad==maxPDError);
     %% 计算总线损
@@ -128,7 +136,12 @@ PD(Loadi)=thesis.MeanPD();
 QD(Loadi)=thesis.MeanQD();
 PG(Balance)=thesis.MeanPG();
 QG(PVi)=thesis.MeanQG();
-thesis.MaxDeviation(PG0(Balance),QG0(PVi),PD0(Loadi),QD0(Loadi));
-thesis.StatDeviation(PG0(Balance),QG0(PVi),PD0(Loadi),QD0(Loadi))
+thesis.MaxDeviation(RealPG(Balance),RealQG(PVi),RealPD(Loadi),RealQD(Loadi))
+thesis.StatDeviation(RealPG(Balance),RealQG(PVi),RealPD(Loadi),RealQD(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)
+thesis.StatBranchDeviation(Linei,Linej,Liner,Linex,Lineb,Transfori,Transforj,Transforr,Transforx,Branchi,Branchg,Branchb,Transfork0,Volt0,UAngel0,Volt,UAngel)
 toc
 
diff --git a/OPF_Init.m b/OPF_Init.m
index 1bc6caf..b9fdfc6 100644
--- a/OPF_Init.m
+++ b/OPF_Init.m
@@ -15,19 +15,20 @@ tPU=sparse(GenU(:,2));% 
 tQU=sparse(PVQU(:,1));% 无功上限
 tPL=sparse(GenL(:,2));% 发电机有功下限
 tQL=sparse(PVQL(:,1));% 无功下限
-PG(PGi)=(tPU+tPL)/2;
-wPG=0*ones(size(PGi,1),1);
-wQG=0*ones(size(PVi,1),0);
+wPG=0;
+wQG=0;
 %randInt=randperm(size(Loadi,1));
 %randPDind=randInt(1:10);
 randPDind=0;
-wPD=20*ones(size(Loadi,1),1);
-wQD=0;
+wPD=0/.05^2*ones(size(Loadi,1),1);
+% wPD([2 7])=0;
+wQD=0/.05^2*ones(size(Loadi,1),1);
+% wQD([2 7])=0;
 %wD(randPDind)=0;%一些负荷不约束
 %wD(7)=0;
 % wD(11)=0;
-PD=1*PD0;
+PD=0.5*PD0;
 %powerFacter=0.98;
 %QD=PD*sqrt((1-powerFacter^2)/powerFacter^2);
-QD=QD0;
+QD=.5*QD0;
 end
\ No newline at end of file
diff --git a/func_deltG.m b/func_deltG.m
index 893ff02..230db27 100644
--- a/func_deltG.m
+++ b/func_deltG.m
@@ -1,36 +1,20 @@
 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,length(Loadi));
-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),length(Loadi));
 dg5_dPD=sparse(size(Loadi,1),Busnum);
 %% 
-dg1_dQD=sparse(size(Loadi,1),size(PGi,1));
-dg2_dQD=sparse(size(Loadi,1),size(PVi,1));
 dg3_dQD=sparse(length(Loadi),length(Loadi));
 dg4_dQD=sparse(1:size(Loadi,1),1:size(Loadi,1),ones(size(Loadi,1),1),size(Loadi,1),size(Loadi,1));
 dg5_dQD=sparse(size(Loadi,1),Busnum);
 
 %% 
-dg1_dx=sparse(2*Busnum,sizePGi);
-dg2_dx=sparse(2*Busnum,sizePVi);
 dg3_dx=sparse(2*Busnum,sizeLoadi);
 dg4_dx=sparse(2*Busnum,length(Loadi));
 dg5_dx=[sparse(1:Busnum,1:Busnum,ones(Busnum,1),Busnum,Busnum);
diff --git a/func_deltH.m b/func_deltH.m
index 366f9ef..e2b8401 100644
--- a/func_deltH.m
+++ b/func_deltH.m
@@ -1,6 +1,4 @@
 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_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);  %形成雅克比矩阵
diff --git a/pf.m b/pf.m
index a1ed312..7527375 100644
--- a/pf.m
+++ b/pf.m
@@ -1,4 +1,4 @@
-function [kmax,Precision,Uangle,U,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(FileName)
+function [kmax,Precision,Uangle,U,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(FileName)
 %**************************************************************************
 %    程序名称：电力系统潮流计算程序
 %    程序算法：极坐标下的牛顿-拉夫逊法