基于广义预测控制的水轮机系统频率调节研究

宋运忠,冉皓

(河南理工大学 电气工程与自动化学院，河南 焦作 454000)

0 引 言

电力系统中，频率是衡量电能质量的一个重要指标，频率波动不仅会影响用户设备，还会影响发电厂设备的安全稳定运行，因此，频率质量的改善对经济效益可产生重要影响。我国城市电网的谐波污染问题比较严重，并呈现扩大化趋势，具有明显的地域性和渗透性等特点[1]。各行业和居民用户都对电能质量提出了更高要求，这需要从发电端解决电能质量问题。

ZHOU X等[2]，HOU G等[3]分别将广义预测控制应用于气轮机，提出一种基于熵的聚类和子空间辨识的模糊建模方法识别燃气轮机系统模型，在成本函数中加入惩罚因子，实现对转子转速的极限保护，为本文将其应用于水轮机提供了一些思路。上述所提的二次调频由部分电厂承担，主要承担负荷增量。火电厂的可调容量不大[4]，水电厂可调容量大，且调节速度快，枯水期最适合调频。由于水电厂调频的各种优点，故本文选择水电厂作为调频对象。李少远等[5]针对水轮机的特点提出一套按能量守恒准则启机的优化控制策略，采用的设计理论观点与预测控制一致。

本文主要研究连续时间线性系统广义预测控制(GPC)问题，TIAN Y等[6]基于模糊模型的广义预测控制(FGPC)方法，得到了变换控制的自回归积分移动平均模型；ZHANG R等[7]采用T-S模型作为预测模型，该方法既继承了广义预测控制的优点，又能有效控制复杂的非线性系统，由于模糊模型过于复杂，导致增加了大量计算和控制时间；李少远[8]讨论了由连续时间系统采样得到的离散时间系统组成的计算机控制系统采样周期的选择，指导了控制系统的设计；M.Alhajeri等[9]分析了预测控制的调整准则，主要是为了设计广义预测控制；S.KAYAIVIZHI等[10]提出了一种应用于电网的模糊自适应模型预测控制，对光伏、风能、燃料电池等分布式能源的负荷频率进行控制，但目前我国分布式能源并未广泛使用，与国情不符；HOU G等[11]则将广义预测控制用于大型水电站汽轮机模型，提高联合循环发电厂燃气轮机的控制性能和降低发电过程中的能耗，以节能为目标限制燃油，导致发电效率降低；M.BEUS等[12]将干扰和不确定性的弹性SISO模型应用于水电厂的线性化，将约束GPC控制器与经典PI控制器进行了响应比较，但由于预测控制实现的复杂性和工业应用的困难，目前还没有实现涡轮调速器应用于预测控制，需要进一步研究这种类型的控制器是否有实际实现的潜力；S.Olaru等[13]从水电厂的基本原理出发，以水轮机的控制为目标，对水电厂进行了建模。

本文提出基于广义预测控制的水轮机系统频率控制策略，首次将广义预测控制应用于水轮机系统。与文献[14]不同的是，不再以单一的调速器为对象，而是以整个水轮机调节系统作为控制对象进行预测控制。首先，分析基于现代控制理论连续系统的传递函数模型；然后，结合GPC控制算法的滚动优化、反馈校正和有限时域步进迭代算法对搭建的数学模型进行仿真求解，得到该系统的最优控制策略；最后，利用仿真算例，验证GPC在水轮机调节系统中的应用效果。

1 水电机组频率-功率调节系统

水电机组的功率和频率是同一个控制回路，图1为微机调速器工作原理图[15]，根据电网频率发生的改变，生成Δf，然后根据控制算法调节导叶开度。如果能量管理系统命令改变，则生成ΔP，也通过控制算法调节导叶开度，改变机组发电功率，最终使机组的发电功率与功率指令相符，实现二次调频。

图1 微机调速器工作原理图Fig.1 Working principle diagram of microcomputer governor

当调节系统进入稳定状态时，输入E=0，得到式(1)，式(1)可转化为式(2)。式(1)中ep是一次调频的不等式表达式。将式(2)转化为式(3)[16]，得到功率指令，根据AGC指令和频率偏差调整机组的功率，实现调频。

epΔP+Δf=0，

(1)

(2)

(3)

为了保证过渡过程中水轮机的转速升高和有压过水系统的压力升高不超过设计值，导叶运动速度应有一定限制，不能超过导叶最大运动速度。另外，导叶的控制指令应设置一定限幅，使机组运行在安全范围内[17]。

2 基于受控自回归积分滑动平均模 型的广义预测控制

目前各种预测控制算法的基本特征为预测模型、滚动优化和反馈校正[18]。

2.1 预测模型

考虑如下线性离散时间系统的状态空间模型：

(4)

式中：u(k)为输入；y(k)为k时刻状态方程的输出；ξ(k)为白噪声，其方差为σ2，均值为0；d=[τ/T]+1，为系统时滞系数，[τ/T]为不超过τ/T的最大整数，τ为下一通信延时预测值，T为离散系统的采样时间；z-1为后移算子，Δ=1-z-1为差分算子[19]。

2.2 输出预测

根据多步Diophantine方程，可求得多步最优预测输出[20-22]。取N1=1,N2=Nu=N，根据输出方程，可以推导出以下方程，

y(k+1)=F1Δu(k)+G1y(k)+E1ξ(k+1)=

f1,0Δu(k)+f1,1Δu(k-1)+…+

f1,nbΔu(k-nb)+G1y(k)+ξ(k+1)，

y(k+2)=F2Δu(k+1)+G2y(k)+E2ξ(k+2)=

f2,0Δu(k+1)+f2,1Δu(k)+f2,2Δu(k-1)+

…+f1,nb+1Δu(k-nb)+G2y(k)+ξ(k+

2)+e2,1ξ(k+1)，

⋮

y(k+N)=FNΔu(k+N-1)+GNy(k)+

ENξ(k+N)=fN,0Δu(k+N-1)+…+

fN,N-1Δu(k)+fN,NΔu(k-1)+…+

fN,nb+N-1Δu(k-nb)+GNy(k)+ξ(k+N)+

eN,1ξ(k+N-1)+…+eN,N-1ξ(k+1)。 (5)

对系统未来预测的输出可以写为如下矩阵形式：

Y=F1ΔU+F2ΔU(k-j)+GY(k)+Eξ=F1ΔU+Y1+Eξ，

(6)

式中，Y1=F2ΔU(k-j)+GY(k)，为基于过去与当前输入输出的输出预测响应，

Y=[y(k+1),y(k+2),…,y(k+N)]T

为未来的预测输出；

ΔU=[Δu(k),Δu(k+1),…,Δu(k+N-1)]T，

为当前和未来的控制增量向量；

ΔU(k-j)=[Δu(k-1),Δu(k-2),…,Δu(k-nb)]T为过去的控制增量向量；

Y(k)=[y(k),y(k-1),…,y(k-na)]T

为当前及过去的实际输出；

ξ=[ξ(k+1),ξ(k+2),…,ξ(k+N)]T

为未来的白噪声向量；

(7)

(8)

(9)

(10)

2.3 性能指标函数

预测控制的本质是优化控制算法，性能指标的优化决定控制效果。性能指标与系统的未来行为有关。目标输出为最小的方差跟踪预期轨迹。本文采用如下性能指标:

(11)

式中：y(k+j)和yr(k+j)分别为系统k+j时刻的实际输出和估计输出；N1为起始时刻，N2为终止时刻；Nu为控制长度；γj为加权。

由式(11)可知，优化不是离线一次完成的，而是在线反复运行完成的。在每个采样时间段，优化后的性能指标只涉及未来的有限时域，在下一个采样时间，优化后的时域向后移动，实现滚动优化[23]。

文献[24-25]给出了广义预测控制的控制律，可由参考信号与输出预测误差、控制输入误差的平方和计算出的性能指标推导。通过添加带有设计参数的新信号，对GPC进行扩展。

实际运行过程中不断更新信息，并始终根据最新信息进行优化，可有效解决模型失配、时变和干扰等引起的变化，使控制具有更强的鲁棒性，在运行中保持最优[26]。

性能指标函数式(11)可表示为矩阵形式，即

J=E{[Y-Yr]T[Y-Yr]+ΔUTΓΔU}，

(12)

式中：Γ=diag(γ1,γ2,…,γNu)；

Y=[y(k+N1),y(k+N1+1),…,y(k+N2)]T；Yr=[yr(k+N1),yr(k+N1+1),…,yr(k+N2)]T；ΔU=[Δu(k),Δu(k+1),…,Δu(k+Nu-1)]T。 以往的控制加权矩阵Γ为Nu阶单位阵，考虑控制与输出的权重，需要弱化控制加权系数在程序运行以后的权重，因此令控制加权系数

(13)

将式(6)代入式(12)中，得

J=E{[F1ΔU+Y1+Eξ-Yr]T[F1ΔU+
Y1+Eξ-Yr]+ΔUTΓΔU}，

(14)

(15)

则当前时刻的控制量为

(16)

此外，本文采用带遗忘因子入的递推最小二乘法估计对象参数，此种辨识方法更适用于参数突变但不频繁和参数缓慢变化两种情况，如式(17)所示：

(17)

综上所述，具体实施步骤可归纳如下[27]：

(2)采样当前实际输出y(k)和参考轨迹输出yr(k+j);

(4)求解Diophantine方程(18)，得到Ej,Gj,Fj。

(5)构造向量Yr,ΔU(k-j),Y(k)及矩阵G,F1,F2。

(6)利用式(16)计算并实施u(k)。

(7)返回(2)，这时k变为k+1，继续循环。

2.4 多步Diophantin方程迭代算法

多步Diophantin方程模型的形式如下。

C(z-1)=A(z-1)Ej(z-1)+z-jGj(z-1),

Fj(z-1)=B(z-1)Ej(z-1),

(18)

式中：

A(z-1)=1+a1z-1+…+anaz-na；

B(z-1)=b0+b1z-1+…+bnbz-nb；

C(z-1)=1+c1z-1+…+cncz-nc；

Ej(z-1)=1+ej,1z-1+…+ej,nejz-nej；

Gj(z-1)=gj,0+gj,1z-1+…+gj,ngjz-ngj；

Fj(z-1)=fj,0+fj,1z-1+…+fj,nfjz-nfj；

degEj=j-1,degGj=na-1；degFj=nb+j-1，

(19)

j=1,2,3,…,N，

j=1时，有

C=AE1+z-1G1，

(20)

从而有

(21)

根据式(18)，可以写出第(j+1)项，

(22)

式(18)与式(22)相减，得

A(Ej+1-Ej)=z-j(Gj-z-1Gj+1)，

(23)

由于等号右边乘上z-j，所以右边最低次幂也从z-j开始，左边前(j-1)项全为零0，Ej+1与Ej的前(j-1)项系数必相等，即

z-1Gj+1=Gj-ej+1,jA，

(24)

展开，得

(25)

可得递推公式为

3 水轮机调节系统模型

如图2所示[15]，为了使调节系统具有良好的动态特性，需要应用自动控制理论分析水轮机系统。建立水轮机系统各部分的数学模型分析系统的动态特性[28-30]。

图2 水轮机转速自动调节系统方框图Fig.2 Block diagram of automatic adjustment system for turbine speed

3.1 调速器液压随动系统

液压随动系统将电子调节器发出的电信号转换并放大成相应的机械信号，从而控制导叶[31]。

配压阀是以机械位移作为输入以连续控制输出的液体压力和流量的装置，也称为液压放大器，起着功率放大的作用。如图3所示[15]，配压阀开口为矩形，阀心位移量ΔS，宽度W，由窗口与接力器油流量平衡，得

图3 配压阀开口示意图Fig.3 Schematic diagram of pressure valve opening

(27)

式中：接力器位移Y=Y0+ΔY；接力器的运动速度

(28)

可得

(29)

将式(29)变化为

(30)

由此可得主接力器传递函数

(31)

3.2 引水系统

由于水击的存在，水轮机系统动态性能会变差。水轮机引水系统如图4所示[15]。

图4 水轮机引水系统Fig.4 Water diversion system of hydropower unit

转速发生变化时，调节器将自动操作导叶开合，从而引起流量和流速变化。在水流的惯性作用下，阀门和水库的压力发生改变，这种现象称为水击。水击压力ΔH和流量Q等变量可以用以下模型表示，

(32)

对式(32)变形，得

(33)

H(s)=-TωsQ(s)，

(34)

于是引水系统传递函数为

(35)

3.3 水轮机数学模型

依据水轮机综合特性曲线求取水轮机的7个传递系数，分别为：力矩对导叶开度传递系数ey；流量对导叶开度传递系数eqy；力矩对转速传递系数ex；流量对转速传递系数eqx；力矩对水头传递系数eh；流量对水位传递系数eqh；发电机负载的调节系数eg。

令eqx=0，en=eg-ex，则可得水轮机的传递函数

(36)

3.4 发电机负载模型

发电机负载模型如图5所示[15]。图5中，G为发电机，Mg1～Mgn为负荷力矩，K1～Kn为断路器。

图5 发电机负荷力矩简化模型Fig.5 Simplified model of generator load moments

(37)

(38)

式中，

(39)

发电机及负载的传递函数为

(40)

综上所述，采用表1可求得GPC控制对象的传递函数模型，为

(41)

对式(41)变形，得

(42)

即

0.5s3y(s)+3.6s2y(s)+5.7sy(s)+y(s)=-su(s)+u(s)。

(43)

对式(43)两边进行拉普拉斯逆变换，得

(44)

对其离散化，令

(45)

把式(45)代入式(44)，可得水轮机系统的差分方程模型为

y(k)=0.102u(k-2)+1.469y(k-1)-0.52y(k-2)+0.051y(k-3)。

(46)

4 仿真结果与分析

为了使输出y(k)平滑地过渡到设定值ω,参考轨迹取为

(47)

式中，α为输出柔化系数，本文取α=0.7，设定值ω则分别取为阶跃与方波信号。结合表1所给参数数值，进行针对水轮机对象模型的GPC仿真。为了验证预测模型的有效性，通过给定不同的设定值，观察输出转速是否满足要求[32]。

表1 水电站机组主要等效模型参数设计Tab.1 Main equivalent model parameter design of hydropower generating units

图6 阶跃信号GPC算法输出波形Fig.6 Output waveforms of GPC algorithm with step signal

图6和图7分别对应不同ω的GPC算法,图8和图9分别对应于各自不同设定值的控制输出序列。由图6～9可以看出，输出值在设定值之间上下波动。系统采用GPC方法进行发电机转速调节时，各典型状态变量的超调量、调节时间、振荡次数、衰减速率等指标均取得了较优的控制效果，证明本文GPC 频率调节方法的有效性及其控制性能的优越性。

图7 方波信号GPC算法输出波形Fig.7 Output waveforms of GPC algorithm with square wave signal

图8 方波信号GPC算法输入控制波形Fig.8 Input control waveforms of GPC algorithm with square wave signal

图9 阶跃信号GPC算法输入控制波形Fig.9 Input control waveforms of GPC algorithm with step signal

系统工况出现大范围变化时，PID控制往往需针对不同工况选取不同的PID参数。虽然GPC中参数N1,Nu,γj,α的取值与系统动态过程相关，但这些参数中，N1取决于系统时延、Nu决定控制量变化倾向、N2决定系统闭环稳定和优化时间、rj用于限制控制能量、α改变期望的响应速度，这些参数均对工况变化的敏感性不高。

当系统矩阵的对象参数未知或变化时，可加入自适应控制策略，首先利用式(4)中的对象参数a1,a2,b0,b1，它们分别为A(z-1)的第二项与第三项系数、B(z-1)的第一项与第二项系数,然后再设计GPC控制律。如图10～11所示，通过辨识得到系统对象参数。

图10 方波信号对象参数估计结果Fig.10 Parameter estimation results of the square wave signal object

图11 阶跃信号对象参数估计结果Fig.11 Parameter estimation results of the step signal object

对比图9与图12 、图8与图13可知，控制输入可调区间更广，控制输入更加剧烈，并对比图6与图14、图7与图15，可知输出波形的各个波峰更小，输出更加平滑，稳定性更强。

图12 阶跃信号系统辨识GPC算法输入控制波形Fig.12 Input control waveforms of system identification of GPC algorithm with step signal

图13 方波信号系统辨识GPC算法输入控制波形Fig.13 Input control waveforms of system identification of GPC algorithm with square wave signal

图14 阶跃信号系统辨识GPC算法输出波形Fig.14 Output waveforms of system identification of GPC algorithm with step signal

图15 方波信号系统辨识GPC算法输出波形Fig.15 Output waveforms of system identification of GPC algorithm with square wave signal

由于GPC的本质是通过在线优化计算出最优控制律，因此，采用GPC进行机组频率调节时，其控制参数具有更强的工况适用性。

从仿真结果可以看出，本文提出的方法使频率调节在水轮机系统模型中具有良好的动态响应性能，从而验证了本文提出控制方法的可行性。

5 结 语

本文对水轮机调速系统广义预测控制进行了研究，针对水电厂微机调节系统，提出了一种基于广义预测控制的水轮机系统频率调节方法，取代传统PID型调速器，以获得更为满意的控制性能。首先，建立了全面反映机组运行过程中调速系统关键状态变量变化过程的水轮机组调速系统模型；其次，根据该模型进行拉普拉斯变换，得出用于GPC的连续系统模型；最后，通过 GPC算法对该模型实施了在线滚动优化。

仿真结果表明，对于不同的设定值，该控制方法都能使水轮机转速调节保持良好的动态响应性能，并且增强了系统的鲁棒性。