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

电力系统仿真作业-三机九节点电力系统 暂态仿真 ）学院：能源与动力工程学院 专业：电力系统及其自动化 学号：姓名：于永生 导师：授课教师：目录 一、概述.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 0 0 0 3 2 0 0 0 2 3 0 0 0 2 4 1 0 0 0 0 0 1 5 1 0 0 0 1 6 1 0 0 0 1 7 1 0 0 0 0 0 1 8 1 0 0 0 1 1 9 1 0 0 0 0 0 1;其中，支路数据如下：%线路数据%首端 末端 电阻 电抗 电纳（1/2)变压器非标准变比 L=4 5 1 4 6 1 5 7 1 6 9 1 7 8 1 8 9 1 1 4 0 0 1 2 7 0 0 1 3 9 0 0 1;发电机数据如下：%发电机 母线 Xd Xd Td0 Xq Xq Tq0 Tj Xf Ge=1 1 0 2 2 3 3 ;系统电路结构拓扑图如下：图 1 WSCC 3 机 9 节点系统（所有参数以 100MVA 为基准值的标幺值）2.潮流计算 首先进行潮流计算，采用牛顿拉夫逊迭代法，电力系统潮流计算是电力系统运行和规划中最基本和最经常的计算，其任务是在已知某些运行参数的情况下，计算出系统中全部 的运行参数，一般来说，各个母线所供负荷的功率是已知的，各个节点电压是未知的（平衡节点除外），可以根据网络结构形成节点导纳矩阵，然后由节点导纳矩阵和网络拓扑结构列写功率方程，由于功率方程里功率是已知的，电压的幅值和相角是未知的，这样潮流计算的问题就转化为求解非线性方程组的问题了。为了便于用迭代法解方程组，需要将上述功率方程改写成功率平衡方程，并对功率平衡方程求偏导，得出对应的雅可比矩阵，给未知节点赋电压初值，一般为额定电压，将初值带入功率平衡方程，得到功率不平衡量，这样由功率不平衡量、雅可比矩阵、节点电压不平衡量（未知的）构成了误差方程，解误差方程，得到节点电压不平衡量，节点电压加上节点电压不平衡量构成新的节点电压初值，将新的初值带入原来的功率平衡方程，并重新形成雅可比矩阵，然后计算新的电压不平衡量，这样不断迭代，不断修正，一般迭代三到五次就能收敛。牛顿拉夫逊算法修正方程 W=-JV 其中 W 是节点不平衡量向量，包括有功，无功，电压；J 是雅克比矩阵；V 是节点电压修正量。令 ijijijiiijBGYjfeV;，则直角坐标形式的功率不平衡量方程为 PQ 节点：0)()(11njnjjijjijjjijjijiisieBfGffBeGePP 0)()(11njnjjijjijjjijjijiisieBfGefBeGfQQ PV 节点：0)()(11njnjjijjijjjijjijiisieBfGffBeGePP)(222222iiisiisifeVVVV 极坐标形式的功率不平衡量方程 0)sincos(1njijijijijjiisiBGVVPP 0)cossin(1njijijijijjiisiBGVVQQ 雅可比矩阵 J 各元素的表达式 LMNHJ 当 ji 时：iijijiijijPHB eG ff iijijiijijPNG eB fe iijijiijijQMB fG ef iijijiijijQLG fB ee 当 j=i 时：iiiii iiiiiiPHB eG fbf iiiii iiiiiiPNG eB fae iiiiiiiiiiiQMG eB faf iiiii iiiiiiQLB eG fbe 其中，nijijjijinijijjijieBfGbfBeGa11)()(。进行牛顿拉夫逊算法迭代后得到电压幅值 V 和相角。3.负荷等效和支路简化 然后求出支路电流，将发电机内电抗 X加入系统导纳矩阵，求出发电机内电势 E。加入发电机内节点后，系统导纳矩阵变成 12*12 阶的矩阵，并将负荷等效成阻抗。然后将支路导纳矩阵分块，如下：AEYCD 其中，A 是 3*3 的方阵，E 是 3*9 的矩阵，C 是 9*3 的矩阵，D 是 9*9 的方阵。经过网络简化得到故障前的 3*3 简化导纳矩阵 Ypre=A-E*(inv(D)*C （1）其中“inv（D）”是 MATLAB 中 D 矩阵的求逆。故障中导纳矩阵的第七行和第七列从矩阵中删除，此时有 AEYCD 此时，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.Simulink-based transient stability analysis of a multimachine power system 中给出的电磁功率计算公式为：nijjjiijijjiiiieiYEEGEP12)cos(本人在此次仿真中用的是公式（2）计算得到的电磁功率。稳态情况下有，机械功率 Pme=Pe 5.求解运动方程 发电机的运动方程可以写成常微分方程组：RinijjjiijijjiiiimijjiRidtdYEEGEPDdtdH12)cos(2 其中Pmi为第i个机组故障前稳态的电磁功率。在本次仿真中Dji为零，即阻尼为零。仿真开始，t=0 时引入故障，后切除故障。求解运动方程后得到 和曲线如下：图 1 曲线 图 2 曲线 2 号发电机 3 号发电机 2 号发电机 1 号发电机 3 号发电机 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 0 0 0 3%2 0 0 0 2%3 0 0 0 2%4 1 0 0 0 0 0 1%5 1 0 0 0 1%6 1 0 0 0 1%7 1 0 0 0 0 0 1%8 1 0 0 0 1 1%9 1 0 0 0 0 0 1;%线路数据%首端 末端 电阻 电抗 电纳（1/2)变压器非标准变比 L=4 5 1 4 6 1 5 7 1 6 9 1 7 8 1 8 9 1 1 4 0 0 1 2 7 0 0 1 3 9 0 0 1;%形成节点导纳矩阵 fb=L(:,1);%起始节点编号 tb=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=9 nl=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）end for 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 end end%节点数据%节点 电压 电压 发电机 发电机 负荷 负荷 节点%号 幅值 相角 有功 无功 有功 无功 类型(1PQ 2PV 3 平衡 N=1 0 0 0 3 2 0 0 0 2 3 0 0 0 2 4 1 0 0 0 0 0 1 5 1 0 0 0 1 6 1 0 0 0 1 7 1 0 0 0 0 0 1 8 1 0 0 0 1 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=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=3 npq=length(pq);%npq=6 Tol=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=Psp-P;%额定有功功率减去流出或流入的有功功率 dQa=Qsp-Q;%额定无功功率减去流出或流入的无功功率 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(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)-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(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);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,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,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);end Vm=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);end Iij=sparse(Iij);%线路数据%首端 末端 电阻 电抗 电纳（1/2)变压器非标准变比 L=4 5 1 4 6 1 5 7 1 6 9 1 7 8 1 8 9 1 1 4 0 0 1 2 7 0 0 1 3 9 0 0 1;Yzengjia=zeros(12,12);Zci=j*j*j*;Yci=1./Zci;for k=1:3 Yzengjia(k,k)=Yci(k);end for k=1:9 for m=1:9 Yzengjia(3+k,3+m)=Y(k,m);end end Yzengjia(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)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);end end E=zeros(3,9);for m=1:3 for n=4:12 E(m,n-3)=Yzengjia(m,n);end end C=zeros(9,3);for m=4:12 for n=1:3 C(m-3,n)=Yzengjia(m,n);end end D=zeros(9,9);for m=4:12 for n=4:12 D(m-3,n-3)=Yzengjia(m,n);end end Ypre=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 Edur(m,n)=Yzengjiadur(m,n+3);%因为 Edur 的第四到第九列是零，而所以这里为了方便，直接取 E 的前六列 end end Cdur=zeros(8,3);for m=1:8 for n=1:3 Cdur(m,n)=Yzengjiadur(m+3,n);%因为 Cdur 的第四到第九行是零，而所以这里为了方便，直接取 C 的前六行 end end Ddur=zeros(8,8);for m=1:8 for n=1:8 Ddur(m,n)=Yzengjiadur(m+3,n+3);end end Ydur=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);end end Eaft=zeros(3,9);for m=1:3 for n=4:12 Eaft(m,n-3)=Yzengjiaaft(m,n);end end Caft=zeros(9,3);for m=4:12 for n=1:3 Caft(m-3,n)=Yzengjiaaft(m,n);end end Daft=zeros(9,9);for m=1:9 for n=1:9 Daft(m,n)=Yzengjiaaft(m+3,n+3);end end Yaft=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)=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=;Pm0=Pe;%机械功率=故障前的电磁功率 X0=ones(3,1).*(2*pi*60),delta;Y_Gen=Ydur;E=abs(Enei);Time_s=0,;%后切除故障 t1,Xout1=ode45(t,X)Gen_fw(t,X,Y_Gen,E,Pm0,Tj),Time_s,X0);Time_s=,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=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_fw(t,X,Y_Gen,E,Pm0,Tj)function Xp=Gen_fw(t,X,Y_Gen,E,Pm0,Tj)Wn=ones(3,1).*2*pi*60;D=zeros(3,1);dW=zeros(3,1);dDeta=zeros(3,1);W=X(1:3);Deta=X(4:6);Gen_U=E.*exp(1i*Deta);Pe=real(Gen_U.*conj(Y_Gen*Gen_U);dW=(Pm0-Pe-D.*W)./Tj.*Wn;dDeta=W-Wn;Xp=dW;dDeta;End 三、课程总结及心得体会 首先感谢任课老师对我仿真的指导与帮助，然后感谢我的导师对我的指导与关怀，还要感谢在程序编写过程中给过我无私帮助的其他同学。虽然我在以前也编写过其他的程序，但这次课程的仿真设计实验让我学到了很多。首 先，对于错误的查找，举个简单的例子，之前我的雅克比矩阵一直有问题，但找了好久没有找到，后来发现是我的节点类型与老师定义的节点类型不一样，将后面的自己的节点类型检索改过后雅克比矩阵就成功了，瞬时感觉豁然开朗。再次，对于某一个参数的计算可能有不同的方法，例如对于电磁功率如果按照文献 Ramnarayan Patel,T.S.Bhatti and D.P.Simulink-based transient stability analysis of a multimachine power system中给出的电磁功率计算公式来计算会变得很繁琐，但经过其他同学的帮助与指点，原来电磁功率就是发电机输出的功率的实部，这样一来计算量就少了很多。最后，完成了这次作业，对于我本人也是一个极大的激励，我之前觉得这个作业是一个很大的挑战，经过这次完成作业我觉得我战胜了一个难题一样，心情无比愉悦。但是，我的仿真还有很多不完善的地方，最明显的就是我的曲线的角度和老师的有很大误差，这是我今后要完善的地方。四、参考文献 1何仰赞.温增银.电力系统分析（上）M.武汉：华中科技大学出版社,2002 2何仰赞.温增银.电力系统分析（下）M.武汉：华中科技大学出版社,2002 3P.M.Anderson.A.A.Fouad.Power System Control and Stability M.(Iowa State University Press,Ames,IA,1977)4P.Kundur.Power System Stability and Control,EPRI Power System Engineering Series(M).(Mc Graw-Hill,New York,1994)(