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