三机九节点电力系统仿真matlab

电力系统仿真作业------------三机九节点电力系统暂态仿真学院：能源与动力工程学院专业：电力系统及其自动化学号：姓名：于永生导师：授课教师：1 / 19目录一、概述 (1)二、课程主要任务 (1)1.系统数据 (1)2.潮流计算 (2)3.负荷等效和支路简化 (4)4.求解电磁功率 (5)5.求解运动方程 (5)6.程序清单 (7)(1).主程序： (7)(2).极坐标转换成直角坐标函数pol2rect(V,del) (16)(3).直角坐标转换成极坐标函数rect2pol(Z) (16)(4).求解微分方程所用的得到微分量的函数Gen_fw(t,X,Y_Gen,E,Pm0,Tj) (16)三、课程总结及心得体会 (16)四、参考文献 (17)一、概述在动态稳定分析中，系统由线性化的微分方程组和代数方程组描写，并用经典的或现代的线性系统理论来进行稳定分析，分析可以在时域或频域进行。

当用计算机和现代线性系统理论分析时，常把系统线性化的微分方程组和代数方程组消去代数变量，化为状态方程形式，并广泛采用特征分析进行稳定分析。

电力系统是由不同类型的发电机组、多种电力负荷、不同电压等级的电力网络等组成的十分庞大复杂的动力学系统。

其暂态过渡过程不仅包括电磁方面的过渡过程，而且还有机电方面的过渡过程。

由此可见，电力系统的数学模型是一个强非线性的高维状态方程组。

在动态稳定仿真中使用简单的电力系统模型，发电机用三阶模型表示。

二、课程主要任务本次课程主要应用P. M. Anderson and A. A. Fouad编写的《Power System Control and Stability》一书中所引用的Western System Coordinated Council (WSCC)三机九节点系统模型。

1.系统数据其中，节点数据如下：%节点数据% 节点电压电压发电机发电机负荷负荷节点% 号幅值相角有功无功有功无功类型(1PQ 2PV 3平衡)N=[ 1 1.04 0 0.7164 0.2705 0 0 32 1.025 0 1.63 0.0665 0 0 23 1.025 0 0.85 -0.1086 0 0 24 1 0 0 0 0 0 15 1 0 0 0 1.25 0.5 16 1 0 0 0 0.9 0.3 17 1 0 0 0 0 0 18 1 0 0 0 1 0.35 19 1 0 0 0 0 0 1];其中，支路数据如下：% 线路数据% 首端末端电阻电抗电纳（1/2) 变压器非标准变比L=[4 5 0.01 0.085 0.088 14 6 0.017 0.092 0.079 15 7 0.032 0.161 0.153 16 9 0.039 0.17 0.179 17 8 0.0085 0.072 0.0745 11 / 198 9 0.0119 0.1008 0.1054 11 4 0 0.0576 0 12 7 0 0.0625 0 13 9 0 0.0586 0 1];发电机数据如下：% 发电机母线Xd Xd' Td0' Xq Xq' Tq0’Tj XfGe=[ 1 1 0.1460 0.0608 8.96 0.0969 0.0969 0 47.28 0.05762 2 0.8958 0.1198 6.00 0.8645 0.1969 0.535 12.80 0.06253 3 1.3125 0.1813 8.59 1.2578 0.2500 0.600 6.02 0.0585];系统电路结构拓扑图如下：图1 WSCC 3机9节点系统（所有参数以100MV A为基准值的标幺值）2.潮流计算首先进行潮流计算，采用牛顿拉夫逊迭代法，电力系统潮流计算是电力系统运行和规划中最基本和最经常的计算，其任务是在已知某些运行参数的情况下，计算出系统中全部3 / 19的运行参数，一般来说，各个母线所供负荷的功率是已知的，各个节点电压是未知的（平衡节点除外），可以根据网络结构形成节点导纳矩阵，然后由节点导纳矩阵和网络拓扑结构列写功率方程，由于功率方程里功率是已知的，电压的幅值和相角是未知的，这样潮流计算的问题就转化为求解非线性方程组的问题了。

为了便于用迭代法解方程组，需要将上述功率方程改写成功率平衡方程，并对功率平衡方程求偏导，得出对应的雅可比矩阵，给未知节点赋电压初值，一般为额定电压，将初值带入功率平衡方程，得到功率不平衡量，这样由功率不平衡量、雅可比矩阵、节点电压不平衡量（未知的）构成了误差方程，解误差方程，得到节点电压不平衡量，节点电压加上节点电压不平衡量构成新的节点电压初值，将新的初值带入原来的功率平衡方程，并重新形成雅可比矩阵，然后计算新的电压不平衡量，这样不断迭代，不断修正，一般迭代三到五次就能收敛。

牛顿拉夫逊算法修正方程W = -J ΔV其中W 是节点不平衡量向量，包括有功，无功，电压；J 是雅克比矩阵；ΔV 是节点电压修正量。

令ijij ij i i i jB G Y jf e V +=+=;，则直角坐标形式的功率不平衡量方程为 PQ 节点：)()(11=+---=∆∑∑==n j nj j ij j ij j j ij j ij i is i e B f G f f B e G e P P)()(11=++--=∆∑∑==nj nj j ij j ij j j ij j ij i is i e B f G e f B e G f Q QPV 节点：)()(11=+---=∆∑∑==n j nj j ij j ij j j ij j ij i is i e B f G f f B e G e P P)(222222i i is i is i f e V V V V +-=-=∆极坐标形式的功率不平衡量方程)sin cos (1=+-=∆∑=nj ij ij ij ij j i is i B G V V P P δδ)cos sin (1=--=∆∑=nj ij ij ij ij j i is i B G V V Q Q δδ雅可比矩阵J 各元素的表达式⎭⎬⎫⎩⎨⎧=L MN H J当j ≠i 时：iij ij i ij ijP H B e G f f ∂∆==-∂iij ij i ij ijP N G e B f e ∂∆==--∂iij ij i ij ijQ M B f G e f ∂∆==+∂iij ij i ij ijQ L G f B e e ∂∆==-+∂当j=i 时：iii ii i ii i iiP H B e G f b f ∂∆==--∂ iii ii i ii i iiP N G e B f a e ∂∆==---∂ iii ii i ii i iiQ M G e B f a f ∂∆==+-∂iii ii i ii i iiQ L B e G f b e ∂∆==-+∂其中，∑∑==+=-=ni j ij j ij i ni j ij j ij i e B f G b f B e G a 11)()(。

进行牛顿拉夫逊算法迭代后得到电压幅值V 和相角θ。

3. 负荷等效和支路简化然后求出支路电流，将发电机内电抗X ’加入系统导纳矩阵，求出发电机内电势E ’。

加入发电机内节点后，系统导纳矩阵变成12*12阶的矩阵，并将负荷等效成阻抗。

然后将支路导纳矩阵分块，如下：A E Y C D ⎛⎫= ⎪⎝⎭其中，A 是3*3的方阵，E 是3*9的矩阵，C 是9*3的矩阵，D 是9*9的方阵。

经过网络简化得到故障前的3*3简化导纳矩阵5 / 19Ypre=A-E*(inv(D))*C （1）其中“inv （D ）”是MA TLAB 中D 矩阵的求逆。

故障中导纳矩阵的第七行和第七列从矩阵中删除，此时有A E Y C D ⎛⎫= ⎪⎝⎭此时，A 是3*3的方阵，E 是3*8的矩阵，C 是8*3的矩阵，D 是8*8的方阵。

简化矩阵求法如同公式（1）。

故障后的Y 矩阵相对于故障前的Y 矩阵只是第5个节点和第7个节点有变化，反映到12*12的矩阵中即为第（10,8），（8,10），（8,8），（10,10）位置的元素有变化，其中（10,8），（8,10）位置的元素变为零，(8,8)，（10,10）节点在故障前的基础上加上（8,10）位置元素的值。

然后简化导纳矩阵的求法同式（1）。

4. 求解电磁功率得到故障前，故障中，故障后三个不同的导纳矩阵后，就开始计算电磁功率和机械功率，机械功率等于稳态的电磁功率中的有功分量。

所以可以有Pe=real(E*I) (2)公式（2）中，E 为发电机内电势，I 为从发电机流出的电流。

但在参考文献Ramnarayan Patel, T. S. Bhatti and D. P. Kothari.MA TLAB/Simulink-based transient stability analysis of a multimachine power system 中给出的电磁功率计算公式为：∑≠=+-+=nij j j i ij ij j i ii i ei Y E E G E P 12)cos(δδθ本人在此次仿真中用的是公式（2）计算得到的电磁功率。

稳态情况下有，机械功率Pme=Pe5. 求解运动方程发电机的运动方程可以写成常微分方程组：⎪⎪⎪⎩⎪⎪⎪⎨⎧-=⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡+-+-=+∑≠=R i n i j j j i ij ij j i ii i mi j j i R i dtd Y E E G E P D dt d H ωωδδδθωωω12)cos(2其中Pmi 为第i 个机组故障前稳态的电磁功率。

在本次仿真中D j ωi 为零，即阻尼为零。

仿真开始，t=0时引入故障，0.083s 后切除故障。

求解运动方程后得到ω和δ曲线如下：2号发电机3号发电机图1 ω曲线2号发电机3号发电机1号发电机图2 δ曲线6.程序清单以下是我编写的仿真程序，除主程序外还包含三个函数：①极坐标转换成直角坐标函数pol2rect(V,del),其中输入参数V为极坐标向量的幅值，del为极坐标向量的角度，函数返回值为一个复数，即为极坐标向量在直角坐标中的复数值；②直角坐标转换成极坐标函数rect2pol(Z),其中输入参数Z为一个复数，函数返回值为一个二维向量，向量的第一个数为幅值，第二个数为相角；③求解微分方程所用的得到微分量的函数Gen_fw(t,X,Y_Gen,E,Pm0,Tj),其中输入参数X为ω和δ的迭代初值，Y_Gen，为求解所用的导纳矩阵，这里是三阶的方阵，对应于故障前，故障中，和故障后的三个Y矩阵，E为发电机内电势，Pm0为机械功率，等于稳态时的有功功率，Tj为运动方程中的2H。