MATLAB实验电力系统暂态稳定分析

实验三 电力系统暂态稳定分析电力系统暂态稳定计算实际上就是求解发电机转子运动方程的初值问题，从而得出δ-t 和ω-t 的关系曲线。

每台发电机的转子运动方程是两个一阶非线性的常微分方程。

因此，首先介绍常微分方程的初值问题的数值解法。

一、常微分方程的初值问题 （一）问题及求解公式的构造方法我们讨论形如式（3-1）的一阶微分方程的初值问题⎩⎨⎧=≤≤='00)(),,()(y x y bx a y x f x y （3-1） 设初值问题（3-1）的解为)(x y ，为了求其数值解而采取离散化方法，在求解区间[b a ,]上取一组节点b x x x x x a n i i =<<<<<<=+ 110称i i i x x h -=+1（1,,1,0-=n i ）为步长。

在等步长的情况下，步长为nab h -=用i y 表示在节点i x 处解的准确值)(i x y 的近似值。

设法构造序列{}i y 所满足的一个方程（称为差分方程）),,(1h y x h y y i i i i ϕ⋅+=+ （3-2）作为求解公式，这是一个递推公式，从（0x ，0y ）出发，采用步进方式，自左相右逐步算出)(x y 在所有节点i x 上的近似值i y （n i ,,2,1 =）。

在公式（3-2）中，为求1+i y 只用到前面一步的值i y ，这种方法称为单步法。

在公式（3-2）中的1+i y 由i y 明显表示出，称为显式公式。

而形如（3-3）),,,(11h y y x h y y i i i i i ++⋅+=ψ （3-3）的公式称为隐式公式，因为其右端ψ中还包括1+i y 。

如果由公式求1+i y 时，不止用到前一个节点的值，则称为多步法。

由式（3-1）可得dy =dx y x f ),( （3-4）两边在[i x ，1+i x ]上积分，得⎰++=+1))(,()()(1i ix x i i dx x y x f x y x y （3-5）由此可以看出，如果想构造求解公式，就要对右端的积分项作某种数值处理。

这种求解公式的构造方法叫做数值积分法。

（二）一般的初值问题的解法1． 欧拉法和改进欧拉法对于初值问题（3-1），采用数值积分法，从而得到（3-5）。

对于（3-5）右端的积分用矩形公式(取左端点)，则得到⎰+⋅≈1))(,())(,(i ix x i i x y x f h dx x y x f进而得到（3-1）的求解公式（3-2）),(1i i i i y x f h y y ⋅+=+ （i =0，1，2，n-1） （3-6）此公式称为欧拉（Euler ）格式。

如果对式（3-5）右端的积分用梯形公式)))(,())(,((2))(,(111++⎰++⋅≈i i x x i i x y x f x y x f hdx x y x f i i则可以得到初值问题（3-1）的梯形求解公式如式（3-7）[]),(),(2111=+++⋅+=i i i i i i y x f y x f hy y （i =0，1，2，n-1） （3-7） 式（3-7）是个隐式公式。

可以采取先用欧拉格式求一个)(1+i x y 的初步近似值，记作1+i y ，称之为预报值，然后用预报值1+i y 替代式（3-7）右端的1+i y ，再计算得到1+i y ，称之为校正值，这样建立起来的预报－校正方法称为改进欧拉格式[]⎪⎩⎪⎨⎧+⋅+=⋅+=++++),(),(2),(1111i i i i i i i i i i y x f y x f h y y y x f h y y （3-8）2． 龙格—库塔方法在单步法中，应用最广泛的是龙格－库塔（Runge-kutta ）法，简称R －K 法。

下面直接给出一种四阶的龙格－库塔法的计算公式（3-9）⎪⎪⎪⎪⎩⎪⎪⎪⎪⎨⎧++⋅=++⋅=++⋅=⋅=++++==),()21,2()21,2(),()22(61342312143211K y h x f h K K y h x f h K K y h x f h K y x f h K K K K K y y i i i i i i i i i i （3-9） 它也称为标准（古典）龙格－库塔法。

例3-1 研究下列微分方程的初值问题⎪⎩⎪⎨⎧=-+='0)0(21122y y xy 解：这是一个特殊的微分方程，其解的解析式可以给出，为21x xy +=应用龙格－库塔法，取h =0.25，根据式（3-9）编写一段程序，由零开始自左相右逐步算出)(x y 在所有节点i x 上的近似值i y 。

计算结果见表3-1。

计算结果表明，四阶龙格－库塔方法的精度是较高的。

表3-1实际上，MATLAB 为常微分方程提供了很好的解题指令，使得求解常微分方程变得很容易，并且能将问题及解答表现在图形上。

因此，我们可以不用根据式（3-9）编写较复杂的程序，而只需应用MATLAB 提供的常微分方程解题器来解决问题。

下面给出用MATLAB 编写的解题程序。

首先编写描述常微分方程的ODE 文件，文件名为ˊmyfun ˊ，便于解题器调用它。

function dy = myfun(x,y) dy = zeros(1,1); dy=1/(1+x^2)-2*y^2;再编写利用解题器指令求解y 的程序。

clear x0=0; for i=1:4 xm=2*i;y0=0;[x,y] = ode45('myfun',[x0 xm],[y0]); format long y(length(y)) endplot(x,y,'-')运行上述程序，在得到几个点的函数值的同时，也得到函数y 的曲线，如图3-1所示。

图3-1 根据运算结果画出y 的曲线二、简单电力系统的暂态稳定性 （一）物理过程分析某简单电力系统如图3-2(a)所示，正常运行时发电机经过变压器和双回线路向无限大系统供电。

发电机用电势E' 作为其等值电势，则电势E ' 与无限大系统间的电抗为 212T LT dx x x x x +++'=I （3-10） 这时发电机发出的电磁功率可表示为δδsin sin M P x UE P I II ='=（3-11） 如果突然在一回输电线路始端发生不对称短路，如图3-2(b)所示。

故障期间发电机电势E' 与无限大系统之间的联系电抗为 ∆II ++'++++'=x x x x x x x x x x T LT dT LT d)2)(()2()(2121 （3-12）在故障情况下发电机输出的电磁功率为δδsin sin M P x UE P II IIII ='=（3-13） 在短路故障发生之后，线路继电保护装置将迅速断开故障线路两端的断路器，如图3-2(c)所示。

此时发电机电势E' 与无限大系统间的联系电抗为21T L T dx x x x x +++'=III （3-14） 发电机输出的功率为δδsin sin M P x UE P III IIIIII ='=（3-15） U =cE~G T 1T 2L~jx jx jx x j 'UE' (c )Ljx 1T jx 2T jx dx j 'U图3-2 简单电力系统及其等值电路（a ）正常运行方式及其等值电路；（b ）故障情况及其等值电路；（c ）故障切除后及其等值电路如果正常时发电机向无限大系统输送的有功功率为0P ，则原动机输出的机械功率T P 等于0P 。

假定不计故障后几秒种之内调速器的作用，即认为机械功率始终保持0P 。

因此，可以得到此简单电力系统正常运行、故障期间及故障切除后的功率特性曲线如图3-3所示。

0δk cm h 0P P T =Pδ图3-3 简单系统正常运行、故障期间及故障切除后的功率特性曲线对于上述简单电力系统，我们可以根据等面积定则求得极限切除角。

但是，实际工作需要知道在多少时间之内切除故障线路，也就是要知道与极限切除角对应的极限切除时间。

要解决这个问题，必须求解发电机的转子运动方程。

（二）求解发电机的转子运动方程求解发电机转子运动方程可以得出δ-t 和ω-t 的关系曲线。

其中δ-t 曲线一般称为摇摆曲线。

在上述简单电力系统中故障期间的转子运动方程为⎪⎪⎩⎪⎪⎨⎧-=-=II )sin (1)1(1δωωωδM T J P P T dt d dtd （3-16） 式中，δ——功率角，其单位为弧度；ω——转子角速度，标幺值；1ω——转子的同步角速度，即1ω=f π2=314.16，其单位为弧度/秒；J T ——发电机的惯性时间常数，其单位为秒；T P 、M P II ――分别为机械和电磁功率，标幺值。

这是两个一阶的非线性常微分方程，它的起始条件是已知的，即t =0t =0； ω=0ω=1.0；δ=0δ=MTP P I -1sin 故障切除后，由于系统参数改变，以致发电机功率特性发生变化，必须开始求解另一组微分方程：⎪⎪⎩⎪⎪⎨⎧-=-=III )sin (1)1(1δωωωδM T J P P T dt d dtd （3-17） 式中变量含义同前述，其中M P III 也为标幺值。

这组方程的起始条件为t =c t ；δ=c δ；ω=c ω其中c t 为给定的切除时间；c δ、c ω为与c t 时刻对应的δ和ω，它们可由故障期间的δ-t 和ω-t 的关系曲线求得（δ和ω都是不突变的）。

一般来说，在计算故障发生后几秒种的过程中，如果δ始终不超过180º，而且振荡幅值越来越小，则系统是暂态稳定的。

当发电机与无限大系统之间发生振荡或失去同步时，在发电机的转子回路中，特别是阻尼绕组中将有感应电流而形成阻尼转矩（也称为异步转矩）。

当作微小振荡时，阻尼功率可表达为：D P =ω∆D =)1(-ωD （3-18）式中，D 称为阻尼功率系数；ω∆为转子角速度的偏移量，标幺值；ω为转子角速度，标幺值。

阻尼功率系数D 除了与发电机的参数有关外，还和原始功角、δ∆的振荡频率有关。

在一般情况下它是正数。

在原始功角较小，或者定子回路中有串联电容使定子回路总电阻相对于总电抗较大时，D 可能为负数。

如果考虑阻尼功率的影响，则故障后的转子运动方程又可表达为⎪⎪⎩⎪⎪⎨⎧---=-=III ]sin )1([1)1(1δωωωωδM T J P D P T dt d dtd （3-19） 电力系统暂态稳定计算包括两类问题，一类是应用数值计算法得出故障期间的曲线后，根据曲线找到与极限切除角对应的极限切除时间，此时只需要求解微分方程（3-16）；另一类是已知故障切除时间，需要求出摇摆曲线来判断系统的稳定性，此时需要分段分别求解微分方程（3-16）和（3-17）。