第一篇：第5章高等數學計算的MATLAB實現講稿

第5章 高等數學計算的MATLAB實現

高等數學是大學數學學習的基本內容。利用MATLAB的符號工具箱，可以解決極限、導數、微分、積分、級數和微分方程等方面的問題。

5.1 函數和極限

5.1.1 函數

使用符號表達式，可以進行復合函數運算和反函數運算，下面分別予以介紹。

⒈ 復合函數運算

在MATLAB中，符號表達式的復合函數運算主要是通過函數compose來實現的。compose函數的調用格式如下：

●compose(f,g)：返回復合函數f(g(y))。在這里f?f(x),g?g(y)。其中，x是findsym定義的f函數的符號變量，y是findsym定義的g函數的符號變量。

●compose(f,g,z)：返回自變量為z的復合函數f(g(z))。在這里，f=f(x),g=g(y)，x，y分別是findsym定義的f函數和g函數的符號變量。

●compose(f,g,x,z)：返回復合函數f(g(z))，并且使x成為f函數的獨立變量。如：f=cos(x/t),compose(f,g,x,z)返回cos(g(z)/t)，而compose(f,g,t,z)返回cos(x/g(z))。

●compose(f,g,x,y,z)：返回復合函數f(g(z))，并且使x與y分別成為函數f與g的獨立變量。

例5－1 將f?x和x?tany復合到一個函數中，指定x和y為它們的獨立變量，自變量為z。

程序：

syms x y t z;g=tan(y);f=x^t;compose(f,g,x,y,z)%求復合函數 運行結果： ans = tan(z)^t ⒉ 符號表達式的反函數運算 在MATLAB中，符號表達式的反函數運算主要是通過函數finverse來實現的。finverse函數的調用格式如下。

5.2 導數

5.2.1 求函數的導數

在MATLAB中，微分和求導都可以由函數diff實現。diff函數可同時處理數值和符號兩種情況下的求導和微分。該函數的調用格式如下所示。

●diff(F)：對findsym函數返回獨立變量求微分，F為符號表達式。●diff(F,'a')：對a變量求微分，F為符號表達式。

●diff(F,n)：對findsym函數返回的獨立變量求n次微分，F為符號表達式。●diff(F,'a',n)或diff(F,n,'a')：對變量a求n次微分，F為符號表達式。

(x?1)5例5－4 求f(x)?的二階導數。

x?1程序： syms x f=(x-1)^5/(x+1);df=diff(f,1);%求導數 d2f=diff(f,2);df=simplify(df)%化簡 d2f=simplify(d2f)運行結果： df = 2*(x-1)^4*(2*x+3)/(x+1)^2 d2f = 4*(x-1)^3*(3*x^2+9*x+8)/(x+1)^3 5.2.2 求隱函數的導數

例5－5 求隱函數

1F(x,y)?x?y?siny

2所確定的導數dy。dx程序：

%求隱函數的導數

f=sym('x-y+1/2*sin(y)');fx=diff(f,'x');fy=diff(f,'y');

5.3.1 漸近線

求函數圖形的水平漸近線，需要求x趨于無窮時f的極限，即

limit(f,inf)limit(f,-inf)求f的垂直漸近線，使分母等于0，用下面的命令進行求解。roots=solve(denom)%返回方程x?x?3?0的解。綜合程序： syms x num=3*x^2+6*x-1;denom=x^2+x-3;f=num/denom;a=limit(f,inf);b=double(a);roots=solve(denom);ezplot(f)%符號函數作圖命令。

hold on%在原有的圖形上面疊加圖形。plot([-2*pi 2*pi],[b b], 'g')%繪水平漸近線

plot(double(roots(1))*[1 1],[-5 10], 'r')%繪垂直漸近線 plot(double(roots(2))*[1 1],[-5 10], 'r')%繪垂直漸近線 title('水平漸近線和垂直漸近線')hold off%取消圖形疊加 運行結果見圖5－1。5.3.2 極值

23x2?6x?1從圖5－1可以看出，函數f(x)?至少有2個極值點，求解

x2?x?3程序：

syms x num=3*x^2+6*x-1;denom=x^2+x-3;f=num/denom;a=limit(f,inf);b=double(a);roots=solve(denom);ezplot(f)%符號函數作圖命令。

hold on%在原有的圖形上面疊加圖形。plot([-2*pi 2*pi],[b b], 'g')%繪水平漸近線

圖5－2 表示函數的漸近線和極值

5.3.3 拐點

求函數的拐點，需要先求函數的2階導數，后面的處理方法與求極值方法相似。

5.4 不定積分和定積分

MATLAB中，用符號工具箱的int函數求函數的不定積分和定積分。int函數的調用格式如下所示。

●int(F)：對findsym函數返回獨立變量求不定積分，F為符號表達式。

●int(F,v)：對v變量求不定積分，F為符號表達式。

●int(F,a,b)：對findsym函數返回獨立變量求從a到b的定積分，F為符號表達式。

●int(F,v,a,b)：對v變量求從a到b的定積分，F為符號表達式。

5.4.4 定積分的應用

例5－11 計算由兩條拋物線y2?x，y?x2所圍成的圖形的面積。程序：

%求曲線的交點

[x1,y1]=solve('y^2=x','y=x^2');x1=double(x1);y1=double(y1);n=numel(x1);%下面尋找實數解 m=1;x0=[];y0=[];for k=1:n

if isreal(x1(k))&&isreal(x1(k))

x0(m)=x1(k);y0(m)=y1(k);

m=m+1;

end end x0=sort(x0);y0=sort(y0);%排序 %下面計算定積分 syms x f=sqrt(x)-x^2;A=int(f,x,x0(1),x0(2))運行結果： A = 1/3 x2y2例5－12 計算由橢圓2?2?1所圍成的圖形繞x軸旋轉而成的旋轉ab體的體積。

V??程序：

syms a b x f=pi*b*b*(a*a-x*x)/a/a;V=int(f,x,-a,a)運行結果： V = 4/3*pi*b^2*a 例5－13 計算由曲線y?a?a?b2a2(a2?x2)dx

23/2x上相應于x從a到b的一段弧的長度。3-9

0 d =

-3 f =

-1 g =

-3 例5－16 已知三點M(1,1,1)、A(2,2,1)和B(2,1,2)，求?AMB。程序：

M=[1 1 1];A=[2 2 1];B=[2 1 2];ma=A-M;mb=B-M;c=dot(ma,mb)/sqrt(dot(ma,ma))/sqrt(dot(mb,mb));amb=acos(c)運行結果： amb =

1.0472 例5－17 已知三角形ABC的頂點是A(1,2,3)、B(3,4,5)和C(2,4,7)求三角形ABC的面積。

程序：

A=[1 2 3];B=[3 4 5];C=[2 4 7];ab=B-A;ac=C-A;S=sqrt(dot(cross(ab,ac),cross(ab,ac)))/2 運行結果： S =

3.7417 5.5.2 曲面及其方程

利用MATLAB提供的繪圖函數，可以繪制給定函數的曲面。相關內容在前面已介紹過，請參見4.2.4小節。

5.6 多元函數的極限和求導

對于函數有多個變量的情況，求極限和導數時需要指定函數對哪個變量進行求取。在MATLAB中仍然使用limit和diff函數求多元函數的極限和導數。5.6.1 求多元函數的極限

例5－18 求極限limsin(x?y)?sin(x)。

y?0y-11●symsum(s)：求符號表達式s相對于符號變量k的和，k由findsym函數確定，取值從0到k-1。

●symsum(s,v)：求符號表達式s相對于符號變量v的和，v從0到v-1。●symsum(s,a,b)和symsum(s,v,a,b)：指定符號表達式s從v=a累加到v=b。

?1例5－21 求級數?k、?2和?xk(x?1)。

i?0k?1kk?0n?1?程序：

syms x k n s1=symsum(n)s2=symsum(1/k^2,1,inf)s3=symsum(x^k,k,0,inf)運行結果： s1 = 1/2*n^2-1/2*n s2 = 1/6*pi^2 s3 =-1/(x-1)5.7.2 泰勒級數展開

用taylor函數進行泰勒級數展開。該函數的調用格式如下：

●taylor(f,n,v)：返回f的n-1階馬克勞林多項式近似。f為表示函數的符號表達式，v指定表達式中的獨立變量。v可以是字符串或符號變量。

●taylor(f,n,v,a)：返回f關于a的n-1階泰勒級數近似。變量a可以是數值、符號或表示數值值或未知值的字符串。n，v和a的順序沒有先后之分。taylor函數根據變量的位置和類型確定它們的用途。還可以忽略n,v,a等變量中的任何一個。如果不確定v，taylor函數用findsym函數確定函數的獨立變量。n的默認值為6。

f(n)(a)泰勒級數：f(x)??(x?a)n。

n!n?01例5－22 求函數f(x)?的泰勒級數展開，取前9項。

5?4cosx?程序： syms x f=1/(5+4*cos(x));t=taylor(f,9)運行結果：

3?d2y?2?cos2x?y?dx例5－24 求解微分方程?yx?0?1

?dy?x?0?0dx?程序：

y=dsolve('D2y=cos(2*x)-y','y(0)=1','Dy(0)=0','x');%求解微分方程 y=simplify(y)%化簡y的形式。運行結果： y = 4/3*cos(x)-2/3*cos(x)^2+1/3

習題五

1、完成實驗指導書中的實驗五的上半部分。

第二篇：計算講稿

如何做好有效的計算教學

淺談如何提高小學生數學計算能力

福州路小學 張成勛

《新課程標準》關于小學數學計算指出“使學生能夠正確地進行整數、小數、分數的四則計算,對于其中一些基本的計算,要達到一定的熟練程度,并逐步做到計算方法合理、靈活”。

小學計算教學貫穿于小學數學教學的全過程, 計算教學直接關系著學生對數學基礎知識與基本技能的掌握，關系著學生觀察，記憶，注意，思維等能力的發展，關系著學生的學習習慣，情感，意志等非智力因素的培養。“計算”在教學中所占的比重相當大，無論是應用題、統計知識，還是幾何題、簡易方程，都離不開計算。計算的準確率和速度如何，將直接影響學生學習的質量，因此，小學階段的計算教學就顯得異常重要。很多數學教師都有同感:每次批閱學生作業或試卷,總會看到不少學生在計算中出現錯寫,多寫,少寫數字或符號,打錯小數點等現象。下面就提高計算能力的實效性，我談一點個人的體會。

一、激發學生的計算興趣,提高學生的計算能力

1.適時列舉中外數學家的典型事例。教學中教師可以適時列舉中外數學家的典型事例,來激發學生對數學的興趣和愛好,提高計算效果。比如:我國著名數學家陳景潤為了攻克歌德巴赫猜想,草稿紙演算了幾麻袋；“會移動的黑板”中上街散步還在馬車的后背箱上專心致志演算的安培,這些生動典型的事例都可以激發學生對計算的興趣。

2.有效的創設數學情境。計算比較枯燥,因此,教學時要根據小學生的心理特點,將童話、游戲、比賽等融入課堂教學中,同時要注意題目的靈活性、注意練習形式的多樣性,從而激發小學生的計算興趣,提高計算能力。

3.成立數學興趣小組。成立數學興趣小組不僅可以豐富小學生的課外生活,而且可以調動學生學習的積極性、主動性和創造性。可以定期或不定期舉辦數學講座、速算、巧算比賽,從而使學生達到算得準、算得巧的目的,增強計算情趣。

二、使學生理解和牢固掌握有關基礎知識

學生的計算離不開數學概念、運算定律、運算性質、運算法則和計算公式等內容。對學生不易理解的某些計算法則，往往成為教學的難點。在教學中教師不能急于求成，應幫助學生以掌握基礎知識為突破口，分散、突破難點。例如教學同分母分數加減法時，首先要讓學生領會分母相同即分數單位不同，而分數單位相同，就能直接相加減，懂得了這個道理，例如教學異分母分數加減法時，首先要讓學生領會分母不同即分數單位不同，而分數單位不同，就不能直接相加減，懂得了這個道理，再引導學生運用通分的知識，化異分母分數為同分母分數，于是問題就轉化為已學過的同分母分數相加減了。

乘法分配律的理解和掌握對我們計算四則混合運算很重要。怎樣 讓學生從實質上理解乘法分配律？

在乘法分配律的教學中，如果只求形式不求實質理解，一方面從認識的角度看是不嚴謹的，另一方面很容易造成學生不求甚解、囫圇吞棗的不良認知習慣。如果滿足于從形式上掌握乘法分配律，對于學生的后續發展也極為不利。因此，在教學時結合生活中的購物問題。先出示了這樣一道例題：一張桌子65元，一張椅子35元，學校買了15套課桌椅，一共要花多少元？學生用了兩種解答方法65×15 +35×15，（65+35）×15，借助對同一實際問題的不同解決方法讓學生體會乘法分配律的合理性，體會分配律中“分配”的含義

三、理解算理,優化算法，提高計算技巧

新課標要求在教學中要關注學生的學習過程，那么我們在計算教學中也應該關注過程教學。在計算過程中，理解算理很重要，他是計算的前提，而算法優化則是計算的關鍵。學生計算錯誤的原因常常是算理在學習的過程中沒有理解到位。在計算教學中根據知識體系之間的聯系可以在遷移中幫助學生理解算理。例如教學除數是小數的除法，學生已經學習了除數是整數的除法，積累了以下的兩點認識：計算時就按整數除法的方法算出結果；商的小數點和被除數的小數點對齊。這些認識是學生學習除數是小數除法的基礎，在實際教學中教師可以先復習除數是整數的除法，如：“38.4÷24”，在學生明確商的小數點是如何確定后，把復習題改成“3.84÷2.4”，在學生嘗試計算中著重引導學生分析怎樣把除數是小數轉化成除數是整數的除法，在學生初步理解算理的基礎上進行“試一試”的教學：

0.12÷3= 0.12÷0.03=

學生在兩組題目的練習比較中發現：先運用商不變規律把除數是小數的轉化成除數是整數的除法，再按除數是整數的除法的方法來計算。如果教師直接通過例題的教學就讓學生嘗試計算，學生將缺少再次理解算理的機會。所以“試一試”的教學為學生提供了自主遷移的機會，對學生更深刻地理解算理是十分必要的。在教學整數加法時要讓學生理解為什么要“相同數位要對齊，相加滿十向前一位進一”的算理。在教學小數加法時要讓學生理解為什么小數點要對齊。

學生在理解算理的基礎上，應該引導學生自主探索和合作交流出算法的多樣化，教師還應該引領學生尋找最優的算法。以這種算法為主進行訓練，從而提高學生的計算能力。

例如：在兩位數乘整十數探索“24×10”的口算方法時，有的學生聯系情境圖，先算9箱有多少瓶：24×9=216，再加1箱的24瓶：216+24=240；先算5箱有多少瓶：24×

5=120再算10箱有多少瓶：120×2=240；把每箱中的24瓶分成20瓶和4瓶，先算10個20瓶是200瓶，再算4個10瓶是40瓶，再用200+40=240；還有利用24×1=24遷移出24×10=240。在發散的基礎上引導學生著重理解最后一種算法“24乘1個十得24個十就是240”，在比較中引領學生進行算法的優化，在練習中重點運用這種算法，從而提高學生的計算能力。

為了計算簡便，解題中要訓練學生合理運用運算定律，靈活解題，增強計算技巧。如計算3.4×0.125 +4.6×0.125，學生一眼就能看出運用乘法分配律可以得出（3.4+4.6）×0.125。教學時，教師不應就此滿足，可進一步深化，充分挖掘學生的潛能，如依次出示：1.25× 0.34 + 4.6 × 0.125 3.4 ÷8+ 4.6 × 0.125 這樣，學生也就不會一遇到稍有變化的題目就不會解，同時學生的思維也得到了訓練。教學中要減少學生計算的錯誤，提高計算的正確率，應根據學生的實際情況，因材施教，因人施教，采取相應的對策，才能提高學生計算的能力。

四、加強練習和基本技能訓練（精講精練）

教師設計練習時最好分層進行，形式多樣。特別是練習的內容要注意有針對性，有層次，有梯度，練習的形式要多樣，學生在進行計算練習時才不會覺得枯燥，才會覺得有興趣。在設計練習題時要注意圍繞重點與難點來設計一些有針對性的練習，盡量讓學生能夠練習有所收獲。但要注意練習的數量要有個度，不能只要量不講質，搞題海戰術，就會適得其反。部分學生本身缺乏勤奮學習的精神，再加上計算本身又枯燥乏味，缺乏情節，學生遇到題量較多時，易產生抵觸情緒，不愿計算，嚴重的可影響學生對學習數學的興趣，教學中，首先要注意對題量的控制，其次計算形式要多樣，除了計算題，可適當增加一些判斷、選擇題，趣味題，如：（）+（）－（）+（）=0。這樣既減輕了學生的學業負擔，又增加了學生的學習興趣。

1、加強口算基本訓練,提高學生的計算能力

任何一道四則混合運算題都是有若干道口算題綜合而成的,口算的正確、迅速與否

直接關系到計算能力的提高。設計口算練習時,要有針對性,由易到難,逐步提高,包括一些簡便運算題,經常進行口算練習,做到每堂課練習5分鐘，有利于提高學生的計算能力。

2、加強簡算訓練，提高計算效率。簡便計算是小學計算教學的重要組成部分，它要求學生充分運用學過的運算定律、性質、公式，合理改變運算的數據及運算順序，使計算盡可能簡便、快捷，提高計算效率。因此，在教學中，必須加強簡算訓練，逐步增強簡算意識，提高簡算能力。

計算中，學生容易套用、濫用一些性質、定律，要讓學生進行一些對比練習，自己診斷錯誤，反思計算出錯的癥結點，防止再次出現同樣的錯誤。如：300－175+25，300－175－25；125×8÷125×8，（125×8）÷（125×8）；??讓學生辨析，什么情況下運用性質、定律可以簡便，明白為什么有些可以用，有些題不能用。

3、加強估算教學,提高學生的計算能力

估算可以培養學生的“數感”,可以引導學生深入理解“運算”,可以幫助學生檢查計算的結果正確與否,運用估算的方法可以對計算的結果做預先定位,快速地確定計算結果的取值范圍,通過計算前的估算和計算后的檢查,可以避免由于粗心大意造成的錯誤。如:個位是2和7的兩個數,計算結果的個位如果相加就肯定是9,相乘就一定是4。

4、重視錯題分析,提高學生的計算能力

針對計算中存在的問題和一些常見而又典型的錯例,教師要與學生一起分析、交流,通過集體“會診”,達到既“治病”又“預防”的目的。至于那些形近又易錯的題目,教師應該組織學生經常進行對比練習,培養學生比較、鑒別的能力。

五、培養學生良好的學習習慣,提高學生的計算能力

學生的計算錯誤,表面上看是“粗心”造成的,可是這些“粗心”其實就是沒有養成良好的學習習慣。所以在平時的計算訓練時,要讓學生養成一對、二審、三算、四查的好習慣。

1.對:就是認真核對。題目都抄錯了,結果又怎么能正確呢?所以,要求學生在抄題和每步計算時,都應當及時與原題或上一步算式進行核對,以免抄錯數或運算符號。

2.審:就是認真審題。引導學生在做計算題時,不應拿起筆來就下手算,必須先審題,弄清這道題應該先算什么,再算什么,有沒有簡便的計算方法,然后才能動筆算。

3.算:就是認真書寫、計算。作業、練習的書寫都要工整,不能潦草,格式一定要規范,對題目中的數字、小數點、運算符號的書寫尤其要符合規范,數字間有適當的間隔,草稿上的豎式也要數位對齊、條理清楚,計算時精力集中,不急不搶。另外,計算必須先求準,再求快。

4.查:就是仔細檢查。計算完,首先要檢查計算方法是不是合理;其次,檢查數字、符號會不會抄錯,小數點會不會錯寫或漏寫;再次,對計算中途得到的每一個得數和最后的結果都要進行檢查和驗算。

總之，在計算教學中也要關注過程，讓學生在理解算理的基礎上掌握方法，掌握計算法則，優化算法。這樣學生經歷了過程，明白了道理，對學生的計算技能和創新能力的培養有好處，有利于學生的后續發展。

第三篇：MATLAB計算24點游戲代碼

clear,closeall clc a=5;b=7;c=10;

d=4;%這里輸入需要計算的四個數字a,b,c,d f=[a b c d];tic;g=perms(f);[m,n]=size(g);h='+-*/';fori=1:24 for k1=1:4 for k2=1:4 for k3=1:4

str11=[num2str(g(i,1)),h(k1),num2str(g(i,2)),h(k2),num2str(g(i,3)),h(k3),num2str(g(i,4))];

str22=['(',num2str(g(i,1)),h(k1),num2str(g(i,2)),')',h(k2),num2str(g(i,3)),h(k3),num2str(g(i,4))];

str33=['(',num2str(g(i,1)),h(k1),num2str(g(i,2)),h(k2),num2str(g(i,3)),')',h(k3),num2str(g(i,4))];

str44=['(',num2str(g(i,1)),h(k1),num2str(g(i,2)),')',h(k2),'(',num2str(g(i,3)),h(k3),num2str(g(i,4)),')',];A=str2num(str11);B=str2num(str22);C=str2num(str33);D=str2num(str44);if A==24||B==24||C==24||D==24 break else end end

if A==24||B==24||C==24||D==24 break else end end

if A==24||B==24||C==24||D==24 break else end end

if A==24||B==24||C==24||D==24 break else end end

if A==24

answer=str11;elseif B==24

answer=str22;elseif C==24

answer=str33;elseif D==24

answer=str44;else

answer='無解';end

disp(['計算方法',num2str(answer)])time=toc;

disp(['計算耗時',num2str(time),'s'])

第四篇：matlab計算AHP層次分析法

用matlab解決層次分析法AHP

1、求矩陣最大特征值及特征向量 用matlab求：

輸入：A=[1 1/2 2 1/4;2 1 1 1/3;1/2 1 1 1/3;4 3 3 1]

[x,y]=eig(A)得出：特征向量x=[0.2688 0.3334 0.2373 0.8720]

最大特征值λmax=4.1964

2、一致性檢驗

CI=（λmax-n）/(n-1)=(4.1964-4)/(4-1)=0.0655 CR=CI/RI=0.0655/0.9=0.0727

(注：維數為4時，RI=0.9)CR=0.0727<0.1,矩陣一致性通過檢驗

3、對最大特征值進行歸一化處理，即可得到各指標權重（歸一化：分項/分項之和）W=[0.157 0.195 0.139 0.510]

第五篇：用matlab電力系統潮流計算

題目：潮流計算與matlab

教學單位 電氣信息學院

姓 名 學 號

年 級

專 業 電氣工程及其自動化

指導教師

職 稱 副教授

摘 要

電力系統穩態分析包括潮流計算和靜態安全分析。本文主要運用的事潮流計算，潮流計算是電力網絡設計與運行中最基本的運算，對電力網絡的各種設計方案及各種運行方式進行潮流計算，可以得到各種電網各節點的電壓，并求得網絡的潮流及網絡中的各元件的電力損耗，進而求得電能損耗。本位就是運用潮流計算具體分析，并有MATLAB仿真。

關鍵詞： 電力系統 潮流計算 MATLAB

Abstract Electric power system steady flow calculation and analysis of the static safety analysis.This paper, by means of the calculation, flow calculation is the trend of the power network design and operation of the most basic operations of electric power network, various design scheme and the operation ways to tide computation, can get all kinds of each node of the power grid voltage and seek the trend of the network and the network of the components of the power loss, and getting electric power.The standard is to use the power flow calculation and analysis, the specific have MATLAB simulation.Key words: Power system;Flow calculation;MATLAB simulation

目 錄 任務提出與方案論證....................................................................................................................................2 2 總體設計........................................................................................................................................................3 2.1潮流計算等值電路.............................................................................................................................3 2.2建立電力系統模型.............................................................................................................................3 2.3模型的調試與運行.............................................................................................................................3 3 詳細設計........................................................................................................................................................4 3.1 計算前提............................................................................................................................................4 3.2手工計算.............................................................................................................................................7 4設計圖及源程序...........................................................................................................................................11 4.1MATLAB仿真.......................................................................................................................................11 4.2潮流計算源程序...............................................................................................................................11 5 總結.............................................................................................................................................................31 參考文獻..........................................................................................................................................................32 任務提出與方案論證

潮流計算是在給定電力系統網絡結構、參數和決定系統運行狀態的邊界條件的情況下確定系統穩態運行狀態的一種基本方法,是電力系統規劃和運營中不可缺少的一個重要組成部分。可以說,它是電力系統分析中最基本、最重要的計算,是系統安全、經濟分析和實時控制與調度的基礎。常規潮流計算的任務是根據給定的運行條件和網路結構確定整個系統的運行狀態，如各母線上的電壓（幅值及相角）、網絡中的功率分布以及功率損耗等。潮流計算的結果是電力系統穩定計算和故障分析的基礎。在電力系統運行方式和規劃方案的研究中，都需要進行潮流計算以比較運行方式或規劃供電方案的可行性、可靠性和經濟性。同時，為了實時監控電力系統的運行狀態，也需要進行大量而快速的潮流計算。因此，潮流計算是電力系統中應用最廣泛、最基本和最重要的一種電氣運算。在系統規劃設計和安排系統的運行方式時，采用離線潮流計算；在電力系統運行狀態的實時監控中，則采用在線潮流計算。是電力系統研究人員長期研究的一個課題。它既是對電力系統規劃設計和運行方式的合理性、可靠性及經濟性進行定量分析的依據，又是電力系統靜態和暫態穩定計算的基礎。

潮流計算經歷了一個由手工到應用數字電子計算機的發展過程，現在的潮流算法都以計算機的應用為前提用計算機進行潮流計算主要步驟在于編制計算機程序，這是一項非常復雜的工作。對系統進行潮流分析，本文利用 MATLAB中的SimpowerSystems工具箱設計電力系統，在simulink 環境下，不僅可以仿真系統的動態過程，還可以對系統進行穩態潮流分析。

總體設計

SimpowerSystems使用Simulink環境，可以將該系統中的發電機、變壓器，線路等模型聯結起來，形成電力系統仿真模擬圖。在加人測量模塊，并對各元件的參數進行設置后，用measurement和sink中的儀器可以觀察各元件的電壓、電流、功率的大小。

2.1潮流計算等值電路

10MWYN,d114?63MWVA15MWGGGG120MW10kV??p0?15.7kW??ps?73kW?I%?0.5?0?Vs%?10.5YN,d1116MWVA4?63MW“?xd?0.134?x2?0.161?x?0.06?0?cos?N?0.8510kV??p0?11kW??ps?50kW?I%?0.55?0?Vs%?10.5YN,d112?10MWVA35kV32km25MW110kV80km25MW110kV70km110kVYN,d112?20MWVA20MWGGG4?15MW”?xd?0.136?x2?0.16?x?0.073?0?cos?N?0.8??p0?18.6kW??ps?89kW?I%?0.530MW?0?Vs%?10.510kV??p0?15.7kW??ps?73kW?I%?0.5GG?0V%?10.5s?YN,d112?16MWVA63MWVA??p0?44kW??ps?121kW10kV?I%?0.35?0?Vs%?10.5GG35MWYN,Y,d112?10MVA10kVGGG3?12MW1?50MW“?xd?0.128?x2?0.154?x?0.054?0?cos?N?0.852?25MW”?xd?0.128?x2?0.157?x?0.0591?0?cos?N?0.8?x?0.136?x2?0.161?x?0.075?0?cos?N?0.8“d

2.2建立電力系統模型

在Simulink中按照電力系統原型選擇元件進行建模。所建立的模型和建立的方法在詳細設計中詳述。

在電力系統模型的建立工程中主要涉及到的是：元器件的選擇及其參數的設置；發電機選型；變壓器選擇；線路的選擇；負荷模型的選擇；母線選擇。

2.3模型的調試與運行

建立系統模型，并設置好參數以后，就可以在Simulink環境下進行仿真運行。運行的具體結果和分析也在詳細設計中詳述。

30km35kV0km31??p0?44kW??ps?121kW?I%?0.35?0?Vs%?10.5??p0?13.2kW??ps?63kW?I0%?0.55?V%?10.5s(1?2)80MW?Vs(2?3)%?6.55MW??Vs(1?3)%?17.53 詳細設計

3.1 計算前提

首先是發電機的參數計算，先對5個發電廠簡化為5臺發電機來計算。發電機G1：

P1?4?15?60MWQ1?60?tan(arccos0.8)?45MVar發電機G2：

P2?4?63?252MWQ2?252?tan(arccos0.85)?156MVar發電機G3：

P3?3?12?36MWQ3?36?tan(arccos0.8)?27MVar發電機G4：

P4?1?50?50MWQ4?50?tan(arccos0.85)?31MVar發電機G5：

P5?2?25?50MWQ5?50?tan(arccos0.8)?37.5MVar

其次是變電站的參數計算，我們還是對7個變電站簡化為7臺變壓器來計算。變壓器T1：

2?ps?VN73?11023RT1??10??103?3.450?232SN(16?10)2Vs%?VN10.5?1102XT1??10??10?79.406?SN16?103?S01??p0?j變壓器T2：（雙并聯）

I0%?SN?(0.0157?j0.0800)MVA 100RT2XT221?ps?VN189?11023???10???103?1.346?2322SN2(20?10)21Vs%?VN110.5?1102???10???10?31.7625? 32SN220?10?S02?2?(?p0?j變壓器T3：（四并聯）

I0%?SN)?(0.0372?j0.2000)MVA100 1?ps?VN21121?11023RT3???10???103?0.092?2324SN4(63?10)XT31Vs%?VN2110.5?1102???10???10?5.042? 4SN463?103I0%?SN)?(0.1760?j0.8820)MVA100?S03?4?(?p0?j變壓器T4：（雙并聯）

1RT1?1.7250?21 XT4?XT1?39.7030?2?S04?2?S01?(0.0314?j0.1600)MVART4?變壓器T5：

RT5?4RT3?0.3680?XT5?4XT3?20.168??S05?1?S03?(0.0440?j0.2205)MVA4163?3523???10?0.386? 322(10?10)

變壓器T6：（兩個三繞組變壓器并聯）

RT6?1?RT6?2?RT6?31?[Vs(1?2)%?Vs(1?3)%?Vs(2?3)%]?10.7521Vs2%??[Vs(1?2)%?Vs(2?3)%?Vs(1?3)%]??0.25

21Vs3%??[Vs(1?3)%?Vs(2?3)%?Vs(1?2)%]?6.752Vs1%?21Vs1%?VNXT6?1???10?6.584?2SNXT6?2XT6?321Vs2%?VN???10??0.153?2SN21Vs3%?VN???10?4.134? 2SN?S06?2?(?P06?j變壓器T7：（雙并聯）

I0%?10)?(0.0264?j0.1100)MVA 100 RT7XT721?ps?VN150?3523???10???103?0.306?2322SN2(10?10)21Vs%?VN110.5?352???10???10?6.431?2SN210?103

?S07?2?(?p0?jI0%?SN)?(0.0220?j0.1100)MVA100再次是傳輸線參數計算，5條傳輸線的具體計算如下。

根據教材查得r0?0.21?/km x0?0.4?/km b0?2.8?10S/km ?6線路L1：

線路L2：

線路L3：（雙回路）

線路L4：

線路L5：（雙回路）RL1?r0?l1?0.21?40?8.4?XL1?x0?l1?0.4?40?16?B?6L1?b0?l1?2.8?10?40?1.12?10?4S ?Q1L1??2BL1V2N??0.6776MVarRL2?r0?l2?0.21?130?27.3?XL2?x0?l2?0.4?130?52?B?6L2?b0?l2?2.8?10?130?3.64?10?4S ?Q1L2??2BL2V2N??2.2022MVarR?12?rl1L30?3?2?0.21?70?7.35?X11L3?2?x0?l3?2?0.4?70?14? BL3?2?b?40?l3?2?2.8?10?6?70?3.92?10S?Q1L3??2B2L3VN??2.3716MVarRL4?r0?l4?0.21?60?12.6?XL4?x0?l4?0.4?60?24?BL4?b0?l4?2.8?10?6?60?1.68?10?4S ?Q12L4??2BL4VN??1.0164MVar

11RL5??r0?l5??0.21?20?2.1?2211XL5??x0?l5??0.4?20?4? 22BL5?2?b0?l5?2?2.8?10?6?20?1.12?10?4S1?QL5??BL3VN2??0.0686MVar23.2手工計算

FLR1：

P2102?ST1?2(RT1?jXT1)?(3.450?j74.406)?(0.0285?j0.6562)MVA2VN110Sa?10MW??ST1??S01?j?QL1?(10.0442?j0.1142)MVAP2?Q210.04422?0.11422?SL1?(RL1?jXL1)?(8.4?j16)?(0.070?j0.1334)MVAVN21102?ST2P2?Q2402?452?(RT2?jXT2)?(1.346?j31.7625)?(0.4032?j9.5156)MVAVN21102FLR2Sb?SG1?20??ST2?60?j45?20?0.4032?j9.5156?(39.5968?j35.4844)MVASc?Sb?Sa?25?jQL1??SL1?(4.4826?j35.9144)MVA：

?ST3P2?Q22522?1562?(RT3?jXT3)?(0.092?j5.042)?(0.6679?j36.6024)MVA22VN110P?Q4.4931?34.1048(R?jX)?(27.3?j52)?(2.67?j5.0854)MVAL2L222VN1102222Sc'?(4.4931?j34.1048)MVA?SL2?FLRSd?SG2?Sc'?120??ST3??S03?jQL2??SL2?(132.9792?j149.229)MVA3：

?ST4P2?Q262?272?(RT4?jXT4)?(1.725?j39.703)?(0.1091?j2.5101)MVAVN21102P2?Q2133.59552?149.99562?(RL3?jXL3)?(7.35?j14)?(24.51?j46.682)MVA22VN110'Sd?(133.5955?j149.9956)MVA?SL3'Se?SG3?Sd?30?25??ST4??S04?jQL3??SL3?(89.945?j130.0151)MVAFLR4： ?ST5P2?Q2502?312?(RT5?jXT5)?(0.368?j20.168)?(0.1052?j5.7687)MVA22VN110P2?Q292.74872?133.99372?(RL4?jXL4)?(12.6?j24)?(27.654?j52.674)MVA22VN110Se'?(92.7481?j133.9937)MVA?SL4Sf?SG4?Se'?80??ST5??S05?jQL4??SL4?(34.9449?j107.3469)MVAFLR5: 152?ST7?2?(0.306?j6.431)?(0.0562?j1.1812)MVA35Sh?15??ST7??S07?j?QL5?(15.0782?j0.3422)MVA15.07822?0.34222?SL5??(2.1?j4)?(0.3899?j0.743)MVA352Si?Sh??SL5??S06?j?QL5?5?(20.4945?j1.1266)MVA 152?37.52?ST6?3??(0.386?j4.34)?(0.514?j5.7793)MVA35220.65052?0.54512?ST6?2??(0.386?j0.153)?(0.1345?j0.0533)MVA23526.3362?98.73692?ST6?1??(0.386?j6.584)?(3.2905?j56.1256)MVA352Sg?Sf??ST6?1?SG5??ST6?2??ST6?3?Si?35?(25.5114?j194.12)MVA計算每一個FLR的功率分布和電壓分布計算如下： FLR1：

?VT2?PR?QX40?1.346?45?31.7625??12.8970kVVN115 Vb?115??VT2?102.1030kVPR?QX10.0442?8.4?0.1442?16?VL1???0.8489kVVb102.1030Va?Vb??VL1?101.2541kVFLR2：

功率分布：

SL2?ZZ?Z*T3*L2*Sd?T3(Vb?VN)Z?ZL2****?VN?(0.092?j5.042)?(132.9792?j149.229)?1418.6727.392?j57.042T3?(4.8812?j13.8097)MVAST3?ZZ?Z*L2*L2*Sd?T3(Vb?VN)Z?ZL2?VN?(27.3?j52)?(132.9792?j149.229)?1418.6727.392?j57.042T3?(108.687?j122.62)MVA 電壓分布：

Sc1?SL2??SL2?(4.8812?j13.8097)?(2.67?j5.0854)?(7.5512?j8.7243)MVA7.5512?27.3?8.7243?52??2.424kV102.1030Vd?Vb??VL2?102.103?(?2.424)?104.527kV?VL2?功率分布：

FLR3：

SL3?ZZ?Z*T4*L3*Se?T4(VG3?Vd)Z?ZL3****?VN?(1.725?j39.703)?(89.945?j130.0151)?1037.9279.075?j53.73T4?(59.444?j16.846)MVAST4?ZZ?Z*L3*L3*Se?T4(Vb?VN)Z?ZL3?VN?(7.35?j14)?(89.945?j130.0151)?1037.9279.075?j53.73T4?(31.811?j60.1256)MVA 電壓分布：

Se1?SL3??SL3?(59.444?j19.846)?(24.51?j46.682)?(83.954?j26.836)MVA83.954?7.35?26.836?14?9.404kV105.5643Ve?Vd??VL3?96.16kV?VL3?功率分布：

FLR4：

SL4?ZZ?Z*T5*L4*Sf?T5(VG3?Vd)Z?ZL4****?VN=(0.368?j20.168)?(34.9449?j107.3469)?1037.92712.968?j44.168T5?(20.843?j19.689)MVAST4?ZZ?Z*L4*L4*Sf?T5(VG3?Vd)Z?ZL4?VN=(12.6?j24)?(34.9449?j107.3469)?1037.92712.968?j44.168T5?(1.398?j44.389)MVA 電壓分布：

Se1?SL3??SL3?(59.444?j16.846)?(24.51?j46.682)?(83.954?j63.528)MVA83.954?12.6?63.528?24?24.464kV105.5643Ve?Vd??VL3?81.10kV?VL4?FLR5: 這里我們先將f點和發電機G5當做電源，經過ZT61和ZT63構成兩端供電網絡以g點作為運算負荷進行計算。ST6?ST4(0.386?j4.134)?(20.2656?j70.9293)?(22.0938?37)?35?(3.900?j25.1175)MVA0.772?j10.718(0.386?j6.584)?(20.2656?j70.9293)?(22.0938?37)?35??(16.5061?j91.7905)MVA0.772?j10.718電壓分布：

ST631?ST63??ST63?(16.6421?j97.5698)MVA16.6421?0.386?97.5698?4.134?10.9186kV37Vg?37??VT63?26.0814V?VT63?20.2656?0.386?70.9293?(?0.153)??0.1162kV

26.0814Vi?Vg??VT62?26.1976?VT62?20.4945?2.1?1.1266?4?1.815kV26.1976Vh?Vi??VL5?24.3826?VL5?

4設計圖及源程序

4.1MATLAB仿真

相關的原始數據輸入格式如下：

1、B1是支路參數矩陣，第一列和第二列是節點編號。節點編號由小到大編寫。

2、對于含有變壓器的支路，第一列為低壓側節點編號，第二列為高壓側節點編號，將變壓器的串聯阻抗置于低壓側處理，第三列為支路的串列阻抗參數，第四列為支路的對地導納參數，第五烈為含變壓器支路的變壓器的變比，第六列為變壓器是否是否含有變壓器的參數，其中“1”為含有變壓器，“0”為不含有變壓器。

3、B2為節點參數矩陣，其中第一列為節點注入發電功率參數；第二列為節點負荷功率參數；第三列為節點電壓參數；第六列為節點類型參數，其中“1”為平衡節點，“2”為PQ節點，“3”為PV節點參數。

4、X為節點號和對地參數矩陣。其中第一列為節點編號，第二列為節點對地參數。

4.2潮流計算源程序

%本程序的功能是用牛頓——拉夫遜法進行11節點潮流計算 clear;n=11;%input('請輸入節點數：n=');nl=11;%input('請輸入支路數：nl=');isb=1;%input('請輸入平衡母線節點號：isb=');pr=0.00001;%input('請輸入誤差精度：pr=');B1=[1

0.03512+0.08306i

0.13455i

0;

0.0068+0.18375i

0

1.02381

1;

0.05620+0.13289i

0.05382i

0;

0.00811+0.24549i

0

1.02381

1;

0.05620+0.13289i

0.05382i

0;

0.04215+0.09967i

0.04037i

0;

0.0068+0.18375i

0

1.02381

1;

0.02810+0.06645i

0.10764i

0;

0.05620+0.13289i

0.05382i

0;0.00811+0.24549i

0

1;

0.03512+0.08306i

0.13455i

0] B2=[0

0

1.1

1.1

0

1;

0

0

0

0

2;

0

0.343+0.21256i

0

0

2;

0

0

0

0

2;

0

0.204+0.12638i

0

0

2;

0

0

0

0

2;

0

0.306+0.18962i

0

0

2;

0

0

0

0

2;

0.5

0

1.1

1.1

0

3;

0

0.343+0.21256i

0

0

2;

0

0

0

0

2];% B1矩陣：

1、支路首端號；

2、末端號；

3、支路阻抗；

4、支路對地電納 %

5、支路的變比；

6、支路首端處于K側為1，1側為0 % B2矩陣：

1、該節點發電機功率；

2、該節點負荷功率；

3、節點電壓初始值 %

4、PV節點電壓V的給定值；

5、節點所接的無功補償設備的容量 %

6、節點分類標號：1為平衡節點（應為1號節點）；2為PQ節點； %

3為PV節點；

%input('請輸入各節點參數形成的矩陣： B2=');Y=zeros(n);e=zeros(1,n);f=zeros(1,n);V=zeros(1,n);sida=zeros(1,n);S1=zeros(nl);% % %--------------------for i=1:nl

%支路數

if B1(i,6)==0

%左節點處于1側

p=B1(i,1);q=B1(i,2);

else

%左節點處于K側

p=B1(i,2);q=B1(i,1);

end

Y(p,q)=Y(p,q)-1./(B1(i,3)*B1(i,5));%非對角元

Y(q,p)=Y(p,q);

%非對角元

Y(q,q)=Y(q,q)+1./(B1(i,3)*B1(i,5)^2)+B1(i,4)./2;%對角元K側

Y(p,p)=Y(p,p)+1./B1(i,3)+B1(i,4)./2;

%對角元1側

end %求導納矩陣

disp('導納矩陣 Y=');disp(Y)%---------------------------G=real(Y);B=imag(Y);

%分解出導納陣的實部和虛部

for i=1:n

%給定各節點初始電壓的實部和虛部

e(i)=real(B2(i,3));

f(i)=imag(B2(i,3));

V(i)=B2(i,4);

%PV節點電壓給定模值

end for i=1:n

%給定各節點注入功率

S(i)=B2(i,1)-B2(i,2);

%i節點注入功率SG-SL

B(i,i)=B(i,i)+B2(i,5);%i節點無功補償量

end %===================== P=real(S);Q=imag(S);

%分解出各節點注入的有功和無功功率 ICT1=0;IT2=1;N0=2*n;N=N0+1;a=0;%迭代次數ICT1、a；不滿足收斂要求的節點數IT2 while IT2~=0

% N0=2*n 雅可比矩陣的階數；N=N0+1擴展列

IT2=0;a=a+1;

for i=1:n

if i~=isb

%非平衡節點

C(i)=0;D(i)=0;

for j1=1: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*fj+Bij*ej)

end

P1=C(i)*e(i)+f(i)*D(i);%節點功率P計算eiΣ(Gij*ej-Bij*fj)+fiΣ(Gij*fj+Bij*ej)

Q1=C(i)*f(i)-e(i)*D(i);%節點功率Q計算fiΣ(Gij*ej-Bij*fj)-eiΣ(Gij*fj+Bij*ej)%求i節點有功和無功功率P',Q'的計算值

V2=e(i)^2+f(i)^2;

%電壓模平方

%========= 以下針對非PV節點來求取功率差及Jacobi矩陣元素 =========

if B2(i,6)~=3

%非PV節點

DP=P(i)-P1;

%節點有功功率差

DQ=Q(i)-Q1;

%節點無功功率差

%=============== 以上為除平衡節點外其它節點的功率計算 ================= %================= 求取Jacobi矩陣 ===================

for j1=1:n

if j1~=isb&j1~=i

%非平衡節點&非對角元

X1=-G(i,j1)*e(i)-B(i,j1)*f(i);% dP/de=-dQ/df

X2=B(i,j1)*e(i)-G(i,j1)*f(i);% dP/df=dQ/de

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;J(p,N)=DQ;m=p+1;

% X3=dQ/de J(p,N)=DQ節點無功功率差

J(m,q)=X1;J(m,N)=DP;q=q+1;

% X1=dP/de J(m,N)=DP節點有功功率差

J(p,q)=X4;J(m,q)=X2;

% X4=dQ/df X2=dp/df

elseif j1==i&j1~=isb %非平衡節點&對角元

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(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;J(p,N)=DQ;%擴展列△Q

m=p+1;

J(m,q)=X1;q=q+1;J(p,q)=X4;J(m,N)=DP;%擴展列△P

J(m,q)=X2;

end

end

else

%=============== 下面是針對PV節點來求取Jacobi矩陣的元素 ===========

DP=P(i)-P1;

% PV節點有功誤差

DV=V(i)^2-V2;

% PV節點電壓誤差

for j1=1:n

if j1~=isb&j1~=i

%非平衡節點&非對角元

X1=-G(i,j1)*e(i)-B(i,j1)*f(i);

% dP/de

X2=B(i,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;J(p,N)=DV;% PV節點電壓誤差

m=p+1;

J(m,q)=X1;J(m,N)=DP;q=q+1;J(p,q)=X6;

% PV節點有功誤差

J(m,q)=X2;

elseif j1==i&j1~=isb %非平衡節點&對角元

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);

p=2*i-1;q=2*j1-1;J(p,q)=X5;J(p,N)=DV;% PV節點電壓誤差

m=p+1;

J(m,q)=X1;J(m,N)=DP;q=q+1;J(p,q)=X6;

% PV節點有功誤差

J(m,q)=X2;

end

end

end

end

end %========= 以上為求雅可比矩陣的各個元素及擴展列的功率差或電壓差 =====================

for k=3:N0

% N0=2*n（從第三行開始，第一、二行是平衡節點）

k1=k+1;N1=N;

% N=N0+1 即 N=2*n+1擴展列△P、△Q 或 △U

for k2=k1:N1

% 從k+1列的Jacobi元素到擴展列的△P、△Q 或 △U

J(k,k2)=J(k,k2)./(J(k,k)+eps);% 用K行K列對角元素去除K行K列后的非對角元素進行規格化

end

J(k,k)=1;

% 對角元規格化K行K列對角元素賦1

%==================== 回代運算

=======================================

if k~=3

% 不是第三行

k > 3

k4=k-1;

for k3=3:k4

% 用k3行從第三行開始到當前行的前一行k4行消去

for k2=k1:N1 %

k3行后各行上三角元素

J(k3,k2)=J(k3,k2)-J(k3,k)*J(k,k2);%消去運算（當前行k列元素消為0）

end

%用當前行K2列元素減去當前行k列元素乘以第k行K14 列元素

J(k3,k)=0;%當前行第k列元素已消為0

end

if k==N0

%若已到最后一行

break;

end

%================== 前代運算

==================================

for k3=k1:N0

% 從k+1行到2*n最后一行

for k2=k1:N1

% 從k+1列到擴展列消去k+1行后各行下三角元素

J(k3,k2)=J(k3,k2)-J(k3,k)*J(k,k2);%消去運算

end

%用當前行K2列元素減去當前行k列元素乘以第k行K2列元素

J(k3,k)=0;%當前行第k列元素已消為0

end

else

%是第三行k=3

%====================== 第三行k=3的前代運算 ========================

for k3=k1:N0

%從第四行到2n行（最后一行）

for k2=k1:N1

%從第四列到2n+1列（即擴展列）

J(k3,k2)=J(k3,k2)-J(k3,k)*J(k,k2);%消去運算（當前行3列元素消為0）

end

%用當前行K2列元素減去當前行3列元素乘以第三行K2列元素

J(k3,k)=0;%當前行第3列元素已消為0

end

end

end %====上面是用線性變換方式高斯消去法將Jacobi矩陣化成單位矩陣=====

for k=3:2:N0-1

L=(k+1)./2;

e(L)=e(L)-J(k,N);

%修改節點電壓實部

k1=k+1;

f(L)=f(L)-J(k1,N);

%修改節點電壓虛部

end

%------修改節點電壓-----------

for k=3:N0

DET=abs(J(k,N));

if DET>=pr

%電壓偏差量是否滿足要求

IT2=IT2+1;%不滿足要求的節點數加1

end

end

ICT2(a)=IT2;

%不滿足要求的節點數

ICT1=ICT1+1;

%迭代次數 end %用高斯消去法解”w=-J*V“

disp('迭代次數：');disp(ICT1);disp('沒有達到精度要求的個數：');disp(ICT2);for k=1:n

V(k)=sqrt(e(k)^2+f(k)^2);

%計算各節點電壓的模值

sida(k)=atan(f(k)./e(k))*180./pi;

%計算各節點電壓的角度

E(k)=e(k)+f(k)*j;

%將各節點電壓用復數表示 end %=============== 計算各輸出量 =========================== disp('各節點的實際電壓標幺值E為(節點號從小到大排列)：');disp(E);

%顯示各節點的實際電壓標幺值E用復數表示 disp('----------------------');disp('各節點的電壓大小V為(節點號從小到大排列)：');disp(V);

%顯示各節點的電壓大小V的模值 disp('----------------------');disp('各節點的電壓相角sida為(節點號從小到大排列)：');disp(sida);

%顯示各節點的電壓相角 for p=1:n

C(p)=0;

for q=1:n

C(p)=C(p)+conj(Y(p,q))*conj(E(q));%計算各節點的注入電流的共軛值

end

S(p)=E(p)*C(p);

%計算各節點的功率 S = 電壓 X 注入電流的共軛值 end disp('各節點的功率S為(節點號從小到大排列)：');disp(S);

%顯示各節點的注入功率

disp('----------------------');disp('各條支路的首端功率Si為(順序同您輸入B1時一致)：');for i=1:nl

p=B1(i,1);q=B1(i,2);

if B1(i,6)==0

Si(p,q)=E(p)*(conj(E(p))*conj(B1(i,4)./2)+(conj(E(p)*B1(i,5))...-conj(E(q)))*conj(1./(B1(i,3)*B1(i,5))));

Siz(i)=Si(p,q);

else

Si(p,q)=E(p)*(conj(E(p))*conj(B1(i,4)./2)+(conj(E(p)./B1(i,5))...-conj(E(q)))*conj(1./(B1(i,3)*B1(i,5))));

Siz(i)=Si(p,q);

end

disp(Si(p,q));

SSi(p,q)=Si(p,q);

ZF=['S(',num2str(p),',',num2str(q),')=',num2str(SSi(p,q))];

disp(ZF);

disp('----------------------');end disp('各條支路的末端功率Sj為(順序同您輸入B1時一致)：');

for i=1:nl

p=B1(i,1);q=B1(i,2);

if B1(i,6)==0

Sj(q,p)=E(q)*(conj(E(q))*conj(B1(i,4)./2)+(conj(E(q)./B1(i,5))...-conj(E(p)))*conj(1./(B1(i,3)*B1(i,5))));

Sjy(i)=Sj(q,p);

else

Sj(q,p)=E(q)*(conj(E(q))*conj(B1(i,4)./2)+(conj(E(q)*B1(i,5))...-conj(E(p)))*conj(1./(B1(i,3)*B1(i,5))));

Sjy(i)=Sj(q,p);

end

disp(Sj(q,p));

SSj(q,p)=Sj(q,p);

ZF=['S(',num2str(q),',',num2str(p),')=',num2str(SSj(q,p))];

disp(ZF);

disp('----------------------');end disp('各條支路的功率損耗DS為(順序同您輸入B1時一致)：');for i=1:nl

p=B1(i,1);q=B1(i,2);

DS(i)=Si(p,q)+Sj(q,p);

disp(DS(i));

DDS(i)=DS(i);

ZF=['DS(',num2str(p),',',num2str(q),')=',num2str(DDS(i))];

disp(ZF);

disp('----------------------');end

%本程序的功能是用牛頓——拉夫遜法進行10節點潮流計算 %本程序的功能是用牛頓——拉夫遜法進行潮流計算 clear;n=10;%input('請輸入節點數：n=');nl=10;%input('請輸入支路數：nl=');isb=1;%input('請輸入平衡母線節點號：isb=');pr=0.00001;%input('請輸入誤差精度：pr=');B1=[1

0.03512+0.08306i

0.13455i

0;

0.0068+0.18375i

0

1.02381

1;

0.05620+0.13289i

0.05382i

0;

0.00811+0.24549i

0

1.02381

1;

0.05620+0.13289i

0.05382i

0;

0.04215+0.09967i

0.04037i

0;

0.0068+0.18375i

0

1.02381

1;

0.02810+0.06645i

0.10764i

0;0.00811+0.24549i

0

1;

0.03512+0.08306i

0.13455i

0] B2=[0

0

1.1

1.1

0

1;

0

0

0

0

2;

0

0.343+0.21256i

0

0

2;

0

0

0

0

2;

0

0.204+0.12638i

0

0

2;

0

0

0

0

2;

0

0.306+0.18962i

0

0

2;

0

0

0

0

2;

0.5

0

1.1

1.1

0

3;

0

0.343+0.21256i

0

0

2];% B1矩陣：

1、支路首端號；

2、末端號；

3、支路阻抗；

4、支路對地電納 %

5、支路的變比；

6、支路首端處于K側為1，1側為0 % B2矩陣：

1、該節點發電機功率；

2、該節點負荷功率；

3、節點電壓初始值 %

4、PV節點電壓V的給定值；

5、節點所接的無功補償設備的容量 %

6、節點分類標號：1為平衡節點（應為1號節點）；2為PQ節點； %

3為PV節點；

%input('請輸入各節點參數形成的矩陣： B2=');Y=zeros(n);e=zeros(1,n);f=zeros(1,n);V=zeros(1,n);sida=zeros(1,n);S1=zeros(nl);% % %--------------------for i=1:nl

%支路數

if B1(i,6)==0

%左節點處于1側

p=B1(i,1);q=B1(i,2);

else

%左節點處于K側

p=B1(i,2);q=B1(i,1);

end

Y(p,q)=Y(p,q)-1./(B1(i,3)*B1(i,5));%非對角元

Y(q,p)=Y(p,q);

%非對角元

Y(q,q)=Y(q,q)+1./(B1(i,3)*B1(i,5)^2)+B1(i,4)./2;%對角元K側

Y(p,p)=Y(p,p)+1./B1(i,3)+B1(i,4)./2;

%對角元1側

end %求導納矩陣

disp('導納矩陣 Y=');disp(Y)%---------------------------G=real(Y);B=imag(Y);

%分解出導納陣的實部和虛部

for i=1:n

%給定各節點初始電壓的實部和虛部

e(i)=real(B2(i,3));

f(i)=imag(B2(i,3));

V(i)=B2(i,4);

%PV節點電壓給定模值

end for i=1:n

%給定各節點注入功率

S(i)=B2(i,1)-B2(i,2);

%i節點注入功率SG-SL

B(i,i)=B(i,i)+B2(i,5);%i節點無功補償量

end %===================== P=real(S);Q=imag(S);

%分解出各節點注入的有功和無功功率

ICT1=0;IT2=1;N0=2*n;N=N0+1;a=0;%迭代次數ICT1、a；不滿足收斂要求的節點數IT2 while IT2~=0

% N0=2*n 雅可比矩陣的階數；N=N0+1擴展列

IT2=0;a=a+1;

for i=1:n

if i~=isb

%非平衡節點

C(i)=0;D(i)=0;

for j1=1: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*fj+Bij*ej)

end

P1=C(i)*e(i)+f(i)*D(i);%節點功率P計算eiΣ(Gij*ej-Bij*fj)+fiΣ(Gij*fj+Bij*ej)

Q1=C(i)*f(i)-e(i)*D(i);%節點功率Q計算fiΣ(Gij*ej-Bij*fj)-eiΣ(Gij*fj+Bij*ej)%求i節點有功和無功功率P',Q'的計算值

V2=e(i)^2+f(i)^2;

%電壓模平方

%========= 以下針對非PV節點來求取功率差及Jacobi矩陣元素 =========

if B2(i,6)~=3

%非PV節點

DP=P(i)-P1;

%節點有功功率差

DQ=Q(i)-Q1;

%節點無功功率差

%=============== 以上為除平衡節點外其它節點的功率計算 ================= %================= 求取Jacobi矩陣 ===================

for j1=1:n

if j1~=isb&j1~=i

%非平衡節點&非對角元

X1=-G(i,j1)*e(i)-B(i,j1)*f(i);% dP/de=-dQ/df

X2=B(i,j1)*e(i)-G(i,j1)*f(i);% dP/df=dQ/de

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;J(p,N)=DQ;m=p+1;

% X3=dQ/de J(p,N)=DQ節點無功功率差

J(m,q)=X1;J(m,N)=DP;q=q+1;

% X1=dP/de J(m,N)=DP節點有功功率差

J(p,q)=X4;J(m,q)=X2;

% X4=dQ/df X2=dp/df

elseif j1==i&j1~=isb %非平衡節點&對角元

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(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;J(p,N)=DQ;%擴展列△Q

m=p+1;

J(m,q)=X1;q=q+1;J(p,q)=X4;J(m,N)=DP;%擴展列△P

J(m,q)=X2;

end

end

else

%=============== 下面是針對PV節點來求取Jacobi矩陣的元素

===========

DP=P(i)-P1;

% PV節點有功誤差

DV=V(i)^2-V2;

% PV節點電壓誤差

for j1=1:n

if j1~=isb&j1~=i

%非平衡節點&非對角元

X1=-G(i,j1)*e(i)-B(i,j1)*f(i);

% dP/de

X2=B(i,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;J(p,N)=DV;% PV節點電壓誤差

m=p+1;

J(m,q)=X1;J(m,N)=DP;q=q+1;J(p,q)=X6;

% PV節點有功誤差

J(m,q)=X2;

elseif j1==i&j1~=isb %非平衡節點&對角元

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);

p=2*i-1;q=2*j1-1;J(p,q)=X5;J(p,N)=DV;% PV節點電壓誤差

m=p+1;

J(m,q)=X1;J(m,N)=DP;q=q+1;J(p,q)=X6;

% PV節點有功誤差

J(m,q)=X2;

end

end

end

end

end %========= 以上為求雅可比矩陣的各個元素及擴展列的功率差或電壓差 =====================

for k=3:N0

% N0=2*n（從第三行開始，第一、二行是平衡節點）

k1=k+1;N1=N;

% N=N0+1 即 N=2*n+1擴展列△P、△Q 或 △U

for k2=k1:N1

% 從k+1列的Jacobi元素到擴展列的△P、△Q 或 △U

J(k,k2)=J(k,k2)./J(k,k);% 用K行K列對角元素去除K行K列后的非對角元素進行規格化

end

J(k,k)=1;

% 對角元規格化K行K列對角元素賦1

%==================== 回代運算

=======================================

if k~=3

% 不是第三行

k > 3

k4=k-1;

for k3=3:k4

% 用k3行從第三行開始到當前行的前一行k4行消

去

for k2=k1:N1 %

k3行后各行上三角元素

J(k3,k2)=J(k3,k2)-J(k3,k)*J(k,k2);%消去運算（當前行k列元素消為0）

end

%用當前行K2列元素減去當前行k列元素乘以第k行K2列元素

J(k3,k)=0;%當前行第k列元素已消為0

end

if k==N0

%若已到最后一行

break;

end

%================== 前代運算

==================================

for k3=k1:N0

% 從k+1行到2*n最后一行

for k2=k1:N1

% 從k+1列到擴展列消去k+1行后各行下三角元素

J(k3,k2)=J(k3,k2)-J(k3,k)*J(k,k2);%消去運算

end

%用當前行K2列元素減去當前行k列元素乘以第k行K2列元素

J(k3,k)=0;%當前行第k列元素已消為0

end

else

%是第三行k=3

%====================== 第三行k=3的前代運算 ========================

for k3=k1:N0

%從第四行到2n行（最后一行）

for k2=k1:N1

%從第四列到2n+1列（即擴展列）

J(k3,k2)=J(k3,k2)-J(k3,k)*J(k,k2);%消去運算（當前行3列元素消為0）

end

%用當前行K2列元素減去當前行3列元素乘以第三行K2列元素

J(k3,k)=0;%當前行第3列元素已消為0

end

end

end %====上面是用線性變換方式高斯消去法將Jacobi矩陣化成單位矩陣=====

for k=3:2:N0-1

L=(k+1)./2;

e(L)=e(L)-J(k,N);

%修改節點電壓實部

k1=k+1;

f(L)=f(L)-J(k1,N);

%修改節點電壓虛部

end

%------修改節點電壓-----------

for k=3:N0

DET=abs(J(k,N));

if DET>=pr

%電壓偏差量是否滿足要求

IT2=IT2+1;%不滿足要求的節點數加1

end

end

ICT2(a)=IT2;

%不滿足要求的節點數

ICT1=ICT1+1;

%迭代次數 end %用高斯消去法解”w=-J*V“ disp('迭代次數：');disp(ICT1);disp('沒有達到精度要求的個數：');disp(ICT2);for k=1:n

V(k)=sqrt(e(k)^2+f(k)^2);

%計算各節點電壓的模值

sida(k)=atan(f(k)./e(k))*180./pi;

%計算各節點電壓的角度

E(k)=e(k)+f(k)*j;

%將各節點電壓用復數表示 end %=============== 計算各輸出量 =========================== disp('各節點的實際電壓標幺值E為(節點號從小到大排列)：');disp(E);

%顯示各節點的實際電壓標幺值E用復數表示 disp('----------------------');disp('各節點的電壓大小V為(節點號從小到大排列)：');disp(V);

%顯示各節點的電壓大小V的模值 disp('----------------------');disp('各節點的電壓相角sida為(節點號從小到大排列)：');disp(sida);

%顯示各節點的電壓相角 for p=1:n

C(p)=0;

for q=1:n

C(p)=C(p)+conj(Y(p,q))*conj(E(q));%計算各節點的注入電流的共軛值

end

S(p)=E(p)*C(p);

%計算各節點的功率 S = 電壓 X 注入電流的共軛值 end disp('各節點的功率S為(節點號從小到大排列)：');disp(S);

%顯示各節點的注入功率

disp('----------------------');disp('各條支路的首端功率Si為(順序同您輸入B1時一致)：');for i=1:nl

p=B1(i,1);q=B1(i,2);

if B1(i,6)==0

Si(p,q)=E(p)*(conj(E(p))*conj(B1(i,4)./2)+(conj(E(p)*B1(i,5))...-conj(E(q)))*conj(1./(B1(i,3)*B1(i,5))));

Siz(i)=Si(p,q);

else

Si(p,q)=E(p)*(conj(E(p))*conj(B1(i,4)./2)+(conj(E(p)./B1(i,5))...-conj(E(q)))*conj(1./(B1(i,3)*B1(i,5))));

Siz(i)=Si(p,q);

end

disp(Si(p,q));

SSi(p,q)=Si(p,q);

ZF=['S(',num2str(p),',',num2str(q),')=',num2str(SSi(p,q))];

disp(ZF);

disp('----------------------');end disp('各條支路的末端功率Sj為(順序同您輸入B1時一致)：');for i=1:nl

p=B1(i,1);q=B1(i,2);

if B1(i,6)==0

Sj(q,p)=E(q)*(conj(E(q))*conj(B1(i,4)./2)+(conj(E(q)./B1(i,5))...-conj(E(p)))*conj(1./(B1(i,3)*B1(i,5))));

Sjy(i)=Sj(q,p);

else

Sj(q,p)=E(q)*(conj(E(q))*conj(B1(i,4)./2)+(conj(E(q)*B1(i,5))...-conj(E(p)))*conj(1./(B1(i,3)*B1(i,5))));

Sjy(i)=Sj(q,p);

end

disp(Sj(q,p));

SSj(q,p)=Sj(q,p);

ZF=['S(',num2str(q),',',num2str(p),')=',num2str(SSj(q,p))];

disp(ZF);

disp('----------------------');end disp('各條支路的功率損耗DS為(順序同您輸入B1時一致)：');for i=1:nl

p=B1(i,1);q=B1(i,2);

DS(i)=Si(p,q)+Sj(q,p);

disp(DS(i));

DDS(i)=DS(i);

ZF=['DS(',num2str(p),',',num2str(q),')=',num2str(DDS(i))];

disp(ZF);

disp('----------------------');end

%本程序的功能是用牛頓——拉夫遜法進行12節點潮流計算 %本程序的功能是用牛頓——拉夫遜法進行潮流計算 clear;n=12;%input('請輸入節點數：n=');nl=12;%input('請輸入支路數：nl=');isb=1;%input('請輸入平衡母線節點號：isb=');pr=0.00001;%input('請輸入誤差精度：pr=');B1=[1

0.03512+0.08306i

0.13455i

0;

0.0068+0.18375i

0

1.02381

1;

0.05620+0.13289i

0.05382i

0;

0.00811+0.24549i

0

1.02381

1;

0.05620+0.13289i

0.05382i

0;

0.04215+0.09967i

0.04037i

0;

0.0068+0.18375i

0

1.02381

1;

0.02810+0.06645i

0.10764i

0;

0.05620+0.13289i

0.05382i

0;0.00811+0.24549i

0

1;

0.03512+0.08306i

0.13455i

0;

0.03512+0.08306i

0.13455i

0] B2=[0

0

1.1

1.1

0

1;

0

0

0

0

2;

0

0.343+0.21256i

0

0

2;

0

0

0

0

2;

0

0.204+0.12638i

0

0

2;

0

0

0

0

2;

0

0.306+0.18962i

0

0

2;

0

0

0

0

2;

0.5

0

1.1

1.1

0

3;

0

0.343+0.21256i

0

0

2;

0

0

0

0

2;

0

0

0

0

2];% B1矩陣：

1、支路首端號；

2、末端號；

3、支路阻抗；

4、支路對地電納 %

5、支路的變比；

6、支路首端處于K側為1，1側為0 % B2矩陣：

1、該節點發電機功率；

2、該節點負荷功率；

3、節點電壓初始值 %

4、PV節點電壓V的給定值；

5、節點所接的無功補償設備的容量 %

6、節點分類標號：1為平衡節點（應為1號節點）；2為PQ節點； %

3為PV節點；

%input('請輸入各節點參數形成的矩陣： B2=');Y=zeros(n);e=zeros(1,n);f=zeros(1,n);V=zeros(1,n);sida=zeros(1,n);S1=zeros(nl);% % %--------------------for i=1:nl

%支路數

if B1(i,6)==0

%左節點處于1側

p=B1(i,1);q=B1(i,2);

else

%左節點處于K側

p=B1(i,2);q=B1(i,1);

end

Y(p,q)=Y(p,q)-1./(B1(i,3)*B1(i,5));%非對角元

Y(q,p)=Y(p,q);

%非對角元

Y(q,q)=Y(q,q)+1./(B1(i,3)*B1(i,5)^2)+B1(i,4)./2;%對角元K側

Y(p,p)=Y(p,p)+1./B1(i,3)+B1(i,4)./2;

%對角元1側

end %求導納矩陣

disp('導納矩陣 Y=');disp(Y)%---------------------------G=real(Y);B=imag(Y);

%分解出導納陣的實部和虛部

for i=1:n

%給定各節點初始電壓的實部和虛部

e(i)=real(B2(i,3));

f(i)=imag(B2(i,3));

V(i)=B2(i,4);

%PV節點電壓給定模值

end for i=1:n

%給定各節點注入功率

S(i)=B2(i,1)-B2(i,2);

%i節點注入功率SG-SL

B(i,i)=B(i,i)+B2(i,5);%i節點無功補償量

end %===================== P=real(S);Q=imag(S);

%分解出各節點注入的有功和無功功率 ICT1=0;IT2=1;N0=2*n;N=N0+1;a=0;%迭代次數ICT1、a；不滿足收斂要求的節點數IT2 while IT2~=0

% N0=2*n 雅可比矩陣的階數；N=N0+1擴展列

IT2=0;a=a+1;

for i=1:n

if i~=isb

%非平衡節點

C(i)=0;D(i)=0;

for j1=1: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*fj+Bij*ej)

end

P1=C(i)*e(i)+f(i)*D(i);%節點功率P計算eiΣ(Gij*ej-Bij*fj)+fiΣ(Gij*fj+Bij*ej)

Q1=C(i)*f(i)-e(i)*D(i);%節點功率Q計算fiΣ(Gij*ej-Bij*fj)-eiΣ(Gij*fj+Bij*ej)%求i節點有功和無功功率P',Q'的計算值

V2=e(i)^2+f(i)^2;

%電壓模平方

%========= 以下針對非PV節點來求取功率差及Jacobi矩陣元素 =========

if B2(i,6)~=3

%非PV節點

DP=P(i)-P1;

%節點有功功率差

DQ=Q(i)-Q1;

%節點無功功率差

%=============== 以上為除平衡節點外其它節點的功率計算 ================= %================= 求取Jacobi矩陣 ===================

for j1=1:n

if j1~=isb&j1~=i

%非平衡節點&非對角元

X1=-G(i,j1)*e(i)-B(i,j1)*f(i);% dP/de=-dQ/df

X2=B(i,j1)*e(i)-G(i,j1)*f(i);% dP/df=dQ/de

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;J(p,N)=DQ;m=p+1;

% X3=dQ/de J(p,N)=DQ節點無功功率差

J(m,q)=X1;J(m,N)=DP;q=q+1;

% X1=dP/de J(m,N)=DP節點有功功率差

J(p,q)=X4;J(m,q)=X2;

% X4=dQ/df X2=dp/df

elseif j1==i&j1~=isb %非平衡節點&對角元

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(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;J(p,N)=DQ;%擴展列△Q

m=p+1;

J(m,q)=X1;q=q+1;J(p,q)=X4;J(m,N)=DP;%擴展列△P

J(m,q)=X2;

end

end

else

%=============== 下面是針對PV節點來求取Jacobi矩陣的元素 ===========

DP=P(i)-P1;

% PV節點有功誤差

DV=V(i)^2-V2;

% PV節點電壓誤差

for j1=1:n

if j1~=isb&j1~=i

%非平衡節點&非對角元

X1=-G(i,j1)*e(i)-B(i,j1)*f(i);

% dP/de

X2=B(i,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;J(p,N)=DV;% PV節點電壓誤差

m=p+1;

J(m,q)=X1;J(m,N)=DP;q=q+1;J(p,q)=X6;

% PV節點有功誤差

J(m,q)=X2;

elseif j1==i&j1~=isb %非平衡節點&對角元

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);

p=2*i-1;q=2*j1-1;J(p,q)=X5;J(p,N)=DV;% PV節點電壓誤差

m=p+1;

J(m,q)=X1;J(m,N)=DP;q=q+1;J(p,q)=X6;

% PV節點有功誤差

J(m,q)=X2;

end

end

end

end

end %========= 以上為求雅可比矩陣的各個元素及擴展列的功率差或電壓差 =====================

for k=3:N0

% N0=2*n（從第三行開始，第一、二行是平衡節點）

k1=k+1;N1=N;

% N=N0+1 即 N=2*n+1擴展列△P、△Q 或 △U

for k2=k1:N1

% 從k+1列的Jacobi元素到擴展列的△P、△Q

或 △U

J(k,k2)=J(k,k2)./(J(k,k)+eps);% 用K行K列對角元素去除K行K列后的非對角元素進行規格化

end

J(k,k)=1;

% 對角元規格化K行K列對角元素賦1

%==================== 回代運算

=======================================

if k~=3

% 不是第三行

k > 3

k4=k-1;

for k3=3:k4

% 用k3行從第三行開始到當前行的前一行k4行消去

for k2=k1:N1 %

k3行后各行上三角元素

J(k3,k2)=J(k3,k2)-J(k3,k)*J(k,k2);%消去運算（當前行k列元素消為0）

end

%用當前行K2列元素減去當前行k列元素乘以第k行K2列元素

J(k3,k)=0;%當前行第k列元素已消為0

end

if k==N0

%若已到最后一行

break;

end

%================== 前代運算

==================================

for k3=k1:N0

% 從k+1行到2*n最后一行

for k2=k1:N1

% 從k+1列到擴展列消去k+1行后各行下三角元素

J(k3,k2)=J(k3,k2)-J(k3,k)*J(k,k2);%消去運算

end

%用當前行K2列元素減去當前行k列元素乘以第k行K2列元素

J(k3,k)=0;%當前行第k列元素已消為0

end

else

%是第三行k=3

%====================== 第三行k=3的前代運算 ========================

for k3=k1:N0

%從第四行到2n行（最后一行）

for k2=k1:N1

%從第四列到2n+1列（即擴展列）

J(k3,k2)=J(k3,k2)-J(k3,k)*J(k,k2);%消去運算（當前行3列元素消為0）

end

%用當前行K2列元素減去當前行3列元素乘以第三行K2列元素

J(k3,k)=0;%當前行第3列元素已消為0

end

end

end %====上面是用線性變換方式高斯消去法將Jacobi矩陣化成單位矩陣=====

for k=3:2:N0-1

L=(k+1)./2;

e(L)=e(L)-J(k,N);

%修改節點電壓實部

k1=k+1;

f(L)=f(L)-J(k1,N);

%修改節點電壓虛部

end

%------修改節點電壓-----------

for k=3:N0

DET=abs(J(k,N));

if DET>=pr

%電壓偏差量是否滿足要求

IT2=IT2+1;%不滿足要求的節點數加1

end

end

ICT2(a)=IT2;

%不滿足要求的節點數

ICT1=ICT1+1;

%迭代次數 end %用高斯消去法解”w=-J*V" disp('迭代次數：');disp(ICT1);disp('沒有達到精度要求的個數：');disp(ICT2);for k=1:n

V(k)=sqrt(e(k)^2+f(k)^2);

%計算各節點電壓的模值

sida(k)=atan(f(k)./e(k))*180./pi;

%計算各節點電壓的角度

E(k)=e(k)+f(k)*j;

%將各節點電壓用復數表示 end %=============== 計算各輸出量 =========================== disp('各節點的實際電壓標幺值E為(節點號從小到大排列)：');disp(E);

%顯示各節點的實際電壓標幺值E用復數表示 disp('----------------------');disp('各節點的電壓大小V為(節點號從小到大排列)：');disp(V);

%顯示各節點的電壓大小V的模值 disp('----------------------');disp('各節點的電壓相角sida為(節點號從小到大排列)：');disp(sida);

%顯示各節點的電壓相角 for p=1:n

C(p)=0;

for q=1:n

C(p)=C(p)+conj(Y(p,q))*conj(E(q));%計算各節點的注入電流的共軛值

end

S(p)=E(p)*C(p);

%計算各節點的功率 S = 電壓 X 注入電流的共軛值 end disp('各節點的功率S為(節點號從小到大排列)：');disp(S);

%顯示各節點的注入功率

disp('----------------------');disp('各條支路的首端功率Si為(順序同您輸入B1時一致)：');for i=1:nl

p=B1(i,1);q=B1(i,2);

if B1(i,6)==0

Si(p,q)=E(p)*(conj(E(p))*conj(B1(i,4)./2)+(conj(E(p)*B1(i,5))...-conj(E(q)))*conj(1./(B1(i,3)*B1(i,5))));

Siz(i)=Si(p,q);

else

Si(p,q)=E(p)*(conj(E(p))*conj(B1(i,4)./2)+(conj(E(p)./B1(i,5))...-conj(E(q)))*conj(1./(B1(i,3)*B1(i,5))));

Siz(i)=Si(p,q);

end

disp(Si(p,q));

SSi(p,q)=Si(p,q);

ZF=['S(',num2str(p),',',num2str(q),')=',num2str(SSi(p,q))];

disp(ZF);

disp('----------------------');end disp('各條支路的末端功率Sj為(順序同您輸入B1時一致)：');for i=1:nl

p=B1(i,1);q=B1(i,2);

if B1(i,6)==0

Sj(q,p)=E(q)*(conj(E(q))*conj(B1(i,4)./2)+(conj(E(q)./B1(i,5))...-conj(E(p)))*conj(1./(B1(i,3)*B1(i,5))));

Sjy(i)=Sj(q,p);

else

Sj(q,p)=E(q)*(conj(E(q))*conj(B1(i,4)./2)+(conj(E(q)*B1(i,5))...-conj(E(p)))*conj(1./(B1(i,3)*B1(i,5))));

Sjy(i)=Sj(q,p);

end

disp(Sj(q,p));

SSj(q,p)=Sj(q,p);

ZF=['S(',num2str(q),',',num2str(p),')=',num2str(SSj(q,p))];

disp(ZF);

disp('----------------------');end disp('各條支路的功率損耗DS為(順序同您輸入B1時一致)：');for i=1:nl

p=B1(i,1);q=B1(i,2);

DS(i)=Si(p,q)+Sj(q,p);

disp(DS(i));

DDS(i)=DS(i);

ZF=['DS(',num2str(p),',',num2str(q),')=',num2str(DDS(i))];

disp(ZF);

disp('----------------------');end

如果源程序的運行結果需要作圖可用下面的程序

figure(1);subplot(1,2,1);plot(V);xlabel('節點號');ylabel('電壓標幺值');grid on;subplot(1,2,2);plot(sida);xlabel('節點號');ylabel('電壓角度');grid on;figure(2);subplot(2,2,1);P=real(S);Q=imag(S);bar(P);xlabel('節點號');ylabel('節點注入有功');grid on;subplot(2,2,2);bar(Q);xlabel('節點號');ylabel('節點注入無功');grid on;subplot(2,2,3);P1=real(Siz);Q1=imag(Siz);bar(P1);xlabel('支路號');ylabel('支路首端注入有功');grid on;subplot(2,2,4);bar(Q1);xlabel('支路號');ylabel('支路首端注入無功');grid on;

總結

通過本次課程設計讓我有復習了一次潮流計算的相關知識，跟家清晰了什么事潮流計算以及潮流計算的在電力系統的重要性。電力系統的穩定運行狀況即是正常運行狀況，是指電力系統在穩定運行條件下電壓、功率的分布，也稱為潮流分布。電力系統分析的潮流計算是電力系統分析的一個重要的部分。通過對電力系統潮流分布的分析和計算，可進一步對系統運行的安全性，經濟性進行分析、評估，提出改進措施。同時潮流分布也是電力系統規劃設計的一項基礎工作。

整個計算過程的模型建立并不是十分復雜，但計算過程十分繁瑣、計算量相當的大，而且由于枝節太多很容易算錯。不過在計算潮流計算的過程中卻對以往學過的電力系統分析的相關知識進行了一次較為深入的復習。而且整個計算對計算量的要求很大，鍛煉了我們的計算能力。而且對細節的把握也得到了鍛煉，做題的精細程度得到了提高。

參考文獻

[1]何仰贊, 溫增銀《電力系統分析》（第三版）［M].華中科技大學，2002 [2]http://baike.baidu.com/view/627420.[3]王守相，劉玉田 電力系統潮流計算研究現狀--《山東電力技術》1996年05期

[4]何仰贊，溫增銀.電力系統分析（上冊）第三版［M].湖北：華中科技大學出版社，2002 [5] 劉同娟.MATLAB在電路分析中的應用.電氣電子教學學報.2002 [6] 西安交通大學等.電力系統計算[M].北京：水利電力出版社，1993.12 [7] 李光琦.電力系統暫態分析[M].北京： 水利電力出版社，2002.5 [8]何仰贊，溫增銀.電力系統分析（下冊）第三版［M].湖北：華中科技大學出版社，2002 [9]韋化，李濱，杭乃善，等．大規模水一火電力系統最優潮流的現代內點算法實現[J]．中國電機工程學報，2003．23(6)：13一l8．

[10]Chen Luo—nan，Suzuki Hideki，Katou Kazuo．Mean fieldtheory for optimal power flow[J]．IEEE Transactions OilPower Systems，1997，12(4)：1481·1486