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

1、PQ分解法计算大电网 潮流程序 -CAL-FENGHAI-(2020YEAR-YICAI) JINGBIAN function PQ %用卩(1分解法计算大电网潮流 % %bus数组1.节点编号2.节点电压3.节点电压角度4.注入有功5.注入无功6.节点类型(1PQ 2PV 3 平衡) %line数组1.始端节点编号2.末端节点编号3.电阻4电抗5电导G 6电纳B 7.变比 %打开数据文件 clear clc bus=load(,!); line=load(); linenum(:/l/2)=line(:,1,2); (nb/v=size(bus); (nl/=size(line); % no

2、denum=(l:nb) bus(: J); %带入子函数数据处理 busJinePCbnPVSWodenum =changel_busline( busjine );%对盯点晅新编匕 Y = admittance(busjine/l);%生成节点导纳矩阵 Yl= admittancefbusjine );%生成化简条件3的矩阵B1 Y2=admittance(bus/line/3 );%生成化简条件 3 的矩阵 B2 %临时添加的测试数据 % nPQ二4; nPV=0;nSW=l;nb=5; % Y= +5i+5i+-5+15i % +5i 10+30i 0 + % +5i 10+30i +

3、 0 % + 0 + 0 % -5+15i+ 0 0 % %Yl=+5i+5i + -5+15i % +5i 10+30i 0 + % +5i 10+30i + 0 % + 0 + 0 % -5+15i+ 0 0 % % bus=l 101 %21 0 1 %31 01 %41 0 1 %5 0 0 0 3; % % line=5 20 0 %23 10 -30 0 0 0 %34 0 0 0 %41 0 0 0 %12 5 0 0 0 %13 5 0 0 0 %15 5 15 0 0 0； % bus_PV0=bus(nPQ+l):end/2),;%*ones(l/nPV+nSW); bus

4、_U=ones(l/nPQ) bus_PVO;%电JE幅值 bus_e=zeros(nb,l); %电压角度 delta_P=zeros(nPQ+nPV/l); delta_Q=zeros(nPCbl)； % delta_e=zerosnb-l/l); %delta_U=zeros( nPQl); c=O；KP=l;KQ=l；%KP KQ用来判断有功、无功是否收敛 G=real(Y);B=imag(Y);B10=imag(Yl);B20=imag(Y2);% BO 是进化简三后的讥导纳矩阵虚部 %形成解耦潮流的系数矩阵B1和B2 Bl=B10(l:nb-l,l:nb-l); B2=B2O(l:

5、nPCbl：nPQ); while c delta_e=-(UP#Bl)error_P; bus_e=bus_e+delta_e;O); c=c+l; KQ 二 1; else KP=O; if KCT二0 else break end end for ii=l:nPQ delta_Q(ii)=bus(ii/5); for jj=l: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(

6、bus_U(l:( nbnPV nSW); erro r_Q 二 U Qd e 11 a_Q; if max(abs(error_Q) delta_U=-B2error_Q; bus_U=bus_U+delta_U;zeros( nPV+nSW)，l); c=c+l; KP=1; else KQ=O; if KP JO else break end end end %至此得到收敛的节点电压值 % % %对计算结果进行数据处理 %将节点结果用原节点编号农示 bus_Ue=zeros(nb/3)； bus_Ue(:/l/2,3)=nodenum(:/2) bus_U bus_e/pi*180; f

7、or ii=l:nb for jj=ii+l:nb if bus_Ue(ii/l)bus_Ue(jj/l) t=bus_Ue(ii,:); bus_Ue(ii,:)二 bus_UG(jj,:); bus_Ue(jj，:)=t; end end end %r_U是收敛的电压衣达成复数的形式 r_U=zeros(nbzl); for k=l:nb r_U(k)=bus_U(k)#(cos(bus_e(k)+li#sin(bus_e(k); end %计算平衡节点功率 sw_s=o； SW_S 二 SW_S+r_U(nb 广 conj(Y( nb,:)广 c on j(r_U); %计算各支路功率S

8、ij line_S=zeros( nb,nb);line_SO 二 zeros( nb，nb); for ii=l:nb for jj=l:nb line_S(ii/jj)=r_U(ii)*(conj(r_U(ii)*conj(Yii/ii)+(conj(r_U(ii)-conj(r_U(jj)*conj(Y(ii/jj); end end % %把线路结果还原成原节点编号对应的结果 for ii=l:nb for jj=l:nb line_S0(nodenum(ii/2)/nodenum(jj/2)=line_S(ii/jj); end end line_P=real(line_SO);li

9、ne_Q=imag(line_SO); %计算各支路损耗 delta_S=zeros(nhl); for k=l:nl a=linenum(k,l);b 二 line nu m(k,2); delta_S(k)=line_SO(a/b)+line_S(b,a); end %计算网络总损耗 SO=sum(delta_S); % %将计算结果输入指定文件 fid=fopen(C:UserslrDesktopmatlab 练习训练题大电网潮流计算twt1)； fprintf(fid； ij点号t节点电Jk幅值t节点电压角度n; for k=l:nb fprintf(fid/!%dt%ft%fn,/k

10、,bus_Ue(kzl)/bus_Ue(k/2); end fprintf(fid|支路首端t支路末端t支路仃功t支路无功t支路损耗W); for k=l:nl fprintf(fid/%dtt%dtt%At%ft%fn,Jinenum(k/l)/linenum(k/2)/line_P(linenum(k/l) ,linenum(k/2)/line_Q (linGnum(kJ)inenum(k,2),delta_S(k); end fprintf(fid；平衡节点功率=%fn：SW_S); fprintf(fid；网络总损耗=%fn : SO); fclose(fid); end functi

11、on busJinePCbnPVSWodenum =changel_busline( busjine ) %此函数用来对原始输入节点、线路数据进行重新编号 % %bus数组1.节点编号2.节点电压3.节点电压角度4.注入有功5.注入无功6.节点类型(1PQ2PV3 平衡) %line数组1 始端节点编号2末端节点编号3电阻4电抗5电导G 6电纳B 7变比 nb/sizbus); nl/=size(line); % nodenum=(l:nb) bus nPQ=O; nPV=0;nSW=0;%PQ=;PQ=zeros(nb/6);PZ=zeros(nb/6);SW=zeros(nb/6);%PQ

12、 PV 平衡节丄 的个数 for k=l:nb switch bus(k,6) case 1 nPQ=nPQ+l; PQ(nPQ:)=bus(匕:); case 2 nPV=nPV+l; PV(nPV/:)=bus(k/:); case 3 nSW=nSW+l; SW(nSW,:)=bus(k/:); otherwise dispC节点数据类型出错！ 7； end end %生成重新编号后的节点数据矩阵 bus=PQ;PV;SW; nodenum=(l:nb)bus(:zl);%第列为新的节点编号，第二列为对应的旧节点编号 bus(:,l)=(l:n; %至此实现广节点数据的重新编号 % %对

13、线路数据重新编号 % nodenum=(l:nb)1 bus(:zl);%第列为新的节点编号，第二列为对应的I口节点编号 for ii=l:nl rl/v=find(nodenum(:/2)=line(ii/l); line(ii/l)=nodenum(rl/l); r2/v=find(nodenum(:/2)=line(ii/2); Iine(ii/2)=nodenum(r2/l); end end function Y = admittance(busjin已c ) %此函数用来形成节点导纳矩阵 % %bus数组1 节点编号2节点电压3节点电压角度4注入有功5注入无功6节点类型(1PQ 2

14、PV 3 平衡) %line数组1.始端节点编号2.末端节点编号3.电阻4电抗5电导6电纳B/2 7.变比 %c是用来控制形成节点导纳矩阵的方式的，c=l形成般的节点导纳矩阵，可以用来确定B2, c=2,形成化简条件3的节点导纳矩阵，确定B1 nbz=size(bus); (nl/AZ=size(line); Y=zeros (n b,nb); zt=zeros( nl,l);yt二zeros( nl，l);ym 二 ZGros(nl,l);l=Zros(nlJ);J=zeros( nl);K二zeros(nl,l); switch c case 1 for k=l:nl zt(k)=line

15、(kz3)+i 产 line(匕 4); yt(k)=l/zt(k); ym(k)=line(k/5)+li*line(k/6); l(k)=line(k/l);J(k)=line(kz2);K(k)=line(k/7); end case 2 for k=l:nl yt(k)=l/(line(k/3)+li*line(k/4); ym=zeros(nbl); l(k)二 line(kJ);J(k)=line(k，2); if line(k,7)性 1 K(k)=l; else K(k)=line(k/7); end end case 3 for k=l:nl yt(k)=l/(li*line(k/4); ym(k)=line(kz5)+li*line(k/6); l(k)=line(k/l);J(k)=line(kz2);K(k)=line(k/7); end %针对不同线路设置节点导纳值 for k=l:nl if (K(k)=l) Y(J(k),J(k)=Y(J(k),J(k)+yt(k)+ym(k); Y(l(k),J(k)=Y(l(k),J(k)-yt(k); Y(J(k),l(k)=Y(l(k),J(k); end if