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

电力系统仿真作业-三机九节点电力系统暂态仿真学院:能源与动力工程学院专业：电力系统及其自动化学号：姓名：于永生导师:授课教师：目录一、概述1二、课程主要任务11.系统数据12.潮流计算23.负荷等效和支路简化44.求解电磁功率55.求解运动方程56.程序清单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。040 0.7164 0.2705 0 0 3 2 1。0250 1。63 0.0665 0 0 2 3 1。0250 0。85 0。1086 0 0 2 4 1 0 0 0 0 0 1 5 1 0 0 0 1.25 0。5 1 6 1 0 0 0 0.9 0。3 1 7 1 0 0 0 0 0 1 8 1 0 0 0 1 0.35 1 9 1 0 0 0 0 0 1;其中,支路数据如下： 线路数据 首端 末端 电阻 电抗 电纳（1/2） 变压器非标准变比L=4 5 0。01 0。085 0.088 1 4 6 0。017 0。092 0.079 1 5 7 0。032 0。161 0。153 1 6 9 0。039 0。17 0。179 1 7 8 0.0085 0。072 0。0745 1 8 9 0.0119 0.1008 0.1054 1 1 4 0 0。0576 0 1 2 7 0 0.0625 0 1 3 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。0576 2 2 0。8958 0.1198 6.00 0。8645 0。1969 0.535 12.80 0.0625 3 3 1.3125 0。1813 8.59 1。2578 0.2500 0.600 6。02 0.0585；系统电路结构拓扑图如下:图1 WSCC 3机9节点系统（所有参数以100MVA为基准值的标幺值）2. 潮流计算首先进行潮流计算，采用牛顿拉夫逊迭代法，电力系统潮流计算是电力系统运行和规划中最基本和最经常的计算，其任务是在已知某些运行参数的情况下，计算出系统中全部的运行参数，一般来说，各个母线所供负荷的功率是已知的,各个节点电压是未知的（平衡节点除外),可以根据网络结构形成节点导纳矩阵，然后由节点导纳矩阵和网络拓扑结构列写功率方程，由于功率方程里功率是已知的，电压的幅值和相角是未知的，这样潮流计算的问题就转化为求解非线性方程组的问题了。为了便于用迭代法解方程组，需要将上述功率方程改写成功率平衡方程，并对功率平衡方程求偏导，得出对应的雅可比矩阵，给未知节点赋电压初值,一般为额定电压，将初值带入功率平衡方程,得到功率不平衡量，这样由功率不平衡量、雅可比矩阵、节点电压不平衡量（未知的）构成了误差方程，解误差方程，得到节点电压不平衡量，节点电压加上节点电压不平衡量构成新的节点电压初值,将新的初值带入原来的功率平衡方程,并重新形成雅可比矩阵,然后计算新的电压不平衡量，这样不断迭代，不断修正，一般迭代三到五次就能收敛.牛顿拉夫逊算法修正方程W = -JV其中W是节点不平衡量向量，包括有功，无功，电压；J是雅克比矩阵;V是节点电压修正量。令,则直角坐标形式的功率不平衡量方程为PQ节点：PV节点：极坐标形式的功率不平衡量方程雅可比矩阵J各元素的表达式当ji时:当j=i时:其中，。进行牛顿拉夫逊算法迭代后得到电压幅值V和相角。3. 负荷等效和支路简化然后求出支路电流，将发电机内电抗X加入系统导纳矩阵，求出发电机内电势E。加入发电机内节点后,系统导纳矩阵变成1212阶的矩阵，并将负荷等效成阻抗。然后将支路导纳矩阵分块，如下：其中,A是33的方阵，E是39的矩阵，C是93的矩阵,D是99的方阵。经过网络简化得到故障前的33简化导纳矩阵 （1）其中“inv(D）”是MATLAB中D矩阵的求逆.故障中导纳矩阵的第七行和第七列从矩阵中删除，此时有此时,A是3*3的方阵，E是38的矩阵，C是8*3的矩阵，D是8*8的方阵。简化矩阵求法如同公式(1）。故障后的Y矩阵相对于故障前的Y矩阵只是第5个节点和第7个节点有变化，反映到1212的矩阵中即为第（10,8），(8，10),（8，8），（10，10）位置的元素有变化,其中（10，8)，（8，10）位置的元素变为零，（8,8）,（10,10）节点在故障前的基础上加上（8,10)位置元素的值。然后简化导纳矩阵的求法同式(1）。4. 求解电磁功率得到故障前，故障中，故障后三个不同的导纳矩阵后，就开始计算电磁功率和机械功率，机械功率等于稳态的电磁功率中的有功分量。所以可以有Pe=real（EI) （2）公式（2）中，E为发电机内电势，I为从发电机流出的电流。但在参考文献Ramnarayan Patel, T. S。 Bhatti and D. P。 Kothari.MATLAB/Simulink-based transient stability analysis of a multimachine power system中给出的电磁功率计算公式为：本人在此次仿真中用的是公式（2）计算得到的电磁功率。稳态情况下有，机械功率Pme=Pe5. 求解运动方程发电机的运动方程可以写成常微分方程组:其中Pmi为第i个机组故障前稳态的电磁功率.在本次仿真中Dji为零，即阻尼为零。仿真开始，t=0时引入故障,0.083s后切除故障.求解运动方程后得到和曲线如下：3号发电机2号发电机图1 曲线1号发电机3号发电机2号发电机图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。(1). 主程序：clc；clear;%节点数据 节点 电压 电压 发电机 发电机 负荷 负荷 节点% 号 幅值 相角 有功 无功 有功 无功 类型(1PQ 2PV 3平衡 N= 1 1。040 0。71640.2705 0 0 3% 2 1。0250 1。63 0.0665 0 0 2% 3 1。0250 0。85 0.1086 0 0 2 4 1 0 0 0 0 0 1 5 1 0 0 0 1.250。5 1% 6 1 0 0 0 0.90.3 1 7 1 0 0 0 0 0 1 8 1 0 0 0 1 0.35 1 9 1 0 0 0 0 0 1; 线路数据% 首端 末端 电阻 电抗 电纳（1/2） 变压器非标准变比L=4 5 0.01 0.085 0.088 1 4 6 0。017 0.092 0。079 1 5 7 0.032 0。161 0.153 1 6 9 0。039 0.17 0.179 1 7 8 0。0085 0。072 0.0745 1 8 9 0。0119 0。1008 0。1054 1 1 4 0 0.0576 0 1 2 7 0 0.0625 0 1 3 9 0 0。0586 0 1；%形成节点导纳矩阵fb=L（：，1）;起始节点编号tb=L（：,2)；%末节点编号r=L（：，3)；电阻x=L（：，4);电抗b=L（:，5）；%电纳a=L(：,6);%变压器非标准变比z=r + ix；z矩阵,复数，9乘1的矩阵y=1。/z；求倒数，9乘1的导纳矩阵b=ib；%虚数nb=max(max(fb），max（tb）);nb=9nl=length（fb）；支路个数Y=zeros（nb，nb）；%Y是9乘9的支路矩阵for k=1:nl Y(fb(k），tb（k）)=Y(fb（k），tb（k)）y(k）/a(k）;%Y方阵的元素,非对角元（4,5)（4，6)（5,7）（6，9）（7，8)（8，9）（1，4）（2,7)（3，9） Y(tb(k），fb（k）)=Y（fb（k），tb（k)）;Y方阵的元素，上一句的对称部分(5,4）（6，4）（7，5)(9,6)（8，7）(9，8）(4，1）（7，2）（9，3）endfor m=1：nb for n=1：nl if fb(n)=m支路中的第n个数的首节点编号等于m Y（m，m）=Y（m,m）+y（n）/(a(n)2）+b（n);Y矩阵中的对角元 elseif tb(n）=m%支路中的第n个数的末节点编号等于m Y(m，m）=Y(m，m）+y（n)+b（n)；因为支路是首端或末端等于m,不能是首端编号等于末端编号 end endend节点数据 节点 电压 电压 发电机 发电机 负荷 负荷 节点% 号 幅值 相角 有功 无功 有功 无功 类型(1PQ 2PV 3平衡N= 1 1.040 0。71640.2705 0 0 3 2 1。0250 1。63 0.0665 0 0 2 3 1.0250 0.85 -0.1086 0 0 2 4 1 0 0 0 0 0 1 5 1 0 0 0 1。250.5 1 6 1 0 0 0 0。90。3 1 7 1 0 0 0 0 0 1 8 1 0 0 0 1 0.35 1 9 1 0 0 0 0 0 1；nbus=9;节点数bus=N（：，1)；节点编号type=N（：，8)；节点类型V=N（:,2）;额定电压del=N（:，3)；%电压相角Pg=N(：,4）；%发电机有功功率Qg=N(：,5)；发电机无功功率Pl=N（：，6）；负荷有功功率Ql=N（:，7）；负荷无功功率P=PgPl；发出减消耗Q=QgQl;发出减消耗Psp=P；额定有功功率Qsp=Q;额定无功功率G=real（Y）；B=imag（Y);pv=find(type=3type=2）;%pv节点编号和平衡节点编号pq=find（type=1)；pq节点编号npv=length（pv）；%npv=3npq=length（pq）;npq=6Tol=1；Iter=1；while（Tol1e-7） P=zeros（nbus，1）;迭代初值P为零 Q=zeros（nbus，1）；%迭代初值Q为零 for i=1：nbus for k=1：nbus P（i）=P(i）+V（i)*V（k）（G（i,k）cos（del(i）-del（k)）+B（i，k）*sin（del（i)del（k)）； Q（i）=Q（i)+V（i)V（k）（G(i，k)sin(del(i）-del(k）B（i,k）cos(del(i)del（k）; end end dPa=PspP;额定有功功率减去流出或流入的有功功率 dQa=QspQ;额定无功功率减去流出或流入的无功功率 k=1； dQ=zeros（npq,1)；pq节点的无功功率，pq乘1的矩阵，6乘1 for i=1:nbus if type(i）=1%pq节点 dQ(k，1）=dQa（i）; k=k+1; end end dP=dPa（2:nbus）； Z=dP;dQ； %下面计算雅克比矩阵 J1=zeros(nbus1，nbus1);有功功率对功角求偏微分 for i=1：（nbus-1） m=i+1; for k=1：(nbus1） n=k+1； if n=m for n=1：nbus J1(i,k）=J1(i，k)+V(m)V（n)*(-G（m,n）*sin(del(m）-del（n)+B（m，n）cos（del（m）del（n)）; end J1（i，k)=J1(i,k)V(m）2*B(m，m）； else J1（i，k）=V（m）V（n)（G(m，n）*sin（del（m）-del（n)）-B（m，n）*cos（del（m）del(n）)； end end end J2=zeros（nbus1，npq)；%有功功率对电压求偏微分 for i=1：（nbus1) m=i+1; for k=1：npq n=pq（k）； if n=m for n=1：nbus J2(i，k）=J2(i，k)+V（n)*(G（m，n）cos（del（m)del（n)）+B（m，n)sin（del(m）-del(n）; end J2（i,k)=J2（i，k)+V（m)G（m,m）； else J2(i，k）=V(m）（G（m，n)*cos（del(m）-del（n)）+B（m,n）*sin（del（m）del(n)）; end end end J3=zeros(npq，nbus1）；%无功功率对相角求偏微分 for i=1:npq m=pq(i)； for k=1:(nbus1） n=k+1; if n=m for n=1：nbus J3（i，k）=J3（i，k)+V（m)*V(n)(G（m，n）cos（del(m)del(n）)+B（m，n）*sin（del（m）-del（n)）; end J3（i，k)=J3（i，k)V(m)2*G(m，m）; else J3(i，k)=V（m）*V(n）（-G（m，n）cos（del(m)-del(n）-B（m,n)sin（del(m）del(n)）; end end end J4=zeros（npq，npq）；无功功率对电压求偏微分 for i=1：npq m=pq（i）; for k=1：npq n=pq（k)； if n=m for n=1：nbus J4（i，k)=J4（i，k）+V（n）（G(m，n)*sin(del（m）del(n）)B（m,n）*cos(del（m）del（n）; end J4(i，k）=J4（i，k）-V（m)*B(m，m)； else J4（i，k）=V（m)*(G（m，n）*sin(del（m）del（n)）B(m，n)*cos（del(m)-del（n)）； end end end J=J1 J2;J3 J4； X=inv（J)Z； X=JZ; dTh=X(1：nbus-1）; dV=X(nbus:end）； del（2：nbus)=dTh+del(2：nbus）； k=1; for i=2:nbus if type(i）=1 V（i）=V(i)+dV（i-3)； k=k+1; end end Iter=Iter+1； Tol=max（abs（Z）；endVm = pol2rect(V,del）； % 极坐标转换成直角坐标Del = 180/pi*del； 弧度转化成角度Iij = zeros（nb,nb）；%支路电流for m = 1:nl%每一条支路都过一遍 p = fb(m）； q = tb(m)； Iij（p，q） = （Vm（p) Vm（q）)Y(p,q）; Y（m,n） = y(m，n)。 Iij（q，p） = -Iij（p，q)；endIij = sparse（Iij)； % 线路数据 首端 末端 电阻 电抗 电纳(1/2) 变压器非标准变比 L=4 5 0。01 0.085 0.088 1 4 6 0。017 0.092 0。079 1 5 7 0.032 0。161 0。153 1 6 9 0。039 0。17 0。179 1 7 8 0.0085 0.072 0。0745 1 8 9 0。0119 0。1008 0。1054 1 1 4 0 0。0576 0 1 2 7 0 0.0625 0 1 3 9 0 0.0586 0 1；Yzengjia=zeros(12，12）；Zci=j*0。0608 j0。1198 j*0.1813;Yci=1./Zci；for k=1：3 Yzengjia（k，k）=Yci（k）;endfor k=1：9 for m=1：9 Yzengjia（3+k，3+m）=Y（k，m)； endendYzengjia（1，4)=Yzengjia(1，1)；Yzengjia（4,1）=Yzengjia(1，4);Yzengjia(2，5)=Yzengjia（2，2）；Yzengjia(5,2）=Yzengjia(2,5）；Yzengjia（3,6）=Yzengjia（3，3)；Yzengjia（6，3）=Yzengjia（3,6）；Yzengjia（4，4)=Yzengjia（4，4）+Yzengjia（1,1）；Yzengjia（5，5)=Yzengjia（5，5）+Yzengjia（2,2）;Yzengjia（6，6）=Yzengjia（6,6)+Yzengjia（3，3）;求负荷等效导纳Y5=N(5,6）/(V（5)2-j（N（5，7)/(V（5）)2);Y6=N(6,6）/(V（6）2j（N(6,7）/(V（5）2）；Y8=N(8，6）/(V(8）2-j*(N（8,7）/（V（5）)2)；Yzengjia（8，8)=Yzengjia（8,8)+Y5;Yzengjia（9，9)=Yzengjia（9，9）+Y6;Yzengjia（11，11）=Yzengjia(11,11）+Y8;A=zeros（3，3);for m=1：3 for n=1：3 A（m,n)=Yzengjia(m，n）； endendE=zeros(3,9)；for m=1：3 for n=4：12 E(m,n3)=Yzengjia(m，n）; endendC=zeros（9,3）;for m=4：12 for n=1：3 C（m-3，n）=Yzengjia(m，n）； endendD=zeros（9，9）；for m=4：12 for n=4：12 D(m3，n3）=Yzengjia(m,n)； endendYpre=AE*（inv（D)*C;%得到故障前的33矩阵Yzengjiadur=Yzengjia;Yzengjiadur(：,10）=；Yzengjiadur（10，：)=;Adur=A；Edur=zeros（3,8)；for m=1：3 for n=1:8 Edur(m,n)=Yzengjiadur（m，n+3）;因为Edur的第四到第九列是零，而所以这里为了方便，直接取E的前六列 endendCdur=zeros(8，3）;for m=1：8 for n=1:3 Cdur（m,n）=Yzengjiadur（m+3,n）；%因为Cdur的第四到第九行是零,而所以这里为了方便，直接取C的前六行 endendDdur=zeros(8，8）；for m=1:8 for n=1：8 Ddur(m，n）=Yzengjiadur（m+3，n+3）； endendYdur=Adur-Edur*（inv（Ddur）)Cdur;%故障中的3*3矩阵Yzengjiaaft=Yzengjia；temp=Yzengjiaaft(8,10）；Yzengjiaaft(8，10）=0；Yzengjiaaft(8，8）=Yzengjiaaft（8，8)+temp;Yzengjiaaft(10,8）=0;Yzengjiaaft（10,10）=Yzengjiaaft（10,10）+temp；Aaft=zeros(3,3）;for m=1：3 for n=1:3 Aaft(m,n）=Yzengjiaaft（m，n)； endendEaft=zeros(3，9)；for m=1:3 for n=4：12 Eaft（m，n3）=Yzengjiaaft（m，n); endendCaft=zeros（9,3）;for m=4：12 for n=1：3 Caft（m-3，n）=Yzengjiaaft（m，n)； endendDaft=zeros(9,9）;for m=1：9 for n=1：9 Daft（m,n)=Yzengjiaaft（m+3,n+3)； endendYaft=AaftEaft（inv（Daft）Caft；故障后的3*3矩阵求电磁功率Enei=ones（3，1);%发电机内电势Pe=ones（3，1）;Enei（2)=Vm（2）+Zci（2）Iij（2，7）;Enei（1）=Vm（1)+Zci(1）(Iij(4，1)）；Enei(3）=Vm(3）+Zci(3）*Iij（3，9）；polEnei1=rect2pol(Enei(1)；直角坐标转换成极坐标后得到角度polEnei2=rect2pol（Enei（2）；polEnei3=rect2pol(Enei(3)）;Pe(1)=real（Enei(1）(Iij（4，1）)）；%故障前的电气功率Pe（2)=real（Enei(2）*(Iij（2，7）；Pe（3）=real（Enei（3)*(Iij(3，9)）)；delta=zeros（3，1）；delta(1）=polEnei1（2）；delta(2）=polEnei2（2)；delta（3）=polEnei3（2）；Tj=47。28;12。8;6。02；Pm0=Pe;机械功率=故障前的电磁功率X0=ones(3，1）.*（2pi60），delta；Y_Gen=Ydur；E=abs(Enei)；Time_s=0,0。083；0。083s后切除故障t1,Xout1=ode45(（t,X） Gen_fw（t，X，Y_Gen，E,Pm0，Tj)，Time_s,X0);Time_s=0。083，2；Y_Gen=Yaft;X0=Xout1（end，：)；t2,Xout2=ode45(（t,X） Gen_fw（t，X,Y_Gen，E，Pm0,Tj)，Time_s，X0）;XOUT=Xout1；Xout2；T=t1；t2；XOUT=XOUT180/pi；figure，plot(T，XOUT（：，4）,r,T,XOUT（:,5），g，T,XOUT（：，6），y)；figure，plot(T，XOUT（：,5）XOUT（：，4），r,T,XOUT（：，6）-XOUT（：，4)，Y)；(2). 极坐标转换成直角坐标函数pol2rect（V，del）function rect = pol2rect（rho，theta）rect = rho。cos（theta） + jrho.sin(theta)；(3). 直角坐标转换成极坐标函数rect2pol（Z)function pol=rect2pol(Z）t=real(Z);k=imag（Z)；pol（1）=(t2+k2）（1/2）；模pol（2）=atan（k/t）；(4). 求解微分方程所用的得到微分量的函数Gen_fw（t,X,Y_Gen,E，Pm0，Tj）function Xp = Gen_fw（t,X，Y_Gen，E,Pm0,Tj)Wn=ones（3,1)。2*pi60；D=zeros（3,1）;dW=zeros（3，1）；dDeta=zeros(3,1)；W=X（1：3）;Deta=X(4：6）；Gen_U=E.exp（1iDeta）;Pe=real（Gen_U.conj(Y_Gen*Gen_U)）；dW=（Pm0-Pe-D。W）。/Tj.*Wn；dDeta=WWn;Xp=dW；dDeta；End三、 课程总结及心得体会首先感谢任课老师对我仿真的指导与帮助，然后感谢我的导师对我的指导与关怀，还要感谢在程序编写过程中给过我无私帮助的其他同学。虽然我在以前也编写过其他的程序,但这次课程的仿真设计实验让我学到了很多。首先,对于错误的查找，举个简单的例子，之前我的雅克比矩阵一直有问题，但找了好久没有找到，后来发现是我的节点类型与老师定义的节点类型不一样,将后面的自己的节点类型检索改过后雅克比矩阵就成功了，瞬时感觉豁然开朗.再次,对于某一个参数的计算可能有不同的方法,例如对于电磁功率如果按照文献Ramnarayan Patel， T. S. Bhatti and D。 P. Kothari.MATLAB/Simulink-based transient stability analysis of a multimachine power