基于MATLAB的电力系统潮流计算分析与仿真

1基于 MATLAB 的电力系统潮流计算分析与仿真摘 要 本文首先对电力系统潮流计算进行了深入的研究，熟练掌握了牛顿拉夫逊极坐标形式和直角坐标形式的两种算法的原理。接着，为电力系统网络中各元件建立了等效的数学模型：输电线路用 型等值电路等效，考虑了导线对地导纳；变压器用 型等值电路等效。进一步，为电力系统网络建立极坐标与直角坐标形式的牛顿拉夫逊潮流计算数学模型。在此基础上，运用 MATLAB 语言编写了基于牛顿拉夫逊方法的潮流计算程序，该程序可以实现原始数据输入、自动形成节点导纳矩阵、形成雅可比矩阵和求解满足精度要求的各节点电压和功率等功能。为验证程序的准确性，将该程序应用于某三节点环形电力网络潮流计算，并与手工算法求得的理论结果进行对比。最后，采用该程序对某复杂环形电力网络进行潮流计算。关键词 潮流计算，牛顿拉夫逊法，MATLAB，电力系统1 前言潮流计算是研究电力系统稳态运行情况的一种重要的分析计算,随着电力系统越来越复杂，使得计算机计算潮流成为主流，目前计算机计算潮流有三种算法：高斯赛德尔、牛顿-拉夫逊、PQ分解方法 1,2,3。而牛顿- 拉夫逊法是电力系统潮流计算的常用算法之一，它收敛性好，迭代次数少,精度高。且 MATLAB 使用方便，有着其他高级语言无法比拟的强大的矩阵处理功能 4，所以本毕设研究的方向为基于 MATLAB 的电力系统牛顿- 拉夫逊法潮流计算分析与仿真。本毕设在实现潮流计算基本功能的同时，还扩展了一些新的功能。具体如下：根据从命令窗口输入 i 的值来选择算法，增加了程序算法可选择性；扩展了根据自己的需要规定相应的精度 eps的功能；根据导纳矩阵形成原理自定义了函数 Y=matrix(y,n,y0,Ltype,KT)来形成导纳矩阵，省去了以往手工算导纳矩阵花费的大量精力；增加了各支路间流动功率，各节点注入电流和各节点间线路功率损耗的计算功能；将已知线路串联的导纳 y，线路并联的导纳 y0，节点的功率 S，线路类型 Ltype（有变压器时置 2，无变压器时置 1）,变压器变比 KT，存入相应矩阵中，并将矩阵保存为 mat 文件，程序运行前通过加载命令 load 调入。避免每次计算时重复输入线路参数。2 潮流计算的数学模型2.1 潮流计算的节点模型根据电力系统中各节点性质的不同，很自然地把节点分成三类 1,2,3,5：（1） PQ 节点基于 MATLAB 的电力系统潮流计算分析与仿真2对这一类点，事先给定的是节点功率(P，Q)，待求的未知量是节点电压向量 (U， )，所以叫PQ 节点。（2） PV 节点这类节点给出的参数是该节点的有功功率 P 及电压幅值 U，待求量为该节点的无功功率 Q 及电压向量的相角 。（3） 平衡节点在潮流计算中，这类节点一般只设一个。对该节点，给定其电压值，并在计算中取该节点电压向量的方向作为参考轴，相当于给定该点电压向量的角度为零。也就是说，对平衡节点给定的运行参数是 U 和 ，因此待求量是该节点的 P 和 Q，整个系统的功率平衡由这一节点承担。以上三类节点 4 个运行参数 P、Q、U 、 中，已知量都是两个，待求量也是两个，只是类型不同而已。2.2 修正方程采用牛顿法计算潮流时，需要对功率方程进行修改。下面将根据在不同坐标内的修改进行讨论：当在极坐标系内时，计算的 PQ 节点功率与给定的 PQ 节点的注入功率 的差值满足下isQP,面方程：（21） 0)cossin(ic1nj ijijjjiiij jjijijiii BGUQP式中： i 与 j 节点电压的相角差。jiij由方程（21）可以写出修正方程式如下:VLKNQPD12（22）3其中：H 是 阶方阵，其元素为 ；N 是 阶矩阵，其元素为)1(njiijPmn)1(；K 是 阶矩阵，其元素为 ；L 是 阶矩阵，其元素为ijijVPNmjiijQK。jjijQL对式 （21）求偏导数，可得雅克比矩阵元素的表达式如下：非对角元素（ ）i（23）)cossin(ic)cssi(ijijjjiijj ijjijijiij iBGVLKNH对角元素（ ）ji（24）iii iii iiijQBVLPGKNH22在直角坐标系内时，计算的 PQ 节点功率给定的 PQ 节点的注入功率 的差值应满足以is,下方程：（25）0)()(11nj nj jijiiijjiiisiisi j j jijiiijjiiisiisi eBfGefBeGfQffP对于 PV 节点，已知节点的注入有功功率及节点电压大小，记作 ，其节点的有功功率应满isUP,方程：基于 MATLAB 的电力系统潮流计算分析与仿真4)1,.21(0(0)()(211 nmifeUeBfGfBeGPiiisnj nj jijiiijjiiisiisi（26）对于平衡节点，因为其电压给定，故不需要迭代求解。由方程（25）和（26）可以得到修正方程如下：UJW（27）式中 TnmmUPPQPW212111 nfefeffeU 12122112 1111 21212121211211 11111 111111 11111 nmnmnn nnmmmm nnmmmm nnnnmm fUefeUfefUe PPPP fefefefe PPPP fQefQeffQe fefeffeJ 5其中：雅克比矩阵的各元素可以对式（25）和式（26）求偏导数获得。对于非对角元素（ ）有ji（28）0)(22jiji ijijjiji ijijjijifUefGeBQfPffe对于对角元素（ 有)ji（29）iiiinj iijijiij iijijiinj iijijiij iijijiiffUefBeGffQBfe efefPGfe2)()(2112.3 线路和变压器等效模型I1. I2.ZU.2Y基于 MATLAB 的电力系统潮流计算分析与仿真6图 2.1 变压器 型等值电路 图 2.2 线路 型等值电路 3 电力系统潮流计算程序设计及算例3.1 电力系统网络潮流计算程序设计总体思路本毕业设计使用了牛顿拉夫逊极坐标和直角坐标的两种方法，可以根据已知电压的不同表达形式选用两种方法中的对应的形式，省去了繁琐的手工转化。本毕业设计实现了电力网络原始数据输入程序设计、形成节点导纳矩阵的程序设计、形成雅克比矩阵程序设计，并将该程序应用于某一简单环形电力网络潮流计算，将计算结果与手工算法的理论结果进行对比，验证程序的准确性，最终采用该程序求解更为复杂环形电力网络潮流计算。从而证明该程序适用于任意复杂电力网络的潮流计算。3.2 潮流计算程序设计3.2.1 牛顿法潮流计算主要流程1) 形成节点导纳矩阵；2) 给各节点电压设初值,极坐标表示的 PQ 节点的电压相角一律为 0，电压幅值一律为 1，直角坐标表示的 PQ 节点的电压实部一律为 1，电压虚部一律为 0；（标幺值情况下），极坐标表示的 PQ 节点的电压相角一律为 0，电压幅值一律为额定电压，直角坐标表示的PQ 节点的电压实部一律为额定电压，电压虚部一律为 0；（真实值情况下）3) 将节点电压初值代入（21）或（25），求出不同坐标形势下修正方程式的常数项向量 ， ；iPiQ4) 将节点电压初值代入（23）、（24）或（28）、（29），求出不同坐标形式下的雅可比矩阵元素；5) 求解修正方程式（22）或（27），求出不同坐标形式下的电压变量的修正向量 ；6) 求出节点电压的新值；7) 检查是否收敛，若收敛，电压趋近于真解时，功率偏移量将趋于零。若不收敛，则以各节点电压的新值作为初值自第 3 步重新开始下一次迭代，否则转入下一步。8) 计算支路功率分布，PV 节点无功功率和平衡节点注入功率，最后输出结果，并结束。7P1 123456 7输 入 原 始 数 据形 成 节 点 导 纳 矩 阵 Y给 定 节 点 电 压 迭 代 初 值置 迭 代 次 数 k= 0计 算 节 点 偏 移 量计 算 结 束收 敛 判 断 ？计 算 J矩 阵 元 素解 修 正 方 程 求 节 点 电 压 修 正 量修 正 节 点 电 压 获 得 新 的 近 似 解计 算 PV 、 QV 、 QPV ； 计 算 线 路 潮 流 ； 计 算 全 网 功 率 损 耗迭 代 次 数k=+1 YesNo准 备2P3 牛顿-拉夫逊潮流计算程序框图 6如图 3.1 所示。图 3.1 程序流程图其中： 是平衡节点的功率， 是 PV 节点的无功功率v, pv3.2.2 电力网络原始输入数据设计潮流计算中用到的线路参数包括节点间的串联导纳、节点对地并联导纳，在本程序运行前将线路参数存入矩阵y和y0中，并将矩阵保存为mat文件，形成导纳矩阵时通过加载命令load调入。同样将已知节点的功率S，线路类型Ltype（有变压器时置2，无变压器时置1）,变压器变比KT，存入相应矩阵中，并将矩阵保存为mat文件，程序运行前通过加载命令load调入。然后根据需要选择主函数输入变量的个数，例如只有2个输入变量PowerFlow(y, S）时，则默认调用极坐标形式的子函数function jizuobiao(y, S, eps),然后主函数会通过jizuobiao(y, S)调用语句将实参变量y、S传递给子函数jizuobiao(y, S, eps)中的形参变量y、S，并利用 nargin函数统计变量个数功能计算出变量个数为2，最后利用if (nargin = 2)语句判断变量个数是否为 2，如果是，则自动给精度eps 赋默认值；3个输入变量PowerFlow(y, S，i ）时，则根据命令窗口输入的i值来依此选择是运行极坐标子函数function jizuobiao(y, S, eps)（当i=1时)，还是直角坐标形式的子函数function zhijiaozuobiao(y, S, eps)(当i=2时)，并按上述传递实参变量给形参变量后统计变量个数的原理给精度eps 赋默认值，如果输入的i不在范围内，则在窗口显示 只有两种方法，你有更好的方法？此外本程序还可以根据自己的需要规定相应的精度，即在命令窗口直接输入eps 的值后，输入4个变量 PowerFlow(y, S, i, eps)，同样根据命令窗口输入的i 值来依此选择是运行极坐标子函数function jizuobiao(y, S, eps)（当i=1 时) ，还是直角坐标子函数function zhijiaozuobiao(y, S, eps)(当i=2时)，如果输入的i不在范围内，则在窗口显示只有两种方法，你有更好的方法？具体程序如下：%潮流计算主函数 function PowerFlow(y, S, i, eps)基于 MATLAB 的电力系统潮流计算分析与仿真8clcif (nargin = 2)jizuobiao(y, S, eps);elseif (nargin = 3)if i = 1jizuobiao(y, S, eps);elseif i = 2zhijiaozuobiao(y, S, eps);elsedisp(只有两种方法，你有更好的方法？);endelse (nargin = 4)if i = 1jizuobiao(y, S, eps);elseif i = 2zhijiaozuobiao(y, S, eps);elsedisp(只有两种方法，你有更好的方法？);endEnd%潮流计算极坐标形式下的子函数function jizuobiao(y, S, eps)if nargin = 2eps = 1.0e-4;end%潮流计算直角坐标形式下的子函数function zhijiaozuobiao(y, S, eps)if nargin = 2eps = 1.0e-4;end3.2.3 形成节点导纳矩阵的程序设计根据命令窗口显示的语句请输入 PQ 节点的个数：和语句请输入 PV 节点的个数：依次给 PQ 节点数（nPQ) 和 PV 节点数(nPV)赋值后，系统会根据输入的节点数和线路参数 串联导纳 y，并联导纳 y0，线路类型 Ltype（有变压器时置 2，无变压器时置 1），变压器变比 KT，来自动调用按照导纳矩阵形成原理定义的函数 Y=matrix(y,n,y0,Ltype,KT)形成导纳矩阵，并通过取精度命令 Y=round(Y*100)/100 来选择合适的精度，最后应用 real 和 imag 函数依次将导纳矩阵实部和虚部分离，赋值给电导矩阵 G 和电纳矩阵 B。具体程序如下：nPQ=input(请输入PQ节点的个数：);nPV=input(请输入PV节点的个数：);n=nPQ+nPV+1;load y0.matload Ltype.matload KT.matY=matrix(y,n,y0,Ltype,KT);Y=round(Y*100)/100;G=real(Y);B=imag(Y);其中导纳矩阵的形成：9求导纳矩阵 Y 中的非对角元元素 Yij，若无变压器,则 Yij 直接为线路导纳取负值，若有变压器，Yij 为线路导纳乘以 KT 分之一再取负值。求导纳矩阵Y中对角元元素 Yii，无变压器时Yii 为Yij 加上线路对地电导的一半乘 j，有变压器时，对角元元素就与所输入的折算到哪一侧有关，如果支路起始端处于高压侧，支路起始节点的自导纳中要加上变压器等值导纳模型的对地支路的(1-KT)/KT2倍，支路终止节点的自导纳要加上变压器等值导纳模型的对地支路的(KT-1)/KT倍，如果支路起始端处于低压侧，情况正好相反支路起始节点的自导纳中要加上变压器等值导纳模型的对地支路的(KT-1)/KT倍，支路终止节点的自导纳要加上变压器等值导纳模型的对地支路的(1-KT)/KT2倍。具体程序如下：function Y=matrix(y,n,y0,Ltype,KT)for i=1:nfor j=1:nswitch Ltype(i,j)case 1y(i,j)=y(j,i);y0(i,j)=y0(j,i);case 2y(i,j)=y(i,j)*1/KT(i,j);y(i,j)=y(j,i); y0(i,j)=(1-KT(i,j)*y(i,j)/(KT(i,j)2); y0(j,i)=(KT(i,j)-1)*y(i,j)/KT(i,j);endend end%求互导纳for i=1:nfor j=1:nif i=jY(i,j)=-y(i,j);endendend%求自导纳for i=1:nY(i,i)=sum(y(i,:)+sum(y0(i,:);endY %Y 为导纳矩阵3.2.4 形成雅克比矩阵的程序设计1、极坐标形式：根据计算雅克比矩阵元素公式，并运用for i=1:n-1循环和for j=1:n-1循环嵌套结构，区分 和 两种情况依次对Hi,j,Ni,j,Ki,j,Li,j 计算并且赋值，最后依据矩阵每两行或两jiji列互换后矩阵不变原理，对矩阵进行分块化简，形成了简化后的雅克比矩阵。具体程序如下：%求雅克比矩阵%当i=j时候求H,N,M,L 如下：for i=1:n-1for j=1:n-1if i=j H(i,j)=-U(i)*U(j)*(G(i,j)*sin(delt(i)-delt(j)-B(i,j)*cos(delt(i)-delt(j);N(i,j)=-U(i)*U(j)*(G(i,j)*cos(delt(i)-基于 MATLAB 的电力系统潮流计算分析与仿真10delt(j)+B(i,j)*sin(delt(i)-delt(j);K(i,j)=-N(i,j);L(i,j)=H(i,j);endendendH,N,K,L%当i=j 时H,N,M,L如下：for i=1:n-1for j=1:nP(j)=U(i)*U(j)*(G(i,j)*cos(delt(i)-delt(j)+B(i,j)*sin(delt(i)-delt(j);Q(j)=U(i)*U(j)*(G(i,j)*sin(delt(i)-delt(j)-B(