一、实验学时：

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。

求解代码示例：

clearx(1)=1.0;x(2)=2.0;k=0; precision=1;k,xwhile 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));end

2. 采用牛顿-拉夫逊法实现电力系统潮流计算编程

网络接线如图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) 节点参数和导纳矩阵相关代码：

clearG(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,qqfor 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));uufor 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 endY JJ S B 	pp qq uu U k Q4 S5 

运行结果：

>> dianli2k =     0delt =     0     0     0     0u =     1     1     1     1pp =    0.4500   -0.3900   -0.4000   -0.6000qq =    0.9500    0.0350   -0.0500   -0.1000uu =   -0.0476    0.0329   -0.0903    0.0041   -0.0967    0.0032   -0.1097k =     1delt =   -0.0476   -0.0903   -0.0967   -0.1097         0u =    1.0329    1.0041    1.0032    1.0000    1.0500pp =   -0.0306   -0.0002    0.0078    0.0101qq =   -0.0751   -0.0217   -0.0095   -0.0326uu =    0.0001   -0.0032    0.0008   -0.0036    0.0008   -0.0033    0.0000k =     2delt =   -0.0475   -0.0896   -0.0959   -0.1097         0u =    1.0297    1.0005    0.9998    1.0000    1.0500pp =   -0.0012    0.0001    0.0002    0.0003qq =   -0.0035    0.0001    0.0005   -0.0696uu =   1.0e-03 *   -0.0025   -0.1127    0.0003   -0.0276    0.0000   -0.0226   -0.0071k =     3delt =   -0.0475   -0.0896   -0.0959   -0.1097         0u =    1.0296    1.0005    0.9998    1.0000    1.0500pp =   1.0e-04 *   -0.3598    0.0441    0.0600    0.0992qq =   -0.0001    0.0000    0.0000   -0.0705uu =   1.0e-05 *   -0.0003   -0.3269   -0.0003   -0.0022   -0.0004   -0.0003   -0.0008k =     4delt =   -0.0475   -0.0896   -0.0959   -0.1097         0u =    1.0296    1.0005    0.9998    1.0000    1.0500Q4 =   -0.0294Y =  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.7000iJJ =  -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.2295S =   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.0000iB =  -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.7000pp =   1.0e-04 *   -0.3598    0.0441    0.0600    0.0992qq =   -0.0001    0.0000    0.0000   -0.0705uu =   1.0e-05 *   -0.0003   -0.3269   -0.0003   -0.0022   -0.0004   -0.0003   -0.0008U =   1.0285 - 0.0489i   0.9965 - 0.0895i   0.9952 - 0.0958i   0.9940 - 0.1095i   1.0500 + 0.0000ik =     4Q4 =   -0.0294S5 =   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)。经过六次迭代就已经收敛。