PQ1-(1).zip_matlab电力系统_power system_潮流计算pq_电力系统_电力系统matlab

% function [iiSij]=PQ(Bus_number Branch_number jing_du Bus Branch Bus_type Bus_U Bus_a Bus_PG Bus_PL Bus_QG Bus_QL Branch_I Branch_J Rij Xij Zij B0 B00 RT XT ZT KT T_flag X) function DP=PQ(X) % [Bus_number]=xlsread('IEEE11.xls','A2:A2'); %节点个数 % [Branch_number]=xlsread('IEEE11.xls','B2:B2'); %支路个数 % jing_du=xlsread('IEEE11.xls','B4:B4'); %精度要求 % [Bus]=xlsread('IEEE11.xls','D3:J100'); %节点信息 % [Branch]=xlsread('IEEE11.xls','K3:T100'); %支路信息 % [Bus_type]=zeros(Bus_number,1); %节点类型 % D=length(X); % for i=1:D-1 % if round(X(1,2))==5 % DP=100; % else Bus_number=11; Branch_number=10; jing_du=0.00001; Bus_U=zeros(Bus_number,1); %节点电压 Bus_Ub=ones(Bus_number,1); Bus_a=zeros(Bus_number,1);%节点电压相角 Bus_PG=zeros(Bus_number,1); %发电机有功 Bus_PL=zeros(Bus_number,1);%负荷有功 Bus_QG=zeros(Branch_number,1);%发电机无功 Bus_QL=zeros(Branch_number,1);%负荷无功 Branch_I=zeros(Branch_number,1); %支路起始节点 Branch_J=zeros(Branch_number,1); %支路终止节点 Rij=zeros(Branch_number,1); Xij=zeros(Branch_number,1); Zij=j*Xij; %支路阻抗 Y=zeros(Bus_number); %节点导纳矩阵 G=zeros(Bus_number); B=zeros(Bus_number); B0=zeros(Branch_number,1); B00=zeros(Branch_number,1); RT=zeros(Branch_number,1); XT=zeros(Branch_number,1); ZT=j*XT; KT=zeros(Branch_number,1); Bus =[ 3.0000 1.0000 0 35 0 12 0 1.0000 1.0000 0 0 2.28 0 1.49 1.0000 1.0000 0 0 1.51 0 0.316 1.0000 1.0000 0 0 1.51 0 0.316 1.0000 1.0000 0 0 1.51 0 0.316 1.0000 1.0000 0 0 0.888 0 0.55 1.0000 1.0000 0 0 2.738 0 1.395 1.0000 1.0000 0 0 2.91 0 1.185 1.0000 1.0000 0 0 1.51 0 0.816 1.0000 1.0000 0 0 1.956 0 0.993 1.0000 1.0000 0 0 2.295 0 1.422]; Branch =[ 1.0000 2.0000 0.00018 0.00240 0 0 0 0 0 0 1.0000 8.0000 0.00044 0.00330 0 0 0 0 0 0 2.0000 3.0000 0.00009 0.00055 0 0 0 0 0 0 3.0000 4.0000 0.00051 0.00320 0 0 0 0 0 0 4.0000 5.0000 0.00033 0.00220 0 0 0 0 0 0 5.0000 6.0000 0.00008 0.00065 0 0 0 0 0 0 5.0000 7.0000 0.00014 0.00120 0 0 0 0 0 0 8.0000 9.0000 0.00006 0.00005 0 0 0 0 0 0 9.0000 10.0000 0.00006 0.00060 0 0 0 0 0 0 10.0000 11.0000 0.00043 0.00130 0 0 0 0 0 0]; Bus_type=Bus(:,1); Bus_U=Bus(:,2); Bus_a=Bus(:,3); Bus_PG=Bus(:,4); Bus_PL=Bus(:,5); Bus_QG=Bus(:,6); Bus_QL=Bus(:,7); Branch_I=Branch(:,1); Branch_J=Branch(:,2); Rij=Branch(:,3); Xij=Branch(:,4); Zij=j*Xij; B0=Branch(:,5); B00=Branch(:,10); RT=Branch(:,6); XT=Branch(:,7); ZT=j*XT; KT=Branch(:,8); T_flag=Branch(:,9); B0=zeros(Branch_number,1); B00=zeros(Branch_number,1); X(1,2)=round(X(1,2)); if X(1,2)<=1 X(1,2)=1; end DG=X(1,2); Bus_PG(DG)=X(1,1); %——————————————计算导纳矩阵————————————% for m=1:Bus_number %（1） for n=1:Branch_number%(2) if(Branch_J(n)~=0&Branch_I(n)==m&Branch_J(n)~=m) if KT(n)==0 %支路中没有变压器 %（4） Y(Branch_I(n),Branch_J(n))=-1/Zij(n); Y(Branch_J(n),Branch_I(n))=Y(Branch_I(n),Branch_J(n)); else %支路中有变压器 %（5） Y(Branch_I(n),Branch_J(n))=-1/(KT(n)*ZT(n)); Y(Branch_J(n),Branch_I(n))=Y(Branch_I(n),Branch_J(n)); end end if (Branch_J(n)~=0&(Branch_I(n)==m|Branch_J(n)==m))%计算自导纳 if KT(n)==0%支路中没有变压器(7) Y(m,m)=Y(m,m)+1/Zij(n); end if KT(n)~=0%支路中有变压器 if T_flag(n)==0%变压器所在线路的起始节点在低压侧，终止节点是高压 if Branch_I(n)==m%(9) Y(m,m)=Y(m,m)+1/(KT(n)*ZT(n))+(1-KT(n))/(KT(n)*KT(n)*ZT(n))%起始节点在低压侧，所以该节点的π型等值电路接地导纳为Yπ2 end if Branch_J(n)==m %(10) Y(m,m)=Y(m,m)+1/(KT(n)*ZT(n))+(KT(n)-1)/(KT(n)*ZT(n)); %终止节点在高压侧，所以该节点π型等值电路接地导纳为Yπ1 end end if T_flag(n)==1%变压器所在线路的起始节点在高压侧，终止节点在低压侧，一下计算过程刚好与上计算过程相反(11) if Branch_I(n)==m%(12) Y(m,m)=Y(m,m)+1/(KT(n)*ZT(n))+(KT(n)-1)/(KT(n)*ZT(n)); end if Branch_J(n)==m Y(m,m)=Y(m,m)+1/(KT(n)*ZT(n))+(1-KT(n))/(KT(n)*KT(n)*ZT(n)); end end end end if(Branch_I(n)==m&Branch_J(n)==0) Y(m,m)=Y(m,m)+j*B00(n); end end end % %——————————————导纳矩阵存档——————————————% % fid=fopen('E:\xwy1（k=5）\Y.txt','wt'); % for m=1:Bus_number % for n=1:Bus_number % fprintf(fid, '%3.4f+i*(%3.4f)\t',real(Y(m,n)),imag(Y(m,n))); % end % fprintf(fid, '\n'); % end % fclose(fid); %——————————————求解B‘——————————————% PHwz=find(Bus_type(:,1)==3);%查找平衡节点所在位置 B1=imag(Y);%取Y矩阵的虚部 B1(:,PHwz)=[];%将平衡节点所在的列置空 B1(PHwz,:)=[];%将平衡节点所在的行置空 %——————————————求解B“——————————————% PHPVwz=find(Bus_type(:,1)==3|Bus_type(:,1)==2);%查找平衡节点和PV节点所在位置 B2=imag(Y); B2(:,PHPVwz)=[]; B2(PHPVwz,:)=[]; %——————————————求解ΔP数据准备阶段——————————————% G=real(Y); B=imag(Y); PQPVwz=find(Bus_type(:,1)==1|Bus_type(:,1)==2);%获取PQ和PV节点所在位置 PQPVnum=0; for m=1:Bus_number%统计PQ和PV节点的个数，确定ΔP的阶数 if (Bus_type(m,1)==1|Bus_type(m,1)==2) PQPVnum=PQPVnum+1; end end iPi=zeros(PQPVnum,1);%ΔP的阶数同PQPV节点个数之和 U_PQPV=Bus_U; U_PQPV(PHwz)=[]; %——————————————求解ΔQ数据准备阶段——————————————% G=real(Y); B=imag(Y); PQwz=find(Bus_type(:,1)==1); PQnum=0; for m=1:Bus_number%统计PQ节点的个数，确定ΔQ的阶数 if Bus_type(m,1)==1 PQnum=PQnum+1; end end iQi=zeros(PQnum,1); U_PQ=Bus_U; U_PQ(PHPVwz)=[]; %——————————————迭代过程——————————————% k=0; Kp=1; Kq=1; b=1; c=1; num1=1;%画图用 while(((Kp~=0)||(Kq~=0))&&k<=7) Kp=1; Kq=1; %——————————————求解ΔP——————————————% for m=1:PQPVnum%累加公式中i是除了平衡节点以外的所有节点 sum1=0;%i不同，ΔPi不同，算完一个ΔPi后，算下一个时要将累加寄存器中的数据置零 for n=1:Bus_number%累加公式中j是所有节点 sum1=sum1+Bus_U(n)*(G(PQPVwz(m),n)*cos(Bus_a(PQPVwz(m))-Bus_a(n))+B(PQPVwz(m),n)*sin(Bus_a(PQPVwz(m))-Bus_a(n))); %ΔPi的计算公式 end iPi(m)=(Bus_PG(PQPVwz(m))-Bus_PL(PQPVwz(m)))-sum1*Bus_U(PQPVwz(m)); end iPi_max=max(abs(iPi)); % e1=max(b); if iPi_max<=jing_du Kp=0; if Kq==0 break else %——————————————求解ΔQ——————————————% for m=1:PQnum sum2=0; for n=1:Bus_numb