matlab電力系統(tǒng)快速解耦法潮流計(jì)算及短路計(jì)算程序

1、電力系統(tǒng)快速解耦法潮流分析及短路計(jì)算 一程序設(shè)計(jì)的基本思想：（1） 由于電力系統(tǒng)潮流分析中要利用到矩陣運(yùn)算，復(fù)數(shù)運(yùn)算，故采用matlab編程。采用文件輸入，將系統(tǒng)的各個(gè)參數(shù)以文件的形式輸入，便于程序的通用化。（2） 本程序共有兩個(gè)輸入文件，分別為線路參數(shù)的文件，和已知的節(jié)點(diǎn)狀態(tài)文件（pq）（3） 為了使程序不僅僅局限于計(jì)算9節(jié)點(diǎn)網(wǎng)絡(luò)，在形成節(jié)點(diǎn)導(dǎo)納的函數(shù)yn（）中，利用循環(huán)，找出線路首節(jié)點(diǎn)中的最大編號(hào)數(shù)，自動(dòng)確定節(jié)點(diǎn)導(dǎo)納矩陣的維數(shù)。故對(duì)于任意n節(jié)點(diǎn)網(wǎng)絡(luò)，均可以計(jì)算出節(jié)點(diǎn)導(dǎo)納矩陣（4） 在（3）的前提下，為了使程序支持系統(tǒng)增加節(jié)點(diǎn)，增加負(fù)荷等造成的pq參數(shù)改變，或者pq表的加長。對(duì)程序做了如下優(yōu)

2、化。首先，程序執(zhí)行的基礎(chǔ)是pq表中平衡節(jié)點(diǎn)在第一行，接下來是pv節(jié)點(diǎn)，最后是pq節(jié)點(diǎn)，如果系統(tǒng)添加節(jié)點(diǎn)，或者刪除節(jié)點(diǎn)，均在pq表的末端操作，會(huì)造成pq表的順序不是平衡節(jié)點(diǎn)、pv節(jié)點(diǎn)、pq節(jié)點(diǎn)的順序。故引入了seqencing（）函數(shù)，其作用就是不論輸入的pq表是什么順序，在程序讀入后均按平衡-pv-pq的順序排列。其次，順序打亂的pq表必須與支路參數(shù)表對(duì)應(yīng)，故在yn（）函數(shù)中加入了兩段循環(huán)體，使之對(duì)應(yīng)（見相應(yīng)函數(shù)體注釋）（5） 在滿足了上述4個(gè)條件后，程序便可以通用化了。當(dāng)然，由于水平有限，且程序未能由大量數(shù)據(jù)測試，故缺陷在所難免，這里僅是做了通用化的嘗試。在本文最末附加了該程序通用化的實(shí)例。

3、二、潮流計(jì)算框圖輸入系統(tǒng)的支路參數(shù)和pq參數(shù)計(jì)算節(jié)點(diǎn)導(dǎo)納矩陣設(shè)置初始電壓值u0設(shè)置收斂條件及最大迭代次數(shù)kmax判斷是否收斂求解修正方程得到u和對(duì)u0進(jìn)行修正，u0=u0+ukkmax ?ny輸出結(jié)果結(jié)束輸出潮流不收斂信息yn開始三定義相應(yīng)的函數(shù)1.形成節(jié)點(diǎn)導(dǎo)納矩陣的函數(shù)yn（）function y=yn(x,y) %定義一名為yn的函數(shù)，其功能是自動(dòng)識(shí)別輸入表中節(jié)點(diǎn)的個(gè)數(shù)，形成相應(yīng)的節(jié)點(diǎn)導(dǎo)納矩陣fid,message=fopen(x,r) ; %從x文件中讀入支路參數(shù)if fid=-1; %判斷文件是否正確打開 error(message);end;headpoint,headnumber,

4、 endpoint,endnumber,r,x,b,k=textread(x,%s %d %s %d %f %f %f %f); %將讀入的參數(shù)處理為以列為向量的數(shù)組fclose(fid);%關(guān)閉文件l=length(headnumber); %確定輸入表的行數(shù) pointstyle,pointnumber,ps,qs,uk,ok=seqencing(y); %調(diào)用seqencing函數(shù)，引入y文件中的pq參數(shù)a=pointnumber;for i=1:l; %通過以下兩循環(huán)體，實(shí)現(xiàn)pq參數(shù)與支路參數(shù)的編號(hào)對(duì)應(yīng) for j=1:l; if headnumber(i)=j; headnumber

5、(i)=a(j); break; end; end;end;for i=1:l; for j=1:l; if endnumber(i)=j; endnumber(i)=a(j); break; end; end;end;y=zeros(l,l); %根據(jù)txt文件中數(shù)據(jù)表的長度建立空的節(jié)點(diǎn)導(dǎo)納矩陣for i=1:l m=headnumber(i);n=endnumber(i); if k(i)=0; %判斷是否何種元件，為輸電線元件 if n=0; y(m,m)=y(m,m)+1j*b(i)+1/(r(i)+1j*x(i); y(n,n)=y(n,n)+1j*b(i)+1/(r(i)+1j*x

6、(i); y(m,n)=y(m,n)-1/(r(i)+1j*x(i); y(n,m)=y(n,m)-1/(r(i)+1j*x(i); else y(m,m)=y(m,m)+r(i)+1j*x(i); end; else %為變壓器元件 if n=0; y(m,m)=y(m,m)+1/(r(i)+1j*x(i); y(m,n)=y(m,n)-1/(k(i)*(r(i)+1j*x(i); y(n,n)=y(n,n)+1/(k(i)*k(i)*(r(i)+1j*x(i); y(n,m)=y(n,m)-1/(k(i)*(r(i)+1j*x(i); else y(m,m)=y(m,m)+r(i)+1j*

7、x(i); end; end; end;maxm=headnumber(1);%通過下面兩個(gè)循環(huán)體，確定輸入表中節(jié)點(diǎn)編號(hào)的最大值，及為節(jié)點(diǎn)導(dǎo)納矩陣的維數(shù)for i=1:l; if maxm=headnumber(i); maxm=headnumber(i); end;end;maxn=endnumber(1);for i=1:l; if maxnpointstyle(j+1); t=pointstyle(j+1); pointstyle(j+1)=pointstyle(j); pointstyle(j)=t; t=pointnumber(j+1); pointnumber(j+1)=point

8、number(j); pointnumber(j)=t; t=ps(j+1); ps(j+1)=ps(j); ps(j)=t; t=qs(j+1); qs(j+1)=qs(j); qs(j)=t; t=uk(j+1); uk(j+1)=uk(j); uk(j)=t; end; end;end; 3、形成解耦算法b矩陣的函數(shù) formb1（）function b1=formb1(x,y)%定義名為b1的函數(shù)形成解耦算法中的b矩陣，得到的b矩陣用b1表示fid,message=fopen(x,r) ; %從x文件中讀入支路參數(shù)if fid=-1; %判斷文件是否正確打開 error(message

9、);end;headpoint,headnumber,endpoint,endnumber,r,x,b,k=textread(x,%s %d %s %d %f %f %f %f); %將讀入的參數(shù)處理為以列為向量的數(shù)組fclose(fid);l=length(headnumber); pointstyle,pointnumber,ps,qs,uk,ok=seqencing(y);%調(diào)用seqencing函數(shù)，引入y文件中的pq參數(shù)a=pointnumber;%通過以下兩循環(huán)體，實(shí)現(xiàn)pq參數(shù)與支路參數(shù)的編號(hào)對(duì)應(yīng)for i=1:l; for j=1:l; if headnumber(i)=j; h

10、eadnumber(i)=a(j); break; end; end;end;for i=1:l; for j=1:l; if endnumber(i)=j; endnumber(i)=a(j); break; end; end;end;b1=zeros(l,l);for i=1:l %以行為單位，通過循環(huán)，用支路參數(shù)對(duì)b1進(jìn)行修改，形成b矩陣 m=headnumber(i);n=endnumber(i); b1(m,m)=b1(m,m)-1/x(i); b1(n,n)=b1(n,n)-1/x(i); b1(m,n)=b1(m,n)+1/x(i); b1(n,m)=b1(n,m)+1/x(i)

11、;endmaxm=headnumber(1);for i=1:l; if maxm=headnumber(i); maxm=headnumber(i); end;end;maxn=endnumber(1);for i=1:l; if maxn=endnumber(i); maxn=endnumber(i); end;end;b1=b1(2:max(maxm,maxn),2:max(maxm,maxn); %形成b矩陣 4、形成解耦算法b矩陣的函數(shù) formb11（）function b11=formb11(x,y)%定義名為b11的函數(shù)形成解耦算法中b矩陣，用b11表示從x文件中讀入支路參數(shù)確

12、定y，從y文件中讀入pq參數(shù)確定b11的維數(shù)，即除去平衡節(jié)點(diǎn)和pv節(jié)點(diǎn)，此處要求pq參數(shù)錄入時(shí)，將平衡節(jié)點(diǎn)和pq節(jié)點(diǎn)放在前排，這一要求在yn函數(shù)中通過seqencing函數(shù)已經(jīng)滿足y=yn(x,y);b=imag(y);pointstyle,pointnumber,ps,qs,uk,ok=seqencing(y);i=1;j=1;while pointstyle(i)=1; i=i+1; j=j+1;end;b11=b(j:end,j:end); %形成b矩陣5、計(jì)算正常情況下系統(tǒng)節(jié)點(diǎn)電壓的函數(shù) powerflow（）function u0,o0=powerflow(x,y)%定義名為powe

13、rflow的函數(shù)，利用快速解耦算法來計(jì)算正常情況下系統(tǒng)內(nèi)各個(gè)節(jié)點(diǎn)的電壓和相角pointstyle,pointnumber,ps,qs,uk,ok=seqencing(y);% 調(diào)用seqencing函數(shù)對(duì)pq參數(shù)表進(jìn)行排序 y=yn(x,y); %形成節(jié)點(diǎn)導(dǎo)納矩陣，yn為n維 b1=formb1(x,y); %形成解耦算法中的b矩陣，b1為n-1維 b11=formb11(x,y); %形成解耦算法中的b矩陣，b為m維 g=real(y); %取y的實(shí)部 b=imag(y); %取y的虛部 u0=uk; o0=ok; l=length(pointnumber); p=zeros(l,1); q

14、=zeros(l,1); dp=zeros(l,1); dq=zeros(l,1); number=1; i=1;k=1; while pointstyle(i)=1; %通過k值確定系統(tǒng)中pq節(jié)點(diǎn)的個(gè)數(shù) i=i+1; k=k+1;end; while number100 %定義迭代次數(shù)上限為100次for i=2:l; sum1=0; for j=1:l; sum1=sum1+u0(j)*(g(i,j)*cos(o0(i)-o0(j)+b(i,j)*sin(o0(i)-o0(j); %潮流方程，n-1維 end; dp(i)=ps(i)-u0(i)*sum1;endfor i=k:l; su

15、m2=0; for j=1:l; sum2=sum2+u0(j)*(g(i,j)*sin(o0(i)-o0(j)-b(i,j)*cos(o0(i)-o0(j); %潮流方程，m維 end; dq(i)=qs(i)-u0(i)*sum2;end dp1=dp(2:l)./u0(2:l); dq1=dq(k:l)./u0(k:l); a=max(norm(dp1,inf); b=max(norm(dq1,inf); if max(a,b)0.00001 %判斷是否收斂 break;disp(迭代)disp(k);disp(次后收斂); else %如不收斂， do=-inv(b1)*dp1; %d

16、o為n-1維 du=-inv(b11)*dq1; %du為m維 zero1=zeros(k-1,1); zero2=0; du=zero1;du; do=zero2;do; u0=u0+du; o0=o0+do; number=number+1; end; if number=100; disp(迭代100次后不收斂，迭代結(jié)束); end;end;四對(duì)相應(yīng)的系統(tǒng)進(jìn)行潮流分析和短路計(jì)算定義完上述函數(shù)之后，可直接調(diào)用函數(shù)形成導(dǎo)納矩陣，計(jì)算正常情況下的節(jié)點(diǎn)電壓，進(jìn)行短路計(jì)算計(jì)算短路電流，短路后各個(gè)節(jié)點(diǎn)電壓以及支路潮流分布。程序的輸入表共有三個(gè)，為network.txt ，network2.txt ，

17、pq參數(shù).txt。其中network.txt存放不含發(fā)電機(jī)節(jié)點(diǎn)導(dǎo)納和負(fù)荷導(dǎo)納的系統(tǒng)支路參數(shù)，network2.txt存放含發(fā)電機(jī)節(jié)點(diǎn)導(dǎo)納和負(fù)荷導(dǎo)納的系統(tǒng)支路參數(shù)，pq參數(shù).txt存放系統(tǒng)的pq參數(shù)，現(xiàn)列表如下第一列：首節(jié)點(diǎn)名稱 第二列：首節(jié)點(diǎn)編號(hào) 第三列：尾節(jié)點(diǎn)名稱 第四列：尾節(jié)點(diǎn)編號(hào)第五列：支路參數(shù)r 第六列：支路參數(shù)x 第七列：支路參數(shù)b/2 第八列：變比k第一列：節(jié)點(diǎn)類型，0-平衡節(jié)點(diǎn)，1-pv節(jié)點(diǎn)，2-pq節(jié)點(diǎn) 第二列：節(jié)點(diǎn)編號(hào)第三列：節(jié)點(diǎn)的有功p 第四列：節(jié)點(diǎn)的無功q 第五列：節(jié)點(diǎn)電壓，pq節(jié)點(diǎn)置為0第六列：節(jié)點(diǎn)電壓的相角，初始值為0程序代碼如下：%正常情況下節(jié)點(diǎn)導(dǎo)納矩陣y=yn(n

18、etwork.txt,pq參數(shù).txt);disp(節(jié)點(diǎn)導(dǎo)納矩陣)disp(sparse(y);%正常情況下潮流計(jì)算結(jié)果兵繪制潮流計(jì)算框圖 u0,o0=powerflow(network.txt,pq參數(shù).txt) %各個(gè)節(jié)點(diǎn)電壓為u0，相角為o0disp(正常情況下系統(tǒng)各個(gè)節(jié)點(diǎn)電壓);disp(u0);disp(正常情況下系統(tǒng)各個(gè)節(jié)點(diǎn)相角);disp(o0);o0=o0*180/pidisp(正常情況下系統(tǒng)各個(gè)節(jié)點(diǎn)相角角度值)%計(jì)算平衡節(jié)點(diǎn)功率和pv節(jié)點(diǎn)無功功率：for m=1:9 p1=p1+u0(1)*u0(m)*(g(1,m)*cos(o0(1)-o0(m)+b(1,m)*sin(o0

19、(1)-o0(m); q1=q1+u0(1)*u0(m)*(g(1,m)*sin(o0(1)-o0(m)-b(1,m)*cos(o0(1)-o0(m);q2=q2+u0(2)*u0(m)*(g(2,m)*sin(o0(2)-o0(m)-b(2,m)*cos(o0(2)-o0(m);q3=q3+u0(3)*u0(m)*(g(3,m)*sin(o0(3)-o0(m)-b(3,m)*cos(o0(3)-o0(m);end;disp(平衡節(jié)點(diǎn)有功、無功，兩個(gè)電壓節(jié)點(diǎn)有功分別為)disp(p1); %有功功率和無功功率disp(q1); %平衡節(jié)點(diǎn)無功功率disp(q2); %pv節(jié)點(diǎn)2無功功率disp

20、(q3); %pv節(jié)點(diǎn)3無功功率%增加發(fā)電機(jī)導(dǎo)納yi和負(fù)荷導(dǎo)納yldi后的發(fā)電機(jī)節(jié)點(diǎn)和負(fù)荷節(jié)點(diǎn)的自導(dǎo)納 y=yn(network2.txt,pq參數(shù).txt) %network2文件中存放的是經(jīng)修改的支路參數(shù)disp(加發(fā)電機(jī)導(dǎo)納yi和負(fù)荷導(dǎo)納yldi后的發(fā)電機(jī)節(jié)點(diǎn)和負(fù)荷節(jié)點(diǎn)的自導(dǎo)納);diag(y)%精確法計(jì)算短路電流，短路后各節(jié)點(diǎn)電壓以及網(wǎng)絡(luò)中各支路的電流分布y=yn(network2.txt,pq參數(shù).txt); %形成包括發(fā)電機(jī)內(nèi)阻抗和符合阻抗節(jié)點(diǎn)導(dǎo)納矩陣 z=inv(y); %形成節(jié)點(diǎn)阻抗矩陣 %題中給出的是4節(jié)點(diǎn)短路zf=z(:,4);u0,o0=powerflow(network

21、.txt,pq參數(shù).txt);u01=u0.*cos(o0)+1j*u0.*sin(o0);if=u01(4)/z(4,4); %計(jì)算短路電流ifu=u01-zf*if; %計(jì)算短路后的各個(gè)節(jié)點(diǎn)電壓fid,message=fopen(network2.txt,r) ; %從x文件中讀入支路參數(shù)if fid=-1; %判斷文件是否正確打開 error(message);end;headpoint,headnumber,endpoint,endnumber,r,x,b,k=textread(network2.txt,%s %d %s %d %f %f %f %f); row=headnumber,

22、endnumber;for i=1:9 %計(jì)算系統(tǒng)中各個(gè)支路的潮流分布 i(row(i,1),row(i,2)=(u(row(i,1)-u(row(i,2)/z(row(i,1),row(i,2);end;disp(阻抗矩陣中的第f列);disp(精確算法短路電流 if);disp(if模值);disp(abs(if);disp(相角為);disp(angle(if);disp(短路后各個(gè)節(jié)點(diǎn)的電壓);disp(u);disp(模值為);disp(短路后各個(gè)支路的電流值 i) ;disp(模值為 i);disp(相角為（度）); %近似計(jì)算短路電流，短路后各節(jié)點(diǎn)電壓以及網(wǎng)絡(luò)中各支路的電流分布i

23、f2=1/z(4,4); %計(jì)算短路電流ifu2=1-zf*if; %計(jì)算短路后的各個(gè)節(jié)點(diǎn)電壓disp(近似算法短路電流 if2);disp(if2模值);disp(abs(if);disp(相角為（度）);disp(angle(if2);disp(短路后各個(gè)節(jié)點(diǎn)的電壓);disp(u2);disp(模值為);程序運(yùn)行結(jié)果如下：節(jié)點(diǎn)導(dǎo)納矩陣y = (1,1) 0 -17.3611i (4,1) 0 +17.3611i (2,2) 0 -16.0000i (7,2) 0 +16.0000i (3,3) 0 -17.0648i (9,3) 0 +17.0648i (1,4) 0 +17.3611i

24、 (4,4) 3.3074 -39.3089i (5,4) -1.3652 +11.6041i (6,4) -1.9422 +10.5107i (4,5) -1.3652 +11.6041i (5,5) 2.5528 -17.3382i (7,5) -1.1876 + 5.9751i (4,6) -1.9422 +10.5107i (6,6) 3.2242 -15.8409i (9,6) -1.2820 + 5.5882i (2,7) 0 +16.0000i (5,7) -1.1876 + 5.9751i (7,7) 2.8047 -35.4456i (8,7) -1.6171 +13.698

25、0i (7,8) -1.6171 +13.6980i (8,8) 2.7722 -23.3032i (9,8) -1.1551 + 9.7843i (3,9) 0 +17.0648i (6,9) -1.2820 + 5.5882i (8,9) -1.1551 + 9.7843i (9,9) 2.4371 -32.1539i迭代 k=9次后收斂正常情況下系統(tǒng)各個(gè)節(jié)點(diǎn)電壓u0 = 1.0400 1.0250 1.0250 1.0258 0.9956 1.0127 1.0258 1.0159 1.0324正常情況下系統(tǒng)各個(gè)節(jié)點(diǎn)相角o0 = 0 0.1620 0.0814 -0.0387 -0.069

26、6 -0.0644 0.0649 0.0127 0.0343正常情況下系統(tǒng)各個(gè)節(jié)點(diǎn)相角角度值o0 = 0 9.2800 4.6647 -2.2168 -3.9888 -3.6875 3.7197 0.7275 1.9667平衡節(jié)點(diǎn)有功、無功，兩個(gè)電壓節(jié)點(diǎn)有功分別為 p1=0.7164 q1= 0.2705q2= 0.0665 q3= -0.1086加發(fā)電機(jī)導(dǎo)納yi和負(fù)荷導(dǎo)納yldi后的發(fā)電機(jī)節(jié)點(diǎn)和負(fù)荷節(jié)點(diǎn)的自導(dǎo)納ans = 0 -20.6944i 0 -19.3333i 0 -20.3982i 3.3074 -39.3089i 1.2917 -16.8338i 2.3466 -15.5484i

27、 2.8047 -35.4456i 1.8033 -22.9641i2.4371 -32.1539i 阻抗矩陣中的第f列zf = -0.0630 + 0.1475i -0.0645 + 0.0860i -0.0634 + 0.0864i -0.0751 + 0.1758i -0.0892 + 0.1534i -0.0854 + 0.1532i -0.0779 + 0.1040i -0.0832 + 0.1032i -0.0758 + 0.1032i精確算法短路電流 if = -2.2963 - 4.8486iif模值為 5.3649 相角為 -115.3422短路后各個(gè)節(jié)點(diǎn)的電壓 u = 0.

28、1801 + 0.0333i 0.4463 + 0.0502i 0.4573 - 0.0257i -0.0000 - 0.0000i 0.0446 - 0.1498i 0.0716 - 0.1271i 0.3406 - 0.0725i 0.3243 - 0.1533i 0.3572 - 0.0950i模值為0.1831 0.4491 0.4580 0.0000 0.1563 0.1459 0.3482 0.3587 0.3696短路后各個(gè)支路的電流值 i(1,4) -0.2501 - 1.1140i (4,5) 0.8560 - 0.2072i (4,6) 0.8320 + 0.0041i (

29、2,7) 0.5086 - 0.9694i (5,7) 0.6818 + 1.8688i (7,8) 0.3965 - 0.3175i (3,9) 0.2141 - 0.7887i (6,9) 0.8241 + 1.7182i (8,9) -0.2355 + 0.3676i模值為 i(1,4) 1.1417 (4,5) 0.8807 (4,6) 0.8320 (2,7) 1.0947 (5,7) 1.9893 (7,8) 0.5079 (3,9) 0.8172 (6,9) 1.9056 (8,9) 0.4365相角為（度）(1,4) -102.6530 (4,5) -13.6077 (4,6)

30、 0.2855 (2,7) -62.3148 (5,7) 69.9555 (7,8) -38.6844 (3,9) -74.8141 (6,9) 64.3752 (8,9) 122.6426近似算法短路電流 if2 = 0.2105 - 5.2357iif2模值 5.2399 相角為（度） -87.6981短路后各個(gè)節(jié)點(diǎn)電壓值u2 = 0.1611 + 0.0000i 0.4785 - 0.0351i 0.4773 - 0.0384i 0.0000 0.1155 - 0.0063i 0.1178 - 0.0098i 0.3698 - 0.0424i 0.3673 - 0.0436i 0.3752 - 0.0459i模值為 0.1611 0.4797 0.4788 0.0000 0.1157 0.1182 0.3722 0.3698 0.3780精確算法與近似算法的比較分析（1） 短路電流誤差（if精-if近）/