電力系統(tǒng)分析專業(yè)課程設(shè)計牛頓拉夫遜潮流計算

電力系統(tǒng)編程潮流計算1設(shè)計任務(wù)及初步分析1.1設(shè)計任務(wù)條件：節(jié)點數(shù)：3支路數(shù)：3計算精度：0.00010支路1：0.0300+j0.09001┠—————□—————┨2支路2：0.0200+j0.09002┠—————□—————┨3支路3：0.0300+j0.09003┠—————□—————┨1節(jié)點1：PQ節(jié)點，S（1）=-0.5000-j0.節(jié)點2：PQ節(jié)點，S（2）=-0.6000-j0.2500節(jié)點3：平衡節(jié)點，U（3）=1.0000∠0.0000規(guī)定：編寫程序計算潮流1.2初步分析潮流計算在數(shù)學(xué)上可歸結(jié)為求解非線性方程組，其數(shù)學(xué)模型簡寫如下：2牛頓-拉夫遜法簡介2.1概述 牛頓－拉夫遜法是當(dāng)前求解非線性方程最佳一種辦法。這種辦法特點就是把對非線性方程求解過程變成重復(fù)對相應(yīng)線性方程求解過程，普通稱為逐次線性化過程，就是牛頓－拉夫遜法核心。 牛頓-拉夫遜法基本原理是在解某一鄰域內(nèi)某一初始點出發(fā)，沿著該點一階偏導(dǎo)數(shù)——雅可比矩陣，朝減小方程殘差方向邁進一步，在新點上再計算殘差和雅可矩陣?yán)^續(xù)邁進，重復(fù)這一過程直到殘差達到收斂原則，即得到了非線性方程組解。由于越接近解，偏導(dǎo)數(shù)方向越準(zhǔn)，收斂速度也越快，因此牛頓法具備二階收斂特性。而所謂“某一鄰域”是指雅可比喻向均指向解范疇，否則也許走向非線性函數(shù)其他極值點，普通來說潮流由平電壓即各母線電壓(相角為0，幅值為1)啟動即在此鄰域內(nèi)。2．2普通概念 對于非線性代數(shù)方程組即(2－1) 在待求量某一種初始計算值附件，將上式展開泰勒級數(shù)并略去二階及以上高階項，得到如下線性化方程組(2－2)上式稱之為牛頓法修正方程式。由此可以求得第一次迭代修正量(2－3) 將和相加，得到變量第一次改進值。接著再從出發(fā)，重復(fù)上述計算過程。因而從一定初值出發(fā)，應(yīng)用牛頓法求解迭代格式為(2－4)(2－5)上兩式中：是函數(shù)對于變量一階偏導(dǎo)數(shù)矩陣，即雅可比矩陣；為迭代次數(shù)。 由式（2－4）和式子（2－5）可見，牛頓法核心便是重復(fù)形成求解修正方程式。牛頓法當(dāng)時始預(yù)計值和方程精準(zhǔn)解足夠接近時，收斂速度非常快，具備平方收斂特性。2.3潮流計算修正方程 運用牛頓－拉夫遜法計算潮流分布時，一方面要找出描述電力系統(tǒng)非線性方程。這里仍從節(jié)點電壓方程入手，設(shè)電力系統(tǒng)導(dǎo)納矩陣已知，則系統(tǒng)中某節(jié)點（節(jié)點）電壓方程為從而得 進而有 （2－6） 式（2－6）中，左邊第一項為給定節(jié)點注入功率，第二項為由節(jié)點電壓求得節(jié)點注入功率。她們兩者之差就是節(jié)點功率不平衡量。當(dāng)前有待解決問題就是各節(jié)點功率不平衡量都趨近于零時，各節(jié)點電壓應(yīng)具備價值。 由此可見，如將式（2－6）作為牛頓－拉夫遜中非線性函數(shù)，其中節(jié)點電壓就相稱于變量。建立了這種相應(yīng)關(guān)系，就可列出修正方程式，并迭代求解。但由于節(jié)點電壓可有兩種表達方式——以直角做表或者極坐標(biāo)表達，因而列出迭代方程相應(yīng)地也有兩種，下面分別討論。2.4直角坐標(biāo)表達修正方程 節(jié)點電壓以直角坐標(biāo)表達時，令、，且將導(dǎo)納矩陣中元素表達為，則式（2－7）變化為 （2－7） 再將實部和虛某些開，可得 （2－8）這就是直角坐標(biāo)下功率方程?？梢?，一種節(jié)點列出了有功和無功兩個方程。 對于節(jié)點（），給定量為節(jié)點注入功率，記為、，則由式（2－8）可得功率不平衡量，作為非線性方程 （2－9）式中、——分別表達第節(jié)點有功功率不平衡量和無功功率不平衡量。 對于節(jié)點（），給定量為節(jié)點注入有功功率及電壓數(shù)值，記為、，因而，可以運用有功功率不平衡量和電壓不平衡量表達出非線性方程，即有 （2－10）式中為電壓不平衡量。 對于平衡節(jié)點（），由于電壓數(shù)值及相位角給定，因此也擬定，不需要參加迭代求節(jié)點電壓。因而，對于個節(jié)點系統(tǒng)只能列出個方程，其中有功功率方程個，無功功率方程個，電壓方程個。將式（2－9）、式（2－10）非線性方程聯(lián)立，稱為個節(jié)點系統(tǒng)非線性方程組，且按泰勒級數(shù)在、（）展開，并略去高次項，得到以矩陣形式表達修正方程如下。（2－11）上式中雅可比矩陣各個元素則分別為 將（2－11）寫成縮寫形式 （2－12）對雅可比矩陣各元素可做如下討論：當(dāng)時，對于特定，只有該特定點和是變量，于是雅可比矩陣中各非對角元素表達為 當(dāng)時，雅可比矩陣中各對角元素表達式為 由上述表達式可知，直角坐標(biāo)雅可比矩陣有如下特點：雅可比矩陣是階方陣，由于、等等，因此它是一種不對稱方陣。雅可比矩陣中諸元素是節(jié)點電壓函數(shù)，在迭代過程中隨電壓變化而不斷地變化。雅可比矩陣非對角元素與節(jié)點導(dǎo)納矩陣中相應(yīng)非對角元素關(guān)于，當(dāng)中為零時，雅可比矩陣中相應(yīng)、、、也都為零，因而，雅可比矩陣也是一種稀疏矩陣。3程序設(shè)計3.1程序流程圖圖3-1程序流程圖3.2源程序n=3n1=3isb=3pr=0.0001B1=[120.03+0.09i010;130.03+0.09i010;230.02+0.09i010]B2=[0-0.5-0.2i1002;0-0.6-0.25i1002;001101]X=[10;20;30]%X=input('節(jié)點號和對地參數(shù):X=')；Y=zeros(n)；Times=1；%置迭代次數(shù)為初始值%創(chuàng)立節(jié)點導(dǎo)納矩陣fori=1:n1ifB1(i,6)==0%不含變壓器支路p=B1(i,1)；q=B1(i,2)；Y(p,q)=Y(p,q)-1/B1(i,3)；Y(q,p)=Y(p,q)；Y(p,p)=Y(p,p)+1/B1(i,3)+0.5*B1(i,4)；Y(q,q)=Y(q,q)+1/B1(i,3)+0.5*B1(i,4)；else%具有變壓器支路p=B1(i,1)；q=B1(i,2)；Y(p,q)=Y(p,q)-1/(B1(i,3)*B1(i,5))；Y(q,p)=Y(p,q)；Y(p,p)=Y(p,p)+1/B1(i,3)；Y(q,q)=Y(q,q)+1/(B1(i,5)^2*B1(i,3))；endendYOrgS=zeros(2*n-2,1)；DetaS=zeros(2*n-2,1)；%將OrgS、DetaS初始化%創(chuàng)立OrgS，用于存儲初始功率參數(shù)h=0；j=0；fori=1:n%對PQ節(jié)點解決ifi~=isb&B2(i,6)==2h=h+1；forj=1:nOrgS(2*h-1,1)=OrgS(2*h-1,1)+real(B2(i,3))*(real(Y(i,j))*real(B2(j,3))-imag(Y(i,j))*imag(B2(j,3)))+imag(B2(i,3))*(real(Y(i,j))*imag(B2(j,3))+imag(Y(i,j))*real(B2(j,3)))；OrgS(2*h,1)=OrgS(2*h,1)+imag(B2(i,3))*(real(Y(i,j))*real(B2(j,3))-imag(Y(i,j))*imag(B2(j,3)))-real(B2(i,3))*(real(Y(i,j))*imag(B2(j,3))+imag(Y(i,j))*real(B2(j,3)))；endendendfori=1:n%對PV節(jié)點解決，注意這時不可再將h初始化為0ifi~=isb&B2(i,6)==3h=h+1；forj=1:nOrgS(2*h-1,1)=OrgS(2*h-1,1)+real(B2(i,3))*(real(Y(i,j))*real(B2(j,3))-imag(Y(i,j))*imag(B2(j,3)))+imag(B2(i,3))*(real(Y(i,j))*imag(B2(j,3))+imag(Y(i,j))*real(B2(j,3)))；OrgS(2*h,1)=OrgS(2*h,1)+imag(B2(i,3))*(real(Y(i,j))*real(B2(j,3))-imag(Y(i,j))*imag(B2(j,3)))-real(B2(i,3))*(real(Y(i,j))*imag(B2(j,3))+imag(Y(i,j))*real(B2(j,3)))；endendendOrgS%創(chuàng)立PVU用于存儲PV節(jié)點初始電壓PVU=zeros(n-h-1,1)；t=0；fori=1:nifB2(i,6)==3t=t+1；PVU(t,1)=B2(i,3)；endendPVU%創(chuàng)立DetaS，用于存儲有功功率、無功功率和電壓幅值不平衡量h=0；fori=1:n%對PQ節(jié)點解決ifi~=isb&B2(i,6)==2h=h+1；DetaS(2*h-1,1)=real(B2(i,2))-OrgS(2*h-1,1)；DetaS(2*h,1)=imag(B2(i,2))-OrgS(2*h,1)；endendt=0；fori=1:n%對PV節(jié)點解決，注意這時不可再將h初始化為0ifi~=isb&B2(i,6)==3h=h+1；t=t+1；DetaS(2*h-1,1)=real(B2(i,2))-OrgS(2*h-1,1)；DetaS(2*h,1)=real(PVU(t,1))^2+imag(PVU(t,1))^2-real(B2(i,3))^2-imag(B2(i,3))^2；endendDetaS%創(chuàng)立I，用于存儲節(jié)點電流參數(shù)i=zeros(n-1,1)；h=0；fori=1:nifi~=isbh=h+1；I(h,1)=(OrgS(2*h-1,1)-OrgS(2*h,1)*sqrt(-1))/conj(B2(i,3))；endendI%創(chuàng)立Jacbi(雅可比矩陣)Jacbi=zeros(2*n-2)；h=0；k=0；fori=1:n%對PQ節(jié)點解決ifB2(i,6)==2h=h+1；forj=1:nifj~=isbk=k+1；ifi==j%對角元素解決Jacbi(2*h-1,2*k-1)=-imag(Y(i,j))*real(B2(i,3))+real(Y(i,j))*imag(B2(i,3))+imag(I(h,1))；Jacbi(2*h-1,2*k)=real(Y(i,j))*real(B2(i,3))+imag(Y(i,j))*imag(B2(i,3))+real(I(h,1))；Jacbi(2*h,2*k-1)=-Jacbi(2*h-1,2*k)+2*real(I(h,1))；Jacbi(2*h,2*k)=Jacbi(2*h-1,2*k-1)-2*imag(I(h,1))；else%非對角元素解決Jacbi(2*h-1,2*k-1)=-imag(Y(i,j))*real(B2(i,3))+real(Y(i,j))*imag(B2(i,3))；Jacbi(2*h-1,2*k)=real(Y(i,j))*real(B2(i,3))+imag(Y(i,j))*imag(B2(i,3))；Jacbi(2*h,2*k-1)=-Jacbi(2*h-1,2*k)；Jacbi(2*h,2*k)=Jacbi(2*h-1,2*k-1)；endifk==(n-1)%將用于內(nèi)循環(huán)指針置于初始值，以保證雅可比矩陣換行k=0；endendendendendk=0；fori=1:n%對PV節(jié)點解決ifB2(i,6)==3h=h+1；forj=1:nifj~=isbk=k+1；ifi==j%對角元素解決Jacbi(2*h-1,2*k-1)=-imag(Y(i,j))*real(B2(i,3))+real(Y(i,j))*imag(B2(i,3))+imag(I(h,1))；Jacbi(2*h-1,2*k)=real(Y(i,j))*real(B2(i,3))+imag(Y(i,j))*imag(B2(i,3))+real(I(h,1))；Jacbi(2*h,2*k-1)=2*imag(B2(i,3))；Jacbi(2*h,2*k)=2*real(B2(i,3))；else%非對角元素解決Jacbi(2*h-1,2*k-1)=-imag(Y(i,j))*real(B2(i,3))+real(Y(i,j))*imag(B2(i,3))；Jacbi(2*h-1,2*k)=real(Y(i,j))*real(B2(i,3))+imag(Y(i,j))*imag(B2(i,3))；Jacbi(2*h,2*k-1)=0；Jacbi(2*h,2*k)=0；endifk==(n-1)%將用于內(nèi)循環(huán)指針置于初始值，以保證雅可比矩陣換行k=0；endendendendendJacbi%求解修正方程，獲取節(jié)點電壓不平衡量DetaU=zeros(2*n-2,1)；DetaU=inv(Jacbi)*DetaS；DetaU%修正節(jié)點電壓j=0；fori=1:n%對PQ節(jié)點解決ifB2(i,6)==2j=j+1；B2(i,3)=B2(i,3)+DetaU(2*j,1)+DetaU(2*j-1,1)*sqrt(-1)；endendfori=1:n%對PV節(jié)點解決ifB2(i,6)==3j=j+1；B2(i,3)=B2(i,3)+DetaU(2*j,1)+DetaU(2*j-1,1)*sqrt(-1)；endendB2whileabs(max(DetaS))>prOrgS=zeros(2*n-2,1)；%!!!初始功率參數(shù)在迭代過程中是不累加，因此在這里必要將其初始化為零矩陣h=0；j=0；fori=1:nifi~=isb&B2(i,6)==2h=h+1；forj=1:nOrgS(2*h-1,1)=OrgS(2*h-1,1)+real(B2(i,3))*(real(Y(i,j))*real(B2(j,3))-imag(Y(i,j))*imag(B2(j,3)))+imag(B2(i,3))*(real(Y(i,j))*imag(B2(j,3))+imag(Y(i,j))*real(B2(j,3)))；OrgS(2*h,1)=OrgS(2*h,1)+imag(B2(i,3))*(real(Y(i,j))*real(B2(j,3))-imag(Y(i,j))*imag(B2(j,3)))-real(B2(i,3))*(real(Y(i,j))*imag(B2(j,3))+imag(Y(i,j))*real(B2(j,3)))；endendendfori=1:nifi~=isb&B2(i,6)==3h=h+1；forj=1:nOrgS(2*h-1,1)=OrgS(2*h-1,1)+real(B2(i,3))*(real(Y(i,j))*real(B2(j,3))-imag(Y(i,j))*imag(B2(j,3)))+imag(B2(i,3))*(real(Y(i,j))*imag(B2(j,3))+imag(Y(i,j))*real(B2(j,3)))；OrgS(2*h,1)=OrgS(2*h,1)+imag(B2(i,3))*(real(Y(i,j))*real(B2(j,3))-imag(Y(i,j))*imag(B2(j,3)))-real(B2(i,3))*(real(Y(i,j))*imag(B2(j,3))+imag(Y(i,j))*real(B2(j,3)))；endendendOrgS%創(chuàng)立DetaSh=0；fori=1:nifi~=isb&B2(i,6)==2h=h+1；DetaS(2*h-1,1)=real(B2(i,2))-OrgS(2*h-1,1)；DetaS(2*h,1)=imag(B2(i,2))-OrgS(2*h,1)；endendt=0；fori=1:nifi~=isb&B2(i,6)==3h=h+1；t=t+1；DetaS(2*h-1,1)=real(B2(i,2))-OrgS(2*h-1,1)；DetaS(2*h,1)=real(PVU(t,1))^2+imag(PVU(t,1))^2-real(B2(i,3))^2-imag(B2(i,3))^2；endendDetaS%創(chuàng)立Ii=zeros(n-1,1)；h=0；fori=1:nifi~=isbh=h+1；I(h,1)=(OrgS(2*h-1,1)-OrgS(2*h,1)*sqrt(-1))/conj(B2(i,3))；endendI%創(chuàng)立JacbiJacbi=zeros(2*n-2)；h=0；k=0；fori=1:nifB2(i,6)==2h=h+1；forj=1:nifj~=isbk=k+1；ifi==jJacbi(2*h-1,2*k-1)=-imag(Y(i,j))*real(B2(i,3))+real(Y(i,j))*imag(B2(i,3))+imag(I(h,1))；Jacbi(2*h-1,2*k)=real(Y(i,j))*real(B2(i,3))+imag(Y(i,j))*imag(B2(i,3))+real(I(h,1))；Jacbi(2*h,2*k-1)=-Jacbi(2*h-1,2*k)+2*real(I(h,1))；Jacbi(2*h,2*k)=Jacbi(2*h-1,2*k-1)-2*imag(I(h,1))；elseJacbi(2*h-1,2*k-1)=-imag(Y(i,j))*real(B2(i,3))+real(Y(i,j))*imag(B2(i,3))；Jacbi(2*h-1,2*k)=real(Y(i,j))*real(B2(i,3))+imag(Y(i,j))*imag(B2(i,3))；Jacbi(2*h,2*k-1)=-Jacbi(2*h-1,2*k)；Jacbi(2*h,2*k)=Jacbi(2*h-1,2*k-1)；endifk==(n-1)k=0；endendendendendk=0；fori=1:nifB2(i,6)==3h=h+1；forj=1:nifj~=isbk=k+1；ifi==jJacbi(2*h-1,2*k-1)=-imag(Y(i,j))*real(B2(i,3))+real(Y(i,j))*imag(B2(i,3))+imag(I(h,1))；Jacbi(2*h-1,2*k)=real(Y(i,j))*real(B2(i,3))+imag(Y(i,j))*imag(B2(i,3))+real(I(h,1))；Jacbi(2*h,2*k-1)=2*imag(B2(i,3))；Jacbi(2*h,2*k)=2*real(B2(i,3))；elseJacbi(2*h-1,2*k-1)=-imag(Y(i,j))*real(B2(i,3))+real(Y(i,j))*imag(B2(i,3))；Jacbi(2*h-1,2*k)=real(Y(i,j))*real(B2(i,3))+imag(Y(i,j))*imag(B2(i,3))；Jacbi(2*h,2*k-1)=0；Jacbi(2*h,2*k)=0；endifk==(n-1)k=0；endendendendendJacbiDetaU=zeros(2*n-2,1)；DetaU=inv(Jacbi)*DetaS；DetaU%修正節(jié)點電壓j=0；fori=1:nifB2(i,6)==2j=j+1；B2(i,3)=B2(i,3)+DetaU(2*j,1)+DetaU(2*j-1,1)*sqrt(-1)；endendfori=1:nifB2(i,6)==3j=j+1；B2(i,3)=B2(i,3)+DetaU(2*j,1)+DetaU(2*j-1,1)*sqrt(-1)；endendB2Times=Times+1；%迭代次數(shù)加1enddisp('迭代次數(shù)為');disp(Times)forno=1:nVn(no)=B2(no,3);str1=sprintf('節(jié)點%d電壓是%d',no);disp(str);disp(Vn(no));end3.3成果及分析節(jié)點導(dǎo)納矩陣為Y=迭代運營了3次，每次功率不平衡量表3-1所示表3-1迭代過程中各節(jié)點功率不平衡量k?P?P0-0.5-j0.2-0.6-j0.251-0.0083-j0.0279-0.0091-j0.03602-0.0000297-j0.0000955-0.0000326-j0.000133230.0000+j0.00000.0000+j0.0000表3-2迭代過程中各節(jié)點電壓kee10.9664-j0.04230.9658-j0.045720.9632-j0.04230.9623-j0.045730.9632-j0.04230.9623-j0.0457