课程设计 4支路,4节点简单闭式网络潮流计算及程序设计

1 设计任务及要求
1.1设计任务
4支路，4节点简单闭式网络潮流计算及程序设计。

本任务要求用MATLAB 编写程序，包括计算原理，计算程序，计算结果，结果分析等。

实验原始数据如下：
节点数：4 支路数：4 计算精度：0.00010
支路1： 3 (0.0200+j0..0800)
支路2： 4 （0.0400+j0.1200)
支路3： 4 （0.0500+j0.1400）
支路4： 4 （0.0400+j0.1200)
节点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.2 数据处理
根据原始数据所画电路简化图如图1
1 3
图1电路简化图
2 理论方法
2.1 PQ 分解法思路
PQ 分解法即有功-无功功率分解法，它的基本思路是：把节点功率表示为电压向量的极坐标形式，然后以有功功率误差作为修正电压向量角度的依据，以无功功率误差作为修正电压幅值的依据，这样，n-1+m 阶的方程式便分解为一个n-1阶和一个m 阶的方程，这两组方程分别进行轮流迭代，这就是所谓的PQ 分解法。

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

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

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

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

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

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

考虑到以上关系后，式(2)中系数矩阵中的元素表达式可以化简为：
（3）
这样，式(5)中系数矩阵可以表示为：
（4）
进一步可以把它们表示为以下矩阵的乘积：
（5）
将它代入(2)中，并利用乘法结合率，可以把修正方程式变为：
（6）
及
（7）
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 ====21111212
1122221212222
1122
n 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 ⎛⎫
⎪
⎪
== ⎪
⎪ ⎪⎝
⎭1112
111221
2212
0000n n n n n n nn B B B V V B B B H L V V B B B ⎛⎫⎛⎫⎛⎫ ⎪ ⎪ ⎪
⎪== ⎪ ⎪
⎪ ⎪ ⎪ ⎪⎝
⎭⎝
⎭
⎝
⎭1111111212222221
221
2
00n 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 θθθ⎛⎫∆ ⎪⎛⎫⎛⎫⎛⎫
⎪ ⎪ ⎪ ⎪∆ ⎪ ⎪ ⎪ ⎪= ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪∆⎝
⎭⎝⎭
⎝⎭ ⎪⎝
⎭
11111
121222221
221
2
00n 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 ⎛⎫∆ ⎪⎛⎫⎛⎫⎛⎫
⎪ ⎪ ⎪ ⎪∆ ⎪ ⎪ ⎪ ⎪= ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪∆⎝
⎭⎝⎭
⎝⎭ ⎪⎝
⎭
将以上两式的左右两侧用以下矩阵左乘
(8)
就可以得到
（9）
（10）
在这两个修正方程式中，系数矩阵都由节点导纳矩阵的虚部构成，只是阶次不同，是对称矩阵，而且在迭代过程中维持不变。

它们与功率误差方程式
()
)
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 =+-=∆∑==θθ （11）
()
)
3,2,1(cos sin 1n i B G V V Q Q ij ij ij ij n
j j j i is i =--=∆∑==θθ （12）
构成了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)。

⎪⎪⎪⎪
⎪⎭⎫ ⎝⎛=⎪⎪⎪⎪⎪⎭⎫ ⎝⎛n n V V V V V V /10...
/10/10 (02121)
1
1
2
21111
12122221221
2n 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 θθθ∆∆∆⎛⎫∆⎛⎫⎛⎫ ⎪ ⎪⎪ ⎪∆ ⎪⎪= ⎪
⎪⎪ ⎪ ⎪⎪ ⎪∆ ⎪⎝⎭⎝⎭⎝
⎭
11221111212221221
2
n n Q V n Q n V Q n n n nn V B B B B V B B V B B B ∆∆∆⎛⎫
∆⎛⎫⎛⎫ ⎪ ⎪⎪ ⎪∆ ⎪⎪= ⎪
⎪⎪ ⎪ ⎪⎪ ⎪∆ ⎪⎝⎭
⎝⎭⎝
⎭
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 MATLAB编程及运行
图2 PQ分解法潮流计算流程框图
3.2运行程序及结果分析：
请输入节点数：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
每次迭代后各节点的电压值如图所示
图3电压迭代次数曲线
从图中可以看出，经过九轮迭代，节点功率不平衡量下降到0.000010以下，迭代结束，电压幅值和相角都能够满足计算精度为0.000010的要求，所以完成任务。

1
2 3 4
5 6 7 8 9
11.1迭代次数
电压
4.小结
在这次课程设计中，我的题目是四节点四支路简单闭式网络潮流计算程序设计。

我对整个设计的思路是先从理论知识着手，我先从图书馆借阅了大量的参考书并详细研读了课本上关于潮流计算以及PQ分解法的知识，确定了做题的基本思路。

接下来就是运用MATLAB进行编程计算，这对我来讲是一个非常难的问题，因为没有这方面的实践，所以一开始根本不知道从哪里入手写程序，后来我与同学们还有几个学长交流之后，终于确定了思路，虽然很麻烦但是终于可以解决问题了。

通过这次课程设计，我明白了许多，学到了许多知识。

首先，MATLAB软件的应用非常重要，尤其是对于我们电气工程及其自动化专业的学生来讲，要学会熟练的运用它来解决各种各样的问题，而且不能停留在理论方面，要多上机操作，并反复思考，孔子云“学而不思则罔，思而不学则殆”告诉我们要熟练的运用必须学习理论并且还要多思考。

其次，做好一门课程设计要从整个全局去规划，比如首先确定理论再实践，只有从全局规划好了，才能少走很多弯路。

我觉得跟同学们以及学长们的交流非常重要，当你遇到问题的时候而自己又解决不了，通过和同学们的讨论交流各自的想法就能很快的解决问题。

在这次课程设计中，我认识了专业知识需要学习的有很多，通过这次课程设计，我知道了PQ分解法就是牛顿-拉夫逊法的简化算法，并且简化也是有条件的，并且熟悉了学过的线性代数中的三角分解。

回顾整个从无到有的学习，各种困难摆在我的面前，虽然也有烦躁，但是经过一系列的学习我终于完成了这次课程设计。

在7天的时间里，我认识到了自己的不足，也认识到了主动学习的重要性，并增长了自己的知识，同时也意识到团队合作的重要性，通过交流我很快的解决了问题。

最后我要感谢在这次课程设计中帮助我的同学们以及几位学长，没有他们的帮助我很难完成，总的来说，这次课程设计我不仅从学习上收获了血多，从生活中也有了很大的收获。

5参考文献
[1] 何仰赞等.电力系统分析[M]. 武汉：华中科技大学出版社，2002.3
[2] 陈珩,电力系统稳态分析[M].北京:中国电力出版社,1999
[3]胡博.基于预处理条件GMRES的不精确牛顿法潮流计算[J].电工技术学
报,2007,22(2):99-104.
[4] 王沫然MATLAB6.0与科学计算[M].北京：电子工业出版社，2001.1
[5] 张伯明，陈寿孙高等电力网络分析[M].北京：清华大学出版社，1996
附录
程序代码：
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
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;
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
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('电压迭代次数曲线');。