matlab三機(jī)九節(jié)點(diǎn)電力系統(tǒng)仿真（帶程序）

matlab三機(jī)九節(jié)點(diǎn)電力系統(tǒng)仿真（帶程序）

三機(jī)九節(jié)點(diǎn)電力系統(tǒng)

暫態(tài)仿真

學(xué)院：

專業(yè)：

學(xué)號：

姓名：

授課教師:

一、摘要

電力系統(tǒng)仿真計(jì)算己經(jīng)成為電力系統(tǒng)設(shè)計(jì)、運(yùn)行與控制中不可缺少的手段。通過設(shè)置

不同故障類型、不同故障地點(diǎn)運(yùn)用仿真技術(shù)可以對電力系統(tǒng)的暫態(tài)穩(wěn)定進(jìn)行分析。本文采

用IEEE3機(jī)9節(jié)點(diǎn)的經(jīng)典多機(jī)模型,基于隱式梯形積分法對系統(tǒng)發(fā)生三相金屬性短路故障

進(jìn)行仿真，分析系統(tǒng)在這種情況下的暫態(tài)穩(wěn)定。發(fā)電機(jī)模型采用經(jīng)典的二階模型;負(fù)荷采

用恒定阻抗負(fù)荷。在Matlab2010上編寫程序進(jìn)行調(diào)試和運(yùn)行。

電力系統(tǒng)是由不同類型的發(fā)電機(jī)組、多種電力負(fù)荷、不同電壓等級的電力網(wǎng)絡(luò)等組成

的十分龐大復(fù)雜的動(dòng)力學(xué)系統(tǒng)。其暫態(tài)過渡過程不僅包括電磁方面的過渡過程，而且還有

機(jī)電方面的過渡過程。由此可見，電力系統(tǒng)的數(shù)學(xué)模型是一個(gè)強(qiáng)非線性的高維狀態(tài)方程組。

在動(dòng)態(tài)穩(wěn)定仿真中使用簡單的電力系統(tǒng)模型，通過仿真計(jì)算分析說明，此仿真方法可以進(jìn)

行簡單的電力系統(tǒng)暫態(tài)分析，對提高電力系統(tǒng)暫態(tài)穩(wěn)定具有重要意義。

二、案例

本次課程主要應(yīng)用P.M.AndersonandA.A.Fouad編寫的《PowerSystem

ControlandStability》一書中所引用的WesternSystemCoordinatedCouncil(WSCC)

三機(jī)九節(jié)點(diǎn)系統(tǒng)模型。系統(tǒng)電路結(jié)構(gòu)拓?fù)鋱D如下：

0.0625

j0.0586

j0.0576

j18/230

230/13.8

16.5/230

18KV

230KV230KV

230KV

13.8KV16.5KV

11

2

2

3

3

4

5

6

7

8

9

A

負(fù)荷負(fù)荷B

負(fù)荷C0.00850.072

j+/20.0745

Bj=/20.1045

Bj=/20.153

Bj=/20.179

Bj=/20.088

Bj=/20.079

Bj=0.01190.1008j+0.0320.161

j+0.0100.085

j+0.0390.170

j+0.0170.092

J十

圖2-13機(jī)9節(jié)點(diǎn)系統(tǒng)

系統(tǒng)數(shù)據(jù)其中，節(jié)點(diǎn)數(shù)據(jù)如下：

節(jié)點(diǎn)號有無負(fù)載類型電壓相角有功負(fù)荷無功負(fù)荷有功出力無功出力電壓基準(zhǔn)期望電壓

N=[l031.04000.000.000.0071.6027.0016.501.040

2021.02500.000.000.00163.006.7018.001.025

3021.02500.000.000.0085.00-10.9013.801.025

4001.00000.000.000.000.000.00230.001,026

5101.00000.00125.0050.000.000.000.000.996

6101.00000.0090.0030.000.000.000.001.013

7001.00000.000.000.000.000.00230.001.026

8101.00000.00100.0035.000.000.000.001.016

9001.00000.000.000.000.000.00230.001.032];%支路數(shù)據(jù)

%從到電阻電抗容納類型變比

B=[140.00.05760.011

270.00.06250.011

390.00.05860.011

450.0100.0850.17600

460.0170.0920.15800

570.0320.1610.30

600

690.0390.1700.35800

780.00850.0720.14900

890.01190.10080.20900];

發(fā)電機(jī)數(shù)據(jù)如下：

%發(fā)電機(jī)母線XdXd'TdO'XqXq'TqO'TjXf

Ge=[110.14600.06088.960.09690.0969047.280.0576

220.89580.11986.000.86450.19690.53512.800.0625

331.31250.18138.591.25780.25000.6006.020.0585];

三、仿真框圖

在仿真之前，首先，應(yīng)明確仿真的所要到達(dá)的結(jié)果，即仿真目標(biāo):本此仿真的結(jié)果主

要是得到發(fā)電機(jī)攻角、轉(zhuǎn)速隨時(shí)間變化的值，包括故障前、故障中、故障后。故障前，系

統(tǒng)處于穩(wěn)定狀態(tài)，發(fā)電機(jī)的攻角、轉(zhuǎn)速基本穩(wěn)定。而當(dāng)系統(tǒng)發(fā)生故障以及故障切除，系統(tǒng)

結(jié)構(gòu)拓?fù)浒l(fā)生變化，系統(tǒng)的狀態(tài)也將隨時(shí)間發(fā)生變化，為了求取系統(tǒng)狀態(tài)的變化，我們通

過對系統(tǒng)進(jìn)行簡化建立數(shù)學(xué)模型，得到相關(guān)的代數(shù)一微分方程組，進(jìn)行數(shù)值計(jì)算，從而得

到系統(tǒng)狀態(tài)的隨時(shí)間的變化值。此次仿真的系統(tǒng)以發(fā)電機(jī)二階經(jīng)典模型來進(jìn)行系統(tǒng)是數(shù)學(xué)

建模，系統(tǒng)的狀態(tài)量為發(fā)電機(jī)攻角、發(fā)電機(jī)轉(zhuǎn)速。

其次，當(dāng)明確仿真目標(biāo)后，我們就得明確大體的仿真框架流程。

仿真框架流程如下：

數(shù)據(jù)準(zhǔn)備

（支路、節(jié)點(diǎn)、發(fā)電機(jī)）

潮流計(jì)算

計(jì)算故障前中后發(fā)電機(jī)

內(nèi)節(jié)點(diǎn)的導(dǎo)納矩陣

發(fā)電機(jī)初值計(jì)算

列寫系統(tǒng)狀態(tài)方程

（轉(zhuǎn)子運(yùn)動(dòng)方程）

調(diào)用ode45計(jì)算發(fā)電攻角、

轉(zhuǎn)速變化情況

后處理

圖3-1仿真流程圖

四、仿真模型

在電力系統(tǒng)的機(jī)電暫態(tài)仿真中，常根據(jù)實(shí)際要求的不同，采用不同時(shí)間尺度的仿真模

型，而仿真算法和采用的模型有直接的關(guān)系，下面就本次仿真實(shí)例機(jī)電暫態(tài)過程的仿真模

型及其仿真算法。

一、潮流計(jì)算

由于本文以三機(jī)九節(jié)點(diǎn)為模型，假定節(jié)點(diǎn)一為參考節(jié)點(diǎn)，這樣就有2兩個(gè)發(fā)電機(jī)的

PV節(jié)點(diǎn)，6個(gè)PQ節(jié)點(diǎn)，未知量為8個(gè)節(jié)點(diǎn)（包括2個(gè)PV節(jié)點(diǎn)和6個(gè)PQ節(jié)點(diǎn)）的電壓相

角，還有6個(gè)節(jié)點(diǎn)（PQ節(jié)點(diǎn)）的電壓幅值。

可以先求出Y矩陣

圖4-1Y矩陣

然后，我們列寫方程，也就是利用各個(gè)節(jié)點(diǎn)的有功、無功功率的平衡關(guān)系,列寫14

個(gè)功率平衡方程。這樣就能使用牛頓一拉夫遜算法來求解這14個(gè)非線性方程。

其中的關(guān)鍵是要計(jì)算出雅克比矩陣

圖4-2雅克比矩陣

然后計(jì)算出修正量。在設(shè)定精度和最大迭代次數(shù)的前提下進(jìn)行迭代，直到滿足要求。

電力網(wǎng)絡(luò)的節(jié)點(diǎn)功率方程可用一般形式表示如下：

牛頓拉夫遜算法修正方程

W=-JAV

其中w是節(jié)點(diǎn)不平衡量向量，包括有功，無功，電壓；J是雅克比矩陣；戰(zhàn)是節(jié)點(diǎn)

電壓修正量。

令

ij

ijijiiijBGYjfeV?=+=;,

則極坐標(biāo)形式的功率不平衡量方程

)sincos(l=+-=?2=n

jijijijijjiisiBGVVPP88

)cossin(1

jijijijijjiisiBGVVQQ65

雅可比矩陣J各元素的表達(dá)式

=LM

NHJ

當(dāng)jHi時(shí):當(dāng)j=i時(shí):

ijijiiji

j

PHBeGff??=

=-?

iiiiin11

i

PHBeGfbf??=

=—?

IJIJIIJI

J

PNGeBfe??=

——7

i

iiiiiiiii

i

PNGeBfae??=

zz--?

i

ijijiiji

QMBfGef??==+?

QMGeBfaf??=

=+-?

i

UUiUi

J

QLGfBee??=

="+?

i

iiiiiiiii

QLBeGfbe??=

=-+?

其中，

ZZ==+=-=n

ijijjijin

ijijjijieBfGbfBeGal

1)

0

(o

進(jìn)行牛頓拉夫遜算法得到潮流結(jié)果

圖4-3潮流結(jié)果

二、故障前中后僅含發(fā)電機(jī)內(nèi)節(jié)點(diǎn)的導(dǎo)納矩陣

圖4-4故障前中后僅含發(fā)電機(jī)內(nèi)節(jié)點(diǎn)的導(dǎo)納矩陣

三、求解電磁功率

得到故障前，故障中，故障后三個(gè)不同的導(dǎo)納矩陣后，就開始計(jì)算電磁功率和機(jī)械功

率，機(jī)械功率等于穩(wěn)態(tài)的電磁功率中的有功分量。所以可以有

Pe=real(E*I)

如上中，E為發(fā)電機(jī)內(nèi)電勢，I為從發(fā)電機(jī)流出的電流。

但在參考文獻(xiàn)RamnarayanPatel,T.S.BhattiandD.RKothari.MATLAB/Simulink-

basedtransientstabilityanalysisofamultimachinepowersystem中給出的電磁功

率計(jì)算公式為：

WH=+-+=n

i

jjjiijijjiiiieiYEEGEPl

2

)cos(660

穩(wěn)態(tài)情況下有，機(jī)械功率Pme=Pe四、求解運(yùn)動(dòng)方程

發(fā)電機(jī)的運(yùn)動(dòng)方程可以寫成常微分方程組：

-=??????????+-+-=+￡w=RinijjjiijijjiiiimijjiRidt

dYEEGEPDdtdHWOJ6660U）（JOU）1

2）cos（2

其中Pmi為第i個(gè)機(jī)組故障前穩(wěn)態(tài)的電磁功率。在本次仿真中Djcoi為零，即阻尼

為零。仿真開始，t=0時(shí)引入故障，0.083s后切除故障。求解運(yùn)動(dòng)方程后得到曲線如下:

五、結(jié)果分析

上圖分別顯示了各臺發(fā)電機(jī)的轉(zhuǎn)子角與時(shí)間的關(guān)系曲線，顯示了發(fā)電機(jī)轉(zhuǎn)速差的曲線，

和2121588=-.3131655二-的曲線，由圖可以看到，最大角差21S為85o,出現(xiàn)

在0.4ts二處,無論是218還是316第二個(gè)搖擺都不大于第一個(gè)搖擺，可見系統(tǒng)是穩(wěn)

定的。

六、程序代碼

主程序：

globalPmYrungenGenE

%節(jié)點(diǎn)數(shù)據(jù)

%節(jié)點(diǎn)號有無負(fù)載類型電壓相角有功負(fù)荷無功負(fù)荷有功負(fù)荷無功負(fù)荷電壓基準(zhǔn)期望電

N=[l031.04000.000.000.0071.6027.0016.501.040

2021.02500.000.000.00163.006.7018.001.025

3021.02500.000.000.0085.00-10.9013.801.025

4001.00000.000.000.000.000.00230.001.026

5101.00000.00125.0050.000.000.000.000.996

6101.00000.0090.0030.000.000.000.001.013

7001.00000.000.000.000.000.00230.001.026

8101.00000.00100.0035.000.000.000.001.016

9001.00000.000.000.000.000.00230.001.032];

%支路數(shù)據(jù)

%從到電阻電抗容納類型變比

B=[140.00.05760.011

270.00.06250.011

390.00.05860.011

450.0100.0850.17600

460.0170.0920.15800

570.0320.1610.30

600

690.0390.1700.35800

780.00850.0720.14900

890.01190.10080.20900];

%發(fā)電機(jī)數(shù)據(jù)

%HMVAxd'*10000nodexdxqxlxadxaqxftdO'rf

gen=[2364247.560810.14600.09690.03360.11240.06330.14838.96

0.0000439640192.0119820.89580.86450.05210.84370.81240.91736.00

0.0004054301128.018133131251.25780.07421.23831.28361.35555.89

0.0006105];

Y=zeros(9,9);%導(dǎo)納矩陣

fori=l:9

a=B(bl);b=B(i,2);

ifB(i,6)==0

Y(a,b)=-l./(B(i,3)+B(i,4)*li);

Y(b,a)=-l./(B(i/3)+B(if4)*li);

Y(a,a)=Y(a,a)+l./(B(if3)+B(i,4)*j)+B(i,5)*li./2;

Y(b,b)=Y(b,b)+l./(B(⑶+B(i,4)*j)+B(i,5)*li./2;

else

Y(a,b)=-L/((B(i,3)+li*B(i,4))*B(i,6));

Y(b,a)=-l./((B(i,3)+li*B(i(4))*B(i6));

Y(afa)=Y(afa)+l./(B(ir3)+B(L4)*li);

Y(b,b)=Y(b,b)+L/(B(i,3)+B(i,4)*j*B(i,6”2)+B(i,5)*li./2;

end

end%導(dǎo)納矩陣

forT=l:100

[dRdQ]=Caoliu(N,Y);%潮流

J=Ykb(N,Y);%雅克比矩陣

U=zeros(6);

fori=4:9

U(i-3J-3)=N(i,4);

end

dAngU=J\[dP;dQ];

dAng=dAngll(l:8,l);

dU=U*(dAngU(9:14/l));

N(4:9,4)=N(4:9,4)-dU;

N(2:9,5)=N(2:9,5)-dAng;

if(max(abs(dU))<0.00001)&8i(max(abs(dAng))<0.00001)

break

end

end

[Yc/Yb,Ya]=Ynew(gen/N/B/Y);

GEgj=zeros(l,3);

GenE=zeros(l,3)；

fori=l:3

GenE(l,i)=abs(N(i,4)*exp(li*N(i,5))+li*gen(i/3)/10000*conj(((N(i,8)/100+li*N(ir

9)/100)/(N(i,4)*exp(li*N(i/5))))));

GEgj(l,i)=angle(N(i/4)*exp(li*N(i/5))+li*gen(i,3)/10000*conj(((N(i,8)/100+li*N

(i,9)/100)/(N(i,4)*exp(li*N(i,5))))));

end

Yrun=zeros(3);

Pm=zeros(l,3)；

fori=l:3

PmQ,i)=N(i,8)./100;

end

options=odeset('REITOL'/le-10);%設(shè)置精度

X0=[GEgj(l)zGEgj(2),GEgj(3),2*pi*60*ones(l,3)]；

t0=0;

tc=0.083;

tspanl=[tO,tc];

Yrun=Yb;

,

[Tl,Yl]=ode45(fzz,,tspanl/X0,options);

Xl=Yl(end/);

tf=2;

tspan2=[tc,tf];

Yrun=Ya;

[T2,Y2]=ode45('fzz',tspan2,Xl/options);

T=[T1;T2];

Y12=[Y1;Y2];

subplot。,1,1);

plot(T,Y12(:,l)*180/pi,T,Y12(:,2)*180/pi/T,Y12(:,3)*180/pi);%發(fā)電順角

subplot(3,l/2);

plol(T,Y12(:,5)-Y12(:,4),T/Y12(:,6)-Y12(:,4));%發(fā)電機(jī)轉(zhuǎn)速差

subplot(3,l/3)；

plot(T,Y12(:,2)*180/pi-Y12(:/l)*180/pi/T/Y12(:,3)*180/pi-Y12(:/l)*180/pi)%發(fā)電

機(jī)攻角差

雅克比矩陣：

functionJ=Ykb(N,Y)

H=zeros(8);Nl=zeros(8/6);K=zeros(6/8);L=zeros(6);

fori=2:9

forj=2:9

ifi~=j

H(i-l,j-l)=-N(i,4)*N(j,4)*(real(Y(i,j))*sin(N(i,5)-N(j,5))-imag(Y(i,j))*cos(N(i.5)-

N(j,5)));

else

H(i?l”)=(N(i⑼?N(i,7))/100+imag(Y(ij))*((N(i,4))八2);

end

end

end%H

fori=2:9

forj=4:9

ifi~=j

Nl(i-lJ-3)=-N(i,4)*N(j/4)*(real(Y(ij))*cos(N(i/5)-N(j/5))+imag(Y(iJ))*sin(N(if5)-

N(j,5)));

else

Nl(i?Lj-3)=?(N(i,8)?N(i,6))/100?real(Y(i,j))*((N(i,4)r2);

end

end

end%N

fori=4:9

forj=2:9

ifi-=j

K(i?3,j?l)=N(i,4)*N(j,4)*(real(Y(i,j))*cos(N(i,5)?N(j,5))+imag(Y(iJ))*sin(N(i,5)?

N(j,5)));

else

K(i-3J?l)=?(N(i,8)?N(i,6))/100+real(Y(iJ))*((N(i,4))A2);

end

end

end%K

fori=4:9

forj=4:9

ifi~=j

L(i-3J-3)=-N(i/4)*N(j/4)*(real(Y(ij))*sin(N(i/5)-N(j,5))-imag(Y(iJ))*cos(N(i/5)-

N(j,5)));

else

L(i-3,j-3)=-(N(i,9)-N(i,7))/100+imag(Y(i,j))*((N(i,4))入2);

end

end

end%L

J=[HN1;KL];

潮流計(jì)算

function[dBdQ]=Caoliu(N,Y)

dP=zeros(8,l);

dQ=zeros(6,l);

fori=2:9

dP(i-l,l)=(N(i,8)-N(i,6))/100;

end

fori=4:9

dQ(i-3,l)=(N(i,9)-N(i/7))/100;

end

fori=2:9

forj=1:9

dP(i-l,l)=dP(i-lJ)-N(i,4)*N(j/4)*(real(Y(ij))*cos(N(i/5)-N(j/5))+imag(Y(ij))

*sin(N(i/5)-N(j,5)));

end

end

fori=4:9

forj=1:9

dQ(i?3,l)=dQ(i?3,l)-N(i,4)*N(j,4)*(real(Y(i,j))*sin(N(i,5)?N(j,5))-imag(Y(iJ))

*cos(N(i,5)-N(j,5)));

end

end

故障前中后僅含發(fā)電機(jī)內(nèi)節(jié)點(diǎn)的導(dǎo)納矩陣：

function[Yp/Yd,Ya]=Ysim(gen/N/B/Y)

Yp=zeros(3);Yd=zeros(3);Ya=zeros(3);

Yl=zeros(12);Y2=zeros(ll);Y3=zeros(12);

Yl=[zeros(3),zeros(3,9);zeros(9,3),Y];%故障前增廣

AA

Y5=N(5,6)/100/((N(5,4)2))-(N(5/7)/100/((N(5.4)2)))*li;

AA

Y6=N(6,6)/100/((N(6,4)2))-(N(6/7)/100/((N(6,4)2)))*li;

Y8=N(8,6)/100/((N(8,4)A2))-(N(8,7)/100/((N(8,4)A2)))*li;%負(fù)載等效導(dǎo)納

fori=l:3

Yl(i,i)=-(l/gen(i,3)*10000)*li;

Yl(i,i+3)=(l/gen(i,3)*10000)*11；

Yl(i+3J)=Yl(ij+3);

end

fori=l:3

Yl(i+3,i+3)=Yl(i+3,i+3)?(l/gen(i,3)*10000)*li;%發(fā)電機(jī)節(jié)點(diǎn)修改

end

Yl(3+5,3+5)=Yl(3+5,3+5)+Y5;

Yl(3+6,3+6)=Yl(3+6,3+6)+Y6;

Yl(3+8,3+8)=Yl(3+8,3+8)+Y8;%負(fù)載節(jié)點(diǎn)修改

Ynn=zeros(3);Ynr=zeros(3/9);Ynrl=zeros(3/8);Yrn=zeros(9,3)；Yrnl=zeros(8,3)；

Yrr=zeros(9);Yrrl=zeros(8);

Ynn=Yl(l:3,l:3);Ynr=Yl(l:3,4:12);Yrn=Yl(4:12,l:3);Yrr=Yl(4:12,4:12);

Yp=Ynn-Ynr*(YrrA-l)*Yrn;

Y3=Y1;

Y3(8/8)=Y3(8/8)-(l./(B(6,3)+B(6f4)*li)+B(6,5)*li/2);Y3(10,10)=Y3(10/10)-

(l./(B(6/3)+B(6f4)*li)+B(6/5)*li/2);Y3(8,10)=0;Y3(10,8)=0;%故障后增廣

Yl(10,;)=[];

Yl(:,10)=[];

Y2=Y1;%故障中增廣

Ynn=Y2(l:3/l:3);Ynrl=Y2(l:3,4:ll);Yrnl=V2(4:ll/l:3);Yrrl=Y2(4:ll,4:ll);

Yd=Ynn-Ynrl*(YrrlA-l)*Yrnl;

Ynn=Y3(l:3/l:3);Ynr=Y3(l:3,4:12);Ym=Y3(4:12,l:3);Yrr=Y3(4:12,4:12);

Ya=Ynn-Ynr*(YrrA-l)*Yrn;

求解運(yùn)動(dòng)方程：

functionxtt=fouad3(T,DX)

gl