潮流计算

武汉理工大学《电力系统分析》课程设计说明书节点数：4 支路数：4 计算精度：0.00010

支路1：0.0200+j0.0800

1┠—————□—————┨3

支路2：0.0400+j0.1200

1┠—————□—————┨4

支路3：0.0500+j0.1400

2┠—————□—————┨4

支路4：0.0400+j0.1200

3┠—————□—————┨4

节点1：PQ节点，S(1)=-0.6000-j0.2500

节点2：PQ节点，S(2)=-0.8000-j0.3500

节点3：PV节点，P(3)=0.4000 V(3)=0.9500

节点4：平衡节点，U(4)=1.0000∠0.0000

运用matlab软件对选定课设题目进行潮流计算。潮流计算是电力系统课程中必须掌握也是非常重要的计算。潮流计算是指对电力系统正常运行状况的分析和计算。在已知系统条件情况下，给定一些初始条件，进而计算出系统运行的电压和功率等；潮流计算方法很多：高斯-塞德尔法、牛顿-拉夫逊法、PQ分解法、直流潮流法等。

通过潮流计算，可以确定各母线的电压幅值和相角,各元件流过的功率和整个系统的功率损耗。潮流计算是实现安全经济发供电的必要手段和重要工作环节。因此潮流计算在电力系统的规划计算，生产运行，调度管理及科学计算中都有广泛的运用。

本课程设计采用PQ分解法进行电力系统分析的潮流计算程序的编制与调试，获得电力系统中各节点电压，为进一步进行电力系统分析作准备。

关键词：matlab 潮流计算PQ分解法

1.题目原始数据及其化简 (1)

2.PQ分解法 (2)

2.1PQ分解法基本思想 (2)

2.2 PQ分解法潮流计算基本步骤 (5)

3编程及运行 (6)

3.1 PQ分解法潮流计算程序框图 (6)

3.2源程序代码 (7)

3.3运行程序及结果分析： (16)

4.小结 (18)

5.参考文献 (19)

1.题目原始数据及其化简

原始数据：

节点数：4 支路数：4 计算精度：0.00010 支路1：0.0200+j0.0800

1┠—————□—————┨3

支路2：0.0400+j0.1200

1┠—————□—————┨4

支路3：0.0500+j0.1400

2┠—————□—————┨4

支路4：0.0400+j0.1200

3┠—————□—————┨4

节点1：PQ节点，S(1)=-0.6000-j0.2500

节点2：PQ节点，S(2)=-0.8000-j0.3500

节点3：PV节点，P(3)=0.4000 V(3)=0.9500

节点4：平衡节点，U(4)=1.0000∠0.0000

根据原始数据所画电路简化图如图1：

1 3

4 2

图1电路简化图

2.PQ 分解法

2.1PQ 分解法基本思想

PQ 分解法是从改进和简化牛顿法潮流程序的基础上提出来的，它的基本思想是：把节点功率表示为电压向量的极坐标形式，以有功功率误差作为修正电压向量角度的依据，以无功功率误差作为修正电压幅值的依据，这样，n-1+m 阶的方程式便分解为一个n-1阶和一个m 阶的方程，这两组方程分别进行轮流迭代，这就是所谓的有功-无功功率分解法。

牛顿法潮流程序的核心是求解修正方程式，当节点功率方程式采取极坐标系统时，修正方程式为：

/P H N Q J L V V δ????????

=-??????????????

或展开为： /P H N V δ?=??+??

V V L J Q V /??+??=?????δ

（1）

以上方程式是从数学上推倒出来的，并没有考虑电力系统这个具体对象的特点。

在交流高压电网中，输电线路的电抗要比电阻大得多，系统中有功功率变化主要受电压相位的影响，无功功率则主要受母线电压幅值变化的影响。在修正方程式的系数矩阵中，偏导数/P V ???和/Q δ???的数值相对于偏导数/P δ???和

/Q V ???是相当小的，所以，矩阵N 及J 中各元素的数值相对是很小的，因此对

牛顿法的第一步简化就是把有功功率和无功功率分开来进行迭代，即将式(1)化简为：

P H δ?=-?

/Q L V V ?=-??

（2）

这样，由于我们把2n 阶的线性方程组变成了二个n 阶的线性方程组，因而大大节省了机器内存和解题时间。但是矩阵H 和L 都是节点电压幅值和相角差的函数，在迭代过程中仍然不断变化，而且又都是不对称矩阵。对牛顿法的第二个化简，也是比较关键的一个化简，即把式(2)中的系数矩阵简化为在迭代过程中不变的对称矩阵，即常数矩阵。

在一般情况下，线路两端电压的相角差是不大的(不超过10~20度)，因此可以认为：

cos 1ij δ≈

cos 1sin ij ij ij ij

G B θθ≈= （3）

此外，与系统各节点无功功率相应的导纳Li B 必定远远小于该节点自导纳的虚部，即：

2i

Li ii i

Q B B V =

=

或

2i i ii Q V B =

（4）

考虑到以上关系后，式(2)中系数矩阵中的元素表达式可以化简为：

22ii i ii ij i j ij ii i ii

ij i j ij

H V B H VV B L V B L VV B ==== （5）

这样，式(5)中系数矩阵可以表示为：

211112121122221212222

1122n n n n n n n n n nn V B VV B VV B V V B V V B V B H L V V B

V V B V B ??

?

?

== ?

? ???L L M M

M L

（6）

进一步可以把它们表示为以下矩阵的乘积：

1112111221221

20000n n n n n n nn B B B V V B B B H L V V B B B ??

???? ? ?

?

?== ? ?

? ? ? ?????

??L L O O M M M L

（7）

将它代入(2)中，并利用乘法结合率，可以把修正方程式变为：

1111111212222221221200n n n n n n n n nn V P V B B B V P B V B B V P V B B B θθθ??? ???????

? ? ? ?? ? ? ?

?= ? ? ? ? ? ? ? ??????

?? ???

V L V L O M M

M M M V L

（8）

及

11111121222221221200n n n n n n n nn V Q V B B B V Q B V B B V Q V B B B ???

???????

? ? ? ?? ? ? ? ?= ? ? ?

? ? ? ? ??????

?? ???

V L V L O M M M M M V L

（9）

将以上两式的左右两侧用以下矩阵左乘

????

??? ??=??????? ??n n V V V V V V /10...

/10/10 (02121)

就可以得到

1

1

22

1111121222212212n n P V n P n V P n n n n nn V V B B B B V B B V B B B θθθ?????????? ? ??? ??

???

= ? ??? ? ??? ?? ???????

L L M

M M

M M L

（10）

112211112122212212n n Q V n Q n V Q n n n nn V V B B B B V B B V B B B ?????

????? ? ??? ??

???= ? ???

? ??? ?? ???

????

L L M M M

M M L

（11）

以上两式就是P-Q 分解法达到修正方程式，其中系数矩阵只不过是系统导纳矩阵的虚部，只是阶次不同，是对称矩阵，而且在迭代过程中维持不变。它们与功率误差方程式

()

)

3,2,1(sin cos 1

n i B G V V P P ij ij ij ij n

j j j i is i Λ=+-=?∑==θθ （12）

()

)

3,2,1(cos sin 1

n i B G V V Q Q ij ij ij ij n j j j i is i Λ=--=?∑==θθ （13）

构成了P-Q 分解法迭代过程中基本计算公式

2.2 PQ 分解法潮流计算基本步骤

1) 形成系数矩阵B B '''、，并求其逆矩阵。

2) 设PQ 节点电压的初值和各节点相角初值为(0)i U (i=1,2，…，n ，i ≠s)和(0)i δ (i=1,2，…，m ，i ≠s)。

3) 通过（12）式计算各节点有功功率误差(0)i

P ?,从而求出

(0)

i (0)

i U P ?

(i=1,2，…，n ，i ≠s)。

4) 解修正方程式，求各节点电压相位角的变量(0)i δ?(i=1,2，…，n ，i ≠s) 5) 求各节点电压相位角的新值(0)i (0)i (1)i δδδ?+=(i=1,2，…，n ，i ≠s)。

6) 通过（13）式计算无功功率的不平衡量(0)

i Q ?，从而求出(0)

i (0)i U Q ?

(i=1,2，…，m ，i ≠s)。

7) 解修正方程式，求各节点电压大小的变量(0)i U ?(i=1,2，…，m ，i ≠s)。 8) 求各节点电压大小的新值(0)i (0)i (1)i U U U ?+=(i=1,2，…，m ，i ≠s)。 9) 运用各节点电压的新值自第三步开始进入下一次迭代。 10) 计算平衡节点功率和线路功率

3编程及运行

图2 PQ分解法潮流计算流程框图

3.2源程序代码

n=input('请输入节点数：n=');

nl=input('请输入支路数：nl=');

isb=input('请输入平衡母线节点号：isb=');

pr=input('请输入误差精度：pr=');

B1=input('请输入由支路参数形成的矩阵：B1='); %输入B1

B2=input('请输入由支路参数形成的矩阵：B2='); %输入B2

X=input('请输入由节点号及其对地阻抗形成的矩阵：X=');%输入X

na=input('请输入PQ节点数na=');

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

for i=1:n

if X(i,2)~=0;

p=X(i,1);

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

end

end

for i=1:nl

if B1(i,6)==0

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

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

end

Y(p,q)=Y(p,q)-1./(B1(i,3)*B1(i,5));

YI(p,q)=YI(p,q)-1./B1(i,3);

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

YI(q,p)=YI(p,q);

Y(q,q)=Y(q,q)+1./(B1(i,3)*B1(i,5)^2)+B1(i,4)./2;

YI(q,q)=YI(q,q)+1./B1(i,3);

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

YI(p,p)=YI(p,p)+1./B1(i,3);

end %求导纳矩阵

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

for i=1:n

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

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

end

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

for i=1:n

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

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

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

end

for i=1:n

if B2(i,6)==2

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

O(i)=atan(f(i)./e(i));

end

end

for i=2:n

if i==n

B(i,i)=1./B(i,i);

else IC1=i+1;

for j1=IC1:n

B(i,j1)=B(i,j1)./B(i,i);

end

B(i,i)=1./B(i,i);

for k=i+1:n

for j1=i+1:n

B(k,j1)=B(k,j1)-B(k,i)*B(i,j1);

end

end

end

end

p=0;q=0;

for i=1:n

if B2(i,6)==2

p=p+1;k=0;

for j1=1:n

if B2(j1,6)==2

k=k+1;

A(p,k)=BI(i,j1);

end

end

end

end

for i=1:na

if i==na

A(i,i)=1./A(i,i);

else k=i+1;

for j1=k:na

A(i,j1)=A(i,j1)./A(i,i);

end

A(i,i)=1./A(i,i);

for k=i+1:na

for j1=i+1:na

A(k,j1)=A(k,j1)-A(k,i)*A(i,j1);

end

end

end

end

ICT2=1;ICT1=0;kp=1;kq=1;K=1;DET=0;ICT3=1;

while ICT2~=0|ICT3~=0

ICT2=0;ICT3=0;

for i=1:n

if i~=isb

C(i)=0;

for k=1:n

C(i)=C(i)+V(k)*(G(i,k)*cos(O(i)-O(k))+BI(i,k)*sin(O(i)-O(k)));

end

DP1(i)=P(i)-V(i)*C(i);

DP(i)=DP1(i)./V(i);

DET=abs(DP1(i));

if DET>=pr

ICT2=ICT2+1;

end

end

end

Np(K)=ICT2;

if ICT2~=0

for i=2:n

DP(i)=B(i,i)*DP(i);

if i~=n

IC1=i+1;

for k=IC1:n

DP(k)=DP(k)-B(k,i)*DP(i);

end

else

for LZ=3:i

L=i+3-LZ;

IC4=L-1;

for MZ=2:IC4

I=IC4+2-MZ;

DP(I)=DP(I)-B(I,L)*DP(L);

end

end

end

end

for i=2:n

O(i)=O(i)-DP(i);

end

kq=1;L=0;

for i=1:n

if B2(i,6)==2

C(i)=0;L=L+1;

for k=1:n

C(i)=C(i)+V(k)*(G(i,k)*sin(O(i)-O(k))-BI(i,k)*cos(O(i)-O(k)));

end

DQ1(i)=Q(i)-V(i)*C(i);

DQ(L)=DQ1(i)./V(i);

DET=abs(DQ1(i));

if DET >=pr

ICT3=ICT3+1;

end

end

end

else kp=0;

if kq~=0;

L=0;

for i=1:n

if B2(i,6)==2

C(i)=0;L=L+1;

for k=1:n

C(i)=C(i)+V(k)*(G(i,k)*sin(O(i)-O(k))-BI(i,k)*cos(O(i)-O(k)));

end

DQ1(i)=Q(i)-V(i)*C(i);

DQ(L)=DQ1(i)./V(i);

DET=abs(DQ1(i));

end

end

end

end

Nq(K)=ICT3;

if ICT3~=0

L=0;

for i=1:na

DQ(i)=A(i,i)*DQ(i);

if i==na

for LZ=2:i

L=i+2-LZ;

IC4=L-1;

for MZ=1:IC4

I=IC4+1-MZ;

DQ(I)=DQ(I)-A(I,L)*DQ(L);

end

end

else

IC1=i+1;

for k=IC1:na

DQ(k)=DQ(k)-A(k,i)*DQ(i);

end

end

end

L=0;

for i=1:n

if B2(i,6)==2

L=L+1;

V(i)=V(i)-DQ(L);

end

end

kp=1;

K=K+1;

else

kq=0;

if kp~=0

K=K+1;

end

end

for i=1:n

Dy(K-1,i)=V(i);

end

end

disp('迭代次数');

disp(K);

disp('每次没有达到精度要求的有功功率个数为'); disp(Np);

disp('每次没有达到精度要求的无功功率个数为'); disp(Nq);

for k=1:n

E(k)=V(k)*cos(O(k))+V(k)*sin(O(k))*j;

O(k)=O(k)*180./pi;

end

disp('各节点的电压标幺值E为（节点号从小到大排）：');

disp(E);

disp('各节点的电压V大小（节点号从小到大排）为：');

disp(V);

disp('各节点的电压相角O（节点号从小到大排）为：');

disp(O);

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

end

disp('各节点的功率S（节点号从小到大排）为：');

disp(S);

disp('各条支路的首端功率Sj（顺序同您输入B1时一样）为：');

for i=1:nl

if B1(i,6)==0

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

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

end

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

disp(Si(p,q));

end

disp('各条支路的末端功率Sj（顺序同您输入B1时一样）为：');

for i=1:nl

if B1(i,6)==0

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

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

end

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

disp(Sj(q,p));

end

disp('各条支路的功率损耗DS（顺序同您输入B1时一样）为：');

for i=1:nl

if B1(i,6)==0

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

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

end

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

disp(DS(i));

end

for i=1:K

Cs(i)=i;

for j=1:n

Dy(K,j)=Dy(K-1,j);

end

end

disp('以下是每次迭代后各节点的电压值（如图所示）');

plot(Cs,Dy),xlabel('迭代次数'),ylabel('电压'),title('电压迭代次数曲线');

3.3运行程序及结果分析：

请输入节点数：n=4

请输入支路数：nl=4

请输入平衡母线节点号：isb=1

请输入误差精度：pr=0.00001

请输入由支路参数形成的矩阵：B1=[1 3 0.02+0.08i 0 1 0;1 4 0.04+0.12i 0 1 0;2 4 0.05+0.14i 0 1 0;3 4 0.04+0.12i 0 1 0]

请输入由支路参数形成的矩阵：B2=[0 -0.6-0.25i 1.0 0 0 2;0 -0.8-0.35i 1.0 0 0 2;0.4 0 0.95 0.95 0 3;0 0 1.0 1.0 0 1]

请输入由节点号及其对地阻抗形成的矩阵：X=[1 0;2 0;3 0;4 0]

请输入PQ节点数na=2

迭代次数9

每次没有达到精度要求的有功功率个数为

3 3 3 3 3 2 2 2 0

每次没有达到精度要求的无功功率个数为

2 2 2 2 1 1 1 0 0

各节点的电压标幺值E为（节点号从小到大排）：

0.9641 1.0658 + 0.1659i 0.9484 + 0.0546i 0.9978 + 0.0666i

各节点的电压V大小（节点号从小到大排）为：

0.9641 1.0787 0.9500 1.0000

各节点的电压相角O（节点号从小到大排）为：

0 8.8456 3.2936 3.8196

各节点的功率S（节点号从小到大排）为：

-1.1374 + 0.2500i 0.8000 + 0.3500i 0.4000 - 0.6290i 0.0000 + 0.2198i

各条支路的首端功率Sj（顺序同您输入B1时一样）为：

-0.5746 + 0.3328i

-0.5628 - 0.0828i

0.8000 + 0.3500i

-0.1841 - 0.3341i

各条支路的末端功率Sj （顺序同您输入B1时一样）为： 0.5841 - 0.2948i 0.5767 + 0.1246i -0.7672 - 0.2583i 0.1905 + 0.3535i

各条支路的功率损耗DS （顺序同您输入B1时一样）为： 0.0095 + 0.0379i 0.0139 + 0.0418i 0.0328 + 0.0917i 0.0065 + 0.0194i

每次迭代后各节点的电压值如图所示

1234

56789

0.94

0.960.9811.021.041.061.081.1

1.12迭代次数

电压

电压迭代次数曲线

图3电压迭代次数曲线

经过九轮迭代，节点功率不平衡量下降到0.000010以下，迭代到此结束，电压幅值和相角都能够满足计算精度的要求。