CN104656448A - 一种基于解耦和扰动观测的超临界机组预测控制方法 - Google Patents

Abstract

本发明公开了一种基于解耦和扰动观测的超临界机组预测控制方法，该方法以超临界机火电机组为被控对象，以燃料量、给水量、汽机调门开度为输入，主蒸汽压力、中间点温度、机组负荷为输出，首先通过前馈补偿器使得被控对象近似解耦，然后分别通过在每个通道设置的扰动观测器对干扰进行观测，以实现对被控对象进行多变量预测控制，解决超临界机组协调控制系统因干扰影响所带来的控制效果差的问题，能够有效地抑制外部不可测量扰动的影响，同时抑制输入变量之间耦合引起的内部扰动，从而提高超临界机组协调控制系统的控制性能。

Description

一种基于解耦和扰动观测的超临界机组预测控制方法

技术领域

本发明属于超临界火电机组热工自动控制领域，尤其涉及一种基于解耦和扰动观测的超临界机组预测控制方法。

背景技术

在超临界火电机组的运行过程中，协调控制系统受到模型失配、变量耦合等引起的内部扰动，以及多种外部扰动的影响，这些扰动的存在明显降低了控制器的控制效果。

很多先进控制算法在处理强扰动时，因其在控制器中没有针对性的抗扰设计，系统的动态特性仍存在不足，影响控制性能。由于协调控制系统中除了可测扰动外还存在大量不可在线测量的扰动，为了改进控制器的控制效果，必须尽量消除不可测扰动的不利影响。不可测扰动的影响无法通过前馈控制的方式进行消除，此时，扰动观测器DOB提供了一种可行的解决办法。

然而，超临界机组协调控制系统具有多个输入和多个输出，将扰动观测器DOB应用到超临界机组协调控制系统时，存在以下几个问题：1、多变量系统输入与输出之间关系复杂，如果对每个通道都加入扰动观测器，则系统结构将变得相当复杂，不易于调试和维护；2、如果只针对主对角线通道设计扰动观测器，则系统存在强耦合的情况下，扰动观测器无法实现扰动的有效补偿。

发明内容

发明目的：为了解决现有技术中超临界机组协调控制系统因干扰所带来的控制效果差的问题，本发明提供了一种基于解耦和扰动观测的超临界机组预测控制方法，能够有效地抑制现场中各种不可测量的扰动的影响，进一步提高协调控制系统的控制性能，该方法只针对主对角线通道设计扰动观测器，降低了控制系统结构的复杂度，有效抑制了系统内因强耦合而带来的内部扰动。

技术方案：本发明提供的基于解耦和扰动观测的超临界机组预测控制方法，该方法包括以下步骤：

(1)利用多变量动态矩阵控制器采用滚动优化方法，对于输出，根据某个稳态工况下k时刻所检测到的实际输出值y_i(k)计算k时刻的初始输出预测值和k+1时刻输出预测值并利用所述输出预测值与给定值w_i(k)比较并计算得到最优控制序列ΔU_j,M(k)，再以所述最优控制序列ΔU_j,M(k)的第一个元素Δu_j(k)作为k时刻的控制增量计算k时刻输入的预输入值μ_j(k)，其中，i＝1,2,3,j＝1,2,3；

(2)利用扰动观测器对扰动进行在线观测，当扰动发生变化时将解扰预输入υ_j与超临界机组预测控制系统的集总扰动D₁(s)作为输入参量通过扰动观测器DOB_j观测出等效的干扰信号再利用所述等效的干扰信号作为补偿对解扰预输入υ_j进行实时更新；

(3)通过前馈补偿器D(s)对所述解扰预输入υ_j进行解耦得到k时刻的实际输入u_j并将所述实际输入u_j作用于传递函数模型G(s)所表征的被控对象；

(4)在k时刻的控制作用之后，检测k+1时刻的实际输出值y_i(k+1)并与所述k+1时刻输出预测值计算输出误差e_i(k+1)，利用所述输出误差e_i(k+1)对k+1时刻的初始输出预测值进行反馈校正，其中，i＝1,2,3；

(5)在之后周期内反复执行步骤(1)至(4)。

其中：在步骤(3)中，所述传递函数模型G(s)通过以下步骤建立：

1)在稳态工况下，将超临界机组协调控制系统切换到手动状态，针对三个输入分别做阶跃响应实验，获取每一个输出相对于每一个输入的阶跃响应曲线；

2)利用最小二乘法对所述阶跃响应曲线进行辨识，得到每个输出与每个输入之间的传递函数；

3)将所述传递函数作为矩阵元素得到传递函数模型G(s)。

其中，步骤(3)中所述前馈补偿器D(s)的表达式为：

$D\left(s\right)=\left[\begin{array}{ccc}1& D_{12}\left(s\right)& D_{13}\left(s\right)\\ D_{21}\left(s\right)& 1& D_{23}\left(s\right)\\ D_{31}\left(s\right)& D_{32}\left(s\right)& 1\end{array}\right]$

式中，

$D_{\mathrm{ij}}\left(s\right)=\frac{M_{\mathrm{ji}}\left(s\right)}{M_{\mathrm{jj}}\left(s\right)},i=\mathrm{1,2,3},j=\mathrm{1,2,3},$ M_ji(s)表示矩阵G(s)第j行第i列元素的代数余子式。

其中，步骤(2)中获得所述等效的干扰信号包括如下步骤：

将所述集总扰动D₁(s)通过Q(s)G_Dii ^-1(s)环节得到扰动估计值

将所述k-1时刻补偿后的预输入υ_j(k-1)通过低通滤波器Q_i(s)得到实际扰动等效值

利用所述扰动估计值减去所述实际扰动等效值得到所述等效的干扰信号

其中，在步骤(1)中，获取所述k时刻输入的预输入值μ_j(k)和所述k+1时刻输出预测值包括以下步骤：

1)设置所述多变量动态矩阵控制器的相关参数，包括采样时间T_s、预测时域P、控制时域M、模型时域N、误差校正矩阵H、输出误差权矩阵Q及控制权矩阵R；

2)采用预测模型对被控对象在预测时域P内的输出进行预测，所述模型的表达式为：

${\tilde{y}}_{\mathrm{PM}}\left(k\right)={\tilde{y}}_{\mathrm{PO}}\left(k\right)+{\mathrm{A\ΔU}}_M\left(k\right)$

式中，

${\tilde{y}}_{\mathrm{PM}}\left(k\right)=\left[\begin{array}{c}{\tilde{y}}_{1,\mathrm{PM}}\left(k\right)\\ .\\ .\\ .\\ {\tilde{y}}_{3,\mathrm{PM}}\left(k\right)\end{array}\right],{\tilde{y}}_{i,\mathrm{PM}}\left(k\right)=\left[\begin{array}{c}{\tilde{y}}_{i,M}(k+1|k)\\ .\\ .\\ .\\ {\tilde{y}}_{i,M}(k+P|k)\end{array}\right]$ ${\tilde{y}}_{P0}\left(k\right)=\left[\begin{array}{c}{\tilde{y}}_{1,P0}\left(k\right)\\ .\\ .\\ .\\ {\tilde{y}}_{3,P0}\left(k\right)\end{array}\right],{\tilde{y}}_{i,P0}\left(k\right)=\left[\begin{array}{c}{\tilde{y}}_{i,0}(k+1|k)\\ .\\ .\\ .\\ {\tilde{y}}_{i,0}(k+P|k)\end{array}\right]$

$A_{\mathrm{ij}}=\left[\begin{array}{ccc}{a}_{\mathrm{ij}}\left(1\right)& ...& ...\\ ...& ...& ...\\ a_{\mathrm{ij}}\left(M\right)& ...& a_{\mathrm{ij}}\left(1\right)\\ ...& ...& ...\\ a_{\mathrm{ij}}\left(P\right)& ...& a_{\mathrm{ij}}(P-M+1)\end{array}\right],A=\left[\begin{array}{ccc}{A}_{11}& ...& A_{13}\\ ...& ...& ...\\ A_{31}& ...& A_{33}\end{array}\right]$

${\mathrm{\ΔU}}_M\left(k\right)=\left[\begin{array}{c}{\mathrm{\ΔU}}_{1,M}\left(k\right)\\ .\\ .\\ .\\ {\mathrm{\ΔU}}_{3,M}\left(k\right)\end{array}\right],{\mathrm{\ΔU}}_{j,M}\left(k\right)=\left[\begin{array}{c}{\mathrm{\Δu}}_j\left(k\right)\\ .\\ .\\ .\\ {\mathrm{\Δu}}_j(k+M-1)\end{array}\right]$

P表示第i个输出在k时刻对未来k+j时刻的输出预测值，P表示第i个输出在k时刻对未来k+j时刻的初始输出预测值，a_ij为输出y_i对输入u_j的阶跃响应系数，其中，i＝1,2,3，j＝1,2,3；

3)采用滚动优化方法求解得到k时刻的最优控制序列ΔU_M(k)；

4)只执行采样时刻k第一个控制增量Δμ(k)，计算所述k时刻的预输入μ_j(k)的表达式如下：

μ_j(k)＝Δμ_j(k)+μ_j(k-1),j＝1,...,3；

5)利用所述预测模型计算所述k+1时刻输出预测值

其中，步骤(4)中，所述对k+1时刻的初始输出预测值进行反馈校正，包括以下步骤：

(1)k+1时刻的实际输出向量y(k+1)与k+1时刻输出的预测向量的输出误差向量e(k+1)为：

$e(k+1)=\left[\begin{array}{c}{e}_1(k+1)\\ .\\ .\\ .\\ e_3(k+1)\end{array}\right]=\left[\begin{array}{c}{y}_1(k+1)-{\tilde{y}}_{\mathrm{1,1}}(k+1|k)\\ .\\ .\\ .\\ y_3(k+1)-{\tilde{y}}_{\mathrm{3,1}}(k+1|k)\end{array}\right];$

(2)利用所述误差向量e(k+1)补偿所述k+1时刻输出的预测向量得到经校正的预测向量为：

${\tilde{y}}_{\mathrm{cor}}(k+1)={\tilde{y}}_{\mathrm{PM}}\left(k\right)+\mathrm{He}(k+1)$

式中，H为误差校正矩阵，表达式为：

(3)将所述经校正的预测向量进行移位得到k+1时刻初始输出预测向量为：

${\tilde{y}}_{P0}(k+1)=S_0{\tilde{y}}_{\mathrm{cor}}(k+1)$

式中，

有益效果：本发明与现有技术相比，具有以下显著优点：超临界机组一般被简化为三入三出的多变量系统，输入与输出之间关系复杂，如果对每个通道都加入扰动观测器，则系统结构将变得相当复杂，不易于调试和维护；本发明方法提供的基于解耦和扰动观测器的预测控制方法，无需对每个通道都加入扰动观测器，只针对主对角线通道设计扰动观测器，系统结构相对简单，易于调试。解耦后的广义对象虽不能实现完全解耦，考虑到扰动观测器并不要求多变量系统动态完全解耦，只要求非主对角线通道的增益较小，相对于主对角线通道可以忽略即可，这样主对角线上的扰动观测器便可将其他输入对本输出的影响视作不可测扰动；本发明方法在系统存在强耦合的情况下，由于前馈解耦的作用，可以有效地抑制输入变量之间耦合引起的内部扰动，还能有效地抑制控制系统的外部可测及不可测扰动、模型失配引起的内部扰动。

附图说明

图1是本发明的系统结构图；

图2是本发明的扰动观测器结构图；

图3是本发明的仿真效果对比图；图3(a)是输出变量主蒸汽压力响应曲线；图3(b)是控制变量汽机调门开度响应曲线；图3(c)是输出变量中间点温度响应曲线；图3(d)是控制变量给水量响应曲线；图3(e)是输出变量机组负荷响应曲线；图3(f)是控制变量燃料量响应曲线。

具体实施方式

下面结合附图和具体实施实例，进一步阐明本发明，应理解这些实施实例仅用于说明本发明而不用于限制本发明的范围，在阅读了本发明之后，本领域技术人员对本发明的各种等价形式的修改均落于本申请所附权利要求所限定的范围。

本实施例中的基于解耦和扰动观测的超临界机组预测控制方法的被控对象为超临界机组，可简化为一个三入三出的多变量系统，输入分别为燃料量B、给水量D、汽机调门开度U，输出分别为主蒸汽压力P、中间点温度T、机组负荷N。如图1所示，在该超临界机组预测控制系统中，包括：多变量动态矩阵控制器MDMC、设置在各自输出通道上的扰动观测器、前馈补偿器D(s)，多变量动态矩阵控制器MDMC用于根据当前检测时刻主蒸汽压力P的实际值y₁、中间点温度T的实际值y₂、机组负荷N的实际值y₃分别与各自的给定值w₁(k)、w₂(k)、w₃(k)得到预输入值μ₁、μ₂、μ₃，扰动观测器用于对扰动进行在线观测，一旦发生扰动则进行实时更新，μ₁、μ₂、μ₃分别经过扰动观测器DOB₁、DOB₂、DOB₃的补偿作用，得到补偿后的解扰预输入υ₁、υ₂、υ₃并输出；前馈补偿器D(s)用于对υ₁、υ₂、υ₃进行转换最终得到汽机调门开度U的实际输入u₁、给水量D的实际输入u₂、燃料量B的实际输出u₃，作用于该被控对象G(s)，图中d₁、d₂、d₃分别代表三个输出通道上的扰动信号。

在进行预测控制之前需要根据抽象模型即传递函数模型G(s)设计相应的前馈补偿器D(s)和扰动观测器DOB，然后再基于前三者所形成的带DOB的广义被控对象设计多变量动态矩阵控制器MDMC，具体如下：

(1)建立表征被控对象输入与输出之间响应关系的传递函数模型G(s)，具体如下：

在稳态工况下，将超临界机组协调控制系统切换到手动状态，针对三个输入变量分别做阶跃响应实验，获取超临界机组协调控制系统被控对象的传递函数模型，即：在燃料量B、给水量D不变的情况下，获取主蒸汽压力P、中间点温度T、机组负荷N三个输出对汽机调门开度U的阶跃响应曲线，同理，获取三个输出对其他两个输入的阶跃响应曲线。

利用最小二乘法辨识阶跃响应曲线得到控制对象的下列传递函数：机组负荷N对燃料量B传递函数G_BN(MW/t)，主蒸汽压力P对燃料量B传递函数G_BP(MPa/t)，中间点温度T对燃料量B的传递函数G_BT(℃/t)；机组负荷N对汽机调门开度U传递函数G_UN(MW/％)，主蒸汽压力P对汽机对调门开度U传递函数G_UP(MPa/％)，中间点温度T对汽机调门开度U传递函数G_UT(℃/％)；主蒸汽压力P对给水量D传递函数G_DP(MPa/t)，中间点温度T对给水量D传递函数G_DT(℃/t)，机组负荷N对给水量D传递函数G_DN(MW/t)，则输入与输出之间的矩阵表达式为：

$\left[\begin{array}{c}P\\ T\\ N\end{array}\right]=\left[\begin{array}{ccc}{G}_{\mathrm{UP}}\left(s\right)& G_{\mathrm{DP}}\left(s\right)& G_{\mathrm{BP}}\left(s\right)\\ G_{\mathrm{UT}}\left(s\right)& G_{\mathrm{DT}}\left(s\right)& G_{\mathrm{BT}}\left(s\right)\\ G_{\mathrm{UN}}\left(s\right)& G_{\mathrm{DN}}\left(s\right)& G_{\mathrm{BN}}\left(s\right)\end{array}\right]\left[\begin{array}{c}U\\ D\\ B\end{array}\right]---\left(1\right)$

辨识得到的传递函数矩阵用G(s)来表示则传递函数矩阵G(s)的表达式为：

$G\left(s\right)=\left[\begin{array}{ccc}{G}_{11}\left(s\right)& G_{12}\left(s\right)& G_{13}\left(s\right)\\ G_{21}\left(s\right)& G_{22}\left(s\right)& G_{23}\left(s\right)\\ G_{31}\left(s\right)& G_{32}\left(s\right)& G_{33}\left(s\right)\end{array}\right]=\left[\begin{array}{ccc}{G}_{\mathrm{UP}}\left(s\right)& G_{\mathrm{DP}}\left(s\right)& G_{\mathrm{BP}}\left(s\right)\\ G_{\mathrm{UT}}\left(s\right)& G_{\mathrm{DT}}\left(s\right)& G_{\mathrm{BT}}\left(s\right)\\ G_{\mathrm{UN}}\left(s\right)& G_{\mathrm{DN}}\left(s\right)& G_{\mathrm{BN}}\left(s\right)\end{array}\right]---\left(2\right)$

(2)在被控对象输入端设计一个前馈补偿器D(s)，使其与原被控对象的传递函数矩阵G(s)组成的广义被控对象的传递函数G_D(s)为对角阵，将被控对象近似为三个独立的SISO系统，从而实现输入输出间的近似解耦。具体步骤如下：

设前馈补偿器D(s)具有如下结构：

$D\left(s\right)=\left[\begin{array}{ccc}1& D_{12}\left(s\right)& D_{13}\left(s\right)\\ D_{21}\left(s\right)& 1& D_{23}\left(s\right)\\ D_{31}\left(s\right)& D_{32}\left(s\right)& 1\end{array}\right]---\left(3\right)$

解耦后的广义被控对象G_D(s)为对角阵，具有如下结构：

$G_D\left(s\right)=\left[\begin{array}{ccc}{G}_{D1}\left(s\right)& & \\ & G_{D2}\left(s\right)& \\ & & G_{D3}\left(s\right)\end{array}\right]---\left(4\right)$

广义被控对象G_D(s)与被控对象的传递函数矩阵G(s)、前馈补偿器D(s)的关系为：

G_D(s)＝G(s)D(s)          (5)

计算可得：

D(s)＝G^-1(s)G_D(s)         (6)

根据D(s)主对角线元素全为1，可知：

$G_{\mathrm{jj}}^{-1}\left(s\right)G_{\mathrm{Dj}}\left(s\right)=1---\left(7\right)$

结合式(5)则有：

$D_{\mathrm{ij}}\left(s\right)=G_{\mathrm{ij}}^{-1}\left(s\right)G_{\mathrm{Dj}}\left(s\right)=\frac{G_{\mathrm{ij}}^{-1}\left(s\right)}{G_{\mathrm{jj}}^{-1}\left(s\right)}---\left(8\right)$

由逆矩阵公式A^-1＝A^*/|A|可知：

$G_{\mathrm{ij}}^{-1}\left(s\right)=\frac{M_{\mathrm{ji}}\left(s\right)}{\left|G\left(s\right)\right|}---\left(9\right)$

其中，M_ji(s)表示矩阵G(s)第j行第i列元素的代数余子式。

因此，前馈解耦补偿器D(s)的解为：

$D_{\mathrm{ij}}\left(s\right)=\frac{M_{\mathrm{ji}}\left(s\right)}{M_{\mathrm{jj}}\left(s\right)}---\left(10\right)$

采用平衡截断降阶方法对计算得到的D(s)进行降阶处理，得到阶次较低的D(s)；根据式(5)计算得到解耦后的广义被控对象G_D(s)；采用平衡截断降阶方法对计算得到的G_D(s)进行降阶处理，得到阶次较低的G_D(s)。

(3)本实施例以第1个输入汽机调门开度U和第1个输出主蒸汽压力作为主控关系描述扰动观测器的具体设计。图2中，将外部扰动D_ex1、对象参数变化造成的实际对象G_n11(s)与抽象模型G_p11(s)的模型失配扰动D_m1(s)、其他两个输入对第1个输出主蒸汽压力的控制影响即G_p12(s)和G_p13(s)对G_p11(s)的耦合扰动D_c1(s)全部等效到输入端，即观测出等效的输入端干扰量并在控制中引入等量的补偿，实现对干扰的完全抑制。

系统的集总扰动D₁(s)包括外部扰动和模型失配及变量之间的耦合引起的内部扰动，即

D₁(s)＝D_ex1(s)+D_m1(s)+D_c1(s)          (11)

其中，D_ex1为外部扰动，D_m1为模型失配扰动，D_c1为耦合扰动。

忽略测量噪声ξ的影响，可以推出扰动估计值为：

${\hat{D}}_{f1}\left(s\right)=Q\left(s\right)G_{n11}^{-1}\left(s\right)D_1\left(s\right)---\left(12\right)$

则集总扰动D₁(s)和扰动估计值的误差为：

${\tilde{D}}_1\left(s\right)=D_1\left(s\right)-G_{n11}\left(s\right){\tilde{D}}_{f1}\left(s\right)=[1-Q\left(s\right)]D_1\left(s\right)---\left(13\right)$

根据终值定理可以推出

$\begin{array}{c}{\tilde{d}}_1(\∞)=\underset{s\→0}{\lim }s{\tilde{D}}_1\left(s\right)\\ =\underset{s\→0}{\lim }[1-Q\left(s\right)]\underset{s\→0}{\lim }s{D}_1\left(s\right)\\ =\underset{s\→0}{\lim }[1-Q\left(s\right)]d_1(\∞)\end{array}---\left(14\right)$

由式(14)可知，只要使得滤波器Q(s)的稳态增益为1，即可得到表明扰动观测器可以在线估计扰动。

因此，可以得到以下扰动观测器设计方法，具体步骤为：

①针对每个输出通道分别设计一个扰动观测器，即设计一个低通滤波器Q(s)和Q(s)G_Dii ^-1(s)，G_Dii(s)为G_D(s)主对角线上的传递函数，其中，i＝1,2,3；

②根据G_Dii ^-1(s)的阶次，计算低通滤波器Q(s)的阶次，来补偿G_Dii ^-1(s)的相对阶，即保证Q(s)G_Dii ^-1(s)物理上可实现；

③分别整定每个低通滤波器Q(s)的时间常数λ，保证能够有效的抑制扰动。

(4)计算带DOB的广义对象的阶跃响应模型，将其作为MDMC控制器的预测模型，得到各个模型的阶跃响应系数a_ij＝[a_ij(1)...a_ij(N)]^T，其中，i＝1,2,3，j＝1,2,3，N为阶跃响应的时域长度，N的选择应该确保使各输出响应值已接近稳态值。

利用上述设计好的前馈补偿器D(s)、扰动观测器DOB以及多变量动态矩阵控制器MDMC进行预测控制的具体步骤如下：

(1)利用多变量动态矩阵控制器采用滚动优化方法，根据某个稳态工况下k时刻所检测到的实际输出值y_i(k)计算k时刻的初始输出预测值和k+1时刻的输出预测值并利用输出预测值与给定值w_i(k)比较并计算得到最优控制序列ΔU_j,M(k)，以最优控制序列ΔU_j,M(k)的第一个元素Δu_j(k)作为k时刻的控制增量计算k时刻输入的预输入值μ_j(k)，其中，i＝1,2,3,j＝1,2,3。

具体步骤如下：

1)设置MDMC控制器的相关参数，包括采样时间T_s、预测时域P、控制时域M、模型时域N、误差校正矩阵H、输出误差权矩阵Q及控制权矩阵R。

2)超临界协调控制系统对象有3个输入和3个输出，假设每一输出y_i对每一输入u_j的阶跃响应系数矩阵a_ij＝[a_ij(1)...a_ij(N)]^T，其中，i＝1,2,3，j＝1,2,3，N为模型阶跃响应长度。根据一定的输入增量，利用这个模型，就可以预测模型未来时刻的输出。

被控变量的模型预测输出可以表示为如下形式：

${\tilde{y}}_{\mathrm{PM}}\left(k\right)={\tilde{y}}_{P0}\left(k\right)+\mathrm{A\Δ}{U}_M\left(k\right)---\left(15\right)$

式中，

${\tilde{y}}_{\mathrm{PM}}\left(k\right)=\left[\begin{array}{c}{\tilde{y}}_{1,\mathrm{PM}}\left(k\right)\\ .\\ .\\ .\\ {\tilde{y}}_{3,\mathrm{PM}}\left(k\right)\end{array}\right],{\tilde{y}}_{i,\mathrm{PM}}\left(k\right)=\left[\begin{array}{c}{\tilde{y}}_{i,M}(k+1|k)\\ .\\ .\\ .\\ {\tilde{y}}_{i,M}(k+P|k)\end{array}\right]$

${\tilde{y}}_{P0}\left(k\right)=\left[\begin{array}{c}{\tilde{y}}_{1,P0}\left(k\right)\\ .\\ .\\ .\\ {\tilde{y}}_{3,P0}\left(k\right)\end{array}\right],{\tilde{y}}_{i,P0}\left(k\right)=\left[\begin{array}{c}{\tilde{y}}_{i,0}(k+1|k)\\ .\\ .\\ .\\ {\tilde{y}}_{i,0}(k+P|k)\end{array}\right]$

$A_{\mathrm{ij}}=\left[\begin{array}{ccc}{a}_{\mathrm{ij}}\left(1\right)& ...& ...\\ ...& ...& ...\\ a_{\mathrm{ij}}\left(M\right)& ...& a_{\mathrm{ij}}\left(1\right)\\ ...& ...& ...\\ a_{\mathrm{ij}}\left(P\right)& ...& a_{\mathrm{ij}}(P-M+1)\end{array}\right],A=\left[\begin{array}{ccc}{A}_{11}& ...& A_{13}\\ ...& ...& ...\\ A_{31}& ...& A_{33}\end{array}\right]$

${\mathrm{\ΔU}}_M\left(k\right)=\left[\begin{array}{c}{\mathrm{\ΔU}}_{1,M}\left(k\right)\\ .\\ .\\ .\\ {\mathrm{\ΔU}}_{3,M}\left(k\right)\end{array}\right],{\mathrm{\ΔU}}_{j,M}\left(k\right)=\left[\begin{array}{c}{\mathrm{\Δu}}_j\left(k\right)\\ .\\ .\\ .\\ {\mathrm{\Δu}}_j(k+M-1)\end{array}\right]$

P表示第i个输出在k时刻对未来k+j时刻的预测值，P表示在k时刻对未来k+i时刻输出预测的初值，其中，i＝1,2,3，j＝1,2,3，表示施加控制作用时预测时域P各未来时刻的模型预测输出向量，表示无控制作用时预测输出初始向量，A是由模型的单位阶跃响应系数组成的动态矩阵，ΔU_M(k)表示各个输入依次有M个增量变化的增量矩阵。

3)MDMC采用滚动优化的策略，其目标函数具有如下形式：

$\underset{{\mathrm{\ΔU}}_M\left(k\right)}{\min }J\left(k\right)=[W\left(k\right)-{\tilde{y}}_{\mathrm{PM}}\left(k\right){]}^{T}Q[W\left(k\right)-{\tilde{y}}_{\mathrm{PM}}\left(k\right)]+{\mathrm{\ΔU}}_M{\left(k\right)}^T{\mathrm{R\ΔU}}_M\left(k\right)---\left(16\right)$

式中，W(k)表示未来P个采样周期输出的设定值，Q为输出误差权矩阵，R为控制权矩阵。

$W\left(k\right)=\left[\begin{array}{c}{w}_1\left(k\right)\\ .\\ .\\ .\\ w_3\left(k\right)\end{array}\right],w_i\left(k\right)=\left[\begin{array}{c}{w}_i(k+1)\\ .\\ .\\ .\\ w_i(k+P)\end{array}\right],i=1,...,3$

在不考虑约束的情况下，可求得最优的控制增量矩阵

${\mathrm{\ΔU}}_M\left(k\right)={(A^T\mathrm{QA}+R)}^{-1}{A}^{T}Q[W\left(k\right)-{\tilde{y}}_{P0}\left(k\right)]---\left(17\right)$

MDMC采用滚动优化算法，每个采样时刻只执行第一个控制增量Δμ(k)，到下一采样时刻，重新计算优化的控制增量矩阵。

4)计算当前时刻各输入的预输入值并输出。

μ_j(k)＝Δμ_j(k)+μ_j(k-1),j＝1,...,3          (18)

需要注意的是，此处计算得到的μ_j(k)并非实际的控制量输出。

5)利用根据式(15)算出预测模型计算所述k+1时刻输出的预测值

(2)利用扰动观测器对扰动进行在线观测，当扰动发生变化时将解扰预输入υ_j与超临界机组预测控制系统的集总扰动D₁(s)作为输入参量通过扰动观测器DOB_j观测出等效的干扰信号再利用所述等效的干扰信号作为补偿对解扰预输入υ_j进行实时更新；

具体步骤如下：

将集总扰动D₁(s)通过Q(s)G_Dii ^-1(s)环节得到扰动估计值

将第k-1时刻的经扰动补偿后的预输入值υ_j通过低通滤波器Q_i(s)得到实际扰动等效值

利用扰动估计值减去扰动等效值得到等效的干扰信号

(3)υ_j(k)再经过前馈补偿器D(s)的转换，最终得到实际的输入作用u_j(k)，并输出作用于对象G(s)得到每个输出的实际输出值y_i(k)；

(4)k时刻施加控制作用后，可根据式(15)算出对象在未来时刻的各输出值，到k+1时刻测得各实际输出y_i(k+1)后，与相应的预测值构成误差向量

$e(k+1)=\left[\begin{array}{c}{e}_1(k+1)\\ .\\ .\\ .\\ e_3(k+1)\end{array}\right]=\left[\begin{array}{c}{y}_1(k+1)-{\tilde{y}}_{\mathrm{1,1}}(k+1|k)\\ .\\ .\\ .\\ y_3(k+1)-{\tilde{y}}_{\mathrm{3,1}}(k+1|k)\end{array}\right]---\left(19\right)$

利用误差信息补偿基于模型的预测，可得到经校正的预测向量

${\tilde{y}}_{\mathrm{cor}}(k+1)={\tilde{y}}_{\mathrm{PM}}\left(k\right)+\mathrm{He}(k+1)---\left(20\right)$

式中，H为误差校正矩阵。

由于时间基点从k时刻移到k+1时刻，故校正后的预测向量可通过移位构成k+1时刻的初始预测值

${\tilde{y}}_{P0}(k+1)=S_0{\tilde{y}}_{\mathrm{cor}}(k+1)---\left(21\right)$

式中，

(5)将修正过的预测值作为k+1时刻的预测初值在之后的周期内，重复执行第(2)步到第(4)步。

针对上述具体实施例，做了如下仿真实验：

通过阶跃响应试验，利用最小二乘法分别得到：机组负荷对燃料量传递函数G_BN(MW/t)，主蒸汽压力对燃料量传递函数G_BP(MPa/t)，中间点温度对燃料量的传递函数G_BT(℃/t)；机组负荷对汽机调门开度传递函数G_UN(MW/％)，主蒸汽压力对汽机对调门开度传递函数G_UP(MPa/％)，中间点温度对汽机调门开度传递函数G_UT(℃/％)；机组负荷对给水量传递函数G_DN(MW/t)，主蒸汽压力对给水量传递函数G_DP(MPa/t)，中间点温度对给水量传递函数G_DT(℃/t)，如式(22)所示：

$\left[\begin{array}{c}P\\ T\\ N\end{array}\right]=\left[\begin{array}{ccc}\frac{-0.3}{26s+1}& \frac{0.009(1223s+1)}{(28s+1)(301s+1)}& \frac{0.17(206s+1)}{(68s+1)(406s+1)}\\ \frac{-7.9s}{{556s}^2+41s+1}& \frac{-4.1}{394s+1}& \frac{2.26}{402s+1}\\ \frac{60s}{(20s+1)(6.5s+1)}& \frac{409s-0.13}{(40s+1)(386s+1)}& \frac{3.7(123s+1)}{(71s+1)(410s+1)}\end{array}\right]\left[\begin{array}{c}U\\ D\\ B\end{array}\right]---\left(22\right)$

(1)针对得到的传递函数模型设计前馈补偿器D(s)，并降阶得到：

$D\left(s\right)=\left[\begin{array}{ccc}1& \frac{0.02s^2-0.0035s+5.5\×{10}^{-6}}{s^2+0.027s+0.00017}& \frac{{0.145s}^2+0.0073s+{3.16\×10}^{-5}}{s^2+0.024s+5.44\×{10}^{-5}}\\ \frac{{-9s}^2-0.09s+0.00015}{s^2+0.19s+0.0048}& 1& \frac{{0.34s}^2+0.029s+0.00074}{s^2+0.058s+0.00135}\\ \frac{{-14s}^2-0.3s+0.00015}{s^2+0.2s+0.01}& \frac{{-1.1s}^2-0.04s-6.3{\×10}^{-7}}{s^2+0.038s+0.0004}& 1\end{array}\right]---\left(23\right)$

(3)根据式(22)及式(23)计算广义对象的传递函数G_D(s)，并降阶得到G_D(s)实现近似解耦，

$G_D\left(s\right)=\left[\begin{array}{ccc}\frac{-{0.04s}^2-0.003s-9.2{\×10}^{-5}}{s^3+{0.23s}^2+0.015s+0.0003}& G_{D12}\left(s\right)& G_{D13}\left(s\right)\\ G_{D21}\left(s\right)& \frac{{-0.016s}^2-0.0005s-3.4\×{10}^{-6}}{s^3+{0.034s}^2+0.0004s+8{\×10}^{-7}}& G_{D23}\left(s\right)\\ G_{D31}\left(s\right)& G_{D32}\left(s\right)& \frac{{0.093s}^2+0.0052s+2.4{\×10}^{-5}}{s^3+{0.19s}^2+0.0032s+6.4{\×10}^{-6}}\end{array}\right]---\left(24\right)$

式中，

$G_{D12}\left(s\right)=\frac{{3\×10}^{-5}{s}^4-1.6{\×10}^{-6}{s}^3+1.6{\×10}^{-8}{s}^2-1.3{\×10}^{-10}s-3.4{\×10}^{-14}}{s^5+0.{035s}^4+{9\×10}^{-4}{s}^3+8{\×10}^{-6}{s}^2+4{\×10}^{-8}s+4.2\×{10}^{-11}}$

$G_{D13}\left(s\right)=\frac{{-3.3\×10}^{-6}{s}^5+7.5\×{10}^{-7}{s}^4-7.9{\×10}^{-9}{s}^3+4.2{\×10}^{-10}{s}^2+5.5{\×10}^{-14}s+7{\×10}^{-16}}{s^6+0.052s^5+{0.002s}^4+3.3{\×10}^{-5}{s}^3+3{\×10}^{-7}{s}^2+{9\×10}^{-10}s+1.1{\×10}^{-12}}$

$G_{D21}\left(s\right)=\frac{{6\×10}^{-4}{s}^4-1.4{\×10}^{-3}{s}^3+1.5{\×10}^{-5}{s}^2-{9\×10}^{-8}s-{4.7\×10}^{-9}}{s^5+{0.34s}^4+{0.05s}^3+0.002s^2+2.4{\×10}^{-5}s+5\×{10}^{-8}}$

$G_{D23}\left(s\right)=\frac{1.{8\×10}^{-5}{s}^4-1.4{\×10}^{-6}{s}^3+5.4\×{10}^{-9}{s}^2-3.2{\×10}^{-10}s+5{\×10}^{-13}}{s^5+{0.03s}^4+{0.001s}^3+1\×{10}^{-5}{s}^2+4{\×10}^{-8}s+4{\×10}^{-11}}$

$G_{D31}\left(s\right)=\frac{{1.3\×10}^{-4}{s}^5-0{.003s}^4-7{\×10}^{-5}{s}^3+4{\×10}^{-6}{s}^2-8.5{\×10}^{-8}s+1.5{\×10}^{-10}}{s^6+{0.5s}^5+{0.07s}^4+{0.006s}^3+1.4{\×10}^{-4}{s}^2+1.5\×{10}^{-6}s+2.8\×{10}^{-9}}$

$G_{D32}\left(s\right)=\frac{{-5\×10}^{-4}{s}^4+4.6{\×10}^{-5}{s}^3-1.7{\×10}^{-6}{s}^2+1.6{\×10}^{-8}s-8.5\×{10}^{-11}}{s^5+{0.07s}^4+0.003s^3+4.2\×{10}^{-5}{s}^2+3.2{\×10}^{-7}s+6{\×10}^{-10}}$

(4)分别设置三个通道扰动观测器的低通滤波器Q(s)的阶次和参数。阶次均取1时，即可满足要求，取参数λ＝5，即

$Q\left(s\right)=\frac{1}{5s+1}---\left(25\right)$

三个通道的低通滤波器Q(s)均为式(24)形式。

(5)在主蒸汽压力通道加入一个单位阶跃扰动信号，扰动通道模型为

$H_1\left(s\right)=\frac{1}{4.8s+1}---\left(26\right)$

(6)在每个控制周期内，计算MDMC控制器的输出，扰动观测器补偿MDMC控制作用并经解耦器D(s)转换后，最终得到实际的控制系统输入变量，并输出作用于对象。k时刻施加控制作用后，预测出对象在未来时刻的各输出值，到k+1时刻测得各实际输出y_i(k+1)后，与相应的预测值构成误差向量，利用误差信息校正模型的预测值。整个控制过程如此循环进行。

①设置MDMC控制器的相关参数，采样时间为T_s＝1s，模型长度N＝3000，预测时域P＝2500，控制时域M＝10；

②利用式(12)计算MDMC控制器的输出μ_j(k)；

$\left[\begin{array}{c}{\mathrm{\μ}}_U\\ {\mathrm{\μ}}_D\\ {\mathrm{\μ}}_B\end{array}\right]=\left[\begin{array}{cccccccc}0& 0.0075& 0.0201& 0.0361& 0.0540& 0.0729& ...& 1.1{\×10}^{-4}\\ 0& {-3.1\×10}^{-5}& {-1.3\×10}^{-4}& -3.9\×{10}^{-4}& {-8.8\×10}^{-4}& -0.0017& ...& {-6.5\×10}^{-7}\\ 0& {8.3\×10}^{-5}& {3.44\×10}^{-4}& 8.8{\×10}^{-4}& 0.0018& 0.0032& ...& {-2.3\×10}^{-5}\end{array}\right]---\left(27\right)$

③经过扰动观测器的补偿，得到补偿后的控制作用υ_j(k)并输出；

$\left[\begin{array}{c}{\mathrm{\υ}}_U\\ {\mathrm{\υ}}_D\\ {\mathrm{\υ}}_B\end{array}\right]=\left[\begin{array}{cccccccc}0& 0& 4{.7\×10}^{-14}& 0.0523& 0.1084& 0.1599& ...& 3.3570\\ 0& 0& {9.93\×10}^{-29}& 7.21{\×10}^{-5}& {7.82\×10}^{-5}& {-1.23\×10}^{-4}& ...& -0.0772\\ 0& 0& {-2.68\×10}^{-30}& {9.69\×10}^{-6}& {1.62\×10}^{-4}& {3.68\×10}^{-4}& ...& -0.0480\end{array}\right]---\left(28\right)$

④υ_j(k)再经过前馈补偿器D(s)的转换，最终得到实际的输入作用u_j(k)，并输出作用于对象。

$\left[\begin{array}{c}{u}_U\\ u_D\\ u_B\end{array}\right]=\left[\begin{array}{cccccccc}0& 0& {4.74\×10}^{-14}& 0.0523& 0.1084& 0.16& ...& 3.3337\\ 0& 0& {-4.26\×10}^{-13}& -0.4492& -0.8972& -1.2558& ...& 0.0011\\ 0& 0& -{6.63\×10}^{-13}& -0.6991& -1.3961& -1.9526& ...& 0.0020\end{array}\right]---\left(29\right)$

为了做对比，设计三种控制器分别为多变量动态矩阵控制器MDMC、扰动观测器与预测控制复合控制方案DOBDMC及本发明所述的基于前馈解耦的扰动观测器与预测控制复合控制方法DeDOBDMC，并比较其控制性能。图3(a)中，在扰动相同情况下，DeDOBDMC的主蒸汽压力输出通道的超调量为0.08，与DOBDMC的相差不大，都比DMC的要小很多；图3(c)中，DeDOBDMC的中间点温度输出通道耦合扰动带来影响几乎为0，而DOBDMC的影响虽然也比较小，但是相比DMC的影响还要稍大；同样地，图3(e)中，DeDOBDMC的机组负荷输出通道耦合扰动带来影响几乎为0，而DOBDMC的影响达到了0.52，甚至比DMC的影响还要大。图3(b)和图3(d)表明DeDOBDMC的汽机调门开度及给水量响应曲线与DOBDMC的相同；在图3(f)中，DOBDMC的燃料量响应曲线要比DeDOBDMC的剧烈很多。上述数据证明本发明提出的带前馈解耦的扰动观测器与预测控制器复合控制方案(DeDOBDMC)能够进一步改善控制效果。

Claims (7)

1.一种基于解耦和扰动观测的超临界机组预测控制方法，其特征在于，该方法包括以下步骤：

(1)利用多变量动态矩阵控制器采用滚动优化方法，对于输出，根据某个稳态工况下k时刻所检测到的实际输出值y_i(k)计算k时刻的初始输出预测值和k+1时刻输出预测值并利用所述输出预测值与给定值w_i(k)比较并计算得到最优控制序列ΔU_j,M(k)，再以所述最优控制序列ΔU_j,M(k)的第一个元素Δu_j(k)作为k时刻的控制增量计算k时刻输入的预输入值μ_j(k)，其中，i＝1,2,3,j＝1,2,3；

(2)利用扰动观测器对扰动进行在线观测，当扰动发生变化时将解扰预输入υ_j与超临界机组预测控制系统的集总扰动D₁(s)作为输入参量通过扰动观测器DOB_j观测出等效的干扰信号再利用所述等效的干扰信号作为补偿对解扰预输入υ_j进行实时更新；

(3)通过前馈补偿器D(s)对所述解扰预输入υ_j进行解耦得到k时刻的实际输入u_j并将所述实际输入u_j作用于传递函数模型G(s)所表征的被控对象；

(4)在k时刻的控制作用之后，检测k+1时刻的实际输出值y_i(k+1)并与所述k+1时刻输出预测值计算输出误差e_i(k+1)，利用所述输出误差e_i(k+1)对k+1时刻的初始输出预测值进行反馈校正，其中，i＝1,2,3；

(5)在之后周期内反复执行步骤(1)至(4)。

2.根据权利要求1所述的基于解耦和扰动观测的超临界机组预测控制方法，其特征在于，在步骤(3)中，所述传递函数模型G(s)通过以下步骤建立：

(1)在稳态工况下，将超临界机组协调控制系统切换到手动状态，针对三个输入分别做阶跃响应实验，获取每一个输出相对于每一个输入的阶跃响应曲线；

(2)利用最小二乘法对所述阶跃响应曲线进行辨识，得到每个输出与每个输入之间的传递函数；

(3)将所述传递函数作为矩阵元素得到传递函数模型G(s)。

3.根据权利要求2所述的基于解耦和扰动观测的超临界机组预测控制方法，其特征在于，所述前馈补偿器D(s)的表达式为：

$D\left(s\right)=\left[\begin{array}{ccc}1& D_{12}\left(s\right)& D_{13}\left(s\right)\\ D_{21}\left(s\right)& 1& D_{23}\left(s\right)\\ D_{31}\left(s\right)& D_{32}\left(s\right)& 1\end{array}\right]$

式中，

i＝1,2,3，j＝1,2,3，M_ji(s)表示矩阵G(s)第j行第i列元素的代数余子式。

4.根据权利要求3所述的基于解耦和扰动观测的超临界机组预测控制方法，其特征在于，所述等效的干扰信号通过以下步骤得到：

(1)将所述集总扰动D₁(s)通过环节得到扰动估计值

(2)将所述k-1时刻补偿后的预输入υ_j(k-1)通过低通滤波器Q_i(s)得到实际扰动等效值

(3)利用所述扰动估计值减去所述实际扰动等效值得到所述等效的干扰信号

5.根据权利要求1所述的基于解耦和扰动观测的超临界机组预测控制方法，其特征在于，在步骤(1)中，获取所述k时刻输入的预输入值μ_j(k)和所述k+1时刻输出预测值包括以下步骤：

(1)设置所述多变量动态矩阵控制器的相关参数，包括采样时间T_s、预测时域P、控制时域M、模型时域N、误差校正矩阵H、输出误差权矩阵Q及控制权矩阵R；

(2)采用预测模型对被控对象在预测时域P内的输出进行预测，所述模型的表达式为：

${\tilde{y}}_{\mathrm{PM}}\left(k\right)={\tilde{y}}_{P0}\left(k\right)+{\mathrm{A\ΔU}}_M\left(k\right)$

式中，

${\tilde{y}}_{\mathrm{PM}}\left(k\right)=\left[\begin{array}{c}{\tilde{y}}_{1,\mathrm{PM}}\left(k\right)\\ .\\ .\\ .\\ {\tilde{y}}_{3,\mathrm{PM}}\left(k\right)\end{array}\right],{\tilde{y}}_{i,\mathrm{PM}}\left(k\right)=\left[\begin{array}{c}{\tilde{y}}_{i,M}(k+1|k)\\ .\\ .\\ .\\ {\tilde{y}}_{i,M}(k+P|k)\end{array}\right]$

${\tilde{y}}_{P0}\left(k\right)=\left[\begin{array}{c}{\tilde{y}}_{1,P0}\left(k\right)\\ .\\ .\\ .\\ {\tilde{y}}_{3,P0}\left(k\right)\end{array}\right],{\tilde{y}}_{i,P0}\left(k\right)=\left[\begin{array}{c}{\tilde{y}}_{i,0}(k+1|k)\\ .\\ .\\ .\\ {\tilde{y}}_{i,0}(k+P|k)\end{array}\right]$

$A_{\mathrm{ij}}=\left[\begin{array}{ccc}{a}_{\mathrm{ij}}\left(1\right)& ...& ...\\ ...& ...& ...\\ a_{\mathrm{ij}}\left(M\right)& ...& a_{\mathrm{ij}}\left(1\right)\\ ...& ...& ...\\ a_{\mathrm{ij}}\left(P\right)& ...& a_{\mathrm{ij}}(P-M+1)\end{array}\right],A=\left[\begin{array}{ccc}{A}_{11}& ...& A_{13}\\ ...& ...& ...\\ A_{31}& ...& A_{33}\end{array}\right]$

${\mathrm{\ΔU}}_M\left(k\right)=\left[\begin{array}{c}{\mathrm{\ΔU}}_{1,M}\left(k\right)\\ .\\ .\\ .\\ {\mathrm{\ΔU}}_{3,M}\left(k\right)\end{array}\right],{\mathrm{\ΔU}}_{j,M}\left(k\right)=\left[\begin{array}{c}{\mathrm{\Δu}}_j\left(k\right)\\ .\\ .\\ .\\ {\mathrm{\Δu}}_j(k+M-1)\end{array}\right]$

i＝1,2,3,j＝1,…,P表示第i个输出在k时刻对未来k+j时刻的输出预测值，i＝1,2,3,j＝1,…,P表示第i个输出在k时刻对未来k+j时刻的初始输出预测值，a_ij为输出y_i对输入u_j的阶跃响应系数，其中，i＝1,2,3，j＝1,2,3；

(3)采用滚动优化方法求解得到k时刻的最优控制序列ΔU_M(k)；

(4)只执行采样时刻k第一个控制增量Δμ(k)，计算所述k时刻的预输入μ_j(k)的表达式如下：

μ_j(k)＝Δμ_j(k)+μ_j(k-1),j＝1,...,3；

(5)利用所述预测模型计算所述k+1时刻输出预测值

6.根据权利要求5所述的基于解耦和扰动观测的超临界机组预测控制方法，其特征在于，所述对k+1时刻的初始输出预测值进行反馈校正，包括以下步骤：

(1)k+1时刻的实际输出向量y(k+1)与k+1时刻输出的预测向量的输出误差向量e(k+1)为：

$e(k+1)=\left[\begin{array}{c}{e}_1(k+1)\\ .\\ .\\ .\\ e_3(k+1)\end{array}\right]=\left[\begin{array}{c}{y}_1(k+1)-{\tilde{y}}_{\mathrm{1,1}}(k+1|k)\\ .\\ .\\ .\\ y_3(k+1)-{\tilde{y}}_{\mathrm{3,1}}(k+1|k)\end{array}\right];$

(2)利用所述误差向量e(k+1)补偿所述k+1时刻输出的预测向量得到经校正的预测向量为：

${\tilde{y}}_{\mathrm{cor}}(k+1)={\tilde{y}}_{\mathrm{PM}}\left(k\right)+\mathrm{He}(k+1)$

式中，H为误差校正矩阵，表达式为：

(3)将所述经校正的预测向量进行移位得到k+1时刻初始输出预测向量为：

${\tilde{y}}_{P0}(k+1)=S_0{\tilde{y}}_{\mathrm{cor}}(k+1)$

式中，

7.根据权利要求1所述的基于解耦和扰动观测的超临界机组预测控制方法，其特征在于，所述输入分别为燃料量B、给水量D、汽机调门开度U，所述输出分别为主蒸汽压力P、中间点温度T、机组负荷N。