一、实验学时：
6学时。
二、实验目的：
通过采用牛顿-拉夫逊法实现电力系统潮流计算的编程和仿真实验，强化学生对复杂电力系统潮流计算相关知识的理解，使学生具备通过MATLAB编程实现数值计算的能力，培养学生解决电力系统中复杂工程问题的能力。
三、实验原理：
电力系统分析的潮流计算是电力系统分析的一个重要的部分。通过对电力系统潮流分布的分析和计算，可进一步对系统运行的安全性，经济性进行分析、评估，提出改进措施。电力系统潮流的计算和分析是电力系统运行和规划工作的基础。潮流计算是指对电力系统正常运行状况的分析和计算。通常需要已知系统参数和条件，给定一些初始条件，从而计算出系统运行的电压和功率等；潮流计算方法很多：高斯-塞德尔法、牛顿-拉夫逊法、P-Q分解法、直流潮流法，以及由高斯-塞德尔法、牛顿-拉夫逊法演变的各种潮流计算方法。
采用牛顿-拉夫逊迭代法实现潮流计算的一般步骤：
(1) 输入原始数据和信息：y、Pis、Qis、Uis、约束条件；
(2) 形成节点导纳矩阵YB；
(3) 设置各节点电压初值ei(0)、fi(0)或Ui(0)、δi(0)；
(4) 将初始值代入直角坐标或极坐标形式的功率方程，求不平衡量ΔPi(0)、ΔQi(0)、ΔUi2(0)；
(5) 计算雅可比矩阵各元素（Hij、Lij、Nij、Jij、Rij、Sij）；
(6) 求解修正方程，解得Δei(k)、Δfi(k)或ΔUi(k)、Δδi(k)；
(7) 求节点电压新值ei(k+1) = ei(k) + Δei(k)、fi(k+1) = fi(k) + Δfi(k)，或Ui(k+1) = Ui(k) + ΔUi(k)、δi(k+1) = δi(k) + Δδi(k)；
(8) 判断是否收敛：Max|Δei(k)| ≤ ε，Max|Δfi(k)| ≤ ε或Max|ΔUi(k)| ≤ ε，Max|Δδi(k)| ≤ ε；
(9) 重复迭代步骤(4)、(5)、(6)、(7)，直到满足步骤(8)的收敛条件；
(10) 求平衡节点的功率和PV节点的Qi及各支路的功率。
四、实验要求：
1. 通过牛顿-拉夫逊法求非线性方程组近似解的MATLAB编程范例，理解实现牛顿-拉夫逊法的基本代码编写方法；
2. 根据给定的电力系统网络接线图、节点类型和具体参数，运用以极坐标形式的牛顿-拉夫逊法计算系统的潮流分布。
五、实验内容：
1. 牛顿-拉夫逊法求非线性方程组近似解范例
求解过程：
令
迭代次数为k。
F(X)的雅克比矩阵为
设初始近似解为
迭代精度取0.0001。
求解代码示例：
clear
x(1)=1.0;x(2)=2.0;
k=0; precision=1;
k,x
while precision>0.0001
f1=3*x(1)^2+2*x(2)^2+x(1)*x(2)-x(2)-10.5;
f2=2*x(1)^2+x(2)^2+2*x(1)*x(2)+x(1)-11.3;
f=[f1 f2]'
k=k+1;
k
J=[6*x(1)+x(2) x(1)+4*x(2)-1
4*x(1)+2*x(2)+1 2*x(1)+2*x(2)];
dx=-J\f;
x(1)=x(1)+dx(1);
x(2)=x(2)+dx(2);
x
precision=max(abs(dx));
end2. 采用牛顿-拉夫逊法实现电力系统潮流计算编程
网络接线如图1.1所示，各支路导纳均以标幺值标于图1.1中。其中：
(1) 节点1、2、3、4为PQ节点，注入功率分别为：
，节点1连接给定功率的发电厂；
(2) 节点5为平衡节点，电压保持为定值，V5 = 1.05；
试运用极坐标形式的牛顿-拉夫逊法计算该系统各节点的电压和各线路的功率。计算精度要求个节点电压修正量不大于10?5。
图1.1 电力系统接线图
程序编写提示：
(1) 注意编写程序时节点编号应与图中对应，特别是平衡节点必须编为5号；
(2) 在MATLAB中i和j是作为虚数单位，所以在编写代码时表示节点导纳矩阵的行号和列号的变量用m和n；
(3) 极坐标形式的牛顿-拉夫逊法潮流计算相关公式：
(4) 节点参数和导纳矩阵相关代码：
clear
G(1,1)=10.2;B(1,1)=-31.5;G(1,2)=-1.2;B(1,2)=4.0;G(1,3)=-1.5;
B(1,3)=5.0;G(1,4)=-2.5;B(1,4)=7.5;G(1,5)=-5.000;B(1,5)=15.000;
G(2,1)=-1.2;B(2,1)=4.0; G(2,2)=10.4;B(2,2)=-31.7;G(2,3)=-8.0;
B(2,3)=24.0;G(2,4)=0;B(2,4)=0;G(2,5)=-1.2;B(2,5)=3.7;
G(3,1)=-1.5;B(3,1)=5.0;G(3,2)=-8.0;B(3,2)=24.0;G(3,3)=10.7;B(3,3)=-32.7; G(3,4)=-1.2;B(3,4)=3.7;G(3,5)=0;B(3,5)=0;
G(4,1)=-2.500;B(4,1)=7.500;G(4,2)=0;B(4,2)=0;G(4,3)=-1.2;B(4,3)=3.7;G(4,4)=3.7;B(4,4)=-11.2;G(4,5)=0;B(4,5)=0;
G(5,1)=-5.0;B(5,1)=15.0;G(5,2)=-1.2;B(5,2)=3.7;G(5,3)=0;B(5,3)=0;G(5,4)=0;B(5,4)=0;G(5,5)=6.2;B(5,5)=-18.7;
Y=G+j*B;
delt(1)=0;delt(2)=0;delt(3)=0;delt(4)=0;
u(1)=1.0;u(2)=1.0;u(3)=1.0;u(4)=1.0;
p(1)=0.20; q(1)=0.20; p(2)=-0.45; q(2)=-0.15;
p(3)=-0.40; q(3)=-0.05; p(4)=-0.60; q(4)=-0.10;
d(1,4)=0; d(4,1)=0; d(1,5)=0;d(5,1)=0;
k=0;precision=1;
k,delt,u
N1=4;
while precision>0.00001
delt(5)=0;u(5)=1.05;
for m=1:N1
for n=1:N1+1
pt(n)=u(m)*u(n)*(G(m,n)*cos(delt(m)-delt(n))+B(m,n)*sin(delt(m)-delt(n)));
qt(n)=u(m)*u(n)*(G(m,n)*sin(delt(m)-delt(n))-B(m,n)*cos(delt(m)-delt(n)));
end
pp(m)=p(m)-sum(pt);qq(m)=q(m)-sum(qt);
end
pp,qq
for m=1:N1
for n=1:N1+1
h0(n)=u(m)*u(n)*(G(m,n)*sin(delt(m)-delt(n))-B(m,n)*cos(delt(m)-delt(n)));
n0(n)=-u(m)*u(n)*(G(m,n)*cos(delt(m)-delt(n))+B(m,n)*sin(delt(m)-delt(n)));
j0(n)=-u(m)*u(n)*(G(m,n)*cos(delt(m)-delt(n))+B(m,n)*sin(delt(m)-delt(n)));
l0(n)=-u(m)*u(n)*(G(m,n)*sin(delt(m)-delt(n))-B(m,n)*cos(delt(m)-delt(n)));
end
H(m,m)=sum(h0)-u(m)^2*(G(m,m)*sin(delt(m)-delt(m))-B(m,m)*cos(delt(m)-delt(m)));
N(m,m)=sum(n0)-2*u(m)^2*G(m,m)+u(m)^2*(G(m,m)*cos(delt(m)-delt(m)))+B(m,m)*sin(delt(m)-delt(m));
J(m,m)=sum(j0)+u(m)^2*(G(m,m)*cos(delt(m)-delt(m))+B(m,m)*sin(delt(m)-delt(m)));
L(m,m)=sum(l0)+2*u(m)^2*B(m,m)+u(m)^2*(G(m,m)*sin(delt(m)-delt(m)) -B(m,m)*cos(delt(m)-delt(m)));
end
for m=1:N1-1
JJ(2*m-1,2*m-1)=H(m,m);JJ(2*m-1,2*m)=N(m,m);
JJ(2*m,2*m-1)=J(m,m);JJ(2*m,2*m)=L(m,m);
end
for m=N1:N1
JJ(2*m-1,2*m-1)=H(m,m);
end
for m=1:N1
for n=1:N1
if m==n
else
H(m,n)=-u(m)*u(n)*(G(m,n)*sin(delt(m)-delt(n))-B(m,n)*cos(delt(m)-delt(n)));
J(m,n)=u(m)*u(n)*(G(m,n)*cos(delt(m)-delt(n))+B(m,n)*sin(delt(m)-delt(n)));
N(m,n)=-J(m,n);L(m,n)=H(m,n);
end
end
end
for m=1:N1-1
for n=1:N1-1
if m==n
else
JJ(2*m-1,2*n-1)=H(m,n);JJ(2*m-1,2*n)=N(m,n);
JJ(2*m,2*n-1)=J(m,n);JJ(2*m,2*n)=L(m,n);
end
end
end
for m=N1
for n=1:N1-1
JJ(2*m-1,2*n-1)=H(m,n);JJ(2*m-1,2*n)=N(m,n);
end
end
for n=N1
for m=1:N1-1
JJ(2*m-1,2*n-1)=H(m,n);JJ(2*m,2*n-1)=J(m,n);
end
end
for m=1:N1-1
PP(2*m-1)=pp(m);PP(2*m)=qq(m);
end
for m=N1
PP(2*m-1)=pp(m);
end
uu=-inv(JJ)*PP';precision=max(abs(uu));uu
for n=1:N1-1
delt(n)=delt(n)+uu(2*n-1);
u(n)=u(n)+uu(2*n);
end
for n=N1
delt(n)=delt(n)+uu(2*n-1);
end
k=k+1;
k,delt,u
end
for n=1:N1+1
U(n)=u(n)*(cos(delt(n))+j*sin(delt(n)));
end
for m=1:N1+1
I(m)=Y(5,m)*U(m);
end
S5=U(5)*sum(conj(I));
for n=1:N1+1
q4(n)=u(4)*u(n)*(G(4,n)*sin(delt(4)-delt(n))-B(4,n)*cos(delt(4)-delt(n)));
end
Q4=sum(q4)
for m=1:N1+1
for n=1:N1+1
S(m,n)=U(m)*(conj(U(m))*conj(d(m,n))+(conj(U(m))-conj(U(n)))*conj(-Y(m,n)));
end
end
Y
JJ
S
B
pp
qq
uu
U
k
Q4
S5
运行结果：
>> dianli2
k =
0
delt =
0 0 0 0
u =
1 1 1 1
pp =
0.4500 -0.3900 -0.4000 -0.6000
qq =
0.9500 0.0350 -0.0500 -0.1000
uu =
-0.0476
0.0329
-0.0903
0.0041
-0.0967
0.0032
-0.1097
k =
1
delt =
-0.0476 -0.0903 -0.0967 -0.1097 0
u =
1.0329 1.0041 1.0032 1.0000 1.0500
pp =
-0.0306 -0.0002 0.0078 0.0101
qq =
-0.0751 -0.0217 -0.0095 -0.0326
uu =
0.0001
-0.0032
0.0008
-0.0036
0.0008
-0.0033
0.0000
k =
2
delt =
-0.0475 -0.0896 -0.0959 -0.1097 0
u =
1.0297 1.0005 0.9998 1.0000 1.0500
pp =
-0.0012 0.0001 0.0002 0.0003
qq =
-0.0035 0.0001 0.0005 -0.0696
uu =
1.0e-03 *
-0.0025
-0.1127
0.0003
-0.0276
0.0000
-0.0226
-0.0071
k =
3
delt =
-0.0475 -0.0896 -0.0959 -0.1097 0
u =
1.0296 1.0005 0.9998 1.0000 1.0500
pp =
1.0e-04 *
-0.3598 0.0441 0.0600 0.0992
qq =
-0.0001 0.0000 0.0000 -0.0705
uu =
1.0e-05 *
-0.0003
-0.3269
-0.0003
-0.0022
-0.0004
-0.0003
-0.0008
k =
4
delt =
-0.0475 -0.0896 -0.0959 -0.1097 0
u =
1.0296 1.0005 0.9998 1.0000 1.0500
Q4 =
-0.0294
Y =
10.2000 -31.5000i -1.2000 + 4.0000i -1.5000 + 5.0000i -2.5000 + 7.5000i -5.0000 +15.0000i
-1.2000 + 4.0000i 10.4000 -31.7000i -8.0000 +24.0000i 0.0000 + 0.0000i -1.2000 + 3.7000i
-1.5000 + 5.0000i -8.0000 +24.0000i 10.7000 -32.7000i -1.2000 + 3.7000i 0.0000 + 0.0000i
-2.5000 + 7.5000i 0.0000 + 0.0000i -1.2000 + 3.7000i 3.7000 -11.2000i 0.0000 + 0.0000i
-5.0000 +15.0000i -1.2000 + 3.7000i 0.0000 + 0.0000i 0.0000 + 0.0000i 6.2000 -18.7000i
JJ =
-33.1934 -11.0132 4.1688 1.0618 5.2158 1.2931 7.8673
10.6131 -33.5937 -1.0618 4.1688 -1.2931 5.2158 -2.0888
4.0649 1.4083 -31.8813 -9.9603 24.0578 7.8488 0
-1.4083 4.0649 10.8603 -31.5813 -7.8488 24.0578 0
5.0663 1.7916 23.9556 8.1557 -32.7373 -10.2958 3.7155
-1.7916 5.0663 -8.1557 23.9556 11.0958 -32.6373 -1.1486
7.5471 3.0493 0 0 3.6824 1.2507 -11.2295
S =
0.0000 + 0.0000i 0.2103 + 0.0716i 0.2971 + 0.0847i 0.5615 + 0.0835i -0.8689 - 0.0399i
-0.2071 - 0.0609i 0.0000 + 0.0000i 0.1591 - 0.0341i 0.0000 + 0.0000i -0.4020 - 0.0549i
-0.2921 - 0.0682i -0.1588 + 0.0351i 0.0000 + 0.0000i 0.0509 - 0.0169i 0.0000 + 0.0000i
-0.5493 - 0.0471i 0.0000 + 0.0000i -0.0507 + 0.0176i 0.0000 + 0.0000i 0.0000 + 0.0000i
0.8831 + 0.0827i 0.4151 + 0.0952i 0.0000 + 0.0000i 0.0000 + 0.0000i 0.0000 + 0.0000i
B =
-31.5000 4.0000 5.0000 7.5000 15.0000
4.0000 -31.7000 24.0000 0 3.7000
5.0000 24.0000 -32.7000 3.7000 0
7.5000 0 3.7000 -11.2000 0
15.0000 3.7000 0 0 -18.7000
pp =
1.0e-04 *
-0.3598 0.0441 0.0600 0.0992
qq =
-0.0001 0.0000 0.0000 -0.0705
uu =
1.0e-05 *
-0.0003
-0.3269
-0.0003
-0.0022
-0.0004
-0.0003
-0.0008
U =
1.0285 - 0.0489i 0.9965 - 0.0895i 0.9952 - 0.0958i 0.9940 - 0.1095i 1.0500 + 0.0000i
k =
4
Q4 =
-0.0294
S5 =
1.2982 + 0.1779i
>> 六、报告要求
(1) 画出程序流程图；
(2) 给出雅克比矩阵参数求解、不平衡量求解、各条线路功率求解的关键代码；
(3) 给出迭代过程中各节点电压的计算值；
(4) 给出各条线路功率的计算结果。
七、报告
《电力系统基础》实验报告
姓名： 学号： 日期： 成绩：
实验名称：复杂电力系统的潮流计算编程
实验学时：6学时
实验内容：
1. 牛顿-拉夫逊法求非线性方程组近似解范例；
2. 采用牛顿-拉夫逊法实现电力系统潮流计算编程。
实验要求：
1. 通过MATLAB编程范例，理解实现牛顿-拉夫逊法的基本代码编写方法；
2. 根据给定的电力系统网络接线图、节点类型和具体参数，运用以极坐标形式的牛顿-拉夫逊法计算系统的潮流分布。
范例程序求解结果：迭代次数k = 6 ，结果x1 = 1.3478 ，x2 = 1.5045 。
潮流计算程序流程图：
部分关键代码：
见上
潮流计算结果：
实验结果分析：
(1)本次课程设计让我从真正认识了牛顿拉夫逊计算潮流的方法,通过自己编程和同学一起讨论我会使用牛拉法求潮流了。首先最重要的就是建立节点导纳短阵,对于变压器和线路导纳的处理。会用编程语言实现。
(2)在编程以前要做好前提工作,分析题目,画好等值电路,求出各节点的导纳阻抗,还要分析好 PQ , PV 节点,没立平衡节点。身将平衡节点编号放到最大。
(3)从本次实例可以看出,牛拉法收敛速度快。结果精确(误差< 0.00001)。经过六次迭代就已经收敛。