PQ分解法計算大電網(wǎng)潮流程序

1、精選優(yōu)質(zhì)文檔-傾情為你奉上function PQ%用PQ分解法計算大電網(wǎng)潮流% %bus數(shù)組 1.節(jié)點編號 2.節(jié)點電壓 3.節(jié)點電壓角度 4.注入有功 5.注入無功 6.節(jié)點類型(1PQ 2PV 3平衡) %line數(shù)組 1.始端節(jié)點編號 2.末端節(jié)點編號 3.電阻 4電抗 5電導(dǎo)G 6電納B 7.變比 %打開數(shù)據(jù)文件clearclcbus=load('ieee14bus.txt');line=load('ieee14line.txt');linenum(:,1,2)=line(:,1,2);nb,=size(bus);nl,=size(line);% nod

2、enum=(1:nb)' bus(:,1);%帶入子函數(shù)數(shù)據(jù)處理bus,line,nPQ,nPV,nSW,nodenum =change1_busline( bus,line );%對節(jié)點重新編號Y = admittance(bus,line,1 );%生成節(jié)點導(dǎo)納矩陣Y1= admittance(bus,line,2 );%生成化簡條件3的矩陣B1Y2=admittance(bus,line,3 );%生成化簡條件3的矩陣B2 %-% %臨時添加的測試數(shù)據(jù)% nPQ=4; nPV=0;nSW=1;nb=5;% Y=10.834-32.5i -1.667+5i -1.667+5i -2

3、.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

4、-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

5、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=

6、real(Y);B=imag(Y);B10=imag(Y1);B20=imag(Y2);%矩陣B0是進(jìn)行化簡三后的節(jié)點導(dǎo)納矩陣虛部 %形成解耦潮流的系數(shù)矩陣B1和B2B1=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(

7、ii)-bus_e(jj); endendUP=diag(bus_U(1:(nb-1);%U矩陣?yán)酶鞴?jié)點電壓形成對角陣，來計算修正方程，對角線上的元素與bus_U列元素一一對應(yīng)error_P=UPdelta_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 endend for ii=1:nPQ delta_Q(ii)=bus(ii,5); for jj=1:nb delta_Q(ii)=delta

8、_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); endend UQ=diag(bus_U(1:(nb-nPV-nSW);error_Q=UQdelta_Q; if max(abs(error_Q)>0.00001 delta_U=-B2error_Q; bus_U=bus_U+delta_U;zeros(nPV+nSW),1); c=c+1; KP=1;else KQ=0; if KP=0 else break endend end%至此得到收斂的

9、節(jié)點電壓值%- %-% 對計算結(jié)果進(jìn)行數(shù)據(jù)處理%將節(jié)點結(jié)果用原節(jié)點編號表示bus_Ue=zeros(nb,3);bus_Ue(:,1,2,3)=nodenum(:,2) bus_U bus_e/pi*180;for ii=1:nb for jj=ii+1:nb if bus_Ue(ii,1)>bus_Ue(jj,1) t=bus_Ue(ii,:); bus_Ue(ii,:)=bus_Ue(jj,:); bus_Ue(jj,:)=t; end endend %r_U是收斂的電壓表達(dá)成復(fù)數(shù)的形式r_U=zeros(nb,1); for k=1:nbr_U(k)=bus_U(k)*(cos(b

10、us_e(k)+1i*sin(bus_e(k);end %計算平衡節(jié)點功率SW_S=0;SW_S=SW_S+r_U(nb)*conj(Y(nb,:)*conj(r_U); %計算各支路功率Sijline_S=zeros(nb,nb);line_S0=zeros(nb,nb);for ii=1:nb for jj=1:nb line_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); endend %-%把線路結(jié)果還原成原節(jié)點編號對應(yīng)的結(jié)果for ii=1:nb for jj

11、=1:nb line_S0(nodenum(ii,2),nodenum(jj,2)=line_S(ii,jj); endend line_P=real(line_S0);line_Q=imag(line_S0);%計算各支路損耗delta_S=zeros(nl,1);for k=1:nl a=linenum(k,1);b=linenum(k,2); delta_S(k)=line_S0(a,b)+line_S(b,a);end%計算網(wǎng)絡(luò)總損耗S0=sum(delta_S);%- %將計算結(jié)果輸入指定文件fid=fopen('C:UserslrDesktopmatlab練習(xí)訓(xùn)練題大電網(wǎng)潮

12、流計算ieee14_out.txt','wt');fprintf(fid,'節(jié)點號t節(jié)點電壓幅值t節(jié)點電壓角度n');for k=1:nb fprintf(fid,'%dt%ft%fn',k,bus_Ue(k,1),bus_Ue(k,2);endfprintf(fid,'支路首端t支路末端t支路有功t支路無功t支路損耗n');for k=1:nl fprintf(fid,'%dtt%dtt%ft%ft%fn',linenum(k,1),linenum(k,2),line_P(linenum(k,1),lin

13、enum(k,2),line_Q(linenum(k,1),linenum(k,2),delta_S(k);endfprintf(fid,'平衡節(jié)點功率=%fn',SW_S);fprintf(fid,'網(wǎng)絡(luò)總損耗=%fn',S0);fclose(fid);end function bus,line,nPQ,nPV,nSW,nodenum =change1_busline( bus,line )%此函數(shù)用來對原始輸入節(jié)點、線路數(shù)據(jù)進(jìn)行重新編號 % %bus數(shù)組 1.節(jié)點編號 2.節(jié)點電壓 3.節(jié)點電壓角度 4.注入有功 5.注入無功 6.節(jié)點類型(1PQ 2PV

14、3平衡) %line數(shù)組 1.始端節(jié)點編號 2.末端節(jié)點編號 3.電阻 4電抗 5電導(dǎo)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 平衡節(jié)點的個數(shù)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,:);

15、case 3 nSW=nSW+1; SW(nSW,:)=bus(k,:); otherwise disp('節(jié)點數(shù)據(jù)類型出錯！');endend %生成重新編號后的節(jié)點數(shù)據(jù)矩陣 bus=PQ;PV;SW; nodenum=(1:nb)' bus(:,1);%第一列為新的節(jié)點編號，第二列為對應(yīng)的舊節(jié)點編號 bus(:,1)=(1:nb)' %至此實現(xiàn)了節(jié)點數(shù)據(jù)的重新編號%- %對線路數(shù)據(jù)重新編號 % nodenum=(1:nb)' bus(:,1);%第一列為新的節(jié)點編號，第二列為對應(yīng)的舊節(jié)點編號 for ii=1:nl r1,=find(nodenum(

16、:,2)=line(ii,1); line(ii,1)=nodenum(r1,1); r2,=find(nodenum(:,2)=line(ii,2); line(ii,2)=nodenum(r2,1); end end function Y = admittance(bus,line,c )%此函數(shù)用來形成節(jié)點導(dǎo)納矩陣 % %bus數(shù)組 1.節(jié)點編號 2.節(jié)點電壓 3.節(jié)點電壓角度 4.注入有功 5.注入無功 6.節(jié)點類型(1PQ 2PV 3平衡) %line數(shù)組 1.始端節(jié)點編號 2.末端節(jié)點編號 3.電阻 4電抗 5電導(dǎo) 6電納B/2 7.變比 %c是用來控制形成節(jié)點導(dǎo)納矩陣的方式的，c

17、=1 形成一般的節(jié)點導(dǎo)納矩陣，可以用來確定B2，c=2,形成化簡條件3的節(jié)點導(dǎo)納矩陣，確定B1 nb,=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 c case 1 for k=1:nl zt(k)=line(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)

18、=line(k,2);K(k)=line(k,7); end case 2 for k=1:nl yt(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)=1 K(k)=1; else K(k)=line(k,7); end end case 3 for k=1:nl yt(k)=1/(1i*line(k,4); ym(k)=line(k,5)+1i*line(k,6); I(k)=line(k,1);J(k)=line(k,2);K(k)=line(k,7); endend%針對不同線路設(shè)置節(jié)點導(dǎo)納值for k=1:nlif (K(k)=1)&&am