电力系统分析课程设计——牛顿拉夫逊潮流计算.doc

武汉理工大学电力系统分析课程设计说明书 1 电力系统编程潮流计算 1 设计任务及初步分析 1。1 设计任务 条件：节点数:3 支路数:3 计算精度:0。00010 支路 1:0.0300+j0.0900 1-2 支路 2:0。0200+j0。0900 2-3 支路 3:0。0300+j0。0900 3-1 节点 1:PQ 节点,S（1）=0.5000j0.2000 节点 2：PQ 节点,S（2）=-0.6000j0。2500 节点 3：平衡节点，U（3）=1。00000。0000 要求:编写程序计算潮流 1。2 初步分析 潮流计算在数学上可归结为求解非线性方程组，其数学模型简写如下：2 牛顿-拉夫逊法简介 2。1 概述 牛顿拉夫逊法是目前求解非线性方程最好的一种方法。这种方法的特点就是把对非线性方程的求解过程变成反复对相应的线性方程求解的过程,通常称为逐次线性化过程，就是牛顿拉夫逊法的核心。牛顿-拉夫逊法的基本原理是在解的某一邻域内的某一初始点出发，沿着该点的一阶偏导数雅可比矩阵，朝减小方程的残差的方向前进一步，在新的点上再计算残差和雅可矩阵继续前进，重复这一过程直到残差达到收敛标准，即得到了非线性方程组的解.因为越靠近解,偏导数的方向越准，收敛速度也越快，所以牛顿法具有二阶收敛特性.而所谓“某一邻域是指雅可比方向均指向解的范围，否则可能走向非线性函数的其它极值点，一般来说潮流由平电压即各母线电压（相角为0,幅值为1)启动即在此邻域内.22 一般概念 对于非线性代数方程组 即 （21)武汉理工大学电力系统分析课程设计说明书 2 在待求量的某一个初始计算值附件,将上式展开泰勒级数并略去二阶及以上的高阶项，得到如下的线性化的方程组 （22)上式称之为牛顿法的修正方程式。由此可以求得第一次迭代的修正量 （23)将和相加，得到变量的第一次改进值。接着再从出发，重复上述计算过程。因此从一定的初值出发，应用牛顿法求解的迭代格式为 (24)(25）上两式中:是函数对于变量的一阶偏导数矩阵，即雅可比矩阵；为迭代次数.由式（24)和式子（25）可见，牛顿法的核心便是反复形成求解修正方程式。牛顿法当初始估计值和方程的精确解足够接近时，收敛速度非常快,具有平方收敛特性.2。3 潮流计算的修正方程 运用牛顿拉夫逊法计算潮流分布时,首先要找出描述电力系统的非线性方程.这里仍从节点电压方程入手，设电力系统导纳矩阵已知，则系统中某节点（节点)电压方程为 从而得 进而有 （26)式（26）中，左边第一项为给定的节点注入功率，第二项为由节点电压求得的节点注入功率。他们二者之差就是节点功率的不平衡量.现在有待解决的问题就是各节点功率的不平衡量都趋近于零时，各节点电压应具有的价值。由此可见，如将式（26）作为牛顿拉夫逊中的非线性函数,其中节点电压就相当于变量。建立了这种对应关系，就可列出修正方程式，并迭代求解。但由于节点电压可有两种表示方式-以直角做表或者极坐标表示，因而列出的迭代方程相应地也有两种,下面分别讨论。2。4 直角坐标表示的修正方程 节点电压以直角坐标表示时,令、,且将导纳矩阵中元素表示为，则式（27）改变为 （27）再将实部和虚部分开，可得 (28）这就是直角坐标下的功率方程。可见,一个节点列出了有功和无功两个方程.对于节点（），给定量为节点注入功率，记为、,则由式（28）可得功率的不平衡量，作为非线性方程 （29)武汉理工大学电力系统分析课程设计说明书 3 式中、-分别表示第节点的有功功率的不平衡量和无功功率的不平衡量。对于节点(），给定量为节点注入有功功率及电压数值，记为、，因此,可以利用有功功率的不平衡量和电压的不平衡量表示出非线性方程，即有 (210）式中为电压的不平衡量.对于平衡节点(),因为电压数值及相位角给定，所以也确定，不需要参加迭代求节点电压。因此,对于个节点的系统只能列出个方程,其中有功功率方程个，无功功率方程个，电压方程个.将式(29)、式（210）非线性方程联立,称为个节点系统的非线性方程组，且按泰勒级数在、（）展开，并略去高次项，得到以矩阵形式表示的修正方程如下.（211）上式中雅可比矩阵的各个元素则分别为 将（211）写成缩写形式 （212)对雅可比矩阵各元素可做如下讨论：当时，对于特定的，只有该特定点的和是变量，于是雅可比矩阵中各非对角元素表示为 当时，雅可比矩阵中各对角元素的表示式为 由上述表达式可知，直角坐标的雅可比矩阵有以下特点：1)雅可比矩阵是阶方阵,由于、等等，所以它是一个不对称的方阵。2)雅可比矩阵中诸元素是节点电压的函数，在迭代过程中随电压的变化而不断地改变.3)雅可比矩阵的非对角元素与节点导纳矩阵中对应的非对角元素有关，当中的为零时，雅可比矩阵中相应的、也都为零，因此，雅可比矩阵也是一个稀疏矩阵。3 程序设计 3.1 程序流程图 图 31 程序流程图 3。2 源程序 n=3 n1=3 isb=3 pr=0。0001 B1=1 2 0。03+0.09i 0 1 0；1 3 0.03+0。09i 0 1 0;2 3 0。02+0。09i 0 1 0 B2=0-0。5-0.2i 1 0 0 2;0-0.6-0。25i 1 0 0 2；0 0 1 1 0 1 X=1 0;2 0;3 0 X=input(节点号和对地参数:X=);Y=zeros（n）;Times=1；%置迭代次数为初始值%创建节点导纳矩阵 武汉理工大学电力系统分析课程设计说明书 4 for i=1:n1 if B1（i,6）=0 不含变压器的支路 p=B1(i,1);q=B1（i，2）;Y(p，q)=Y(p，q)-1/B1(i,3);Y(q,p）=Y（p,q）;Y(p,p)=Y(p,p)+1/B1（i,3)+0。5B1(i,4）；Y(q,q）=Y(q，q)+1/B1(i，3)+0.5B1(i,4）；else 含有变压器的支路 p=B1(i,1）；q=B1（i,2）；Y（p，q)=Y（p，q)1/(B1(i,3）B1（i，5)；Y(q，p)=Y（p,q）；Y(p，p）=Y（p，p）+1/B1(i,3)；Y（q，q)=Y(q,q）+1/（B1(i，5)2B1（i,3);end end Y OrgS=zeros(2n2，1）;DetaS=zeros(2*n2，1）；将 OrgS、DetaS 初始化 创建 OrgS，用于存储初始功率参数 h=0；j=0；for i=1:n%对 PQ 节点的处理 if i=isbB2(i,6)=2 h=h+1；for j=1：n OrgS（2h1，1）=OrgS(2h1,1)+real(B2(i,3)(real（Y(i，j)）real（B2(j，3）)imag(Y（i，j）imag(B2(j，3）+imag（B2（i，3）*（real（Y（i，j)）imag（B2（j,3）+imag（Y（i,j)）*real(B2（j,3）)；OrgS(2*h，1)=OrgS(2*h,1）+imag（B2（i,3）*（real(Y（i,j)）*real（B2(j,3）-imag(Y(i，j)）*imag（B2(j,3)）)real（B2(i,3)(real(Y(i，j）)imag（B2(j,3)）+imag(Y(i，j)）real（B2（j,3);end end end for i=1:n 对 PV 节点的处理，注意这时不可再将 h 初始化为 0 if i=isbB2（i，6)=3 h=h+1；for j=1:n OrgS（2h1,1）=OrgS（2h1,1)+real（B2（i，3）(real(Y（i，j)）*real(B2（j,3)）imag(Y（i,j)*imag(B2（j,3）)+imag（B2(i,3)*（real（Y(i,j)imag（B2(j，3)）+imag(Y(i，j）real(B2（j,3）)；OrgS（2*h,1)=OrgS(2h，1)+imag（B2(i，3)）*(real(Y（i,j）*real（B2(j，3)）武汉理工大学电力系统分析课程设计说明书 5 -imag（Y(i，j）)imag(B2(j,3）)-real（B2(i,3）*(real（Y（i，j）*imag（B2（j，3）+imag（Y(i,j)）real（B2(j,3）)；end end end OrgS 创建 PVU 用于存储 PV 节点的初始电压 PVU=zeros（nh-1,1)；t=0；for i=1:n if B2(i，6）=3 t=t+1;PVU(t，1)=B2(i,3）;end end PVU 创建 DetaS，用于存储有功功率、无功功率和电压幅值的不平衡量 h=0;for i=1：n 对 PQ 节点的处理 if i=isbB2（i,6）=2 h=h+1；DetaS（2h-1，1)=real(B2（i，2)）OrgS(2*h1,1）；DetaS(2*h，1)=imag(B2(i,2）OrgS（2h，1);end end t=0；for i=1:n 对 PV 节点的处理,注意这时不可再将 h 初始化为 0 if i=isb&B2(i，6)=3 h=h+1；t=t+1；DetaS(2h1，1）=real（B2(i,2）OrgS(2*h1，1)；DetaS（2h,1）=real（PVU(t，1）2+imag（PVU（t,1)2real（B2(i，3)）2imag(B2(i,3）2；end end DetaS%创建 I，用于存储节点电流参数 i=zeros(n-1,1）；h=0;for i=1:n if i=isb h=h+1；I(h,1)=（OrgS(2*h1,1）-OrgS(2*h，1)*sqrt(1）/conj（B2（i,3）；end 武汉理工大学电力系统分析课程设计说明书 6 end I%创建 Jacbi(雅可比矩阵）Jacbi=zeros(2n2）;h=0；k=0；for i=1：n 对 PQ 节点的处理 if B2（i，6)=2 h=h+1;for j=1:n if j=isb k=k+1；if i=j%对角元素的处理 Jacbi（2*h-1,2k-1)=imag(Y（i，j）real（B2（i,3）+real（Y（i,j）imag(B2（i,3）)+imag（I（h，1）)；Jacbi（2*h-1，2*k）=real(Y（i，j)real(B2(i，3）+imag(Y(i，j)）imag(B2(i，3）+real(I（h,1）);Jacbi(2*h,2*k1)=Jacbi(2h1，2*k)+2real（I（h,1）；Jacbi（2h，2k）=Jacbi(2*h1，2k-1）-2*imag(I(h，1）);else 非对角元素的处理 Jacbi(2*h-1,2k-1）=-imag（Y(i,j)）*real（B2（i，3)）+real(Y（i,j）)imag(B2（i,3);Jacbi(2h1,2k）=real（Y（i,j)real（B2(i，3)）+imag（Y(i,j）)imag（B2（i,3);Jacbi(2h,2k-1）=-Jacbi(2*h-1，2*k);Jacbi(2*h，2k）=Jacbi（2*h-1,2*k-1）;end if k=(n-1)将用于内循环的指针置于初始值，以确保雅可比矩阵换行 k=0;end end end end end k=0;for i=1：n 对 PV 节点的处理 if B2(i,6）=3 h=h+1；for j=1:n if j=isb k=k+1；if i=j%对角元素的处理 Jacbi（2h1，2*k-1)=-imag（Y(i，j）)*real（B2(i,3）)+real(Y(i,j）imag(B2(i,3）)+imag(I（h，1）；武汉理工大学电力系统分析课程设计说明书 7 Jacbi(2h1,2k）=real(Y（i，j）*real（B2（i，3)+imag（Y(i，j)imag（B2（i,3）)+real(I（h，1)）;Jacbi（2h,2*k-1)=2*imag(B2(i，3）);Jacbi(2h,2k)=2real（B2(i,3）；else%非对角元素的处理 Jacbi（2*h1,2*k1)=imag(Y（i，j）)real（B2(i,3)）+real(Y(i,j)imag(B2（i，3）;Jacbi（2*h1,2k)=real（Y(i，j）)real(B2（i,3）)+imag(Y（i,j）)imag(B2(i，3）)；Jacbi(2*h，2k1)=0；Jacbi(2h,2k）=0；end if k=（n1)将用于内循环的指针置于初始值，以确保雅可比矩阵换行 k=0;end end end end end Jacbi 求解修正方程，获取节点电压的不平衡量 DetaU=zeros（2*n2，1）；DetaU=inv(Jacbi)DetaS;DetaU%修正节点电压 j=0；for i=1:n%对 PQ 节点处理 if B2（i，6)=2 j=j+1;B2(i，3)=B2(i，3)+DetaU（2j，1)+DetaU（2j-1,1）*sqrt(-1)；end end for i=1:n 对 PV 节点的处理 if B2（i，6)=3 j=j+1;B2（i，3）=B2（i,3）+DetaU(2j,1）+DetaU（2j-1，1）*sqrt(1）；end end B2 while abs（max(DetaS)）pr OrgS=zeros(2*n-2,1）;！!初始功率参数在迭代过程中是不累加的，所以在这里必须将其初始化为零矩阵 h=0；j=0；武汉理工大学电力系统分析课程设计说明书 8 for i=1：n if i=isbB2（i，6）=2 h=h+1；for j=1：n OrgS(2*h1,1）=OrgS（2h1,1)+real（B2（i，3)）(real(Y（i，j)real（B2(j,3)imag（Y（i,j)）*imag(B2(j，3）)+imag(B2(i，3）)*（real（Y（i，j)）imag(B2(j,3)+imag（Y(i,j)real(B2（j，3)）;OrgS(2*h,1）=OrgS(2*h，1）+imag（B2(i，3）*（real(Y（i,j)*real（B2（j,3）)imag(Y(i，j)imag(B2(j，3)）)real(B2（i,3）)*(real（Y(i,j)）imag（B2(j,3)+imag（Y（i，j）*real（B2（j，3)；end end end for i=1：n if i=isb&B2(i,6)=3 h=h+1；for j=1:n OrgS（2*h-1,1)=OrgS(2h-1，1）+real（B2(i,3)）(real(Y(i,j)real（B2（j，3）-imag(Y(i，j)）imag(B2(j,3)）)+imag（B2（i,3)）*(real（Y（i,j)*imag（B2(j,3)）+imag(Y（i，j)）*real(B2(j，3）)；OrgS（2h，1)=OrgS(2*h,1）+imag(B2(i,3)*(real（Y(i,j)real（B2(j，3）-imag（Y（i，j)）imag（B2(j，3）)）real（B2(i,3)(real（Y（i，j）*imag(B2(j,3）)+imag（Y(i,j）*real（B2（j，3）；end end end OrgS 创建 DetaS h=0;for i=1：n if i=isb&B2（i，6)=2 h=h+1;DetaS（2h1,1）=real（B2(i,2）)OrgS(2h-1，1)；DetaS(2h,1）=imag(B2(i，2)）OrgS（2h,1）；end end t=0；for i=1：n if i=isb&B2(i,6)=3 h=h+1；t=t+1；DetaS(2h1，1）=real（B2(i,2）-OrgS(2h1,1);DetaS(2*h，1)=real（PVU（t,1）)2+imag（PVU(t，1)2-real（B2(i,3）2-imag(B2(i,3）2；武汉理工大学电力系统分析课程设计说明书 9 end end DetaS%创建 I i=zeros(n-1，1）；h=0;for i=1:n if i=isb h=h+1;I（h，1)=（OrgS（2h1，1）OrgS(2h,1）*sqrt(-1）)/conj(B2（i,3）;end end I 创建 Jacbi Jacbi=zeros(2n-2）;h=0;k=0;for i=1:n if B2（i,6)=2 h=h+1;for j=1：n if j=isb k=k+1;if i=j Jacbi(2*h-1，2k1）=imag（Y（i,j）real（B2(i，3)）+real(Y（i,j）)imag(B2(i,3）)+imag（I（h,1)；Jacbi(2 h1,2*k）=real（Y(i，j）real(B2（i,3）+imag(Y(i,j）)*imag(B2(i,3）)+real(I(h,1）);Jacbi（2h,2k-1)=Jacbi（2*h1，2*k)+2real（I(h,1)；Jacbi(2h,2k）=Jacbi(2h-1，2*k1）-2*imag(I(h,1）)；else Jacbi（2h1，2k-1)=-imag（Y(i,j）)*real(B2(i,3)）+real（Y(i，j）)*imag（B2（i，3)）;Jacbi(2*h1，2k）=real(Y（i,j)*real（B2(i，3）+imag（Y(i，j)*imag（B2（i,3)；Jacbi（2h,2k1）=Jacbi（2*h1,2k)；Jacbi(2h，2*k)=Jacbi(2*h1,2*k1）;end if k=(n1)k=0;end end end end 武汉理工大学电力系统分析课程设计说明书 10 end k=0；for i=1:n if B2（i,6）=3 h=h+1;for j=1:n if j=isb k=k+1;if i=j Jacbi（2h-1，2*k-1）=-imag（Y（i,j）*real(B2（i，3)+real(Y(i，j）)imag（B2(i，3)）+imag(I(h,1)；Jacbi（2*h1，2*k）=real(Y(i，j)）*real（B2（i，3）+imag(Y（i，j)）*imag(B2（i,3)）+real(I(h，1)；Jacbi(2*h，2k-1）=2*imag(B2（i，3）;Jacbi(2h，2*k)=2*real（B2(i,3)）；else Jacbi(2h-1,2*k1)=imag（Y（i,j)real(B2（i，3）+real(Y(i,j）)imag（B2(i,3）)；Jacbi(2h1,2*k）=real(Y(i,j）real(B2(i，3)）+imag（Y(i，j)*imag(B2（i，3)；Jacbi（2*h,2*k-1)=0;Jacbi(2*h,2k)=0；end if k=(n-1)k=0;end end end end end Jacbi DetaU=zeros(2*n2,1);DetaU=inv（Jacbi)DetaS；DetaU 修正节点电压 j=0；for i=1：n if B2（i，6)=2 j=j+1;B2(i，3)=B2（i，3)+DetaU（2*j,1）+DetaU（2j1，1）*sqrt（-1）;end end for i=1：n if B2(i，6)=3 武汉理工大学电力系统分析课程设计说明书 11 j=j+1；B2（i，3)=B2（i，3）+DetaU（2*j，1）+DetaU（2j1,1）sqrt（1）;end end B2 Times=Times+1；迭代次数加 1 end disp（迭代次数为);disp(Times)for no=1:n Vn(no)=B2（no,3)；str1=sprintf（节点d 的电压是d,no)；disp(str)；disp（Vn（no);end 3。3 结果及分析 节点导纳矩阵为 迭代运行了 3 次，每次功率的不平衡量表 31 所示 表 31 迭代过程中各节点的功率不平衡量 k 0-0。5-j0。2-0.6j0.25 1-0。0083-j0。0279-0.0091-j0。0360 2-0。0000297-j0.0000955-0。0000326j0。0001332 3 0.0000+j0。0000 0。0000+j0.0000 表 32 迭代过程中各节点电压 k 1 0。9664 j0。0423 0。9658 j0。0457 2 0.9632 j0。0423 0。9623-j0.0457 3 0。9632-j0.0423 0。9623-j0.0457 平衡节点功率为 下面计算网络中的功率分布，根据所有线路功率计算如下:武汉理工大学电力系统分析课程设计说明书 12 同样可以得出 =0.5646 j0.2543 因此各支路功率损耗为 S12=4.1166e005+j1.2350e004 S13=0。0105+j0.0315 S23=0.0083+j0.0372 4 小结 此次课程设计学习了不少复杂潮流计算的方法，包括牛顿拉夫逊法、PQ分解法等等，合理运用这些方法使得复杂网络的潮流分布计算有迹可循，相对而言也更为简便。尤其这些方法，可以很好地通过程序来模拟，将传统的算法变为计算机算法.这次课程设计要求编程计算潮流，经过多天的设计，我明白了以下道理：（1）要培养多学科关联的意识,构建自己的理论.电力系统分析涉及到很多的其他学科知识,尤其是数学，地位尤为突出。电力系统分析最关键地要把实际模型转换为实际模型，因此数学功底好，电力系统分析也就游刃有余.此次课程设计的题目需要将其写成非线性方程组,而牛顿拉夫逊法正好可以解非线性方程组，因此建模时电力系统分析与数学之间沟通的桥梁。（2）电力系统与软件相结合，可以大大节省人工工作量,电力系统是个有机整体，相对而言比较复杂，各种潮流计算对于人工而言不仅仅是大量重复的计算分析，更是导致出错的重要原因，因此在要求比较高的场合下往往都不仅仅是靠人工去直接控制操作的,而是借助计算机软件去控制。本次设计可以明显看出程序大大简化的计算，减小了工作量。（3）本次通过 matlab 来编写程序。Matlab 是一款非常强大的软件,涵盖了多门学科，包括数学、电路等学科，用它来做仿真分析，可以设计比较合理稳定的程序、电路等，因此掌握好软件的用法显得尤为重要，不过电力系统分析的软件不止这一款.（4）编写程序是一件考研耐心与细心的事情,稍有不慎，可能就会出现意想不到的结果，甚至不能出现结果，因此这是一件非常头疼的事情，经常因为一个小小的 bug 要去反复查看整个程序,因此细心一点可以大大提高工作效率.（5）电力系统分析是一门与电力息息相关的一门课，学好它并且运用好它显得尤为重要.在今后地学习中，我还必须不断更新自己的知识，坚持终身学习。5 参考文献 1 何仰赞，温增银.电力系统分析.武汉：华中科技大学出版社，2006。武汉理工大学电力系统分析课程设计说明书 13 2 陈珩.电力系统稳态分析（第三版).北京：中国电力出版社，2007。3 于永源，杨绮雯。电力系统分析（第二版）。北京:中国电力出版社，2004。