【正文】 程的振蕩或者不收斂，從數(shù)學(xué)上講，非線性的潮流計(jì)算方程組本來就是無解的。另外，為了優(yōu)化系統(tǒng)的運(yùn)行，從所有以上的可行潮流解中挑選出滿足一定指標(biāo)要求的一個(gè)最佳方案就是最優(yōu)潮流問題。四、本文主要工作 本文主要工作主要是詳細(xì)介紹牛頓——拉夫遜法的原理、算法設(shè)計(jì)和Matlab程序的編寫。步驟是：（1）設(shè)一組結(jié)點(diǎn)電壓；（2）求功率和電壓的不平衡量；（3）求雅克比矩陣的各個(gè)元素；（4）解修正方程式。但是牛頓——拉夫遜法對于初值有比較高的要求，當(dāng)給定的初值與精確值相差較大時(shí)，計(jì)算結(jié)果會產(chǎn)生很大的誤差，甚至不能收斂。在利用計(jì)算機(jī)的復(fù)雜電力系統(tǒng)的潮流計(jì)算中用的最多的是節(jié)點(diǎn)電壓方程。將其用矩陣形式表式為： ()式（）可以展開為 ()在式（）中是節(jié)點(diǎn)注入電流的列向量。是一個(gè)階的節(jié)點(diǎn)導(dǎo)納矩陣，其就等于網(wǎng)絡(luò)中出參考地之外的節(jié)點(diǎn)數(shù)。（一）輸電線路的等值電路。當(dāng)線路長度在100公里到300公里之間時(shí)，要考慮線路的導(dǎo)納，導(dǎo)納用集中參數(shù)表示。理想變壓器的變比是實(shí)際的變壓器變比，而阻抗一般是折算到低壓側(cè)的變壓器短路阻抗，串聯(lián)在理想變壓器的低壓側(cè)。節(jié)點(diǎn)導(dǎo)納矩陣求取時(shí)，注意一下幾點(diǎn)。與節(jié)點(diǎn)2對應(yīng)的對角元。但在某些情況下，迭代法也適用于線性方程，下面以一個(gè)線性方程為例介紹迭代法。然后再將，代入上式，解得又一組解，經(jīng)過三次迭代后可得，誤差已經(jīng)減小到了千分之一。但是由于在工程實(shí)踐中，通常我們已知的既不是節(jié)點(diǎn)電壓，也不是節(jié)點(diǎn)電流，而是節(jié)點(diǎn)功率，實(shí)際計(jì)算中幾乎無一例外的要迭代解非線性的節(jié)點(diǎn)電壓方程（*表示復(fù)數(shù)的共軛）。由式（）可見，在功率方程中，母線電壓的相位角是以差=的形式出現(xiàn)的，亦即決定功率大小的是相對相位角或相對功率角，而不是絕對相位角或絕對功率角。它們就稱為不可控變量或擾動變量?？刂谱兞坑昧邢蛄勘硎?。狀態(tài)變量一般用列向量表示。看來似乎將變量做如上分類后，只要已知或給定擾動變量和控制變量就可以運(yùn)用功率方程式（）解出狀態(tài)變量。這一對控制變量、將使系統(tǒng)的功率保持平衡。這樣，原則上可以從個(gè)方程中解出個(gè)未知量。對于某些狀態(tài)變量還有如下的約束條件這條件主要是保證系統(tǒng)運(yùn)行的穩(wěn)定性所要求的。這樣，系統(tǒng)的節(jié)點(diǎn)就因給定變量的不同而分為三類。屬于這一類的節(jié)點(diǎn)又按給定有功、無功功率的發(fā)電廠母線和沒有其它電源的變電所母線。待求的則是等值電源的無功功率，也就是要求注入無功功率和節(jié)點(diǎn)電壓的相位角。對這個(gè)節(jié)點(diǎn)，等值負(fù)荷功率、是給定的，節(jié)點(diǎn)電壓的大小和相位角是給定的。（二）高斯——賽德爾潮流計(jì)算自1956年成功地運(yùn)用數(shù)字計(jì)算機(jī)計(jì)算潮流分布以來，曾先后出現(xiàn)過許多種結(jié)算方法。高斯——賽爾德法比較簡單，是由于它可以直接迭代解節(jié)點(diǎn)電壓方程。迭代解式（）的步驟是：先假設(shè)一組，一般可設(shè)，將它們代入第一式，可解得。這就是第一次迭代。然后將、…代入第二式，可解得。第二次迭代結(jié)束時(shí)，解得了所有的。設(shè)節(jié)點(diǎn)是PV節(jié)點(diǎn)，則由于已經(jīng)給定，在每次迭代求得后，應(yīng)首先將求得的修正為，即將求得的電壓大小由改為，而求得的相位則不改動。換言之，這時(shí)只能滿足約束條件而不能滿足=給定值。作如下說明。（2）框1和框2的步驟在調(diào)用高斯——塞爾德子程序之前完成，在調(diào)用高斯塞爾德程序時(shí)輸入的數(shù)據(jù)有：平衡節(jié)點(diǎn)的電壓；PQ節(jié)點(diǎn)的視在功率；PV節(jié)點(diǎn)的無功功率和節(jié)點(diǎn)電壓的大小，；PV節(jié)點(diǎn)的無功功率限額，；（3）對于框4，迭代次數(shù)k的設(shè)置，是為了防止迭代不收斂時(shí)來終止迭代的。 三節(jié)點(diǎn)等值網(wǎng)絡(luò)【解】，我們先建立一個(gè)矩陣，是一個(gè)階的矩陣，其中第3節(jié)點(diǎn)分別對應(yīng)的第3行（列），參考地對應(yīng)第4行（列），的每一個(gè)元素表示連接第個(gè)節(jié)點(diǎn)和第個(gè)節(jié)點(diǎn)的導(dǎo)納，或者表示連接第個(gè)節(jié)點(diǎn)和地之間的導(dǎo)納的值。 第三次迭代后的節(jié)點(diǎn)電壓節(jié)點(diǎn)編號123節(jié)點(diǎn)電壓++最大電壓變化量dumax為：dumax=。計(jì)算結(jié)果如下： 節(jié)點(diǎn)注入功率節(jié)點(diǎn)編號123節(jié)點(diǎn)注入功率++ 各線路上的流動功率 ij12310+2003+0 各線路上的損耗功率 ij12310+2+030三、牛頓——拉夫遜潮流計(jì)算牛頓——拉夫遜是廣泛采用的解非線性方程式組的方法，也是當(dāng)前廣泛采用的計(jì)算潮流分布的方法。以第一式為例式中、…分別表示以、…代入這些偏導(dǎo)數(shù)的表示式時(shí)計(jì)算所得；則是一包含、…的高次方與的高階偏導(dǎo)數(shù)乘積的函數(shù)，如果初值與精確解相差不大，則可以略去。然后運(yùn)用解線性矩陣方程的解法，可求的，從而可以求得經(jīng)第一次迭代后的新值=。正因如此，在運(yùn)用牛頓——拉夫遜法的潮流計(jì)算中，第一、二次的迭代先采用高斯——賽德爾法得到一個(gè)比較精確的初值，然后在用牛頓——拉夫遜法進(jìn)行迭代。有待解決的問題就是各節(jié)點(diǎn)功率不平衡量都趨近于零時(shí)，各節(jié)點(diǎn)電壓應(yīng)具有何值。節(jié)點(diǎn)電壓以直角坐標(biāo)表示時(shí)，令，并將式（）改寫為 ()將實(shí)部和虛部分開寫，可得 () ()鑒于系統(tǒng)中還有電壓大小給定的PV節(jié)點(diǎn)，還應(yīng)由列出 ()式中，第一部分為給定的PV節(jié)點(diǎn)的節(jié)點(diǎn)電壓的平方，第二部分為求得的節(jié)點(diǎn)電壓的平方，它們二者的差可看做為節(jié)點(diǎn)電壓大小的不平衡量。這樣就可建立類似式（）的修正方程式如下 () 式中的、分別就是如式組（）所示。（3）迭代次數(shù)k設(shè)為1000，每次迭代后減1，直至為零，如果迭代不收斂，則當(dāng)k減為零時(shí)，程序自動終止。圖中，節(jié)點(diǎn)1為平衡節(jié)點(diǎn)，給定；其它節(jié)點(diǎn)都是PQ節(jié)點(diǎn)，給定。在進(jìn)行牛頓——拉夫遜潮流計(jì)算之前，我們先調(diào)用高斯——賽德爾潮流計(jì)算程序計(jì)算初值。根據(jù)式（）和（）計(jì)算PV節(jié)點(diǎn)的無功功率和平衡節(jié)點(diǎn)的視在功率，根據(jù)（）計(jì)算線路上的視在功率，根據(jù)式（）計(jì)算線路上的功率損耗。各節(jié)點(diǎn)之間輸電線路的長度，輸電線路單位長度（1千米）的阻抗、導(dǎo)納以及變壓器的變比和阻抗（折算到低壓側(cè)）。一、輸入數(shù)據(jù)和輸出數(shù)據(jù)輸入數(shù)據(jù)和輸出數(shù)據(jù)都以數(shù)組的形式在Excel文檔中讀寫，組織好輸入數(shù)據(jù)和輸出數(shù)據(jù)是編寫程序的關(guān)鍵一步。（一）支路導(dǎo)納矩陣的計(jì)算 ，第一步數(shù)要計(jì)算支路導(dǎo)納矩陣。 主程序流程圖在步驟1中，需要調(diào)用輸電線路和變壓器數(shù)學(xué)模型的計(jì)算子程序，在步驟3中需要調(diào)用高斯——斯德爾潮流計(jì)算子程序。寫文件子程序，主要完成計(jì)算結(jié)果的保存，計(jì)算結(jié)果保存在excel文件中。for i=1:n for j=1:n if(i==j) Y(i,j)=sum(J(i,:))。)以上是所有子程序的設(shè)計(jì)方法，每一個(gè)子程序分別完成了不同的工作，從而將整個(gè)潮流計(jì)算完成。牛頓——拉夫遜法具有很好的收斂性，而且計(jì)算速度快。為了加快計(jì)算速度，可以對節(jié)點(diǎn)導(dǎo)納矩陣進(jìn)行部分修改，而不必進(jìn)行重新計(jì)算。% Qmax；PV節(jié)點(diǎn)無功功率最大值% Qmin：PV節(jié)點(diǎn)無功功率最小值% epxl: 最大電壓變量dumax的最大值，即規(guī)定的潮流計(jì)算的精度，實(shí)數(shù)。%節(jié)點(diǎn)電壓的初始值；U(s)=us。%PV節(jié)點(diǎn)在計(jì)算過程中轉(zhuǎn)化為PQ節(jié)點(diǎn)的標(biāo)志，PV轉(zhuǎn)化為PQ節(jié)點(diǎn)時(shí)，相應(yīng)的標(biāo)志位置1Q_PV=ones(1,nm)。%第a個(gè)節(jié)點(diǎn)是PV節(jié)點(diǎn)，暫時(shí)將PV節(jié)點(diǎn)電壓給定值賦給節(jié)點(diǎn)電壓初始%值,其相角暫時(shí)為0endwhile(k=1)for i=1:n u=U(i)。 for a=1:n if(a~=i) A=A+Y(i,a)*U(a)。 flag(im)==0) %第i各節(jié)點(diǎn)是否為PV節(jié)點(diǎn)？ Q_PV(im)=imag(U(i)39。amp。 else if(Q_PV(im)Qmax) Q_PV(im)=Qmax。,i)。PV節(jié)點(diǎn)%gQ小于Qmin轉(zhuǎn)化為PQ節(jié)點(diǎn)\n39。 end else if(i=m)%第i個(gè)節(jié)點(diǎn)是PQ節(jié)點(diǎn) U(i)=(S_PQ(i)39。amp。 if(dudumax) dumax=du。,1000k+1) break。高斯——塞爾德潮流計(jì)算不能收斂！?。n39。%平衡節(jié)點(diǎn)的視在功率 for i=1:n%計(jì)算各個(gè)節(jié)點(diǎn)的視在功率 if(i==s) s_node(i)=Ss。%直角坐標(biāo)系%輸出值：s_node：節(jié)點(diǎn)視在功率% u_node：節(jié)點(diǎn)電壓%輸入值：Y： 導(dǎo)納矩陣% m： PQ節(jié)點(diǎn)數(shù),包含平衡節(jié)點(diǎn)% n： 總結(jié)點(diǎn)數(shù)% s： 平衡節(jié)點(diǎn)編號，一般定為1號編號% S_PQ: PQ節(jié)點(diǎn)的視在功率[1:m],其中平衡節(jié)點(diǎn)的視在功率暫定為0；% P_PV：PV節(jié)點(diǎn)的有功功率[m+1:n]% U_PV: PV節(jié)點(diǎn)的電壓給定值,只有幅值，沒有相位角，是實(shí)數(shù)[m+1:n]% us: 平衡節(jié)點(diǎn)的電壓給定值，復(fù)數(shù)。 fprintf(39。%雅克比矩陣；I=zeros(1,n)。%設(shè)置迭代次數(shù)Ss=0。i=1。*U(a)39。 %*************************************** if(i=s1)%如果i號節(jié)點(diǎn)是平衡節(jié)點(diǎn)之前的PQ節(jié)點(diǎn) dPQV2(2*i1)=real(S_PQ(i))real(A)。 i=m)%如果i號節(jié)點(diǎn)是平衡節(jié)點(diǎn)之后的PQ節(jié)點(diǎn) dPQV2(2*(i1)1)=real(S_PQ(i))real(A)。 i=n)%如果i號節(jié)點(diǎn)是PV節(jié)點(diǎn) dPQV2(2*(i1)1)=P_PV(im)real(A)。 end a=1。 for i=1:n %計(jì)算雅克比矩陣 if(i=s1)%如果i號節(jié)點(diǎn)是平衡節(jié)點(diǎn)之前的PQ節(jié)點(diǎn) for a=1:n if(a~=i) if(a=s1)%如果a號節(jié)點(diǎn)是平衡節(jié)點(diǎn)之前的PQ節(jié)點(diǎn) Ykb(2*i1,2*a1)=imag(Y(i,a))*real(U(i))real(Y(i,a))*imag(U(i))。%PQ節(jié)點(diǎn)的Lia參數(shù)=Hia else if(a=s+1 amp。%PQ節(jié)點(diǎn)Hia參數(shù) Ykb(2*i1,2*(a1))=real(Y(i,a))*real(U(i))imag(Y(i,a))*imag(U(i))。amp。%PQ節(jié)點(diǎn)Nia參數(shù)Ykb(2*i,2*(a1)1)=Ykb(2*i1,2*(a1))。%PQ節(jié)點(diǎn)Nia參數(shù)Ykb(2*i,2*a1)=Ykb(2*i1,2*a)2*real(I(i))。 a=m)%如果a號節(jié)點(diǎn)是平衡節(jié)點(diǎn)之后的PQ節(jié)點(diǎn) Ykb(2*i1,2*(a1)1)=imag(Y(i,a))*real(U(i))real(Y(i,a))*imag(U(i))imag(I(i))。%PQ節(jié)點(diǎn)的Lia參數(shù)=Hia+2imag(I(i))。%PQ節(jié)點(diǎn)Hia參數(shù) Ykb(2*i1,2*(a1))=real(Y(i,a))*real(U(i))imag(Y(i,a))*imag(U(i))real(I(i))。 end end end end a=1。Ybk中的ii1。%PQ節(jié)點(diǎn)Hia參數(shù) Ykb(2*(i1)1,2*a)=real(Y(i,a))*real(U(i))imag(Y(i,a))*imag(U(i))。amp。 Ykb(2*(i1)1,2*(a1)1)=imag(Y(i,a))*real(U(i))real(Y(i,a))*imag(U(i))。%PQ節(jié)點(diǎn)的Lia參數(shù)=Hia elseif(a=m+1 amp。%PQ節(jié)點(diǎn)Nia參數(shù) Ykb(2*(i1),2*(a1)1)=Ykb(2*(i1)1,2*(a1))。a不變 Ykb(2*(i1)1,2*a1)=imag(Y(i,a))*real(U(i))real(Y(i,a))*imag(U(i))imag(I(i))。%PQ節(jié)點(diǎn)的Lia參數(shù) else if(a=s+1 amp。%PQ節(jié)點(diǎn)Nia參數(shù) Ykb(2*(i1),2*(a1)1)=Ykb(2*(i1)1,2*(a1))2*real(I(i))。 a=n