ipmthreephasese/OPF.m

tic
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');
%% <EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>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);
%<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>
%ƽ<EFBFBD><EFBFBD><EFBFBD>ڵ<EFBFBD><EFBFBD><EFBFBD>0<EFBFBD><EFBFBD>1
fsY2(Balance,:)=0;
fsY2(:,Balance)=0;
fsY2=fsY2+sparse(Balance,Balance,ones(length(Balance),1),busNum,busNum);
%ƽ<EFBFBD><EFBFBD><EFBFBD>ڵ<EFBFBD><EFBFBD><EFBFBD>0<EFBFBD><EFBFBD>1
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));
    %<EFBFBD>Ȳ<EFBFBD>Ҫ<EFBFBD><EFBFBD><EFBFBD><EFBFBD>
%     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);
    [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))))
%% <EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>end


busNum=length(nodeNum);
Busnum=busNum;
PQi=setxor(nodeNum,Balance);
QGi=[Balance];
% fprintf('<EFBFBD><EFBFBD>ʼţ<EFBFBD>ٷ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>\n');
[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;
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;
PGA=sum(PD3P(1:3:end))*1.03;
PGB=sum(PD3P(2:3:end))*1.03;
PGC=sum(PD3P(3:3:end))*1.03;
QGA=sum(QD3P(1:3:end))*1.03;
QGB=sum(QD3P(2:3:end))*1.03;
QGC=sum(QD3P(3:3:end))*1.03;
PG3P=sparse( (Balance-1)*3+1:(Balance-1)*3+3,1,[PGA,PGB,PGC],Busnum*3,1);
QG3P=sparse( (Balance-1)*3+1:(Balance-1)*3+3,1,[QGA,QGB,QGC],Busnum*3,1);
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;
Loadi=PQi3P;
wVolt=ones(busNum*3,1);
wPD=ones(length(Loadi),1);
wQD=ones(length(Loadi),1);
mVolt=Vp3m;
PD0=PD3P;
QD0=QD3P;
%%
RestraintCount=size(Loadi,1)*2; %Լ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>,<EFBFBD>ſ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>QD
Init_Z=sparse(ones(1,RestraintCount));
Init_W=sparse(-1*ones(1,RestraintCount));
Init_L=1*sparse(ones(1,RestraintCount));
Init_U=1*sparse(ones(1,RestraintCount));
Init_Y=sparse(1,2*Busnum*3);%<EFBFBD><EFBFBD>ѧ<EFBFBD><EFBFBD>һ<EFBFBD><EFBFBD>
Gap=(Init_L*Init_Z'-Init_U*Init_W');
KK=0;
ContrlCount=size(Loadi,1)*2+Busnum*6;
kmax=30;
%%
Precision=1e-5;
CenterA=0.1;
%% <EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
Volt=Vp3m;
UAngel=Vp3a;
while(abs(Gap)>Precision)
    if KK>kmax
        break;
    end
    Init_u=Gap/2/RestraintCount*CenterA;
    AngleIJMat=0;
    %% <EFBFBD><EFBFBD>ʼ<EFBFBD><EFBFBD><EFBFBD><EFBFBD>OPF
    %% <EFBFBD>γɵ<EFBFBD>ʽԼ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ſ˱<EFBFBD>
    deltH=func_deltH(Busnum,Volt,Y,UAngel,r,c,Yangle,Loadi);
    %% <EFBFBD>γɲ<EFBFBD><EFBFBD><EFBFBD>ʽԼ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ſ˱<EFBFBD>
    deltG=func_deltG(Busnum,Loadi);
    %%
    L_1Z=diag(Init_Z./Init_L);
    U_1W=diag(Init_W./Init_U);
    %% <EFBFBD>γɺ<EFBFBD>ɭ<EFBFBD><EFBFBD>
    deltdeltF=func_deltdeltF(wVolt,wPD,wQD,ContrlCount);
%     deltdeltF=0;
    %% <EFBFBD>γ<EFBFBD>ddHy
    ddh=func_ddh(Volt,Init_Y,Busnum,Y,UAngel,r,c,Yangle,Loadi,ContrlCount);
    %%  <EFBFBD><EFBFBD>ʼ<EFBFBD><EFBFBD><EFBFBD><EFBFBD>ddg
%     ddg=func_ddg(PGi,PVi,Busnum,RestraintCount,Loadi,PD,QD);
    ddg=0;
    %%  <EFBFBD><EFBFBD>ʼ<EFBFBD><EFBFBD><EFBFBD><EFBFBD>deltF
    deltF=func_deltF(wVolt,wPD,wQD,PD0,PD3P,QD3P,QD0,Volt,mVolt,Busnum,Loadi);
%     deltF=0;
    %% <EFBFBD>γɷ<EFBFBD><EFBFBD>̾<EFBFBD><EFBFBD><EFBFBD>
    Luu=Init_U'.*Init_W'+Init_u*ones(RestraintCount,1);
    Lul=Init_L'.*Init_Z'-Init_u*ones(RestraintCount,1);
    Mat_G=FormG(PD3P,QD3P,Loadi);
    Mat_H=FormH(Busnum,Volt,PG3P,PD3P,QG3P,QD3P,Y,UAngel,r,c,Yangle);
    Ly=Mat_H;
    Lz=FormLz(Mat_G,Init_L,Loadi);
    Lw=FormLw(Mat_G,Init_U,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>
%     fprintf('<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD> %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);
    %%ȡ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
    [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,PD3P,QD3P]=Modification(Init_Z,Init_L,Init_W,Init_U,Init_Y,deltZ,deltL,deltW,deltU,deltX,deltY,PG3P,QG3P,Volt,UAngel,ContrlCount,Balance,Busnum,PD3P,QD3P,Loadi);
    Gap=(Init_L*Init_Z'-Init_U*Init_W');
    fprintf('%f\n',full(Gap));
    KK=KK+1;
end
fprintf('<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>%d\n',KK);
toc