文档库 最新最全的文档下载
当前位置:文档库 › PQ分解法计算大电网潮流程序

PQ分解法计算大电网潮流程序

PQ分解法计算大电网潮流程序
PQ分解法计算大电网潮流程序

function PQ

%用PQ分解法计算大电网潮流

% %bus数组 1.节点编号 2.节点电压 3.节点电压角度 4.注入有功 5.注入无功 6.节点类型(1PQ 2PV 3平衡)

%line数组 1.始端节点编号 2.末端节点编号 3.电阻 4电抗 5电导G 6电纳B 7.变比

%打开数据文件

clear

clc

bus=load('ieee14bus.txt');

line=load('ieee14line.txt');

linenum(:,[1,2])=line(:,[1,2]);

[nb,~]=size(bus);

[nl,~]=size(line);

% nodenum=[(1:nb)' bus(:,1)];

%带入子函数数据处理

[bus,line,nPQ,nPV,nSW,nodenum] =change1_busline( bus,line );%对节点重新编号Y = admittance(bus,line,1 );%生成节点导纳矩阵

Y1= admittance(bus,line,2 );%生成化简条件3的矩阵B1

Y2=admittance(bus,line,3 );%生成化简条件3的矩阵B2

%-----------------------------------------------------

% %临时添加的测试数据

% nPQ=4; nPV=0;nSW=1;nb=5;

% Y=[10.834-32.5i -1.667+5i -1.667+5i -2.5+7.5i -5+15i

% -1.667+5i 12.917-38.75i -10+30i 0 -1.25+3.75i

% -1.667+5i -10+30i 12.917-38.75i -1.25+3.75i 0

% -2.5+7.5i 0 -1.25+3.75i 3.75-11.25i 0

% -5+15i -1.25+3.75i 0 0 6.25-18.75i];

%

% Y1=[10.834-32.5i -1.667+5i -1.667+5i -2.5+7.5i -5+15i

% -1.667+5i 12.917-38.75i -10+30i 0 -1.25+3.75i

% -1.667+5i -10+30i 12.917-38.75i -1.25+3.75i 0

% -2.5+7.5i 0 -1.25+3.75i 3.75-11.25i 0

% -5+15i -1.25+3.75i 0 0 6.25-18.75i];

%

% bus=[1 1 0 0.2 0.2 1

% 2 1 0 -0.45 -0.15 1

% 3 1 0 -0.4 -0.05 1

% 4 1 0 -0.6 -0.1 1

% 5 1.06 0 0 0 3];

%

% line=[5 2 1.25 -3.75 0 0 0

% 2 3 10 -30 0 0 0

% 3 4 1.25 -3.75 0 0 0

% 4 1 2.5 -7.5 0 0 0

% 1 2 1.667 -5 0 0 0

% 1 3 1.667 -5 0 0 0

% 1 5 5 -15 0 0 0];

%-------------------------------------------------------

bus_PV0=bus((nPQ+1):end,2)';%1.05*ones(1,nPV+nSW);

bus_U=[ones(1,nPQ) bus_PV0]';%电压幅值

bus_e=zeros(nb,1); %电压角度

delta_P=zeros(nPQ+nPV,1);

delta_Q=zeros(nPQ,1);

% delta_e=zeros(nb-1,1);

%delta_U=zeros(nPQ,1);

c=0;KP=1;KQ=1;%KP KQ用来判断有功、无功是否收敛

G=real(Y);B=imag(Y);B10=imag(Y1);B20=imag(Y2);%矩阵B0是进行化简三后的节点导纳矩阵虚部

%形成解耦潮流的系数矩阵B1和B2

B1=B10(1:nb-1,1:nb-1);

B2=B20(1:nPQ,1:nPQ);

while c<80

%求解P Q的不平衡量

for ii=1:nPQ+nPV

delta_P(ii)=bus(ii,4);

for jj=1:nb

delta_P(ii)=delta_P(ii)-bus_U(ii)*bus_U(jj)*(G(ii,jj)*cos(bus_e(ii)-bus _e(jj))+B(ii,jj)*sin(bus_e(ii)-bus_e(jj)));

end

end

UP=diag(bus_U(1:(nb-1)));%U矩阵利用各节点电压形成对角阵,来计算修正方程,对角线上的元素与bus_U列元素一一对应

error_P=UP\delta_P;

if max(abs(error_P))>0.00001

delta_e=-(UP*B1)\error_P;

bus_e=bus_e+[delta_e;0];

c=c+1;

KQ=1;

else KP=0;

if KQ~=0

else

break

end

end

for ii=1:nPQ

delta_Q(ii)=bus(ii,5);

for jj=1:nb

delta_Q(ii)=delta_Q(ii)-bus_U(ii)*bus_U(jj)*(G(ii,jj)*sin(bus_e(ii)-bus _e(jj))-B(ii,jj)*cos(bus_e(ii)-bus_e(jj)));

end

end

UQ=diag(bus_U(1:(nb-nPV-nSW)));

error_Q=UQ\delta_Q;

if max(abs(error_Q))>0.00001

delta_U=-B2\error_Q;

bus_U=bus_U+[delta_U;zeros((nPV+nSW),1)];

c=c+1;

KP=1;

else KQ=0;

if KP~=0

else

break

end

end

end

%至此得到收敛的节点电压值

%----------------------------------------------

%--------------------------------------------

% 对计算结果进行数据处理

%将节点结果用原节点编号表示

bus_Ue=zeros(nb,3);

bus_Ue(:,[1,2,3])=[nodenum(:,2) bus_U bus_e/pi*180];

for ii=1:nb

for jj=ii+1:nb

if bus_Ue(ii,1)>bus_Ue(jj,1)

t=bus_Ue(ii,:);

bus_Ue(ii,:)=bus_Ue(jj,:);

bus_Ue(jj,:)=t;

end

end

end

%r_U是收敛的电压表达成复数的形式

r_U=zeros(nb,1);

for k=1:nb

r_U(k)=bus_U(k)*(cos(bus_e(k))+1i*sin(bus_e(k)));

end

%计算平衡节点功率

SW_S=0;

SW_S=SW_S+r_U(nb)*conj(Y(nb,:))*conj(r_U);

%计算各支路功率Sij

line_S=zeros(nb,nb);line_S0=zeros(nb,nb);

for ii=1:nb

for jj=1:nb

line_S(ii,jj)=r_U(ii)*(conj(r_U(ii))*conj(Y(ii,ii))+(conj(r_U(ii))-conj (r_U(jj)))*conj(Y(ii,jj)));

end

end

%------------------------------------------

%把线路结果还原成原节点编号对应的结果

for ii=1:nb

for jj=1:nb

line_S0(nodenum(ii,2),nodenum(jj,2))=line_S(ii,jj);

end

end

line_P=real(line_S0);line_Q=imag(line_S0);

%计算各支路损耗

delta_S=zeros(nl,1);

for k=1:nl

a=linenum(k,1);b=linenum(k,2);

delta_S(k)=line_S0(a,b)+line_S(b,a);

end

%计算网络总损耗

S0=sum(delta_S);

%-------------------------------------------------

%将计算结果输入指定文件

fid=fopen('C:\Users\lr\Desktop\matlab练习\训练题\大电网潮流计算

\ieee14_out.txt','wt');

fprintf(fid,'节点号\t节点电压幅值\t节点电压角度\n');

for k=1:nb

fprintf(fid,'%d\t%f\t%f\n',k,bus_Ue(k,1),bus_Ue(k,2));

end

fprintf(fid,'支路首端\t支路末端\t支路有功\t支路无功\t支路损耗\n');

for k=1:nl

fprintf(fid,'%d\t\t%d\t\t%f\t%f\t%f\n',linenum(k,1),linenum(k,2),line_P (linenum(k,1),linenum(k,2)),line_Q(linenum(k,1),linenum(k,2)),delta_S(k ));

end

fprintf(fid,'平衡节点功率=%f\n',SW_S);

fprintf(fid,'网络总损耗=%f\n',S0);

fclose(fid);

end

function [bus,line,nPQ,nPV,nSW,nodenum] =change1_busline( bus,line )

%此函数用来对原始输入节点、线路数据进行重新编号

% %bus数组 1.节点编号 2.节点电压 3.节点电压角度 4.注入有功 5.注入无功 6.节点类型(1PQ 2PV 3平衡)

%line数组 1.始端节点编号 2.末端节点编号 3.电阻 4电抗 5电导G 6电纳B 7.变比

[nb,~]=size(bus);

[nl,~]=size(line);

% nodenum=[(1:nb)' bus(:,1)];

nPQ=0;

nPV=0;nSW=0;%PQ=[];PQ=zeros(nb,6);PV=zeros(nb,6);SW=zeros(nb,6);%PQ PV 平衡节点的个数

for k=1:nb

switch bus(k,6)

case 1

nPQ=nPQ+1;

PQ(nPQ,:)=bus(k,:);

case 2

nPV=nPV+1;

PV(nPV,:)=bus(k,:);

case 3

nSW=nSW+1;

SW(nSW,:)=bus(k,:);

otherwise

disp('节点数据类型出错!');

end

end

%生成重新编号后的节点数据矩阵

bus=[PQ;PV;SW];

nodenum=[(1:nb)' bus(:,1)];%第一列为新的节点编号,第二列为对应的旧节点编号

bus(:,1)=(1:nb)'; %至此实现了节点数据的重新编号

%------------------------------------------------------

%对线路数据重新编号

% nodenum=[(1:nb)' bus(:,1)];%第一列为新的节点编号,第二列为对应的旧节点编号

for ii=1:nl

[r1,~]=find(nodenum(:,2)==line(ii,1));

line(ii,1)=nodenum(r1,1);

[r2,~]=find(nodenum(:,2)==line(ii,2));

line(ii,2)=nodenum(r2,1);

end

end

function Y = admittance(bus,line,c )

%此函数用来形成节点导纳矩阵

% %bus数组 1.节点编号 2.节点电压 3.节点电压角度 4.注入有功 5.注入无功 6.节点类型(1PQ 2PV 3平衡)

%line数组 1.始端节点编号 2.末端节点编号 3.电阻 4电抗 5电导 6电纳B/2 7.变比

%c是用来控制形成节点导纳矩阵的方式的,c=1 形成一般的节点导纳矩阵,可以用来确定B2,c=2,形成化简条件3的节点导纳矩阵,确定B1

[nb,~]=size(bus);

[nl,~]=size(line);

Y=zeros(nb,nb);

zt=zeros(nl,1);yt=zeros(nl,1);ym=zeros(nl,1);I=zeros(nl,1);J=zeros(nl,1 );K=zeros(nl,1);

switch c

case 1

for k=1:nl

zt(k)=line(k,3)+1i*line(k,4);

yt(k)=1/zt(k);

ym(k)=line(k,5)+1i*line(k,6);

I(k)=line(k,1);J(k)=line(k,2);K(k)=line(k,7);

end

case 2

for k=1:nl

yt(k)=1/(line(k,3)+1i*line(k,4));

ym=zeros(nl,1);

I(k)=line(k,1);J(k)=line(k,2);

if line(k,7)~=1

K(k)=1;

else K(k)=line(k,7);

end

end

case 3

for k=1:nl

yt(k)=1/(1i*line(k,4));

ym(k)=line(k,5)+1i*line(k,6);

I(k)=line(k,1);J(k)=line(k,2);K(k)=line(k,7);

end

end

%针对不同线路设置节点导纳值

for k=1:nl

if (K(k)==1)&&(J(k)~=0) %普通线路

Y(I(k),I(k))=Y(I(k),I(k))+yt(k)+ym(k);

Y(J(k),J(k))=Y(J(k),J(k))+yt(k)+ym(k);

Y(I(k),J(k))=Y(I(k),J(k))-yt(k);

Y(J(k),I(k))=Y(I(k),J(k));

end

if (K(k)==1)&&(J(k)==0) %对地支路

Y(I(k),I(k))=Y(I(k),I(k))+ym(k);

end

if (K(k)>1)&&(J(k)~=0) %非标准变比在j侧的变压器支路,zt ym为折算到i侧的值 Y(I(k),I(k))=Y(I(k),I(k))+yt(k)+ym(k);

Y(J(k),J(k))=Y(J(k),J(k))+yt(k)/K(k)^2;

Y(I(k),J(k))=Y(I(k),J(k))-yt(k)/K(k);

Y(J(k),I(k))=Y(I(k),J(k));

end

if (K(k)<1)&&(J(k)~=0)

Y(I(k),I(k))=Y(I(k),I(k))+yt(k)+ym(k);

Y(J(k),J(k))=Y(J(k),J(k))+yt(k)*K(k)^2;

Y(I(k),J(k))=Y(I(k),J(k))-yt(k)*K(k);

Y(J(k),I(k))=Y(I(k),J(k));

end

end

end

相关文档