facat@lab.com
parent
fbffecb341
commit
4a6fa42f5b
|
|
@ -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
|
||||
|
|
@ -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
|
||||
142
OPF.m
142
OPF.m
|
|
@ -5,6 +5,148 @@ 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);
|
||||
% %全部转换为负荷电流
|
||||
% CurpA=conj((iterPhaseASpotLoadP+1j*iterPhaseASpotLoadQ)./VoltpA);
|
||||
% CurpB=conj((iterPhaseBSpotLoadP+1j*iterPhaseBSpotLoadQ)./VoltpB);
|
||||
% CurpC=conj((iterPhaseCSpotLoadP+1j*iterPhaseCSpotLoadQ)./VoltpC);
|
||||
% %转换为序电流
|
||||
% f012=Tp2f*conj([CurpA';CurpB';CurpC']);
|
||||
% %把三序电流分离出来
|
||||
% If0=conj(f012(1,:)');
|
||||
% If1=conj(f012(2,:)');
|
||||
% If2=conj(f012(3,:)');
|
||||
% %试着算一下正序电流
|
||||
% fsY11*V1;
|
||||
% %形成负荷序电流的测量值
|
||||
% mIf0=If0;
|
||||
% mIf1=If1;
|
||||
% mIf1(3)=-mIf1(2);
|
||||
% mIf2=If2;
|
||||
% %计算
|
||||
% fsY11=fsY11+sparse(Balance,Balance,ones(length(Balance),1),busNum,bus
|
||||
% Num);%这里要置0，置1，否则是奇异的
|
||||
%%做最小二乘法
|
||||
[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);
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
|
@ -0,0 +1,18 @@
|
|||
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)
|
||||
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));
|
||||
dP(Balance)=0;
|
||||
dQ(QGi)=0;
|
||||
end
|
||||
|
|
@ -0,0 +1,25 @@
|
|||
function [ output_args ] = checkSSatisfied(Balance,phaseABCY,VoltpABC,phaseASpotLoadP,phaseBSpotLoadP,phaseCSpotLoadP,phaseASpotLoadQ,phaseBSpotLoadQ,phaseCSpotLoadQ )
|
||||
%%利用Fortiscue方法计算得到了相电压和相角后，反推负荷功率，检查是否和给定的功率一致。
|
||||
busNum=size(phaseABCY,1)/3;
|
||||
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);
|
||||
pLoadABC=sparse(length(phaseASpotLoadP)*3,1);
|
||||
pLoadABC(1:3:end)=phaseASpotLoadP+1j*phaseASpotLoadQ;
|
||||
pLoadABC(2:3:end)=phaseBSpotLoadP+1j*phaseBSpotLoadQ;
|
||||
pLoadABC(3:3:end)=phaseCSpotLoadP+1j*phaseCSpotLoadQ;
|
||||
ck=Sp3+pLoadABC;
|
||||
ck(3*(Balance-1)+1:3*(Balance-1)+3)=0;
|
||||
if any(abs(ck)>1e-5)
|
||||
fprintf('反推回的功率不匹配。潮流方程约束不满足。\n');
|
||||
end
|
||||
output_args=ck;
|
||||
end
|
||||
|
||||
Loading…
Reference in New Issue