diff --git a/AssignXX.m b/AssignXX.m
index 9353ff3..4d50403 100644
--- a/AssignXX.m
+++ b/AssignXX.m
@@ -1,18 +1,13 @@
 function [deltX,deltY]=AssignXX(XX,ContrlCount,Loadi,Balance,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);
-t=size(Loadi,1)*2;
+deltX=XX(1:Busnum*3*2);
+t=0;
 deltX(t+(Balance-1)*3+1)=0;
 deltX(t+(Balance-1)*3+2)=0;
 deltX(t+(Balance-1)*3+3)=0;
-t=size(Loadi,1)*2+Busnum*3;
+t=Busnum*3;
 deltX(t+(Balance-1)*3+1)=0;
 deltX(t+(Balance-1)*3+2)=0;
 deltX(t+(Balance-1)*3+3)=0;
-deltY=XX(ContrlCount+1:end);
+t=Busnum*3*2;
+deltY=XX(t+1:end);
 end
\ No newline at end of file
diff --git a/FormH.m b/FormH.m
index bd64484..480256e 100644
--- a/FormH.m
+++ b/FormH.m
@@ -6,8 +6,10 @@ function Mat_H=FormH(Busnum,Volt,PG,PD,QG,QD,Y,UAngel,r,c,Angle,Loadi)
 %% 
 %PD(Loadi)=QD(Loadi)./tan(acos(0.98));
 AngleIJ=sparse(r,c,UAngel(r)-UAngel(c)-Angle,Busnum*3,Busnum*3);
-dP=PG-sparse(Loadi,1,PD,Busnum*3,1)-diag(Volt)*Y.*cos(AngleIJ)*Volt;
-dQ=QG-sparse(Loadi,1,QD,Busnum*3,1)-diag(Volt)*Y.*sin(AngleIJ)*Volt;
+% dP=-sparse(Loadi,1,PD,Busnum*3,1)-diag(Volt)*Y.*cos(AngleIJ)*Volt;
+% dQ=-sparse(Loadi,1,QD,Busnum*3,1)-diag(Volt)*Y.*sin(AngleIJ)*Volt;
+dP=diag(Volt)*Y.*cos(AngleIJ)*Volt;
+dQ=diag(Volt)*Y.*sin(AngleIJ)*Volt;
 
 Mat_H=[dP;dQ;];
 Mat_H=sparse(Mat_H);
diff --git a/Modification.m b/Modification.m
index 0a8b7f4..1151509 100644
--- a/Modification.m
+++ b/Modification.m
@@ -1,12 +1,9 @@
 function [Init_Y,PG,QG,Volt,UAngel,PD,QD]=Modification(Init_Y,deltX,deltY,PG,QG,Volt,UAngel,ContrlCount,Balance,Busnum,PD,QD,Loadi)
-Init_Y=Init_Y+deltY';
+Init_Y=Init_Y+deltY;
 t=deltX(1:size(Loadi,1)*2);
 PD=PD+t(1:length(Loadi));
 QD=QD+t(length(Loadi)+1:length(Loadi)*2);
-t=deltX(size(Loadi,1)*2+1:ContrlCount)';
-t(Busnum*3+(Balance-1)*3+1)=0;
-t(Busnum*3+(Balance-1)*3+2)=0;
-t(Busnum*3+(Balance-1)*3+3)=0;
+t=deltX(1:Busnum*3*2)';
 %Volt=Volt+AlphaP*t(2:2:2*Busnum);暂时改一下20111227
 %UAngel=UAngel+AlphaP*t(1:2:2*Busnum);暂时改一下20111227
 % balVolt=Volt( (Balance-1)*3+1 );
@@ -14,5 +11,8 @@ Volt=Volt+t(1:Busnum*3)';
 Volt( (Balance-1)*3+1 )=1;
 Volt( (Balance-1)*3+2 )=1;
 Volt( (Balance-1)*3+3 )=1;
+t(Busnum*3+(Balance-1)*3+1)=0;
+t(Busnum*3+(Balance-1)*3+2)=0;
+t(Busnum*3+(Balance-1)*3+3)=0;
 UAngel=UAngel+t(Busnum*3+1:2*Busnum*3)';
 end
\ No newline at end of file
diff --git a/OPF.m b/OPF.m
index 199933a..382e8ea 100644
--- a/OPF.m
+++ b/OPF.m
@@ -3,10 +3,10 @@ function [JMeasurement,AME_Volt,AME_VAngle,AME_PD,AME_QD,AME_mVolt,AME_mPD,AME_m
 tic
 clc
 clear
-lineZ=readLineZ('./../DistributionNetwork-Power2Current/modified-feeder69\lineParameter.txt');
+lineZ=readLineZ('feeder13\lineParameter.txt');
 [ fsY0, fsY1, fsY2,phaseASpotLoadP,phaseBSpotLoadP,phaseCSpotLoadP ...
     phaseASpotLoadQ,phaseBSpotLoadQ,phaseCSpotLoadQ,setIJ,nodeNum,Balance,phaseABCY ...
-    cap]=dataRead(lineZ,'./../DistributionNetwork-Power2Current/modified-feeder69\data-某相没负荷.txt');
+    cap]=dataRead(lineZ,'feeder13\data1.txt');
 % phaseASpotLoadP(phaseASpotLoadP==0)=0.002;
 % phaseBSpotLoadP(phaseBSpotLoadP==0)=0.002;
 % phaseCSpotLoadP(phaseCSpotLoadP==0)=0.002;
@@ -15,36 +15,36 @@ lineZ=readLineZ('./../DistributionNetwork-Power2Current/modified-feeder69\linePa
 % phaseCSpotLoadQ(phaseCSpotLoadQ==0)=0.002;
 
 %负荷不平衡
-for I=1:length(phaseASpotLoadP)
-    roll=mod(round(rand()*10),3);
-    if roll==0
-        phaseBSpotLoadP(I)=phaseBSpotLoadP(I)+phaseASpotLoadP(I)*.45;
-        phaseCSpotLoadP(I)=phaseCSpotLoadP(I)+phaseASpotLoadP(I)*.45;
-        phaseASpotLoadP(I)=phaseASpotLoadP(I)*.10;
-        
-        phaseBSpotLoadQ(I)=phaseBSpotLoadQ(I)+phaseASpotLoadQ(I)*.45;
-        phaseCSpotLoadQ(I)=phaseCSpotLoadQ(I)+phaseASpotLoadQ(I)*.45;
-        phaseASpotLoadQ(I)=phaseASpotLoadQ(I)*.10;
-    end
-    if roll==1
-        phaseASpotLoadP(I)=phaseASpotLoadP(I)+phaseBSpotLoadP(I)*.45;
-        phaseCSpotLoadP(I)=phaseCSpotLoadP(I)+phaseBSpotLoadP(I)*.45;
-        phaseBSpotLoadP(I)=phaseBSpotLoadP(I)*.10;
-        
-        phaseASpotLoadQ(I)=phaseASpotLoadQ(I)+phaseBSpotLoadQ(I)*.45;
-        phaseCSpotLoadQ(I)=phaseCSpotLoadQ(I)+phaseBSpotLoadQ(I)*.45;
-        phaseBSpotLoadQ(I)=phaseBSpotLoadQ(I)*.10;
-    end
-    if roll==2
-        phaseASpotLoadP(I)=phaseASpotLoadP(I)+phaseCSpotLoadP(I)*.45;
-        phaseBSpotLoadP(I)=phaseBSpotLoadP(I)+phaseCSpotLoadP(I)*.45;
-        phaseCSpotLoadP(I)=phaseCSpotLoadP(I)*.10;
-        
-        phaseASpotLoadQ(I)=phaseASpotLoadQ(I)+phaseCSpotLoadQ(I)*.45;
-        phaseBSpotLoadQ(I)=phaseBSpotLoadQ(I)+phaseCSpotLoadQ(I)*.45;
-        phaseCSpotLoadQ(I)=phaseCSpotLoadQ(I)*.10;
-    end
-end
+% for I=1:length(phaseASpotLoadP)
+%     roll=mod(round(rand()*10),3);
+%     if roll==0
+%         phaseBSpotLoadP(I)=phaseBSpotLoadP(I)+phaseASpotLoadP(I)*.45;
+%         phaseCSpotLoadP(I)=phaseCSpotLoadP(I)+phaseASpotLoadP(I)*.45;
+%         phaseASpotLoadP(I)=phaseASpotLoadP(I)*.10;
+%         
+%         phaseBSpotLoadQ(I)=phaseBSpotLoadQ(I)+phaseASpotLoadQ(I)*.45;
+%         phaseCSpotLoadQ(I)=phaseCSpotLoadQ(I)+phaseASpotLoadQ(I)*.45;
+%         phaseASpotLoadQ(I)=phaseASpotLoadQ(I)*.10;
+%     end
+%     if roll==1
+%         phaseASpotLoadP(I)=phaseASpotLoadP(I)+phaseBSpotLoadP(I)*.45;
+%         phaseCSpotLoadP(I)=phaseCSpotLoadP(I)+phaseBSpotLoadP(I)*.45;
+%         phaseBSpotLoadP(I)=phaseBSpotLoadP(I)*.10;
+%         
+%         phaseASpotLoadQ(I)=phaseASpotLoadQ(I)+phaseBSpotLoadQ(I)*.45;
+%         phaseCSpotLoadQ(I)=phaseCSpotLoadQ(I)+phaseBSpotLoadQ(I)*.45;
+%         phaseBSpotLoadQ(I)=phaseBSpotLoadQ(I)*.10;
+%     end
+%     if roll==2
+%         phaseASpotLoadP(I)=phaseASpotLoadP(I)+phaseCSpotLoadP(I)*.45;
+%         phaseBSpotLoadP(I)=phaseBSpotLoadP(I)+phaseCSpotLoadP(I)*.45;
+%         phaseCSpotLoadP(I)=phaseCSpotLoadP(I)*.10;
+%         
+%         phaseASpotLoadQ(I)=phaseASpotLoadQ(I)+phaseCSpotLoadQ(I)*.45;
+%         phaseBSpotLoadQ(I)=phaseBSpotLoadQ(I)+phaseCSpotLoadQ(I)*.45;
+%         phaseCSpotLoadQ(I)=phaseCSpotLoadQ(I)*.10;
+%     end
+% end
 
 
 %% 潮流计算begin
@@ -143,7 +143,6 @@ while(k<=kmax+10 && maxD> EPS)
     %Vf2=fsY2\If2;
     Vf2=fsY2Q*(fsY2U\(fsY2L\(fsY2P*(fsY2R\If2))));
     fprintf('迭代时间%f\n',toc);
-    
     %
 end
 FortiscueToc=toc;
@@ -181,14 +180,22 @@ PQi3P(1:3:end)=(Loadi-1)*3+1;
 PQi3P(2:3:end)=(Loadi-1)*3+2;
 PQi3P(3:3:end)=(Loadi-1)*3+3;
 Loadi=PQi3P;
+%把发电机也作为负荷
+Loadi=[Loadi;(Balance-1)*3+1;(Balance-1)*3+2;(Balance-1)*3+3];
 PD3P=sparse(Busnum*3,1);
 QD3P=sparse(Busnum*3,1);
-PD3P(1:3:end)=phaseASpotLoadP*0.9;
-PD3P(2:3:end)=phaseBSpotLoadP*0.9;
-PD3P(3:3:end)=phaseCSpotLoadP*0.9;
-QD3P(1:3:end)=phaseASpotLoadQ*0.9;
-QD3P(2:3:end)=phaseBSpotLoadQ*0.9;
-QD3P(3:3:end)=phaseCSpotLoadQ*0.9;
+phaseASpotLoadP(Balance)=-real(PGQG((Balance-1)*3+1));
+phaseBSpotLoadP(Balance)=-real(PGQG((Balance-1)*3+2));
+phaseCSpotLoadP(Balance)=-real(PGQG((Balance-1)*3+3));
+phaseASpotLoadQ(Balance)=-imag(PGQG((Balance-1)*3+1));
+phaseBSpotLoadQ(Balance)=-imag(PGQG((Balance-1)*3+2));
+phaseCSpotLoadQ(Balance)=-imag(PGQG((Balance-1)*3+3));
+% PD3P(1:3:end)=phaseASpotLoadP*0.9;
+% PD3P(2:3:end)=phaseBSpotLoadP*0.9;
+% PD3P(3:3:end)=phaseCSpotLoadP*0.9;
+% QD3P(1:3:end)=phaseASpotLoadQ*0.9;
+% QD3P(2:3:end)=phaseBSpotLoadQ*0.9;
+% QD3P(3:3:end)=phaseCSpotLoadQ*0.9;
 PD3P=PD3P(Loadi);
 QD3P=QD3P(Loadi);
 QGi3P=zeros(length(QGi)*3,1);
@@ -274,7 +281,7 @@ wQD(ismember( Loadi,noPQi3P))=1./(abs(mQD3P(ismember( Loadi,noPQi3P))*0.15).^2);
 %%
 % RestraintCount=size(Loadi,1)*2+length(rVoltABCV); %约束条件数
 % RestraintCount=size(Loadi,1)*2; %约束条件数
-Init_Y=sparse(1,2*Busnum*3);%与学姐一致
+Init_Y=sparse(length(setdiff(1:Busnum*3,Loadi))*2,1);%与学姐一致
 KK=0;
 ContrlCount=size(Loadi,1)*2+Busnum*6;
 kmax=200;
@@ -294,7 +301,7 @@ while(maxD>Precision)
     end
     %% 开始计算OPF
     %% 形成等式约束的雅克比
-    deltH=func_deltH(Busnum,Volt,Y,UAngel,r,c,Yangle,Loadi);
+    dH_dx = jacobian_M(Busnum,Volt,Y,Yangle,UAngel,r,c);  %形成雅克比矩阵
     %% 形成不等式约束的雅克比
     %%  开始构建ddg
     %%  开始构建deltF
@@ -302,8 +309,8 @@ while(maxD>Precision)
     %% 形成方程矩阵
     Mat_H=FormH(Busnum,Volt,PG3P,PD3P,QG3P,QD3P,Y,UAngel,r,c,Yangle,Loadi);
     %% 开始解方程
-    %     fprintf('迭代次数 %d Gap %f\n',KK+1,plotGap(KK+1));
-    XX=SolveIt(deltH,ContrlCount,Balance,Busnum,Loadi,deltF,mPD3P,PD3P,mQD3P,QD3P,mVoltABCV,Volt,Mat_H,wVolt,wPD,wQD);
+%         fprintf('迭代次数 %d Gap %f\n',KK+1,plotGap(KK+1));
+    XX=SolveIt(ContrlCount,Balance,Busnum,Loadi,deltF,mPD3P,PD3P,mQD3P,QD3P,mVoltABCV,Volt,Mat_H,wVolt,wPD,wQD,dH_dx);
     %%取各分量
     [deltX,deltY]=AssignXX(XX,ContrlCount,Loadi,Balance,Busnum);
     [Init_Y,PG,QG,Volt,UAngel,PD3P,QD3P]=Modification(Init_Y,deltX,deltY,PG,QG,Volt,UAngel,ContrlCount,Balance,Busnum,PD3P,QD3P,Loadi);
@@ -313,6 +320,8 @@ while(maxD>Precision)
 end
 totalTime=toc
 (rVoltABCV-Volt)./rVoltABCV*100;
+PD3P=-Mat_H(Loadi);
+QD3P=-Mat_H(Loadi+Busnum*3);
 %% 计算统计量
 %目标函数均值
 JMeasurement=sum(((mVoltABCV-Volt)./mVoltABCV./sigma).^2)+sum(((mPD3P-PD3P)./mPD3P./sigma).^2)+sum(((mQD3P-QD3P)./mQD3P./sigma).^2);
diff --git a/SolveIt.m b/SolveIt.m
index 8483710..e552455 100644
--- a/SolveIt.m
+++ b/SolveIt.m
@@ -1,59 +1,66 @@
-function XX=SolveIt(deltH,ContrlCount,Balance,Busnum,Loadi,deltF,PD0,PD,QD0,QD,mVolt,Volt,Mat_H,wVolt,wPD,wQD)
-% LxComa=FormLxComa(Lx);%Lx=deltF-deltH*Init_Y';
-% Lx=LxComa;
-% H=-deltdeltF+ddh;
-% aa=[
-%     H,deltH;
-%     deltH',zeros(size(Init_Y,2));
-%     ];
-aa=[
-    deltF'*diag([wPD ;wQD ;wVolt])*deltF deltH'
-    deltH zeros(size(deltH,1),size(deltH',2))
-];
-% yy=[Lx;-Ly];
-% t3=2*wPD.*(PD-PD0);
-% t4=2*wQD.*(QD-QD0);
-% t5=2*wVolt.*(Volt-mVolt);
+function XX=SolveIt(ContrlCount,Balance,Busnum,Loadi,deltF,PD0,PD,QD0,QD,mVolt,Volt,Mat_H,wVolt,wPD,wQD,dH_dx)
+% dH_dx
+dP_x=-dH_dx(Loadi,:);
+dQ_x=-dH_dx(Loadi+Busnum*3,:);
+dV_x=[sparse(1:Busnum*3,1:Busnum*3,1,Busnum*3,Busnum*3),sparse(Busnum*3,Busnum*3)];
+H=[dP_x;dQ_x;dV_x];
+c=dH_dx(setdiff(1:Busnum*3,Loadi),:);
+c=[c;dH_dx(setdiff(1:Busnum*3,Loadi)+Busnum*3,:)];
+aa=[H'*diag([wPD;wQD;wVolt])*H c'
+    c zeros(size(c,1),size(c',2))];
+% aa=H'*eye(length(H))*H;
+PD=-Mat_H(Loadi,:);
+QD=-Mat_H(Loadi+Busnum*3,:);
 deltZ=[(PD0-PD);
     (QD0-QD);
     (mVolt-Volt);
     ];
 yy=[
-    deltF'*diag([wPD; wQD; wVolt])*deltZ;
--Mat_H    ;
+    H'*diag([wPD;wQD;wVolt])*deltZ;
+    -Mat_H(setdiff(1:Busnum*3,Loadi));
+    -Mat_H(setdiff(1:Busnum*3,Loadi)+Busnum*3);
 ];
+% yy=[
+%     H'*eye(length(H))*deltZ;
+% ];
 %% 平衡节点电压不变 
-t=size(Loadi,1)*2;
+t=0;
 aa(t+(Balance-1)*3+1,:)=0;
 aa(t+(Balance-1)*3+2,:)=0;
 aa(t+(Balance-1)*3+3,:)=0;
 aa(:,t+(Balance-1)*3+1)=0;
 aa(:,t+(Balance-1)*3+2)=0;
 aa(:,t+(Balance-1)*3+3)=0;
-%aa(t+Balance,t+Balance)=1;
-aa=aa+sparse(t+(Balance-1)*3+1,t+(Balance-1)*3+1,ones(length(Balance),1),ContrlCount+2*Busnum*3,ContrlCount+2*Busnum*3);
-aa=aa+sparse(t+(Balance-1)*3+2,t+(Balance-1)*3+2,ones(length(Balance),1),ContrlCount+2*Busnum*3,ContrlCount+2*Busnum*3);
-aa=aa+sparse(t+(Balance-1)*3+3,t+(Balance-1)*3+3,ones(length(Balance),1),ContrlCount+2*Busnum*3,ContrlCount+2*Busnum*3);
+aa=aa+sparse(t+(Balance-1)*3+1,t+(Balance-1)*3+1,ones(length(Balance),1),size(c,1)+2*Busnum*3,size(c,1)+2*Busnum*3);
+aa=aa+sparse(t+(Balance-1)*3+2,t+(Balance-1)*3+2,ones(length(Balance),1),size(c,1)+2*Busnum*3,size(c,1)+2*Busnum*3);
+aa=aa+sparse(t+(Balance-1)*3+3,t+(Balance-1)*3+3,ones(length(Balance),1),size(c,1)+2*Busnum*3,size(c,1)+2*Busnum*3);
+
+% aa=aa+sparse(t+(Balance-1)*3+1,t+(Balance-1)*3+1,ones(length(Balance),1),2*Busnum*3,2*Busnum*3);
+% aa=aa+sparse(t+(Balance-1)*3+2,t+(Balance-1)*3+2,ones(length(Balance),1),2*Busnum*3,2*Busnum*3);
+% aa=aa+sparse(t+(Balance-1)*3+3,t+(Balance-1)*3+3,ones(length(Balance),1),2*Busnum*3,2*Busnum*3);
+
 %%
-t=size(Loadi,1)*2+Busnum*3;
+t=Busnum*3;
 aa(t+(Balance-1)*3+1,:)=0;
 aa(t+(Balance-1)*3+2,:)=0;
 aa(t+(Balance-1)*3+3,:)=0;
 aa(:,t+(Balance-1)*3+1)=0;
 aa(:,t+(Balance-1)*3+2)=0;
 aa(:,t+(Balance-1)*3+3)=0;
-%aa(t+Balance,t+Balance)=1;
-aa=aa+sparse(t+(Balance-1)*3+1,t+(Balance-1)*3+1,ones(length(Balance),1),ContrlCount+2*Busnum*3,ContrlCount+2*Busnum*3);
-aa=aa+sparse(t+(Balance-1)*3+2,t+(Balance-1)*3+2,ones(length(Balance),1),ContrlCount+2*Busnum*3,ContrlCount+2*Busnum*3);
-aa=aa+sparse(t+(Balance-1)*3+3,t+(Balance-1)*3+3,ones(length(Balance),1),ContrlCount+2*Busnum*3,ContrlCount+2*Busnum*3);
+aa=aa+sparse(t+(Balance-1)*3+1,t+(Balance-1)*3+1,ones(length(Balance),1),size(c,1)+2*Busnum*3,size(c,1)+2*Busnum*3);
+aa=aa+sparse(t+(Balance-1)*3+2,t+(Balance-1)*3+2,ones(length(Balance),1),size(c,1)+2*Busnum*3,size(c,1)+2*Busnum*3);
+aa=aa+sparse(t+(Balance-1)*3+3,t+(Balance-1)*3+3,ones(length(Balance),1),size(c,1)+2*Busnum*3,size(c,1)+2*Busnum*3);
+% aa=aa+sparse(t+(Balance-1)*3+1,t+(Balance-1)*3+1,ones(length(Balance),1),2*Busnum*3,2*Busnum*3);
+% aa=aa+sparse(t+(Balance-1)*3+2,t+(Balance-1)*3+2,ones(length(Balance),1),2*Busnum*3,2*Busnum*3);
+% aa=aa+sparse(t+(Balance-1)*3+3,t+(Balance-1)*3+3,ones(length(Balance),1),2*Busnum*3,2*Busnum*3);
 %%
 dxdy=aa\yy;
 %% KLU
 %spy(aa)
 %dxdy = klu(aa,'\',full(yy));
 %%
-dX=dxdy(1:ContrlCount);
-dY=dxdy(ContrlCount+1:end);
+dX=dxdy(1:Busnum*3*2);
+dY=dxdy(Busnum*3*2+1:end);
 XX=[
     dX;
     dY;    
diff --git a/feeder13/data3.txt b/feeder13/data3.txt
new file mode 100644
index 0000000..0e6873d
--- /dev/null
+++ b/feeder13/data3.txt
@@ -0,0 +1,60 @@
+650 4.16
+0
+600	632	645	50
+602	632	633	50
+601	633	634	60
+600	645	646	30
+601	650	632	200
+602	684	652	80
+601	632	671	200
+600	671	684	30
+600	671	692	30
+602	684	611	30
+601	692	675	50
+0
+692 0 0 0
+0
+684	120	90	120	90	160	110
+633	120	90	120	90	160	110
+632	120	90	120	90	160	110
+634	120	90	120	90	160	110
+645	40	30	50	40	60	50
+646	100	50	70	40	60	40
+652	20	10	60	30	40	40
+671	385	220	385	220	385	220
+675	68	60	290	212	285	190
+692	260	200	40	40	70	40
+611	80	20	60	40	30	20
+0
+600	632	645	500
+602	632	633	500
+601	633	634	600
+600	645	646	300
+601	650	632	2000
+602	684	652	800
+601	632	671	2000
+600	671	684	300
+602	671	680	1000
+600	671	692	300
+602	684	611	300
+601	692	675	500
+0
+632 0 0 0
+0
+634 160 110 120 90 120 90
+645 123 102 170 125 133 110
+646 172 150 230 132 0 211 160
+652 128 86 30 10 80 70
+671 385 220 385 220 385 220
+675 485 190 68 60 290 212
+692 142 140 153 121 170 151
+611 156 142 185 60 170 80
+0
+634 160 110 120 90 120 90
+645 123 102 170 125 133 110
+646 17 15 23 13 21 16
+652 128 86 30 10 80 70
+675 85 10 68 60 20 12
+611 52 42 85 60 70 40
+680 12 4 10 6 7 4
+
diff --git a/jacobian_M.m b/jacobian_M.m
index fb999de..37aa10a 100644
--- a/jacobian_M.m
+++ b/jacobian_M.m
@@ -8,13 +8,12 @@ function [Jacob]=jacobian_M3(Busnum,Volt,Y,Angle,UAngel,r,c)
 %% 分别求雅克比矩阵的子阵H,L,N,J及有功无功分量P,Q
 AngleIJ=UAngel(r)-UAngel(c)-Angle;
 mat_AngleIJ=sparse(r,c,AngleIJ,Busnum*3,Busnum*3);
-mat_IvAngleIJ=mat_AngleIJ';
-H=sparse(1:Busnum*3,1:Busnum*3,Volt,Busnum*3,Busnum*3)*Y.*sin(mat_IvAngleIJ)*sparse(1:Busnum*3,1:Busnum*3,Volt,Busnum*3,Busnum*3)-sparse(1:Busnum*3,1:Busnum*3,Y.*sin(mat_AngleIJ)*Volt,Busnum*3,Busnum*3)*sparse(1:Busnum*3,1:Busnum*3,Volt,Busnum*3,Busnum*3);
-N=-sparse(1:Busnum*3,1:Busnum*3,Volt,Busnum*3,Busnum*3)*Y.*cos(mat_IvAngleIJ)*sparse(1:Busnum*3,1:Busnum*3,Volt,Busnum*3,Busnum*3)+sparse(1:Busnum*3,1:Busnum*3,Y.*cos(mat_AngleIJ)*Volt,Busnum*3,Busnum*3)*sparse(1:Busnum*3,1:Busnum*3,Volt,Busnum*3,Busnum*3);
-J=sparse(1:Busnum*3,1:Busnum*3,Y.*cos(mat_AngleIJ)*Volt,Busnum*3,Busnum*3)+Y.*cos(mat_IvAngleIJ)*sparse(1:Busnum*3,1:Busnum*3,Volt,Busnum*3,Busnum*3);
-L=sparse(1:Busnum*3,1:Busnum*3,Y.*sin(mat_AngleIJ)*Volt,Busnum*3,Busnum*3)+Y.*sin(mat_IvAngleIJ)*sparse(1:Busnum*3,1:Busnum*3,Volt,Busnum*3,Busnum*3);
-t1=[J,L;
-    H,N;
-]';
-Jacob=-t1;
+H=sparse(1:Busnum*3,1:Busnum*3,Volt,Busnum*3,Busnum*3)*Y.*sin(mat_AngleIJ)*sparse(1:Busnum*3,1:Busnum*3,Volt,Busnum*3,Busnum*3)-sparse(1:Busnum*3,1:Busnum*3,Y.*sin(mat_AngleIJ)*Volt,Busnum*3,Busnum*3)*sparse(1:Busnum*3,1:Busnum*3,Volt,Busnum*3,Busnum*3);
+N=-sparse(1:Busnum*3,1:Busnum*3,Volt,Busnum*3,Busnum*3)*Y.*cos(mat_AngleIJ)*sparse(1:Busnum*3,1:Busnum*3,Volt,Busnum*3,Busnum*3)+sparse(1:Busnum*3,1:Busnum*3,Y.*cos(mat_AngleIJ)*Volt,Busnum*3,Busnum*3)*sparse(1:Busnum*3,1:Busnum*3,Volt,Busnum*3,Busnum*3);
+J=sparse(1:Busnum*3,1:Busnum*3,Y.*cos(mat_AngleIJ)*Volt,Busnum*3,Busnum*3)+Y.*cos(mat_AngleIJ)*sparse(1:Busnum*3,1:Busnum*3,Volt,Busnum*3,Busnum*3);
+L=sparse(1:Busnum*3,1:Busnum*3,Y.*sin(mat_AngleIJ)*Volt,Busnum*3,Busnum*3)+Y.*sin(mat_AngleIJ)*sparse(1:Busnum*3,1:Busnum*3,Volt,Busnum*3,Busnum*3);
+t1=[J,H;
+    L,N;
+];
+Jacob=t1;
 end
\ No newline at end of file