From f55585e444b98b579d6b63e2eaf282564a92111c Mon Sep 17 00:00:00 2001 From: facat Date: Thu, 17 Apr 2014 00:32:05 +0800 Subject: [PATCH] =?UTF-8?q?=E6=94=B6=E6=95=9B=E4=BA=86=E3=80=82=E4=BD=86?= =?UTF-8?q?=E6=98=AF=E7=B2=BE=E5=BA=A6=E8=BF=98=E6=9C=89=E7=82=B9=E9=97=AE?= =?UTF-8?q?=E9=A2=98=EF=BC=8C=E9=9C=80=E8=A6=81=E8=B0=83=E4=B8=80=E4=B8=8B?= =?UTF-8?q?=E3=80=82=E6=84=9F=E8=A7=89=E6=98=AF=E6=9C=89=E4=B8=80=E6=AD=A5?= =?UTF-8?q?=E6=B2=A1=E6=9B=B4=E6=96=B0=E3=80=82?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Signed-off-by: facat --- .gitignore | 3 + Jacobi.m | 25 +++++ Modify.m | 5 + Solv.m | 16 +++ Unbalance.m | 27 +++++ YMatrix.m | 52 ++++++++++ Yf2p.m | 23 ++++ dataRead.m | 43 ++++++++ lineWithConfig.m | 89 ++++++++++++++++ numberNode.m | 8 ++ readLineZ.m | 4 + run.m | 265 +++++++++++++++++++++++++++++++++++++++++++++++ 12 files changed, 560 insertions(+) create mode 100644 .gitignore create mode 100644 Jacobi.m create mode 100644 Modify.m create mode 100644 Solv.m create mode 100644 Unbalance.m create mode 100644 YMatrix.m create mode 100644 Yf2p.m create mode 100644 dataRead.m create mode 100644 lineWithConfig.m create mode 100644 numberNode.m create mode 100644 readLineZ.m create mode 100644 run.m 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 +%% 平衡节点相角不变 +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/3，IEEE上的公式错了。 +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