数值分析应用

1、西安理工大学研究生课程论文/研究报告课程名称： 数值分析 任课教师： 闵涛 论文/研究报告题目： 数值分析在潮流计算 方面的应用 完成日期： 2012 年 12 月 26 日学 科： 电力系统及其自动化 学 号： 1208080915 姓 名： 黄亭 成 绩： 数值分析在潮流计算方面的应用1 潮流计算的计算机算法简介潮流计算的计算机算法是以电网络理论为基础的，应用数值计算方法求解一组描述电力系统稳态特性的方程。从数学上讲是一组多元的非线性方程式的求解问题，这类方程的求解过程都离不开迭代。由于电力系统结构及参数的一些特点，同时随着电力系统不断扩大，潮流问题的方程式的阶数也越来越高，这样的非线性方

2、程式并不是任何数学方法都能保证给出正确答案的。这种情况就成为促使电力系统计算人员不断寻求新的且更可靠方法的一个重要因素。电网潮流计算的性能优劣一般依据的是能否可靠收敛，计算速度的快慢，内存占有多少，使用是否方便灵活，调整和修改是否容易，是否满足工程需要等来判别，其中以是否可靠收敛作为评价的主要标准。常用的分析法包括高斯-塞德尔法、牛顿-拉夫逊潮流算法、快速解耦算法（PQ 分解法）等。2 潮流计算的约束条件电力系统运行必须满足一定技术和经济上的要求。这些要求够成了潮流问题中某些变量的约束条件，常用的约束条件如下：2.1节点电压应满足: （2.1） 2.2节点的有功功率和无功功率应满足： （2.2

3、） 2.3节点之间电压的相位差应满足： （2.3） 3 节点导纳矩阵的形成与修改3.1 节点电压方程 （1）自、互导纳的物理意义自导纳在数值上等于与该节点I直接连接的所有支路导纳的总和。如。互导纳在数值上等于连接节点、支路导纳的负值，即。如。（2）节点导纳矩阵YB为对称方阵。（3）节点导纳矩阵YB为稀疏矩阵。（4）节点导纳矩阵具有对角优势。3.2 节点导纳矩阵的形成 用直接形成法形成节点导纳矩阵YB。节点导纳矩阵即可根据自导纳和互导纳的定义直接形成，也可用支路节点关联矩阵计算。 3.3 节点导纳矩阵的修改（1）从原有网络引出一支路，同时增加一节点，节点导纳矩阵将增加一阶。新增的对角元，；新增的

4、非对角元，；原有矩阵中的对角元将增加 ，。（2）在原有网络的节点、之间增加一支路。 ，（3）在原有网络的节点，之间切除一支路，（4）原有网络的节点、之间的导纳由改变为：，（5）原有网络节点i、j之间变压器的变比由改变为 ；4 牛顿-拉夫逊法（直角坐标）4.1概述1. 牛顿-拉夫逊法的意义和推导过程把按泰勒级数在点展开 修正方程 （4.1）2牛顿拉夫逊法的特点(1)牛顿-拉夫逊法是迭代法，逐渐逼近的方法；(2)修正方程是线性化方程，它的线性化过程体现在把非线性方程在按泰勒级数展开，并略去高阶小量；(3)用牛顿拉夫逊法解题时，其初始值要求严格(较接近真解)，否则迭代不收敛。3多变量非线性方程的解牛

5、顿拉夫逊法修正方程： （4.2）缩写为 (4.3) 4.2潮流计算时的修正方程（直角坐标）PQ节点 (4.4)PV节点 (4.5) 平衡节点只设一个,电压为已知,不参见迭代,其电压为 (4.6) 修正方程 (4.7) (4.8) 4.3雅可比矩阵各元素当时, 雅可比矩阵中非对角元素为 (4.9)当时,雅可比矩阵中对角元素为 (4.10)图1 牛顿-拉夫逊法计算步骤5 用MATLAB进行编程 牛顿-拉夫逊法（直角坐标）5.1 某电网接线图及给定的参数其中，1,2,3,4为PQ节点，5为平衡节点各支路阻抗：Z12=Z21=0.06+j0.18 Z13=Z31=0.06+j0.18 Z14=Z41=

6、0.04+j0.12Z15=Z51=0.02+j0.06 Z23=Z32=0.01+j0.03 Z25=Z52=0.08+j0.24 Z34=Z43=0.08+j0.24 各节点输出功率1：-0.2-j0.22: 0.45+0.15 3: 0.4+j0.054: 0.6+j0.1 5: 05.2 潮流计算计算机算法流程图开始形成节点导纳矩阵输入原始数据设节点电压，i=1,2,n,is置迭代次数置节点号i=1按式（4.9），（4.10）计算雅克比矩阵元素按式（4.4）计算节点的，节点的，求解修正方程式，得，雅克比矩阵是否已全部形成？计算平衡节点及PV节点功率求，迭代次数 k=k+1i=i+1？潮

7、流计算完成计算各节点电压的新值：结 论本文采用牛顿-拉夫逊法（直角坐标）进行潮流计算。牛顿迭代法是求方程根的重要方法之一，是常用的解非线性方程组的方法，也是当前广泛采用的计算潮流的方法。其最大优点是在方程的单根附近具有平方收敛，收敛性较好，而且该法还可以用来求方程的重根、复根。但是对初始值的要求比较严格。用MATLAB编程并仿真，MATLAB语言允许用户以数学形式的语言编写程序, 其比BASIC 语言和FORTRAN 等更为接近书写的数学表达格式, 且程序易于调试。在计算要求相同的情况下, 使用MATLAB 编程, 工作量将会大为减少。其程序的编写也因MATLAB提供了许多功能函数而变得简单易

8、行。参 考 文 献1 陈衍.电力系统稳态分析M.北京：中国电力出版社，2007. 2 李维波. MATLAB 在电气工程中的应用M. 北京: 中国电力出版社，2007.3 武晓朦，张飞廷. 电力系统的P- Q 分解法潮流计算J.现代电子技术,2002 ,142(11): 105-106 4 谢威，彭志炜，张朝纲，马春生.一种基于牛顿拉夫逊的潮流就算方法J.许昌学院学报. 2006.3,25(2):27-305 杨帆.电力系统潮流计算程序设计J.山西冶金，2007，106（2）：42-446 刘军,刘学军. MATLAB 在电力系统分析中的应用J. 电力系统及其自动化学报, 2000.4,12(

9、2): 23-257 李家坤,刘姣姣. MATLAB仿真技术在电力电子教学中的应用J. 长江工程职业技术学院学院2009.3,26(1): 75-778 连小洲，焦洁. MATLAB在电力系统计算中的应用J. 江西电力职业技术学院学报，2005.12,18(4): 10-119 张宁,江红梅,张渭. 基于MATLAB的电力系统潮流计算J.西北农林科技大学学报(自然科学版)，2004.12,32(12): 124-12610 周卫星，张颖.基于Matlab的电力系统潮流计算J.科技咨询导报，2007，10：70-7111 徐一哲，沈瑞寒.基于Matlab的电力系统潮流分析J.中外企业家，2009

10、,5(下)：206-20812 赵芳,陈士方,张圣集.保留非线性的电力系统概率潮流计算J.电力情报，2000,03：22-25附录A 基于MATLAB的牛顿拉夫逊法（直角坐标）潮流计算程序清单clear;clcy=0;%输入原始数据，求节点导纳矩阵y (1,2)=1/(0.06+0.18i); y (1,3)=1/(0.06+0.18i); y (1,4)=1/(0.04+0.12i);y(1,5)=1/(0.02+0.06i);y(2,3)=1/(0.01+0.03i);y(2,5)=1/(0.08+0.24i);y(3,4)=1/(0.08+0.24i);y(4,5)=0;for i=1:

11、5for j=i:5y(j,i)=y(i,j);endendY=0;%求互导纳for i=1:5for j=1:5if i=jY(i,j)=-y(i,j);endendend%求自导纳for i=1:5Y(i,i)=sum(y(i,:);endY %Y 为导纳矩阵G=real(Y);B=imag(Y);%原始节点功率S(1)=0.2+0.2i;S(2)=-0.45-0.15i;S(3)=-0.4-0.05i;S(4)=-0.6-0.1i;S(5)=0;P=real(S);Q=imag(S);%赋初值U=ones(1,5);U(5)=1.06;e=zeros(1,5);ox=ones(8,1);

12、fx=ones(8,1);count=0 %计算迭代次数while max(fx)1e-5for i=1:4 for j=1:4 H(i,j)=0;N(i,j)=0;M(i,j)=0;L(i,j)=0;oP(i)=0;oQ(i)=0; endendfor i=1:4for j=1:5 oP(i)=oP(i)-U(i)*U(j)*(G(i,j)*cos(e(i)-e(j)+B(i,j)*sin(e(i)-e(j);oQ(i)=oQ(i)-U(i)*U(j)*(G(i,j)*sin(e(i)-e(j)-B(i,j)*cos(e(i)-e(j);endoP(i)=oP(i)+P(i); oQ(i)=

13、oQ(i)+Q(i);endfx=oP,oQ;%求雅克比矩阵%当i=j时候求H,N,M,L 如下：for i=1:4for j=1:4if i=j H(i,j)=-U(i)*U(j)*(G(i,j)*sin(e(i)-e(j)-B(i,j)*cos(e(i)-e(j);N(i,j)=-U(i)*U(j)*(G(i,j)*cos(e(i)-e(j)+B(i,j)*sin(e(i)-e(j);L(i,j)=H(i,j);M(i,j)=-N(i,j);endendendH,N,M,L%当i=j 时H,N,M,L如下：for i=1:4for j=1:5if i=jH(i,i)=H(i,i)+U(i)

14、*U(j)*(G(i,j)*sin(e(i)-e(j)-B(i, j)*cos (e(i)-e(j); N(i,i)=N(i,i)-U(i)*U(j)*(G(i, j)*cos(e(i)-e(j)+B(i,j)*sin(e(i)-e(j);M(i,i)=M(i,i)-U(i)*U(j)*(G(i,j)*cos(e(i)-e(j)+B(i,j)*sin(e(i)-e(j); L(i,i)=L(i,i)-U(i)*U(j)*(G(i,j)*sin(e(i)-e(j)-B(i,j)*cos(e(i)-e(j);endendN(i,i)=N(i,i)-2*(U(i)2*G(i,i);L(i,i)=L(i,i)+2*(U(i)2*B(i,i);endJ=H,N;M,L %J 为雅克比矩阵ox=-(inv(J)*fx);for i=1:4oe(i)=ox(i); oU(i)=ox(i+4)*U(i);endfor i=1:4e(i)=e(i)+oe(i); U(i)=U(i)+