matlab流体力学计算代码：基于MATLAB牛顿-拉夫逊算法潮流计算算例

一、实验学时：

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. 根据给定的电力系统网络接线图、节点类型和具体参数，运用以极坐标形式的牛顿-拉夫逊法计算系统的潮流分布。

求解过程：

令

迭代次数为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));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) 节点参数和导纳矩阵相关代码：

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)。经过六次迭代就已经收敛。