diff --git a/Jacobi.m b/Jacobi.m
new file mode 100644
index 0000000..73f31fd
--- /dev/null
+++ b/Jacobi.m
@@ -0,0 +1,25 @@
+function [jaco]=Jacobi(Balance,busNum,QGi,Volt,YdotSinVolt,YdotCosVolt,diag_Volt_YdotSin,diag_Volt_YdotCos)
+diag_YdotSinVolt_=diag(YdotSinVolt);
+diag_YdotCosVolt_=diag(YdotCosVolt);
+dPdTyta=diag_Volt_YdotSin*diag(Volt)-diag_YdotSinVolt_*diag(Volt); % 简化第三次
+dQdTyta=-diag_Volt_YdotCos*diag(Volt)+diag_YdotCosVolt_*diag(Volt);%dQ/dThyta
+dPdV=diag_YdotCosVolt_+diag_Volt_YdotCos;%dP/dV
+dQdV=diag_YdotSinVolt_+diag_Volt_YdotSin;%dQ/dV
+%% 平衡节点相角不变
+dPdTyta(Balance,:)=0;
+dPdTyta(:,Balance)=0;
+dPdTyta=dPdTyta+sparse(Balance,Balance,ones(length(Balance),1),busNum,busNum);
+dQdTyta(:,Balance)=0;
+dPdV(Balance,:)=0;
+%% PV节点电压不变
+dQdV(QGi,:)=0;
+dQdV(:,QGi)=0;
+dQdV=dQdV+sparse(QGi,QGi,ones(length(QGi),1),busNum,busNum);
+dQdTyta(QGi,:)=0;
+dPdV(:,QGi)=0;
+% jaco=[    
+%     dPdV,dPdTyta;
+%     dQdV,dQdTyta;
+% ];
+jaco=cat(1,cat(2,dPdV,dPdTyta),cat(2,dQdV,dQdTyta));
+end
\ No newline at end of file
diff --git a/Modify.m b/Modify.m
new file mode 100644
index 0000000..7524f42
--- /dev/null
+++ b/Modify.m
@@ -0,0 +1,5 @@
+function [Volt, Vangle]=Modify(Volt,Vangle,dV,dVangle,optimalu)
+%Volt=Volt+Volt.*dV;
+Volt=Volt+optimalu*dV;
+Vangle=Vangle+optimalu*dVangle;
+end
\ No newline at end of file
diff --git a/OPF.m b/OPF.m
index cae9621..a42e3f1 100644
--- a/OPF.m
+++ b/OPF.m
@@ -5,6 +5,148 @@ lineZ=readLineZ('feeder13\lineParameter.txt');
 [ fsY0, fsY1, fsY2,phaseASpotLoadP,phaseBSpotLoadP,phaseCSpotLoadP ...
     phaseASpotLoadQ,phaseBSpotLoadQ,phaseCSpotLoadQ,setIJ,nodeNum,Balance,phaseABCY ...
     cap]=dataRead(lineZ,'feeder13\data1.txt');
+%% 潮流计算begin
+a=exp(1j*2*pi/3);
+Tp2f=1/3*[1 1 1;
+    1 a a^2;
+    1 a^2 a];
+Tp2f=sparse(Tp2f);
+Tf2p=inv(Tp2f);
+fsY1amp=abs(fsY1);
+[r,c,fsY1ang]=find(fsY1);
+fsY1ang=angle(fsY1ang);
+Pabc=phaseASpotLoadP+phaseBSpotLoadP+phaseCSpotLoadP;
+Qabc=phaseASpotLoadQ+phaseBSpotLoadQ+phaseCSpotLoadQ;
+busNum=length(phaseASpotLoadP);
+%给序电压赋初值
+Vmf1=sparse(ones(busNum,1));
+Vaf1=sparse(zeros(busNum,1));
+%先求解正序的
+PQi=nodeNum;
+PG=sparse(busNum,1);
+QG=sparse(busNum,1);
+QGi=[Balance];
+PD=Pabc/3;
+QD=Qabc/3;
+Loadi=find(PD~=0);
+maxD=100000;% 最大不平衡量
+EPS=1e-5;
+k=0;
+kmax=20;
+fsY11=fsY1;
+fsY00=fsY0;
+fsY22=fsY2;
+Vf2=sparse(busNum,1);
+If2=sparse(busNum,1);
+Vf0=sparse(busNum,1);
+If0=sparse(busNum,1);
+%准备序矩阵
+%平衡节点置0置1
+fsY2(Balance,:)=0;
+fsY2(:,Balance)=0;
+fsY2=fsY2+sparse(Balance,Balance,ones(length(Balance),1),busNum,busNum);
+%平衡节点置0置1
+fsY0(Balance,:)=0;
+fsY0(:,Balance)=0;
+fsY0=fsY0+sparse(Balance,Balance,ones(length(Balance),1),busNum,busNum);
+%%LU分解
+[fsY0L,fsY0U,fsY0P,fsY0Q,fsY0R]=lu(fsY0);
+[fsY2L,fsY2U,fsY2P,fsY2Q,fsY2R]=lu(fsY2);
+%算初始补偿功率
+tic
+VoltpA=sparse(ones(busNum,1));
+VoltpB=sparse(ones(busNum,1)).*exp(1j*-120/180*pi);
+VoltpC=sparse(ones(busNum,1)).*exp(1j*+120/180*pi);
+while(k<=kmax && maxD> EPS)
+    k=k+1;
+    %把补偿电容看作负荷
+    SA=VoltpA.*conj(VoltpA.*sparse(cap.capNode,1,1j*cap.capB(:,1),busNum,1));
+    SB=VoltpB.*conj(VoltpA.*sparse(cap.capNode,1,1j*cap.capB(:,2),busNum,1));
+    SC=VoltpC.*conj(VoltpA.*sparse(cap.capNode,1,1j*cap.capB(:,3),busNum,1));
+    %先不要电容
+%     SA=0;
+%     SB=0;
+%     SC=0;
+    iterPD=PD+real(SA+SB+SC)/3;
+    iterQD=QD+imag(SA+SB+SC)/3;
+    iterPhaseASpotLoadP=phaseASpotLoadP+real(SA);
+    iterPhaseBSpotLoadP=phaseBSpotLoadP+real(SB);
+    iterPhaseCSpotLoadP=phaseCSpotLoadP+real(SC);
+    iterPhaseASpotLoadQ=phaseASpotLoadQ+imag(SA);
+    iterPhaseBSpotLoadQ=phaseBSpotLoadQ+imag(SB);
+    iterPhaseCSpotLoadQ=phaseCSpotLoadQ+imag(SC);
+%    %全部转换为负荷电流
+%    CurpA=conj((iterPhaseASpotLoadP+1j*iterPhaseASpotLoadQ)./VoltpA);
+%    CurpB=conj((iterPhaseBSpotLoadP+1j*iterPhaseBSpotLoadQ)./VoltpB);
+%    CurpC=conj((iterPhaseCSpotLoadP+1j*iterPhaseCSpotLoadQ)./VoltpC);
+%    %转换为序电流
+%    f012=Tp2f*conj([CurpA';CurpB';CurpC']);
+%    %把三序电流分离出来
+%    If0=conj(f012(1,:)');
+%    If1=conj(f012(2,:)');
+%    If2=conj(f012(3,:)');
+%    %试着算一下正序电流
+%     fsY11*V1;
+%    %形成负荷序电流的测量值
+%    mIf0=If0;
+%    mIf1=If1;
+%    mIf1(3)=-mIf1(2);
+%    mIf2=If2;
+%    %计算
+%     fsY11=fsY11+sparse(Balance,Balance,ones(length(Balance),1),busNum,bus
+%     Num);%这里要置0，置1，否则是奇异的
+    %%做最小二乘法
+    [dP, dQ, YdotSinVolt, YdotCosVolt, diag_Volt_YdotSin, diag_Volt_YdotCos]=Unbalance(Balance,busNum, ...
+        PQi,PG,QG,QGi,iterPD,iterQD,Vmf1,Vaf1,fsY1amp,fsY1ang,r,c,Vf2,If2,Vf0,If0);%不平衡量
+    maxD=max(abs([dP;dQ;]));
+    jaco=Jacobi(Balance,busNum,QGi,Vmf1,YdotSinVolt,YdotCosVolt,diag_Volt_YdotSin,diag_Volt_YdotCos);%雅克比矩阵
+    [dV, dVangle]=Solv(busNum,jaco,dP,dQ);%解出修正量
+    [Vmf1, Vaf1]=Modify(Vmf1,Vaf1,dV,dVangle,1);
+    fprintf('第 %d 次迭代， 最大不平衡量为 %f\n',k,full(maxD));
+    %转换为三相电压
+    VoltpABC=Tp2f\conj([ Vf0'; (Vmf1.*exp(1j*Vaf1))'; Vf2']);%用Tp2f\ 代替Tf2p*
+    VoltpA=conj(VoltpABC(1,:)');
+    CurpA=-conj((iterPhaseASpotLoadP+1j*iterPhaseASpotLoadQ)./VoltpA);
+    VoltpB=conj(VoltpABC(2,:)');
+    CurpB=-conj((iterPhaseBSpotLoadP+1j*iterPhaseBSpotLoadQ)./VoltpB);
+    VoltpC=conj(VoltpABC(3,:)');
+    CurpC=-conj((iterPhaseCSpotLoadP+1j*iterPhaseCSpotLoadQ)./VoltpC);
+    f012=Tp2f*conj([CurpA';CurpB';CurpC']);
+    If0=conj(f012(1,:)');
+    If1=conj(f012(2,:)');
+    If2=conj(f012(3,:)');
+    
+    
+    If0(Balance)=0;
+    If2(Balance)=0;
+    %Vf0=fsY0\If0;
+    Vf0=fsY0Q*(fsY0U\(fsY0L\(fsY0P*(fsY0R\If0))));
+    %Vf2=fsY2\If2;
+    Vf2=fsY2Q*(fsY2U\(fsY2L\(fsY2P*(fsY2R\If2))));
+    fprintf('迭代时间%f\n',toc);
+    
+    %
+end
+FortiscueToc=toc;
+fprintf('Fortiscue法计算时间 %f\n',FortiscueToc);
+Vf1=Vmf1.*exp(1j*Vaf1);
+(Vf0.*conj(fsY00*Vf0)+Vf1.*conj(fsY11*Vf1)+Vf2.*conj(fsY22*Vf2))*3;%包含补偿电容的功率
+conj(Tf2p*[If0(2);If1(2);If2(2)]).*(Tf2p*[Vf0(2);Vf1(2);Vf2(2)]);
+IpABC=Tf2p*conj([If0';If1';If2']);
+%转换回三相电压
+VoltpABC=Tf2p*conj([ Vf0'; Vf1'; Vf2']);
+disp(['                  A                        B                          C'])
+full(abs(VoltpABC'))
+fprintf('节点号对应\n');
+disp([setIJ,nodeNum ])
+%%检查反推回去的功率是否满足
+ub=checkSSatisfied(Balance,phaseABCY,VoltpABC, ...
+    phaseASpotLoadP,phaseBSpotLoadP,phaseCSpotLoadP, ...
+    phaseASpotLoadQ,phaseBSpotLoadQ,phaseCSpotLoadQ );
+fprintf('最大不平衡量为%f\n\n',full(max(abs(ub))))
+%% 潮流计算end
+
+
 busNum=length(nodeNum);
 Busnum=busNum;
 PQi=setxor(nodeNum,Balance);
diff --git a/Solv.m b/Solv.m
new file mode 100644
index 0000000..f658297
--- /dev/null
+++ b/Solv.m
@@ -0,0 +1,16 @@
+function [dV, dVangle]=Solv(busNum,jaco,dP,dQ)
+
+%y=klu (jaco, '\', full(-[dP;dQ]));
+y=jaco\(-[dP;dQ]);
+% [L,U] = luinc(jaco,1e-3); %luinc(A,'0')也可以试一下，是一种完全不同的ILU
+% tol=1e-5; %残量的精度要求
+% restart=30; % 30-50之间吧，不要过小
+% maxit=100; %看情况，如果不收敛就适当调大
+% [x,flag]=gmres(jaco,-[dP;dQ],restart,tol,maxit);
+% y=sparse(x);
+
+y=sparse(y);
+dV=y(1:busNum);
+dVangle=y(busNum+1:end);
+
+end
\ No newline at end of file
diff --git a/Unbalance.m b/Unbalance.m
new file mode 100644
index 0000000..e8bfcf6
--- /dev/null
+++ b/Unbalance.m
@@ -0,0 +1,18 @@
+function [dP, dQ, YdotSinVmf1, YdotCosVmf1, diag_Vmf1_YdotSin, diag_Vmf1_YdotCos]=Unbalance(Balance,busNum, ... 
+    PQi,PG,QG,QGi,PD,QD,Vmf1,Vaf1,Y,Yangle,r,c,Vf2,If2,Vf0,If0)
+t1=sparse(r,c,Vaf1(r)-Vaf1(c) -Yangle,busNum,busNum) ;
+YdotSin=Y.* ( spfun(@sin,t1) );
+YdotCos=Y.* ( spfun (@cos, t1 ) );
+orderedPDPG=sparse(PQi,1,PD-PG,busNum,1); %节点功率按顺序排列
+orderedQDQG=sparse(PQi,1,QD-QG,busNum,1); %节点功率按顺序排列
+diag_Vmf1_YdotCos=diag(Vmf1)*YdotCos;
+diag_Vmf1_YdotSin=diag(Vmf1)*YdotSin;
+YdotCosVmf1=YdotCos*Vmf1;
+YdotSinVmf1=YdotSin*Vmf1;
+diag_Vmf1_YdotCosVmf1=diag_Vmf1_YdotCos*Vmf1;
+diag_Vmf1_YdotSinVmf1=diag_Vmf1_YdotSin*Vmf1;
+dP=diag_Vmf1_YdotCosVmf1+orderedPDPG+real(Vf2.*conj(If2)+Vf0.*conj(If0));
+dQ=diag_Vmf1_YdotSinVmf1+orderedQDQG+imag(Vf2.*conj(If2)+Vf0.*conj(If0));
+dP(Balance)=0; 
+dQ(QGi)=0;  
+end
\ No newline at end of file
diff --git a/checkSSatisfied.m b/checkSSatisfied.m
new file mode 100644
index 0000000..d92f537
--- /dev/null
+++ b/checkSSatisfied.m
@@ -0,0 +1,25 @@
+function [ output_args ] = checkSSatisfied(Balance,phaseABCY,VoltpABC,phaseASpotLoadP,phaseBSpotLoadP,phaseCSpotLoadP,phaseASpotLoadQ,phaseBSpotLoadQ,phaseCSpotLoadQ )
+%%利用Fortiscue方法计算得到了相电压和相角后，反推负荷功率，检查是否和给定的功率一致。
+busNum=size(phaseABCY,1)/3;
+VoltpA=VoltpABC(1,:);
+VoltpB=VoltpABC(2,:);
+VoltpC=VoltpABC(3,:);
+Vp3=sparse(busNum*3,1);
+Vp3(1:3:end)=VoltpA;
+Vp3(2:3:end)=VoltpB;
+Vp3(3:3:end)=VoltpC;
+Ip3=phaseABCY*Vp3;
+Sp3=Vp3.*conj(Ip3);
+%VoltpABC.*conj(IpABC);
+pLoadABC=sparse(length(phaseASpotLoadP)*3,1);
+pLoadABC(1:3:end)=phaseASpotLoadP+1j*phaseASpotLoadQ;
+pLoadABC(2:3:end)=phaseBSpotLoadP+1j*phaseBSpotLoadQ;
+pLoadABC(3:3:end)=phaseCSpotLoadP+1j*phaseCSpotLoadQ;
+ck=Sp3+pLoadABC;
+ck(3*(Balance-1)+1:3*(Balance-1)+3)=0;
+if any(abs(ck)>1e-5)
+    fprintf('反推回的功率不匹配。潮流方程约束不满足。\n');
+end
+output_args=ck;
+end
+