distributionnetwork-power2current/run.m

clc
clear
lineZ=readLineZ('feeder13\lineParameter.txt');
[ fsY0, fsY1, fsY2,phaseASpotLoadP,phaseBSpotLoadP,phaseCSpotLoadP ...
    phaseASpotLoadQ,phaseBSpotLoadQ,phaseCSpotLoadQ,setIJ,nodeNum,Balance,phaseABCY ...
    cap]=dataRead(lineZ,'feeder13\data1.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);
%<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ѹ<EFBFBD><EFBFBD><EFBFBD><EFBFBD>ֵ
Vmf1=sparse(ones(busNum,1));
Vaf1=sparse(zeros(busNum,1));
%<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
PQi=nodeNum;
PG=sparse(busNum,1);
QG=sparse(busNum,1);
QGi=[Balance];
PD=Pabc/3;
QD=Qabc/3;
Loadi=find(PD~=0);
maxD=100000;% <EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ƽ<EFBFBD><EFBFBD><EFBFBD><EFBFBD>
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);
%׼<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
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<EFBFBD>ֽ<EFBFBD>
[fsY0L,fsY0U,fsY0P,fsY0Q,fsY0R]=lu(fsY0);
[fsY2L,fsY2U,fsY2P,fsY2Q,fsY2R]=lu(fsY2);
%<EFBFBD><EFBFBD><EFBFBD><EFBFBD>ʼ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
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;
    %<EFBFBD>Ѳ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ݿ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
    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));
    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);
    [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);%<EFBFBD><EFBFBD>ƽ<EFBFBD><EFBFBD><EFBFBD><EFBFBD>
    maxD=max(abs([dP;dQ;]));
    jaco=Jacobi(Balance,busNum,QGi,Vmf1,YdotSinVolt,YdotCosVolt,diag_Volt_YdotSin,diag_Volt_YdotCos);%<EFBFBD>ſ˱Ⱦ<EFBFBD><EFBFBD><EFBFBD>
    [dV, dVangle]=Solv(busNum,jaco,dP,dQ);%<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
    [Vmf1, Vaf1]=Modify(Vmf1,Vaf1,dV,dVangle,1);
    fprintf('<EFBFBD><EFBFBD> %d <EFBFBD>ε<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD> <EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ƽ<EFBFBD><EFBFBD><EFBFBD><EFBFBD>Ϊ %f\n',k,full(maxD));
    %ת<EFBFBD><EFBFBD>Ϊ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ѹ
    VoltpABC=Tp2f\conj([ Vf0'; (Vmf1.*exp(1j*Vaf1))'; Vf2']);%<EFBFBD><EFBFBD>Tp2f\ <EFBFBD><EFBFBD><EFBFBD><EFBFBD>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('<EFBFBD><EFBFBD><EFBFBD><EFBFBD>ʱ<EFBFBD><EFBFBD>%f\n',toc);
    
    %
end
FortiscueToc=toc;
fprintf('Fortiscue<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ʱ<EFBFBD><EFBFBD> %f\n',FortiscueToc);
Vf1=Vmf1.*exp(1j*Vaf1);
%%
(Vf0.*conj(fsY00*Vf0)+Vf1.*conj(fsY11*Vf1)+Vf2.*conj(fsY22*Vf2))*3;%<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ݵĹ<EFBFBD><EFBFBD><EFBFBD>
conj(Tf2p*[If0(2);If1(2);If2(2)]).*(Tf2p*[Vf0(2);Vf1(2);Vf2(2)]);
IpABC=Tf2p*conj([If0';If1';If2']);
%ת<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ѹ
VoltpABC=Tf2p*conj([ Vf0'; Vf1'; Vf2']);
disp(['                  A                        B                          C'])
full(abs(VoltpABC'))
fprintf('<EFBFBD>ڵ<EFBFBD><EFBFBD>Ŷ<EFBFBD>Ӧ\n');
disp([setIJ,nodeNum ])
%%<EFBFBD><EFBFBD><EFBFBD>鷴<EFBFBD>ƻ<EFBFBD>ȥ<EFBFBD>Ĺ<EFBFBD><EFBFBD><EFBFBD><EFBFBD>Ƿ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
ub=checkSSatisfied(Balance,phaseABCY,VoltpABC, ...
    phaseASpotLoadP,phaseBSpotLoadP,phaseCSpotLoadP, ...
    phaseASpotLoadQ,phaseBSpotLoadQ,phaseCSpotLoadQ );
fprintf('<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ƽ<EFBFBD><EFBFBD><EFBFBD><EFBFBD>Ϊ%f\n\n',full(max(abs(ub))))
fprintf('<EFBFBD><EFBFBD>ʼţ<EFBFBD>ٷ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>\n');
%% <EFBFBD><EFBFBD>ţ<EFBFBD>ٷ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>begin
% [r,c,GB]=find(phaseABCY);
% Y=abs(phaseABCY);
% Yangle=angle(GB);
% Vp3=sparse(ones(busNum*3,1));%<EFBFBD><EFBFBD><EFBFBD><EFBFBD>ѹ<EFBFBD><EFBFBD><EFBFBD><EFBFBD>ֵ
% 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);%<EFBFBD>ſ˱Ⱦ<EFBFBD><EFBFBD><EFBFBD>
%     [dV, dVangle]=Solv(busNum*3,jaco,dP,dQ);%<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
%     [Vp3m, Vp3a]=Modify(Vp3m,Vp3a,dV,dVangle,1);
%     fprintf('<EFBFBD><EFBFBD> %d <EFBFBD>ε<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD> <EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ƽ<EFBFBD><EFBFBD><EFBFBD><EFBFBD>Ϊ %f\n',k,full(maxD));
%     fprintf('<EFBFBD><EFBFBD><EFBFBD><EFBFBD>ʱ<EFBFBD><EFBFBD>%f\n',toc);
% end
% NewtonToc=toc;
% fprintf('ţ<EFBFBD>ٷ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ʱ<EFBFBD><EFBFBD> %f\n',NewtonToc);
% fprintf('<EFBFBD><EFBFBD><EFBFBD>ٱ<EFBFBD>Ϊ%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));
%% <EFBFBD><EFBFBD>ţ<EFBFBD>ٷ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>end
%% <EFBFBD><EFBFBD>ʼ<EFBFBD><EFBFBD><EFBFBD><EFBFBD>״̬<EFBFBD><EFBFBD><EFBFBD><EFBFBD>
%״̬<EFBFBD><EFBFBD>
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);
%<EFBFBD><EFBFBD>ʼ<EFBFBD><EFBFBD>
%״̬<EFBFBD><EFBFBD>Ϊ 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);%<EFBFBD><EFBFBD>ʽԼ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
Gap=(Init_L'*Init_Z-Init_U'*Init_W);
PG=sparse(busNum,1);
PG(Balance)=0.1105;
QD=sparse(busNum,1);
QG(Balance)=0.0984;
while Gap>1e-5 && KK<60
    KK=KK+1;
    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);
    %% <EFBFBD><EFBFBD>ʼ<EFBFBD>ⷽ<EFBFBD><EFBFBD>
    plotGap(KK)=Gap;
    fprintf('<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD> %d Gap %f\n',KK,plotGap(KK));
    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);
end