电力系统稳态潮流计算课程设计

1、课程设计(论文题目名称电力系统潮流计算课程名称电力系统稳态分析学生姓名学号系、专业电气工程系电气工程及其自动化电力方向指导教师2009年1月6日前言在如今的社会,电力已经成为人们必不可少的需求,而建立结构合理的大型电力系统不仅便于电能生产与消费的集中管理、统一调度和分配,减少总装机容量,节省动力设施投资,且有利于地区能源资源的合理开发利用,更大限度地满足地区国民经济日益增长的用电需要。电力系统建设往往是国家及地区国民经济发展规划的重要组成部分。电力系统的出现,使高效、无污染、使用方便、易于调控的电能得到广泛应用,推动了社会生产各个领域的变化,开创了电力时代,发生了第二次技术革命。电力系统的规模

2、和技术水准已成为一个国家经济发展水平的标志之一。电力系统稳态分析包括潮流计算(或潮流分析和静态安全分析。潮流计算针对电力革统各正常运行方式,而静态安全分析则要研究各种运行方式下个别系统元件退出运行后系统的状况。其目的是校验系统是否能安全运行,即是否有过负荷的元件或电压过低的母线等。原则上讲,静态安全分析也可U用潮流计算来代替。但是一般静态安全分析需要校验的状态数非常多,用严格的潮流计算来分析这些状态往往计算量过大,因此不得不寻求一些特殊的算法以满足要求。牛顿法是数学中解决非线性方程式的典型方法,有较好的收敛性。解决电力系统潮流计算问题是以导纳距阵为基础的,因此,只要在迭代过程中尽可能保持方程式

3、系数距阵的稀疏性,就可以大大提高牛顿法潮流程序的放率。自从20 世纪60 年代中期利用了最佳顺序消去法以后,牛顿法在收敛性、内存要求、速度方面都超过了阻抗法,成为直到目前仍在广泛采用的优秀方法。目录第一章系统概述 (41.1 设计目的与要求 (41.2 设计题目 (41.3 设计内容 (4第二章潮流计算设计题目 (52.1 潮流计算题目 (52.2 对课题的分析及求解思路 (6第三章潮流计算算法及手工计算 (63.1 潮流计算算法 (63.2 关于电力系统潮流计算手工计算 (8第四章Matlab概述 (104.1 Matlab简介 (104.2 矩阵的运算 (11第五章潮流计算流程图及源程序

4、(135.1 潮流计算流程图 (135.2 潮流计算源程序 (145.3 运行计算结果 (19总结 (20参考文献 (21第一章系统概述1.1 设计目的与要求1. 掌握电力系统潮流计算的基本原理;2. 掌握并能熟练运用一门计算机语言(MATLAB语言或FORTRAN或C语言或C+语言;3. 采用计算机语言对潮流计算进行计算机编程计算。1. 程序源代码;2. 给定题目的输入,输出文件;3. 程序说明;4. 给定系统的程序计算过程;5. 给定系统的手算过程(至少迭代2次。1.2 设计题目电力系统潮流计算(牛顿-拉夫逊法、P-Q 分解法1.3 设计内容1.根据电力系统网络推导电力网络数学模型,写出节

5、点导纳矩阵;2.赋予各节点电压变量(直角坐标系形式初值后,求解不平衡量;3.形成雅可比矩阵;4.求解修正量后,重新修改初值,从2开始重新循环计算;5.求解的电压变量达到所要求的精度时,再计算各支路功率分布、功率损耗和平衡节点功率;6.上机编程调试;连调;7.计算分析给定系统潮流分析并与手工计算结果作比较分析。8.准备计算机演示答辩,书写该课程设计说明书(必须计算机打印。第二章潮流计算设计题目2.1 潮流计算题目 图2-1 电力系统接线图2.2 对课题的分析及求解思路此电力系统是一个5节点,4支路的电力网络。其中包含3个PQ 节点,一个PV 节点,和一个平衡节点。综合比较牛顿拉夫逊法(直角坐标、

6、极坐标、PQ 分解法等多种求解方法的特点,最后确定采用牛顿拉夫逊法(极坐标。因为此方法所需解的方程组最少。第三章 潮流计算算法及手工计算3.1 潮流计算算法本题采用了题目要求的牛顿-拉夫逊潮流计算的方法。牛顿-拉夫逊法潮流计算的公式。把牛顿法用于潮流计算,采用直角坐标形式表示的如式(1-3所示的形式。其中电压和支路导纳可表示为:i i iij ij ijj j j ij ij ijU e jf Y G jB U e jf Y G jB *=+=+=-=- (1-2将上述表示式(1-2代入(1-1式的右端,展开并分出实部和虚部,便得:1111(n ni i ij j ij j i i j j i

7、j j j i nni i ij j ij j i i j j ij j j j P e G e B f f G f B e Q f G e B f e G f B e =-+=-+ (1-3按照以上的分类,PQ 节点的输出有功功率和无功功率是给定的,则第i 节点的给定功率设为is P 和is Q (称为注入功率。假定系统中的第1、2、m 节点为PQ 节点,对其中每一个节点的N-R 法表达式F(x=0如0i S =、0i P =、0i Q =形式有些下列方程:1111(0(0n ni is i is i ij j ij j i ij j ij j j j n ni is i is i ij j

8、 ij j i ij j ij j j j P P P P e G e B f f G f B e Q Q Q Q f G e B f e G f B e =-=-+=-=-+=(1-4i =(1、2、m PV 节点的有功功率和节点电压幅值是给定的。假定系统中的第m+1、m+2、n-1节点为PV 节点,则对其中每一PV 节点可以列写方程:11222222(0(n ni is i is i ij j ij j i ij j ij j j j is i is i i P P P P e G e B f f G f B e U U U U e f =-=-+=-=-+ (1-5 i =(m+1、m+

9、2、n-1(6形成雅可比矩阵。N-R 法的思想是(0F x F x x +=;本例(P j Q F x +=;对F(x求偏导的式(1-6、式(1-7,即式(1-4、式(1-5中的0i P =、0i Q =、U 是多维变量的函数,对多维变量求偏导(i i P e 、i j P e 、i i P f 、i j P f 、i i Q e 、ijQ e 、iiP e 、,并以矩阵的形式表达称为雅可比矩阵。当j=i 时,对角元素为111122(22ni ij j ij j ii i ii i ii j i n iij j ij j ii i ii i ii j i n iij j ij j ii i i

10、i i ii j i niij j ij j ii i ii ii ii j i i iiii i P G e B f G e B f N e P G f B e B e G f H f Q G f B e B e G f L e Q G e B f G e B f J f U e e U f f =-=-+-=-+-=-+=-=-(1-6当j i 时,矩阵非对角元素为:22(0i iij i ij i ij ij i j i i ij i ij i ij ij j ji i j j P Q G e B f N J e f P Q B e G f H L f e U U e f =-=-+=-

11、=-=(1-7由上式不难看出,雅可比矩阵有以下特点。 雅可比矩阵中的诸元素都是节点电压的函数,因此在迭代过程中,它们将随着节点电压的变化而不断的变化。 雅可比矩阵具有结构对称性,数据不对称。如非对角i j j i H H ,ij ij i ij i H B e G f =-,ji ij j ij j H B e G f =-。 由式(1-7可以看出,当导纳矩阵中非对角元素ij Y 为零时,。雅可比矩阵中相应的元素也为零,即矩阵是非常稀疏的。因此,修正方程的求解同样可以应用稀疏矩阵的求解技巧。正是由于这一点才使N-R 法获得广泛的应用。3.2 关于电力系统潮流计算手工计算Y ii =ij y j

12、i y +0y Y ikik -=各节点的导纳值如下:Y j=-+;4200Y j=+;Y j=+-;431.250 3.750Y j=-+;341.250 3.750Y j=-+;Y j=+-;3500Y j=+;4500Y j=+;Y j=-+;521.250 3.750Y j=-+;5300Y j=+;5400Y j=+;Y j=+-.形成有功迭代和无功迭代的简化雅可比矩阵B/和B/B/= B/= 将B/ 和B/进行三角分解: 给定PQ 节点初值和各节点电压相角初值 V 1=1.050。 ,V 2(0=V 3(0=1.0,V 4=1.1 2(0=3(0=0, 4(0=01 作第一次有功迭

13、代,按公式计算节点有功功率不平衡量 P 2(0=-0.55-(-0.024037=-0.525963 P 3(0=-0.30-(-0.022695=-0.277305 P 4(0=0.500000P 1(0/V 1(0=0.454545 P 2(0/ V 2(0=-0.525963 P 3(0/V 3(0=-0.2773092做第一次无功迭代,按公式计算无功功率不平衡量,计算时电压相角最新的修正值。Q 2(0=-0.13-(-0.001550=-0.039594Q 3(0=-0.18-(-0.14406=-0.039588 Q 2(0/ V 2(0=-0.131553 Q 3(0/V 3(0=

14、-0.039588解修正方程式,可得各节点电压幅值的修正量为035224.00(2-=V V 3(0=-0.014855 于是有:V 2(1 = V 2(0+V 2(1=0.964776 V 3(1 = V 3(0+V 3(1=0.985145 到这里为止,第一轮有功迭代和无功迭代便做完了。 3 按公式计算平衡节点功率,得:P 1+jQ 1=0.367885+j0.264696经过四轮迭代,节点不平衡功率也下降到10-5以下,迭代到此结束。第四章 Matlab 概述4.1 Matlab 简介目前电子计算机已广泛应用于电力系统的分析计算,潮流计算是其基本应用软件之一。现有很多潮流计算方法。对潮流

15、计算方法有五方面的要求:(1计算速度快(2内存需要少(3计算结果有良好的可靠性和可信性(4适应性好,亦即能处理变压器变比调整、系统元件的不同描述和与其它程序配合的能力强(5简单。 MATLAB 是一种交互式、面向对象的程序设计语言,广泛应用于工业界与学术界,主要用于矩阵运算,同时在数值分析、自动控制模拟、数字信号处理、动态分析、绘图等方面也具有强大的功能。MATLAB程序设计语言结构完整,且具有优良的移植性,它的基本数据元素是不需要定义的数组。它可以高效率地解决工业计算问题,特别是关于矩阵和矢量的计算。MATLAB与C语言和FORTRAN语言相比更容易被掌握。通过M语言,可以用类似数学公式的方

16、式来编写算法,大大降低了程序所需的难度并节省了时间,从而可把主要的精力集中在算法的构思而不是编程上。另外,MATLAB提供了一种特殊的工具:工具箱(TOOLBOXES.这些工具箱主要包括:信号处理(SIGNAL PROCESSING、控制系统(CONTROL SYSTEMS、神经网络(NEURAL NETWORKS、模糊逻辑(FUZZY LOGIC、小波(WA VELETS和模拟(SIMULATION等等。不同领域、不同层次的用户通过相应工具的学习和应用,可以方便地进行计算、分析及设计工作。MATLAB设计中,原始数据的填写格式是很关键的一个环节,它与程序使用的方便性和灵活性有着直接的关系。原

17、始数据输入格式的设计,主要应从使用的角度出发,原则是简单明了,便于修改。常数与矩阵的运算即是同该矩阵的每一元素进行运算。但需注意进行数除时,常数通常只能做除数。基本函数运算中,矩阵的函数运算是矩阵运算中最实用的部分,常用的主要有以下几个:det(a 求矩阵a的行列式eig(a 求矩阵a的特征值inv(a或a (-1 求矩阵a的逆矩阵rank(a 求矩阵a的秩trace(a 求矩阵a的迹(对角线元素之和我们在进行工程计算时常常遇到矩阵对应元素之间的运算。这种运算不同于前面讲的数学运算,为有所区别,我们称之为数组运算。数组的加、减与矩阵的加、减运算完全相同。而乘除法运算有相当大的区别,数组的乘除法

18、是指两同维数组对应元素之间的乘除法,它们的运算符为“.*”和“./”或“.”。前面讲过常数与矩阵的除法运算中常数只能做除数。在数组运算中有了“对应关系”的规定,数组与常数之间的除法运算没有任何限制。另外,矩阵的数组运算中还有幂运算(运算符为. 、指数运算(exp、对数运算(log、和开方运算(sqrt等。有了“对应元素”的规定,数组的运算实质上就是针对数组内部的每个元素进行的。矩阵的幂运算与数组的幂运算有很大的区别。逻辑运算是MATLAB中数组运算所特有的一种运算形式,也是几乎所有的高级语言普遍适用的一种运算。第五章潮流计算流程图及源程序5.1 潮流计算流程图 图5-1 潮流计算流程图5.2

19、潮流计算源程序据课题题目,本程序把节点1设为平衡节点,节点2、3、4为PQ节点,节点5为PV节点。G(1,1=10.834;B(1,1=-32.500;G(1,2=-1.667;B(1,2=5.000;G(1,3=-1.667;B(1,3=5.000;G(1,4=-2.500;B(1,4=7.500;G(1,5=-5.000;B(1,5=15.000;G(2,1=-1.667;B(2,1=5.000;G(2,2=12.917;B(2,2=-38.750;G(2,3=-10.000;B(2,3=30.000;G(2,4=0;B(2,4=0;G(2,5=-1.250;B(2,5=3.750;G(3

20、,1=-1.667;B(3,1=5.000;G(3,2=-10.000;B(3,2=30.000;G(3,3=12.917;B(3,3=-38.750;G(3,4=-1.250;B(3,4=3.750;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.250; B(4,3=3.750; G(4,4=3.750; B(4,4=-11.250; G(4,5=0;B(4,5=0;G(5,1=-5.000; B(5,1=15.000; G(5,2=-1.250; B(5,2=3.750; G(5,3=0;B(5,3=

21、0;G(5,4=0;B(5,4=0;G(5,5=6.250; B(5,5=-18.750; 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;k=0;precision=1;N1=4; %the N1 is the amount of the PQ buswhile precision>0.00001delt(5=0;u(5=1.06

22、;for m=1:N1for n=1:N1+1pt(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;endpp(m=p(m-sum(pt;qq(m=q(m-sum(qt;endfor m=1:N1for n=1:N1+1h0(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

23、-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;endH(m,m=sum(h0-u(m2*(G(m,m*sin(delt(m-delt(m-B(m,m*cos(delt(m-delt(m;N(m,m=sum(n0-2*u(m2*G(m,m+u(m2*(G(m,m*cos(delt(m-delt(m+B(m,m*sin(delt(m

24、-delt(m;J(m,m=sum(j0+u(m2*(G(m,m*cos(delt(m-delt(m+B(m,m*sin(delt(m-delt(m;L(m,m=sum(L0+2*u(m2*B(m,m+u(m2*(G(m,m*sin(delt(m-delt(m-B(m,m*cos(delt(m-delt (m;endfor m=1:N1JJ(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;endfor m=1:N1for n=1:N1if m=nelseH(m,n=-u(m*u(n*(G(m,n*

25、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;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;endendendfor m=1:N1PP(2*m-1=pp(m;PP(2*m=qq(m;enduu=-inv(JJ*PP'precision=max(abs(uu;for n=1

26、:N1delt(n=delt(n+uu(2*n-1;u(n=u(n+uu(2*n;endk=k+1;endK=k-1,delt,u'%the following program is used to calculate the S5 and S(m,n for n=1:N1+1U(n=u(n*(cos(delt(n+j*sin(delt(n;endfor m=1:N1+1I(m=Y(5,m*U(m;endS5=U(5*sum(conj(Ifor m=1:N1+1for n=1:N1+1S(m,n=U(m*(conj(U(m-conj(U(n*conj(-Y(m,n;endendYS5.

27、3 运行计算结果K =4;delt = -0.0461 -0.0839 -0.0896 -0.1044 0;U = 1.0365 1.0087 1.0073 1.0016 1.0600S5 =1.2982 + 0.2445iY =10.8340 -32.5000i -1.6670 + 5.0000i -1.6670 + 5.0000i -2.5000 + 7.5000i -5.0000 +15.0000i-1.6670 + 5.0000i 12.9170 -38.7500i -10.0000 +30.0000i 0 -1.2500 + 3.7500i-1.6670 + 5.0000i -10.0000 +30.0000i 12.9170 -38.7500i -1.2500 + 3.75