潮流计算

节点数：4 支路数：4 计算精度：
支路1：+
1┠—————□—————┨3
支路2：+
1┠—————□—————┨4
支路3：+
2┠—————□—————┨4
支路4：+
3┠—————□—————┨4
节点1：PQ节点，S(1)=节点2：PQ节点，S(2)=节点3：PV节点，P(3)= V(3)= 节点4：平衡节点，U(4)=∠
摘要
运用matlab软件对选定课设题目进行潮流计算。

潮流计算是电力系统课程中必须掌握也是非常重要的计算。

潮流计算是指对电力系统正常运行状况的分析和计算。

在已知系统条件情况下，给定一些初始条件，进而计算出系统运行的电压和功率等；潮流计算方法很多：高斯-塞德尔法、牛顿-拉夫逊法、PQ分解法、直流潮流法等。

通过潮流计算，可以确定各母线的电压幅值和相角,各元件流过的功率和整个系统的功率损耗。

潮流计算是实现安全经济发供电的必要手段和重要工作环节。

因此潮流计算在电力系统的规划计算，生产运行，调度管理及科学计算中都有广泛的运用。

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

关键词：matlab 潮流计算 PQ分解法
目录
6
8
9
1.题目原始数据及其化简 原始数据：
节点数：4 支路数：4 计算精度： 支路 1： +
1┠—————□—————┨3 支路 2： +
1┠—————□—————┨4 支路 3： +
2┠—————□—————┨4 支路 4： +
3┠—————□—————┨4
节点1：PQ 节点，S(1)= 节点2：PQ 节点，S(2)= 节点3：PV 节点，P(3)= V(3)=
节点4：平衡节点，U(4)=∠
根据原始数据所画电路简化图如图1：
1
3
4
2
图1电路简化图
分解法
分解法基本思想
PQ 分解法是从改进和简化牛顿法潮流程序的基础上提出来的，它的基本思想是：把节点功率表示为电压向量的极坐标形式，以有功功率误差作为修正电压向量角度的依据，以无功功率误差作为修正电压幅值的依据，这样，n-1+m 阶的方程式便分解为一个n-1阶和一个m 阶的方程，这两组方程分别进行轮流迭代，这就是所谓的有功-无功功率分解法。

牛顿法潮流程序的核心是求解修正方程式，当节点功率方程式采取极坐标系统时，修正方程式为：
或展开为： /P H N V δ∆=⋅∆+⋅∆
V V L J Q V /∆⋅+∆⋅=∆∆⋅∆⋅δ
（1）
以上方程式是从数学上推倒出来的，并没有考虑电力系统这个具体对象的特点。

在交流高压电网中，输电线路的电抗要比电阻大得多，系统中有功功率变化主要受电压相位的影响，无功功率则主要受母线电压幅值变化的影响。

在修正方程式的系数矩阵中，偏导数/P V ∂∆∂和/Q δ∂∆∂的数值相对于偏导数/P δ
∂∆∂和/Q V ∂∆∂是相当小的，所以，矩阵N 及J 中各元素的数值相对是很小的，因
此对牛顿法的第一步简化就是把有功功率和无功功率分开来进行迭代，即将式(1)化简为：
/Q L V V ∆=-⋅∆
（2）
这样，由于我们把2n 阶的线性方程组变成了二个n 阶的线性方程组，因而大大节省了机器内存和解题时间。

但是矩阵H 和L 都是节点电压幅值和相角差的函数，在迭代过程中仍然不断变化，而且又都是不对称矩阵。

对牛顿法的第二个化简，也是比较关键的一个化简，即把式(2)中的系数矩阵简化为在迭代过程中不变的对称矩阵，即常数矩阵。

在一般情况下，线路两端电压的相角差是不大的(不超过10~20度)，因此可以认为：
cos 1sin ij ij ij ij
G B θθ≈=
（3）
此外，与系统各节点无功功率相应的导纳Li B 必定远远小于该节点自导纳的虚部，即：
或
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）
进一步可以把它们表示为以下矩阵的乘积：
1112111221
2212000
0n 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）
将以上两式的左右两侧用以下矩阵左乘 就可以得到
11
2
21111121222212212n 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 分解法迭代过程中基本计算公式
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编程及运行
PQ
是
否
是
否
否
图2 PQ 分解法潮流计算流程框图
源程序代码
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=z eros(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('电压迭代次数曲线');
运行程序及结果分析：
请输入节点数：n=4
请输入支路数：nl=4
请输入平衡母线节点号：isb=1
请输入误差精度：pr=
请输入由支路参数形成的矩阵：B1=[1 3 + 0 1 0;1 4 + 0 1 0;2 4 + 0 1 0;3 4 + 0 1 0]
请输入由支路参数形成的矩阵：B2=[0 0 0 2;0 0 0 2; 0 0 3;0 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为（节点号从小到大排）：
+ + +
各节点的电压V大小（节点号从小到大排）为：
各节点的电压相角O（节点号从小到大排）为：
各节点的功率S（节点号从小到大排）为：
+ + - +
各条支路的首端功率Sj（顺序同您输入B1时一样）为：
+
-
+
-
各条支路的末端功率Sj（顺序同您输入B1时一样）为：
-
+
-
+
各条支路的功率损耗DS （顺序同您输入B1时一样）为：
+
+
+
+
每次迭代后各节点的电压值如图所示
12345
6789
0.940.96
0.98
1
1.02
1.04
1.06
1.08
1.1
1.12
迭代次数电压电压迭代次数曲线
图3电压迭代次数曲线
经过九轮迭代，节点功率不平衡量下降到以下，迭代到此结束，电压幅值和相角都能够满足计算精度的要求。

4.小结
通过本次电力系统分析实验，我更加熟悉地掌握了matlab的用法，在之前也用过matlab进行过实验的分析，都没能十分熟练地运用该软件，这次实践过后，加深了对此软件的理解以及更多的应用。

PQ潮流计算法在运用matlab时会显得比较简单，相比别的软件。

就自动化专业来讲，matlab会是一个很方便去解决问题的软件，因此在之后的学习过程中应该更加重视它，更近一步的掌握对以后的专业学习会有很大帮助。

电力系统分析是今年学习的一门比较重要的专业课，电力系统潮流计算又是这门课非常重要的内容，因此此次课设意义重大。

PQ法是进行潮流计算的一种方法，通过给定条件确定各母线的电压幅值和相角,各元件流过的功率和整个系统的功率损耗。

在做课设的过程中也遇到过很多问题，比如对课本知识掌握不透没能完全理解实验内容耽误了一些时间和精力。

但在同学的帮助以及查阅资料的情况下慢慢的克服了个个问题，完成了实验报告。

感悟到无论做什么事都应该平心静气，认认真真。

5.参考文献
[1] 何仰赞等.电力系统分析[M]. 武汉：华中理工大学出版社，
[2] 西安交通大学等.电力系统计算[M].北京：水利电力出版社，
[3] 陈衍.电力系统稳态分析[M].北京：水利电力出版社，
[4] 李光琦.电力系统暂态分析[M].北京：水利电力出版社，
[5] 于永源，杨绮雯. 电力系统分析（第二版）[M]. 北京：中国电力出版社，
[6] 华智明.张瑞林.电力系统分析.北京.重庆大学出版社.
[7] 弋东方.电力工程电气设计手册.北京.中国电力出版社.
本科生课程设计成绩评定表
指导教师签字：。