clc
clear
lineZ=readLineZ('D:\share\feeder123\lineParameter.txt');
[ fsY0, fsY1, fsY2,phaseASpotLoadP,phaseBSpotLoadP,phaseCSpotLoadP ...
    phaseASpotLoadQ,phaseBSpotLoadQ,phaseCSpotLoadQ,setIJ,nodeNum,Balance,phaseABCY]=dataRead(lineZ,'D:\share\feeder123\data.txt');
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));
% Vf1=Vmf1.*exp(1j*Vaf1);
%Vaf1=sparse(zeros(busNum,1));
%先求解正序的
PQi=nodeNum;
PG=sparse(busNum,1);
QG=sparse(busNum,1);
QGi=[Balance];
PD=Pabc/3;
QD=Qabc/3;
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);
Vf0=sparse(busNum,1);
If0=sparse(busNum,1);
%准备序矩阵
fsY2(Balance,:)=0;
fsY2(:,Balance)=0;
fsY2=fsY2+sparse(Balance,Balance,ones(length(Balance),1),busNum,busNum);
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
while(k<=kmax && maxD> EPS)
    k=k+1;
    [dP, dQ, YdotSinVolt, YdotCosVolt, diag_Volt_YdotSin, diag_Volt_YdotCos]=Unbalance(Balance,busNum, ...
        PQi,PG,QG,QGi,PD,QD,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((phaseASpotLoadP+1j*phaseASpotLoadQ)./VoltpA);
    VoltpB=conj(VoltpABC(2,:)');
    CurpB=-conj((phaseBSpotLoadP+1j*phaseBSpotLoadQ)./VoltpB);
    VoltpC=conj(VoltpABC(3,:)');
    CurpC=-conj((phaseCSpotLoadP+1j*phaseCSpotLoadQ)./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))))
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);