復(fù)雜網(wǎng)絡(luò)NR法潮流分析與計(jì)算的設(shè)計(jì)

1、電氣工程及其自動(dòng)化專業(yè)課程設(shè)計(jì)復(fù)雜網(wǎng)絡(luò)N-R法潮流分析與計(jì)算的設(shè)計(jì)學(xué)生學(xué)號(hào)： 學(xué)生姓名： 班 級(jí)： 指導(dǎo)教師： 起止日期： 哈爾濱工程大學(xué)自動(dòng)化學(xué)院課程設(shè)計(jì)報(bào)告撰寫內(nèi)容一、設(shè)計(jì)要求（宋體，小四號(hào)字，加黑）用matlab編程，N_R法計(jì)算潮流分布具體要求為：（1）給出程序，并給出注釋（2）輸出迭代次數(shù)，各節(jié)點(diǎn)電壓，各支路電流（3）在圖中標(biāo)明功率流向節(jié)點(diǎn)數(shù)據(jù)如下表所示（標(biāo)幺值）123456P31.80.63.55Q10.50.81.3V11.05支路及變壓器數(shù)據(jù)線路T1T2L2L3L4L5阻抗j0.04j0.020.06+j0.0250.01+j0.20.06+j0.50.05+j0.3導(dǎo)納/2j

2、0.25j0.25j0.25j0.25變比1.05:11.05:1精度要求：0.0001二、設(shè)計(jì)方案（要求給出詳細(xì)的設(shè)計(jì)思路及其必要的論證）（1.）潮流計(jì)算的方法（1）高斯雅克比迭代法（2）高斯-塞得爾法（對(duì)初值要求底，迭代次數(shù)多）（3）牛頓-拉夫遜法（使用廣泛）（4）PQ 快速分解法（提升運(yùn)算速度）目前廣泛應(yīng)用的潮流計(jì)算方法都是基于節(jié)點(diǎn)電壓法的，以節(jié)點(diǎn)導(dǎo)納矩陣Y作為電力網(wǎng)絡(luò)的數(shù)學(xué)模型。節(jié)點(diǎn)電壓Ui 和節(jié)點(diǎn)注入電流Ii 由節(jié)點(diǎn)電壓方程YV=I （1）根據(jù)S=VI(I為I 的共軛)可得非線性的節(jié)點(diǎn)方程YV=I=(S/V) （2）在實(shí)際的電力系統(tǒng)中，已知的運(yùn)行條件不是節(jié)點(diǎn)的注入電流，而是負(fù)荷和發(fā)電

3、機(jī)的功率，而且這些功率一般不隨節(jié)點(diǎn)電壓的變化而變化。由于各節(jié)點(diǎn)注入功率與注入電流的關(guān)系為SiPijQi =ViIi，因此可將式（2）改寫為Ii=Si/Vi=Pi+jQi/Vi (i= 1, 2,3 n) （3）式中，Pi 和Qi 分別為節(jié)點(diǎn)i 向網(wǎng)絡(luò)注入的有功功率和無功功率，當(dāng)i 為發(fā)電機(jī)節(jié)點(diǎn)時(shí)Pi0；當(dāng)i 為負(fù)荷節(jié)點(diǎn)時(shí)Pi0；當(dāng)i 為無源節(jié)點(diǎn)Pi0，Qi0；Vi 和Ii 分別為節(jié)點(diǎn)電壓相量Vi 和節(jié)點(diǎn)注入電流相量Ii 的共軛。式（3）亦即潮流計(jì)算的基本方程式,它可以在直角坐標(biāo)也可以在極坐標(biāo)上建立2n 個(gè)實(shí)數(shù)形式功率方程式。發(fā)電機(jī)Pi、Qi 為正，負(fù)荷Pi、Qi 為負(fù)。展開YV=I 為Ii=Y

4、ijVj=YiiVi+YijVi( i=1 2 3 n) (4)將式（4）代入式（3），得n 維的非線性復(fù)數(shù)的電壓方程組潮流計(jì)算的基本方程為 (Pi-jQi)/ Vi= YiiVi+YijVi （i=1, 2, n） (5)（2.）變量的分類假設(shè)系統(tǒng)中有n 個(gè)節(jié)點(diǎn)，構(gòu)成n 個(gè)復(fù)數(shù)方程，2n 個(gè)實(shí)數(shù)方程，變量總數(shù)為6n 個(gè)。a)不可控變量（2n 個(gè)）：負(fù)荷消耗的有功功率Li P 和無功功率Li Q .由于該類變量無法控制，取決于用戶，而且出現(xiàn)事先沒有預(yù)計(jì)的變動(dòng)，使系統(tǒng)偏離原始運(yùn)行狀態(tài)，因此又稱為不可控變量或擾動(dòng)變量。b)控制變量（2n 個(gè)）:發(fā)電機(jī)發(fā)出的有功功率Gi P 和無功功率Gi Q ，因

5、為該類變量可控。也稱獨(dú)立變量。c)狀態(tài)變量（2n 個(gè)）：母線電壓或節(jié)點(diǎn)電壓的幅值大小i V 與相角大小i ，又稱依從變量或因變量。并且i V 受Gi P 控制， i 受Gi Q 控制。其中2n 個(gè)擾動(dòng)變量是給定的，2n 個(gè)控制變量和2n 個(gè)狀態(tài)變量中給定兩個(gè)，求另外兩個(gè)。（3.）變量的約束條件a)擾動(dòng)變量沒有約束條件。b)控制變量約束條件：為滿足發(fā)電機(jī)的技術(shù)經(jīng)濟(jì)特性指標(biāo)。c)狀態(tài)變量的i V 的約束條件：保證良好的電能質(zhì)量。d) 狀態(tài)變量的i 的約束條件：保證系統(tǒng)的穩(wěn)定運(yùn)行。(4.)系統(tǒng)節(jié)點(diǎn)的分類，根據(jù)給定的控制變量和狀態(tài)變量進(jìn)行分類如下：（1）PQ 節(jié)點(diǎn)（即負(fù)荷節(jié)點(diǎn)）：給定Pi、Qi，求Vi

6、 和i （ i i e , f ）。通常變電所都是這一類型的節(jié)點(diǎn)，由于沒有發(fā)電設(shè)備，因而發(fā)電功率為零電力系統(tǒng)中的絕大多數(shù)節(jié)點(diǎn)屬于這一節(jié)點(diǎn)。其包含變電站節(jié)點(diǎn)（即聯(lián)絡(luò)節(jié)點(diǎn)或浮游節(jié)點(diǎn)）。（2）PV 節(jié)點(diǎn)（即調(diào)節(jié)節(jié)點(diǎn)、電壓控制節(jié)點(diǎn)）：給定Pi 和Vi，求Qi 和i （ i i e , f ）。這類節(jié)點(diǎn)必須有足夠的可調(diào)無功容量，用以維持給定的電壓幅值。一般時(shí)選擇有一頂武功儲(chǔ)備的發(fā)電廠和具有可調(diào)無功電源設(shè)備的變電所作為PV 節(jié)點(diǎn)。在電力系統(tǒng)中，這類節(jié)點(diǎn)數(shù)很少。（3）平衡節(jié)點(diǎn)（即松弛節(jié)點(diǎn)、參考節(jié)點(diǎn)、基準(zhǔn)節(jié)點(diǎn)）：給定Vi 和i （ i =0），求Pi 和Qi。（只有一個(gè)）有功功率不能給定，這個(gè)節(jié)點(diǎn)承擔(dān)了系統(tǒng)的

7、有功功率平衡。同時(shí)其電壓幅值也是給定的，相位為零。(5. )P-Q 分解法是從改進(jìn)和簡(jiǎn)化牛頓法潮流程序的基礎(chǔ)上提出來的，它的基本思想是：把節(jié)點(diǎn)功率表示為電壓向量的極坐標(biāo)方程式，抓住主要矛盾，以有功功率誤差作為修正電壓向量角度的依據(jù)，以無功功率誤差作為修正電壓幅值的依據(jù)，把有功功率和無功功率迭代分開來進(jìn)行。牛頓法潮流程序的核心是求解修正方程式，當(dāng)節(jié)點(diǎn)功率方程式采取極坐標(biāo)系統(tǒng)時(shí)，修正方程式展開為：P = H + NV/ VQ = J + LV /V以上方程式是從數(shù)學(xué)上推倒出來的，并沒有考慮電力系統(tǒng)這個(gè)具體對(duì)象的特點(diǎn)。電力系統(tǒng)中有功功率主要與各節(jié)點(diǎn)電壓向量的角度有關(guān)，無功功率則主要受各節(jié)點(diǎn)電壓幅值的

8、影響。大量運(yùn)算經(jīng)驗(yàn)也告訴我們，矩陣N 及J 中各元素的數(shù)值相對(duì)是很小的，因此對(duì)牛頓法的第一步簡(jiǎn)化就是把有功功率和無功功率分開來進(jìn)行迭代，即將式(4)化簡(jiǎn)為：P = HQ = LV /V (5)這樣，由于我們把2n 階的線性方程組變成了二個(gè)n 階的線性方程組，對(duì)牛頓法的第二個(gè)化簡(jiǎn)，也是比較關(guān)鍵的一個(gè)化簡(jiǎn)，即把式(5)中的系數(shù)矩陣簡(jiǎn)化為在迭代過程中不變的對(duì)稱矩陣。眾所周知，一般線路兩端電壓的相角差是不大的(通常不超過1020 度)，因此可以認(rèn)為： (6)此外，與系統(tǒng)各節(jié)點(diǎn)無功功率相應(yīng)的導(dǎo)納Li B 必定遠(yuǎn)遠(yuǎn)小于該節(jié)點(diǎn)自導(dǎo)納的虛部，即：因此， (7)考慮到以上關(guān)系后，式(5)中系數(shù)矩陣中的元素表達(dá)式

9、可以化簡(jiǎn)為： （8）這樣，式(5)中系數(shù)矩陣可以表示為： （9）進(jìn)一步可以把它們表示為以下矩陣的乘積： （10）將它代入(5)中，并利用乘法結(jié)合率，我們可以把修正方程式變?yōu)椋簩⒁陨蟽墒降淖笥覂蓚?cè)用以下矩陣左乘就可得到以上兩式就是P-Q 分解法達(dá)到修正方程式，其中系數(shù)矩陣只不過是系統(tǒng)導(dǎo)納矩陣的虛部，因而是對(duì)稱矩陣，而且在迭代過程中維持不變。它們與功率誤差方程式構(gòu)成了P-Q 分解法迭代過程中基本計(jì)算公式，其迭代步驟大致是：(1)給定各節(jié)點(diǎn)電壓向量的電壓初值V i (0) ,i (0);(2)根據(jù)(12)計(jì)算各節(jié)點(diǎn)有功功率誤差Pi ，并求出;Pi/Vi(3)解修正方程式(11)，并進(jìn)而計(jì)算各節(jié)點(diǎn)電壓

10、向量角度的修正量i (4)修正各節(jié)點(diǎn)電壓向量角度i ；(5)根據(jù)式(16)計(jì)算各節(jié)點(diǎn)無功功率誤差i Q ，并求出 / ; i i Q V(6)解修正方程式(11)，求出各節(jié)點(diǎn)電壓幅值的修正量i V(7)修正各節(jié)點(diǎn)電壓幅值i V (i ) (i 1) (i 1)i i i V =V V (18)(8)返回(2)進(jìn)行迭代，直到各節(jié)點(diǎn)功率誤差及電壓誤差都滿足收斂條件。P-Q 分解法與牛頓法潮流程序的主要差別表現(xiàn)在它們的修正方程式上。P-Q分解法通過對(duì)電力系統(tǒng)具體特點(diǎn)的分析，對(duì)牛頓法修正方程式的雅可比矩陣進(jìn)行了有效的簡(jiǎn)化和改進(jìn)，有以下三個(gè)特點(diǎn)：（1）在提高計(jì)算速度和減少內(nèi)存方面的作用是明顯的，不再敘述

11、。（2）使我們得到以下好處。首先，因?yàn)樾拚匠淌降南禂?shù)矩陣就是導(dǎo)納矩陣的虛部，因此在迭代過程中不必象牛頓法那樣進(jìn)行形成雅可比矩陣的計(jì)算，這樣不僅是僅減少了運(yùn)算量，而且也大大簡(jiǎn)化了程序。其次，由于系數(shù)矩陣在迭代過程中維持不變，因此在求解修正方程式時(shí)，可以迅速求得修正量，從而顯著提高了迭代速度。（3）可以使我們減少形成因子表時(shí)的運(yùn)算量，而且由于對(duì)稱矩陣三角分解后，其上三角矩陣和下三角矩陣有非常簡(jiǎn)單的關(guān)系，所以在計(jì)算機(jī)中可以只儲(chǔ)存上三角矩陣或下三角矩陣，從而也進(jìn)一步節(jié)約了內(nèi)存。三、設(shè)計(jì)內(nèi)容%本程序的功能是用牛頓拉夫遜法進(jìn)行潮流計(jì)算% B1矩陣：1、支路首端號(hào)； 2、末端號(hào)； 3、支路阻抗； 4、線路

12、對(duì)地電納 (或變壓器導(dǎo)納)；% 5、支路的變比； 6、支路首端處于K側(cè)為1，1側(cè)為0；% 7、線路/變壓器標(biāo)識(shí)（0/1）變壓器參數(shù)當(dāng)支路首端處于K側(cè)標(biāo)識(shí)為1時(shí)歸算至末端側(cè)，0歸算至首端側(cè)% B2矩陣：1、該節(jié)點(diǎn)發(fā)電機(jī)功率； 2、該節(jié)點(diǎn)負(fù)荷功率；% 3、PQ節(jié)點(diǎn)電壓初始值，或平衡節(jié)點(diǎn)及PV節(jié)點(diǎn)電壓的給定值% 4、節(jié)點(diǎn)所接無功補(bǔ)償并聯(lián)電容（感）的電納% 5、節(jié)點(diǎn)分類標(biāo)號(hào)：1為平衡節(jié)點(diǎn)（應(yīng)為1號(hào)節(jié)點(diǎn)）；2為PQ節(jié)點(diǎn)；3為PV節(jié)點(diǎn)；clear;isb=1; %input('請(qǐng)輸入平衡母線節(jié)點(diǎn)號(hào)：isb=');pr=0.0001; %input('請(qǐng)輸入誤差精度：pr='

13、);%-n=6;%input('請(qǐng)輸入節(jié)點(diǎn)數(shù)：n=');nl=6;%input('請(qǐng)輸入支路數(shù)：nl=');B1=1 2 0+0.04i 0 1.05 1 1; 2 3 0.06+0.025i 0+0.5i 1 0 0; 2 5 0.01+0.2i 0+0.5i 1 0 0; 3 4 0.06+0.50i 0+0.5i 1 0 0; 4 5 0.05+0.3i 0 1 0 0; 6 5 0+0.02i 0 1.05 1 1B2=0 0 1 0 1; 0 3+1i 1.00 0 2; 0 1.8+0.50i 1.00 0 2; 0 0.6+0.8i 1.00 0

14、2; 0 3.5+1.3i 1.00 0 2; 0 -5+0i 1.05 0 3 %input('請(qǐng)輸入各節(jié)點(diǎn)參數(shù)形成的矩陣： B2=');%X=1 0;2 0;3 0;4 0;5 0;6 0%-%n=4;%input('請(qǐng)輸入節(jié)點(diǎn)數(shù)：n=');nl=4;%input('請(qǐng)輸入支路數(shù)：nl=');%B1=1 2 4+16i 0 1 0 0;1 3 4+16i 0 1 0 0;2 3 2+8i 0 1 0 0;2 4 1.49+48.02i 0 11/110 0 1 %input('請(qǐng)輸入由支路參數(shù)形成的矩陣： B1=');%B2=

15、0 0 115 0 1;0 0 110 0 2;0 20+4i 110 0 2;0 10+6i 10 0 2 %input('請(qǐng)輸入各節(jié)點(diǎn)參數(shù)形成的矩陣： B2=');%-Y=zeros(n);e=zeros(1,n);f=zeros(1,n);V=zeros(1,n);sida=zeros(1,n);S1=zeros(nl);% % %-求導(dǎo)納矩陣-%for i=1:n % if X(i,2)=0; % p=X(i,1); % Y(p,p)=1/X(i,2); %end%endfor i=1:nl%從1到n1（總支路數(shù)）if B1(i,7)=1 %-如果是變壓器支路- if

16、B1(i,6)=0%左節(jié)點(diǎn)(首端)處于1側(cè) p=B1(i,1);q=B1(i,2); else %左節(jié)點(diǎn)(首端)處于K側(cè) p=B1(i,2);q=B1(i,1); end Y(p,q)=Y(p,q)-1./(B1(i,3)*B1(i,5);%非對(duì)角元 Y(q,p)=Y(p,q); %非對(duì)角元 Y(q,q)=Y(q,q)+1./(B1(i,3)*B1(i,5)2);%對(duì)角元K側(cè) Y(p,p)=Y(p,p)+1./B1(i,3)+B1(i,4);%對(duì)角元1側(cè)+勵(lì)磁導(dǎo)納 else %-否則為線路支路-p=B1(i,1);q=B1(i,2);Y(p,q)=Y(p,q)-1./B1(i,3);%非對(duì)角元

17、 Y(q,p)=Y(p,q); %非對(duì)角元 Y(q,q)=Y(q,q)+1./B1(i,3)+B1(i,4)./2.0000;%對(duì)角元j側(cè)+線路電納的一半 Y(p,p)=Y(p,p)+1./B1(i,3)+B1(i,4)./2.0000;%對(duì)角元i側(cè)+線路電納的一半endenddisp('導(dǎo)納矩陣 Y=');disp(Y);%-給定各節(jié)點(diǎn)初始電壓及給定各節(jié)點(diǎn)注入功率-G=real(Y);B=imag(Y);%分解出導(dǎo)納陣的實(shí)部和虛部 for i=1:n%給定各節(jié)點(diǎn)初始電壓的實(shí)部和虛部 e(i)=real(B2(i,3); f(i)=imag(B2(i,3); V(i)=abs(

18、B2(i,3);%PV、平衡節(jié)點(diǎn)及PQ節(jié)點(diǎn)電壓模值 endfor i=1:n%給定各節(jié)點(diǎn)注入功率 S(i)=B2(i,1)-B2(i,2); %i節(jié)點(diǎn)注入功率SG-SL B(i,i)=B(i,i)+B2(i,4);%i節(jié)點(diǎn)無功補(bǔ)償量(電納值)end%=用牛頓-拉夫遜法迭代求解非線性代數(shù)方程（功率方程）=P=real(S);Q=imag(S); %分解出各節(jié)點(diǎn)注入的有功和無功功率ICT1=0;IT2=1;N0=2*n;N1=N0+1;a=0; %迭代次數(shù)ICT1、a；不滿足收斂要求的節(jié)點(diǎn)數(shù)IT2while IT2=0 % N0=2*n 雅可比矩陣的階數(shù)；N1=N0+1擴(kuò)展列 IT2=0;a=a+

19、1; JZ='Jacobi矩陣第(',num2str(a),')次消去運(yùn)算'JZ1='Jacobi矩陣第(',num2str(a),')次回代運(yùn)算'JZ0='功率方程第(',num2str(a),')次差值：' %-求取各個(gè)節(jié)點(diǎn)的功率及功率偏差及PV節(jié)點(diǎn)的電壓偏差- for i=1:n %n個(gè)節(jié)點(diǎn)2n行(每節(jié)點(diǎn)兩個(gè)方程P和Q或U) p=2*i-1;m=p+1;C(i)=0;D(i)=0; for j1=1:n %第i行共n列(n個(gè)節(jié)點(diǎn)間互導(dǎo)納及節(jié)點(diǎn)電壓相乘即電流) C(i)=C(i)+G(i,j1)

20、*e(j1)-B(i,j1)*f(j1);%(Gij*ej-Bij*fj) D(i)=D(i)+G(i,j1)*f(j1)+B(i,j1)*e(j1);%(Gij*fj+Bij*ej) end %求i節(jié)點(diǎn)有功和無功功率P',Q'的計(jì)算值 P1=C(i)*e(i)+f(i)*D(i);%節(jié)點(diǎn)功率P計(jì)算ei(Gij*ej-Bij*fj)+fi(Gij*fj+Bij*ej) Q1=C(i)*f(i)-e(i)*D(i);%節(jié)點(diǎn)功率Q計(jì)算fi(Gij*ej-Bij*fj)-ei(Gij*fj+Bij*ej) V2=e(i)2+f(i)2;%電壓模平方 %=求取功率差及PV節(jié)點(diǎn)電壓模平方

21、差 = if i=isb%非平衡節(jié)點(diǎn)(PQ或PV節(jié)點(diǎn)) if B2(i,5)=3%非PV節(jié)點(diǎn)(只能是PQ節(jié)點(diǎn)) J(m,N1)=P(i)-P1;%PQ節(jié)點(diǎn)有功功率差J(m,N1)擴(kuò)展列P J(p,N1)=Q(i)-Q1; %PQ節(jié)點(diǎn)無功功率差J(p,N1)擴(kuò)展列Q else %PV節(jié)點(diǎn)= J(m,N1)=P(i)-P1;%PV節(jié)點(diǎn)有功功率差J(m,N1)擴(kuò)展列P J(p,N1)=V(i)2-V2;%PV節(jié)點(diǎn)電壓模平方差J(p,N1)擴(kuò)展列U end end %(if i=isb)非平衡節(jié)點(diǎn)(PQ或PV節(jié)點(diǎn)) end %(for i=1:n) n個(gè)節(jié)點(diǎn)2n行(每節(jié)點(diǎn)兩個(gè)方程P和Q或U) for

22、m=1:N0 JJN1(m)=J(m,N1); end disp(JZ0);disp(JJN1); %-判斷功率偏差量及PV節(jié)點(diǎn)的電壓偏差量是否滿足要求- for k=3:N0 %除去平衡節(jié)點(diǎn)1、2號(hào)以外的所有節(jié)點(diǎn) DET=abs(J(k,N1); if DET>=pr; %PQ節(jié)點(diǎn)的功率偏差量及PV節(jié)點(diǎn)的電壓偏差量是否滿足要求 IT2=IT2+1; %不滿足要求的節(jié)點(diǎn)數(shù)加1 end end ICT2(a)=IT2; %不滿足要求的節(jié)點(diǎn)數(shù)；a為迭代次數(shù) ICT1=ICT1+1; %迭代次數(shù) if ICT2(a)=0; %當(dāng)前不滿足要求的節(jié)點(diǎn)數(shù)為零 break %退出迭代運(yùn)算 end %-

23、以上為求取各個(gè)節(jié)點(diǎn)的功率及功率偏差及PV節(jié)點(diǎn)的電壓偏差- %= 求取Jacobi矩陣形成修正方程 = for i=2:n %n個(gè)節(jié)點(diǎn)2n行(每節(jié)點(diǎn)兩個(gè)方程P和Q或U) if i=isb%非平衡節(jié)點(diǎn)(PQ或PV節(jié)點(diǎn)) if B2(i,5)=3 %下面是針對(duì)PQ節(jié)點(diǎn)來求取Jacobi矩陣的元素 = C(i)=0;D(i)=0; for j1=1:n %第i行共n列(n個(gè)節(jié)點(diǎn)間互導(dǎo)納及節(jié)點(diǎn)電壓相乘即電流) C(i)=C(i)+G(i,j1)*e(j1)-B(i,j1)*f(j1);%(Gij*ej-Bij*fj) D(i)=D(i)+G(i,j1)*f(j1)+B(i,j1)*e(j1);%(Gij

24、*fj+Bij*ej) end for j1=2:n %第i行共n列(2n個(gè)Jacobi矩陣元素dP/de及dP/df或dQ/de及dQ/df) if j1=isb&j1=i%非平衡節(jié)點(diǎn)&非對(duì)角元 X1=-G(i,j1)*e(i)-B(i,j1)*f(i);% X1=dP/de=-dQ/df=-X4 X2=B(i,j1)*e(i)-G(i,j1)*f(i);% X2=dP/df=dQ/de=X3 X3=X2; % X2=dp/df X3=dQ/de X4=-X1; % X1=dP/de X4=dQ/df p=2*i-1;q=2*j1-1; J(p,q)=X3;m=p+1; %

25、X3=dQ/de J(p,N)=DQ節(jié)點(diǎn)無功功率差 J(p,N)=DQ; J(m,q)=X1;q=q+1; % X1=dP/de J(m,N)=DP節(jié)點(diǎn)有功功率差 J(m,N)=DP; J(p,q)=X4;J(m,q)=X2; % X4=dQ/df X2=dp/df elseif j1=i&j1=isb%非平衡節(jié)點(diǎn)&對(duì)角元 X1=-C(i)-G(i,i)*e(i)-B(i,i)*f(i);% dP/de X2=-D(i)+B(i,i)*e(i)-G(i,i)*f(i);% dP/df X3=D(i)+B(i,i)*e(i)-G(i,i)*f(i); % dQ/de X4=-C(

26、i)+G(i,i)*e(i)+B(i,i)*f(i);% dQ/df p=2*i-1;q=2*j1-1;J(p,q)=X3;%擴(kuò)展列Q J(p,N)=DQ; m=p+1; J(m,q)=X1;q=q+1;J(p,q)=X4;%擴(kuò)展列P J(m,N)=DP; J(m,q)=X2; end end else %if B2(i,5)=3 % 否則（即為PV節(jié)點(diǎn)） %= 下面是針對(duì)PV節(jié)點(diǎn)來求取Jacobi矩陣的元素 = for j1=1:n if j1=isb&j1=i%非平衡節(jié)點(diǎn)&非對(duì)角元 X1=-G(i,j1)*e(i)-B(i,j1)*f(i); % dP/de X2=B(i,

27、j1)*e(i)-G(i,j1)*f(i); % dP/df X5=0;X6=0; p=2*i-1;q=2*j1-1;J(p,q)=X5; % PV節(jié)點(diǎn)電壓誤差J(p,N)=DV; m=p+1; J(m,q)=X1;q=q+1;J(p,q)=X6; % PV節(jié)點(diǎn)有功誤差J(m,N)=DP; J(m,q)=X2; elseif j1=i&j1=isb %非平衡節(jié)點(diǎn)&對(duì)角元 X1=-C(i)-G(i,i)*e(i)-B(i,i)*f(i);% dP/de X2=-D(i)+B(i,i)*e(i)-G(i,i)*f(i);% dP/df X5=-2*e(i); X6=-2*f(i);

28、 p=2*i-1;q=2*j1-1;J(p,q)=X5; % PV節(jié)點(diǎn)電壓誤差J(p,N)=DV; m=p+1; J(m,q)=X1;q=q+1;J(p,q)=X6; % PV節(jié)點(diǎn)有功誤差J(m,N)=DP; J(m,q)=X2; end end end %(if B2(i,5)=3 else) end %(if i=isb) end %(for i=1:n)n個(gè)節(jié)點(diǎn)2n行(每節(jié)點(diǎn)兩個(gè)方程P和Q或U) JZ0='形成的第(',num2str(a),')次Jacobi矩陣：' disp(JZ0);disp(J); %= 以上為形成完整的Jacobi矩陣 = %=下

29、面用高斯消去法對(duì)由Jacobi矩陣形成的修正方程進(jìn)行求解(按列消去、回代) = for k=3:N0 % N0=2*n （從第三行開始，第一、二行是平衡節(jié)點(diǎn)） for k1=k+1:N1% 從k+1列的Jacobi元素到擴(kuò)展列的P、Q 或 U J(k,k1)=J(k,k1)./J(k,k);% 用K行K列對(duì)角元素去除K行K列后的非對(duì)角元素進(jìn)行規(guī)格化 end J(k,k)=1; % 對(duì)角元規(guī)格化K行K列對(duì)角元素賦1 %= 按列消去運(yùn)算 = for k2=k+1:N0 % 從k+1行到2*n最后一行 for k3=k+1:N1 %從k2+1列到擴(kuò)展列消去k+1行后各行下三角元素 J(k2,k3)=

30、J(k2,k3)-J(k2,k)*J(k,k3);%消去運(yùn)算 end %用當(dāng)前行K3列元素減去當(dāng)前行K列元素乘以第k行K3列元素 J(k2,k)=0; %當(dāng)前行第k列元素已消為0 end end JZ='Jacobi矩陣第(',num2str(a),')次消去運(yùn)算'JZ1='Jacobi矩陣第(',num2str(a),')次回代運(yùn)算' disp(JZ);disp(J); %= 按列回代運(yùn)算 = for k=N0:-1:3 for k1=k-1:-1:3 J(k1,N1)=J(k1,N1)-J(k1,k)*J(k,N1); J(k

31、1,k)=0; end end for m=1:N0 JJN1(m)=J(m,N1); end disp(JZ1);disp(JJN1);%disp(J); %-修改節(jié)點(diǎn)電壓- for k=3:2:N0-1 L=(k+1)./2; e(L)=e(L)-J(k,N1); %修改節(jié)點(diǎn)電壓實(shí)部 k1=k+1; f(L)=f(L)-J(k1,N1); %修改節(jié)點(diǎn)電壓虛部 U(L)=sqrt(e(L)2+f(L)2); end disp('各個(gè)節(jié)點(diǎn)電壓模');disp(U); %= 結(jié)束一次迭代 =end%* 下面為迭代計(jì)算結(jié)束后的有關(guān)輸出過程 *disp('迭代次數(shù)：'

32、);disp(ICT1-1);disp('沒有達(dá)到精度要求的個(gè)數(shù)：');disp(ICT2);for k=1:n V(k)=sqrt(e(k)2+f(k)2); %計(jì)算各節(jié)點(diǎn)電壓的模值 sida(k)=atan(f(k)./e(k)*180./pi; %計(jì)算各節(jié)點(diǎn)電壓的角度 E(k)=e(k)+f(k)*j; %將各節(jié)點(diǎn)電壓用復(fù)數(shù)表示end%= 計(jì)算各輸出量 =disp('各節(jié)點(diǎn)的電壓復(fù)數(shù)值E為(節(jié)點(diǎn)號(hào)從小到大排列)：');disp(E); %顯示各節(jié)點(diǎn)的實(shí)際電壓值E用復(fù)數(shù)表示disp('-');disp('各節(jié)點(diǎn)的電壓模值大小V為(節(jié)點(diǎn)

33、號(hào)從小到大排列)：');disp(V); %顯示各節(jié)點(diǎn)的電壓大小V的模值disp('-');disp('各節(jié)點(diǎn)的電壓相角sida為(節(jié)點(diǎn)號(hào)從小到大排列)：');disp(sida); %顯示各節(jié)點(diǎn)的電壓相角for p=1:n C(p)=0; for q=1:n C(p)=C(p)+conj(Y(p,q)*conj(E(q); %計(jì)算各節(jié)點(diǎn)注入電流的共軛值 end S(p)=E(p)*C(p); %計(jì)算各節(jié)點(diǎn)的功率 S = 電壓 X 注入電流的共軛值enddisp('各節(jié)點(diǎn)的功率S為(節(jié)點(diǎn)號(hào)從小到大排列)：');disp(S); %顯示各節(jié)

34、點(diǎn)的注入功率disp('-');disp('各條支路的首端功率Si為(順序同您輸入B1時(shí)一致)：');for i=1:nl p=B1(i,1);q=B1(i,2); if B1(i,7)=0 Si(p,q)=E(p)*conj(E(p)*B1(i,4)./2+(E(p)-E(q)./B1(i,3); Siz(i)=Si(p,q); else if B1(i,6)=0 Si(p,q)=E(p)*(conj(E(p)*B1(i,4). +(E(p)*B1(i,5)-E(q)*(1./(B1(i,3)*B1(i,5); Siz(i)=Si(p,q); else Si(

35、p,q)=E(p)*conj(E(p)-E(q)*B1(i,5)*(1./(B1(i,3)*B1(i,5)2); Siz(i)=Si(p,q); end end ZF='S(',num2str(p),',',num2str(q),')=',num2str(Si(p,q); disp(ZF); disp('-');enddisp('各條支路的末端功率Sj為(順序同您輸入B1時(shí)一致)：');for i=1:nl p=B1(i,1);q=B1(i,2); if B1(i,7)=0 Sj(q,p)=E(q)*conj(E(

36、q)*B1(i,4)./2+(E(q)-E(p)./B1(i,3); Sjy(i)=Sj(q,p); else if B1(i,6)=0 Sj(q,p)=E(q)*conj(E(q)-E(p)*B1(i,5)*(1./(B1(i,3)*B1(i,5)2); Sjy(i)=Sj(q,p); else Sj(q,p)=E(q)*(conj(E(q)*B1(i,4). +(E(q)*B1(i,5)-E(p)*(1./(B1(i,3)*B1(i,5); Sjy(i)=Sj(q,p); end end ZF='S(',num2str(q),',',num2str(p),&

37、#39;)=',num2str(Sj(q,p); disp(ZF); disp('-');enddisp('各條支路的功率損耗DS為(順序同您輸入B1時(shí)一致)：');for i=1:nl p=B1(i,1);q=B1(i,2); DS(i)=Si(p,q)+Sj(q,p); ZF='DS(',num2str(p),',',num2str(q),')=',num2str(DS(i); disp(ZF); disp('-'); end zws=0;JDDY=;JDP=;JDQ=;JDDYJD=;

38、 for i=1:n %總網(wǎng)損為所有節(jié)點(diǎn)注入功率的代數(shù)和 zws=zws+S(i); JDDYJD=strcat(JDDYJD,num2str(i),'(',num2str(sida(i),'),'); JDDY=strcat(JDDY,num2str(i),'(',num2str(V(i),'),'); JDP=strcat(JDP,num2str(i),'(',num2str(real(S(i),'),'); JDQ=strcat(JDQ,num2str(i),'(',num2s

39、tr(imag(S(i),'),');endJDDYJD=strcat('節(jié)點(diǎn)電壓角度：',JDDYJD);JDDY=strcat('節(jié)點(diǎn)電壓模值：',JDDY);JDP=strcat('節(jié)點(diǎn)有功：',JDP);JDQ=strcat('節(jié)點(diǎn)無功：',JDQ); ZF subplot(4,1,2);plot(sida);title(JDDYJD);ylabel('節(jié)點(diǎn)電壓角度');grid on;subplot(4,1,3);P=real(S);Q=imag(S);bar(P);title(JDP);

40、ylabel('節(jié)點(diǎn)注入有功');grid on;subplot(4,1,4);bar(Q);title(JDQ);xlabel(ZF);ylabel('節(jié)點(diǎn)注入無功');grid on;%*figure(2);subplot(3,2,1);JDH=;JDH1=;for i=1:nl JDH=strcat(JDH,num2str(i),'(',num2str(B1(i,1),',',num2str(B1(i,2),'), '); JDH1=strcat(JDH1,num2str(i),'(',num2str(B1(i,2),',',num2str(B1(i,1),'), ');endP1=real(Siz);Q1=imag(Siz);bar(P1);title(JDH);ylabel('支路首端注入有功');%xlabel('支路號(hào)');grid on;subplot(3,2,2);bar(Q1);title(JDH);ylabel('支路首端注入無功');grid on;subplot(3,