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

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

2、X =+-=- (3 可见, 最优控制的本质就是如何选取反馈矩阵 K , 以使它在给定控制规律下 达到特定条件下的最优。 1.2 二次型性能指标假定 ( y t 为系统的实际响应, ( t 为系统预期的响应。最优控制性能指标是 使两者的偏差最小,即2min 0( ( J t y t dt J =-= (4式中, J 是随函数 ( y t 而变化的一个泛函数,是在 0 时间内求取偏差平方的积分,称为二次型指标。如果以 ( X t 为实际状态向量,以 ( Xt 为预期状态向量,则要求状态向量偏 差最小的二次型性能指标为:( ( ( ( TJ X t X t X t X t dt =- (5 在二次

3、型性能指标中, 也需要引入对控制量的约束。 如果没有这个约束, 所 设计的控制器中, 控制量的变化范围可能会很大, 难以实现。 引入控制量约束的 二次型指标为:( ( ( ( ( ( TT J Xt X t Q X t X t U t RU t dt =-+(6 式中 Q 和 R 表示状态向量和控制向量的权矩阵。为了便于分析,通常把系统平衡点置于状态空间的原点,即 ( 0Xt =,则式 (6可以变换为:0( ( ( ( T TJ X t QX t U t RU t dt =+ (7 以上式作为性能指标设计最优控制系统, 可以证明这个最优控制规律是存在 且唯一的,表达式为:1( ( ( T U

4、t R B PX t KX t -=-=- (8式中, K 为最优反馈增益矩阵。 P 为 n n 维对称常矩阵,矩阵 P 是黎卡梯方 程:10T T A P PA PBR B P Q -+-+=的解。2 单机无穷大系统数学模型与系统并联的同步发电机基本方程组决定了运行中的发电机的稳态与过渡 过程的全部动力学特性。转子运动方程的形式如下:222m e D H d P P P f dt =- (9式中, H 用秒, 用弧度, t 用秒, m P 、 e P 、 D P 均为标幺值。 在实用计算中, 一般采用 02D D d P f dt=,这里的 D 是机组的阻尼系数,用标幺值,单位用弧度。2.2

5、发电机电势、定子电压和电流之间的联系方程将观察发电机的坐标系从 abc 坐标系转换到 dq0坐标系,可以极大地简化发 电机的数学模型。 这是因为 abc 坐标系固定在三相对称的定子上, 由于转子的旋 转,对于每一坐标轴方向上的磁路磁导率都是时间 t 的周期函数,不便于计算和 分析。 而 dq0坐标系是固定在转子上, 坐标轴与磁路相对静止, 这样磁导率就是 一常数,因而惦记的各个参数也会成为与时间 t 无关的常数。由发电机电势与电流矢量图可以直接得到如下关系式组:(cos sin q sq d dqq d d d qtq d d sq q d d sq q d d sd q q sq s sd

6、s s E V I x E E I x x E V I x V E I x V E I x V I x V V V V V =+=+-=+=-=-= (10 由上述的发电机定子电势、电压和电流的关系式组,可以解出定子电流为:cos q s d d E V I x -=(11sin s q q V I x =(12 cos qs d d E V I x -= (13 由矢量图又可以得到:cos tq s d s V V I x =+ (14 sin tq s q s V V I x =+ (15t V = (16 代入式(11和(12 ,可以解得:t V =(17 2.3 功率方程与无穷大系统并

7、联的发电机有功功率的表达式为:2sin sin 22q s d q s e q sq d sd d d q E V x x V P I V I V x x x -=+=+ (18以 qE 和 为自变量的电磁功率 e P 表达式如下: 2sin sin 22q s dq s e d dq E V x x V P x x x -=+ (192.4 转子绕组动态方程在列写转子绕组的实用动态方程时,忽略阻尼绕组,只列写励磁绕组动态方 程。 另外假设定子绕组中磁链式可以突变的, 励磁绕组磁链守恒, 这是因为定子 绕组中暂态电流的非周期分量比励磁绕组非周期分量衰减快得多。带负载的发电机在转子 d 轴上的总

8、磁链为:0fd d fl ad f ad f fl d ad i L i L I L =+-=+- (20式中, 0d 是有用磁通; fl 是转子漏磁通; ad 是电枢反应磁通; ad L 是励磁绕 组与定子 d 轴绕组的互感; fl L 是励磁绕组的漏感。 假定过渡过程中转速标幺值为 *1=,式(20可以写为:fd f ad f fl d ad i x i x I x =+- (21可以求得:fd d adf fI x i x +=(22此处, f ad fl x x x =+为励磁绕组的自感抗。假设定子 d 轴总磁链与端电压 q 轴分量标幺值相等。 类似于转子绕组, 可以 得到:2( ad

9、 adtq fdd d f fx x V I x x x =- (23 这里 d ad l x x x =+为发电机 d 轴同步电抗,其中 l x 是定子漏电抗; 2add d fx x x x =-为d 轴暂态电抗; adqfd fx E x =称为 d 轴暂态电抗后电势。 由励磁绕组等值电路可得回路方程式:f f f f fdi V i R x dt=+ (24以adfx R 乘以上式,可得在空载情况下励磁绕组的动态方程: 0q f q d dE E E T dt=+ (25上式中, 0f d fx T R =为定子绕组开路时励磁绕组的时间常数。若以adfx R 乘以式(24 ,可得: 0

10、q f q d dE E E T dt=+ (26以上就是得到了主发电机励磁绕组的动态方程。 2.5 基本方程组的偏差化与线性化由式(9可以得到偏差方程:22m e D m e H D P P P P P f f =-=- (27 根据式(18和式(19 ,将电磁功率方程进行偏差化和线性化,可得:e E E q P S R E =+ (28e E E qP S R E =+ (29 式中, 2cos cos 2q s d q E sd d q E V x x S Vx x x -=-; 2cos cos 2q s d q E s d dq E V x x S V x x x -=-;sin s

11、 E d V R x =; sin s E dVR x =。将式(28和(29代入式(27 ,得到线性化的转子运动方程:0022m E E q H DP S R E f f =- (300022m E E q H DP S R E f f =- (31由 t t t q qV VV E E =+,可得到运行角偏差 和机端电压偏差 t V 表示 的电磁功率偏差:e V V t P S R V =+ (32其中, tV E VV S S R =-, E V t q R R V E =。于是转子的运动偏差方程有:0022m V V t H D P S R V f f =- (33励磁绕组的动态偏差方

12、程为:0f q d qE E T E =+ (34 3 最优励磁控制系统设计由式(33和式(34可以得到系统状态方程为:0000010000011001100E E E E R q q d E d d f f e e DS R HH H S SV E E T R T T E E T T -=+-(35 由于 qE 不便于测量,所以采用机端电压 t V 来代替 q E ,即 E E V qt E E S S RE V R R -=+ (36又由于在实际系统中采用的是 e P 、 、 t V 和 f E 的一组状态量,通过 状态变量置换,再结合式(35 ,可以得到 e P 、 、 t V 和 f

13、E 所表示的状态 方程:00E VV E E E d V d V d e eRE V E VE E t td V V V d V d V f f e e S S R S R S T ST S T P P D H H V S S S S S R V V T R S R T S T R E E T T -=+-(37 如果发电机采用可控硅励磁系统,则可以忽略励磁时间常数 0e T ,那么控制量 U 就是励磁绕组电压 f E ,系统状态方程可以表示为 XAX BU =+ 形式的 三阶形式:00000E VV E E E d Vd V d e e f t t E E VE VE d V d V VV

14、d V S S R S S R T S T S T P P D E H H VV R S S S S S T R T R S R T S -=-+-(38 最优控制量 U kX =-4 最优控制设计算例发电机电气参数:隐极机纵轴电抗 2.5d X =, 纵轴暂态电抗 0.5d X =, 纵轴开环时间常数 010d T s =, 8H s =,综合阻尼系数 5.0D =; 变压器参数:0.01T X =; 线路参数:1L X =;系统参数:1s V =。所研究的发输电系统的示意图见图2。s 图 2 系统示意图设计最优运行点为 065=, 00.55P =时的最优控制器。 4.1 计算最优控制参数

15、1.9784sin 1sin 65e d q s P x E V +=由式(19可得:2sin d q d sq e d q s x x x VE P x x V -=-= 由式 (17可得:01.0487t V = 数据代入到相应公式,可以求得-0.16950.6564-0.1422-39.2699-0.6250 0-0.1064-0.0545-0.0893A =; 0.071900.0452B =选取 (1,100,5000Q diag =和 1R =作为权矩阵,进而可以得到状态反馈矩阵:46.5370-6.156068.76(, , 1, 7K lqr A B Q R =可知最优控制量

16、58.74805.176730.0717f e t u E KX P V =-=-+-。 上述计 算过程使用 MATLAB 完成,见报告附录中的程序 1。 4.2 验证最优励磁控制系统性能由于励磁控制系统参数是采用线性化状态方程而得到的,而实际系统的状 态方程是非线性的, 所以使用 MATLAB 的 ode45算法, 对单机系统中发电机的 运动微分方程进行计算, 通过观察发电机的转子角和输出电磁功率的变化情况,得到最优励磁控制系统的控制性能。单机系统的发电机运动方程如下:0=-2000( (sin sin 2 2q s dq s m d dq E V x x D V P H H H x x x

17、 -=-+00111cos d dq f qs d d d dx x E E E V T T T x -=-+ 此运动方程以转子角 、角速度 和 q 轴暂态电势 qE 作为三个状态量。微 分方程中, 02100f =; 0.55m P =; f E 作为输入量,引入最优状态反馈矩 阵,即用三个状态量进行表示,为0000000( (58.74805.176746.5370( 6. 30.0711560( 68.7617(7 (f q f q q e t e t q t E E E E KX E P P P V V V E =+=+-+-=+-+-=+-其中, e P 和 t V 作为与 qE 和

18、 有关的变量,分别用式(19和式(17表示。 4.2.1 静态稳定性系统在运行至 5s 时刻,加入一个时间长达 1秒的电磁功率的扰动,模拟负 荷的突然变化,即 5s 到 6s 之间 e P 从原来的 0.55增加到 0.6, 6s 之后 e P 回到最 初设计的稳态运行点。进行小扰动的仿真。仿真程序见报告附录中的程序 27和 13。仿真的功角和电磁功率波形图如下: 图 3 静态稳定功角仿真图从波形图可以看出,在发生小扰动后,采用最优控制理论所设计的励磁系统,能够快速地平息系统的振荡,转子角和电磁功率均能够迅速地回到最初设 计的最优运行点,提高了电力系统的静态稳定性。可见,发电机加入最优励磁 控

19、制器可以有效地改善系统稳定性。4.2.2 暂态稳定性系统运行至 5s 的时刻,在其中一条线路首端即图 2中的 d 处,发生三相短 路故障,故障 0.1s 后保护动作切除故障线路。程序见报告附录中的程序 2、 3和 812。该过程中的功角变化波形图如下: 图 4 暂态稳定功角仿真图由图 4可以看出, 系统在 5s 时刻发生三相短路后功角和电磁功率急剧上升, 5.1s 切除故障线路后系统迅速地恢复稳定,整个过程不超过 5秒钟。因此,最 优励磁控制器也可以很好地维持电力系统的暂态稳定性。5 结论本报告中先介绍了最优控制理论的相关原理,然后通过对与电网并联的发 电机基本方程组的研究从而建立了发电机与励

20、磁系统的线性化状态方程,采用 MATLAB 计算最优状态反馈参数并代入到系统的非线性微分方程组, 然后对系 统的非线性微分方程组进行了时域仿真。仿真结果表明,最优励磁控制器可以 很好地提高系统的静态稳定性和暂态稳定性。可见,最优励磁控制器有效地增 加了发电机励磁系统的阻尼,抑制了系统受到扰动时的振荡,加快系统恢复稳 定的速度。附 录程序 1:%输入发电机与系统的相关参数Xd=2.5;Xdz=0.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+

21、XL;Xdzs=Xdz+XT+XL;Xqs=Xq+XT+XL;Eq=Pe*Xds/(Vs*sin(angEqz=(Pe-0.5*Vs2*(Xdzs-Xqs/(Xdzs*Xqs*sin(2*ang*Xdzs/(Vs*sin(angTd=Td0*Xdzs/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+(

22、Vs*cos(ang*Xd2+2*Xe*Xd*Eq*Vs*cos(ang/Xds2+(Vs*sin(ang* Xq/Xqs2dVtEq=(Eq*Xe2+Xe*Xd*Vs*cos(ang/(Xds2*(Eq2*Xe2+Vs2*(cos(ang2*Xd2+2 *Xe*Xd*Eq*Vs*cos(ang/(Xds2+(Vs*sin(ang*Xq/Xqs2(-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*Xd2+2*Xd*Xe*Eq*Vs

23、*cos(ang/Xds2+(Vs*sin(an g*Xq/Xqs2(-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/H -D/H 0;(SE-SV/(Td*RV*SV (SEz-SV/RV -SE/(Td*SVB=REz/Td0;0;REz/(Td0*RVK=lqr(A,B,Q,R程序 2:%故障前或扰动前的系统状态function dy=gz1(t,y;XT=0.01;XL=0.75;Xe=

24、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(2-2*pi*50-2*pi*50/8*Pe;dy(3=(

25、-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;程序 3:%计算故障前或扰动前的电磁功率function Pe1=dcgl1(yPe1=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/Xdz

26、s+0.5*Vs2*(Xdzs-Xds/(Xdzs*Xds*sin(2*y(i,1;End程序 4:%某线路首端三相短路故障中的系统状态function dy=gz2(t,y;XT=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.5

27、5/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/Td 0*(Xd-Xdz/Xdzs*Vs*cos(y(1;dy=dy(1;dy(2;dy(3;程序 5:%计算故障中的电磁功率function Pe4=dcgl4(yPe4=zeros(length(y(:,1,1;for i=1:length(y(:,1;Pe4(i=0;End程序 6:%切除故障线路后的系统状态function dy=gz3(t,y;XT=0.01;XL

28、=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-2*pi*50-2*pi*50/8*Pe;dy(3=(-46.537*(Pe-0.55+6

29、.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(yPe5=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/(Xd

30、zs*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*(Xdzs-Xds/(Xdzs*Xds*cos(y(1;Vt=1/Xds*sqrt(Eq2*Xe2+Vs2*Xd2+2*Xe*Xd*Eq*Vs*cos(y(1;Td0=1

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

32、e=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/(Xdzs*Xds*sin(2*y

33、(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/Td+1/Td 0*(Xd-Xdz/Xdzs*Vs*cos(y(1;dy=dy(1;dy(2;dy(3;程序 11:%计算扰动后的电磁功率function Pe3=dcgl3(yPe3=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