发电机最优励磁控制系统设计讲解

1、发电机最优励磁控制系统设计1 线性最优控制理论1.1 线性最优控制原理 线性最优控制理论所研究的核心问题是选取最优的控制规律， 使控制系统在特定指标下的性能为最优。控制系统框如图 1 所示。图1 最优控制框图 线性定常系统状态空间方程的一般形式为：X(t) AX BU(1)式中， A为状态系数矩阵； B为控制系数矩阵； X为n维状态向量； U 为r 维控制向量。如图 1 所示，如果要改善系统性能， 可以引入状态反馈构成闭环系统。 反馈 系统的状态向量为：U KX( 2)式中， K 为状态反馈增益矩阵。将式( 1)和式( 2)合并，可以得到：X AX B( KX) (A BK)X(3)可见，最优

2、控制的本质就是如何选取反馈矩阵 K ，以使它在给定控制规律下 达到特定条件下的最优。1.2 二次型性能指标假定 y(t) 为系统的实际响应， (t) 为系统预期的响应。最优控制性能指标是 使两者的偏差最小，即J 0 (t) y(t) 2dt Jmin( 4)式中， J 是随函数 y(t) 而变化的一个泛函数，是在 0 时间内求取偏差平方的积分，称为二次型指标如果以 X (t )为实际状态向量，以 X?(t )为预期状态向量，则要求状态向量偏 差最小的二次型性能指标为：J 0 X?(t) X(t) X?(t) X(t) dt(5)在二次型性能指标中， 也需要引入对控制量的约束。 如果没有这个约束

3、， 所 设计的控制器中， 控制量的变化范围可能会很大， 难以实现。 引入控制量约束的 二次型指标为：TTJ 0 X?(t) X(t) Q X?(t) X(t) UT(t)RU(t) dt(6)式中Q 和 R表示状态向量和控制向量的权矩阵为了便于分析，通常把系统平衡点置于状态空间的原点，即X?(t) 0 ，则式(6)可以变换为：J 0 XT (t)QX(t) UT (t)RU(t) dt(7)以上式作为性能指标设计最优控制系统， 可以证明这个最优控制规律是存在 且唯一的，表达式为：U(t) R 1BT PX (t ) KX(t)(8)式中，K 为最优反馈增益矩阵。 P为n n维对称常矩阵，矩阵

4、P是黎卡梯方 程： ATP PA PBR 1BTP Q 0的解。2 单机无穷大系统数学模型与无穷大系统并联运行的同步发电机的基本数学模型包括 4 个部分：与系统 并联的同步发电机基本方程组，发电机电势、定子电压和电流之间的联系方程， 功率方程以及转子绕组动态方程。2.1 与系统并联的同步发电机基本方程组与系统并联的同步发电机基本方程组决定了运行中的发电机的稳态与过渡 过程的全部动力学特性。转子运动方程的形式如下：H d 22 f dt2Pm Pe PD9)式中，H 用秒， 用弧度， t用秒， Pm 、 Pe 、PD均为标幺值。 在实用计算中，一般采用 PDD d ，这里的 D 是机组的阻尼系数

5、，用标幺值，单位用弧度。D 2 f0 dt2.2 发电机电势、定子电压和电流之间的联系方程将观察发电机的坐标系从 abc 坐标系转换到 dq0 坐标系，可以极大地简化发 电机的数学模型。 这是因为 abc 坐标系固定在三相对称的定子上， 由于转子的旋 转，对于每一坐标轴方向上的磁路磁导率都是时间 t 的周期函数，不便于计算和 分析。而 dq0 坐标系是固定在转子上， 坐标轴与磁路相对静止， 这样磁导率就是 一常数，因而惦记的各个参数也会成为与时间 t 无关的常数。由发电机电势与电流矢量图可以直接得到如下关系式组：Eq Vsq I d xdEq Eq Id（ xd xd）Eq Vtq Id xd

6、Vsq Eq I d xdVsq Eq Id xd（10）Vsd I qxqVsq Vs cosVsd Vs sinVs Vs2q Vs2d由上述的发电机定子电势、电压和电流的关系式组，可以解出定子电流为：Eq Vs cos Id xd（11）IVs sinIqxq（12）Eq Vs cosqsIdxd（13）由矢量图又可以得到：Vtq Vs cos I dxs（14）Vtq Vs sin Iqxs（15）VtVtq2 Vtd2（16）2.3代入式（ 11）和（ 12），可以解得：Vt2 2 2 2 2 2 2Eqxs Vs cos xd 2xsxdEqVs cosVs sinxq22xd2x

7、q217)功率方程与无穷大系统并联的发电机有功功率的表达式为：Pe IqVsq IdVsd ExqdVssinV2s xdxq sin218)xd xq以 Eq 和 为自变量的电磁功率 Pe 表达式如下：Pe EqVs sin Vs2 xd xq sin2（19）xd2 xd xq2.4 转子绕组动态方程 在列写转子绕组的实用动态方程时，忽略阻尼绕组，只列写励磁绕组动态方 程。另外假设定子绕组中磁链式可以突变的， 励磁绕组磁链守恒， 这是因为定子 绕组中暂态电流的非周期分量比励磁绕组非周期分量衰减快得多。带负载的发电机在转子 d 轴上的总磁链为：fd d0 fl ad if Lad if Lf

8、l I dLad（20）式中， d0是有用磁通； fl是转子漏磁通； ad 是电枢反应磁通； Lad是励磁绕21)22)组与定子 d轴绕组的互感； Lfl 是励磁绕组的漏感。 假定过渡过程中转速标幺值为 * 1，式（ 20）可以写为： fd if xad if xfl Idxad 可以求得：fd I dxadifxf此处， xf xad xfl 为励磁绕组的自感抗。假设定子 d轴总磁链与端电压 q轴分量标幺值相等。 类似于转子绕组，可以 得到：223)Vtqfd xadId(xdxad )xf xf这里 xd xad xl 为发电机 d 轴同步电抗，其中 xl 是定子漏电抗；xd xd2xad

9、xfd 轴暂态电抗； Eq fd xad 称为 d 轴暂态电抗后电势。xf由励磁绕组等值电路可得回路方程式：difVf i f Rf xfdt24)以xRafd乘以上式，可得在空载情况下励磁绕组的动态方程：Ef Eq Td0 ddEtqdt25)上式中， Tdxf 为定子绕组开路时励磁绕组的时间常数Rf若以 xad 乘以式（ 24），可得：Rf26)Ef Eq Td0 ddEtqdt以上就是得到了主发电机励磁绕组的动态方程。2.5 基本方程组的偏差化与线性化由式（ 9）可以得到偏差方程：HPm Pe PD Pm Pe2fD2 f027)28)29)根据式（ 18）和式（ 19），将电磁功率方程

10、进行偏差化和线性化，可得：Pe SERE EqPe SERE Eq式中， SE EqVs cos Vs2 xd xq cos2 ； SE EqVs cos Vs2 xd xq cos2 ；xd xd xq xd xd xqREs sin ； REsxdsin 。 xd将式（ 28）和（ 29）代入式（ 27），得到线性化的转子运动方程：HPm SERE Eq2 f0D2 f030)由 VtVtHPmSEREEq D2 f0 m EE q2 f031)EVtq Eq ，可得到运行角偏差和机端电压偏差Vt 表示32)RE。其中， SV SE RV，RVVt Eq的电磁功率偏差：Pe SVRV Vt

11、于是转子的运动偏差方程有：HD33)Pm SVRV Vt2 f02 f0励磁绕组的动态偏差方程为：34)EfEq Td0 Eq3 最优励磁控制系统设计由式（ 33）和式（ 34）可以得到系统状态方程为：010000 SED00 RE0HEHHEEqSESE011EqqTd0RETdTd0qEf1Ef000000 VR135)由于 Eq 不便于测量，所以采用机端电压Vt来代替 Eq，EqSE SERV VtqREREt36)又由于在实际系统中采用的是 Pe、 、 Vt和 Ef 的一组状态量，通过状态变量置换，再结合式（35），可以得到 Pe、 、 Vt和 Ef 所表示的状态方程：SESVPeVt

12、EfTdSV0HSE SVTd0RV SVSERVSETdSVRETd0DHSE SVRVSETdSVREPeVt如果发电机采用可控硅励磁系统，Td0RV1TeEfTe000 VR137)则可以忽略励磁时间常数 Te 0 ，那么控形式的制量U 就是励磁绕组电压 Ef ，系统状态方程可以表示为 X AX BU三阶形式：SE SVPeVtTdSV0HSE SVSEDHSE SVTd0RV SVRVRVSETdSVSETdSVREPeVtTd00RETd0RVEf38)最优控制量 U kX4 最优控制设计算例发电机电气参数： 隐极机纵轴电抗 Xd2.5，纵轴暂态电抗 Xd 0.5 ，纵轴开环时间常数

13、 Td0 10s，H 8s，综合阻尼系数 D 5.0 ；变压器参数： XT 0.01；1线路参数： XL 1 1.5；系统参数： Vs 1。所研究的发输电系统的示意图见图22。Vs发电机XT 0.01XL 1 1.5L2图2 系统示意图 设计最优运行点为 0 65 ， P0 0.55时的最优控制器。由式18）可得：Eq0Pexd0.55 (2.5 0.01 0.75) 1.9784Vs sin1 sin65由式19）可得：Eq02(Pe Vs xd xq sin2 )Vsxsdin1.02392 xd xq由式17）可得：Eq2xs Vs2 cos2 xd2 2xs xd EqVs cosVs

14、2 siVt022sin xq2 q 1.0487xq2数据代入到相应公式，可以求得-0.16950.6564-0.14220.0719-39.2699-0.62500；B0-0.1064-0.0545-0.08930.0452A2 xd4.1 计算最优控制参数 通过式（ 18）和式（ 19），根据最优运行点 0 65 和 P0 0.55 ，计算得到 最优运行点的相关参数如下：选取 Q diag （1,100,5000） 和 R 1作为权矩阵，进而可以得到状态反馈矩阵：K lqr(A,B,Q,R) 46.5370 -6.1560 68.7617可知最优控制量 u EfKX 58.7480 Pe

15、 5.176730.0717 Vt 。上述计 算过程使用 MATLAB 完成，见报告附录中的程序 1。4.2 验证最优励磁控制系统性能 由于励磁控制系统参数是采用线性化状态方程而得到的，而实际系统的状 态方程是非线性的， 所以使用 MATLAB 的 ode45 算法，对单机系统中发电机的 运动微分方程进行计算， 通过观察发电机的转子角和输出电磁功率的变化情况，得到最优励磁控制系统的控制性能。单机系统的发电机运动方程如下：0H0PmD (0)0 ( EqVs sinVs xdxq sin2 )H H xd2 xd xqEq1 EfTd01TdEq1 xd xd Vs cosTd0 xd此运动方程

16、以转子角 、角速度 和 q 轴暂态电势 Eq作为三个状态量。微分方程中， 0 2 f 100 ； Pm 0.55；Ef 作为输入量，引入最优状态反馈矩 阵，即用三个状态量进行表示，为EfEq0EfEq0( KX)Eq0( 58.7480Pe5.1767 30.0717Vt)Eq0 ( 46.5370 (Pe P0) 6.1560 ( 0) 68.7617 (Vt Vt0)其中， Pe 和Vt作为与 Eq和 有关的变量，分别用式（ 19）和式（ 17）表示。4.2.1 静态稳定性系统在运行至 5s时刻，加入一个时间长达 1 秒的电磁功率的扰动，模拟负 荷的突然变化，即 5s到6s之间Pe从原来的

17、 0.55增加到 0.6，6s之后Pe回到最27初设计的稳态运行点。进行小扰动的仿真。仿真程序见报告附录中的程序和 13。仿真的功角和电磁功率波形图如下：图3 静态稳定功角仿真图从波形图可以看出，在发生小扰动后，采用最优控制理论所设计的励磁系统，能够快速地平息系统的振荡，转子角和电磁功率均能够迅速地回到最初设计的最优运行点，提高了电力系统的静态稳定性。可见，发电机加入最优励磁 控制器可以有效地改善系统稳定性。4.2.2 暂态稳定性系统运行至 5s的时刻，在其中一条线路首端即图 2中的 d处，发生三相短 路故障，故障 0.1s 后保护动作切除故障线路。程序见报告附录中的程序2、3和 812。该过

18、程中的功角变化波形图如下：图4 暂态稳定功角仿真图由图 4 可以看出，系统在 5s时刻发生三相短路后功角和电磁功率急剧上升， 5.1s切除故障线路后系统迅速地恢复稳定，整个过程不超过5 秒钟。因此，最优励磁控制器也可以很好地维持电力系统的暂态稳定性。5 结论本报告中先介绍了最优控制理论的相关原理，然后通过对与电网并联的发 电机基本方程组的研究从而建立了发电机与励磁系统的线性化状态方程，采用 MATLAB 计算最优状态反馈参数并代入到系统的非线性微分方程组， 然后对系 统的非线性微分方程组进行了时域仿真。仿真结果表明，最优励磁控制器可以 很好地提高系统的静态稳定性和暂态稳定性。可见，最优励磁控制

19、器有效地增 加了发电机励磁系统的阻尼，抑制了系统受到扰动时的振荡，加快系统恢复稳 定的速度。附录程序 1：%输入发电机与系统的相关参数Xd=2.5;Xdz=0.5;XT=0.01;Xq=Xd;Td0=10;Vs=1;D=5;H=8;Pe=0.55;ang=65*pi/180;f0=50;%中间计算过程Xe=XT+XL;Xds=Xd+XT+XL;Xdzs=Xdz+XT+XL;Xqs=Xq+XT+XL;Eq=Pe*Xds/(Vs*sin(ang) Eqz=(Pe-0.5*Vs2*(Xdzs-Xqs)/(Xdzs*Xqs)*sin(2*ang)*Xdzs/(Vs*sin(ang)Td=Td0*Xdz

20、s/Xds; SE=Eq*Vs*cos(ang)/Xds+Vs*Vs*(Xds-Xqs)/(Xds*Xqs)*cos(2*ang);SEz=Eqz*Vs*cos(ang)/Xdzs+Vs*Vs*(Xdzs-Xqs)/(Xdzs*Xqs)*cos(2*ang);RE=Vs*sin(ang)/Xds;REz=Vs*sin(ang)/Xdzs; Vt=sqrt(Eq2*Xe2+(Vs*cos(ang)*Xd)2+2*Xe*Xd*Eq*Vs*cos(ang)/Xds2+(Vs*sin(ang)* Xq/Xqs)2) dVtEq=(Eq*Xe2+Xe*Xd*Vs*cos(ang)/(Xds2)*(Eq2

21、*Xe2+Vs2*(cos(ang)2*Xd2+2 *Xe*Xd*Eq*Vs*cos(ang)/(Xds2)+(Vs*sin(ang)*Xq/Xqs)2)(-0.5);dVtang=0.5*(Vs2*Xq2*sin(2*ang)/Xqs2-(Vs2*Xd2*sin(2*ang)+2*Xe*Xd*Eq*Vs*sin(an g)/Xds2)*(Eq2*Xe2+(Vs*cos(ang)*Xd)2+2*Xd*Xe*Eq*Vs*cos(ang)/Xds2+(Vs*sin(an g)*Xq/Xqs)2)(-0.5);RV=RE/dVtEq;SV=SE-RV*dVtang;w0=2*pi*f0;%计算得到

22、A,B,K 矩阵Q=1 0 0;0 100 0;0 0 5000;R=1;A=(SE-SV)/(Td*SV) SEz -RV*SE/(Td*SV);-w0/H -D/H 0;(SE-SV)/(Td*RV*SV) (SEz-SV)/RV -SE/(Td*SV)B=REz/Td0;0;REz/(Td0*RV)K=lqr(A,B,Q,R)程序 2：%故障前或扰动前的系统状态function dy=gz1(t,y);Xd=2.5;Xdz=0.5XT=0.01;XL=0.75;Xe=XT+XL;Xds=Xd+Xe;Xdzs=Xdz+Xe;Vs=1;Eq0=1.9784;Pe=y(3)*Vs*sin(y(

23、1)/Xdzs+0.5*Vs2*(Xdzs-Xds)/(Xdzs*Xds)*sin(2*y(1);Eq=Xds*(y(3)/Xdzs+Vs*(Xdzs-Xds)/(Xdzs*Xds)*cos(y(1);Vt=1/Xds*sqrt(Eq2*Xe2+Vs2*Xd2+2*Xe*Xd*Eq*Vs*cos(y(1);Td0=10;Td=Td0*Xdzs/Xds;dy(1)=y(2)-100*pi;dy(2)=2*pi*50*0.55/8-5/8*(y(2)-2*pi*50)-2*pi*50/8*Pe;dy(3)=(-46.537*(Pe-0.55)+6.156*(y(2)-2*pi*50)-68.761

24、7*(Vt-1.0487)+Eq0)/Td0-y(3)/Td+1/Td 0*(Xd-Xdz)/Xdzs*Vs*cos(y(1);dy=dy(1);dy(2);dy(3);程序 3： %计算故障前或扰动前的电磁功率 function Pe1=dcgl1(y) Pe1=zeros(length(y(:,1),1);Xd=2.5;Xdz=0.5;XT=0.01;XL=0.75;Xe=XT+XL;Xds=Xd+Xe;Xdzs=Xdz+Xe;Vs=1;for i=1:length(y(:,1);Pe1(i)=y(i,3)*Vs*sin(y(i,1)/Xdzs+0.5*Vs2*(Xdzs-Xds)/(Xd

25、zs*Xds)*sin(2*y(i,1);End程序 4：%某线路首端三相短路故障中的系统状态function dy=gz2(t,y);Xd=2.5;Xdz=0.5XT=0.01;Xe=XT;Xds=Xd+Xe;Xdzs=Xdz+Xe;Vs=1;Pe=0;Eq=Xds*(y(3)/Xdzs+Vs*(Xdzs-Xds)/(Xdzs*Xds)*cos(y(1);Vt=1/Xds*sqrt(Eq2*Xe2+Vs2*Xd2+2*Xe*Xd*Eq*Vs*cos(y(1);Eq0=1.9784;Td0=10;Td=Td0*Xdzs/Xds;dy(1)=y(2)-100*pi;dy(2)=2*pi*50*0

26、.55/8-5/8*(y(2)-2*pi*50)-2*pi*50/8*Pe;dy(3)=(-46.537*(Pe-0.55)+6.156*(y(2)-2*pi*50)-68.7617*(Vt-1.0487)+Eq0)/Td0-y(3)/Td+1/Td0*(Xd-Xdz)/Xdzs*Vs*cos(y(1);dy=dy(1);dy(2);dy(3);程序 5：%计算故障中的电磁功率function Pe4=dcgl4(y)Pe4=zeros(length(y(:,1),1);for i=1:length(y(:,1);Pe4(i)=0;End程序 6：%切除故障线路后的系统状态function d

27、y=gz3(t,y);Xd=2.5;Xdz=0.5XT=0.01;XL=1.5;Xe=XT+XL;Xds=Xd+Xe;Xdzs=Xdz+Xe;Vs=1;Pe=y(3)*Vs*sin(y(1)/Xdzs+0.5*Vs2*(Xdzs-Xds)/(Xdzs*Xds)*sin(2*y(1);Eq=Pe*Xds/(Vs*sin(y(1);Vt=1/Xds*sqrt(Eq2*Xe2+Vs2*Xd2+2*Xe*Xd*Eq*Vs*cos(y(1);Eq0=1.9784;Td0=10;Td=Td0*Xdzs/Xds;dy(1)=y(2)-100*pi;dy(2)=2*pi*50*0.55/8-5/8*(y(2)

28、-2*pi*50)-2*pi*50/8*Pe;dy(3)=(-46.537*(Pe-0.55)+6.156*(y(2)-2*pi*50)-68.7617*(Vt-1.0487)+Eq0)/Td0-y(3)/Td+1/Td 0*(Xd-Xdz)/Xdzs*Vs*cos(y(1);dy=dy(1);dy(2);dy(3);程序 7： %计算故障后的电磁功率 function Pe5=dcgl5(y)Pe5=zeros(length(y(:,1),1);Xd=2.5;Xdz=0.5;XT=0.01;XL=1.5;Xe=XT+XL;Xds=Xd+Xe;Xdzs=Xdz+Xe;Vs=1;for i=1:

29、length(y(:,1);Pe5(i)=y(i,3)*Vs*sin(y(i,1)/Xdzs+0.5*Vs2*(Xdzs-Xds)/(Xdzs*Xds)*sin(2*y(i,1); end程序 8： %发生扰动的系统状态 function dy=rd(t,y);Xd=2.5;Xdz=0.5;XT=0.01;XL=0.75;Xe=XT+XL;Xds=Xd+Xe;Xdzs=Xdz+Xe;Vs=1;Eq0=1.9784;Pe=y(3)*Vs*sin(y(1)/Xdzs+0.5*Vs2*(Xdzs-Xds)/(Xdzs*Xds)*sin(2*y(1);Eq=Xds*(y(3)/Xdzs+Vs*(Xdz

30、s-Xds)/(Xdzs*Xds)*cos(y(1);Vt=1/Xds*sqrt(Eq2*Xe2+Vs2*Xd2+2*Xe*Xd*Eq*Vs*cos(y(1);Td0=10;Td=Td0*Xdzs/Xds; dy(1)=y(2)-100*pi;dy(2)=2*pi*50*0.6/8-5/8*(y(2)-2*pi*50)-2*pi*50/8*Pe;dy(3)=(-46.537*(Pe-0.6)+6.156*(y(2)-2*pi*50)-68.7617*(Vt-1.0487)+Eq0)/Td0-y(3)/Td+1/Td0 *(Xd-Xdz)/Xdzs*Vs*cos(y(1);dy=dy(1);dy

31、(2);dy(3);程序 9： %计算扰动中的电磁功率 function Pe2=dcgl2(y)Pe2=zeros(length(y(:,1),1);Xd=2.5;Xdz=0.5;XT=0.01;XL=0.75;Xe=XT+XL;Xds=Xd+Xe;Xdzs=Xdz+Xe;Vs=1;for i=1:length(y(:,1);Pe2(i)=y(i,3)*Vs*sin(y(i,1)/Xdzs+0.5*Vs2*(Xdzs-Xds)/(Xdzs*Xds)*sin(2*y(i,1); end程序 10：%扰动后的系统状态 function dy=rdh(t,y);Xd=2.5;Xdz=0.5;XT=

32、0.01;XL=0.75;Xe=XT+XL;Xds=Xd+Xe;Xdzs=Xdz+Xe;Vs=1;Eq0=1.9784;Pe=y(3)*Vs*sin(y(1)/Xdzs+0.5*Vs2*(Xdzs-Xds)/(Xdzs*Xds)*sin(2*y(1);Eq=Xds*(y(3)/Xdzs+Vs*(Xdzs-Xds)/(Xdzs*Xds)*cos(y(1);Vt=1/Xds*sqrt(Eq2*Xe2+Vs2*Xd2+2*Xe*Xd*Eq*Vs*cos(y(1);Td0=10;Td=Td0*Xdzs/Xds;dy(1)=y(2)-100*pi;dy(2)=2*pi*50*0.55/8-5/8*(y(

33、2)-2*pi*50)-2*pi*50/8*Pe;dy(3)=(-46.537*(Pe-0.55)+6.156*(y(2)-2*pi*50)-68.7617*(Vt-1.0487)+Eq0)/Td0-y(3)/Td+1/Td 0*(Xd-Xdz)/Xdzs*Vs*cos(y(1);dy=dy(1);dy(2);dy(3);程序 11：%计算扰动后的电磁功率function Pe3=dcgl3(y)Pe3=zeros(length(y(:,1),1);Xd=2.5;Xdz=0.5;XT=0.01;XL=0.75;Xe=XT+XL;Xds=Xd+Xe;Xdzs=Xdz+Xe;Vs=1;for i=1:length(y(:,1);Pe3(i)=y(i,3)*Vs*sin(y(i,1)/Xdzs+0.5*Vs2*(Xdzs-Xds)/(Xdzs*Xds)*sin(2*y(i,1); end程序 12：%设定仿真时间为 40 秒,5 秒时电磁功率受到扰动 ,6 秒时刻扰动结束 t0=0;t1=5;t2=6;t3=40;%扰动前，以选定的最优运行点为初值 T1,Y1=ode45(gz1,t0 t1,65*pi/180 100*pi 1.0487)