基于半隐式龙格库塔法的电力系统暂态稳定计算方法与流程

本发明属于电力系统暂态稳定仿真与分析
技术领域：
，具体涉及一种基于半隐式龙格库塔法的电力系统暂态稳定计算方法。
背景技术：
：随着智能电网和交直流混合电网的发展，电力系统运行情况变得更加复杂多变。电力系统时常遭受到不同程度的故障和扰动，为了保障电力系统能够持续的安全稳定运行，需要提前判断电力系统在不同的故障扰动下能否稳定运行，因此需要进行暂态稳定仿真和分析。暂态稳定分析一般指系统受到大扰动之后的机电暂态过程，持续时间为几秒到十几秒，快速精准的进行暂态稳定仿真，能够辅助电网工作人员提前掌控电网运行趋势，及时制定有效的安全控制策略从而防止更大规模电力事故的发生，提升电网的安全稳定运行水平。目前的暂态稳定计算中广泛使用的数值积分算法有改进欧拉法，隐式梯形积分法，四阶显式龙格库塔法等，这些算法也广泛使用于暂态稳定仿真的商业化软件中，如电力系统综合分析程序psasp，中国版的bpa，美国pti公司的pss/e等。这些积分算法有各自独特的优势，但也存在不同的缺陷，比如显式积分算法计算简单，计算效率高，但数值稳定性差，不具备a-稳定，其仿真步长受到稳定性的极大限制，也很难适用于刚性系统的暂态稳定仿真；隐式积分算法数值稳定性好，具备a-稳定，但需要进行迭代计算，计算较复杂，计算速度慢。此外，目前使用最广泛的隐式梯形积分法，虽然具备a-稳定，其仿真步长可以不受稳定性要求的限制，但由于在暂态稳定计算中使用了简单迭代法，其仿真步长要受到迭代收敛性的限制，导致其仍然不能使用较大步长，在暂态稳定仿真中一般取步长为0.01s，随着步长的增大会出现数值振荡的问题。因此，为了快速精准的进行暂态稳定仿真和分析，亟需提出一种数值稳定性好，计算精度高，计算速度快，能够使用较大步长的数值积分算法。技术实现要素：本发明正是针对现有技术中的问题，提供了基于半隐式龙格库塔法的电力系统暂态稳定计算方法，建立电力系统暂态稳定仿真计算的动态元件模型，形成全系统微分代数方程组，输入原始数据和信息，进行潮流计算，得到系统扰动前的运行参量初值，并计算状态变量初始值，生成各动态元件的雅克比矩阵和网络导纳矩阵，随后经过系统扰动判断后，生成全系统微分方程组的雅克比矩阵，求解得到下一时刻各状态变量的值，机网交替迭代计算，直至满足收敛条件，得到下一时刻的运行参量值，以下一时刻各状态变量的值和运行参量值为初值，进行下一步长的计算，既具备隐式积分算法的数值稳定性，又无需进行迭代计算，兼具了显式积分算法计算简单，计算量小的优势，同时相比目前使用最广泛的隐式梯形积分法能够使用更大的步长，解决了暂态稳定仿真计算中算法数值稳定性和计算效率难以兼顾的问题，能够更加高效精确稳定的进行电力系统暂态稳定计算。为了实现上述目的，本发明采用的技术方案是：基于半隐式龙格库塔法的电力系统暂态稳定计算方法，包括如下步骤：s1，建立电力系统暂态稳定仿真计算的动态元件模型，形成系统微分方程式和代数方程式,所述系统微分方程式和代数方程式构成全系统微分代数方程组，具体为：其中，x为各动态元件的状态变量；v为动态元件接入点的网络节点电压；i为动态元件注入网络的电流；v,i均为复数；s2，输入原始数据和信息，进行潮流计算，得到系统扰动前的运行参量初值；s3，利用步骤s2获得的运行参量初值求解计算状态变量的初始值；s4，根据步骤s1中动态元件模型的微分方程组，计算动态元件的雅克比矩阵，根据网络结构形成网络的节点导纳矩阵y，所述网络模型用网络导纳矩阵来描述，yv＝i(x,v)；s5，系统扰动判断，根据系统扰动的来源，修改微分方程式、代数方程式或者网络的节点导纳矩阵；s6，利用半隐式龙格库塔法求解步骤s1中的微分方程组，生成全系统微分方程组的雅克比矩阵j，求解得到下一时刻各状态变量的值，所述矩阵为：s7，机网交替迭代计算，直至满足收敛条件，得到下一时刻的运行参量值；s8，根据步骤s6获取的下一时刻各状态变量的值和步骤s7获得的下一时刻的运行参量值为初值，进行下一步长的计算。作为本发明的一种改进，所述步骤s2中运行参量初值包括系统各节点的电压和各动态元件注入网络的功率其中，vx(0),vy(0)为xy坐标系下节点电压的实部和虚部，p(0),q(0)为节点注入的有功和无功功率；通过计算得到各动态元件注入网络的电流初值所述电流初值为：作为本发明的另一种改进，所述步骤s5进一步包括：s51，判断系统是否发生故障或操作，若是，继续步骤；否则进入步骤s6；s52，判断故障来源，若扰动发生在动态元件内部，修改微分方程组和代数方程组，否则继续步骤s53；所述修改计算过程具体为：计算发生扰动之后各状态变量的初值x(tf)，扰动前后状态变量不能突变，当t<tf<t+△t时,仿真步长变为hf＝tf-t，重新求解微分方程组x(tf)＝x(t+hf)；当tf＝t+△t时，x(tf)＝x(t+△t)；其中，tf为扰动发生时刻，t为当前仿真时刻，△t为当前仿真步长；s53，修改网络的导纳矩阵，利用修改后新的导纳矩阵重新计算网络方程。作为本发明的又一种改进，所述步骤s6进一步包括：s61，生成全系统微分方程组的雅克比矩阵j；s62，计算得到稀疏线性方程组，所述稀疏线性方程组为：其中wi,aij,bi,cij为常数系数，xn+1＝x(tn+1),xn＝x(tn)，h为仿真步长，tn+1＝tn+h；依次求解n维稀疏线性方程组，得到系数矩阵k1,k2,…ks，求解得到下一时刻各状态变量值x(t+△t)；s63，根据下一时刻各状态变量值计算动态元件注入网络的电流值i，所述电流值i为：i＝f(x(t+△t),v(t))。作为本发明的又一种改进，所述步骤s7中进一步包括：s71，根据步骤s63中动态元件注入网络的电流值i(0)＝f(x(t+△t),v(0))求解网络模型yv＝i(0)＝f(x(t+△t),v(0))，得到v(1)，其中v(0)＝v(t)；s72，继续求解i(1)＝f(x(t+△t),v(1))，进而得到v(2)；s73，依次进行迭代计算，当收敛条件时，停止迭代，得到y(t+△t)＝y(k+1)，获得下一时刻的运行参量值。作为本发明的更进一步改进，所述步骤s73中收敛条件为：|v(k+1)-v(k)|＜ε其中，k为迭代次数；ε为预设误差限。与现有技术相比，本发明所提出的基于半隐式龙格库塔法的电力系统暂态稳定计算方法，针对电力系统暂态稳定仿真传统数值积分算法的缺点，将rosenbrock线性化策略运用到传统龙格库塔积分格式中，获得了一种既具备隐式积分算法a-稳定的数值稳定性，又无需进行隐式积分算法的迭代计算，具备显式积分算法计算简单，计算效率高的特性，同时也能很好地适用于刚性电力系统的暂态仿真，并且能够使用大步长，数值稳定性好，计算精度更高，计算效率更快，解决了暂态稳定仿真计算中算法数值稳定性和计算效率难以兼顾的问题，能够更加高效精确稳定的进行电力系统暂态稳定计算。附图说明图1是本发明基于半隐式龙格库塔法的电力系统暂态稳定计算方法流程图；图2是本发明实施例2中3机9节点电力系统的结构接线图；图3是本发明三相短路故障节点7电压在半隐式龙格库塔法和隐式梯形积分法的仿真结果图；图4是本发明三相短路故障发电机g2转速在半隐式龙格库塔法和隐式梯形积分法仿真结果图。具体实施方式以下将结合附图和实施例，对本发明进行较为详细的说明。实施例1基于半隐式龙格库塔法的电力系统暂态稳定计算方法，如图1所示：s1，建立电力系统暂态稳定仿真计算所需的动态元件模型，包括发电机模型，励磁器模型，pss模型，原动机模型，负荷模型等；形成系统微分方程式和代数方程式，构成全系统动态仿真数学模型，由下式所示的微分代数方程组来描述：其中，x为描述动态元件特性的状态变量，y为代数方程组中系统的运行参量；s2，输入原始数据和信息，进行潮流计算，得到系统扰动前的运行参量的初值y(0)，包括系统各节点的电压和各动态元件注入网络的功率进而计算得到各动态元件注入网络的电流初值其中，s3，利用y(0)计算各状态变量的初值x(0)；s4，计算各动态元件模型的雅克比矩阵jd，根据网络结构生成网络的节点导纳矩阵y；s5，系统扰动判断，根据系统扰动的来源，修改微分方程式、代数方程式或者网络的节点导纳矩阵；s51，判断系统是否发生故障或操作，若是，继续步骤；否则进入步骤s6；s52，判断故障来源，若扰动发生在动态元件内部，修改微分方程组和代数方程组，否则继续步骤s53；所述修改计算过程具体为：计算发生扰动之后各状态变量的初值x(tf)，扰动前后状态变量不能突变，当t<tf<t+△t时,仿真步长变为hf＝tf-t，重新求解微分方程组x(tf)＝x(t+hf)；当tf＝t+△t时，x(tf)＝x(t+△t)；其中，tf为扰动发生时刻，t为当前仿真时刻，△t为当前仿真步长；s53，修改网络的导纳矩阵，利用修改后新的导纳矩阵重新计算网络方程；s6，运用半隐式龙格库塔法求解式(1)中的微分方程组，得到下一个积分时刻各状态变量的值x(t+△t).对于自治系统，其中，x＝(x1,x2,...xn)t,f＝(f1,f2,...,fn,)t，s级半隐式龙格库塔法的积分格式如下式所示：其中为系统微分方程组的雅克比矩阵，wi,aij,bi,cij为常数系数，xn+1＝x(tn+1),xn＝x(tn)，h为仿真步长，tn+1＝tn+h；观察式(3)，ki可由kj(j<i)显式表示出来，因此，式(3)可以等效写成，记则式(4)可写为其中，矩阵j和m均为稀疏矩阵，且m为非奇异矩阵，式(5)为n维稀疏线性方程组，求解该线性方程组可以通过矩阵求逆或lu分解的方法，由式(5)计算得到下一时刻各状态变量的值x(t+△t)。利用上一步中得到的各状态变量的值x(t+△t)计算动态元件注入网络的电流i(t+△t)。动态元件注入网络电流是描述其动态行为的状态变量和相应节点电压的函数，如式(6)所示。针对不同的动态元件模型，其注入网络的电流有不同的表达式。再进一步的，利用式(7)求解网络方程yv(t+δt)＝i(x(t+δt),v(t+δt))式(7)其中，y为网络导纳矩阵，v(t+△t)为待求的各节点电压。s7，式(7)为隐式方程，需要进行迭代求解，即以当前时刻的节点电压为迭代初值v(0)，首先利用式(6)计算得到i(0)，然后利用式(7)计算得到v(1)，依次继续进行计算，直到满足收敛条件(8)停止迭代，得到下一时刻的节点电压v(t+△t)。|v(k+1)-v(k)|＜ε式(8)其中，k为迭代次数，ε为预设误差限。s8，根据步骤s6获取的下一时刻各状态变量的值和步骤s7获得的下一时刻的运行参量值为初值，进行下一步长的计算。上述方法既具备隐式积分算法的数值稳定性，又无需进行迭代计算，兼具了显式积分算法计算简单，计算量小的优势，同时相比目前使用最广泛的隐式梯形积分法能够使用更大的步长，解决了暂态稳定仿真计算中算法数值稳定性和计算效率难以兼顾的问题，能够更加高效精确稳定的进行电力系统暂态稳定计算。实施例2以3机9节点电力系统为例，其单线结构图如图2所示，该系统中有3台发电机、3个负荷以及9条支路，其中发电机采用经典二阶模型，负荷采用恒定阻抗模拟，电力网络采用导纳矩阵描述，微分方程采用二级三阶半隐式龙格库塔法求解，网络方程采用直接法求解。采用上述模型的整个暂态稳定计算的流程如下：建立动态模型，形成系统的微分-代数方程组，所述系统中的动态元件模型只有发电机，其经典二阶发电机模型为其中，tj为转子转动惯量，d为阻尼系数系统的微分-代数方程组如下式：其中，i为发电机编号，i＝1,2,3,j为网络节点编号，j＝1,2,…,9潮流计算初始化：首先由潮流计算得到各发电机的端电压和各发电机注入网络中的功率进而计算得到发电机注入网络中的电流进而求得虚构电势即转子角度的初值δ(0)＝arctg(eqy(0)/eqx(0))转子转速的初值ω(0)＝1利用坐标变换公式，得到发电机定子电压和电流的d,q轴分量暂态电势的初值e′q(0)＝vq(0)+raiq(0)+x′did(0)原动机机械功率初值采用二级三阶半隐式龙格库塔法求解微分方程：所述二级三阶半隐式龙格库塔的积分格式如下式：其中，w1＝-0.41315432,w2＝1.41315432,a11＝c11＝0.17378663b1＝1.40824829,b2＝0.59175171由式(10)建立系统状态变量矩阵x和函数矩阵f建立全系统微分方程组的雅克比矩阵将式(12)和(13)代入到式(11)中，求解得到系数矩阵k1,k2，进而求得状态变量xn+1，即x(t+△t)；随后计算发电机注入网络的电流交替迭代求解网络方程yv＝i，得到y(t+△t).以(x(t+△t),y(t+△t))为下一时刻的初值，进行下一步长的计算。采用上述仿真计算流程，设置以下仿真场景：t＝0.5s时，在线路5-7靠近母线7处发生三相短路故障，t＝0.6时，故障清除，仿真总时长设为10s，并将本发明所述的算法与隐式梯形积分法进行对比，以辅助说明本发明所提方法的技术优势。仿真在配置intel(r)core(tm)i7-4550cpu@1.50ghz8g内存的笔记本电脑上进行。下表1给出了半隐式龙格库塔法和隐式梯形积分法仿真用时和仿真步长的对比。表1半隐式龙格库塔和隐式梯形积分法仿真对比所用算法仿真步长仿真用时仿真精度半隐式龙格库塔0.05s1.69s基本一致隐式梯形积分法0.01s5.32s基本一致图3和图4也给出了仿真图形结果对比，从图上和表1均可看出采用半隐式龙格库塔算法和隐式梯形积分法的仿真精度基本一致，半隐式龙格库塔算法的精度可以得到保障。除此之外，本发明还做了其他的仿真测试，增大隐式梯形积分法的仿真步长时，例如0.05s，仿真对比结果如表2所示。表2半隐式龙格库塔和隐式梯形积分法仿真对比所用算法仿真步长仿真用时仿真精度半隐式龙格库塔0.05s1.69s精确隐式梯形积分法0.05s/数值振荡下表中可以看出，当增大隐式梯形积分法的仿真步长时，仿真结果出现振荡，无法正常完成整个的暂态仿真计算，而半隐式龙格库塔算法的仿真结果依旧精确。仿真结果表明，本发明提出的半隐式龙格库塔算法能够使用大步长，数值稳定性好，计算精度更高，计算效率更快，能够更加高效精确稳定的进行电力系统暂态稳定计算。以上显示和描述了本发明的基本原理、主要特征和本发明的优点。本行业的技术人员应该了解，本发明不受上述实例的限制，上述实例和说明书中描述的只是说明本发明的原理，在不脱离本发明精神和范围的前提下本发明还会有各种变化和改进，这些变化和改进都落入要求保护的本发明范围内。本发明要求保护范围由所附的权利要求书及其等同物界定。当前第1页12