diff --git a/.gitignore b/.gitignore
index 8f1013f..d1b71e6 100644
--- a/.gitignore
+++ b/.gitignore
@@ -1,4 +1,5 @@
 *.asv
 /feeder13
 *.mexw32
-/DistributionNetwork-OnlyFortescue
\ No newline at end of file
+#/DistributionNetwork-OnlyFortescue
+*.7z
\ No newline at end of file
diff --git a/AssignXX.m b/AssignXX.m
new file mode 100644
index 0000000..2ff0f36
--- /dev/null
+++ b/AssignXX.m
@@ -0,0 +1,23 @@
+function [deltZ,deltL,deltW,deltU,deltX,deltY]=AssignXX1(XX,ContrlCount,RestraintCount,Busnum)
+% deltX=XX(1:14);
+% deltY=XX(15:24);
+% deltZ=XX(25:38);
+% deltW=XX(39:52);
+% deltL=XX(53:66);
+% deltU=XX(67:80);
+deltX=XX(1:ContrlCount);
+k1=ContrlCount+2*Busnum;
+deltY=XX(ContrlCount+1:k1);
+k2=k1;
+k1=k2+RestraintCount;
+deltZ=XX(k2+1:k1);
+k2=k1;
+k1=k2+RestraintCount;
+deltW=XX(k2+1:k1);
+k2=k1;
+k1=k2+RestraintCount;
+deltL=XX(k2+1:k1);
+k2=k1;
+k1=k2+RestraintCount;
+deltU=XX(k2+1:k1);
+end
\ No newline at end of file
diff --git a/DistributionNetwork-OnlyFortescue b/DistributionNetwork-OnlyFortescue
new file mode 160000
index 0000000..ecaed83
--- /dev/null
+++ b/DistributionNetwork-OnlyFortescue
@@ -0,0 +1 @@
+Subproject commit ecaed830b547d025369bd6137707f6a96856556a
diff --git a/FormAlphaD.m b/FormAlphaD.m
new file mode 100644
index 0000000..bd264d2
--- /dev/null
+++ b/FormAlphaD.m
@@ -0,0 +1,9 @@
+function AlphaD=FormAlphaD(Init_Z,deltZ,Init_W,deltW)
+tdeltZinx=find(deltZ<0);
+tdeltWinx=find(deltW>0);
+t1=-Init_Z(tdeltZinx)./deltZ(tdeltZinx);
+t2=-Init_W(tdeltWinx)./deltW(tdeltWinx);
+t3=[t1,t2];
+t4=min(t3);
+AlphaD=0.9995*min([t4 1]);
+end
\ No newline at end of file
diff --git a/FormAlphaP.m b/FormAlphaP.m
new file mode 100644
index 0000000..a3cf2be
--- /dev/null
+++ b/FormAlphaP.m
@@ -0,0 +1,9 @@
+function AlphaP=FormAlphaP(Init_L,deltL,Init_U,deltU)
+tdeltLinx=find(deltL<0);
+tdeltUinx=find(deltU<0);
+t1=-Init_L(tdeltLinx)./deltL(tdeltLinx);
+t2=-Init_U(tdeltUinx)./deltU(tdeltUinx);
+t3=[t1,t2];
+t4=min(t3);
+AlphaP=0.9995*min([t4 1]);
+end
\ No newline at end of file
diff --git a/FormG.m b/FormG.m
new file mode 100644
index 0000000..fb8bbbf
--- /dev/null
+++ b/FormG.m
@@ -0,0 +1,7 @@
+function Mat_G=FormG(Volt,PD,QD)
+Mat_G=[
+    sparse(PD);
+    sparse(QD);
+    Volt;
+    ];
+end
\ No newline at end of file
diff --git a/FormH.m b/FormH.m
new file mode 100644
index 0000000..7ce3fba
--- /dev/null
+++ b/FormH.m
@@ -0,0 +1,14 @@
+function Mat_H=FormH(Busnum,Volt,PG,PD,QG,QD,Y,UAngel,r,c,Angle)
+%% 
+%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
\ No newline at end of file
diff --git a/FormLw.m b/FormLw.m
new file mode 100644
index 0000000..4685222
--- /dev/null
+++ b/FormLw.m
@@ -0,0 +1,10 @@
+function Lw=FormLw(Mat_G,Init_U,Busnum,Loadi)
+VoltU=sparse(1.1*ones(Busnum,1));
+PDU=sparse(1.5*ones(length(Loadi),1));
+QDU=sparse(1.5*ones(length(Loadi),1));
+t1=[PDU;
+    QDU;
+    VoltU];
+t2=Mat_G+Init_U-t1;
+Lw=t2;
+end
\ No newline at end of file
diff --git a/FormLx.m b/FormLx.m
new file mode 100644
index 0000000..826652b
--- /dev/null
+++ b/FormLx.m
@@ -0,0 +1,4 @@
+function Lx=FormLx(deltF,deltH,Init_Y,deltG,Init_Z,Init_W)
+t1=deltF-deltH*Init_Y-deltG*(Init_Z+Init_W);
+Lx=t1;
+end
\ No newline at end of file
diff --git a/FormLxComa.m b/FormLxComa.m
new file mode 100644
index 0000000..9d25a12
--- /dev/null
+++ b/FormLxComa.m
@@ -0,0 +1,11 @@
+function LxComa=FormLxComa(deltF,deltG,deltH,Init_L,Luu,Lul,Init_Z,Init_Y,Lz,Init_U,Init_W,Lw,Lx)
+%t1=deltG*(Init_Z'+Init_W');%%
+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;%
+t8=deltG*(t4+t7);%%
+LxComa=Lx+t8;
+end
\ No newline at end of file
diff --git a/FormLz.m b/FormLz.m
new file mode 100644
index 0000000..5c373fb
--- /dev/null
+++ b/FormLz.m
@@ -0,0 +1,12 @@
+function Lz=FormLz(Mat_G,Init_L,Busnum,Loadi)
+VoltL=sparse(0.9*ones(Busnum,1));
+PDL=sparse(0*ones(length(Loadi),1));
+QDL=sparse(0*ones(length(Loadi),1));
+t1=[PDL;
+    QDL;
+    VoltL
+    ];
+t2=Mat_G-Init_L-t1;
+Lz=t2;
+
+end
\ No newline at end of file
diff --git a/FormYY.m b/FormYY.m
new file mode 100644
index 0000000..cbad106
--- /dev/null
+++ b/FormYY.m
@@ -0,0 +1,10 @@
+function YY=FormYY(Lul,Lz,Ly,Luu,Lw,Lx)
+YY=[
+  Lx;
+  -Ly;
+  -Lz;
+  -Lw;
+  -Lul;
+  -Luu;  
+];
+end
\ No newline at end of file
diff --git a/Modification.m b/Modification.m
new file mode 100644
index 0000000..d55962c
--- /dev/null
+++ b/Modification.m
@@ -0,0 +1,18 @@
+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,ContrlCount,Balance,Busnum,PD,QD,Loadi)
+AlphaP=FormAlphaP(Init_L,deltL,Init_U,deltU);
+AlphaD=FormAlphaD(Init_Z,deltZ,Init_W,deltW);
+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;
+t=deltX(1:size(Loadi,1)*2);
+PD(Loadi)=PD(Loadi)+AlphaP*t(1:length(Loadi));
+QD(Loadi)=QD(Loadi)+AlphaP*t(length(Loadi)+1:length(Loadi)*2);
+t=deltX(size(Loadi,1)*2+1:ContrlCount);
+t(Busnum+Balance)=0;
+balVolt=Volt(Balance);
+Volt=Volt+AlphaP*t(1:Busnum);
+Volt(Balance)=balVolt;
+UAngel=UAngel+AlphaP*t(Busnum+1:2*Busnum);
+end
\ No newline at end of file
diff --git a/SolveIt.m b/SolveIt.m
new file mode 100644
index 0000000..039d223
--- /dev/null
+++ b/SolveIt.m
@@ -0,0 +1,46 @@
+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,Lx,Balance,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;
+t1=diag(Init_Z./Init_L-Init_W./Init_U);
+t2=-deltG*( t1 )*deltG';
+aa=[
+    (H+t2),deltH;
+    deltH',zeros(length(Init_Y));
+    ];
+yy=[LxComa;-Ly];
+%% 平衡节点电压不变 
+t=size(Loadi,1)*2;
+aa(t+Balance,:)=0;
+aa(:,t+Balance)=0;
+%aa(t+Balance,t+Balance)=1;
+aa=aa+sparse(t+Balance,t+Balance,ones(length(Balance),1),ContrlCount+2*Busnum,ContrlCount+2*Busnum);
+deltG(t+Balance,:)=0;
+%%
+t=size(Loadi,1)*2+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),1),ContrlCount+2*Busnum,ContrlCount+2*Busnum);
+deltG(t+Balance,:)=0;
+%%
+dxdy=aa\yy;
+%% KLU
+%spy(aa)
+%dxdy = klu(aa,'\',full(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
\ No newline at end of file
diff --git a/func_ddg.m b/func_ddg.m
new file mode 100644
index 0000000..7406e6d
--- /dev/null
+++ b/func_ddg.m
@@ -0,0 +1,4 @@
+function ddg=func_ddg()
+
+ddg=0;
+end
\ No newline at end of file
diff --git a/func_ddh.m b/func_ddh.m
new file mode 100644
index 0000000..92f4f86
--- /dev/null
+++ b/func_ddh.m
@@ -0,0 +1,48 @@
+function ddh=func_ddh(Volt,Init_Y,Busnum,Y,UAngel,r,c,Angle,Loadi,ContrlCount)
+mat_AngleIJ=sparse(r,c,UAngel(r)-UAngel(c)-Angle,Busnum,Busnum);
+mat_INV_AngleIJ=mat_AngleIJ';
+yP=Init_Y(1:length(Init_Y)/2);%暂时改这里  20111227
+yQ=Init_Y(length(Init_Y)/2+1:length(Init_Y));%暂时改这里  20111227
+t1=-sparse(1:Busnum,1:Busnum,Y.*cos(mat_INV_AngleIJ)*sparse(1:Busnum,1:Busnum,Volt,Busnum,Busnum)*yP,Busnum,Busnum);
+t2=sparse(1:Busnum,1:Busnum,sparse(1:Busnum,1:Busnum,Volt,Busnum,Busnum)*yP,Busnum,Busnum)*Y.*cos(mat_AngleIJ);
+t3=(t1+t2)*sparse(1:Busnum,1:Busnum,Volt,Busnum,Busnum);
+t4=-(sparse(1:Busnum,1:Busnum,Y.*cos(mat_AngleIJ)*Volt,Busnum,Busnum) -sparse(1:Busnum,1:Busnum,Volt,Busnum,Busnum)*Y.*cos(mat_INV_AngleIJ) )*sparse(1:Busnum,1:Busnum,sparse(1:Busnum,1:Busnum,Volt,Busnum,Busnum)*yP,Busnum,Busnum);
+ddPdTdT=t3+t4;%ok1
+t1=(-sparse(1:Busnum,1:Busnum,Y.*sin(mat_AngleIJ)*Volt,Busnum,Busnum)+sparse(1:Busnum,1:Busnum,Volt,Busnum,Busnum)*Y.*sin(mat_INV_AngleIJ) )*sparse(1:Busnum,1:Busnum,yP,Busnum,Busnum);
+t2= -sparse(1:Busnum,1:Busnum, sparse(1:Busnum,1:Busnum,Volt,Busnum,Busnum)*yP ,Busnum,Busnum)*Y.*sin(mat_AngleIJ)+sparse(1:Busnum,1:Busnum,Y.*sin(mat_INV_AngleIJ)*sparse(1:Busnum,1:Busnum,Volt,Busnum,Busnum)*yP,Busnum,Busnum);
+ddPdVdT=t1+t2;%ok1
+t1=diag( Y.*sin(mat_INV_AngleIJ)*diag(Volt)*yP);
+t2=diag(yP)*Y.*sin(mat_AngleIJ)*diag(Volt);
+t3=-diag(yP)*diag(Y.*sin(mat_AngleIJ)*Volt);
+t4=-Y.*sin(mat_INV_AngleIJ)*diag( diag(Volt)*yP );
+ddPdTdV=t1+t2+t3+t4;%存疑与我的不一样
+t1=Y.*cos(mat_INV_AngleIJ)*diag(yP);
+t2=diag(yP)*Y.*cos(mat_AngleIJ);
+ddPdVdV=t1+t2;
+t1=-diag(Y.*sin(mat_AngleIJ)*Volt);
+t2=diag(Volt)*Y.*sin(mat_INV_AngleIJ);
+t3=(t1+t2)*diag( diag(Volt)*yQ );
+t4=-diag( diag(Volt)*yQ )*Y.*sin(mat_AngleIJ);
+t5=diag(Y.*sin(mat_INV_AngleIJ)*diag(Volt)*yQ);
+t6=-(t4+t5)*diag(Volt);
+ddQdTdT=t3+t6;%ok1
+t1=(diag(Y.*cos(mat_AngleIJ)*Volt)-diag(Volt)*Y.*cos(mat_INV_AngleIJ) )*diag(yQ);
+t2=+diag( diag(Volt)*yQ )*Y.*cos(mat_AngleIJ)-diag(Y.*cos(mat_INV_AngleIJ)*diag(Volt)*yQ);
+ddQdVdT=t1+t2;
+t1=Y.*cos(mat_INV_AngleIJ)*diag(diag(Volt)*yQ);
+t2=diag(yQ)*diag(Y.*cos(mat_AngleIJ)*Volt);
+t3=-diag(Y.*cos(mat_INV_AngleIJ)*diag(Volt)*yQ);
+t4=-diag(yQ)*Y.*cos(mat_AngleIJ)*diag(Volt);
+ddQdTdV=t1+t2+t3+t4;
+t1=Y.*sin(mat_INV_AngleIJ)*diag(yQ);
+t2=diag(yQ)*Y.*sin(mat_AngleIJ);
+ddQdVdV=t1+t2;
+t=[ddPdVdV+ddQdVdV,ddPdVdT+ddQdVdT;
+   ddPdTdV+ddQdTdV,ddPdTdT+ddQdTdT;
+];
+sizeLoadi=size(Loadi,1)*2;
+ddh=[
+    sparse(sizeLoadi,ContrlCount);
+    sparse(2*Busnum,sizeLoadi),-t;
+    ];
+end
\ No newline at end of file
diff --git a/func_deltF.m b/func_deltF.m
new file mode 100644
index 0000000..7bfb3c6
--- /dev/null
+++ b/func_deltF.m
@@ -0,0 +1,8 @@
+function deltF=func_deltF(SEPD,ContrlCount)
+t=sparse(ones(length(SEPD),1));
+deltF=[
+    t;
+    sparse(ContrlCount-length(t),1);
+];
+
+end
\ No newline at end of file
diff --git a/func_deltG.m b/func_deltG.m
new file mode 100644
index 0000000..0522de4
--- /dev/null
+++ b/func_deltG.m
@@ -0,0 +1,22 @@
+function deltG=func_deltG(busNum,Loadi)
+dgV_dPD=sparse(length(Loadi),busNum);
+dgV_dQD=sparse(length(Loadi),busNum);
+dgV_dx=[sparse(1:busNum,1:busNum,ones(busNum,1),busNum,busNum);
+    sparse(busNum,busNum);
+    ];
+%
+dgPD_dPD=sparse(eye(length(Loadi)));
+dgPD_dQD=sparse(length(Loadi),length(Loadi));
+dgPD_dx=sparse(busNum*2,length(Loadi));
+%
+dgQD_dPD=sparse(length(Loadi),length(Loadi));
+dgQD_dQD=sparse(eye(length(Loadi)));
+dgQD_dx=sparse(busNum*2,length(Loadi));
+%% 
+
+deltG=[
+    dgPD_dPD,dgQD_dPD,dgV_dPD;
+    dgPD_dQD,dgQD_dQD,dgV_dQD;
+    dgPD_dx,dgQD_dx,dgV_dx;
+];
+end
\ No newline at end of file
diff --git a/func_deltH.m b/func_deltH.m
new file mode 100644
index 0000000..fa65c0a
--- /dev/null
+++ b/func_deltH.m
@@ -0,0 +1,7 @@
+function deltH=func_deltH(Busnum,Volt,Y,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(Busnum,Volt,Y,Angle,UAngel,r,c);  %形成雅克比矩阵
+deltH=[dH_dPD;dH_dQD;dH_dx'];
+
+end
\ No newline at end of file
diff --git a/func_deltdeltF.m b/func_deltdeltF.m
new file mode 100644
index 0000000..17fdbb2
--- /dev/null
+++ b/func_deltdeltF.m
@@ -0,0 +1,5 @@
+function deltdeltF=func_deltdeltF(SEPD,ContrlCount)
+deltdeltF=[
+    0*zeros(ContrlCount);
+    ];
+end
\ No newline at end of file
diff --git a/jacobian.m b/jacobian.m
new file mode 100644
index 0000000..45a7307
--- /dev/null
+++ b/jacobian.m
@@ -0,0 +1,20 @@
+function [Jacob]=jacobian(Busnum,Volt,Y,Angle,UAngel,r,c)   
+%**************************************************************************
+%  程序功能 : 子函数——形成雅可比矩阵Jacobian
+%  编    者： 
+%  编制时间：2010.12
+%**************************************************************************
+%%参照图书馆6楼的书编写
+%% 分别求雅克比矩阵的子阵H,L,N,J及有功无功分量P,Q
+AngleIJ=UAngel(r)-UAngel(c)-Angle;
+mat_AngleIJ=sparse(r,c,AngleIJ,Busnum,Busnum);
+mat_IvAngleIJ=mat_AngleIJ';
+H=sparse(1:Busnum,1:Busnum,Volt,Busnum,Busnum)*Y.*sin(mat_IvAngleIJ)*sparse(1:Busnum,1:Busnum,Volt,Busnum,Busnum)-sparse(1:Busnum,1:Busnum,Y.*sin(mat_AngleIJ)*Volt,Busnum,Busnum)*sparse(1:Busnum,1:Busnum,Volt,Busnum,Busnum);
+N=-sparse(1:Busnum,1:Busnum,Volt,Busnum,Busnum)*Y.*cos(mat_IvAngleIJ)*sparse(1:Busnum,1:Busnum,Volt,Busnum,Busnum)+sparse(1:Busnum,1:Busnum,Y.*cos(mat_AngleIJ)*Volt,Busnum,Busnum)*sparse(1:Busnum,1:Busnum,Volt,Busnum,Busnum);
+J=sparse(1:Busnum,1:Busnum,Y.*cos(mat_AngleIJ)*Volt,Busnum,Busnum)+Y.*cos(mat_IvAngleIJ)*sparse(1:Busnum,1:Busnum,Volt,Busnum,Busnum);
+L=sparse(1:Busnum,1:Busnum,Y.*sin(mat_AngleIJ)*Volt,Busnum,Busnum)+Y.*sin(mat_IvAngleIJ)*sparse(1:Busnum,1:Busnum,Volt,Busnum,Busnum);
+t1=[J,L;
+    H,N;
+]';
+Jacob=-t1;
+end
\ No newline at end of file
diff --git a/run.m b/run.m
index 40a0f0a..4c671ea 100644
--- a/run.m
+++ b/run.m
@@ -1,9 +1,9 @@
 clc
 clear
-lineZ=readLineZ('D:\share\feeder123\lineParameter.txt');
+lineZ=readLineZ('feeder13\lineParameter.txt');
 [ fsY0, fsY1, fsY2,phaseASpotLoadP,phaseBSpotLoadP,phaseCSpotLoadP ...
     phaseASpotLoadQ,phaseBSpotLoadQ,phaseCSpotLoadQ,setIJ,nodeNum,Balance,phaseABCY ...
-    cap]=dataRead(lineZ,'D:\share\feeder123\data.txt');
+    cap]=dataRead(lineZ,'feeder13\data1.txt');
 a=exp(1j*2*pi/3);
 Tp2f=1/3*[1 1 1;
     1 a a^2;
@@ -26,6 +26,7 @@ QG=sparse(busNum,1);
 QGi=[Balance];
 PD=Pabc/3;
 QD=Qabc/3;
+Loadi=find(PD~=0);
 maxD=100000;% 最大不平衡量
 EPS=1e-5;
 k=0;
@@ -35,7 +36,6 @@ fsY00=fsY0;
 fsY22=fsY2;
 Vf2=sparse(busNum,1);
 If2=sparse(busNum,1);
-% Vf0=sparse(busNum,1);
 Vf0=sparse(busNum,1);
 If0=sparse(busNum,1);
 %准备序矩阵
@@ -115,58 +115,110 @@ ub=checkSSatisfied(Balance,phaseABCY,VoltpABC, ...
     phaseASpotLoadQ,phaseBSpotLoadQ,phaseCSpotLoadQ );
 fprintf('最大不平衡量为%f\n\n',full(max(abs(ub))))
 fprintf('开始牛顿法迭代\n');
-%%用牛顿法求解
-[r,c,GB]=find(phaseABCY);
-Y=abs(phaseABCY);
-Yangle=angle(GB);
-Vp3=sparse(ones(busNum*3,1));%给电压赋初值
-Vp3(2:3:end)=Vp3(2:3:end)*exp(1j*-120/180*pi);
-Vp3(3:3:end)=Vp3(3:3:end)*exp(1j*+120/180*pi);
-PQi3P=zeros(length(PQi)*3,1);
-PQi3P(1:3:end)=(PQi-1)*3+1;
-PQi3P(2:3:end)=(PQi-1)*3+2;
-PQi3P(3:3:end)=(PQi-1)*3+3;
-PG=0;
-QG=0;
-PD3P=sparse(busNum*3,1);
-QD3P=sparse(busNum*3,1);
-PD3P(1:3:end)=phaseASpotLoadP;
-PD3P(2:3:end)=phaseBSpotLoadP;
-PD3P(3:3:end)=phaseCSpotLoadP;
-QD3P(1:3:end)=phaseASpotLoadQ;
-QD3P(2:3:end)=phaseBSpotLoadQ;
-QD3P(3:3:end)=phaseCSpotLoadQ;
-QGi3P=zeros(length(QGi)*3,1);
-QGi3P(1:3:end)=(QGi-1)*3+1;
-QGi3P(2:3:end)=(QGi-1)*3+2;
-QGi3P(3:3:end)=(QGi-1)*3+3;
-Vp3m=abs(Vp3);
-Vp3a=angle(Vp3);
-Balance3P=zeros(length(Balance)*3,1);
-Balance3P(1:3:end)=(Balance-1)*3+1;
-Balance3P(2:3:end)=(Balance-1)*3+2;
-Balance3P(3:3:end)=(Balance-1)*3+3;
-Vp3a((Balance-1)*3+1)=0;
-Vp3a((Balance-1)*3+2)=-120/180*pi;
-Vp3a((Balance-1)*3+3)=+120/180*pi;
-k=0;
-maxD=10000;
-tic
-while(k<=kmax && maxD> EPS)
-    k=k+1;
-    [dP, dQ, YdotSinVolt, YdotCosVolt, diag_Volt_YdotSin, diag_Volt_YdotCos]=Unbalance(Balance3P,busNum*3, ...
-        PQi3P,PG,QG,QGi3P,PD3P,QD3P,Vp3m,Vp3a,Y,Yangle,r,c,0,0,0,0);
-    maxD=max(abs([dP;dQ;]));
-    jaco=Jacobi(Balance3P,busNum*3,QGi3P,Vp3m,YdotSinVolt,YdotCosVolt,diag_Volt_YdotSin,diag_Volt_YdotCos);%雅克比矩阵
-    [dV, dVangle]=Solv(busNum*3,jaco,dP,dQ);%解出修正量
-    [Vp3m, Vp3a]=Modify(Vp3m,Vp3a,dV,dVangle,1);
-    fprintf('第 %d 次迭代， 最大不平衡量为 %f\n',k,full(maxD));
-    fprintf('迭代时间%f\n',toc);
-end
-NewtonToc=toc;
-fprintf('牛顿法计算时间 %f\n',NewtonToc);
-fprintf('加速比为%f\n',NewtonToc/FortiscueToc);
-VoltpA=Vp3m(1:3:end).*exp(1j*Vp3a(1:3:end));
-VoltpB=Vp3m(2:3:end).*exp(1j*Vp3a(2:3:end));
-VoltpC=Vp3m(3:3:end).*exp(1j*Vp3a(3:3:end));
-
+%% 用牛顿法求解begin
+% [r,c,GB]=find(phaseABCY);
+% Y=abs(phaseABCY);
+% Yangle=angle(GB);
+% Vp3=sparse(ones(busNum*3,1));%给电压赋初值
+% Vp3(2:3:end)=Vp3(2:3:end)*exp(1j*-120/180*pi);
+% Vp3(3:3:end)=Vp3(3:3:end)*exp(1j*+120/180*pi);
+% PQi3P=zeros(length(PQi)*3,1);
+% PQi3P(1:3:end)=(PQi-1)*3+1;
+% PQi3P(2:3:end)=(PQi-1)*3+2;
+% PQi3P(3:3:end)=(PQi-1)*3+3;
+% PG=0;
+% QG=0;
+% PD3P=sparse(busNum*3,1);
+% QD3P=sparse(busNum*3,1);
+% PD3P(1:3:end)=phaseASpotLoadP;
+% PD3P(2:3:end)=phaseBSpotLoadP;
+% PD3P(3:3:end)=phaseCSpotLoadP;
+% QD3P(1:3:end)=phaseASpotLoadQ;
+% QD3P(2:3:end)=phaseBSpotLoadQ;
+% QD3P(3:3:end)=phaseCSpotLoadQ;
+% QGi3P=zeros(length(QGi)*3,1);
+% QGi3P(1:3:end)=(QGi-1)*3+1;
+% QGi3P(2:3:end)=(QGi-1)*3+2;
+% QGi3P(3:3:end)=(QGi-1)*3+3;
+% Vp3m=abs(Vp3);
+% Vp3a=angle(Vp3);
+% Balance3P=zeros(length(Balance)*3,1);
+% Balance3P(1:3:end)=(Balance-1)*3+1;
+% Balance3P(2:3:end)=(Balance-1)*3+2;
+% Balance3P(3:3:end)=(Balance-1)*3+3;
+% Vp3a((Balance-1)*3+1)=0;
+% Vp3a((Balance-1)*3+2)=-120/180*pi;
+% Vp3a((Balance-1)*3+3)=+120/180*pi;
+% k=0;
+% maxD=10000;
+% tic
+% while(k<=kmax && maxD> EPS)
+%     k=k+1;
+%     [dP, dQ, YdotSinVolt, YdotCosVolt, diag_Volt_YdotSin, diag_Volt_YdotCos]=Unbalance(Balance3P,busNum*3, ...
+%         PQi3P,PG,QG,QGi3P,PD3P,QD3P,Vp3m,Vp3a,Y,Yangle,r,c,0,0,0,0);
+%     maxD=max(abs([dP;dQ;]));
+%     jaco=Jacobi(Balance3P,busNum*3,QGi3P,Vp3m,YdotSinVolt,YdotCosVolt,diag_Volt_YdotSin,diag_Volt_YdotCos);%雅克比矩阵
+%     [dV, dVangle]=Solv(busNum*3,jaco,dP,dQ);%解出修正量
+%     [Vp3m, Vp3a]=Modify(Vp3m,Vp3a,dV,dVangle,1);
+%     fprintf('第 %d 次迭代， 最大不平衡量为 %f\n',k,full(maxD));
+%     fprintf('迭代时间%f\n',toc);
+% end
+% NewtonToc=toc;
+% fprintf('牛顿法计算时间 %f\n',NewtonToc);
+% fprintf('加速比为%f\n',NewtonToc/FortiscueToc);
+% VoltpA=Vp3m(1:3:end).*exp(1j*Vp3a(1:3:end));
+% VoltpB=Vp3m(2:3:end).*exp(1j*Vp3a(2:3:end));
+% VoltpC=Vp3m(3:3:end).*exp(1j*Vp3a(3:3:end));
+%% 用牛顿法求解end
+%% 开始进入状态估计
+%状态量
+SEVoltpA=sparse(ones(busNum,1));
+SEVoltpB=sparse(ones(busNum,1)).*exp(1j*-120/180*pi);
+SEVoltpC=sparse(ones(busNum,1)).*exp(1j*+120/180*pi);
+SEphaseASpotLoadP=zeros(length(phaseASpotLoadP),1);
+SEphaseBSpotLoadP=zeros(length(phaseBSpotLoadP),1);
+SEphaseCSpotLoadP=zeros(length(phaseCSpotLoadP),1);
+SEphaseASpotLoadQ=zeros(length(phaseASpotLoadQ),1);
+SEphaseBSpotLoadQ=zeros(length(phaseBSpotLoadQ),1);
+SEphaseCSpotLoadQ=zeros(length(phaseCSpotLoadQ),1);
+SEVmf1=sparse(ones(busNum,1));
+SEVaf1=sparse(zeros(busNum,1));
+SEPD=sparse(length(Loadi),1);
+SEQD=sparse(length(Loadi),1);
+KK=0;
+plotGap=zeros(1,60);
+%初始化
+%状态量为 SEPD SEQD SEVmf1 SEVaf1 
+RestraintCount=length(SEVmf1)+length(Loadi)*2;
+ContrlCount=length(SEVmf1)*2+length(Loadi)*2;
+CenterA=0.1;
+Init_Z=sparse(ones(RestraintCount,1));
+Init_W=sparse(-1*ones(RestraintCount,1));
+Init_L=1*sparse(ones(RestraintCount,1));
+Init_U=1*sparse(ones(RestraintCount,1));
+Init_Y=sparse(2*busNum,1);%等式约束乘子
+Gap=(Init_L'*Init_Z-Init_U'*Init_W);
+Init_u=Gap/2/RestraintCount*CenterA;
+deltH=func_deltH(busNum,SEVmf1,fsY1amp,SEVaf1,r,c,fsY1ang,Loadi);
+deltG=func_deltG(busNum,Loadi);
+L_1Z=diag(Init_Z./Init_L);
+U_1W=diag(Init_W./Init_U);
+deltdeltF=func_deltdeltF(SEPD,ContrlCount);
+ddh=func_ddh(SEVmf1,Init_Y,busNum,fsY1amp,SEVaf1,r,c,fsY1ang,Loadi,ContrlCount);
+ddg=func_ddg();
+deltF=func_deltF(SEPD,ContrlCount);
+Luu=Init_U.*Init_W+Init_u*ones(RestraintCount,1);
+Lul=Init_L.*Init_Z-Init_u*ones(RestraintCount,1);
+Mat_G=FormG(SEVaf1,SEPD,SEQD);
+Mat_H=FormH(busNum,SEVmf1,PG,SEPD,QG,SEQD,fsY1amp,SEVaf1,r,c,fsY1ang);
+Ly=Mat_H;
+Lz=FormLz(Mat_G,Init_L,busNum,Loadi);
+Lw=FormLw(Mat_G,Init_U,busNum,Loadi);
+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,Lx,Balance,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,SEVmf1,SEVaf1,ContrlCount,Balance,busNum,PD,QD,Loadi);
+Gap=(Init_L'*Init_Z-Init_U'*Init_W);
\ No newline at end of file