第一篇：基于MATLAB的電力系統潮流計算_張寧

DOI :10.13207

/j.cnki.jnwafu.20

4.12.028

第 32卷 第 12期 402 年 12月

西北農林科技大學學報(自然科學版)

ruoJ.of

No

rthw

est

Sci-Tech

Univ.o

f

Ag

ri.and

F

or.(Nat.Sci.Ed.)loV.32 No.12

Dec.204

基于 MATLAB的電力系統潮流計算

1張 寧,江紅梅,張 渭2

(1西北農林科技大學水利與建筑工程學院,陜西楊凌 71210

;2哈密市高級中學計算機室 ,新疆哈密 83900)

[摘 要] 隨著計算機語言技術的不斷發展和成熟 ,基于 MAT

LA

B的潮流計算研究近年來得到了長足的發 展。針對這一現狀,以 P-Q分解法為例,分析了 BA

SIC , FO RT RA

N和 MA

T

LAB高級語言潮流計算的異同 ,指出了 其優缺點,并針對潮流計算模型結構的特點 ,提出了基于 MA

T LAB的潮流算法。

[關鍵詞] 電力系統;潮流計算;MA

T LAB [中圖分類號] T

M715;T M744

[文獻標識碼] A

[文章編號] 1671-9387(204)12-0124-03

潮流計算是電力系統規劃、運行的基本研究方

法,隨著現代電力系統大系統、強非線性與多元件特 點的日益突出,其計算量與計算復雜度急劇增

[ 1]

加。在處理潮流計算時,其計算機軟件的速度已

[ 2] 無法滿足大電網模擬和實時控制的仿真要求 ,而 高效的潮流問題相關軟件的研究已成為大規模電力

3,4] 系統仿真計算的關鍵[。隨著計算機技術的不斷 發展和成熟,對 MA

T

LA B潮流計算的研究為快速、詳細地解決大電網仿真技術問題開辟了新思 路。針對這一現狀,本文以某電力網絡為例,分 析了 BASIC、FO

RT RA

N和 M A

T

LAB 高級語言潮 流計算的異同,并提出了基于 MA

T LAB 的潮流算 法,以期為電力系統潮流計算提供參考。[ 5,6] 的有效方法。該方法把非線性方程線性化,由于線

性方程的系數矩陣結構上是稀疏的非對稱矩陣,結 合稀疏矩陣技術可使計算機內存占用量大大減少 , 計算速度大大加快;P-Q分解法[

8]是在 N

ew

ton jQ ,(i =1 , 2 ,3 ,…, n)(1)

ji j∑ =P j=1 U＊

i nj=1

計算機潮流計算的基本要求是:計算方法具有

[ 11 , 12 ]

∑·

ji jZ I = 式中,P iQ 分別為節點 i向網絡注入的有功功率，·i和無功功率;jU為節點 j的電壓相量 ·＊ P-Q

i,(i =1 ,2 ,3 ,…, n)(2)

ji·

＊

Ui

一定的可靠性和收斂性;盡量選用占計算機內

存較少的存儲方式;在可靠收斂的前提下,選用計算 速度較快的方法;人機交互環境便于數據輸入、校核 和修改,且具有一定的靈活性。

BAS

IC是一種解決數學問題的語言,由于其取

;Uii

為節點

消了編譯工作及連接過程,語言的執行時間相應較

·

AS

IC語句的標識符可為常數和變量,矩陣也 的電壓共軛相量;jj的電流相量;I為節點 Y為節 慢。Bji點導納矩陣;Z為節點阻抗矩陣。式(1)和式(2)各

ji

[ 13 ]

可規定為變量進行計算。

有 n個非線性復數方程,對其作不同的應用和處 FO RT RAN 是 1957年發明的用于科學計算的

語言理,就形成了不同的潮流計算方法。其中, N ew

ton個工程數學運算函數,可實現潮流計算中的矩陣求.70141E

+38

[ 15 ]。要達到題目所要求的計算精積、求逆、稀疏矩陣形成、復數運算以及初等數學運 算等[ 1]。MA

T

LA B語言允許用戶以數學形式的語 度，且兼顧矩陣程序設計的難易程度 ,MA

T LA

B則成為

言編寫程序 ,其比 BASIC 語言和 FO RT RA

N等更

首選潮流計算的計算機語言。對于大矩陣的潮流計 為接近書寫的數學表達格式 ,且程序易調試。在計

算 ,M A

T LAB提供的 M文件可將輸入矩陣按格式

算要求相同的情況下 ,使用 MA

T LAB編程 ,工作量

先寫入一個文本文件,在編程時 ,按文本文件名在命

將會大為減少。

令窗調用 ,大矩陣就被輸入到內存中了。MA

T LAB 3

基于 MA

T LAB的電力系統潮流計 也能直接創建復數矩陣,這兩點也是 BA

SIC

,(2)求系數矩陣 B

′, B ″的逆陣,計算各節點電壓

FO RT

RA

N 算 的相位角 語言不可比擬的。Δδ。相位角 iΔiδ計算公式為

:

P-Q分解法計算電力系統的潮流分布 ,其

試用

-(k)

i步驟為:

-(B ′)(ΔP

(k)/ U

(k))=(U

式中 , ΔP

(k), U

(k)為已知量,相位角 Δδ的計算主(1)形成系數矩陣 B ′, B ″(設該矩陣為

4×4矩

i要(k)Δδ)是矩陣的求積運陣)。

算

。AMA

T L

B和 BA

SIC語言都用命令 inv

(B)或 a

1a

2…

a 14 IN

V

(B)來實現矩陣的求逆功能,B

′=B

″= a

a

…

MA

T LA

B和 … … … …

a 24 BAS

IC語言求積運算均與數學書寫格式一致 ,但 a

a

…

BAS

IC語言必須對矩陣進行 MA

T說明 ,而且執行

a 44

在程序設計上 ,用 MA

T LA

B語言編寫程序形

時間相應較長 ,計算精度較低。FO RT

RAN 語言雖

成的系數矩陣 B

′和B

″,遠比 BASIC

,FO

RT RA

N簡

然在參考資料上提供矩陣求逆、求積程序[ 1618

] ,但

單,矩陣輸入、輸出與數學書寫格式相似。即在命令

其按線性代數的矩陣求逆、求積步驟編程[ ] ,編寫窗口輸入

: 的程序至少需用 1個三重循環語句 , 所以與

>>B

… a

(3)相比計算平衡節點功率和線路功率MA

T LA

B ,工作量大而效率低。

=[

a

1a

… a 2 2a

a 1平衡節點功率為 …

a

4]

再在命令窗口輸入 >>B

=,則窗口將顯示出 B

′, B ″矩陣。

… a

… … … a

4

2￣S

s

=

+Q

線路功率為

·

＊·

U s is

i =P 1 ＊

Y

U

i =1

n

∑

￣S

ij

= U

iI

=P

ji

ji

+Q

ij

平衡節點功率和線路功率的計算屬復數運算。MA

T LA

B和 FO RT

RAN 語言都提供復數功能語 在 BASIC

, FO RT RAN 語言編程時 , 必須進行

句[ 15 ],MA

T LAB 以數學上的復數書寫格式編寫程 矩陣說明 , BASIC 用 DIM B(4, 4)語句說明矩陣，序 ,方便而且不容易出錯;而 FO RT

RA

N必須對復

F ORT

RA

N用 DOU

BLE

PRECISION

B(1 ∶4

, 1 ∶4

數變量進行)

CMP LX語句的復數類型說明[

],顯

語句同時說明矩陣和矩陣精度;在建立矩陣時, 得較為繁瑣。BA

SIC不提供復數功能語句 ,一般將實

BASIC

,F

ORT RA

N語言用 REA

D語句和 DA

TA 語 數、虛數分開計算 ,最后用輸出語句寫在一起 ,所以 句輸入矩陣;FO RT

RAN 語言用 WRIT E語句輸出 一般不用 BASIC 語言編寫復數運算程序

。矩陣,而 BASIC 語言只能用二重循環語句輸出二維

另外 ,MA

T LA

B可以提供潮流計算中稀疏矩陣

· ·＊

126

西北農林科技大學學報(自然科學版)

第 32卷 的建立命令 spconvetr,從而將外部數據轉化為稀疏 矩陣,而且 MA

T

LA B函數 lu也可直接實現 L

R三 角分解[了

方便。

一步。矩陣輸入、輸出格式簡單 ,與數學書寫格式相

似;以雙精度類型進行數據的存儲和運算 ,數據精確

;潮流計算中復雜矩陣的輸入問題可通過創建 M文件來解決;MA

T LA B稱為矩陣實驗室 ,其能進 行潮流計算中的各種矩陣運算,包括求逆、求積和矩 19 , 20

],從而為現代電力系統潮流計算提供度高 討

論

陣 L

R分解等 ,其程序的編寫也因 M A

T LAB提供 了許多功能函數而變得簡單易行。另外 , MA

T LAB

稀疏矩陣技術的引入 ,使電力系統潮流計算由傳統 方法轉變為優化算法成為可能。通過與 BA

SIC, FO RT RA

N語言的比較,基于 M A

T LAB的電力系統潮流計算使計算機在計算、分 析、研究復雜的電力系統潮流分布問題上又前進了

[參考文獻] [ 1 ]

薛

巍 ,舒繼武 ,王心豐 ,等.電力系統并行算法的研究進展 [ J ]

.清華大學學報(自然科學版), 2002

, 42

(9):1192

-195

.[ 2 ]

F

ernaod

L A

.Computalnoi

complexit

y in

pow

er

systems

[ J ]

.IE

E

T

rans

on

PAS, 1976 , 95(4):1028

-1037

.[ 3 ]

周濟 ,羅應立 ,張建華 ,等.基于 MAT

LA B的次同步諧振的特征根分析 [ J ]

.現代電力 , 1999 ,(4):67-71.[ 4 ]

王克英 ,穆

鋼 ,韓學山 ,等.使潮流方程直接可解的 PMV配置方案研究 [ J ]

.中國電機工程學報 , 1999 , 10

(2):67-69.[ 5]

盧天成 ,韓

或 ,李東風.MA

T LAB在鳴禽語圖分析中的應用 [

J ].遼寧師范大學學報(自然科學版), 2001

, 24

(4):56-58.[ 6]

宋維君.基于 MA

T LAB的時域分析法實現 [ J ]

.遼寧師專學報(自然科學版), 2002, 1(4):72

-74.[ 9]

邱關源.電路(下冊)[ M]

.北京 :高等教育出版[ ,7] 0 周孝信等..國電力百科全書 ·電力系統卷 [ M].第 2版.北京 :中國電力出版社 , 2001

.31

.社 199

.67-78 [ 10

]

趙晉泉 ,侯志儉 ,吳際舜.牛頓最優潮流算法中離散控制量的新處理方法 [

J ].電力系統自動化 , 1999 , 23

(3):32

[ 8]

陳

珩.電力系統穩態分析 [ M]

.北京 :高等教育出版社 ,20

.139

-193

.-34

.[ 11

]

[英]勞頓 M A,塞 M G.電氣工程師技術手冊 [

M]

.北京 :機械工業出版社 , 1992

.826

-827

.[ 12

]

樓順天 ,陳生潭 ,雷虎明.MA

T LAB5

.X程序設計語言 [ M]

.西安 :西安電子科技大學出版社 , 2001

.57-78;163

-162

.[ 14

]

潭浩強 ,田淑清.F

ORT

RA

N語言 ———F ORT

RA

N 77結構化程序設計 [

M].北京 :清華大學出版[ 13

]

王沫然.MA

T LAB

6.0與科學計算[

M]

.北京 :電子工業出版社 , 2001

.1

-5;42

.社 , 1995.58;98.[ 15 ]

劉文遠 ,李興成 ,趙

莉.BA

SIC語言與農業實用程序 [ 041

.[ 16 ]

潭浩強 ,崔武子 ,田淑清.FO RT

RA

N程序設計(二級)教程[

社 , 1998.105

-123

.[ 17 ]

鮑有文 ,周海燕 ,趙重敏 ,等.FO RT

RAN

77程序設計試題匯編 [ M]

.北京 :清華大學出版[ 20

]

鄧建中 ,葛仁杰 ,程正興.計算方法[ M]

.西安 :西安交通大學出版社 , 1998

.21

社 ,.68;107

.-31

1998;201

.[ 18 ]

劉萬春 ,朱玉文 ,龔圓明.FO RT

RA

N程序設計[

M].北京 :國防工業出版社 , 2003

.102

-124

.[ 19 ]

同濟大學教研室.線性代數[

M].北京 :高等教育出版社 , 1990

.35

;70-71

.M].北京 :清華大學出版

M]

.沈陽 :遼寧科學技術出版社 , 1987.139

-Pow

er

flow

computaino

of

the

electir

pow

er

system

based

on

MA

TLA

B

1GNAHZ

Ning, JIANG

Hong

-mei

,ZHANG

Wei

(1

Colge

of

Water

Resocru

and

A

rchlarueti

E

ngiren

, Northwest

A

&F University, Yanglin, Shanxi

71210

, China

;

Computer

Rom

of

High

Schol

in

Hami, Ham i, Xinj

iang

83900, China)

Acartsb

:Wit

h

the

devolpment

of

the

computer

langseu

in

rectn

years

,t

he

reshcra

of

M A

T LAB

power wolf

algohmtirs

got

substanli

devolpment

.T

he

paper

, using

P

-Q

method

, analy

zed

the

dif

erensc

andimsilarseti

of

BA

SIC

, FO RT RAN

and

MA

T LAB

advaecn

langseu, and

pointed

out

t hei

r

advatn

ages

and M A

T LAB

pow

er

flow

algoir

thm

s.yeK

words

:electir

pow

er

sy

stem

;power

f low

computaino

;MA

T LAB w

eaknes

, based

on

the

charsietc

of

the

model

const

ructnoi

in

pow

er

flow

algohmtirs

, bring

up

t he

第二篇：用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

第三篇：基于MATLAB的電力系統潮流計算

基于MATLAB的電力系統潮流計算

%簡單潮流計算的小程序，相關的原始數據數據數據輸入格式如下：

%B1是支路參數矩陣，第一列和第二列是節點編號。節點編號由小到大編寫 %對于含有變壓器的支路，第一列為低壓側節點編號，第二列為高壓側節點 %編號，將變壓器的串聯阻抗置于低壓側處理。%第三列為支路的串列阻抗參數。%第四列為支路的對地導納參數。

%第五烈為含變壓器支路的變壓器的變比

%第六列為變壓器是否是否含有變壓器的參數，其中“1”為含有變壓器，%“0”為不含有變壓器。

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

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

n=input('請輸入節點數:n=');n1=input('請輸入支路數:n1=');isb=input('請輸入平衡節點號:isb=');pr=input('請輸入誤差精度:pr=');B1=input('請輸入支路參數:B1=');B2=input('請輸入節點參數:B2=');X=input('節點號和對地參數:X=');Y=zeros(n);Times=1;%置迭代次數為初始值 %創建節點導納矩陣 for i=1:n1 if B1(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));end end Y OrgS=zeros(2*n-2,1);DetaS=zeros(2*n-2,1);%將OrgS、DetaS初始化

%創建OrgS，用于存儲初始功率參數 h=0;j=0;for i=1:n %對PQ節點的處理 if i~=isb&B2(i,6)==2 h=h+1;for j=1:n OrgS(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)));end end end for i=1:n %對PV節點的處理，注意這時不可再將h初始化為0 if i~=isb&B2(i,6)==3 h=h+1;for j=1:n OrgS(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)));end end end OrgS %創建PVU 用于存儲PV節點的初始電壓 PVU=zeros(n-h-1,1);t=0;for i=1:n if B2(i,6)==3 t=t+1;PVU(t,1)=B2(i,3);end end PVU %創建DetaS，用于存儲有功功率、無功功率和電壓幅值的不平衡量 h=0;for i=1:n %對PQ節點的處理 if i~=isb&B2(i,6)==2 h=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);end end t=0;for i=1:n %對PV節點的處理，注意這時不可再將h初始化為0 if i~=isb&B2(i,6)==3 h=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;end end DetaS %創建I，用于存儲節點電流參數 i=zeros(n-1,1);h=0;for i=1:n if i~=isb h=h+1;I(h,1)=(OrgS(2*h-1,1)-OrgS(2*h,1)*sqrt(-1))/conj(B2(i,3));end end I %創建Jacbi(雅可比矩陣)Jacbi=zeros(2*n-2);h=0;k=0;for i=1:n %對PQ節點的處理 if B2(i,6)==2 h=h+1;for j=1:n if j~=isb k=k+1;if i==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);end if k==(n-1)%將用于內循環的指針置于初始值，以確保雅可比矩陣換行

k=0;end end end end end k=0;for i=1:n %對PV節點的處理 if B2(i,6)==3 h=h+1;for j=1:n if j~=isb k=k+1;if i==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;end if k==(n-1)%將用于內循環的指針置于初始值，以確保雅可比矩陣換行

k=0;end end end end end Jacbi %求解修正方程，獲取節點電壓的不平衡量 DetaU=zeros(2*n-2,1);DetaU=inv(Jacbi)*DetaS;DetaU %修正節點電壓 j=0;for i=1:n %對PQ節點處理 if B2(i,6)==2 j=j+1;B2(i,3)=B2(i,3)+DetaU(2*j,1)+DetaU(2*j-1,1)*sqrt(-1);end end for i=1:n %對PV節點的處理 if B2(i,6)==3 j=j+1;B2(i,3)=B2(i,3)+DetaU(2*j,1)+DetaU(2*j-1,1)*sqrt(-1);end end B2 %開始循環********************************************************************** while abs(max(DetaU))>pr OrgS=zeros(2*n-2,1);%！!初始功率參數在迭代過程中是不累加的，所以在這里必須將其初始化為零矩陣 h=0;j=0;for i=1:n if i~=isb&B2(i,6)==2 h=h+1;for j=1:n OrgS(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)));end end end for i=1:n if i~=isb&B2(i,6)==3 h=h+1;for j=1:n OrgS(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)));end end end OrgS %創建DetaS h=0;for i=1:n if i~=isb&B2(i,6)==2 h=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);end end t=0;for i=1:n if i~=isb&B2(i,6)==3 h=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;end end DetaS %創建I i=zeros(n-1,1);h=0;for i=1:n if i~=isb h=h+1;I(h,1)=(OrgS(2*h-1,1)-OrgS(2*h,1)*sqrt(-1))/conj(B2(i,3));end end I %創建Jacbi Jacbi=zeros(2*n-2);h=0;k=0;for i=1:n if B2(i,6)==2 h=h+1;for j=1:n if j~=isb k=k+1;if i==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);end if k==(n-1)k=0;end end end end end k=0;for i=1:n if B2(i,6)==3 h=h+1;for j=1:n if j~=isb k=k+1;if i==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;end if k==(n-1)k=0;end end end end end Jacbi DetaU=zeros(2*n-2,1);DetaU=inv(Jacbi)*DetaS;DetaU %修正節點電壓 j=0;for i=1:n if B2(i,6)==2 j=j+1;B2(i,3)=B2(i,3)+DetaU(2*j,1)+DetaU(2*j-1,1)*sqrt(-1);end end for i=1:n if B2(i,6)==3 j=j+1;B2(i,3)=B2(i,3)+DetaU(2*j,1)+DetaU(2*j-1,1)*sqrt(-1);end end B2 Times=Times+1;%迭代次數加1 end Times 一個原始數據的例子 節點數 5 支路數 5平衡節點編號 5 精度pr 0.000001 B1(支路參數矩陣)[1 2 0.04+0.25i 0.5i 1 0;1 3 0.1+0.35i 0 1 0;2 3 0.08+0.30i 0.5i 1 0;4 2 0.015i 0 1.05 1;5 3 0.03i 0 1.05 1] B2(節點參數矩陣)[0-1.6-0.8i 1 0 0 2;0-2-1i 1 0 0 2;0-3.7-1.3i 1 0 0 2;0 5+0i 1.05 1.05 0 3;0 0 1.05 1.05 0 1] X(節點號和對地參數)[1 0;2 0;3 0;4 0;5 0]

電力系統潮流計算

——9結點算例-PQ法

原始數據錄入data.txt文檔:

標號,起始結點,終止結點,支路電阻參數,支路電抗參數,支路對地導納參數 1,2,5,0.0,0.063,0.0, 2,5,9,0.019,0.072,0.075, 3,6,9,0.012,0.101,0.105, 4,3,6,0.0,0.059,0.0, 5,6,8,0.039,0.17,0.179, 6,4,8,0.017,0.092,0.079, 7,5,7,0.032,0.161,0.153, 8,4,7,0.01,0.085,0.088, 9,1,4,0.0,0.058,0.0, 潮流程序chaoliu.txt文檔: #include #include #define N 9 /*總結點數*/ #define M 6 /*PQ結點數*/ #define K 9 /*線路數*/ #define eps 1e-4 void guass(int n,int m,float c[],float b[][N],float x[])/*高斯函數*/ { float a[N][N],y[N];int i,j,k;for(i=0;i=0;i--){ for(j=i+1;jHeadnode-1;j=t->Endnode-1;lr=t->R;lx=t->X;lb1=t->b;lg=lr/(lr*lr+lx*lx);lb=-lx/(lr*lr+lx*lx);g[i][i]+=lg;g[j][j]+=lg;b[i][i]+=lb+lb1;b[j][j]+=lb+lb1;h[i][j]=h[j][i]=-lb1;g[i][j]-=lg;g[j][i]-=lg;b[i][j]-=lb;b[j][i]-=lb;} getch();printf(“n=====jie dian dao na ju zhen=====n”);for(i=0;i

;

float x1[N-1],x2

;for(i=1;imax)max=fabs(dpu[i]);} if(max>=eps){ for(i=0;i

guass(1,N-1,y1,B,x1);for(i=1;imax)max=fabs(dqu[i]);} if(max>=eps){ for(i=0;i

guass(N-M,M,y2,B,x2);for(i=N-M;i

else { kp=0;if(kq==0)val(u,g,b,r,ku,kr,h);else goto top;} } } void val(float u[N],float g[N][N],float b[N][N],float r[N],int ku, int kr,float h[N][N]){ float ps=0,pv1=0,pv2=0;float qs=0,qv1=0,qv2=0;float p[N][N]={0};float q[N][N]={0};float s[N][N];float dp[N][N]={0};float dq[N][N]={0};float ds[N][N];float dSp=0,dSq=0;int i,j;FILE *fp1;printf(”n=====ping heng jie dian gong lv =====n“);getch();for(i=0;i

printf(”n=======shu ju bao cun=====n“);fp1=fopen(”jieguo.txt“,”w+“);{ fprintf(fp1,”xian lu cao liu:n“);for(i=0;i

第四篇：基于MATLAB的電力系統潮流計算設計

關鍵詞：電力系統分析；潮流計算；matlab仿真

中圖分類號：tm744 文獻標識碼：a 文章編號：1006-4311（2016）21-0185-03

0 引言

潮流計算是電力系統穩態運行中的基本計算方法中的一種計算方法，也是電力系統穩態運行中最重要的運算。潮流計算是保證電力系統安全、經濟運行的根本。在新電網建設的初期規劃中，有了潮流計算，可規劃出電源的容量及其接入點，可計算出無功補償的容量，選擇合適的補償方式，以滿足在電網潮流的控制、調壓、調相、調峰的交換要求。潮流計算可以選擇電力系統的運行方式，便于定期對電力系統中的元件進行檢修。

潮流計算的過程

1.1 原始資料

①系統圖：兩個發電廠分別通過變壓器和輸電線路與四個變電所相連。（圖1）

②發電廠資料：

③變電所資料：

1）變電所1、2、3、4低壓母線的電壓等級分別為：10kv，35kv，10kv，35kv。

3）每個變電所的功率因數均為cosφ=0.9。

④輸電線路資料：

發電廠和變電所之間的輸電線路的電壓等級及長度標于圖中，單位長度的電阻為0.17ω，單位長度的電抗為0.402ω，單位長度的電納為2.78*10-6s。

1.2 基本要求

①對給定的網絡查找潮流計算所需的各元件等值參數，畫出等值電路圖。

②輸入各支路數據，各節點數據，利用simulink搭建仿真模型等方法，進行在變電所的某一負荷情況下的潮流計算及仿真，并對計算結果進行分析。

③如果各母線電壓不滿足要求，進行電壓的調整。（變電所低壓母線電壓10kv要求調整范圍在9.5-10.5之間；電壓35kv要求調整范圍在35-36之間）。

④利用matlab軟件，進行上述各種情況潮流的計算及仿真。

1.3 節點設置及分析

由上述系統圖可知，該系統圖為雙端供電網絡。將母線1，2設為節點1，10，將變電所1、2、3、4的高低壓側分別設為節點2、3、4、5、6、7、8、9。并且，將節點1設為平衡節點，將節點10設為pv節點，其余節點設為pq節點。

1.4 參數求取

將參數整理如表

1、表2所示。

1.5 進行潮流計算

圖2為仿真模型圖。

從潮流計算的結果可得到，系統的各個節點電壓的標幺值可歸納為表3。

由matlab編程調節后，可得到表4的發電廠電壓和變壓器分接頭電壓得標幺值。

在得到了上述調節后的電壓標幺值，對電機模型和變壓器模型進行更改。表5為調節前后各節點的電壓標幺值。

由題意可知，變電所低壓母線電壓10kv要求調整范圍在9.5-10.5之間；電壓35kv要求調整范圍在35-36之間。因此我們可以看出，經過調節后，節點3、5、7、9點電壓已經滿足了系統的要求。表6是電壓調節前后對線路損耗進行分析的記錄。

由表6的電壓調節前后功率損耗對比，可以看出有功功率隨著變壓器分接頭變比的增大而逐漸增大，使得變壓器的低壓側的電壓處于允許范圍內，符合其要求。

表7為調節后的各支路電壓首末端的功率整理表

表8為各節點功率s的標幺值。

1.6 對比

由上面的三種方法簡單地比較，我們可以看出，在同一個電力系統中，用不同的方法進行潮流計算，所得到的結果是大致相同的。

結束語

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

南 京 理 工 大 學

《電力系統穩態分析》

課程報告

姓名

XX

學 號： 5*** 自動化學院 電氣工程

基于牛頓-拉夫遜法的潮流計算例題編程報學院(系): 專

業: 題

目: 任課教師 碩士導師 告

楊偉 XX

2015年6月10號

基于牛頓-拉夫遜法的潮流計算例題編程報告

摘要：電力系統潮流計算的目的在于：確定電力系統的運行方式、檢查系統中各元件是否過壓或者過載、為電力系統繼電保護的整定提供依據、為電力系統的穩定計算提供初值、為電力系統規劃和經濟運行提供分析的基礎。潮流計算的計算機算法包含高斯—賽德爾迭代法、牛頓-拉夫遜法和P—Q分解法等，其中牛拉法計算原理較簡單、計算過程也不復雜，而且由于人們引入泰勒級數和非線性代數方程等在算法里從而進一步提高了算法的收斂性和計算速度。同時基于MATLAB的計算機算法以雙精度類型進行數據的存儲和運算, 數據精確度高，能進行潮流計算中的各種矩陣運算，使得傳統潮流計算方法更加優化。

一 研究內容

通過一道例題來認真分析牛頓-拉夫遜法的原理和方法（采用極坐標形式的牛拉法），同時掌握潮流計算計算機算法的相關知識，能看懂并初步使用MATLAB軟件進行編程，培養自己電力系統潮流計算機算法編程能力。

例題如下：用牛頓-拉夫遜法計算下圖所示系統的潮流分布，其中系統中5為平衡節點，節點5電壓保持U=1.05為定值，其他四個節點分別為PQ節點，給定的注入功率如圖所示。計算精度要求各節點電壓修正量不大于10-6。

二 牛頓-拉夫遜法潮流計算 1 基本原理

牛頓法是取近似解x(k)之后，在這個基礎上，找到比x(k)更接近的方程的根，一步步地迭代，找到盡可能接近方程根的近似根。牛頓迭代法其最大優點是在方程f(x)=0的單根附近時誤差將呈平方減少，而且該法還可以用來求方程的重根、復根。電力系統潮流計算，一般來說，各個母線所供負荷的功率是已知的，各個節點的電壓是未知的（平衡節點外）可以根據網絡結構形成節點導納矩陣，然后由節點導納矩陣列寫功率方程，由于功率方程里功率是已知的，電壓的幅值和相角是未知的，這樣潮流計算的問題就轉化為求解非線性方程組的問題了。為了便于用迭代法解方程組，需要將上述功率方程改寫成功率平衡方程，并對功率平衡方程求偏導，得出對應的雅可比矩陣，給未知節點賦電壓初值，將初值帶入功率平衡方程，得到功率不平衡量，這樣由功率不平衡量、雅可比矩陣、節點電壓不平衡量（未知的）構成了誤差方程，解誤差方程，得到節點電壓不平衡量，節點電壓加上節點電壓不平衡量構成節點電壓新的初值，將新的初值帶入原來的功率平衡方程，并重新形成雅可比矩陣，然后計算新的電壓不平衡量，這樣不斷迭代，不斷修正，一般迭代三到五次就能收斂。2 基本步驟和設計流程圖

形成了雅克比矩陣并建立了修正方程式，運用牛頓-拉夫遜法計算潮流的核心問題已經解決，已有可能列出基本計算步驟并編制流程圖。由課本總結基本步驟如下：

1)形成節點導納矩陣Y；

2)設各節點電壓的初值，如果是直角坐標的話設電壓的實部e和虛部f；如果是極坐標的話則設電壓的幅值U和相角a；

3)將各個節點電壓的初值代入公式求修正方程中的不平衡量以及修正方程的系數矩陣的雅克比矩陣；

4)解修正方程式，求各節點電壓的變化量，即修正量； 5)計算各個節點電壓的新值，即修正后的值；

6)利用新值從第（3）步開始進入下一次迭代，直至達到精度退出循環； 7)計算平衡節點的功率和線路功率，輸出最后計算結果； ① 公式推導

② 流程圖

三

matlab編程代碼

clear;

% 如圖所示1，2，3，4為PQ節點,5為平衡節點

y=0;

% 輸入原始數據，求節點導納矩陣

y(1,2)=1/(0.07+0.21j);

y(4,5)=0;y(1,3)=1/(0.06+0.18j);

y(1,4)=1/(0.05+0.10j);

y(1,5)=1/(0.04+0.12j);

y(2,3)=1/(0.05+0.10j);

y(2,5)=1/(0.08+0.24j);

y(3,4)=1/(0.06+0.18j);

for i=1:5

for j=i:5

y(j,i)=y(i,j);

end

end

Y=0;

% 求節點導納矩陣中互導納

for i=1:5

for j=1:5

if i~=j

Y(i,j)=-y(i,j);

end

end

end

% 求節點導納矩陣中自導納

for i=1:5

Y(i,i)=sum(y(i,:));

end

Y

% Y為導納矩陣

G=real(Y);

B=imag(Y);% 輸入原始節點的給定注入功率

S(1)=0.3+0.3j;

S(2)=-0.5-0.15j;

S(3)=-0.6-0.25j;

S(4)=-0.7-0.2j;

S(5)=0;

P=real(S);

Q=imag(S);

% 賦初值，U為節點電壓的幅值，a為節點電壓的相位角

U=ones(1,5);

U(5)=1.05;

a=zeros(1,5);

x1=ones(8,1);

x2=ones(8,1);

k=0;

while max(x2)>1e-6

for i=1:4

for j=1:4

H(i,j)=0;

N(i,j)=0;

M(i,j)=0;

L(i,j)=0;

oP(i)=0;

oQ(i)=0;

end

end

% 求有功、無功功率不平衡量

for i=1:4

for j=1:5

oP(i)=oP(i)-U(i)*U(j)*(G(i,j)*cos(a(i)-a(j))+B(i,j)*sin(a(i)-a(j)));

oQ(i)=oQ(i)-U(i)*U(j)*(G(i,j)*sin(a(i)-a(j))-B(i,j)*cos(a(i)-a(j)));

end

oP(i)=oP(i)+P(i);

oQ(i)=oQ(i)+Q(i);

end

x2=[oP,oQ]';

% x2為不平衡量列向量

% 求雅克比矩陣

% 當i~=j時，求H,N,M,L

for i=1:4

for j=1:4

if i~=j

H(i,j)=-U(i)*U(j)*(G(i,j)*sin(a(i)-a(j))-B(i,j)*cos(a(i)-a(j)));

N(i,j)=-U(i)*U(j)*(G(i,j)*cos(a(i)-a(j))+B(i,j)*sin(a(i)-a(j)));

L(i,j)=H(i,j);

M(i,j)=-N(i,j);

end

end

end

% 當i=j時，求H,N,M,L

for i=1:4

for j=1:5

if i~=j H(i,i)=H(i,i)+U(i)*U(j)*(G(i,j)*sin(a(i)-a(j))-B(i,j)*cos(a(i)-a(j)));N(i,i)=N(i,i)-U(i)*U(j)*(G(i,j)*cos(a(i)-a(j))+B(i,j)*sin(a(i)-a(j)));

M(i,i)=M(i,i)-U(i)*U(j)*(G(i,j)*cos(a(i)-a(j))+B(i,j)*sin(a(i)-a(j)));

L(i,i)=L(i,i)-U(i)*U(j)*(G(i,j)*sin(a(i)-a(j))-B(i,j)*cos(a(i)-a(j)))

end

end

N(i,i)=N(i,i)-2*(U(i))^2*G(i,i);

L(i,i)=L(i,i)+2*(U(i))^2*B(i,i);

end

J=[H,N;M,L]

% J為雅克比矩陣

x1=-((inv(J))*x2);

% x1為所求△x的列向量

% 求節點電壓新值，準備下一次迭代

for i=1:4

oa(i)=x1(i);

oU(i)=x1(i+4)*U(i);

end

for i=1:4

a(i)=a(i)+oa(i);

U(i)=U(i)+oU(i);

end

k=k+1;

end

k,U,a

% 求節點注入功率

i=5;

for j=1:5

P(i)=U(i)*U(j)*(G(i,j)*cos(a(i)-a(j))+B(i,j)*sin(a(i)-a(j)))+P(i);

Q(i)=U(i)*U(j)*(G(i,j)*sin(a(i)-a(j))-B(i,j)*cos(a(i)-a(j)))+Q(i);

end

S(5)=P(5)+Q(5)*sqrt(-1);

S

% 求節點注入電流

I=Y*U'

四

運行結果

節點導納矩陣

經過五次迭代后的雅克比矩陣

迭代次數以及節點電壓的幅值和相角（弧度數）

節點注入功率和電流

五 結果分析

在這次學習和實際操作過程里：首先，對電力系統分析中潮流計算的部分特別是潮流計算的計算機算法中的牛頓-拉夫遜法進行深入的研讀，弄明白了其原理、計算過程、公式推導以及設計流程。牛頓-拉夫遜法是求解非線性方程的迭代過程，其計算公式為?F?J?X，式中J為所求函數的雅可比矩陣；?X為需要求的修正值；?F為不平衡的列向量。利用x(*)=x(k+1)+?X(k+1)進行多次迭代，通過迭代判據得到所需要的精度值即準確值x(*)。六 結論

通過這個任務，自己在matlab編程，潮流計算，word文檔的編輯功能等方面均有提高，但也暴漏出一些問題：理論知識儲備不足，對matlab的性能和特點還不能有一個全面的把握，對word軟件也不是很熟練，相信通過以后的學習能彌補這些不足，達到一個新的層次。