9节点系统小干扰稳定分析

MACROBUTTONMTEditEquationSection2SEQMTEqn\r\hSEQMTSec\r1\hSEQMTChap\r2\h9节点系统小干扰稳定分析摘要：关键词：引言电力系统在运行过程中无时不遭受到一些小的干扰，例如负荷的随机变化及随后的发电机组调节；因风吹引起架空线路线间距离变化从而导致线路等值电抗的变化，等等。这些现象随时都在发生。互联电力系统可以提高发输电的经济性和可靠性，但是多个地区之间电力系统的多重互联又诱发出许多新的动态问题，使系统失稳的可能性增大。小干扰稳定性和暂态稳定性、电压稳定性、频率稳定性一起作为电力系统保持安全运行的基本要求，其特殊重要的地位逐渐为人们所认识。自从上世纪七十年代以来，世界上多个国家的电力系统都曾发生过因失去小干扰稳定性而导致系统振荡或电压崩溃的严重事故，造成了巨大的经济损失，所以对小干扰稳定性问题的研究逐渐为人们所关注。如今，小干扰稳定性分析已经成为在电力系统的规划和运行中必须进行的重要内容之一。小干扰的发生一般不会引起系统结构的变化。电力系统小干扰稳定分析研究遭受小干扰后电力系统的稳定性。小干扰稳定性是指正常运行的电力系统受到微小的、瞬时出现但又立即消失的扰动后，恢复到它原来运行状况的能力，或者这种扰动虽不消失，但可用原有的运行状态近似地表示可能的新运行状况。系统丧失小干扰稳定性表现为以下几种形式：(1)次同步振荡：电力系统中一种特殊的运行状态，此时系统的输电网络和发电机间以一个或多个略低于系统同步频率的自然频率进行着能量的交换。(2)低频振荡：由于系统中发电机阻尼转矩不足而引起的增幅转子振荡。电力系统中，低频振荡的频率一般在0.2~3Hz范围内。(3)单调不稳定性：在特定的稳态运行条件下，由于缺乏同步转矩而造成的系统不稳定。在实际的大型互联系统中，小干扰稳定性问题通常是由于阻尼不足而引起的低频振荡问题，而低频振荡问题又分为局部和全局两种：(1)局部问题可以分为以下几种：①只与一台发电机或一个单独的电厂相对于系统剩余部分的转子角振荡有关，这种振荡称之为局部电厂模式振荡。②与几台临近的发电机转子之间的振荡有关，这种振荡称之为机间或站间模式振荡。③其他可能的局部问题包括与装置的控制相关的模式的不确定。(2)全局问题：是由大量发电机组之间的相互影响而造成的，因而有广泛的影响，它们表现为一个区域里的一组发电机相对另一区域的一组发电机发生摆动的振荡，这种振荡称之为区域模式振荡。对电力系统小干扰稳定性进行分析的方法大致可以分为以下几种：数值仿真方法、建立在线性模型基础上的分析方法、小干扰稳定域分析方法、非线性理论分析方法和计及模型不确定性的分析方法。对小干扰稳定性问题可以采用线性模型进行研究，这种线性模型是将描述系统动态行为的微分方程和代数方程在稳态运行点处线性化后得到的。目前，建立在线性模型基础上的电力系统小干扰稳定性分析方法主要有两种：以状态空间模型描述为基础的特征值分析法和以传递函数矩阵为基础的频域分析法。本文采用全部特征值分析法（QR法）对9节点系统在稳态处进行小干扰稳定分析。小干扰稳定分析步骤对含有FACTS的交直流系统进行小干扰稳定分析时，其基本流程可总结如下：(1)对给定的西永稳定运行情况进行潮流计算，求出系统各节点电压、电流和功率；(2)将网络方程写成分块矩阵形式，并写出在坐标下节点注入电流偏差与节点电压偏差之间的线性化方程，形成网络导纳矩阵；(3)已知各负荷的功率及负荷节点电压的稳态值为、、、。根据负荷的恒阻抗模型求出负荷模型的导纳矩阵元素、、、，用它们修改导纳矩阵中对应于各负荷节点的对角子块；(4)首先计算出各发电机组中所有变量的初值，然后分别形成矩阵、、、、及、、。然后求出矩阵、、、，从而得到各发电机组的线性化方程。对于其他动态元件，用同样的方法可得到其线性化方程中的系数矩阵。直至得出系统中所有动态元件的线性化方程；(5)形成矩阵、、、，进一步计算出系统的状态矩阵；(6)应用QR法计算矩阵A的全部特征值，从而判断系统在所给定的稳态运行情况下的小干扰稳定性。现对上述步骤进行详细描述。特征值分析法MACROBUTTONMTEditEquationSection2SEQMTEqn\r\hSEQMTSec\r1\hSEQMTChap\r2\h特征值分析法是以李雅普诺夫第一定律为基础，将电力系统模型线性化，以状态空间法将其描述为一般的线性系统，把描述系统动态的一阶微分方程及描述网络的代数方程线性化得到GOTOBUTTONZEqnNum294441REFZEqnNum294441\*Charformat\!(2.1)、GOTOBUTTONZEqnNum125948REFZEqnNum125948\*Charformat\!(2.2)的形式，根据GOTOBUTTONZEqnNum788133REFZEqnNum788133\*Charformat\!(2.3)求解状态矩阵的特征值，然后根据电力系统固有模式和特征值的对应关系，从特征值得到模式的阻尼和频率，从特征向量得到模式和系统中各状态变量的关系，从而给出表示系统稳定性的定性和定量的信息。 上式中，为系统受扰动后的运行状态相对于原平衡状态的偏离量；为运行变量相对于原平衡状态的偏离量。特征值分析法主要分为全部特征值分析法和部分特征值分析法，其优点包括：(1)由于特征值可以准确反映出振荡模式的频率和阻尼以及非振荡模式的衰减率，因此分析人员可以对系统的小扰动情况有一个比较全面的了解，并能够及时发现系统中可能存在的全部弱阻尼模式或不稳定模式。(2)特征值分析法可以提供有关特征值与系统参数间的灵敏度，这将有助于揭示出系统小干扰稳定性的实质，进而寻求改善系统稳定性的对策。(3)特征值分析法不仅可以用于各种电力系统小干扰稳定性问题的分析，而且还可以用于电力系统中各种控制器参数的调整。利用QR法分析电力系统小干扰稳定性的方法已经发展的十分成熟，这种方法数值解稳定，适应性强，并且收敛快，是计算矩阵全部特征值的一种最有效的方法，具有鲁棒性强、收敛速度快等特点。然而，QR法也有一定的局限性，一般只适合于中小规模电力系统的小干扰稳定性分析，这是因为使用QR法用于计算大型电力系统特征值时会遇到以下3个难以克服的困难：(1)QR法不能够利用矩阵的稀疏特性，它需要将矩阵满阵存储，所需的存储容量庞大；(2)使用QR法计算大型矩阵特征值时非常费时；(3)当矩阵的维数过高时还会导致QR法失效而无法得出正确的结果。对于给定的和正交阵，考虑如下迭代： 其中，每个都是正交阵，每个都是上三角阵。由归纳法可得： 这样，每个都与相似。由于有复特征值，因此不会收敛到严格的、“特征值暴露”的三角阵，而只能满足于计算称为实Schur分解的另一种分解。 对角均为块或块的分块上三角阵称为拟上三角阵。实Schur分解相当于将矩阵实归为一个拟上三角阵。若，则存在一个正交阵使得： 其中每个是矩阵或矩阵。若是的，其元素就是的特征值；若是的，的特征值是的一对共轭的复特征值。线性化网络方程为了叙述方便，将网络方程写成分块矩阵形式，并注意到网络方程本身是线性的，因而可以直接写出在坐标系下节点注入电流偏差与节点电压偏差之间的线性化方程： 式中： 对各负荷节点，把负荷的恒阻抗模型给出的注入电流偏差与节点电压偏差关系代入上式，即可消去负荷节点的电流偏差。设负荷接在节点，则消去该负荷后的网络方程仅是对原网络方程GOTOBUTTONZEqnNum267812REFZEqnNum267812\*Charformat\!(2.7)的简单修正：节点的电流偏差变为0，导纳矩阵中的第个对角块变为，而其他内容不变。全系统线性化微分方程各发电机组微分方程的线性化方程可写成如下矩阵形式： 而定子电压方程式的线性化方程表示为： 上式中的和为各发电机本身、轴电压和电流分量的偏差，因此必须把它们转换成统一的同步旋转坐标参考轴下的相应分量，以便将它们和电力网络联系起来。对于发电机电压、电流可推导如下关系： 其中： 将GOTOBUTTONZEqnNum265872REFZEqnNum265872\*Charformat\!(2.11)、GOTOBUTTONZEqnNum800204REFZEqnNum800204\*Charformat\!(2.12)带入GOTOBUTTONZEqnNum123506REFZEqnNum123506\*Charformat\!(2.9)、GOTOBUTTONZEqnNum681572REFZEqnNum681572\*Charformat\!(2.10)消去电流偏差量可得： 其中： 对于各发电机状态变量的个数，取决于发电机选用模型的复杂程度，根据需要确定其励磁系统、调速器等模型。9节点系统描述问题描述9节点电力系统的单线图、支路数据、发电机参数、正常运行情况下的系统潮流分别如REF_Ref471506621\h图31、REF_Ref471506631\h表31、REF_Ref471506633\h表32、REF_Ref471506634\h表33所示。系统频率为60HZ。各负荷均用恒定阻抗模拟。发电机1采用经典模型，发电机2和3采用双轴模型。发电机2和3均装有自并励静止励磁系统，其参数如下：另外，各发电机的阻尼系数均取为1.0.图STYLEREF1\s3SEQ图\*ARABIC\s119节点系统单线结构图表STYLEREF1\s3SEQ表\*ARABIC\s11支路数据首端母线名末端母线名电阻(标幺值)电抗(标幺值)容纳之半(标幺值)变压器非标准变比450.0100.0850.088460.0170.0920.079570.0320.1610.153690.0390.1700.179780.00850.0720.0745890.01190.10080.1045140.00.05761.0270.00.06251.0390.00.05861.0表STYLEREF1\s3SEQ表\*ARABIC\s12发电机数据发电机母线名1147.280.00.14600.06080.06080.06088.961.02212.800.00.89580.11980.86450.19696.000.5351.0336.020.01.31250.18131.25780.25005.890.6001.0表STYLEREF1\s3SEQ表\*ARABIC\s13正常运行情况下的系统潮流母线名电压发电机负荷幅值相角/°有功功率无功功率有功功率无功功率11.0400.00000.71640.270521.02509.28001.63000.066531.02504.66480.8500-0.108641.0258-2.216850.9956-3.98881.25000.500061.0127-3.68740.90000.300071.02583.719781.01590.72751.00000.350091.03241.9667小干扰分析利用潮流结果，计算各发电机变量的初值，如REF_Ref471507017\h表34所示：表STYLEREF1\s3SEQ表\*ARABIC\s14各发电机变量的初值发电机12.271651.039180.041220.678010.287161.05664261.098440.633610.805710.931991.290151.789320.788170.62220354.136620.666070.779090.619410.561471.402900.767860.62424对于发电机1采用经典模型，因而微分方程只有转子运动方程的两阶，励磁系统假定励磁电压保持恒定；对于发电机2和3采用双轴模型，且均装有自并励静止励磁系统，因而其微分方程均为7阶。自并励静止励磁系统传递函数框图如下所示：图STYLEREF1\s3SEQ图\*ARABIC\s12自并励静止励磁系统传递函数框图 发电机1各矩阵计算结果如下所示： 发电机2各矩阵计算结果如下所示： 发电机2各矩阵计算结果如下所示：进一步，系统的线性化方程GOTOBUTTONZEqnNum294441REFZEqnNum294441\*Charformat\!(2.1)中的矩阵为： 计算得、、、如下图所示：图STYLEREF1\s3SEQ图\*ARABIC\s13矩阵A图STYLEREF1\s3SEQ图\*ARABIC\s14矩阵B图STYLEREF1\s3SEQ图\*ARABIC\s15矩阵C图STYLEREF1\s3SEQ图\*ARABIC\s16矩阵D 据此计算状态矩阵如下：图STYLEREF1\s3SEQ图\*ARABIC\s17状态矩阵 运用QR法计算得矩阵如下：图STYLEREF1\s3SEQ图\*ARABIC\s18特征根矩阵 从REF_Ref471513401\h图38中的对角元即可得出9节点系统在稳态下小干扰稳定特征根。小干扰稳定特征值分析上述特征根矩阵中每个对角元为，是矩阵或矩阵。若是的，其元素就是的特征值；若是的，的特征值是的一对共轭的复特征值。因此可求得系统特征根为： 其中可近似为0，即我们已知的零特征值，系统其他所有特征值都具有负实部，因此系统在给定的运行方式下是小干扰稳定的。心得体会附录Small_Signal_Analysis.mclear;clc;Data1=[1400.05760;2700.06250;3900.05860;450.010.0850.088;460.0170.0920.079;570.0320.1610.153;690.0390.170.179;780.00850.0720.0745;890.01190.10080.1045];%支路数据首端母线名末端母线名电阻电抗容纳之半Data2=[11.0400.71640.2705;21.0259.281.630.0665;31.0254.66480.85-0.1086;41.0258-2.216800;50.9956-3.98881.250.5;61.0127-3.68740.90.3;71.02583.719700;81.01590.727510.35;91.03241.966700];%%正常运行情况下系统潮流母线名电压幅值电压相角/°有功无功Data3=[47.2800.1460.06080.06080.06088.9601;12.800.89580.11980.86450.196960.5351;6.0201.31250.18131.25780.255.890.61];%发电机数据TjRaXdXd'XqXq'Td0'Tq0'DData4=zeros(3,8);fort=1:3V=Data2(t,2)*cos(Data2(t,3)*pi/180)+Data2(t,2)*sin(Data2(t,3)*pi/180)*j;I=(Data2(t,4)-Data2(t,5)*j)/conj(V);EQ=V+I*(Data3(t,2)+Data3(t,5)*j);delta=atan(imag(EQ)/real(EQ));Data4(t,1)=delta*180/pi;Data4(t,2)=real(V)*cos(delta)+imag(V)*sin(delta);Data4(t,3)=real(V)*sin(delta)-imag(V)*cos(delta);Data4(t,4)=real(I)*cos(delta)+imag(I)*sin(delta);Data4(t,5)=real(I)*sin(delta)-imag(I)*cos(delta);Data4(t,6)=Data4(t,2)+Data3(t,2)*Data4(t,4)+Data3(t,3)*Data4(t,5);Data4(t,7)=Data4(t,2)+Data3(t,2)*Data4(t,4)+Data3(t,4)*Data4(t,5);Data4(t,8)=Data4(t,3)+Data3(t,2)*Data4(t,5)-Data3(t,6)*Data4(t,4);end%各发电机变量初值δVq0Vd0Iq0Id0Efq0Eq0'Ed0'[M,N]=size(Data1);n=9;nl=[568];Y=zeros(n,n);fort=1:MY(Data1(t,1),Data1(t,2))=-1/(Data1(t,3)+Data1(t,4)*j);endY=Y+Y.';fort=1:nY(t,t)=-sum(Y(t,:));fors=1:Mift==Data1(s,1)||t==Data1(s,2)Y(t,t)=Y(t,t)+Data1(s,5)*j;endend%初始网络导纳矩阵fors=1:length(nl)ift==nl(s)Y(t,t)=Y(t,t)+(Data2(t,4)-Data2(t,5)*j)/(Data2(t,2))^2;break;%恒阻抗负荷模型修正后的网络导纳矩阵endendendGG=real(Y);BB=imag(Y);YY=zeros(2*n,2*n);fort=1:nfors=1:nYY(2*t-1,2*s-1)=GG(t,s);YY(2*t-1,2*s)=-BB(t,s);YY(2*t,2*s-1)=BB(t,s);YY(2*t,2*s)=GG(t,s);endendAg=cell(1,3);BIg=cell(1,3);BVg=cell(1,3);Pg=cell(1,3);Zg=cell(1,3);Zgg=cell(1,3);RVg=cell(1,3);RIg=cell(1,3);Tg=cell(1,3);AG=cell(1,3);BG=cell(1,3);CG=cell(1,3);DG=cell(1,3);wb=120*pi;TB=10;TC=1;TA=0.02;KA=200;TR=0.03;XC=0;fort=1:3ift==1Ag{t}=[0wb;0-Data3(t,9)/Data3(t,1)];BIg{t}=[00;-(Data3(t,5)-Data3(t,4))*Data4(t,4)/Data3(t,1)-(Data4(t,7)+(Data3(t,5)-Data3(t,4))*Data4(t,5))/Data3(t,1)];BVg{t}=[00;00];Pg{t}=[00;00];Zg{t}=[-Data3(t,2)Data3(t,5);-Data3(t,4)-Data3(t,2)];Zgg{t}=pinv(Zg{t});RVg{t}=[Data4(t,2)0;-Data4(t,3)0];RIg{t}=[Data4(t,4)0;-Data4(t,5)0];Tg{t}=[sin(Data4(t,1)*pi/180)-cos(Data4(t,1)*pi/180);cos(Data4(t,1)*pi/180)sin(Data4(t,1)*pi/180)];elseAg{t}=[0wb00000;0-Data3(t,9)/Data3(t,1)-Data4(t,4)/Data3(t,1)-Data4(t,5)/Data3(t,1)000;00-1/Data3(t,7)01/Data3(t,7)00;000-1/Data3(t,8)000;0000-1/TB-(TC/TA-1)/TB-KA*TC/TA/TB;00000-1/TA-KA/TA;000000-1/TR];BIg{t}=[00;-(Data4(t,8)+(Data3(t,6)-Data3(t,4))*Data4(t,4))/Data3(t,1)-(Data4(t,7)+(Data3(t,6)-Data3(t,4))*Data4(t,5))/Data3(t,1);-(Data3(t,3)-Data3(t,4))/Data3(t,7)0;0(Data3(t,5)-Data3(t,6))/Data3(t,8);00;00;00];BVg{t}=[00;00;00;00;00;00;Data4(t,3)/Data2(t,2)/TRData4(t,2)/Data2(t,2)/TR];Pg{t}=[0001000;0010000];Zg{t}=[-Data3(t,2)Data3(t,6);