中國(guó)石油大學(xué)電力系統(tǒng)分析大作業(yè)

1、 電力系統(tǒng)分析大作業(yè) 編程實(shí)現(xiàn)N-R迭代和P-Q迭代潮流計(jì)算 院 系： 信息與控制工程學(xué)院 班 級(jí)： 電氣12-04班 學(xué) 號(hào)： 120534 姓 名： 日期：2015年 06 月 17 日要求：利用NR迭代和P-Q迭代計(jì)算所給系統(tǒng)的潮流，設(shè)發(fā)電機(jī)G1的端電壓為1 p.u.，發(fā)出的有功、無(wú)功可調(diào)；發(fā)電機(jī)G2的端電壓為1 p.u.，按指定的有功P=0.5 p.u.發(fā)電，取=10xe-4。1:牛頓-拉夫遜法潮流算法1.1牛頓-拉夫遜法潮流算法的基本原理 電力系統(tǒng)潮流計(jì)算是電力系統(tǒng)運(yùn)行和規(guī)劃中最基本和最經(jīng)常的計(jì)算，其任務(wù)是在已知某些運(yùn)行參數(shù)的情況下，計(jì)算出系統(tǒng)中全部的運(yùn)行參數(shù)，一般來(lái)說(shuō)，各個(gè)母線所供

2、負(fù)荷的功率是已知的，各個(gè)節(jié)點(diǎn)電壓是未知的（平衡節(jié)點(diǎn)除外），可以根據(jù)網(wǎng)絡(luò)結(jié)構(gòu)形成節(jié)點(diǎn)導(dǎo)納矩陣，然后由節(jié)點(diǎn)導(dǎo)納矩陣和網(wǎng)絡(luò)拓?fù)浣Y(jié)構(gòu)列寫功率方程，由于功率方程里功率是已知的，電壓的幅值和相角是未知的，這樣潮流計(jì)算的問(wèn)題就轉(zhuǎn)化為求解非線性方程組的問(wèn)題了。為了便于用迭代法解方程組，需要將上述功率方程改寫成功率平衡方程，并對(duì)功率平衡方程求偏導(dǎo)，得出對(duì)應(yīng)的雅可比矩陣，給未知節(jié)點(diǎn)賦電壓初值，一般為額定電壓，將初值帶入功率平衡方程，得到功率不平衡量，這樣由功率不平衡量、雅可比矩陣、節(jié)點(diǎn)電壓不平衡量（未知的）構(gòu)成了誤差方程，解誤差方程，得到節(jié)點(diǎn)電壓不平衡量，節(jié)點(diǎn)電壓加上節(jié)點(diǎn)電壓不平衡量構(gòu)成新的節(jié)點(diǎn)電壓初值，將新的初

3、值帶入原來(lái)的功率平衡方程，并重新形成雅可比矩陣，然后計(jì)算新的電壓不平衡量，這樣不斷迭代，不斷修正，一般迭代三到五次就能收斂。1.2節(jié)點(diǎn)導(dǎo)納矩陣的算法 根據(jù)定義直接求取節(jié)點(diǎn)導(dǎo)納矩陣時(shí)，注意以下幾點(diǎn)： 1) 節(jié)點(diǎn)導(dǎo)納矩陣是方陣，其階數(shù)就等于網(wǎng)絡(luò)中除去參考節(jié)點(diǎn)外的節(jié)點(diǎn)數(shù)。參考節(jié)點(diǎn)一般取大地，編號(hào)為零。 2) 節(jié)點(diǎn)導(dǎo)納矩陣是稀疏矩陣，其各行非零非對(duì)角元素就等于與該行相對(duì)應(yīng)節(jié)點(diǎn)所連接的不接地支路數(shù)。 3) 節(jié)點(diǎn)導(dǎo)納矩陣的對(duì)角元素就等于該節(jié)點(diǎn)所連接導(dǎo)納的總和。因此，與沒(méi)有接地支路的節(jié)點(diǎn)對(duì)應(yīng)的行或列中，對(duì)角元素為非對(duì)角元素之和的負(fù)值。 4) 節(jié)點(diǎn)導(dǎo)納矩陣的非對(duì)角元素等于連接節(jié)點(diǎn)i，j支路導(dǎo)納的負(fù)值。因此，在

4、一般情況下，節(jié)點(diǎn)導(dǎo)納矩陣的對(duì)角元素往往大于非對(duì)角元素的負(fù)值。 5) 節(jié)點(diǎn)導(dǎo)納矩陣一般是對(duì)稱矩陣，這是網(wǎng)絡(luò)的互易特性所決定的。從而，一般只要求取這個(gè)矩陣的上三角或下三角部分。 6) 網(wǎng)絡(luò)中的變壓器。1.3 牛頓拉夫遜迭代法的步驟及算法流程圖：1）形成各節(jié)點(diǎn)導(dǎo)納矩陣Y。2）設(shè)個(gè)節(jié)點(diǎn)電壓的初始值U和相角初始值e 還有迭代次數(shù)初值為0。3）計(jì)算各個(gè)節(jié)點(diǎn)的功率不平衡量。4）根據(jù)收斂條件判斷是否滿足，若不滿足則向下進(jìn)行。5）計(jì)算雅可比矩陣中的各元素。6）修正方程式個(gè)節(jié)點(diǎn)電壓7）利用新值自第3）步開始進(jìn)入下一次迭代，直至達(dá)到精度退出循環(huán)。8）計(jì)算平衡節(jié)點(diǎn)輸出功率和各線路功率流程圖如下：形成節(jié)點(diǎn)導(dǎo)納矩陣設(shè)非平

5、衡節(jié)點(diǎn)電壓初值、e輸入原始數(shù)據(jù)啟動(dòng)對(duì)PQ節(jié)點(diǎn)計(jì)算(對(duì)PV節(jié)點(diǎn)計(jì)算) 計(jì)算雅可比矩陣各元素Hij(k)、Nij(k)、Jij(k)、Lij(k)解修正方程，由及雅可比矩陣用牛頓-拉夫遜法求各節(jié)點(diǎn)的ei（k）、Ui(k)/Ui計(jì)算節(jié)點(diǎn)的新電壓增加迭代次數(shù)k=k+1令迭代次數(shù)k=0計(jì)算平衡節(jié)點(diǎn)的功率及線路功率輸出是否1.4程序代碼及詳細(xì)注釋：clear;disp('【輸出結(jié)果如下】：')m=5;m1=5;isb=5;pr=0.0001;A1=1 2 0.025+0.08i 0.07i 1 0; 1 3 0.03+0.1i 0.09i 1 0; 2 3 0.02+0.06i 0.05i

6、 1 0; 4 2 0.1905i 0 1.0522 1; 5 3 0.1905i 0 1.0522 1;%A1是支路參數(shù)矩陣，第一列和第二列是節(jié)點(diǎn)編號(hào).%第三列為之路串聯(lián)阻抗參數(shù).%第四列為含支路的對(duì)地導(dǎo)納矩陣.%第五列為含變壓器之路的變壓器變比.%第六列表式是否含有變壓器的參數(shù)，1表示含有變壓器，0表示不含有.%變壓器側(cè)為1，否則為0.A2=0 -0.8055-0.5320i 1 0 0 2; 0 -0.18-0.12i 1 0 0 2; 0 0 1 0 0 2; 0 0.5 1 1.0522 0 3; 0 0 1 1.0522 0 1;% A2為節(jié)點(diǎn)參數(shù)矩陣.%第二列為節(jié)點(diǎn)負(fù)荷功率參數(shù).

7、%第三列為節(jié)點(diǎn)電壓參數(shù).%第六列為節(jié)點(diǎn)類型參數(shù)其中1為平衡節(jié)點(diǎn)，2為PQ節(jié)點(diǎn)，3為PV節(jié)點(diǎn).disp('-【 已知量 ：】-')disp('【支路參數(shù)矩陣】')disp(A1)disp('【節(jié)點(diǎn)參數(shù)矩陣】')disp(A2)%=計(jì)算導(dǎo)納矩陣=X=1 0;2 0;3 0;4 0;5 0;Y=zeros(m);U=zeros(1,m);cta=zeros(1,m);V=zeros(1,m);O=zeros(1,m);S1=zeros(m1);for i=1:m if X(i,2)=0; p=X(i,1); Y(p,p)=X(i,2); endendf

8、or i=1:m1 if A1(i,6)=0 p=A1(i,1);q=A1(i,2); else p=A1(i,2);q=A1(i,1); end Y(p,q)=Y(p,q)-1./(A1(i,3)*A1(i,5); Y(q,p)=Y(p,q); Y(q,q)=Y(q,q)+1./A1(i,3)+A1(i,4)./2; Y(p,p)=Y(p,p)+1./(A1(i,3)* A1(i,5)2) +A1(i,4)./2;enddisp('【導(dǎo)納矩陣】 Y=');disp(Y)%輸出導(dǎo)納矩陣%=給定各節(jié)點(diǎn)初始電壓及給定各節(jié)點(diǎn)注入功率=G=real(Y);B=imag(Y);%分解出導(dǎo)

9、納陣的實(shí)部和虛部for i=1:m cta(i)=angle(A2(i,3); U(i)=abs(A2(i,3); V(i)=A2(i,4);%PV、平衡節(jié)點(diǎn)及PQ節(jié)點(diǎn)電壓模值endfor i=1:m %給定各節(jié)點(diǎn)注入功率 S(i)=A2(i,1)-A2(i,2);%i節(jié)點(diǎn)注入功率SG-SL B(i,i)=B(i,i)+A2(i,5);%i節(jié)點(diǎn)無(wú)功補(bǔ)償量(電納值)end%=用牛頓-拉夫遜法迭代求解非線性代數(shù)方程（功率方程）= P=real(S);Q=imag(S);%分解出各節(jié)點(diǎn)注入的有功和無(wú)功功率ICT1=0;IT2=1;%迭代次數(shù)ICT1；不滿足收斂要求的節(jié)點(diǎn)數(shù)IT2while IT2=0

10、 IT2=0;t1=1;t2=1; for i=1:m if i=isb C(i)=0; D(i)=0; for j1=1:m %第i行共m列(m個(gè)節(jié)點(diǎn)間互導(dǎo)納及節(jié)點(diǎn)電壓相乘即電流) C(i)=C(i)+U(i)*U(j1)*(G(i,j1)*cos(cta(i)-cta(j1)+B(i,j1)*sin(cta(i)-cta(j1); D(i)=D(i)+U(i)*U(j1)*(G(i,j1)*sin(cta(i)-cta(j1)-B(i,j1)*cos(cta(i)-cta(j1); end DP(t1)=P(i)-C(i); t1=t1+1; if A2(i,6)=2 DQ(t2)=Q(i

11、)-D(i); t2=t2+1; %計(jì)算失配功率 end endenddisp('失配功率: ') %輸出失配功率 disp(DP); disp(DQ); t1=t1-1;t2=t2-1; DPQ=DP'DQ' %求DP,DQ %= 以下為形成完整的雅克比矩陣 = for i=1:t1+t2 if abs(DPQ(i)>pr IT2=IT2+1; end end H=zeros(t1,t1);N=zeros(t1,t2);K=zeros(t2,t1);L=zeros(t2,t2); for i=1:t1 for j1=1:t1 if j1=isb&

12、j1=i H(i,j1)=0-U(i)*U(j1)*(G(i,j1)*sin(cta(i)-cta(j1)-B(i,j1)*cos(cta(i)-cta(j1); elseif j1=isb&j1=i H(i,j1)=U(i)2*B(i,j1)+D(i); end end end %計(jì)算H矩陣 for i=1:t1 for j1=1:t2 if j1=isb&j1=i N(i,j1)=0-U(i)*U(j1)*(G(i,j1)*cos(cta(i)-cta(j1)+B(i,j1)*sin(cta(i)-cta(j1); elseif j1=isb&j1=i N(i,j1

13、)=0-U(i)2*G(i,j1)-C(i); end end end %計(jì)算N矩陣 for i=1:t2 for j1=1:t1 if j1=isb&j1=i K(i,j1)= U(i)*U(j1)*(G(i,j1)*cos(cta(i)-cta(j1)+B(i,j1)*sin(cta(i)-cta(j1); elseif j1=isb&j1=i K(i,j1)=U(i)2*G(i,j1)-C(i); end end end %計(jì)算K矩陣 for i=1:t2 for j1=1:t2 if j1=isb&j1=i L(i,j1)=0-U(i)*U(j1)*(G(i,j

14、1)*sin(cta(i)-cta(j1)-B(i,j1)*cos(cta(i)-cta(j1); elseif j1=isb&j1=i L(i,j1)=U(i)2*B(i,j1)-D(i); end end end %計(jì)算L矩陣 J=H,N;K,L;%求雅可比矩陣 disp('【 雅可比矩陣 ：】'); disp(J); modify=-JDPQ; Dcta=modify(1:t1,:); t3=U(:,1:t2); DU=diag(t3,0)*modify(t1+1:t1+t2,:); t4=1; for i=1:t1 if A2(i,6)=1 cta(1,i)=c

15、ta(1,i)+Dcta(t4,1); t4=t4+1; end end t5=1; for i=1:t2 if A2(i,6)=2 U(1,i)=U(1,i)+DU(t5,1); t5=t5+1; end end ICT1=ICT1+1;end %修正原值for i=1:m UU(i)=U(i)*cos(cta(i)+1i*U(i)*sin(cta(i);endfor p=1:m c(p)=0; for q=1:m c(p)=c(p)+conj(Y(p,q)*conj(UU(q); end s(p)=UU(p)*c(p);enddisp('*');disp('【 各節(jié)

16、點(diǎn)電壓U為（節(jié)點(diǎn)從小到大排列）:】');disp(UU);disp('*');disp('【 各節(jié)點(diǎn)電壓相角為（節(jié)點(diǎn)從小到大排列） ：】');disp(180*angle(UU)/pi);disp('*');disp('【 按公式計(jì)算全部線路功率，結(jié)果如下:】');for i=1:m1 if A1(i,6)=0 p=A1(i,1);q=A1(i,2); else p=A1(i,2);q=A1(i,1); end Si(p,q)=UU(p)*(conj(UU(p)*conj(A1(i,4)./2)+(conj(UU(p)*A

17、1(i,5)-conj(UU(q)*conj(1./(A1(i,3)*A1(i,5);%各條支路首端功率Si f=p,q,Si(p,q); disp(f);endfor i=1:m1 if A1(i,6)=0 p=A1(i,1);q=A1(i,2); else p=A1(i,2);q=A1(i,1); end Sj(q,p)=UU(q)*(conj(UU(q)*conj(A1(i,4)./2)+(conj(UU(q)./A1(i,5)-conj(UU(p)*conj(1./(A1(i,3)*A1(i,5);%各條支路末端功率Sj f=q,p,Sj(q,p); disp(f);enddisp(&

18、#39;*');disp('【 各條支路的功率損耗DS為(順序同輸入A1時(shí)一樣):】');for i=1:m1 if A1(i,6)=0 p=A1(i,1);q=A1(i,2); else p=A1(i,2);q=A1(i,1); end DS(i)=Si(p,q)+Sj(q,p);%各條支路功率損耗DS disp(DS(i);endSp=0;for i=1:m Sp=Sp+UU(isb)*conj(Y(isb,i)*conj(UU(i);enddisp('【 平衡節(jié)點(diǎn)的功率：】');disp(Sp); 1.6輸出結(jié)果：>> dlxt【輸出結(jié)

19、果如下】：-【 已知量 ：】-【支路參數(shù)矩陣】1.0000 2.0000 0.0250 + 0.0800i 0 + 0.0700i 1.0000 0 1.0000 3.0000 0.0300 + 0.1000i 0 + 0.0900i 1.0000 0 2.0000 3.0000 0.0200 + 0.0600i 0 + 0.0500i 1.0000 0 4.0000 2.0000 0 + 0.1905i 0 1.0522 1.0000 5.0000 3.0000 0 + 0.1905i 0 1.0522 1.0000 【節(jié)點(diǎn)參數(shù)矩陣】0 -0.8055 - 0.5320i 1.0000 0

20、0 2.0000 0 -0.1800 - 0.1200i 1.0000 0 0 2.0000 0 0 1.0000 0 0 2.0000 0 0.5000 1.0000 1.0522 0 3.0000 0 0 1.0000 1.0522 0 1.0000 【導(dǎo)納矩陣】 Y=6.3110 -20.4822i -3.5587 +11.3879i -2.7523 + 9.1743i 0 0 -3.5587 +11.3879i 8.5587 -31.0693i -5.0000 +15.0000i 0 + 4.9889i 0 -2.7523 + 9.1743i -5.0000 +15.0000i 7.7

21、523 -28.8457i 0 0 + 4.9889i 0 0 + 4.9889i 0 0 - 5.2493i 0 0 0 0 + 4.9889i 0 0 - 5.2493i失配功率: 0.8055 0.1800 0 -0.5000 0.6120 0.4275 0.3175【 雅可比矩陣 ：】 -20.5622 11.3879 9.1743 0 -6.3110 3.5587 2.7523 11.3879 -31.3768 15.0000 4.9889 3.5587 -8.5587 5.0000 9.1743 15.0000 -29.1632 0 2.7523 5.0000 -7.7523 0

22、4.9889 0 -4.9889 0 0 0 6.3110 -3.5587 -2.7523 0 -20.4022 11.3879 9.1743 -3.5587 8.5587 -5.0000 0 11.3879 -30.7618 15.0000 -2.7523 -5.0000 7.7523 0 9.1743 15.0000 -28.5282失配功率: -0.2695 0.0231 0.0748 0.0799 -0.1030 -0.1277 -0.1137【 雅可比矩陣 ：】 -28.8103 15.9814 12.8289 0 -10.1478 4.5599 3.4378 15.7341 -41

23、.6855 20.1846 5.7668 5.3512 -11.7083 6.6232 12.6028 20.1216 -38.4668 0 4.1917 6.8122 -10.2937 0 5.7668 0 -5.7668 0 0.5799 0 7.9977 -4.5599 -3.4378 0 -30.0803 15.9814 12.8289 -5.3512 11.3945 -6.6232 0.5799 15.7341 -42.1809 20.1846 -4.1917 -6.8122 10.4433 0 12.6028 20.1216 -38.6942失配功率: -0.0107 0.0024

24、 0.0044 0.0015 -0.0018 -0.0041 -0.0042【 雅可比矩陣 ：】 -27.3240 15.1545 12.1695 0 -9.3998 4.4171 3.3502 14.9730 -39.8614 19.2510 5.6374 4.9977 -11.1924 6.3412 12.0040 19.2055 -36.8266 0 3.9019 6.4777 -9.8938 0 5.6374 0 -5.6374 0 0.5015 0 7.7673 -4.4171 -3.3502 0 -28.3915 15.1545 12.1695 -4.9977 10.8373 -6

25、.3412 0.5015 14.9730 -40.1097 19.2510 -3.9019 -6.4777 9.9027 0 12.0040 19.2055 -36.8350失配功率: 1.0e-004 * -0.1524 0.0448 0.0746 0.0085 1.0e-005 * -0.0652 -0.4685 -0.5571【 雅可比矩陣 ：】 -27.2778 15.1288 12.1490 0 -9.3743 4.4143 3.3490 14.9504 -39.8057 19.2220 5.6333 4.9855 -11.1784 6.3329 11.9862 19.1774 -3

26、6.7767 0 3.8915 6.4669 -9.8837 0 5.6333 0 -5.6333 0 0.5000 0 7.7633 -4.4143 -3.3490 0 -28.3418 15.1288 12.1490 -4.9855 10.8184 -6.3329 0.5000 14.9504 -40.0457 19.2220 -3.8915 -6.4669 9.8838 0 11.9862 19.1774 -36.7767*【 各節(jié)點(diǎn)電壓U為（節(jié)點(diǎn)從小到大排列）:】 1.1586 + 0.1243i 1.1292 + 0.0995i 1.1251 + 0.0951i 1.0000 - 0

27、.0007i 1.0000 *【 各節(jié)點(diǎn)電壓相角為（節(jié)點(diǎn)從小到大排列） ：】 6.1212 5.0334 4.8335 -0.0388 0*【 按公式計(jì)算全部線路功率，結(jié)果如下:】 1.0000 2.0000 0.4176 + 0.2856i 1.0000 3.0000 0.3879 + 0.2464i 2.0000 3.0000 0.0923 + 0.0216i 2.0000 4.0000 0.5000 + 1.1124i 3.0000 5.0000 0.4747 + 1.0795i 2.0000 1.0000 -0.4123 - 0.3613i 3.0000 1.0000 -0.3825

28、- 0.3468i 3.0000 2.0000 -0.0921 - 0.0851i 4.0000 2.0000 -0.5000 - 0.8919i 5.0000 3.0000 -0.4747 - 0.8717i*【 各條支路的功率損耗DS為(順序同輸入A1時(shí)一樣):】 0.0053 - 0.0757i 0.0054 - 0.1004i 0.0002 - 0.0635i 0 + 0.2205i -0.0000 + 0.2078i【 平衡節(jié)點(diǎn)的功率：】 -0.4747 - 0.3638i>> 7、結(jié)果截圖：2:P-Q解耦迭代法2.1 程序代碼cleardisp('【電力系統(tǒng)混合

29、坐標(biāo)下的PQ解耦迭代潮流計(jì)算】:');n=5;m=3;%= 參數(shù)初始化= Y=zeros(n,n); %導(dǎo)納矩陣U=ones(n,1); %電壓矢量a=zeros(n,1); %相角矢量Ps=zeros(n,1);Qs=zeros(n,1);P=zeros(n,1);Q=zeros(n,1);p=zeros(n-1,1);q=zeros(m,1);aa=zeros(n-1,1);u=zeros(m,1);k=0;%=節(jié)點(diǎn)導(dǎo)納矩陣的生成=I=sqrt(-1);%由于后面用到ij作為循環(huán)變量，故重定義虛數(shù)單位Y=6.311-I*20.3722 -2.7523+I*9.1743 -3.558

30、7+I*11.3879 0 0; -3.5587+I*11.3879 8.5587-I*30.9993 -5+I*15 I*4.9889 0; -2.7523+I*9.1743 -5+I*15 7.7523-I*28.7557 0 I*4.9889; 0 I*4.9889 0 -I*5.2493 0; 0 0 I*4.9889 0 -I*5.2493;%=初值的設(shè)定=U=1;1;1;1;1;a=0;0;0;0;0;Ps=-0.8055;-0.18;0;0.5;0;Qs=-0.532;-0.12;0;0;0;%求失配功率P(k)和Q(k)for i=1:n-1 s=0; for j=1:n s=

31、s+U(j)*(real(Y(i,j)*cos(a(i)-a(j)+imag(Y(i,j)*sin(a(i)-a(j); end P(i)=U(i)*s;endfor i=1:n-1 s=0; for j=1:n s=s+U(j)*(real(Y(i,j)*sin(a(i)-a(j)-imag(Y(i,j)*cos(a(i)-a(j); end Q(i)=U(i)*s;endfor i=1:n-1 detp(i)=Ps(i)-P(i);endfor i=1:m detq(i)=Qs(i)-Q(i);end %生成B'和B"矩陣 Bp=zeros(n-1,n-1); Bpp=z

32、eros(m,m); for i=1:n-1; for j=1:n-1; Bp(i,j)=imag(Y(i,j); end; end; for i=1:m; for j=1:m; Bpp(i,j)=imag(Y(i,j); end; end; k=0;while (max(detp)>=0.0001 | max(detq)>=0.0001) %設(shè)定精確度 aa=(-1)*inv(Bp)*detp(1)/U(1);detp(2)/U(2);detp(3)/U(3);detp(4)/U(4); u=(-1)*inv(Bpp)*detq(1)/U(1);detq(2)/U(2);detq(3)/U(3); %求解修正方程，得到修正量 aa,u for i=1:n-1 %計(jì)算 a(電壓相角),U(電壓幅值) a(i)=a(i)+aa(i); end for i=1:m U(i)=U(i)+u(i); end% detp'detq'迭代過(guò)程中失配功率變化情況 k=k+1;%(記錄迭代次數(shù)) for i=1:n-1 %計(jì)算功率誤