電力系統(tǒng)分析潮流計算matlab

1、目錄：1、 軟件需求說明書.32、 概要設(shè)計說明書.41、編寫潮流計算程序.4 2、數(shù)據(jù)的輸入測試.4 3、運行得出結(jié)果.4 4、進(jìn)行實驗結(jié)果驗證.43、 詳細(xì)設(shè)計說明書.51、數(shù)據(jù)導(dǎo)入模塊.52、節(jié)點導(dǎo)納矩陣模塊.53、編號判斷模塊.54、收斂條件判定模塊.55、雅可比矩陣模塊.56、迭代計算模塊.57、計算輸出參數(shù)模塊.54、 程序代碼.65、 最測試?yán)?151、輸入結(jié)果.152、輸出結(jié)果.153、結(jié)果驗證.151、 軟件需求說明書 本次設(shè)計利用MATLAB/C+/C（使用MATLAB）編程工具編寫潮流計算，實現(xiàn)對節(jié)點電壓和功率分布的求取。潮流方程的求解基本方法是迭代，包括牛頓-拉夫遜法，

2、以及P-Q分解法，本次設(shè)計采用牛頓迭代法。牛頓迭代法（Newton's method）又稱為牛頓-拉夫遜方法，它是牛頓在17世紀(jì)提出的一種在實數(shù)域和復(fù)數(shù)域上近似求解方程的方法。多數(shù)方程不存在求根公式，因此求精確根非常困難，甚至不可能，從而尋找方程的近似根就顯得特別重要。方法使用函數(shù)f(x)的泰勒級數(shù)的前面幾項來尋找方程f(x) = 0的根。牛頓迭代法是求方程根的重要方法之一，其最大優(yōu)點是在方程f(x) = 0的單根附近具有平方收斂，而且該法還可以用來求方程的重根、復(fù)根。 牛頓迭代法是取x0 之后，在這個基礎(chǔ)上，找到比x0 更接

3、近的方程的跟，一步一步迭代，從而找到更接近方程根的近似跟。牛頓迭代法是求方程根的重要方法之一，其最大優(yōu)點是在方程f(x) = 0 的單根附近具有平方收斂，而且該法還可以用來求方程的重根、復(fù)根。電力系統(tǒng)潮流計算，一般來說，各個母線所供負(fù)荷的功率是已知的，各個節(jié)點電壓是未知的（平衡節(jié)點外）可以根據(jù)網(wǎng)絡(luò)結(jié)構(gòu)形成節(jié)點導(dǎo)納矩陣，然后由節(jié)點導(dǎo)納矩陣列寫功率方程，由于功率方程里功率是已知的，電壓的幅值和相角是未知的，這樣潮流計算的問題就轉(zhuǎn)化為求解非線性方程組的問題了。為了便于用迭代法解方程組，需要將上述功率方程改寫成功率平衡方程，并對功率平衡方程求偏導(dǎo)，得出對應(yīng)的雅可比矩陣，給未知節(jié)點賦電壓初值，一般為額定

4、電壓，將初值帶入功率平衡方程，得到功率不平衡量，這樣由功率不平衡量、雅可比矩陣、節(jié)點電壓不平衡量（未知的）構(gòu)成了誤差方程，解誤差方程，得到節(jié)點電壓不平衡量，節(jié)點電壓加上節(jié)點電壓不平衡量構(gòu)成新的節(jié)點電壓初值，將新的初值帶入原來的功率平衡方程，并重新形成雅可比矩陣，然后計算新的電壓不平衡量，這樣不斷迭代，不斷修正，一般迭代三到五次就能收斂。2、 概要設(shè)計說明書1、編寫潮流計算程序本程序主要分為七個模塊：數(shù)據(jù)導(dǎo)入模塊、節(jié)點導(dǎo)納矩陣模塊、編號判斷模塊、收斂條件判定模塊、雅可比矩陣模塊、迭代計算模塊、計算輸出參數(shù)模塊。下圖為潮流迭代框圖。2、 數(shù)據(jù)的輸入測試 本次設(shè)計是將電力系統(tǒng)分析課本P88例題進(jìn)行潮

5、流計算。3、 運行得出結(jié)果 得出電壓、有功功率、無功功率、角度的數(shù)據(jù)。4、進(jìn)行實驗結(jié)果驗證3、 詳細(xì)設(shè)計說明書 1、數(shù)據(jù)導(dǎo)入模塊： 利用Excel輸入已知節(jié)點、支路數(shù)據(jù)，通過“讀取”將數(shù)據(jù)導(dǎo)入MATLAB中。 2、節(jié)點導(dǎo)納矩陣模塊： 利用已知的電阻電抗及導(dǎo)納的值，根據(jù)導(dǎo)納的計算公式，計算出節(jié)點的自導(dǎo)納及節(jié)點間的互導(dǎo)納的值，按照節(jié)點編號組成導(dǎo)納矩陣，利用MATLAB“real”和“imag”調(diào)用導(dǎo)納矩陣中的實部和虛部，分別形成電導(dǎo)和電納的矩陣。 3、編號判斷模塊： 當(dāng)首節(jié)點在變壓器左側(cè)，設(shè)為1，位于變壓器右側(cè)，設(shè)為2，既非1也非0為不含變壓器； 節(jié)點類型為PQ時，為1，節(jié)點類型為PV時，為2，節(jié)

6、點類型為V時，為3。 4、收斂條件判定模塊： 根據(jù)節(jié)點的類型賦初值，并進(jìn)行失配功率的初步計算，判斷是否符合收斂條件maxPi，Qi<。如不符合，則進(jìn)行后續(xù)的計算。 5、雅可比矩陣模塊： 根據(jù)節(jié)點類型確定雅克比矩陣的階數(shù)，然后根據(jù)n維非線性方程組的修正方程求出雅克比矩陣。 6、迭代計算模塊： 解修正方程，并進(jìn)行收斂條件判斷，如不符合條件則進(jìn)行下一次迭代，以一直到符合條件為止。 7、計算輸出參數(shù)模塊： 當(dāng)滿足收斂條件maxPi，Qi<時，結(jié)束迭代計算，計算電壓、有功功率、無功功率、電壓角度等參數(shù)。4、 程序代碼clear %清除變量；filename='E:MATLABzhil

7、u.xlsx'a=xlsread(filename); %讀取支路信息；filename='E:MATLABjiedian.xlsx' c=xlsread(filename); %讀取節(jié)點信息 b=zeros(a(1,7); %定義節(jié)點導(dǎo)納矩陣G=zeros(a(1,7);B=zeros(a(1,7);for i1=1:(a(1,7) if (a(i1,5)=0) %首節(jié)點在變壓器左側(cè) b(a(i1,1),a(i1,1)=b(a(i1,1),a(i1,1)+1/(a(i1,3)+1j*a(i1,4); b(a(i1,2),a(i1,2)=b(a(i1,2),a(i1,2

8、)+1/(a(i1,3)+1j*a(i1,4)*a(i1,6)2); b(a(i1,1),a(i1,2)=b(a(i1,1),a(i1,2)+(-1/(a(i1,3)+1j*a(i1,4)*a(i1,6); b(a(i1,2),a(i1,1)=b(a(i1,2),a(i1,1)+(-1/(a(i1,3)+1j*a(i1,4)*a(i1,6);%進(jìn)行導(dǎo)納計算 G(a(i1,1),a(i1,1)=real(b(a(i1,1),a(i1,1); G(a(i1,2),a(i1,2)=real(b(a(i1,2),a(i1,2); G(a(i1,1),a(i1,2)=real(b(a(i1,1),a(i

9、1,2); G(a(i1,2),a(i1,1)=real(b(a(i1,2),a(i1,1); B(a(i1,1),a(i1,1)=imag(b(a(i1,1),a(i1,1); B(a(i1,2),a(i1,2)=imag(b(a(i1,2),a(i1,2); B(a(i1,1),a(i1,2)=imag(b(a(i1,1),a(i1,2); B(a(i1,2),a(i1,1)=imag(b(a(i1,2),a(i1,1); elseif (a(i1,5)=1) %首節(jié)點在變壓器右側(cè) b(a(i1,1),a(i1,1)=b(a(i1,1),a(i1,1)+1/(a(i1,3)+1j*a(i1

10、,4)*a(i1,6)2); b(a(i1,2),a(i1,2)=b(a(i1,2),a(i1,2)+1/(a(i1,3)+1j*a(i1,4); b(a(i1,1),a(i1,2)=b(a(i1,1),a(i1,2)+(-1/(a(i1,3)+1j*a(i1,4)*a(i1,6); b(a(i1,2),a(i1,1)=b(a(i1,2),a(i1,1)+(-1/(a(i1,3)+1j*a(i1,4)*a(i1,6); %進(jìn)行導(dǎo)納計算 G(a(i1,1),a(i1,1)=real(b(a(i1,1),a(i1,1); G(a(i1,2),a(i1,2)=real(b(a(i1,2),a(i1,

11、2); G(a(i1,1),a(i1,2)=real(b(a(i1,1),a(i1,2); G(a(i1,2),a(i1,1)=real(b(a(i1,2),a(i1,1); B(a(i1,1),a(i1,1)=imag(b(a(i1,1),a(i1,1); B(a(i1,2),a(i1,2)=imag(b(a(i1,2),a(i1,2); B(a(i1,1),a(i1,2)=imag(b(a(i1,1),a(i1,2); B(a(i1,2),a(i1,1)=imag(b(a(i1,2),a(i1,1); else %不含變壓器支路 b(a(i1,1),a(i1,1)=b(a(i1,1),a(

12、i1,1)+1/(a(i1,3)+1j*a(i1,4)+1j*a(i1,5); b(a(i1,2),a(i1,2)=b(a(i1,2),a(i1,2)+1/(a(i1,3)+1j*a(i1,4)+1j*a(i1,5); b(a(i1,1),a(i1,2)=b(a(i1,1),a(i1,2)+(-1/(a(i1,3)+1j*a(i1,4); b(a(i1,2),a(i1,1)=b(a(i1,2),a(i1,1)+(-1/(a(i1,3)+1j*a(i1,4); %進(jìn)行導(dǎo)納計算 G(a(i1,1),a(i1,1)=real(b(a(i1,1),a(i1,1); G(a(i1,2),a(i1,2)=

13、real(b(a(i1,2),a(i1,2); G(a(i1,1),a(i1,2)=real(b(a(i1,1),a(i1,2); G(a(i1,2),a(i1,1)=real(b(a(i1,2),a(i1,1); B(a(i1,1),a(i1,1)=imag(b(a(i1,1),a(i1,1); B(a(i1,2),a(i1,2)=imag(b(a(i1,2),a(i1,2); B(a(i1,1),a(i1,2)=imag(b(a(i1,1),a(i1,2); B(a(i1,2),a(i1,1)=imag(b(a(i1,2),a(i1,1); endend %計數(shù)各個節(jié)點個數(shù)PQjd=0;P

14、Vjd=0;VOjd=0;for i2=1:a(1,7) if c(i2,6)=1 PQjd=PQjd+1; end if c(i2,6)=2 PVjd=PVjd+1; end if c(i2,6)=3 VOjd=VOjd+1; end end%注入功率的計算zhuru_PQ=zeros(PQjd*2+PVjd,1);k=1;for i3=1:a(1,7) if c(i3,6)=1 zhuru_PQ(k,1)=c(i3,2)-c(i3,4); zhuru_PQ(k+a(1,7)-1,1)=0; for i10=1:a(1,7) zhuru_PQ(k+a(1,7)-1,1)=zhuru_PQ(k

15、+a(1,7)-1,1)+B(i3,i10); end zhuru_PQ(k+a(1,7)-1,1)=zhuru_PQ(k+a(1,7)-1,1)-c(i3,5) k=k+1; end if c(i3,6)=2 zhuru_PQ(k,1)=c(i3,2)-c(i3,4); k=k+1; endendn=a(1,7);%后續(xù)矩陣變量定義angle_u=zeros(PQjd*2+PVjd,1); %迭代角度、電壓矩陣angle_u(n:PQjd*2+PVjd,1)=1; %電壓初值設(shè)為1，角度初值設(shè)為0delta_au=zeros(PQjd*2+PVjd,1);delta_au(n:PQjd*2+

16、PVjd,1)=0;%角度、電壓修正量向量%失配功率的計算 k=1; k_=PQjd+1; delta_PQ= zeros(PQjd*2+PVjd,1); for i4=1:n%PQ節(jié)點的計算 if c(i4,6)=1 delta_PQ(k,1)=zhuru_PQ(k,1); delta_PQ(k+n-1,1)=zhuru_PQ(k+n-1,1); i6=1; for i5=1:n if c(i5,6)=1 delta_PQ(k,1)=delta_PQ(k,1)-angle_u(n-1+k)*angle_u(n-1+i6)*(G(i4,i5)*cos(angle_u(k,1)-angle_u(

17、i6,1)+B(i4,i5)*sin(angle_u(k,1)-angle_u(i6,1); delta_PQ(k+n-1,1)=delta_PQ(k+n-1,1)-angle_u(n-1+k)*angle_u(n-1+i6)*(G(i4,i5)*sin(angle_u(k,1)-angle_u(i6,1)-B(i4,i5)*cos(angle_u(k,1)-angle_u(i6,1); i6=i6+1; end if c(i5,6)=2 delta_PQ(k,1)=delta_PQ(k,1)-angle_u(n-1+k)*c(i5,1)*(G(i4,i5)*cos(angle_u(k,1)-

18、angle_u(i6,1)+B(i4,i5)*sin(angle_u(k,1)-angle_u(i6,1); delta_PQ(k+n-1,1)=delta_PQ(k+n-1,1)-angle_u(n-1+k)*c(i5,1)*(G(i4,i5)*sin(angle_u(k,1)-angle_u(i6,1)-B(i4,i5)*cos(angle_u(k,1)-angle_u(i6,1); i6=i6+1; end if c(i5,6)=3 delta_PQ(k,1)=delta_PQ(k,1)-angle_u(n-1+k)*c(i5,1)*(G(i4,i5)*cos(angle_u(k,1)+

19、B(i4,i5)*sin(angle_u(k,1); delta_PQ(k+n-1,1)=delta_PQ(k+n-1,1)-angle_u(n-1+k)*c(i5,1)*(G(i4,i5)*sin(angle_u(k,1)-B(i4,i5)*cos(angle_u(k,1); end end k=k+1; end %PV節(jié)點的計算 if c(i4,6)=2 delta_PQ(k_,1)=zhuru_PQ(k_,1); i6=1; for i5=1:n if c(i5,6)=1 delta_PQ(k_,1)=delta_PQ(k_,1)-c(i4,1)*angle_u(n-1+i6)*(G(i

20、4,i5)*cos(angle_u(k_,1)-angle_u(i6,1)+B(i4,i5)*sin(angle_u(k_,1)-angle_u(i6,1); i6=i6+1; end if c(i5,6)=2 delta_PQ(k_,1)=delta_PQ(k_,1)-c(i4,1)*c(i5,1)*(G(i4,i5)*cos(angle_u(k_,1)-angle_u(i6,1)+B(i4,i5)*sin(angle_u(k_,1)-angle_u(i6,1); i6=i6+1; end if c(i5,6)=3 delta_PQ(k_,1)=delta_PQ(k_,1)-c(i4,1)*

21、c(i5,1)*(G(i4,i5)*cos(angle_u(k_,1)+B(i4,i5)*sin(angle_u(k_,1); end end k_=k_+1; end endzhongjian_delta_PQ=abs(delta_PQ); %取絕對值jingdu=max(zhongjian_delta_PQ); %為迭代精度的判斷做賦值準(zhǔn)備 while (jingdu>0.0001) %雅各比矩陣的形成 J=zeros(2*PQjd+PVjd); %矩陣定義 U=ones(n,1); %為方便運算將電壓取出放入一個新的矩陣 for i=1:PQjd U(i,1)=angle_u(i+

22、n-1,1); end angle=zeros(n,1); %為方便運算將角度取出放入一個新的矩陣 for i=1:n-1 angle(i,1)=angle_u(i,1); end %H% %H部分生成 for i1=1:n-1 for i2=1:n-1 if i1=i2 J(i1,i2)=-U(i1,1)*U(i2,1)*(G(i1,i2)*sin(angle(i1,1)-angle(i2,1)-B(i1,i2)*cos(angle(i1,1)-angle(i2,1); end if i1=i2 for i3=1:n if i3=i1 J(i1,i2)=J(i1,i2)+U(i1,1)*U(

23、i3,1)*(G(i1,i3)*sin(angle(i1,1)-angle(i3,1)-B(i1,i3)*cos(angle(i1,1)-angle(i3,1); end end end end end %N% %N部分生成 for i1=1:n-1 for i2=1:PQjd if i1=i2 J(i1,i2+n-1)=-U(i1,1)*(G(i1,i2)*cos(angle(i1,1)-angle(i2,1)-B(i1,i2)*sin(angle(i1,1)-angle(i2,1); end if i1=i2 J(i1,i2+n-1)=(-2)*U(i1,1)*U(i1,1)*G(i1,i

24、1); for i3=1:n if i3=i1 J(i1,i2+n-1)=J(i1,i2+n-1)-U(i3,1)*(G(i1,i3)*cos(angle(i1,1)-angle(i3,1)+B(i1,i3)*sin(angle(i1,1)-angle(i3,1); end end end end end %K% %K部分生成 for i1=1:PQjd for i2=1:n-1 if i1=i2 J(i1+n-1,i2)=U(i1,1)*U(i2,1)*(G(i1,i2)*cos(angle(i1,1)-angle(i2,1)-B(i1,i2)*sin(angle(i1,1)-angle(i

25、2,1); end if i1=i2 for i3=1:n if i3=i1 J(i1+n-1,i2)=J(i1+n-1,i2)-U(i1,1)*U(i3,1)*(G(i1,i3)*cos(angle(i1,1)-angle(i3,1)+B(i1,i3)*sin(angle(i1,1)-angle(i3,1); end end end end end %L% %L部分生成 for i1=1:PQjd for i2=1:PQjd if i1=i2 J(i1+n-1,i2+n-1)=-U(i1,1)*(G(i1,i2)*sin(angle(i1,1)-angle(i2,1)-B(i1,i2)*co

26、s(angle(i1,1)-angle(i2,1); end if i1=i2 J(i1+n-1,i2+n-1)=2*U(i1,1)*U(i1,1)*B(i1,i1); for i3=1:n if i3=i1 J(i1+n-1,i2+n-1)=J(i1+n-1,i2+n-1)-U(i3,1)*(G(i1,i3)*sin(angle(i1,1)-angle(i3,1)-B(i1,i3)*cos(angle(i1,1)-angle(i3,1); end end end end end delta_au=-inv(J)*delta_PQ; %解修正方程angle_u=angle_u+delta_au

27、; %形成新的初始值%再次計算失配功率 k=1; k_=PQjd+1; delta_PQ= zeros(PQjd*2+PVjd,1); for i4=1:n% if c(i4,6)=1 delta_PQ(k,1)=zhuru_PQ(k,1); delta_PQ(k+n-1,1)=zhuru_PQ(k+n-1,1); i6=1; for i5=1:n if c(i5,6)=1 delta_PQ(k,1)=delta_PQ(k,1)-angle_u(n-1+k)*angle_u(n-1+i6)*(G(i4,i5)*cos(angle_u(k,1)-angle_u(i6,1)+B(i4,i5)*si

28、n(angle_u(k,1)-angle_u(i6,1); delta_PQ(k+n-1,1)=delta_PQ(k+n-1,1)-angle_u(n-1+k)*angle_u(n-1+i6)*(G(i4,i5)*sin(angle_u(k,1)-angle_u(i6,1)-B(i4,i5)*cos(angle_u(k,1)-angle_u(i6,1); i6=i6+1; end if c(i5,6)=2 delta_PQ(k,1)=delta_PQ(k,1)-angle_u(n-1+k)*c(i5,1)*(G(i4,i5)*cos(angle_u(k,1)-angle_u(i6,1)+B(i

29、4,i5)*sin(angle_u(k,1)-angle_u(i6,1); delta_PQ(k+n-1,1)=delta_PQ(k+n-1,1)-angle_u(n-1+k)*c(i5,1)*(G(i4,i5)*sin(angle_u(k,1)-angle_u(i6,1)-B(i4,i5)*cos(angle_u(k,1)-angle_u(i6,1); i6=i6+1; end if c(i5,6)=3 delta_PQ(k,1)=delta_PQ(k,1)-angle_u(n-1+k)*c(i5,1)*(G(i4,i5)*cos(angle_u(k,1)+B(i4,i5)*sin(angl

30、e_u(k,1); delta_PQ(k+n-1,1)=delta_PQ(k+n-1,1)-angle_u(n-1+k)*c(i5,1)*(G(i4,i5)*sin(angle_u(k,1)-B(i4,i5)*cos(angle_u(k,1); end end k=k+1; end % if c(i4,6)=2 delta_PQ(k_,1)=zhuru_PQ(k_,1); i6=1; for i5=1:n if c(i5,6)=1 delta_PQ(k_,1)=delta_PQ(k_,1)-c(i4,1)*angle_u(n-1+i6)*(G(i4,i5)*cos(angle_u(k_,1)-angle_u(i6,1)+B(i4,i5)*sin(angle_u(k_,1)-angl