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

1、电力系统仿真作业-三机九节点电力系统暂态仿真学院：能源与动力工程学院专业：电力系统及其自动化学号：姓名：于永生导师：授课教师：目录一、概述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于永生电力系统仿真作业一、 概述在动态稳定分析中，系统由线

2、性化的微分方程组和代数方程组描写，并用经典的或现代的线性系统理论来进行稳定分析，分析可以在时域或频域进行。当用计算机和现代线性系统理论分析时，常把系统线性化的微分方程组和代数方程组消去代数变量，化为状态方程形式，并广泛采用特征分析进行稳定分析。电力系统是由不同类型的发电机组、多种电力负荷、不同电压等级的电力网络等组成的十分庞大复杂的动力学系统。其暂态过渡过程不仅包括电磁方面的过渡过程，而且还有机电方面的过渡过程。由此可见，电力系统的数学模型是一个强非线性的高维状态方程组。在动态稳定仿真中使用简单的电力系统模型，发电机用三阶模型表示。二、 课程主要任务本次课程主要应用P. M. Anderson

3、 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

4、 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

5、 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

6、为基准值的标幺值）2. 潮流计算首先进行潮流计算，采用牛顿拉夫逊迭代法，电力系统潮流计算是电力系统运行和规划中最基本和最经常的计算，其任务是在已知某些运行参数的情况下，计算出系统中全部的运行参数，一般来说，各个母线所供负荷的功率是已知的，各个节点电压是未知的（平衡节点除外），可以根据网络结构形成节点导纳矩阵，然后由节点导纳矩阵和网络拓扑结构列写功率方程，由于功率方程里功率是已知的，电压的幅值和相角是未知的，这样潮流计算的问题就转化为求解非线性方程组的问题了。为了便于用迭代法解方程组，需要将上述功率方程改写成功率平衡方程，并对功率平衡方程求偏导，得出对应的雅可比矩阵，给未知节点赋电压初值，一般为

7、额定电压，将初值带入功率平衡方程，得到功率不平衡量，这样由功率不平衡量、雅可比矩阵、节点电压不平衡量（未知的）构成了误差方程，解误差方程，得到节点电压不平衡量，节点电压加上节点电压不平衡量构成新的节点电压初值，将新的初值带入原来的功率平衡方程，并重新形成雅可比矩阵，然后计算新的电压不平衡量，这样不断迭代，不断修正，一般迭代三到五次就能收敛。牛顿拉夫逊算法修正方程W = -JV其中W是节点不平衡量向量，包括有功，无功，电压；J是雅克比矩阵；V是节点电压修正量。令，则直角坐标形式的功率不平衡量方程为PQ节点：PV节点：极坐标形式的功率不平衡量方程雅可比矩阵J各元素的表达式当ji时：当j=i时：其中

8、，。进行牛顿拉夫逊算法迭代后得到电压幅值V和相角。3. 负荷等效和支路简化然后求出支路电流，将发电机内电抗X加入系统导纳矩阵，求出发电机内电势E。加入发电机内节点后，系统导纳矩阵变成12*12阶的矩阵，并将负荷等效成阻抗。然后将支路导纳矩阵分块，如下：其中，A是3*3的方阵，E是3*9的矩阵，C是9*3的矩阵，D是9*9的方阵。经过网络简化得到故障前的3*3简化导纳矩阵 （1）其中“inv（D）”是MATLAB中D矩阵的求逆。故障中导纳矩阵的第七行和第七列从矩阵中删除，此时有此时，A是3*3的方阵，E是3*8的矩阵，C是8*3的矩阵，D是8*8的方阵。简化矩阵求法如同公式（1）。故障后的Y矩阵

9、相对于故障前的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.

10、 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号发电

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

12、三个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

13、 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=

14、L(:,2);%末节点编号r=L(:,3);%电阻x=L(:,4);%电抗b=L(:,5);%电纳a=L(:,6);%变压器非标准变比z=r + i*x;%z矩阵，复数，9乘1的矩阵y=1./z;%求倒数，9乘1的导纳矩阵b=i*b;%虚数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）

15、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%节点数据% 节点 电压 电压 发电机

16、 发电机 负荷 负荷 节点% 号 幅值 相角 有功 无功 有功 无功 类型(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);%额定电压de

17、l=N(:,3);%电压相角Pg=N(:,4);%发电机有功功率Qg=N(:,5);%发电机无功功率Pl=N(:,6);%负荷有功功率Ql=N(:,7);%负荷无功功率P=Pg-Pl;%发出减消耗Q=Qg-Ql;%发出减消耗Psp=P;%额定有功功率Qsp=Q;%额定无功功率G=real(Y);B=imag(Y);pv=find(type=3|type=2);%pv节点编号和平衡节点编号pq=find(type=1);%pq节点编号npv=length(pv);%npv=3npq=length(pq);%npq=6Tol=1;Iter=1;while(Tol>1e-7) P=zeros(

18、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=Psp-P;%额定有功功率减去流出或流入的有功功率 dQa=Qsp-Q;%额定无功功率减去流出或流入的无功功率 k=1; dQ=zeros(npq,1)

19、;%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(nbus-1,nbus-1);%有功功率对功角求偏微分 for i=1:(nbus-1) m=i+1; for k=1:(nbus-1) 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)-de

20、l(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(nbus-1,npq);%有功功率对电压求偏微分 for i=1:(nbus-1) 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); e

21、nd 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,nbus-1);%无功功率对相角求偏微分 for i=1:npq m=pq(i); for k=1:(nbus-1) 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,

22、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

23、(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(

24、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

25、 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 j*0.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); ende

26、ndYzengjia(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);%

27、求负荷等效导纳Y5=N(5,6)/(V(5)2-j*(N(5,7)/(V(5)2);Y6=N(6,6)/(V(6)2-j*(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

28、:12 E(m,n-3)=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(m-3,n-3)=Yzengjia(m,n); endendYpre=A-E*(inv(D)*C;%得到故障前的3*3矩阵Yzengjiadur=Yzengjia;Yzengjiadur(:,10)=;Yzengjiadur(10,:)=;Adur=A;Edur=zeros(3,8);for m=1:3 for n=1:8 Ed

29、ur(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矩阵Yzengjiaa

30、ft=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,n-3)=Yzengjiaaft(m,n); endend

31、Caft=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=Aaft-Eaft*(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)=V

32、m(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

33、(2);Tj=47.28;12.8;6.02;Pm0=Pe;%机械功率=故障前的电磁功率X0=ones(3,1).*(2*pi*60),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

34、=t1;t2;XOUT=XOUT*180/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) + j*rho.*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_f