243 lines
7.7 KiB
Matlab
243 lines
7.7 KiB
Matlab
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');
|
|
%% 潮流计算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);
|
|
%给序电压赋初值
|
|
Vmf1=sparse(ones(busNum,1));
|
|
Vaf1=sparse(zeros(busNum,1));
|
|
%先求解正序的
|
|
PQi=nodeNum;
|
|
PG=sparse(busNum,1);
|
|
QG=sparse(busNum,1);
|
|
QGi=[Balance];
|
|
PD=Pabc/3;
|
|
QD=Qabc/3;
|
|
Loadi=find(PD~=0);
|
|
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);
|
|
If0=sparse(busNum,1);
|
|
%准备序矩阵
|
|
%平衡节点置0置1
|
|
fsY2(Balance,:)=0;
|
|
fsY2(:,Balance)=0;
|
|
fsY2=fsY2+sparse(Balance,Balance,ones(length(Balance),1),busNum,busNum);
|
|
%平衡节点置0置1
|
|
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
|
|
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;
|
|
%把补偿电容看作负荷
|
|
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));
|
|
%先不要电容
|
|
% 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);%不平衡量
|
|
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((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('迭代时间%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))))
|
|
%% 潮流计算end
|
|
|
|
|
|
busNum=length(nodeNum);
|
|
Busnum=busNum;
|
|
PQi=setxor(nodeNum,Balance);
|
|
QGi=[Balance];
|
|
% 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;
|
|
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;
|
|
wPD=ones(length(Loadi),1);
|
|
wQD=ones(length(Loadi),1);
|
|
PD0=PD3P;
|
|
QD0=QD3P;
|
|
%%
|
|
RestraintCount=size(Loadi,1)*2; %约束条件数,放开所有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);%与学姐一致
|
|
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;
|
|
%% 加误差
|
|
Volt=Vp3m;
|
|
UAngel=Vp3a;
|
|
while(abs(Gap)>Precision)
|
|
if KK>kmax
|
|
break;
|
|
end
|
|
Init_u=Gap/2/RestraintCount*CenterA;
|
|
AngleIJMat=0;
|
|
%% 开始计算OPF
|
|
%% 形成等式约束的雅克比
|
|
deltH=func_deltH(Busnum,Volt,Y,UAngel,r,c,Yangle,Loadi);
|
|
%% 形成不等式约束的雅克比
|
|
deltG=func_deltG(Busnum,Loadi);
|
|
%%
|
|
L_1Z=diag(Init_Z./Init_L);
|
|
U_1W=diag(Init_W./Init_U);
|
|
%% 形成海森阵
|
|
deltdeltF=func_deltdeltF(wPD,wQD,ContrlCount);
|
|
% deltdeltF=0;
|
|
%% 形成ddHy
|
|
ddh=func_ddh(Volt,Init_Y,Busnum,Y,UAngel,r,c,Yangle,Loadi,ContrlCount);
|
|
%% 开始构建ddg
|
|
% ddg=func_ddg(PGi,PVi,Busnum,RestraintCount,Loadi,PD,QD);
|
|
ddg=0;
|
|
%% 开始构建deltF
|
|
deltF=func_deltF(wPD,wQD,PD0,PD3P,QD3P,QD0,Busnum,Loadi);
|
|
% deltF=0;
|
|
%% 形成方程矩阵
|
|
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);
|
|
%% 开始解方程
|
|
% fprintf('迭代次数 %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);
|
|
%%取各分量
|
|
[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('迭代次数%d\n',KK);
|
|
toc
|
|
|