commit f55585e444b98b579d6b63e2eaf282564a92111c Author: facat Date: Thu Apr 17 00:32:05 2014 +0800 收敛了。但是精度还有点问题，需要调一下。感觉是有一步没更新。 Signed-off-by: facat diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..b6392e9 --- /dev/null +++ b/.gitignore @@ -0,0 +1,3 @@ +*.asv +/feeder13 +*.mexw32 \ No newline at end of file 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 +%% ƽڵDz +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/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..ca0f61c --- /dev/null +++ b/Unbalance.m @@ -0,0 +1,27 @@ +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,fsY1,Vf1) +%real(diag(Vmf1.*exp(1j*Vaf1))*(fsY1*(Vmf1.*exp(1j*Vaf1)))); +% Y=abs(fsY1); +% [r,c,Yangle]=find(fsY1); +% Yangle=angle(Yangle); +% Vmf1=abs(Vf1); +% Vaf1=angle(Vf1); +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)); +%orderedPDPG+real(Vf2.*conj(If2)+Vf0.*conj(If0))+real(diag(Vf1)*conj(fsY1*Vf1)); +%dP=diag_Vmf1_YdotCosVmf1+orderedPDPG; +%dQ=diag_Vmf1_YdotSinVmf1+orderedQDQG; +dP(Balance)=0; +dQ(QGi)=0; +end \ No newline at end of file diff --git a/YMatrix.m b/YMatrix.m new file mode 100644 index 0000000..d59cc9a --- /dev/null +++ b/YMatrix.m @@ -0,0 +1,52 @@ +function [Y,Yangle,r,c,GB]=YMatrix(Busnum,lineI,lineJ,lineR,lineX,lineB2,groundbranchI,groundbranchB,transI,transJ,transR,transX,transK) +%% · +% t1=-lineR./(lineR.^2+lineX.^2);%絼 +% t2=lineX./(lineR.^2+lineX.^2);% +% G = +sparse(lineI,lineJ,t1,Busnum,Busnum) + sparse(lineJ,lineI,t1,Busnum,Busnum); +% B = sparse(lineI,lineJ,t2,Busnum,Busnum)+sparse(lineJ,lineI,t2,Busnum,Busnum); +% G = G - sparse(1:Busnum,1:Busnum,sum(G,2)'); +% B = B - sparse(1:Busnum,1:Busnum,sum(B,2)'); +t1=lineR+1j*lineX; +t2=1./t1; +% realT2=real(t2); +% imagT2=imag(t2); + +GB=sparse(lineI,lineJ,-t2,Busnum,Busnum) + sparse(lineJ,lineI,-t2,Busnum,Busnum); +GB=GB-sparse(1:Busnum,1:Busnum,sum(GB,2)); +% G = +sparse(lineI,lineJ,-realT2,Busnum,Busnum) + sparse(lineJ,lineI,-realT2,Busnum,Busnum); +% B = sparse(lineI,lineJ,-imagT2,Busnum,Busnum)+sparse(lineJ,lineI,-imagT2,Busnum,Busnum); +% G = G - sparse(1:Busnum,1:Busnum,sum(G,2)'); +% B = B - sparse(1:Busnum,1:Busnum,sum(B,2)'); +% G=real(GB); +% B=imag(GB); +%% ӵص +% t3=sparse(lineI,lineI,lineB2,Busnum,Busnum)+sparse(lineJ,lineJ,lineB2,Busnum,Busnum);%ӵ֧· +%B=B+t3; +t1=1j*lineB2; +GB=GB+sparse(lineI,lineI,t1,Busnum,Busnum)+sparse(lineJ,lineJ,t1,Busnum,Busnum); +% B=imag(GB); +%% ӵ֧· +if isempty(groundbranchI)==0 %ǿ +% B=B+sparse(groundbranchI,groundbranchI,groundbranchB,Busnum,Busnum); + + GB=GB+sparse(groundbranchI,groundbranchI,1j*groundbranchB,Busnum,Busnum); +% B=imag(GB); +end +%% ѹ +if isempty(transI)==0 %ǿ +% t1 = -transR./(transR.^2+transX.^2); +% t2 = transX./(transR.^2+transX.^2); + t1=1./(transR+1j*transX); +% G = G+sparse(transI,transJ,t1./transK,Busnum,Busnum)+sparse(transJ,transI,t1./transK,Busnum,Busnum)-sparse(transI,transI,t1./transK./transK,Busnum,Busnum)-sparse(transJ,transJ,t1,Busnum,Busnum); +% B = B+sparse(transI,transJ,t2./transK,Busnum,Busnum)+sparse(transJ,transI,t2./transK,Busnum,Busnum)-sparse(transI,transI,t2./transK./transK,Busnum,Busnum)-sparse(transJ,transJ,t2,Busnum,Busnum); + GB=GB+sparse(transI,transJ,-t1./transK,Busnum,Busnum)+sparse(transJ,transI,-t1./transK,Busnum,Busnum)+sparse(transI,transI,t1./transK./transK,Busnum,Busnum)+sparse(transJ,transJ,t1,Busnum,Busnum); +% G=real(GB); +% B=imag(GB); +end +%GB=G+1j*B; +Y = abs(GB); +[r,c,Yangle] = find(GB); +Yangle=angle(Yangle); +%Yangle=angle(GB); +%Angle = angle(GB(GB~=0)); +end \ No newline at end of file diff --git a/Yf2p.m b/Yf2p.m new file mode 100644 index 0000000..f197983 --- /dev/null +++ b/Yf2p.m @@ -0,0 +1,23 @@ +function [ phaseABCY ] = Yf2p( fs0,fs1,fs2 ) +[r,c]=find(fs1); +a=exp(1j*2*pi/3); +Tp2f=1/3*[1 1 1; + 1 a a^2; + 1 a^2 a]; +Tf2p=inv(Tp2f); +% +% Tf2p=[1 1 1; +% 1 a^2 a; +% 1 a a^2]; +busNum=size(fs0,1); +phaseABCY=sparse(3*busNum,3*busNum); +for I=1:length(r) + v0=fs0(r(I),c(I)); + v1=fs1(r(I),c(I)); + v2=fs2(r(I),c(I)); + pabc=Tf2p*diag([v0,v1,v2])*Tp2f; + [pr,pc,pv]=find(pabc); + phaseABCY=phaseABCY+sparse(pr+3*(r(I)-1),pc+3*(c(I)-1),pv,3*busNum,3*busNum); +end +end + diff --git a/dataRead.m b/dataRead.m new file mode 100644 index 0000000..02994f2 --- /dev/null +++ b/dataRead.m @@ -0,0 +1,43 @@ +function [ fsY0, fsY1, fsY2,phaseASpotLoadP,phaseBSpotLoadP,phaseCSpotLoadP ... + phaseASpotLoadQ,phaseBSpotLoadQ,phaseCSpotLoadQ,nodeNum,Balance,phaseABCY] = dataRead(lineZ,dataFile ) +data=dlmread(dataFile); +zeroEntries=find(data(:,1)==0); +lines=data(zeroEntries(1)+1:zeroEntries(2)-1,:); +[setIJ,nodeNum]=numberNode(lines); +Balance=data(1,1); +Balance=nodeNum(setIJ==Balance); +%%three-phase +phaseABCY=sparse(3*length(nodeNum),3*length(nodeNum)); +%% 601 begin +[fs30,fs31,fs32,retphaseABCY]=lineWithConfig(setIJ,nodeNum,lineZ,lines,601); +%phaseABCY ĵɾ3n x 3n ά +fsY0=fs30; +fsY1=fs31; +fsY2=fs32; +phaseABCY=phaseABCY+retphaseABCY; +% 601 end +%% 602 begin +[fs30,fs31,fs32,retphaseABCY]=lineWithConfig(setIJ,nodeNum,lineZ,lines,602); +fsY0=fsY0+fs30; +fsY1=fsY1+fs31; +fsY2=fsY2+fs32; +phaseABCY=phaseABCY+retphaseABCY; +% 602 end +%% spot load +busNum=length(nodeNum); +spotloads=data(zeroEntries(3)+1:zeroEntries(4)-1,:); +spotloads(:,2:end)=spotloads(:,2:end)/1000; +phaseASpotLoadP=sparse(busNum,1); +phaseBSpotLoadP=sparse(busNum,1); +phaseCSpotLoadP=sparse(busNum,1); +phaseASpotLoadQ=sparse(length(nodeNum),1); +phaseBSpotLoadQ=sparse(length(nodeNum),1); +phaseCSpotLoadQ=sparse(length(nodeNum),1); +phaseASpotLoadP( ismember(setIJ,spotloads(:,1) ) )=spotloads(:,2); +phaseBSpotLoadP( ismember(setIJ,spotloads(:,1) ) )=spotloads(:,4); +phaseCSpotLoadP( ismember(setIJ,spotloads(:,1) ) )=spotloads(:,6); +phaseASpotLoadQ( ismember(setIJ,spotloads(:,1) ) )=spotloads(:,3); +phaseBSpotLoadQ( ismember(setIJ,spotloads(:,1) ) )=spotloads(:,5); +phaseCSpotLoadQ( ismember(setIJ,spotloads(:,1) ) )=spotloads(:,7); +end + diff --git a/lineWithConfig.m b/lineWithConfig.m new file mode 100644 index 0000000..49d78b6 --- /dev/null +++ b/lineWithConfig.m @@ -0,0 +1,89 @@ +function [ fs30,fs31,fs32,phaseABCY ] = lineWithConfig(setIJ,nodeNum,lineZ,lines,config ) +phase3Entry=find(lines(:,1)==config); +phase3Line=lines(phase3Entry,:); +phase3Line(:,4)=phase3Line(:,4)/1000; +entry=find(lineZ(:,1)==config); +phase3R=lineZ(entry+1:entry+3,:); +phase3X=lineZ(entry+5:entry+7,:); +phase3B2=lineZ(entry+9:entry+11,:); +phase3Y=1./(phase3R+1j*phase3X);%ԭʼർ +%phase3Y=1./(1j*phase3X);%ԭʼർ +%öԽǺͷǶԽ +phase3Y(1,3)=phase3Y(1,2); +phase3Y(3,1)=phase3Y(2,1); +phase3Y(3,2)=phase3Y(3,1); +phase3Y(2,3)=phase3Y(1,3); +phase3Y(2,2)=phase3Y(1,1); +phase3Y(3,3)=phase3Y(1,1); +%ֵ +baseHighU=4.16; +baseLowU=0.48; +baseS=1; +baseY=baseS/(baseHighU^2); +phase3Y=phase3Y/baseY; +%Fortescue任 +%ֱд +a=exp(1j*2*pi/3); +Tp2f=1/3*[1 1 1; + 1 a a^2; + 1 a^2 a]; +Tf2p=inv(Tp2f); +% +% Tf2p=[1 1 1; +% 1 a^2 a; +% 1 a a^2]; +fs3Y=Tp2f*phase3Y*Tf2p; +Zl=phase3Y(1,1); +Zm=phase3Y(2,1); +fs3Y=3/3*diag([Zl+2*Zm,Zl-Zm,Zl-Zm]);%ﲻ1/3IEEEϵĹʽˡ +fs3Y(abs(fs3Y)<1e-5)=0; +fs3Y=sparse(fs3Y); +lineNodeI=zeros(length(phase3Line(:,2) ) ,1); +lineNodeJ=zeros(length(phase3Line(:,2) ) ,1); +for I=1:length(lineNodeI) + lineNodeI(I)=nodeNum(setIJ==phase3Line(I,2) ); + lineNodeJ(I)=nodeNum(setIJ==phase3Line(I,3) ); +end +busNum=length(nodeNum); +%ർɾ +phaseABCY=sparse(busNum*3,busNum*3); +for I=1:length(phase3Line(:,4)) + [phase3Yr,phase3Yc,phase3Yv]=find(phase3Y./phase3Line(I,4)); + offsetI=phase3Yr+3*(lineNodeI(I)-1); + offsetJ=phase3Yc+3*(lineNodeJ(I)-1); + phaseABCY=phaseABCY+sparse(offsetI,offsetJ,-phase3Yv,busNum*3,busNum*3); + phaseABCY=phaseABCY+sparse(offsetJ,offsetI,-phase3Yv,busNum*3,busNum*3); + %1 + diagY=sum(phaseABCY(:,1:3:end),2); + offsetI=1:busNum*3; + offsetJ=zeros(busNum*3,1); + for J=0:busNum*3-1 + offsetJ(J+1)=floor(J/3)*3+1; + end + phaseABCY=phaseABCY-sparse(offsetI,offsetJ,diagY,busNum*3,busNum*3); + %2 + diagY=sum(phaseABCY(:,2:3:end),2); + offsetJ=zeros(busNum*3,1); + for J=0:busNum*3-1 + offsetJ(J+1)=floor(J/3)*3+2; + end + phaseABCY=phaseABCY-sparse(offsetI,offsetJ,diagY,busNum*3,busNum*3); + %3 + diagY=sum(phaseABCY(:,3:3:end),2); + offsetJ=zeros(busNum*3,1); + for J=0:busNum*3-1 + offsetJ(J+1)=floor(J/3)*3+3; + end + phaseABCY=phaseABCY-sparse(offsetI,offsetJ,diagY,busNum*3,busNum*3); +end +fs30=sparse(lineNodeI,lineNodeJ,-fs3Y(1,1)./phase3Line(:,4),busNum,busNum); +fs30=fs30+sparse(lineNodeJ,lineNodeI,-fs3Y(1,1)./phase3Line(:,4),busNum,busNum); +fs30=fs30-sparse(1:busNum,1:busNum,sum(fs30,2)); +fs31=sparse(lineNodeI,lineNodeJ,-fs3Y(2,2)./phase3Line(:,4),busNum,busNum); +fs31=fs31+sparse(lineNodeJ,lineNodeI,-fs3Y(2,2)./phase3Line(:,4),busNum,busNum); +fs31=fs31-sparse(1:busNum,1:busNum,sum(fs31,2)); +fs32=sparse(lineNodeI,lineNodeJ,-fs3Y(3,3)./phase3Line(:,4),busNum,busNum); +fs32=fs32+sparse(lineNodeJ,lineNodeI,-fs3Y(3,3)./phase3Line(:,4),busNum,busNum); +fs32=fs32-sparse(1:busNum,1:busNum,sum(fs32,2)); +end + diff --git a/numberNode.m b/numberNode.m new file mode 100644 index 0000000..2cc497c --- /dev/null +++ b/numberNode.m @@ -0,0 +1,8 @@ +function [ setIJ,nodeNum ] = numberNode( lines ) +%% ڵ +setI=lines(:,2); +setJ=lines(:,3); +setIJ=union(setI,setJ); +nodeNum=(1:length(setIJ))'; +end + diff --git a/readLineZ.m b/readLineZ.m new file mode 100644 index 0000000..5f71d9c --- /dev/null +++ b/readLineZ.m @@ -0,0 +1,4 @@ +function [ data ] = readLineZ( lineParaFile ) +data=dlmread(lineParaFile); +end + diff --git a/run.m b/run.m new file mode 100644 index 0000000..d03be14 --- /dev/null +++ b/run.m @@ -0,0 +1,265 @@ +clc +clear +lineZ=readLineZ('.\feeder13\lineParameter.txt'); +[ fsY0, fsY1, fsY2,phaseASpotLoadP,phaseBSpotLoadP,phaseCSpotLoadP ... + phaseASpotLoadQ,phaseBSpotLoadQ,phaseCSpotLoadQ,nodeNum,Balance,phaseABCY]=dataRead(lineZ,'.\feeder13\data1.txt'); +%fsY1(1,1)=fsY1(1,1)+1j*1e-10; +a=exp(1j*2*pi/3); +Tp2f=1/3*[1 1 1; + 1 a a^2; + 1 a^2 a]; +Tf2p=inv(Tp2f); +fsY1amp=abs(fsY1); +[r,c,fsY1ang]=find(fsY1); +fsY1ang=angle(fsY1ang); +% fsY2amp=abs(fsY2); +% fsY0ang=angle(fsY0); +% fsY1ang=angle(fsY1); +% fsY2ang=angle(fsY2); +Pabc=phaseASpotLoadP+phaseBSpotLoadP+phaseCSpotLoadP; +Qabc=phaseASpotLoadQ+phaseBSpotLoadQ+phaseCSpotLoadQ; +busNum=length(phaseASpotLoadP); +%ѹֵ +Vmf1=sparse(ones(busNum,1)); +Vmf2=sparse(ones(busNum,1)); +Vmf0=sparse(busNum,1); +Vaf1=sparse(busNum,1); +Vaf2=sparse(busNum,1); +Vaf0=sparse(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; +%%Fתᵽ󣬿 +phaseABCYalt=Yf2p( fsY0,fsY1,fsY2 ); +%% + Vf1=-ones(busNum,1); + fsY11=fsY1; + fsY00=fsY0; + fsY22=fsY2; +while(k<=kmax && maxD> EPS) + k=k+1; + % %㸺 + % Vf2=Vmf2.*exp(1j*Vaf2) ; + % If2=fsY2*(Vf2); + % % + % Vf0=Vmf0.*exp(1j*Vaf0) ; + % If0=fsY0*(Vmf0.*exp(1j*Vaf0) ); + %תΪѹ + VoltpABC=Tf2p*conj([ (Vmf0.*exp(1j*Vaf0))'; (Vmf1.*exp(1j*Vaf1))'; (Vmf2.*exp(1j*Vaf2))']); + 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,:)'); + %ƽڵ㸺ѹΪ0 + fsY2(Balance,:)=0; + fsY2(:,Balance)=0; + fsY2=fsY2+sparse(Balance,Balance,ones(length(Balance),1),busNum,busNum); + If2(Balance)=0; + Vf2=fsY2\If2; + Vmf2=abs(Vf2); + Vaf2=angle(Vf2); + %If1(Balance)=1; + %fsY1(Balance,:)=0; + %fsY1(:,Balance)=0; + %fsY1=fsY1+sparse(Balance,Balance,ones(length(Balance),1),busNum,busNum); + %Vf1=fsY1\If1;%һ£Ȼá + %ƽڵѹΪ0 + fsY0(Balance,:)=0; + fsY0(:,Balance)=0; + fsY0=fsY0+sparse(Balance,Balance,ones(length(Balance),1),busNum,busNum); + If0(Balance)=0; + Vf0=fsY0\If0; + Vmf0=abs(Vf0); + Vaf0=angle(Vf0); + + [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,fsY11,Vf1);%ƽ + maxD=max([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\n',k,full(maxD)); + %% + %תΪѹ + VoltpABC=Tf2p*conj([ (Vmf0.*exp(1j*Vaf0))'; (Vmf1.*exp(1j*Vaf1))'; (Vmf2.*exp(1j*Vaf2))']); + 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,:)'); + Vf0=fsY0\If0; + Vmf0=abs(Vf0); + Vaf0=angle(Vf0); + Vf2=fsY2\If2; + Vmf2=abs(Vf2); + Vaf2=angle(Vf2); + Vf1=Vmf1.*exp(1j*Vaf1); + %% +% VoltpABC=Tf2p*conj([ (Vmf0.*exp(1j*Vaf0))'; (Vmf1.*exp(1j*Vaf1))'; (Vmf2.*exp(1j*Vaf2))']); +% 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)=-sum(If0); +% If2(Balance)=-sum(If2); +% If1(Balance)=-sum(If1); + %% +end +%% +Vtest=[1.02*exp(1j*0);1.01*exp(1j*0.21);1.00*exp(1j*-0.13)];%ѹ +%Itest=[0.12*exp(1j*0.15);0.09*exp(1j*1.21);0.81*exp(1j*-0.43)];% +Ytest=diag([0.75*exp(1j*1.7);33*exp(1j*2.1);12*exp(1j*-3)]); +Itest=Ytest*Vtest; +Stest=Vtest.*conj(Itest);%๦ +sum(Stest)/3; +sum((Tp2f*Vtest).*conj(Tp2f*Itest)); +YYtest=Ytest; +Ytest(2,:)=0; +Ytest(:,2)=0; +Ytest(2,2)=1; +tttt=Itest(2); +Itest(2)=1.01*exp(1j*0.21); +tV=Ytest\Itest; +Itest(2)=tttt; +YYtest*tV; +%% +(Vf0.*conj(fsY00*Vf0)+Vf1.*conj(fsY11*Vf1)+Vf2.*conj(fsY22*Vf2))*3; +(Vf1.*conj(fsY11*Vf1))/3; +conj(Tf2p*[If0(2);If1(2);If2(2)]).*(Tf2p*[Vf0(2);Vf1(2);Vf2(2)]); +%% +If1=fsY11*Vf1; +If1(Balance)=Vf1(Balance); +fsY1(Balance,:)=0; +fsY1(:,Balance)=0; +% fsY1=fsY1+sparse(Balance,Balance,ones(length(Balance),1),busNum,busNum); +fsY1(3,3)=1; +%VV1=fsY1(1:2,1:2)\If1(1:2);%һ£Ȼá +VV1=fsY1\If1;%һ£Ȼá +%Vf1=Vmf1.*exp(1j*Vaf1); +If0(Balance)=-sum(If0); +If2(Balance)=-sum(If2); +If1(Balance)=-sum(If1); +fsY11*VV1; +%% +IpABC=Tf2p*conj([If0';If1';If2']); +%תѹ +VoltpABC=Tf2p*conj([ (Vmf0.*exp(1j*Vaf0))'; (Vmf1.*exp(1j*Vaf1))'; (Vmf2.*exp(1j*Vaf2))']); +VoltpABC(:,1); +Tp2f*VoltpABC(:,1); +Tf2p*Tp2f*VoltpABC(:,1); +disp([' A B C']) +abs(VoltpABC') +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); +I2inj=IpABC(:,2); +aa=phaseABCY(1:3,1:3)*[VoltpA(1);VoltpB(1);VoltpC(1)]; +a1=Tf2p*diag([fsY0(2,1),fsY1(2,1),fsY2(2,1)])*[Vf0(1);Vf1(1);Vf2(1)]; +a2=Tf2p*diag([fsY0(2,2),fsY1(2,2),fsY2(2,2)])*[Vf0(2);Vf1(2);Vf2(2)]; +a3=Tf2p*diag([fsY0(2,3),fsY1(2,3),fsY2(2,3)])*[Vf0(3);Vf1(3);Vf2(3)]; +bb=0; +a1=diag([fsY0(2,2)])*[Vf0(2)]; +a2=diag([fsY1(2,2)])*[Vf1(2)]; +a3=diag([fsY2(2,2)])*[Vf2(2)]; +bb=bb+1*[a1;a2;a3]; +a1=diag([fsY0(2,1)])*[Vf0(1)]; +a2=diag([fsY1(2,1)])*[Vf1(1)]; +a3=diag([fsY2(2,1)])*[Vf2(1)]; +bb=bb+1*[a1;a2;a3]; +a1=diag([fsY0(2,3)])*[Vf0(3)]; +a2=diag([fsY1(2,3)])*[Vf1(3)]; +a3=diag([fsY2(2,3)])*[Vf2(3)]; +bb=bb+1*[a1;a2;a3]; +b=0; +aa=phaseABCY(4:6,4:6)*[VoltpA(2);VoltpB(2);VoltpC(2)]; +b=b+aa; +aa=phaseABCY(4:6,1:3)*[VoltpA(1);VoltpB(1);VoltpC(1)]; +b=b+aa; +aa=phaseABCY(4:6,7:9)*[VoltpA(3);VoltpB(3);VoltpC(3)]; +b=b+aa; +aa=phaseABCY(4:6,:)*[VoltpA(1);VoltpB(1);VoltpC(1);VoltpA(2);VoltpB(2);VoltpC(2);VoltpA(3);VoltpB(3);VoltpC(3)]; +conj(aa).*[VoltpA(2);VoltpB(2);VoltpC(2)]; +conj(a1+a2+a3).*[VoltpA(2);VoltpB(2);VoltpC(2)]; +aa=[15.3502321158364 - 26.9339348391836i,30.6661479563297 - 84.3093119155620i,30.6661479563297 - 84.3093119155620i;30.6661479563297 - 84.3093119155620i,15.3502321158364 - 26.9339348391836i,30.6661479563297 - 84.3093119155620i;30.6661479563297 - 84.3093119155620i,30.6661479563297 - 84.3093119155620i,15.3502321158364 - 26.9339348391836i;]; +aa*Vp3(1:3); +fs=Tp2f*aa*Tf2p; +phaseABCY(7,:)=0; +phaseABCY(:,7)=0; +phaseABCY(8,:)=0; +phaseABCY(:,8)=0; +phaseABCY(9,:)=0; +phaseABCY(:,9)=0; +phaseABCY(7,7)=1; +phaseABCY(8,8)=1; +phaseABCY(9,9)=1; +Ip3(7)=1; +Ip3(8)=-0.5 - 0.866025403784439i; +Ip3(9)=-0.5 + 0.866025403784439i; +inv(phaseABCY)*Ip3; +IppABC=sparse(9,1); +IppABC(1:3:end)=IpABC(1,:); +IppABC(2:3:end)=IpABC(2,:); +IppABC(3:3:end)=IpABC(3,:); +IppABC(7)=1; +IppABC(8)=-0.5 - 0.866025403784439i; +IppABC(9)=-0.5 + 0.866025403784439i; +inv(phaseABCY)*IppABC; +aa*Tf2p*[Vf0(1);Vf1(1);Vf2(1)]-Tf2p*[If0(1);If1(1);If2(1)]; +diag([fsY0(1),fsY1(1),fsY2(1)])*[Vf0(1);Vf1(1);Vf2(1)]-[If0(1);If1(1);If2(1)]; +Tp2f*diag([fsY0(1),fsY1(1),fsY2(1)])*Tf2p; +altY012=sparse(3,3); +altY012(1,1)=fsY0(1,1); +altY012(2,2)=fsY1(1,1); +altY012(3,3)=fsY2(1,1); +V11I=altY012*[Vf0(1);Vf1(1);Vf2(1)]; +V11S=[VoltpA(1);VoltpB(1);VoltpC(1)].*conj(Tf2p*altY012*Tp2f*[VoltpA(1);VoltpB(1);VoltpC(1)]); +altY012(1,1)=fsY0(1,2); +altY012(2,2)=fsY1(1,2); +altY012(3,3)=fsY2(1,2); +V12I=altY012*[Vf0(2);Vf1(2);Vf2(2)]; +V12S=[VoltpA(1);VoltpB(1);VoltpC(1)].*conj(Tf2p*altY012*Tp2f*[VoltpA(1);VoltpB(1);VoltpC(1)]); +altY012(1,1)=fsY0(1,3); +altY012(2,2)=fsY1(1,3); +altY012(3,3)=fsY2(1,3); +V13I=altY012*[Vf0(3);Vf1(3);Vf2(3)]; +V12S+V11S; +[VoltpA(1);VoltpB(1);VoltpC(1)].*(Tf2p*(V11I+V12I)); +Tf2p*V11I; +altY012(1,1)=fsY0(1,1); +altY012(2,2)=fsY1(1,1); +altY012(3,3)=fsY2(1,1); +Tf2p*altY012*Tp2f*[VoltpA(1);VoltpB(1);VoltpC(1)]; +Tf2p*[Vf0(1);Vf1(1);Vf2(1)]; \ No newline at end of file