CN103322553B - 一种火电机组过热汽温多模型扰动估计预测控制方法 - Google Patents

Abstract

本发明公开了一种火电机组过热汽温多模型扰动估计预测控制方法，它以过热汽温系统为被控对象，减温喷水流量为系统控制输入量，过热汽温为系统输出，基于被控对象状态空间模型，建立多种扰动模型分别扩增系统状态，构成扰动模型集。扰动模型集中的各子模型通过加权平均，拟合不可测扰动对系统输出的影响。运用预测控制思想，每一步进行一次优化求解获得最佳系统控制输入量。本发明采用的多模型扰动估计预测控制方法，相比一般的预测控制算法能够更加有效地抑制不可测扰动，维持汽温稳定在设定值附近。当无不可测扰动时，本发明与普通预测控制算法等效，具有较好的跟踪调节性能。

Description

一种火电机组过热汽温多模型扰动估计预测控制方法

技术领域

本发明涉及一种过热汽温系统的优化控制方法，尤其涉及一种火电机组过热汽温抗不可测扰动的多模型扰动估计预测控制方法，属于热工自动控制技术。

背景技术

过热汽温直接关系到火电机组运行的安全性和经济性，是机组运行过程中需要重点监控的过程参数之一。

由于过热器管道较长，汽温被控对象存在较大的惯性和迟延，导致传统控制方法往往难以取得满意的控制效果。近年来预测控制算法在电站汽温控制系统应用中取得一定成效。然而基于模型的预测控制算法在应对诸如煤质变化引起的烟气侧不可测扰动时，由于缺乏对不可测扰动的建模，控制效果并不理想。目前大多数多模型预测控制器主要为了解决机组大范围工况变动的非线性，未考虑通过扰动建模主动抗干扰。

多模型扰动估计预测控制算法基于过热汽温对象状态空间模型建立多类扰动模型扩增系统状态，构成多个扩增系统模型。模型簇中的各子系统模型通过贝叶斯概率加权平均得到输出，能够拟合不可测扰动对系统输出的影响。

本发明充分运用预测控制思想，每一步进行一次优化求解获得最佳减温喷水控制输入量。仿真结果表明，本发明中的算法相比一般的预测控制算法能够更加有效地抑制不可测扰动，维持过热汽温稳定在设定值附近。当无不可测扰动时，该算法与普通预测控制算法等效，具有较好的设定值跟踪和调节性能。

发明内容

发明目的：为了克服现有技术中存在的不足，本发明提供一种火电机组过热汽温多模型扰动估计预测控制方法，能够有效抑制过程不可测扰动，提高过热汽温的调节品质。

技术方案：为实现上述目的，本发明采用的技术方案为：

一种火电机组过热汽温多模型扰动估计预测控制方法，包括如下步骤：

(1)在稳态工况下，将过热汽温系统切换到手动状态，以减温喷水流量为输入进行过热汽温的开环阶跃响应试验，利用两点法求出过热汽温对减温喷水流量的传递函数W(s)(℃/(kg/s))；

(2)根据线性系统理论的系统实现方法，将过热汽温对减温喷水流量的传递函数W(s)(℃/(kg/s))转换为系统状态空间模型 $\left\{\begin{array}{c}\stackrel{\·}{x}=\mathrm{Ax}+\mathrm{Bu}\\ y=\mathrm{Cx}+\mathrm{Du}\end{array}\right.,$ 其中[A,B,C,D]为系统矩阵，x为过热汽温系统的状态，y为过热汽温系统的输出，u为预测控制器的输出；

(3)建立n种类型扰动模型分别扩增过热汽温系统状态，得到n个扩增系统状态模型[ⁱA^a，ⁱB^a，ⁱC^a，ⁱD^a]，其中ⁱD^a＝0，i＝1,2,…,n，所述扩增系统状态空间模型的结构如式(1)和式(2)所示：

其中，表示k时刻扩增的状态向量，包括两个子向量ⁱx_k和 ⁱx_k表示k时刻过热汽温系统原有的状态向量，表示k时刻第i个扰动模型的扰动、扰动变化率或周期的状态；ω_k-1表示k-1时刻过程噪声，v_k表示k时刻输出噪声；

采用式(3)和式(4)简洁表达式(1)和式(2)为：

$x_k^a{}_i=A_a^{i}x_{k-1}^a{}_i+B_a^i{u}_{k-1}+\mathrm{\Ω}_a^i{\mathrm{\ω}}_{k-1}---\left(3\right)$

$y_k^i=C_a^{i}x_k^a{}_i+v_k---\left(4\right)$

其中代表 $\left[\begin{array}{c}x_k^i\\ d_k^a{}_i\end{array}\right],$ ⁱA^a代表 $\left[\begin{array}{cc}A_1^i& A_2^i\\ A_3^i& A_4^i\end{array}\right],$ ⁱB^a代表 $\left[\begin{array}{c}B_1^i\\ B_2^i\end{array}\right],$ ⁱC^a代表[ⁱC₁ ⁱC₂]；

(4)采用卡尔曼滤波方法对n个扩增系统状态模型分别进行状态估计，如式(5)、式(6)和式(7)：

$\hat{x}_{k|k-1}^a{}_i=A_a^i\hat{x}_{k-1|k+1}^a{}_i+B_a^i{u}_{k-1}---\left(5\right)$

$\hat{x}_{k|k}^a{}_i=\hat{x}_{k|k-1}^a{}_i+L_k^i(y_k-C_a^i\hat{x}_{k|k-1}^a{}_i)---\left(6\right)$

$\hat{y}_{k|k-1}^i=C_a^i\hat{x}_{k|k-1}^a{}_i---\left(7\right)$

在状态未知情况下，首先使用式(5)和式(7)预估状态和输出再通过式(6)由实际输出y_k和预估输出间的误差来矫正过热汽温系统状态得到新的预估状态和输出状态估计增益矩阵ⁱL_k通过求解Riccati方程获得，Riccati方程中的Qi代表输入噪声项的协方差，Ri代表输出噪声项的协方差，卡尔曼滤波参数Qi/Ri为预测控制器的一个调节参数；

(5)采用贝叶斯概率加权方法对n个扩增系统状态模型进行加权平均，获得过热汽温的P步平均预测值，应用预测控制算法计算最佳减温喷水流量，将计算得到的最佳减温喷水流量作为过热汽温系统的控制输入量，具体包括如下步骤：

(51)根据预估输出通过式(8)计算实际输出y_k与预估输出间的偏差ⁱε_k，并根据偏差ⁱε_k计算各个扩增系统状态模型的权重ⁱw_k：

$\mathrm{\ϵ}_k^i=y_k-\hat{y}_{k|k-1}^i---\left(8\right)$

(52)通过递推贝叶斯公式计算各扰动模型与过热汽温系统的匹配概率，通过上一时刻的匹配概率ⁱρ_k-1递推获得当前时刻的匹配概率ⁱρ_k，如式(9)：

$\mathrm{\ρ}_k^i=\frac{\exp (-\frac{1}{2}\mathrm{\ϵ}_k^T{}_i\mathrm{\Λ}_i\mathrm{\ϵ}_k^i)\mathrm{\ρ}_{k-1}^i}{{\mathrm{\Σ}}_{j=1}^4\exp (-\frac{1}{2}\mathrm{\ϵ}_k^T{}_j\mathrm{\Λ}_i\mathrm{\ϵ}_k^j)\mathrm{\ρ}_{k-1}^j}---\left(9\right)$

设定下限值δ，当上一时刻的匹配概率ⁱρ_k-1低于下限值δ时，以下限值δ置换ⁱρ_k-1后带入式(9)进行ⁱρ_k计算；贝叶斯概率加权矩阵ⁱΛ为预测控制器的一个调节参数；

(53)对各扰动模型与热汽温系统的匹配概率进行规格化处理，如式(10)：

$w_k^i=\left\{\begin{array}{cc}\frac{\mathrm{\&rho;}_k^i}{{\mathrm{\&Sigma;}}_{j=1}^4\mathrm{\&rho;}_k^j}& \mathrm{\&rho;}_k^i>\mathrm{\&delta;}\\ 0& \mathrm{\&rho;}_k^i<\mathrm{\&delta;}\end{array}\right.---\left(10\right)$

(55)根据式(11)计算过热汽温的P步平均预测值：

${\stackrel{\&OverBar;}{y}}_{k+j|k}={\mathrm{\&Sigma;}}_{i=1}^{4}w_k^i\hat{y}_{k+j|k}^i---\left(11\right)$

其中表示第i个扩增系统状态模型的预测输出值，ⁱw_k表示第i个扩增系统状态模型的权值，为所有扩增系统状态模型的平均预测输出值，j＝1,2,…,P；

(56)计算减温喷水流量控制量增量Δu_k，取式(12)性能指标函数式：

其中，Y_sp为过热汽温设定值，为过热汽温的P步平均预测矩阵， $\stackrel{\&OverBar;}{Y}={[{\stackrel{\&OverBar;}{y}}_{k+1|k},{\stackrel{\&OverBar;}{y}}_{k+2|k},\&CenterDot;\&CenterDot;\&CenterDot;{,\stackrel{\&OverBar;}{y}}_{k+P|k}]}^T,\mathrm{\&Delta;U}={[{\mathrm{\&Delta;u}}_k,{\mathrm{\&Delta;u}}_{k+1},\&CenterDot;\&CenterDot;\&CenterDot;,{\mathrm{\&Delta;u}}_{k+m-1}]}^T,$ W_y为误差权矩阵和W_u为控制量权矩阵；

每一采样时刻进行一次式(12)的优化计算，获得m时域的最佳系统控制输入增量ΔU_opt＝[Δu_k，Δu_k+1，…，Δu_k+m-1]^T；只以当前时刻的控制输入增量Δu_k对当前过热汽温系统实施控制，计算减温喷水流量控制量u_k，如式(13)：

u_k＝u_k-1+Δu_k   (13)

(6)设置预测控制器的调节参数，包括采样时间Ts、预测步数P、减温喷水流量控制步数M、输出误差权矩阵Wy、控制权矩阵Wu、n个扩增系统状态模型对应的卡尔曼滤波参数Qi/Ri和贝叶斯概率加权矩阵Λi，i＝1,2,…,n；

(7)向过热汽温系统输入减温喷水流量控制量u_k，获得热汽温系统实际输出y_k。

重复步骤(4)～(7)实现连续控制。

所述步骤(6)中，Ts可以用经验规则T95/Ts＝5～15来选取，其中，T95为过渡过程上升到95％的调节时间；P选为近似等于过热汽温阶跃响应的上升时间；Qi/Ri一般根据扰动模型形式和对象特性选取，通常阶跃输入扰动模型的Q/R值介于周期性扰动与斜坡扰动之间，斜坡扰动的Q/R值一般最小。

有益效果：本发明提供的火电机组过热汽温多模型扰动估计预测控制方法，与现有技术相比具有以下优点：具有多模型建立方便、建模工作量小、控制算法计算量小；能够延续预测控制，便于处理大惯性、纯滞后对象，在此基础上进行的过程扰动建模，获得了比一般预测控制算法更优的抗干扰能力，能够在过程不可测干扰作用情况下，快速估计并抑制扰动对系统的影响；应用于电站过热汽温系统能够有效抑制干扰，确保汽温稳定在设定值附近；同时该方法具有一般预测控制器同等的设定值跟踪调节能力，总体提高了过热汽温的调节品质。

附图说明

图1为本发明的系统结构框图；

图2为本发明的对象多模型扰动估计预测控制算法框图；

图3为本发明与普通测测控制算法在无干扰情况下，汽温设定值变化时的控制效果对比图；

图4为本发明与普通测测控制算法在存在负荷侧斜坡输入不可测扰动时的控制效果对比图；

图5为本发明与普通测测控制算法在过热汽温受负荷侧周期性输入不可测扰动时的控制效果对比图；

图6为本发明与普通测测控制算法在过热汽温同时受负荷侧和燃料侧两种类型不可测扰动时的控制效果对比图。

具体实施方式

下面结合附图对本发明作更进一步的说明。

一种火电机组过热汽温多模型扰动估计预测控制方法，通过设计多个包含不同扰动模型的扩增系统模型，采用卡尔曼滤波方法估计扩增系统状态，通过贝叶斯概率加权方法计算平均预测输出值，每一采样时刻由预测控制优化算法计算最佳减温喷水量，使系统能够及时有效抑制不可测干扰，同时算法在无扰情况下，与普通预测控制器具有同等的设定值跟踪和调节能力，本发明总体提高了过热汽温的调节品质。该方法包括如下步骤：

一种火电机组过热汽温多模型扰动估计预测控制方法，包括如下步骤：

(1)在稳态工况下，将过热汽温系统切换到手动状态，以减温喷水流量为输入进行过热汽温的开环阶跃响应试验，利用两点法求出过热汽温对减温喷水流量的传递函数W(s)(℃/(kg/s))；

(2)根据线性系统理论的系统实现方法，将过热汽温对减温喷水流量的传递函数W(s)(℃/(kg/s))转换为系统状态空间模型 $\left\{\begin{array}{c}\stackrel{\&CenterDot;}{x}=\mathrm{Ax}+\mathrm{Bu}\\ y=\mathrm{Cx}+\mathrm{Du}\end{array}\right.,$ 其中[A,B,C,D]为系统矩阵，x为过热汽温系统的状态，y为过热汽温系统的输出，u为预测控制器的输出；

(3)建立n种类型扰动模型分别扩增过热汽温系统状态，得到n个扩增系统状态模型[ⁱA^a，ⁱB^a，ⁱC^a，ⁱD^a]，其中ⁱD^a＝0，i＝1,2,…,n，所述扩增系统状态空间模型的结构如式(1)和式(2)所示：

其中，表示k时刻扩增的状态向量，包括两个子向量ⁱx_k和 ⁱx_k表示k时刻过热汽温系统原有的状态向量，表示k时刻第i个扰动模型的扰动、扰动变化率或周期的状态；ω_k-1表示k-1时刻过程噪声，v_k表示k时刻输出噪声；

采用式(3)和式(4)简洁表达式(1)和式(2)为：

$x_k^a{}_i=A_a^{i}x_{k-1}^a{}_i+B_a^i{u}_{k-1}+\mathrm{\&Omega;}_a^i{\mathrm{\&omega;}}_{k-1}---\left(3\right)$

$y_k^i=C_a^{i}x_k^a{}_i+v_k---\left(4\right)$

其中代表 $\left[\begin{array}{c}x_k^i\\ d_k^a{}_i\end{array}\right],$ ⁱA^a代表 $\left[\begin{array}{cc}A_1^i& A_2^i\\ A_3^i& A_4^i\end{array}\right],$ ⁱB^a代表 $\left[\begin{array}{c}B_1^i\\ B_2^i\end{array}\right],$ ⁱC^a代表[ⁱC₁ ⁱC₂]；

(4)采用卡尔曼滤波方法对n个扩增系统状态模型分别进行状态估计，如式(5)、式(6)和式(7)：

$\hat{x}_{k|k-1}^a{}_i=A_a^i\hat{x}_{k-1|k+1}^a{}_i+B_a^i{u}_{k-1}---\left(5\right)$

$\hat{x}_{k|k}^a{}_i=\hat{x}_{k|k-1}^a{}_i+L_k^i(y_k-C_a^i\hat{x}_{k|k-1}^a{}_i)---\left(6\right)$

$\hat{y}_{k|k-1}^i=C_a^i\hat{x}_{k|k-1}^a{}_i---\left(7\right)$

在状态未知情况下，首先使用式(5)和式(7)预估状态和输出再通过式(6)由实际输出y_k和预估输出间的误差来矫正过热汽温系统状态得到新的预估状态和输出状态估计增益矩阵ⁱL_k通过求解Riccati方程获得，Riccati方程中的Qi代表输入噪声项的协方差，Ri代表输出噪声项的协方差，卡尔曼滤波参数Qi/Ri为预测控制器的一个调节参数；

(5)采用贝叶斯概率加权方法对n个扩增系统状态模型进行加权平均，获得过热汽温的P步平均预测值，应用预测控制算法计算最佳减温喷水流量，将计算得到的最佳减温喷水流量作为过热汽温系统的控制输入量，具体包括如下步骤：

(51)根据预估输出通过式(8)计算实际输出y_k与预估输出间的偏差ⁱε_k，并根据偏差ⁱε_k计算各个扩增系统状态模型的权重ⁱw_k：

$\mathrm{\&epsiv;}_k^i=y_k-\hat{y}_{k|k-1}^i---\left(8\right)$

(52)通过递推贝叶斯公式计算各扰动模型与过热汽温系统的匹配概率，通过上一时刻的匹配概率ⁱρ_k-1递推获得当前时刻的匹配概率ⁱρ_k，如式(9)：

$\mathrm{\&rho;}_k^i=\frac{\exp (-\frac{1}{2}\mathrm{\&epsiv;}_k^T{}_i\mathrm{\&Lambda;}_i\mathrm{\&epsiv;}_k^i)\mathrm{\&rho;}_{k-1}^i}{{\mathrm{\&Sigma;}}_{j=1}^4\exp (-\frac{1}{2}\mathrm{\&epsiv;}_k^T{}_j\mathrm{\&Lambda;}_i\mathrm{\&epsiv;}_k^j)\mathrm{\&rho;}_{k-1}^j}---\left(9\right)$

设定下限值δ，当上一时刻的匹配概率ⁱρ_k-1低于下限值δ时，以下限值δ置换ⁱρ_k-1后带入式(9)进行ⁱρ_k计算；贝叶斯概率加权矩阵ⁱΛ为预测控制器的一个调节参数；

(53)对各扰动模型与热汽温系统的匹配概率进行规格化处理，如式(10)：

$w_k^i=\left\{\begin{array}{cc}\frac{\mathrm{\&rho;}_k^i}{{\mathrm{\&Sigma;}}_{j=1}^4\mathrm{\&rho;}_k^j}& \mathrm{\&rho;}_k^i>\mathrm{\&delta;}\\ 0& \mathrm{\&rho;}_k^i<\mathrm{\&delta;}\end{array}\right.---\left(10\right)$

(55)根据式(11)计算过热汽温的P步平均预测值：

${\stackrel{\&OverBar;}{y}}_{k+j|k}={\mathrm{\&Sigma;}}_{i=1}^{4}w_k^i\hat{y}_{k+j|k}^i---\left(11\right)$

其中表示第i个扩增系统状态模型的预测输出值，ⁱw_k表示第i个扩增系统状态模型的权值，为所有扩增系统状态模型的平均预测输出值，j＝1,2,…,P；

(56)计算减温喷水流量控制量增量Δu_k，取式(12)性能指标函数式：

其中，Y_sp为过热汽温设定值，为过热汽温的P步平均预测矩阵， $\stackrel{\&OverBar;}{Y}={[{\stackrel{\&OverBar;}{y}}_{k+1|k},{\stackrel{\&OverBar;}{y}}_{k+2|k},\&CenterDot;\&CenterDot;\&CenterDot;{,\stackrel{\&OverBar;}{y}}_{k+P|k}]}^T,\mathrm{\&Delta;U}={[{\mathrm{\&Delta;u}}_k,{\mathrm{\&Delta;u}}_{k+1},\&CenterDot;\&CenterDot;\&CenterDot;,{\mathrm{\&Delta;u}}_{k+m-1}]}^T,$ W_y为误差权矩阵和W_u为控制量权矩阵；

每一采样时刻进行一次式(12)的优化计算，获得m时域的最佳系统控制输入增量ΔU_opt＝[Δu_k，Δu_k+1，…，Δu_k+m-1]^T；只以当前时刻的控制输入增量Δu_k对当前过热汽温系统实施控制，计算减温喷水流量控制量u_k，如式(13)：

u_k＝u_k-1+Δu_k   (13)

(6)设置预测控制器的调节参数，包括采样时间Ts、预测步数P、减温喷水流量控制步数M、输出误差权矩阵Wy、控制权矩阵Wu、n个扩增系统状态模型对应的卡尔曼滤波参数Qi/Ri和贝叶斯概率加权矩阵Λi，i＝1,2,…,n；Ts可以用经验规则T95/Ts＝5～15来选取，其中，T95为过渡过程上升到95％的调节时间；P选为近似等于过热汽温阶跃响应的上升时间；Qi/Ri一般根据扰动模型形式和对象特性选取，通常阶跃输入扰动模型的Q/R值介于周期性扰动与斜坡扰动之间，斜坡扰动的Q/R值一般最小；

(7)向过热汽温系统输入减温喷水流量控制量u_k，获得热汽温系统实际输出y_k；

(8)重复步骤(4)～(7)实现连续控制。

由于两级减温喷水的控制都采用相同的控制策略，两者的设计过程基本相同，因此下面重点以二级减温控制为例进行说明。具体步骤为：

(1)在稳态工况下，获得以二级减温喷水量为输入，二级过热汽温为输出的单入单出系统状态空间模型 $\left\{\begin{array}{c}{x}_k={\mathrm{Ax}}_{k-1}+{\mathrm{Bu}}_{k-1}+{\mathrm{\&omega;}}_{k-1}\\ y_k={\mathrm{Cx}}_k+v_k\end{array}\right..$ 其系数矩阵[A,B,C]分别：A＝[0.9858,0.0001；0.0002,0.9994]；B＝[0.0009；-0.0020]；C＝[0,1]；

(2)选取工业过程控制中四类常见扰动，分别建立扰动模型估计实际过程中的不可测扰动。四种扰动模型的表达式如下：

扰动模型一：输出扰动模型，如式(14-16)

x_k＝Ax_k-1+Bu_k-1   (14)

d_k＝d_k-1+ω_k-1   (15)

y_k＝Cx_k+d_k+v_k   (16)

扰动模型二：阶跃输入扰动模型，如式(17-19)

x_k＝Ax_k-1+Bu_k-1+B^dd_k-1   (17)

d_k＝d_k-1+ω_k-1   (18)

y_k＝Cx_k+v_k   (19)

扰动模型三：斜坡输入扰动模型，如式(20-23)

x_k＝Ax_k-1+Bu_k-1+B^dd_k-1   (20)

d_k＝d_k-1+Δd_k-1   (21)

Δd_k＝Δd_k-1+ω_k-1(22)

y_k＝Cx_k+v_k   (23)

说明：d_k为待估计的扰动量，Δd_k为扰动变化速率的估计值，表示斜坡扰动是以固定增速率进入系统的。

扰动模型四：周期性输入扰动模型，周期性扰动通常以连续函数模型表示，如式(24)

$\stackrel{\&CenterDot;\&CenterDot;}{d}+a^{2}d=d^{*}---\left(24\right)$

其中a表示周期性扰动的频率，d^*表示扰动源。上述微分方程表示成状态空间方程的形式如式(25)

在取定采样时间Ts后，将上述模型离散化得到式(26)

$\left[\begin{array}{c}{d}_k\\ {\mathrm{\&Delta;d}}_k\end{array}\right]=d_k^a=A^d{d}_{k-1}^a+B^{*}{d}^{*}---\left(26\right)$

扰动源d^*可根据具体应用过程进行设置。

综上，四种扰动模型可以由如式(1-2)相同的模型结构表示

扩增的状态向量由两个子向量组成：ⁱx_k系统模型原有的状态向量，表示扰动模型簇中第i类模型的扰动、扰动变化率或周期的状态。

上述模型可以写成如式(3-4)所示简洁的表达形式

$x_k^a{}_i=A_a^{i}x_{k-1}^a{}_i+B_a^i{u}_{k-1}+\mathrm{\&Omega;}_a^i{\mathrm{\&omega;}}_{k-1}---\left(3\right)$

$y_k^i=C_a^{i}x_k^a{}_i+v_k---\left(4\right)$

其中代表 $\left[\begin{array}{c}x_k^i\\ d_k^a{}_i\end{array}\right],$ ⁱA^a代表 $\left[\begin{array}{cc}A_1^i& A_2^i\\ A_3^i& A_4^i\end{array}\right],$ ⁱB^a代表 $\left[\begin{array}{c}B_1^i\\ B_2^i\end{array}\right],$ ⁱC^a代表[ⁱC₁ ⁱC₂]。

故扩增系统参数矩阵分别为：

输出扰动扩增模型系数矩阵为

¹A^a＝[0.9858,0.0001,0；0.0002,0.9994,0；0,0,1.0000]；

¹B^a＝[0.0009；-0.0020；0]；

¹C^a＝[0,1,1]；

阶跃输入扰动扩增模型系数矩阵为

²A^a＝[0.9858,0.0001,0.0009；0.0002,0.9994,-0.0020；0,0,1.0000]；

²B^a＝[0.0009；-0.0020；0]；

²C^a＝[0,1,0]；

斜坡输入扰动扩增模型系数矩阵为

³A^a＝[0.9858,0.0001,0.0009,0；0.0002,0.9994,-0.0020,0；0,0,1.0000,1.0000；0,0,0,1.0000]；

³B^a＝[0.0009；-0.0020；0；0]；

³C^a＝[0,1,0,0]；

周期性输入扰动扩增模型系数矩阵为

^aA^a＝[0.9858,0.0001,0.0009,0,0；0.0002,0.9994,-0.0020,0,0；0,0,0.9921,0.1253,0.0079；0,0,-0.1253,0.9921,0.1253；0,0,0,0,1.0000]；

⁴B^a＝[0.0009；-0.0020；0；0；0]；

⁴C^a＝[0,1,0,0,0]；

(3)采用卡尔曼滤波方法对n个扩增模型分别进行状态估计，如式(5-7)：

$\hat{x}_{k|k-1}^a{}_i=A_a^i\hat{x}_{k-1|k+1}^a{}_i+B_a^i{u}_{k-1}---\left(5\right)$

$\hat{x}_{k|k}^a{}_i=\hat{x}_{k|k-1}^a{}_i+L_k^i(y_k-C_a^i\hat{x}_{k|k-1}^a{}_i)---\left(6\right)$

$\hat{y}_{k|k-1}^i=C_a^i\hat{x}_{k|k-1}^a{}_i---\left(7\right)$

在未知状态测量值的的情况下首先预估扩增系统状态和输出获得系统实际输出y_k后，再由实际输出和输出估计值间的误差矫正系统状态得到新的状态估计值和输出估计值四个扩增系统状态估计的关键在于获得观测器增益矩阵ⁱL_k。对于输出扰动模型，ⁱL_k等价于非周期观测器，取为L1＝[0，0，1]。其余3个模型的ⁱL_k均可通过求解Riccati方程获得。Riccati方程中的Qi和Ri分别代表输入输出噪声项的协方差。这里取阶跃、斜坡、周期性输入扰动模型的Q/R值分别为100，2，20。

(4)对4个模型采用贝叶斯概率加权方法进行加权平均，获得平均预测输出值，应用预测控制算法计算最佳减温喷水量，其算法过程如下：

得到各子模型的输出估计值如式(8-10)分别计算实际输出值y_k与估计值间的偏差ⁱε_k，以此为基础计算各子模型的权重ⁱw_k。

$\mathrm{\&epsiv;}_k^i=y_k-\hat{y}_{k|k-1}^i---\left(8\right)$

扰动模型簇中各子模型与实际对象的匹配概率可通过递推贝叶斯公式进行计算。由上一时刻的模型匹配概率ⁱρ_k-1递推获得当前时刻的概率ⁱρ_k，如式(9)。

$\mathrm{\&rho;}_k^i=\frac{\exp (-\frac{1}{2}\mathrm{\&epsiv;}_k^T{}_i\mathrm{\&Lambda;}_i\mathrm{\&epsiv;}_k^i)\mathrm{\&rho;}_{k-1}^i}{{\mathrm{\&Sigma;}}_{j=1}^4\exp (-\frac{1}{2}\mathrm{\&epsiv;}_k^T{}_j\mathrm{\&Lambda;}_i\mathrm{\&epsiv;}_k^j)\mathrm{\&rho;}_{k-1}^j}---\left(9\right)$

实际应用中协方差矩阵ⁱΛ未知，因此将ⁱΛ作为调节参数。由于上述概率计算是递推形式，当前时刻的值由上一时刻获得。若某一时刻概率计算值为0，则之后时刻的概率计算无效。因此设定一个下限δ，所有低于该限值的概率均被置为δ。对扰动模型簇中的概率进行规格化处理，便于后续预测控制算法的实现。规格化公式如式(10)

$w_k^i=\left\{\begin{array}{cc}\frac{\mathrm{\&rho;}_k^i}{{\mathrm{\&Sigma;}}_{j=1}^4\mathrm{\&rho;}_k^j}& \mathrm{\&rho;}_k^i>\mathrm{\&delta;}\\ 0& \mathrm{\&rho;}_k^i<\mathrm{\&delta;}\end{array}\right.---\left(10\right)$

(5)设置控制器相关参数，取采样时间T_s为0.5，预测步数P为50，减温水调门控制步数M为3，输出误差权矩阵Wy＝I₅₀，控制权矩阵Wu＝I₃。

(6)系统的P步平均预测输出值由式(11)线性平均模型获得

${\stackrel{\&OverBar;}{y}}_{k+j|k}={\mathrm{\&Sigma;}}_{i=1}^{4}w_k^i\hat{y}_{k+j|k}^i---\left(11\right)$

其中和ⁱw_k分别表示扰动模型簇中第i个子模型的预测输出值和权值，为扰动模型簇的平均预测输出值。

(7)计算减温喷水阀门控制量增量Δu_k，取式(12)性能指标函数式：

每一采样时刻进行一次优化计算，获得m时域的最优控制输出增量[Δu_k，Δu_k+1，…，Δu_k+m-1]^T。只有当前时刻的控制输入增量Δu_k对系统实施控制。计算减温水调门控制量，如式(13)。

u_k＝u_k-1+Δu_k   (13)

(8)输出减温水调门控制量u_k，获得汽温控制系统实际汽温输出值y_k，重复(3)到(8)步。

如图3所示，过热汽温设定值阶跃增加5℃情况下，MDMPC代表本发明对过热汽温系统的优化控制效果曲线，MPC代表普通预测控制器的控制效果曲线。从图3可以看出，两条控制曲线基本重合，说明本发明具有普通预测控制器同等的设定值跟踪和调节能力。如图4所示，当负荷侧有斜坡扰动进入过热汽温系统时，本发明能够快速识别扰动类型为斜坡型，提前动作减温喷水，得到比普通预测控制器更小的过热汽温动态和稳态偏差。如图5所示，当负荷侧有周期性扰动进入过热汽温系统时，本发明能够快速识别扰动类型为周期型。当与普通预测控制器获得基本一致的控制效果时，本发明能够预测扰动提前动作减温喷水，在调节过程中的控制量波动小，能够有效延长阀门设备使用寿命。如图6所示，当存在负荷侧周期性扰动和燃料侧斜坡扰动同时进入过热汽温系统时，本发明能够快速识别两种扰动类型，得到比普通预测控制器更小的过热汽温动态和稳态偏差，提高了过热汽温的调节品质。

以上所述仅是本发明的优选实施方式，应当指出：对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可以做出若干改进和润饰，这些改进和润饰也应视为本发明的保护范围。

Claims (1)

1.一种火电机组过热汽温多模型扰动估计预测控制方法，其特征在于：包括如下步骤：

(1)在稳态工况下，将过热汽温系统切换到手动状态，以减温喷水流量为输入进行过热汽温的开环阶跃响应试验，利用两点法求出过热汽温对减温喷水流量的传递函数W(s)(℃/(kg/s))；

(2)根据线性系统理论的系统实现方法，将过热汽温对减温喷水流量的传递函数W(s)(℃/(kg/s))转换为系统状态空间模型 $\left\{\begin{array}{c}\stackrel{\&CenterDot;}{x}=\mathrm{Ax}+\mathrm{Bu}\\ y=\mathrm{Cx}+\mathrm{Du}\end{array}\right.,$ 其中[A,B,C,D]为系统矩阵，x为过热汽温系统的状态，y为过热汽温系统的输出，u为预测控制器的输出；

(3)建立n种类型扰动模型分别扩增过热汽温系统状态，得到n个扩增系统状态模型[ⁱA^a，ⁱB^a，ⁱC^a，ⁱD^a]，其中ⁱD^a＝0，i＝1,2,…,n，所述扩增系统状态空间模型的结构如式(1)和式(2)所示：

其中，表示k时刻扩增的状态向量，包括两个子向量ⁱx_k和 ⁱx_k表示k时刻过热汽温系统原有的状态向量，表示k时刻第i个扰动模型的扰动、扰动变化率或周期的状态；ω_k-1表示k-1时刻过程噪声，v_k表示k时刻输出噪声；

采用式(3)和式(4)简洁表达式(1)和式(2)为：

$x_k^a{}_i=A_a^{i}x_{k-1}^a{}_i+B_a^i{u}_{k-1}+\mathrm{\&Omega;}_a^i{\mathrm{\&omega;}}_{k-1}---\left(3\right)$

$y_k^i=C_a^{i}x_k^a{}_i+v_k---\left(4\right)$

其中代表 $\left[\begin{array}{c}x_k^i\\ d_k^a{}_i\end{array}\right],$ ⁱA^a代表 $\left[\begin{array}{cc}A_1^i& A_2^i\\ A_3^i& A_4^i\end{array}\right],$ ⁱB^a代表 $\left[\begin{array}{c}B_1^i\\ B_2^i\end{array}\right],$ ⁱC^a代表[ⁱC₁ ⁱC₂]；

(4)采用卡尔曼滤波方法对n个扩增系统状态模型分别进行状态估计，如式(5)、式(6)和式(7)：

$\hat{x}_{k|k-1}^a{}_i=A_a^i\hat{x}_{k-1|k+1}^a{}_i+B_a^i{u}_{k-1}---\left(5\right)$

$\hat{x}_{k|k}^a{}_i=\hat{x}_{k|k-1}^a{}_i+L_k^i(y_k-C_a^i\hat{x}_{k|k-1}^a{}_i)---\left(6\right)$

$\hat{y}_{k|k-1}^i=C_a^i\hat{x}_{k|k-1}^a{}_i---\left(7\right)$

在状态未知情况下，首先使用式(5)和式(7)预估状态和输出再通过式(6)由实际输出y_k和预估输出间的误差来矫正过热汽温系统状态得到新的预估状态和输出状态估计增益矩阵ⁱL_k通过求解Riccati方程获得，Riccati方程中的Qi代表输入噪声项的协方差，Ri代表输出噪声项的协方差，卡尔曼滤波参数Qi/Ri为预测控制器的一个调节参数；

(5)采用贝叶斯概率加权方法对n个扩增系统状态模型进行加权平均，获得过热汽温的P步平均预测值，应用预测控制算法计算最佳减温喷水流量，将计算得到的最佳减温喷水流量作为过热汽温系统的控制输入量，具体包括如下步骤：

(51)根据预估输出，通过式(8)计算实际输出y_k与预估输出间的偏差ⁱε_k，并根据偏差ⁱε_k计算各个扩增系统状态模型的权重ⁱw_k：

$\mathrm{\&epsiv;}_k^i=y_k-\hat{y}_{k|k-1}^i---\left(8\right)$

(52)通过递推贝叶斯公式计算各扰动模型与过热汽温系统的匹配概率，通过上一时刻的匹配概率ⁱρ_k-1递推获得当前时刻的匹配概率ⁱρ_k，如式(9)：

$\mathrm{\&rho;}_k^i=\frac{\exp (-\frac{1}{2}\mathrm{\&epsiv;}_k^T{}_i\mathrm{\&Lambda;}_i\mathrm{\&epsiv;}_k^i)\mathrm{\&rho;}_{k-1}^i}{{\mathrm{\&Sigma;}}_{j=1}^4\exp (-\frac{1}{2}\mathrm{\&epsiv;}_k^T{}_j\mathrm{\&Lambda;}_i\mathrm{\&epsiv;}_k^j)\mathrm{\&rho;}_{k-1}^j}---\left(9\right)$

设定下限值δ，当上一时刻的匹配概率ⁱρ_k-1低于下限值δ时，以下限值δ置换ⁱρ_k-1后带入式(9)进行ⁱρ_k计算；贝叶斯概率加权矩阵ⁱΛ为预测控制器的一个调节参数；

(53)对各扰动模型与热汽温系统的匹配概率进行规格化处理，如式(10)：

$w_k^i=\left\{\begin{array}{cc}\frac{\mathrm{\&rho;}_k^i}{{\mathrm{\&Sigma;}}_{j=1}^4\mathrm{\&rho;}_k^j}& \mathrm{\&rho;}_k^i>\mathrm{\&delta;}\\ 0& \mathrm{\&rho;}_k^i<\mathrm{\&delta;}\end{array}\right.---\left(10\right)$

(55)根据式(11)计算过热汽温的P步平均预测值：

${\stackrel{\&OverBar;}{y}}_{k+j|k}={\mathrm{\&Sigma;}}_{i=1}^{4}w_k^i\hat{y}_{k+j|k}^i---\left(11\right)$

其中表示第i个扩增系统状态模型的预测输出值，ⁱw_k表示第i个扩增系统状态模型的权值，为所有扩增系统状态模型的平均预测输出值，j＝1,2,…,P；

(56)计算减温喷水流量控制量增量Δu_k，取式(12)性能指标函数式：

其中，Y_sp为过热汽温设定值，为过热汽温的P步平均预测矩阵， $\stackrel{\&OverBar;}{Y}={[{\stackrel{\&OverBar;}{y}}_{k+1|k},{\stackrel{\&OverBar;}{y}}_{k+2|k},\&CenterDot;\&CenterDot;\&CenterDot;{,\stackrel{\&OverBar;}{y}}_{k+P|k}]}^T,\mathrm{\&Delta;U}={[{\mathrm{\&Delta;u}}_k,{\mathrm{\&Delta;u}}_{k+1},\&CenterDot;\&CenterDot;\&CenterDot;,{\mathrm{\&Delta;u}}_{k+m-1}]}^T,$ W_y为误差权矩阵和W_u为控制量权矩阵；

每一采样时刻进行一次式(12)的优化计算，获得m时域的最佳系统控制输入增量ΔU_opt＝[Δu_k，Δu_k+1，…，Δu_h+m-1]^T；只以当前时刻的控制输入增量Δu_k对当前过热汽温系统实施控制，计算减温喷水流量控制量u_k，如式(13)：

u_k＝u_k-1+Δu_k   (13)

(6)设置预测控制器的调节参数，包括采样时间Ts、预测步数P、减温喷水流量控制步数M、输出误差权矩阵Wy、控制权矩阵Wu、n个扩增系统状态模型对应的卡尔曼滤波参数Qi/Ri和贝叶斯概率加权矩阵Λi，i＝1,2,…,n；

(7)向过热汽温系统输入减温喷水流量控制量u_k，获得热汽温系统实际输出y_k。