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

1、function PQ%用 PQ 分解法计算大电网潮流% %bus 数组 1. 节点编号 2. 节点电压 3. 节点电压角度 4. 注入有功 5. 注入无功 6. 节点类型 (1PQ 2PV 3 平衡 )%line 数组 1. 始端节点编号 2. 末端节点编号 3. 电阻 4 电抗 5 电导 G 6 电纳 B 7. 变比%打开数据文件 clearclcbus=load( ieee14bus.txt ); line=load( ieee14line.txt ); linenum(:,1,2)=line(:,1,2);nb,=size(bus);nl,=size(line);% nodenum=(

2、1:nb) bus(:,1);%带入子函数数据处理 bus,line,nPQ,nPV,nSW,nodenumY = admittance(bus,line,1 );Y1= admittance(bus,line,2 );Y2=admittance(bus,line,3 );=change1_busline( bus,line%生成节点导纳矩阵%生成化简条件 3 的矩阵 B1%生成化简条件 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

3、 -5+15i%-1.667+5i12.917-38.75i -10+30i0 -1.25+3.75i%-1.667+5i-10+30i 12.917-38.75i-1.25+3.75i 0%-2.5+7.5i0 -1.25+3.75i 3.75-11.25i 0%-5+15i-1.25+3.75i 0 06.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.75

4、i -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 310 -30 000% 3 41.25 -3.7500 0% 4 12.5 -7.5000% 1 21.667 -5000% 1 31.667 -5000%

5、 1 55 -15 000;%bus_PV0=bus(nPQ+1):end,2); 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);%1.05*ones(1,nPV+nSW);%电压幅值%电压角度B0 是进行化简三后的节点导c=0;KP=1;KQ=1;%KP KQ 用来判断有功、无功是否收敛G=real(Y);B=imag(Y);B10=imag(Y1);B20=im

6、ag(Y2);%矩阵纳矩阵虚部%形成解耦潮流的系数矩阵 B1 和 B2B1=B10(1:nb-1,1:nb-1);B2=B20(1:nPQ,1:nPQ);while c0.00001delta_e=-(UP*B1)error_P;bus_e=bus_e+delta_e;0;c=c+1;KQ=1;else KP=0;if KQ=0elsebreakendendfor ii=1:nPQdelta_Q(ii)=bus(ii,5);for jj=1:nbdelta_Q(ii)=delta_Q(ii)-bus_U(ii)*bus_U(jj)*(G(ii,jj)*sin(bus_e(ii)-bus _e(

7、jj)-B(ii,jj)*cos(bus_e(ii)-bus_e(jj);endendUQ=diag(bus_U(1:(nb-nPV-nSW);error_Q=UQdelta_Q;if max(abs(error_Q)0.00001delta_U=-B2error_Q;bus_U=bus_U+delta_U;zeros(nPV+nSW),1);c=c+1;KP=1;else KQ=0;if KP=0elsebreakend endend %至此得到收敛的节点电压值 % 对计算结果进行数据处理 %将节点结果用原节点编号表示 bus_Ue=zeros(nb,3);bus_Ue(:,1,2,3)=n

8、odenum(:,2) bus_U bus_e/pi*180;for ii=1:nbfor jj=ii+1:nbif bus_Ue(ii,1)bus_Ue(jj,1) t=bus_Ue(ii,:); bus_Ue(ii,:)=bus_Ue(jj,:); bus_Ue(jj,:)=t;endendend%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(

9、r_U);%计算各支路功率 Sij line_S=zeros(nb,nb);line_S0=zeros(nb,nb);for ii=1:nbfor jj=1:nbline_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:nbfor jj=1:nb line_S0(nodenum(ii,2),nodenum(jj,2)=line_S(ii,jj); endendline_P=real(line_S0);

10、line_Q=imag(line_S0); %计算各支路损耗delta_S=zeros(nl,1);for k=1:nla=linenum(k,1);b=linenum(k,2); delta_S(k)=line_S0(a,b)+line_S(b,a);end %计算网络总损耗 S0=sum(delta_S);%将计算结果输入指定文件ieee14_out.txt, wt );fprintf(fid,for k=1:nb 节点号 t 节点电压幅值 t 节点电压角度 n );fprintf(fid,end%dt%ft%fn ,k,bus_Ue(k,1),bus_Ue(k,2);fprintf(fi

11、d,for k=1:nl 支路首端 t 支路末端 t 支路有功 t 支路无功 t 支路损耗 nfid=fopen( C:UserslrDesktopmatlab练习 训练题 大电网潮流计算fprintf(fid, %dtt%dtt%ft%ft%fn ,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);endfprintf(fid, 平衡节点功率 =%fn ,SW_S);fprintf(fid, 网络总损耗 =%fn ,S0);fclose(

12、fid);end%此函数用来对原始输入节点、线路数据进行重新编号% %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 平 衡

13、节点的个数for k=1:nbswitch 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,:);otherwisedisp( 节点数据类型出错！ );endend%生成重新编号后的节点数据矩阵bus=PQ;PV;SW;nodenum=(1:nb) bus(:,1);%第一列为新的节点编号，第二列为对应的旧节点编号bus(:,1)=(1:nb);%至此实现了节点数据的重新编号%对线路数据重新编号% nodenum=(

14、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);endendfunctionY = admittance(bus,line,c )%此函数用来形成节点导纳矩阵% %bus 数组 1. 节点编号 2. 节点电压 3. 节点电压角度 4. 注入有功 5. 注入无功 6. 节点类型 (1PQ 2PV 3 平衡 )%l

15、ine 数组 1. 始端节点编号 2. 末端节点编号 3. 电阻 4 电抗 5 电导 6 电纳 B/2 7. 变比%c 是用来控制形成节点导纳矩阵的方式的， c=1 形成一般的节点导纳矩阵， 可以用来确定 B2， c=2, 形成化简条件 3 的节点导纳矩阵，确定 B1nb,=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 ccase 1for k=1:nlzt(k)=line

16、(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);endcase 2for k=1:nlyt(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)=1K(k)=1;else K(k)=line(k,7);endendcase 3for k=1:nlyt(k)=1/(1i*line(k,4);ym(k)=line(k

17、,5)+1i*line(k,6);I(k)=line(k,1);J(k)=line(k,2);K(k)=line(k,7);endend%针对不同线路设置节点导纳值for k=1:nlif (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);endif (K(k)=1)&(J(k)=0)%对地支路Y(I(k),I(k)=Y(I(k),I(k)+ym(k);endif (K(k)1)&(J(k)=0)%非标准变比在 j 侧的变压器支路 ,zt ym 为折算到 i 侧的值Y(I(k),I(k)=Y(I(k),I(k)