薛程雙案

課程名稱：電力系統分析實驗

任課教師：_________________________

所屬院部：_____________

教學班級：____________

教學時間：_________

p1

實驗一、節點導納矩陣的形成與修改(設計性)

一、本次課主要內容

(1)MATLAB環境;

(2)MATLAB常用命令;

(3)節點導納矩陣的形成；

(4)節點導納矩陣的修改。

二、教學目的與要求

(1)熟悉MATLAB開發環境；

(2)掌握節點導納矩陣的形成方法和編程；

(3)掌握節點導納矩陣的修改方法。

三、教學重點難點

(1)MATLAB環境;

(2)節點導納矩陣的形成方法和編程。

四、教學方法和手段

課堂講授、提問、討論、演示、實際操作等。

五、作業與習題布置

撰寫實驗報告

P2

實驗一、節點導納矩陣的形成與修改（設計性）

實驗時間：_______________________________

一、實驗目的

（1）熟悉MATLAB開發環境;

（2）掌握節點導納矩陣的形成方法和編程；

（3）掌握節點導納矩陣的修改方法。

二、實驗內容或原理

（1）MATLAB開發環境和MATLAB常用命令的復習；

（2）形成節點導納矩陣的原始數據；

1）節點數n；2）支路數nl；

3）支路參數矩陣B

它包括六個數據［i,j為支路兩端節點號,z為支路的阻抗,

b為線路電納,t為變比，it為高低壓側標志（高為1,低為0）。

4）節點對地阻抗矩陣X（由節點號與接地阻抗構成）。

（3）編寫程序流程圖；

P4

2、教材P122例4-1的驗證：

1）原始數據：

2）節點導納矩陣：

3、修改節點導納矩陣：

（1）□□增加一條線路：

1）原始數據：

2）節點導納矩陣：

（2）□□之間變比變為0.95；

1）原始數據：

2）節點導納矩陣：

（3）□□節點之間的阻抗變為j0.5o

1）原始數據：

2）節點導納矩陣：

4、程序流程圖：

七、實驗總結

附形成節點導納矩陣程序

n=input(*n-);

nl=input('nl=')；

isb=input(*isb-);

P5

pr=input('pr=');

Bl=input('Bl=')；

B2=input(,B2-);

X=input(fX=');

Y=zeros(n);e=zeros(l,n);f=zeros(1,n);V=zeros(1,n);O=zeros(l,n);S1=zeros(

nl)；

fori=l:n

ifX(i,2)f

p=X(i,l);

Y(p,p)=l./X(i,2);

end

end

fbri=l:nl

ifBl(i,6)==0

p=Bl(i,l);q=Bl(i,2);

elsep=Bl(i,2);q=Bl(i,l);

end

Y(p,q)=Y(p,q)-l./(Bl(i,3)*Bl(i,5));

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

Y(q,q)=Y(q,q)+l./(Bl(i,3)*B1(i,5)人2)+Bl(i,4)./2;

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

end

disp('Y=');

disp(Y)

P6

實驗二牛頓一拉夫遜法潮流計算實驗

一、本次課主要內容

1、牛頓一拉夫遜法的原理；

2、編寫牛頓一拉夫遜法的程序流程圖；

3、編寫牛頓一拉夫遜法的M程序的方法。

二、教學目的與要求

1、掌握牛頓一拉夫遜法的原理；

2、掌握編寫牛頓一拉夫遜法的程序流程圖；

3、掌握編寫牛頓一拉夫遜法的M程序的方法。

三、教學重點難點

1、重點是牛頓一拉夫遜法的程序流程圖；

2、難點是編寫牛頓一拉夫遜法的M程序的方法。

四、教學方法和手段

課堂講授、提問、討論、演示、實際操作等。

五、作業與習題布置

撰寫實驗報告

P7

實驗二牛頓一拉夫遜法潮流計算實驗

實驗時間：_______________________________

一.實驗目的

1、掌握牛頓一拉夫遜法的原理；

2、掌握編寫牛頓一拉夫遜法的程序流程圖；

3、掌握編寫牛頓一拉夫遜法的M程序的方法。

二,實驗內容與原理

1、牛頓一拉夫遜法基本原理

牛拉法作為求解非線性方程式的一種極其有效的方法，在數學上的地位很高。

其中最重要的一點就是將求解非線性方程問題轉變成為與之對應的線性方程的求解

過程，也就是我們俗稱的逐次線性化。

對于一個非線性方程組：

/U）=0即￡（不々,…，七）=。（/=1,2,（1）

x作為待求量，在它的任意一個初始估計值的周圍，將式子（1）按泰勒級數展

開，展開過程中注意略去高階項（二階及以上），從而可以得出以下經過線性化的

方程組：

/（x（o））+/（x（o））zkr（o）=o（2）

式（2）就是牛拉法中的修正方程式，根據這個我們就可以求出第一次迭代的修

正量：

堤（。）=_"'（。。））『/（。。））（3）

我們將和兩個量相加，從而就出現了，變量的第一次修改值。接下來我們拿著

的值，重復上面的計算過程。所以，從出發，運用牛拉法求解的迭代格式為：

/'（?幻）囚￡）=一/（系幻）（4）

/+d=/）+△?"）（5）

2、極坐標形式潮流方程式

、

△4=Q一巨UUj（G）cos%+Brsin0..}=P.-P；=Q

j=l

A

△2=2—巨UUj（G）sin%—Bucos%）=Q,—Q；=0

P8

3、設系統有n個節點，其中系統中有m個PQ節點，而除了PQ節點

和一個平衡節點外其余的都是PV節點。顯然，PV點數目為n-m-l。

在潮流計算中實際上需要求解的非線性方程組中包含n-1個有功功率

方程和m個無功功率方程，總共有n+m-1個。未知量有(n-1)個電

壓相角(i=l,2,3,…,m)和m個電壓幅值Ui(i=l,2,3,…,m),總數

為m+n-lo方程數和未知數相等，方程有解。全部節點的修正方程如

下：

■△八?△心

65

*

△后???:△a

?

?必/"

4P,一見一、.0U\

???=一?????????:????????????=0

弘QiaA<2,.MQ\64Q\△5

?Ou、?

????*

[△a,」

外(2,〃二2：S△以^Q.n

dU34」

LSO-x

4、迭代方程式:

一，”+】)--ef一△￡/)一

?

**

心)4")

???—???+???=o

U\k)△u*

?

*?

??

Z/(4+1)

m_s?

P9

5、收斂判據為:

max{|M，A0}YE

6、雅克比矩陣各分塊矩陣的非對角元素（iWj）分別為:

Hu=翳"=—UUj(G..sinBjjcos吟

Nij=~dULUj=~UiUjGcos4+B.sin6>.)

M“==LUj(G..cos%+Busin%)

LiJ=喘務。，=~~UUj(G..sin6>..

—B(.cosq,)

7、各分塊矩陣的對角元素(尸i)分別為：

AAry

Hn=''=-Ui￡Uj(Gjjsin%-與cosq)=+Q-

N"=答4=—U，W4Gcosq+嗎sin%)-2U：G,

jj*i

=-U；G「K

ML蜉=UgUi(dcosW+Bjjsin%)=u；G.「P；

0%尸i

4=箸4=MZUj(G,sine..-B.cos%)+2UH

=U；B「Q：

P10

8、應用計算機來進行潮流計算其求解修正方程收斂速度快、計算精度

高，可以得到滿意的計算結果。利用牛頓-拉夫遜法求解潮流的計算流

程如圖所示。

p11

三.實驗步驟

1.輸入原始數據

1）節點數n；

2）支路數nl；

3）支路參數矩陣Bl

它包括六個數據數j,z,b,t,it],i,j為支路兩端節點號,z為支路的阻抗,b為

線路電納,t為變比,it為高低壓側標志（高為L低為0）0

4）初始給定值矩陣B2

□節點所接發電機功率SG

□節點負荷功率SL

□節點電壓的初始值

□PV節點電壓的給定值

□節點所接無功補償設備的容量

□節點編號igl

平衡節點為1,PQ節點為2,PV節點為3

5）節點對地阻抗矩陣X（由節點號與接地阻抗構成）。

6）收斂精度pr

2.形成節點導納矩陣

3.根據極坐標形式潮流方程式，求出修正方程式的常數項向量

4.求J各元素值

5.解修正方程式得到修正量

6.修正各節點電壓。

7.按以上步驟迭代

8.求各支路潮流

四、實驗要求

對照設計內容和任務，設計相應的MATLAB的M文件，并調試運行得出正確

結果。

五、實驗儀器設備與器材

計算機（安裝有MATLAB軟件平臺，包含SIMULINK仿真模塊）。

P12

六、實驗結果（或數據）與分析

1、實驗一圖的潮流計算結果

1）原始數據：

2）計算結果：

2、（1）□□增加一條線路:

1）原始數據：

2）計算結果：

3、□□之間變比變為0.95；

1）原始數據：

2）計算結果：

4、□□節點之間的阻抗變為jO.5o

1）原始數據：

2）計算結果:

附程序

%本程序是用牛拉法進行潮流計算

n=input('n=,);

nl=input('nl-);

isb=input('isb=');

pr=input('pr=');

Bl=input('BI=')；

B2=input('B2=,);

X=input(,X=,);

Y=zeros(n);e=zeros(I,n);f=zeros(I,n);V=zeros(1,n);O=zeros(l,n);S1=zeros(nl);

fbri=l:n

P13

ifX(i,2)-=0;

P=X(i,I);

Y(p,p)=l./X(i,2);

end

end

fbri=l:nl

ifBl(i,6)=0

p=Bl(i,l);q=Bl(i,2);

elsep=Bl(i,2);q=Bl(i,l);

end

Y(p,q)=Y(p,q)-l./(Bl(i,3)*Bl(i,5));

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

Y(q,q尸Y(q,q)+l./(Bl(i,3)*Bl(i,5)A2)+Bl(i,4)./2;

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

end

G=real(Y);B=imag(Y);

fbri=l:n

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

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

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

end

fbri=l:n

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

B(i,i)=B(i,i)+B2(i,5);

end

P=rcal(S);Q=iinag(S);

ICTl=0;IT2=l;N0=2*n;N=N0+l;a=0;

whileIT2?=0

IT2=0;a=a+l;

fbri=l:n

ifi-isb

C(i)=0;

D(i)=0;

fbrjl=l:n

C(i)=C(i)+G(i,jl)*e(jl)-B(i,j

D(i)=D(i)+G(i,jl)*f(jl)+B(i,j

end

Pl=C(i)*e(i)+f(i)*D(i);

Ql=f(i)*C(i)-D(i)*e(i);

V2=e(i)A2+f{i}A2;

ifB2(i,6)-=3

DP=P(i)-Pl;

DQ=Q(i)-Ql;

fbrjl=l:n

ifjl-=isb&jl-=i

P14

X2=B(i,jl)*e(i)-G(i,jl)*ai);

X3=X2;

X4=-X1;

p=2*i-l;q=2*j1-1;J(p,q)=X3;J(p,N)=DQ;m=p+1;

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

elseifj1=i&jl-=isb

Xl=-C(i)-G(i,i)*e(i)-B(i,i)*f(i);

X2=-D(i)+B(i,i)*e(i)-G(i,i)*f(i);

X3=D(i)+B(i,i)*c(i)-G(i,i)*f(i);

X4=-C(i)+G(i,i)*e(i)+B(i,i)*f(i);

p=2*i-l；q=2*j1-1;J(p,q)=X3;J(p,N尸DQ;m=p+1;

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

end

end

else

DP=P(i)-Pl;

DV=V(i)A2-V2;

fbrjl=l:n

ifjl~=isb&jl~=i

Xl=-G(ijl)*c(i)-B(i,jl)*f(i)；

X2=B(ijl)*e(i)-G(ijl)*f(i);

X5=0;

X6=0;

p=2*i-l;q=2*j1-1;J(p,q)=X5;J(p,N)=DV;m=p+1;

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

clscifj1—i&j1Vsb

Xl=-C(i)-G(i,i)*e(i)-B(i,i)*f(i)；

X2=-D(i)+B(i,i)*e(i)-G(i,i)*f(i);

X5=-2*e(i);

X6=-2*f(i);

p=2*i-1;q=2*j1-1;J(p,q)=X5;J(p,N)=DV;m=p+1;

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

end

end

end

end

end

fork=3:N0

kl=k+l;Nl=N;

fork2=kl:Nl

J(k,k2)=J(k,k2)./J(k,k);

end

J(k,k)=l;

ifk~=3

k4=k-l;

fork3=3:k4

P15

Ibrk2=kl:NI

J(k3,k2)=J(k3,k2)-J(k3,k)*J(k,k2);

end

J(k3,k)=0;

end

ifk=NO,break;end

fbrk3=kl:N0

fork2=kl:NI

J(k3,k2)=J(k3,k2)-J(k3,k)*J(k,k2);

end

J(k3,k)=0;

end

else

fbrk3=kl:N0

fbrk2=kl:Nl

J(k3,k2)=J(k3,k2)-J(k3,k)*J(k,k2);

end

J(k3,k)=0;

end

end

end

fork=3:2:N0-l

L=(k+l)./2;

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

kl=k+l;

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

end

fork=3:N0

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

ifDET>=pr

IT2=IT2+1;

end

end

ICT2(a)=IT2;

ICT1=ICT1+1;

fbrk=l:n

dy(k)=sqrt(e(k)A2+fi(k)A2);

end

fbri=l:n

Dy(ICTl,i)=dy(i);

end

end

dispC迭代次數上

disp(ICTl);

dispC沒有達到精度要求的個數)；

disp(ICT2);

P16

fork=l:n

V(k)=sqrt(e(k)A2+nk)A2);

O(k)=atan(f(k)./c(k))*1XO./pi;

end

E=e+f*j;

disp。各節點的實際電壓標幺值E為(節點號從小到大排列)：,);

disp(E);

dis