内容发布更新时间 : 2024/10/25 6:25:14星期一 下面是文章的全部内容请认真阅读。

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); Tlta_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