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

1、发电机最优励磁控制系统设计1 线性最优控制理论1.1 线性最优控制原理线性最优控制理论所研究的核心问题是选取最优的控制规律，使控制系统在特定指标下的性能为最优。控制系统框如图1所示。图1 最优控制框图线性定常系统状态空间方程的一般形式为： （1）式中，为状态系数矩阵；为控制系数矩阵；为n维状态向量；为r维控制向量。如图1所示，如果要改善系统性能，可以引入状态反馈构成闭环系统。反馈系统的状态向量为： （2）式中，为状态反馈增益矩阵。将式（1）和式（2）合并，可以得到： （3）可见，最优控制的本质就是如何选取反馈矩阵，以使它在给定控制规律下达到特定条件下的最优。1.2 二次型性能指标假定为系统的实

2、际响应，为系统预期的响应。最优控制性能指标是使两者的偏差最小，即 （4）式中，是随函数而变化的一个泛函数，是在时间内求取偏差平方的积分，称为二次型指标。如果以为实际状态向量，以为预期状态向量，则要求状态向量偏差最小的二次型性能指标为： （5）在二次型性能指标中，也需要引入对控制量的约束。如果没有这个约束，所设计的控制器中，控制量的变化范围可能会很大，难以实现。引入控制量约束的二次型指标为： （6）式中和表示状态向量和控制向量的权矩阵。为了便于分析，通常把系统平衡点置于状态空间的原点，即，则式（6）可以变换为： （7）以上式作为性能指标设计最优控制系统，可以证明这个最优控制规律是存在且唯一的，表

3、达式为： （8）式中，为最优反馈增益矩阵。为维对称常矩阵，矩阵是黎卡梯方程：的解。2 单机无穷大系统数学模型与无穷大系统并联运行的同步发电机的基本数学模型包括4个部分：与系统并联的同步发电机基本方程组，发电机电势、定子电压和电流之间的联系方程，功率方程以及转子绕组动态方程。2.1与系统并联的同步发电机基本方程组与系统并联的同步发电机基本方程组决定了运行中的发电机的稳态与过渡过程的全部动力学特性。转子运动方程的形式如下： （9）式中，用秒，用弧度，用秒，、均为标幺值。在实用计算中，一般采用，这里的是机组的阻尼系数，用标幺值，单位用弧度。2.2发电机电势、定子电压和电流之间的联系方程将观察发电机的

4、坐标系从abc坐标系转换到dq0坐标系，可以极大地简化发电机的数学模型。这是因为abc坐标系固定在三相对称的定子上，由于转子的旋转，对于每一坐标轴方向上的磁路磁导率都是时间t的周期函数，不便于计算和分析。而dq0坐标系是固定在转子上，坐标轴与磁路相对静止，这样磁导率就是一常数，因而惦记的各个参数也会成为与时间t无关的常数。由发电机电势与电流矢量图可以直接得到如下关系式组： （10）由上述的发电机定子电势、电压和电流的关系式组，可以解出定子电流为： （11） （12） （13）由矢量图又可以得到： （14） （15） （16）代入式（11）和（12），可以解得： （17）2.3 功率方程与无穷大

5、系统并联的发电机有功功率的表达式为： （18）以和为自变量的电磁功率表达式如下： （19）2.4 转子绕组动态方程在列写转子绕组的实用动态方程时，忽略阻尼绕组，只列写励磁绕组动态方程。另外假设定子绕组中磁链式可以突变的，励磁绕组磁链守恒，这是因为定子绕组中暂态电流的非周期分量比励磁绕组非周期分量衰减快得多。带负载的发电机在转子d轴上的总磁链为： （20）式中，是有用磁通；是转子漏磁通；是电枢反应磁通；是励磁绕组与定子d轴绕组的互感；是励磁绕组的漏感。假定过渡过程中转速标幺值为，式（20）可以写为： （21）可以求得： （22）此处，为励磁绕组的自感抗。假设定子d轴总磁链与端电压q轴分量标幺值相

6、等。类似于转子绕组，可以得到： （23）这里为发电机d轴同步电抗，其中是定子漏电抗；为d轴暂态电抗；称为d轴暂态电抗后电势。由励磁绕组等值电路可得回路方程式： （24）以乘以上式，可得在空载情况下励磁绕组的动态方程： （25）上式中，为定子绕组开路时励磁绕组的时间常数。若以乘以式（24），可得： （26）以上就是得到了主发电机励磁绕组的动态方程。2.5 基本方程组的偏差化与线性化由式（9）可以得到偏差方程： （27）根据式（18）和式（19），将电磁功率方程进行偏差化和线性化，可得： （28） （29）式中，；。将式（28）和（29）代入式（27），得到线性化的转子运动方程： （30） （31

7、）由，可得到运行角偏差和机端电压偏差表示的电磁功率偏差： （32）其中，。于是转子的运动偏差方程有： （33）励磁绕组的动态偏差方程为： （34）3 最优励磁控制系统设计 由式（33）和式（34）可以得到系统状态方程为： （35） 由于不便于测量，所以采用机端电压来代替，即 （36）又由于在实际系统中采用的是、和的一组状态量，通过状态变量置换，再结合式（35），可以得到、和所表示的状态方程： （37）如果发电机采用可控硅励磁系统，则可以忽略励磁时间常数，那么控制量就是励磁绕组电压，系统状态方程可以表示为形式的三阶形式： （38）最优控制量4 最优控制设计算例发电机电气参数：隐极机纵轴电抗，纵轴

8、暂态电抗，纵轴开环时间常数，综合阻尼系数；变压器参数：；线路参数：；系统参数：。所研究的发输电系统的示意图见图2。图2 系统示意图设计最优运行点为，时的最优控制器。4.1 计算最优控制参数通过式（18）和式（19），根据最优运行点和，计算得到最优运行点的相关参数如下：由式（18）可得：由式（19）可得：由式（17）可得：数据代入到相应公式，可以求得；选取和作为权矩阵，进而可以得到状态反馈矩阵：可知最优控制量。上述计算过程使用MATLAB完成，见报告附录中的程序1。4.2 验证最优励磁控制系统性能由于励磁控制系统参数是采用线性化状态方程而得到的，而实际系统的状态方程是非线性的，所以使用MATLA

9、B的ode45算法，对单机系统中发电机的运动微分方程进行计算，通过观察发电机的转子角和输出电磁功率的变化情况，得到最优励磁控制系统的控制性能。单机系统的发电机运动方程如下：此运动方程以转子角、角速度和q轴暂态电势作为三个状态量。微分方程中，；作为输入量，引入最优状态反馈矩阵，即用三个状态量进行表示，为其中， 和作为与和有关的变量，分别用式（19）和式（17）表示。4.2.1 静态稳定性系统在运行至5s时刻，加入一个时间长达1秒的电磁功率的扰动，模拟负荷的突然变化，即5s到6s之间从原来的0.55增加到0.6，6s之后回到最初设计的稳态运行点。进行小扰动的仿真。仿真程序见报告附录中的程序27和1

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

11、线路后系统迅速地恢复稳定，整个过程不超过5秒钟。因此，最优励磁控制器也可以很好地维持电力系统的暂态稳定性。5 结论 本报告中先介绍了最优控制理论的相关原理，然后通过对与电网并联的发电机基本方程组的研究从而建立了发电机与励磁系统的线性化状态方程，采用MATLAB计算最优状态反馈参数并代入到系统的非线性微分方程组，然后对系统的非线性微分方程组进行了时域仿真。仿真结果表明，最优励磁控制器可以很好地提高系统的静态稳定性和暂态稳定性。可见，最优励磁控制器有效地增加了发电机励磁系统的阻尼，抑制了系统受到扰动时的振荡，加快系统恢复稳定的速度。附 录程序1：%输入发电机与系统的相关参数Xd=2.5;Xdz=0

12、.5;XT=0.01;Xq=Xd;XL=0.5*1.5;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*Xdzs/Xds;SE=Eq*Vs*cos(ang)/Xds+Vs*Vs*(Xds-Xqs)/(Xds*Xqs)*cos(2*ang);SEz=

13、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*Xe2+Vs2*(cos(ang)2*Xd2+2*Xe*Xd*Eq*Vs*cos(ang)/(Xds2)+(Vs*sin(ang)*Xq/Xqs

14、)2)(-0.5);dVtang=0.5*(Vs2*Xq2*sin(2*ang)/Xqs2-(Vs2*Xd2*sin(2*ang)+2*Xe*Xd*Eq*Vs*sin(ang)/Xds2)*(Eq2*Xe2+(Vs*cos(ang)*Xd)2+2*Xd*Xe*Eq*Vs*cos(ang)/Xds2+(Vs*sin(ang)*Xq/Xqs)2)(-0.5);RV=RE/dVtEq;SV=SE-RV*dVtang;w0=2*pi*f0;%计算得到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/

15、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(1)/Xdzs+0.5*Vs2*(Xdzs-Xds)/(Xdzs*Xds)*sin(2*y(1);Eq=Xds*(y(3)/Xdzs+Vs*(Xdzs-Xds

16、)/(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.7617*(Vt-1.0487)+Eq0)/Td0-y(3)/Td+1/Td0*(Xd-Xdz)/Xdzs*Vs*cos(y(1);dy=dy(1);dy(2);d

17、y(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)/(Xdzs*Xds)*sin(2*y(i,1);End程序4：%某线路首端三相短路故障中的系统状态function dy=gz2(t,y);Xd=2.5;Xdz=0.5XT=

18、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.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)-6

19、8.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 dy=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*

20、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)-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

21、)+Eq0)/Td0-y(3)/Td+1/Td0*(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: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);e

22、nd程序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*(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*X

23、dzs/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(2);dy(3);程序9：%计算扰动中的电磁功率function Pe2=dcgl2(y)Pe2=zeros(length(y(:,1),1);Xd=2.5;Xdz=0.5;XT=0.01;X

24、L=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=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)

25、/(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.7617*(Vt-1.0487)+Eq0)/Td0-y(3

26、)/Td+1/Td0*(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：%设定仿真

27、时间为40秒,5秒时电磁功率受到扰动,6秒时刻扰动结束t0=0;t1=5;t2=6;t3=40;%扰动前，以选定的最优运行点为初值T1,Y1=ode45(gz1,t0 t1,65*pi/180 100*pi 1.0487);Pe1=dcgl1(Y1);ac=Y1(length(Y1),1);wc=Y1(length(Y1),2);ec=Y1(length(Y1),3);%扰动中，以扰动前状态量的最后一组解为扰动中的初值T2,Y2=ode45(rd,t1 t2,ac wc ec);Pe2=dcgl2(Y2);ac2=Y2(length(Y2),1);wc2=Y2(length(Y2),2);ec2=Y2(length(Y2)