CN1758026A - 基于新型广义预测控制的在线关键参数测量方法 - Google Patents

Abstract

本发明公开了一种基于广义预测控制的在线关键参数测量方法。关键参量的实时测量是过程控制中至关重要的问题，使用广义预测控制能够使得关键参数测量具有良好的精度和趋势。本发明提出的实时在线测量模型很好地利用了广义预测控制中的思想，实现一些时延长或难以用普通传感器测量的关键参量的实时测量，为在线关键参数测量方法提供理论基础。实现了工业过程关键参量实时在线测量，解决了锅炉控制中测量与控制的时域适配问题，提出的基于广义预测的含氧量软测量技术也为锅炉燃烧过程中其它关键参量如烟气温度的准确、实时、在线软测量开辟了新的途径，对实现燃烧系统的闭环控制和优化运行具有重要的意义。本发明同样适用于毛坯加热炉、煤气生产炉等生产过程中烟气含氧量、CO烟气含量的测量。

Description

基于新型广义预测控制的在线关键参数测量方法

技术领域

本发明涉及流程工业测控领域，是一种基于新型广义预测控制的在线关键参数测量方法。

背景技术

在流程工业领域，特别是在热工行业中，存在大量的难以测量的关键参数，一般只能靠人工手段或者是人工经验判断的情况，实现在线关键参数测量是必须的。故目前主要采用软测量模型进行实时的测量，即通过测量易于检测的辅助变量，基于估计模型来估算出难以甚至无法测量的变量。

目前热工行业中测量设备精度不高，投资大，使用寿命短，而且测量滞后较大，很难实现燃烧过程的在线实时测控。近年来逐步采用了软测量技术来解决此类问题。采用间接测量的思路，利用易于获取的其它参数，通过计算来实现被检测量的估计，是近来在过程控制和检测领域涌现出的一种新技术—软测量技术(Soft-Sensing Technique)。软测量技术也称为软仪表技术，就是利用易测过程变量(如工业过程中容易获取的压力、温度等过程参数)，依据这些易测过程变量与难以直接测量的待测过程变量(如锅炉燃烧过程中的烟气含氧量等)之间的软测量模型，通过各种计算和估计方法，从而实现对待测过程变量的测量。

关键参量是指那些在控制中必不可少的且用常规测量方法难以得到实时准确数据的参数，如炉膛温度、烟气含氧量(O₂)、碳氧化合物含量(CO₂、CO)、氮氧化合物含量(NOx)、水蒸气含量(H₂O)、粉尘浓度以及燃烧过程时滞常数(τ)等。其中，炉膛温度是保证燃烧和锅炉运行的关键参量；含氧量是关系到风煤配比、评价燃烧优劣重要过程参数，也是重要的目标参数；而SO₂、CO₂、CO、H₂O、NOx、粉尘浓度等则是热损计算和评估锅炉效率必不可少的参数。这些参数采用常规的检测方法很难保证实时、准确的要求。而对于那些虽然是控制必需但可容易获取的参数，如炉膛负压、出水压力、烟气温度等参数，则不列入关键参量中。

目前，工程上烟气含氧量主要依赖氧化锆传感器实现测量。SO₂、CO₂等参数由安装在烟囱上烟气分析仪提供。由于反应过程的时滞、安装位置限制、仪表自身特性等原因，使得这些参数不能满足测量与控制时域适配的要求。目前，之所以一些高水平的理论研究不能应用到工程实际，其中一个重要的原因就是忽略了这一因素的影响。

氧化锆测量烟气含氧量的弊病主要有，一是中毒现象时有发生，造成测量失准；二是由于燃烧反应的过渡过程、烟气流通的时延造成测量滞后；三是氧化锆自身时滞性造成测量滞后。事实上，现在工程上的氧化锆测量氧含量只有在稳态情况下才可能反应工况的真正结果，类似的，SO₂的测量也存在同样问题。

关键参量的实时测量是锅炉燃烧过程控制中至关重要的问题。对于解决那些具有大时滞特性信号，或不能用直接方法测量的关键参量，通常采用软测量的方法。

广义预测控制是具有预测模型、滚动优化和反馈校正和自适应特点的一种先进算法，论文研究的新型算法更具有快速收敛性的特点，可以把广义预测技术应用到烟气含氧量测量之中，实现锅炉燃烧过程关键参量的软测量。

发明内容

本发明的目的是提供一种基于新型广义预测控制的在线关键参数测量方法，它利用广义预测控制的优势，对软测量关键参数模型进行滚动优化和反馈校正，使得软测量关键参数模型输出具有良好的精度和趋势。

本发明解决其技术问题所采用的技术方案是：根据过程控制的时滞特性和关键参数的测量方法，利用新型广义预测控制器，建立了基于广义预测算法的软测量预测模型，实现了过程控制关键参量的实时在线测量，具体包括测量模型建立、初始化、新型预测控制方法、反馈校正和基于新型广义预测控制的在线关键参数测量方法四个步骤：

(1)建立测量模型：根据广义预测控制技术特征和过程控制时滞特性建立测量模型：

(2)初始化及求解最优控制律：根据测量模型设定控制参数，利用新型广义预测控制方法求解最优控制律；

(3)模型反馈校正：采用时变遗忘因子两段参数估计方法进行模型反馈校正，优化控制参数；

(4)基于广义预测控制在线关键参数测量方法：根据新型广义预测控制器计算、输出关键参数测量值。

所述建立的测量模型，具体如下：

确定系统的过程控制的目标函数为

$J=E\{{(\mathrm{\ψ}-w)}^T(\mathrm{\ψ}-w)+\mathrm{\λ}{\mathrm{\ψ}}^T\mathrm{\ψ}\}$

$=E\{{(G\stackrel{\→}{u}+\stackrel{\→}{f}+\stackrel{\→}{e}-\stackrel{\→}{w})}^T(G\stackrel{\→}{u}+\stackrel{\→}{f}+\stackrel{\→}{e}-\stackrel{\→}{w})+\mathrm{\λ}{(P\stackrel{\→}{u}+\stackrel{\→}{\mathrm{\δ}})}^T(P\stackrel{\→}{u}+\stackrel{\→}{\mathrm{\δ}})\}$

为了求得最优控制律 由目标函数式可得

$(G^{T}G+\mathrm{\λ}{P}^{T}P)\stackrel{\→}{u}=G^T(\stackrel{\→}{w}-\stackrel{\→}{f}-\stackrel{\→}{e})-\mathrm{\λ}{P}^T\stackrel{\→}{\mathrm{\δ}}$

由数学分析可知由方程所求得的 即为使得目标函数式中的J取得最小值的最优解。

所述初始化及求解最优控制律，具体如下：

由于新型控制算法采用矩阵递推迭代方法，避免了传统方法求解最优控制律中逆矩阵的计算，运算复杂度下降了M(矩阵阶数)的平方阶次，节省的计算量和存储量这样就可以扩大预测时域和控制时域的选取范围，增加测控系统的稳定性，可以实现锅炉关键参数实时在线软测量；

递推优化算法如下：

记 $G^{T}G+\mathrm{\λ}{P}^{T}P=\left[\begin{array}{cccc}{a}_{11}& a_{12}& \·\·\·& a_{1M}\\ a_{21}& a_{22}& \·\·\·& a_{2M}\\ \·\·\·& \·\·\·& \·\·\·& \·\·\·\\ a_{M1}& a_{M2}& \·\·\·& a_{\mathrm{MM}}\end{array}\right]=\left[\begin{array}{c}{{\stackrel{\→}{a}}_1}^t\\ {{\stackrel{\→}{a}}_2}^t\\ \·\\ \·\\ \·\\ {{\stackrel{\→}{a}}_M}^t\end{array}\right]$

$G^T(\stackrel{\→}{w}-\stackrel{\→}{f})-\mathrm{\λ}{P}^T\stackrel{\→}{\mathrm{\δ}}={\stackrel{\→}{b}}^t=[{\stackrel{\→}{b}}_1,\·\·\·,{\stackrel{\→}{b}}_M]$

(1)读入M，t， ${{\stackrel{\→}{a}}_i}^t=(a_{\mathrm{il}},\·\·\·,a_{\mathrm{iM}}),$ i＝1，2，…，M

e₁＝(1，0，…，0)，…，e_j＝(0，…，0，1，0，…，0)，…，e_M＝(0，0，…，1)， ${\stackrel{\→}{b}}^t=[{\stackrel{\→}{b}}_1,\·\·\·,{\stackrel{\→}{b}}_M],$

${\stackrel{\→}{u}}^T=(u_1,\·\·\·,u_M),$ u_i＝Δu_f(t+i-1)，i＝1，2，…，M

(2)任取 ${\stackrel{\→}{u}}^1\∈E^M,$ 这里不妨取 ${\stackrel{\→}{u}}^1=(0,\·\·\·,0)\∈E^M,$ 取H₁∈E^M×M为任意非奇异对称阵，这里不妨取 $H_1^T=(e_1,\·\·\·,e_M),$ 置i＝1，iflag＝0；

(3)计算在

处残差向量的第i个分量， ${\mathrm{\τ}}_i={\stackrel{\→}{a}}_i^T{\stackrel{\→}{u}}^i-{\stackrel{\→}{b}}_i,$ 计算搜索向量 $S^i=H_i{\stackrel{\→}{a}}_i;$

(4)若 $S^i\≠\stackrel{\→}{0},$ 转步(5)；若 $S^i=\stackrel{\→}{0},$ 同时 ${\mathrm{\τ}}_i=\stackrel{\→}{0},$ 则置 ${\stackrel{\→}{u}}^{i+1}={\stackrel{\→}{u}}^i,$ H_i+1＝H_i，iflag＝iflag+1，若i＜M转步(7)；否则停止计算；此时

为方程式的解。若 $S^i=\stackrel{\→}{0},$ 但 ${\mathrm{\τ}}^i\≠\stackrel{\→}{0},$ 则置iflag＝-i停止计算，此时方程组式不相容，也即目标函数式无最优解；

(5)修正解的近似值， ${\stackrel{\→}{u}}^{i+1}={\stackrel{\→}{u}}^i-{\mathrm{\λ}}_i{S}^i,$ 其中步长 ${\mathrm{\λ}}_i={\mathrm{\τ}}_i/{\stackrel{\→}{a}}_i^T{S}^i,$

若i＝M，停止计算，

为方程式的解；

(6)修正矩阵H_i， $H_{i+1}=H_i-S^i{\left(S^i\right)}^T/{\stackrel{\→}{a}}_i^T{S}^i;$

(7)置i＝i+1，转步(3)；

$\stackrel{\→}{u}=u^{M+1}$ 停止，此时得到方程式的解，也即得到了使得目标函数式中J取得最小值的最优解。从而

$\mathrm{\Δ}{u}_f\left(t\right)=[\mathrm{1,0},\·\·\·,0]\stackrel{\→}{u}$

则在t时刻的控制量为

u(t)＝T(Z^-1)[u_f(t-1)+Δu_f(t)]

所述初始化及求解最优控制律，模型的反馈校正，是利用此采样时刻的广义预测模型计算最优值和给定值之间的偏差实时在线修正模型参数，这种基于广义预测控制模型，根据系统的实际输出不断对预测控制输出值进行滚动优化修正，而且利用了反馈信息，构成闭环优化校正；反馈校正参数估计采用时变遗忘因子两段参数估计性能指标加权的广义预测自适应的优化控制算法，具体如下：

预测模型采用

A(z^-1)Δy_f(t)＝B(z^-1)Δu_f(t^-1)+D(z^-1)Δv_f(t-1)+C(z^-1)ξ_f(t)

定理若A(z^-1)和C(z^-1)是稳定的，则预测模型式可用如下高阶CAR模型近似代替

${\mathrm{\ξ}}_f\left(t\right)=\underset{j=0}{\overset{n_0}{\mathrm{\Σ}}}{\mathrm{\α}}_j\mathrm{\Δ}{y}_f(t-j)-\underset{j=0}{\overset{n_0}{\mathrm{\Σ}}}{\mathrm{\β}}_j\mathrm{\Δ}{u}_f(t-j)-\underset{j=0}{\overset{n_0}{\mathrm{\Σ}}}{\mathrm{\γ}}_j\mathrm{\Δ}{v}_f(t-j)$

其中α₀＝1，β₀＝γ₀＝0

它有LS结构

$\mathrm{\Δ}{y}_f\left(t\right)={\mathrm{\Φ}}_1^T(t-1){\mathrm{\θ}}_1+{\mathrm{\ξ}}_f\left(t\right)$

其中 ${\mathrm{\Φ}}_1^T(t-1)=[-\mathrm{\Δ}{y}_f(t-1),\·\·\·,-\mathrm{\Δ}{y}_f(t-n_0),\mathrm{\Δ}{u}_f(t-1),\·\·\·,\mathrm{\Δ}{u}_f(t-n_0)$

$\mathrm{\Δ}{v}_f(t-1),\·\·\·,-\mathrm{\Δ}{v}_f(t-n_0)]$

${\mathrm{\θ}}_1={[{\mathrm{\α}}_1,\·\·\·,{\mathrm{\α}}_{n_0},{\mathrm{\β}}_1,\·\·\·,{\mathrm{\β}}_{n_0},{\mathrm{\γ}}_1,\·\·\·,{\mathrm{\γ}}_{n_0}]}^T$

于是用下面改进的RLS算法可得θ₁的LS估值

${\widehat{\mathrm{\θ}}}_1\left(t\right)={\widehat{\mathrm{\θ}}}_1(t-1)+\frac{P(t-1){\mathrm{\Φ}}_1(t-1)[\mathrm{\Δ}{y}_f\left(t\right)-{\mathrm{\Φ}}_1^T(t-1){\widehat{\mathrm{\θ}}}_1(t-1)]}{\mathrm{\β}\left(t\right)+{\mathrm{\Φ}}_1^T(t-1)P(t-1){\mathrm{\Φ}}_1(t-1)}$

$\mathrm{\β}\left(t\right)=1-\frac{1}{N\left(t\right)}$

$N\left(t\right)=1+\frac{1+{\mathrm{\Φ}}_1(t-1)P(t-1){\mathrm{\Φ}}_1^T(t-1)}{{[\mathrm{\Δ}{y}_f\left(t\right)-{\mathrm{\Φ}}_1^T(t-1){\widehat{\mathrm{\θ}}}_1(t-1)]}^2}\mathrm{\η}\left(t\right)$

记 $\mathrm{\ψ}\left(t\right)=P(t-1)-\frac{P(t-1){\mathrm{\Φ}}_1(t-1){\mathrm{\Φ}}_1^T(t-1)P^T(t-1)}{1+{\mathrm{\Φ}}_1^T(t-1)P(t-1){\mathrm{\Φ}}_1(t-1)}$

如果                            trace[ψ(t)]/β(t)≤ρ

则                              P(t)＝ψ(t)/β(t)

否则                            P(t)＝ψ(t)

参数η(t)与ρ的选择如下

η(t)＝mη_σ(t)

ρ＞α²

m是记忆长度，一般可取为1000及以上，η_σ是η(t)的均方差，α²是初始协方差矩阵P(0)＝α²I中的取值(I为单位阵)， $\widehat{\mathrm{\θ}}\left(0\right)=\mathrm{\ϵ}$ (α为尽可能大的数，一般可取100以上，ε为充分小的实向量)

一般取 $\left\{\begin{array}{c}\widehat{\mathrm{\θ}}\left(0\right)=\stackrel{\→}{0}\\ P\left(0\right)={\mathrm{\α}}^{2}I\end{array}\right.$

trace［ψ(t)]是矩阵ψ(t)的迹。

从而可得ξ_f(j)的平滑估值。

${\widehat{\mathrm{\ξ}}}_f\left(j\right)=\mathrm{\Δ}{y}_f\left(j\right)-{\mathrm{\Φ}}_1^T\left(j\right){\widehat{\mathrm{\θ}}}_1(t+1),j=t-1,\·\·\·,t-n.$

${\widehat{\mathrm{\σ}}}_{{\mathrm{\ξ}}_f}^2\left(t\right)={\widehat{\mathrm{\σ}}}_{{\mathrm{\ξ}}_f}^2(t-1)+\frac{1}{t}[{{\widehat{\mathrm{\ξ}}}_f}^2\left(t\right)-{\widehat{\mathrm{\σ}}}_{{\mathrm{\ξ}}_f}^2(t-1)]$

将平滑估值 ${\widehat{\mathrm{\ξ}}}_f(t-1),\·\·\·,{\widehat{\mathrm{\ξ}}}_f(t-n_c)$ 代入Φ^T(t-1)中，得到用时变遗忘因子两段参数估计法给出的参数估值

时变遗忘因子两段参数估计解决了控制参数与噪声干扰紧耦合时参数估值变慢的问题，算法随着系统动态特征的变化调整遗忘因子，增大或减小记忆长度，提高了系统的稳定性。

所述基于广义预测控制在线关键参数测量方法，就是利用带滤波器的广义预测控制自校正算法，计算、输出关键参数测量值，具体如下：

(1)读入参数N、M、λ、P(z^-1)、Q(z^-1)、T(z^-1)及参数估计算法中初始值α、η_σ(t)、m单位阵I。

(2)t＝0，P(t)＝α²I， $\widehat{\mathrm{\θ}}\left(t\right)=0;$

(3)用公式

${\widehat{\mathrm{\θ}}}_1\left(t\right)={\widehat{\mathrm{\θ}}}_1(t-1)+\frac{P(t-1){\mathrm{\Φ}}_1^T\mathrm{\Φ}(t-1)[\mathrm{\Δ}{y}_f\left(t\right)-{\mathrm{\Φ}}_1^T(t-1){\widehat{\mathrm{\θ}}}_1(t-1)]}{\mathrm{\β}\left(t\right)+{\mathrm{\Φ}}_1(t-1)P(t-1){\mathrm{\Φ}}_1^T(t-1)}$

$\mathrm{\β}\left(t\right)=1-\frac{1}{N\left(t\right)}$

$N\left(t\right)=1+\frac{1+{\mathrm{\Φ}}_1^T(t-1)P(t-1){\mathrm{\Φ}}_1(t-1)}{{[\mathrm{\Δ}{y}_f\left(t\right)-{\mathrm{\Φ}}_1^T(t-1){\widehat{\mathrm{\θ}}}_1(t-1)]}^2}\mathrm{\η}\left(t\right)$

记 $\mathrm{\ψ}\left(t\right)=P(t-1)-\frac{P(t-1){\mathrm{\Φ}}_1(t-1){\mathrm{\Φ}}_1^T(t-1)P(t-1)}{1+{\mathrm{\Φ}}_1^T(t-1)P(t-1){\mathrm{\Φ}}_1(t-1)}$

如果                                 trace[ψ(t)]/β(t)≤ρ

则                                   P(t)＝ψ(t)/β(t)

否则                                 P(t)＝ψ(t)

参数η(t)与ρ的选择同上，初值P(0)及

选法同上，从而可得ξ_f(j)的平滑估值。

${\widehat{\mathrm{\ξ}}}_f\left(j\right)=\mathrm{\Δ}{y}_f\left(j\right)-{\mathrm{\Φ}}_1^T\left(j\right){\widehat{\mathrm{\θ}}}_1(t+1)$

j＝t-1，…，t-n_c

给出

(4)将

代入公式

Φ^T(t-1)＝[-Δy_f(t-1)，…，-Δy_f(t-n_a)，Δu_f(t-1)，Δu_f(t-2)，…，

Δu_f(t-n_b-1)，Δv_f(t-1)，…，Δv_f(t-n_d-1)， ${\widehat{\mathrm{\ξ}}}_f(t-1),\·\·\·,{\widehat{\mathrm{\ξ}}}_f(t-n_c)]$ 的Φ^T(t-1)中。

利用公式 $N\left(t\right)=1+\frac{1+{\mathrm{\Φ}}^T(t-1)p(t-1)\mathrm{\Φ}(t-1)}{{[\mathrm{\Δ}{y}_f\left(t\right)-{\mathrm{\Φ}}^T(t-1)\widehat{\mathrm{\θ}}(t-1)]}^2}\mathrm{\η}\left(t\right)$

$\mathrm{\β}\left(t\right)=1-\frac{1}{N\left(t\right)}$

$\widehat{\mathrm{\θ}}\left(t\right)=\widehat{\mathrm{\θ}}(t-1)+\frac{p(t-1)\mathrm{\Φ}(t-1)[\mathrm{\Δ}{y}_f\left(t\right)-{\mathrm{\Φ}}^T(t-1)\widehat{\mathrm{\θ}}(t-1)]}{\mathrm{\β}\left(t\right)+{\mathrm{\Φ}}^T(t-1)p(t-1)\mathrm{\Φ}(t-1)}$

引进 $\mathrm{\ψ}\left(t\right)=p(t-1)-\frac{p(t-1)\mathrm{\Φ}(t-1){\mathrm{\Φ}}^T(t-1)p(t-1)}{1+{\mathrm{\Φ}}^T(t-1)p(t-1)\mathrm{\Φ}(t-1)}$

如果                               trace[ψ(t)]/β(t)≤ρ

则                                 p(t)＝ψ(t)/β(t)

估计参数 $\hat{A}\left(z^{-1}\right),\hat{B}\left(z^{-1}\right),\hat{C}\left(z^{-1}\right),\hat{D}\left(z^{-1}\right);$

(5)利用递推优化算法得到最优值

(6)利用公式

$\mathrm{\Δ}{u}_f\left(t\right)=[1,0,...,0]\stackrel{\→}{u}$

则在t时刻的控制量为：

u(t)＝T(Z^-1)[u_f(t-1)+Δu_f(t)]

计算t时刻的控制量(鼓风量)u(t)，把u(t)作为历史数据存储；

利用t时刻u(t)值计算被测量(烟气含氧量)y(t)的值；

输出烟气含氧量软测量值；

(7)置t＝t+1，返回步(2)。

本发明的有益效果是：本发明根据广义预测控制和软测量技术的各自优点进行有机结合，以充分发挥各自的优点；提出的实时在线测量模型很好的利用了广义预测控制中的思想，实现一些时延长或难以用普通传感器测量的关键参量的实时测量，为在线关键参数测量方法提供理论基础，也进一步保证了在线参数测量技术能在实际过程中得到成功应用；本发明提出的基于广义预测的含氧量软测量技术为锅炉燃烧过程中其它关键参量如烟气温度的准确、实时、在线软测量开辟了新的途径，对实现燃烧系统的闭环控制和优化运行具有重要的意义；实现了工业过程关键参量实时在线测量，解决了锅炉控制中测量与控制的时域适配问题。

附图说明

下面结合附图和实施例对本发明进一步说明；

图1为烟气含氧量测控图；

图2为模型校正参数估计法框图；

图3为广义预测自校正测量程序流程图；

图4为新型算法含氧量测试仿真图；

图5为传统算法含氧量测试仿真图。

具体实施方案

实施例：基于新型广义预测控制的锅炉燃烧过程中关键参数的在线测量方法。

1、烟气含氧量测控模型的建立

保持锅炉燃烧过程的最佳状态和经济性是锅炉燃烧过程自动控制的重要任务。在锅炉运行中，必须控制好进入锅炉的空气量与燃料量两者的比例。若空气量与燃料量之比相对较小，则化学未完全燃烧损失大，反之，则排烟热损失增大。为了使锅炉保持最佳燃烧工况，必须使空气量与燃烧量的比例合适，这个比例称为过剩空气系数，数值应在1.20～1.30之间。烟气中含氧量与过剩空气系数有确定(单值)的函数关系

$\mathrm{\α}=\frac{{\mathrm{\β}}_0}{{\mathrm{\β}}_0-\mathrm{\β}}$

式中β₀＝20.9％，为标准空气中氧的体积分数。

烟气的氧含量是锅炉控制的关键参量，它表征燃烧工况优劣，也是控制进入炉膛的空气量、维持最佳风煤比、达到优化燃烧的根本依据。

在计算机优化各种锅炉燃烧过程中，通常根据烟气含氧量来控制调节风量，同时参考调节燃料量m，噪声干扰信号，用以控制空气过剩系数α或最佳氧含量β，烟气含氧量测控如图1所示。

下面给出基于广义预测CARIMA预测模型构建的烟气含氧量测控模型：

A(z^-1)y(t)＝B(z^-1)Δu(t-1)+D(z^-1)Δv(t-1)+C(z^-1)ξ(t)

式中y、u、v、ξ分别表示输出量(被测量-烟气含氧量)、输入量(控制量-鼓风量)、前馈(可测扰动-进煤量)和正态白噪声干扰。

引入辅助输出(t)＝P(z^-1)y(t)

考虑Diophantine(丢番图方程)给出目标函数如下：

$J=E\{{(\mathrm{\ψ}-w)}^T(\mathrm{\ψ}-w)+\mathrm{\λ}{\mathrm{\ψ}}^T\mathrm{\ψ}\}$

$=E\{{(G\stackrel{\→}{u}+\stackrel{\→}{f}+\stackrel{\→}{e}-\stackrel{\→}{w})}^T(G\stackrel{\→}{u}+\stackrel{\→}{f}+\stackrel{\→}{e}-\stackrel{\→}{w})$

$+\mathrm{\λ}{(P\stackrel{\→}{u}+\stackrel{\→}{\mathrm{\δ}})}^T(P\stackrel{\→}{u}+\stackrel{\→}{\mathrm{\δ}})\}$

2、广义预测软仪表的设计

广义预测软仪表设计如下，参数辨识设计及软仪表测量如图2、3所示。

(1)读入参数N＝20、M＝20、λ＝12、P(z^-1)、Q(z^-1)、T(z^-1)及参数估计算法中初始值α＝100、m＝100单位阵I、η_σ(t)。

(2)t＝0，P(t)＝α²I， $\widehat{\mathrm{\θ}}\left(t\right)=0;$

(3)用公式

${\widehat{\mathrm{\θ}}}_1\left(t\right)={\widehat{\mathrm{\θ}}}_1(t-1)+\frac{P(t-1){\mathrm{\Φ}}_1^T(t-1)[\mathrm{\Δ}{y}_f\left(t\right)-{\mathrm{\Φ}}_1^T(t-1){\widehat{\mathrm{\θ}}}_1(t-1)]}{\mathrm{\β}\left(t\right)+{\mathrm{\Φ}}_1(t-1)P(t-1){\mathrm{\Φ}}_1^T\mathrm{\Φ}(t-1)}$

$\mathrm{\β}\left(t\right)=1-\frac{1}{N\left(t\right)}$

$N\left(t\right)=1+\frac{1+{\mathrm{\Φ}}_1^T(t-1)P(t-1){\mathrm{\Φ}}_1(t-1)}{{[\mathrm{\Δ}{y}_f\left(t\right)-{\mathrm{\Φ}}_1^T(t-1){\widehat{\mathrm{\θ}}}_1(t-1)]}^2}\mathrm{\η}\left(t\right)$

记 $\mathrm{\ψ}\left(t\right)=P(t-1)-\frac{P(t-1){\mathrm{\Φ}}_1(t-1){\mathrm{\Φ}}_1^T(t-1)P(t-1)}{1+{\mathrm{\Φ}}_1^T(t-1)P(t-1){\mathrm{\Φ}}_1(t-1)}$

如果                             trace[ψ(t)]/β(t)≤ρ

则                               P(t)＝ψ(t)/β(t)

否则                             P(t)＝ψ(t)

参数η(t)与ρ的选择同上，初值P(0)及 选法同上，从而可得ξ_f(j)的平滑估值。

${\widehat{\mathrm{\ξ}}}_f\left(j\right)=\mathrm{\Δ}{y}_f\left(j\right)-{\mathrm{\Φ}}_1^T\left(j\right){\widehat{\mathrm{\θ}}}_1(t+1),j=t-1,\·\·\·,t-n_c$

给出

(4)将 代入公式

Φ^T(t-1)＝[-Δy_f(t-1)，…，-Δy_f(t-n_a)，Δu_f(t-1)，Δu_f(t-2)，…，

Δu_f(t-n_b-1)，Δv_f(t-1)，…，Δv_f(t-n_d-1)， ${\widehat{\mathrm{\ξ}}}_f(t-1),\·\·\·,{\widehat{\mathrm{\ξ}}}_f(t-n_c)]$ 的Φ^T(t-1)中，利用公式

$N\left(t\right)=1+\frac{1+{\mathrm{\Φ}}^T(t-1)p(t-1)\mathrm{\Φ}(t-1)}{{[\mathrm{\Δ}{y}_f\left(t\right)-{\mathrm{\Φ}}^T(t-1)\widehat{\mathrm{\θ}}(t-1)]}^2}\mathrm{\η}\left(t\right)$

$\mathrm{\β}\left(t\right)=1-\frac{1}{N\left(t\right)}$

$\widehat{\mathrm{\θ}}\left(t\right)=\widehat{\mathrm{\θ}}(t-1)+\frac{p(t-1)\mathrm{\Φ}(t-1)[\mathrm{\Δ}{y}_f\left(t\right)-{\mathrm{\Φ}}^T(t-1)\widehat{\mathrm{\θ}}(t-1)]}{\mathrm{\β}\left(t\right)+{\mathrm{\Φ}}^T(t-1)p(t-1)\mathrm{\Φ}(t-1)}$

引进 $\mathrm{\ψ}\left(t\right)=p(t-1)-\frac{p(t-1)\mathrm{\Φ}(t-1){\mathrm{\Φ}}^T(t-1)p(t-1)}{1+{\mathrm{\Φ}}^T(t-1)p(t-1)\mathrm{\Φ}(t-1)}$

如果                             trace[ψ(t)]/β(t)≤ρ

则                               p(t)＝ψ(t)/β(t)

估计参数 $\hat{A}\left(z^{-1}\right),\hat{B}\left(z^{-1}\right),\hat{C}\left(z^{-1}\right),\hat{D}\left(z^{-1}\right);$

(5)利用递推优化算法得到最优值

记 $G^{T}G+\mathrm{\λ}{P}^{T}P=\left[\begin{array}{cccc}{a}_{11}& a_{12}& \·\·\·& a_{1M}\\ a_{21}& a_{22}& \·\·\·& a_{2M}\\ \·\·\·& \·\·\·& \·\·\·& \·\·\·\\ a_{M1}& a_{M2}& \·\·\·& a_{\mathrm{MM}}\end{array}\right]=\left[\begin{array}{c}{{\stackrel{\→}{a}}_1}^t\\ {{\stackrel{\→}{a}}_2}^t\\ \·\\ \·\\ \·\\ {{\stackrel{\→}{a}}_M}^t\end{array}\right]$

$G^T(\stackrel{\→}{w}-\stackrel{\→}{f})-\mathrm{\λ}{P}^T\stackrel{\→}{\mathrm{\δ}}={\stackrel{\→}{b}}^t=[{\stackrel{\→}{b}}_1,\·\·\·,{\stackrel{\→}{b}}_M]$

①读入M＝20，t＝1， ${\stackrel{\→}{a}}_i^t=(a_{i1},\·\·\·,a_{\mathrm{iM}}),$ i＝1，2，…，M

e₁＝(1，0，…，0)，…，e_j＝(0，…，0，1，0，…，0)，…，e_M＝(0，0，…，1)， ${\stackrel{\→}{b}}^t=[{\stackrel{\→}{b}}_1,\·\·\·,{\stackrel{\→}{b}}_M],$ ${\stackrel{\→}{u}}^T=(u_1,\·\·\·,u_M),$ u_i＝Δu_f(t+i-1)，i＝1，2，…，M

②任取 ${\stackrel{\→}{u}}^1\∈E^M,$ 这里不妨取 ${\stackrel{\→}{u}}^1=(0,\·\·\·,0)\∈E^M,$ 取H₁∈E^M×M为任意非奇异对称阵，这里不妨取 $H_1^T=(e_1,\·\·\·,e_M),$ 置i＝1，iflag＝0；

③计算在

处残差向量的第i个分量， ${\mathrm{\τ}}_i={\stackrel{\→}{a}}_i^T{\stackrel{\→}{u}}^i-{\stackrel{\→}{b}}_i,$ 计算搜索向量 $S^i=H_i{\stackrel{\→}{a}}_i;$

④若 $S^i\≠\stackrel{\→}{0},$ 转步⑤；若 $S^i=\stackrel{\→}{0},$ 同时 ${\mathrm{\τ}}_i=\stackrel{\→}{0},$ 则置 ${\stackrel{\→}{u}}^{i+1}={\stackrel{\→}{u}}^i,$ H_i+1＝H_i，iflag＝iflag+1，若i＜M，转步⑦；否则停止计算；此时

为方程式的解。若 $S^i=\stackrel{\→}{0},$ 但 ${\mathrm{\τ}}^i\≠\stackrel{\→}{0},$ 则置iflag＝-i停止计算，此时方程组式不相容，也即目标函数式无最优解；

⑤修正解的近似值， ${\stackrel{\→}{u}}^{i+1}={\stackrel{\→}{u}}^i-{\mathrm{\λ}}_i{S}^i,$ 其中步长 ${\mathrm{\λ}}_i={\mathrm{\τ}}_i/{\stackrel{\→}{a}}_i^T{S}^i,$

若i＝M，停止计算，

为方程式的解；

⑥修正矩阵H_i， $H_{i+1}=H_i-S^i{\left(S^i\right)}^T/{\stackrel{\→}{a}}_i^T{S}^i;$

⑦置i＝i+1，转步③；

$\stackrel{\→}{u}=u^{M+1}$ 停止，此时得到方程式的解，也即得到了使得目标函数式中J取得最小值的最优解。

(6)利用公式 ${\mathrm{\Δu}}_f\left(t\right)=[\mathrm{1,0},\·\·\·,0]\stackrel{\→}{u}$

则在t时刻的控制量为u(t)＝T(Z^-1)[u_f(t-1)+Δu_f(t)]

计算t时刻的控制量(鼓风量)u(t)，把u(t)作为历史数据存储；利用t时刻u(t)值计算被测量(烟气含氧量)y(t)的值；输出烟气含氧量软测量值；

烟气含氧量控制区间为3.5％～4.0％。

(7)置t＝t+1，返回步(2)。

3、仿真结果

应用基于新型广义预测算法和传统广义预测算法的软仪表对烟气含氧量进行测量，情况如图4、图5所示。在同样的计算机环境下进行研究的结果如下：M取20时，新旧算法的计算时间分别0.02s和0.209s。显然，系统的维数越大，即预测时域和控制时域加大，新算法节省的时间越多，稳定性更好。围绕给定值，图5所示传统算法测量烟气含氧量明显有滞后，滞后约100步左右，且波动范围较大，仿真结果表明利用新型广义预测算法构建的软仪表具有快速跟踪和输出稳定的优点。

如果采用氧化锆传感器测量，那么系统滞后一般接近60s以上。0～60s内测量时，由于控制量(风量)的改变和给煤的扰动，以及炉内燃烧过程的时滞特点，使得氧化锆测量值不能真实反映炉内燃烧情况，依据氧化锆实测值来实施控制，对于阶跃扰动，超调量高达50％～300％，时滞时间长达60s以上；而采用广义预测软测量方法测得数据超调量很小，能够反映出含氧量变化的实际情况。60s～120s内测量时，由于炉内燃烧已经达到稳定状态，软硬仪表的实测值是相近的，其误差是可以忽略不计的。

因此，运用新型广义预测控制算法构建的软仪表能够客观反映锅炉燃烧过程中任何时段的真实情况，提高了烟气含氧量的测量精度和实时性，优化了锅炉燃烧控制系统的性能。随着广义预测控制理论的发展和预测算法的改进，测量精度还可进一步提高，基于广义预测的含氧量软测量技术为工业锅炉的氧量测量提供了新的手段，也为锅炉燃烧过程中其它关键参量如烟气温度的准确、实时、在线软测量开辟了新的途径，对实现燃烧系统的闭环控制和优化运行具有重要意义。

本发明同样适用于毛坯加热炉、煤气生产炉等生产过程中烟气含氧量、CO烟气含量的测量。

Claims (5)

1、一种基于新型广义预测控制的在线关键参数测量方法，其特征是：根据过程控制的时滞特性和关键参数的测量方法，利用新型广义预测控制器，建立了基于广义预测算法的软测量预测模型，实现了过程控制关键参量的实时在线测量，具体包括测量模型建立、初始化、新型预测控制方法、反馈校正和基于新型广义预测控制的在线关键参数测量方法四个步骤：

(1)建立测量模型：根据广义预测控制技术特征和过程控制时滞特性建立测量模型；

(2)初始化及求解最优控制律：根据测量模型设定控制参数，利用新型广义预测控制方法求解最优控制律；

(3)模型反馈校正：采用时变遗忘因子两段参数估计方法进行模型反馈校正，优化控制参数；

(4)基于广义预测控制在线关键参数测量方法：根据新型广义预测控制器计算、输出关键参数测量值。

2、根据权利要求1所述的基于新型广义预测控制的在线关键参数测量方法，其特征是：所述建立的测量模型，具体如下：

确定系统的过程控制的目标函数为：

$\begin{array}{c}J=E\{{(\mathrm{\ψ}-w)}^T(\mathrm{\ψ}-w)+{\mathrm{\λ\ψ}}^T\mathrm{\ψ}\}\\ =E\{{(G\stackrel{\→}{u}+\stackrel{\→}{f}+\stackrel{\→}{e}-\stackrel{\→}{w})}^T(G\stackrel{\→}{u}+\stackrel{\→}{f}+\stackrel{\→}{e}-\stackrel{\→}{w})+\mathrm{\λ}{(P\stackrel{\→}{u}+\stackrel{\→}{\mathrm{\δ}})}^T(P\stackrel{\→}{u}+\stackrel{\→}{\mathrm{\δ}})\}\end{array}$

为了求得最优控制律 由目标函数式可得

$(G^{T}G+{\mathrm{\λP}}^{T}P)\stackrel{\→}{u}=G^T(\stackrel{\→}{w}-\stackrel{\→}{f}-\stackrel{\→}{e})-{\mathrm{\λP}}^T\stackrel{\→}{\mathrm{\δ}}$

由数学分析可知由方程所求得的

即为使得目标函数式中的J取得最小值的最优解。

3、根据权利要求1所述的基于新型广义预测控制的在线关键参数测量方法，其特征是：初始化及求解最优控制律，具体如下：

由于新型控制算法采用矩阵递推迭代方法，避免了传统方法求解最优控制律中逆矩阵的计算，运算复杂度下降了M(矩阵阶数)的平方阶次，节省的计算量和存储量这样就可以扩大预测时域和控制时域的选取范围，增加测控系统的稳定性，可以实现锅炉关键参数实时在线软测量；

递推优化算法如下：

记 $G^{T}G+{\mathrm{\λP}}^{T}P=\left[\begin{array}{cccc}{a}_{11}& a_{12}& \·\·\·& a_{1M}\\ a_{21}& a_{22}& \·\·\·& a_{2M}\\ \·\·\·& \·\·\·& \·\·\·& \·\·\·\\ a_{M1}& a_{M2}& \·\·\·& a_{\mathrm{MM}}\end{array}\right]=\left[\begin{array}{c}{{\stackrel{\→}{a}}_1}^t\\ {{\stackrel{\→}{a}}_2}^t\\ \·\\ \·\\ \·\\ {{\stackrel{\→}{a}}_M}^t\end{array}\right]$

$G^T(\stackrel{\→}{w}-\stackrel{\→}{f})-{\mathrm{\λP}}^T\stackrel{\→}{\mathrm{\δ}}={\stackrel{\→}{b}}^t=[{\stackrel{\→}{b}}_1,\·\·\·,{\stackrel{\→}{b}}_M]$

(1)读入 $M,t,{\stackrel{\→}{a}}_i^t=(a_{i1},\·\·\·,a_{\mathrm{iM}}),i=\mathrm{1,2},\·\·\·,M$

e₁＝(1，0，...，0)，...，e_j＝(0，...，0，1，0，...，0)，...，e_M＝(0，0，...，1)， ${\stackrel{\→}{b}}^t=[{\stackrel{\→}{b}}_1,\·\·\·,{\stackrel{\→}{b}}_M],$ ${\stackrel{\→}{u}}^T=(u_1,\·\·\·,u_M),$ u_i＝Δu_f(t+i-1)，i＝1，2，...，M

(2)任取 ${\stackrel{\→}{u}}^1\∈E^M,$ 这里不妨取 ${\stackrel{\→}{u}}^1=(0,\·\·\·,0)\∈E^M,$ 取H₁∈E^M×M为任意非奇异对称阵，这里不妨取 $H_1^T=(e_1,\·\·\·,e_M),$ 置i＝1，iflag＝0；

(3)计算在

处残差向量的第i个分量， ${\mathrm{\τ}}_i={\stackrel{\→}{a}}_i^T{\stackrel{\→}{u}}^i-{\stackrel{\→}{b}}_i,$ 计算搜索向量 $S^i=H_i{\stackrel{\→}{a}}_i;$

(4)若 $S^i\≠\stackrel{\→}{0},$ 转步(5)；若 $S^i=\stackrel{\→}{0},$ 同时 ${\mathrm{\τ}}_i=\stackrel{\→}{0},$ 则置 ${\stackrel{\→}{u}}^{i+1}={\stackrel{\→}{u}}^i,$ H_i+1＝H_i，iflag＝iflag+1，若i＜M，转步(7)；否则停止计算；此时

为方程式的解。若 $S^i=\stackrel{\→}{0},$ 但 ${\mathrm{\τ}}^i\≠\stackrel{\→}{0},$ 则置iflag＝-i停止计算，此时方程组式不相容，也即目标函数式无最优解；

(5)修正解的近似值， ${\stackrel{\→}{u}}^{i+1}={\stackrel{\→}{u}}^i-{\mathrm{\λ}}_i{S}^i,$ 其中步长 ${\mathrm{\λ}}_i={\mathrm{\τ}}_i/{\stackrel{\→}{a}}_i^T{S}^i,$ 若i＝M，停止计算，

为方程式的解；

(6)修正矩阵H_i， $H_{i+1}=H_i-S^i{\left(S^i\right)}^T/{\stackrel{\→}{a}}_i^T{S}^i;$

(7)置i＝i+1，转步(3)；

$\stackrel{\→}{u}=u^{M+1}$ 停止，此时得到方程式的解，也即得到了使得目标函数式中J取得最小值的最优解。从而

${\mathrm{\Δu}}_f\left(t\right)=[\mathrm{1,0},\·\·\·,0]\stackrel{\→}{u}$

则在t时刻的控制量为

u(t)＝T(Z^-1)[u_f(t-1)+Δu_f(t)]

4、根据权利要求1所述的基于广义预测控制的在线关键参数测量方法，其特征是：初始化及求解最优控制律，模型的反馈校正，是利用此采样时刻的广义预测模型计算最优值和给定值之间的偏差实时在线修正模型参数，这种基于广义预测控制模型，根据系统的实际输出不断对预测控制输出值进行滚动优化修正，而且利用了反馈信息，构成闭环优化校正；反馈校正参数估计采用时变遗忘因子两段参数估计性能指标加权的广义预测自适应的优化控制算法，具体如下：

预测模型采用

A(z^-1)Δy_f(t)＝B(z^-1)Δu_f(t-1)+D(z^-1)Δv_f(t-1)+C(z^-1)ξ_f(t)

定理若A(z^-1)和C(z^-1)是稳定的，则预测模型式可用如下高阶CAR模型近似代替

${\mathrm{\ξ}}_f\left(t\right)=\underset{j=0}{\overset{n_0}{\mathrm{\Σ}}}{\mathrm{\α}}_j{\mathrm{\Δy}}_f(t-j)-\underset{j=0}{\overset{n_0}{\mathrm{\Σ}}}{\mathrm{\β}}_j{\mathrm{\Δu}}_f(t-j)-\underset{j=0}{\overset{n_0}{\mathrm{\Σ}}}{\mathrm{\γ}}_j{\mathrm{\Δv}}_f(t-j)$

其中α₀＝1，β₀＝γ₀＝0

它有LS结构

${\mathrm{\Δy}}_f\left(t\right)={{\mathrm{\Φ}}_1}^T(t-1){\mathrm{\θ}}_1+{\mathrm{\ξ}}_f\left(t\right)$

其中 ${{\mathrm{\Φ}}_1}^T(t-1)=[-\mathrm{\Δ}{y}_f(t-1),\·\·\·,-\mathrm{\Δ}{y}_f(t-n_0),$ Δu_f(t-1)，...，Δu_f(t-n₀)

                   Δv_f(t-1)，...，-Δv_f(t-n₀)]

${\mathrm{\θ}}_1={[{\mathrm{\α}}_1,\·\·\·,{\mathrm{\α}}_{n_0},{\mathrm{\β}}_1,\·\·\·,{\mathrm{\β}}_{n_0},{\mathrm{\γ}}_1,\·\·\·,{\mathrm{\γ}}_{n_0}]}^T$

于是用下面改进的RLS算法可得θ₁的LS估值

${\widehat{\mathrm{\θ}}}_1\left(t\right)={\widehat{\mathrm{\θ}}}_1(t-1)+\frac{P(t-1){\mathrm{\Φ}}_1(t-1)[\mathrm{\Δ}{y}_f\left(t\right)-{\mathrm{\Φ}}_1^T(t-1){\widehat{\mathrm{\θ}}}_1(t-1)]}{\mathrm{\β}\left(t\right)+{{\mathrm{\Φ}}_1}^T(t-1)P(t-1){\mathrm{\Φ}}_1(t-1)}$

$\mathrm{\β}\left(t\right)=1-\frac{1}{N\left(t\right)}$

$N\left(t\right)=1+\frac{1+{\mathrm{\Φ}}_1(t-1)P(t-1){\mathrm{\Φ}}_1^T(t-1)}{{[{\mathrm{\Δy}}_f\left(t\right)-{\mathrm{\Φ}}_1^T(t-1){\widehat{\mathrm{\θ}}}_1(t-1)]}^2}\mathrm{\η}\left(t\right)$

记 $\mathrm{\ψ}\left(t\right)=P(t-1)-\frac{P(t-1){\mathrm{\Φ}}_1(t-1){\mathrm{\Φ}}_1^T(t-1)P^T(t-1)}{1+{\mathrm{\Φ}}_1^T(t-1)P(t-1){\mathrm{\Φ}}_1(t-1)}$

如果                        trace[ψ(t)]/β(t)≤ρ

则                          P(t)＝ψ(t)/β(t)

否则                        P(t)＝ψ(t)

参数η(t)与ρ的选择如下：

η(t)＝mη_σ(t)

ρ＞α²

m是记忆长度，一般可取为1000及以上，η_σ是η(t)的均方差，α²是初始协方差矩阵p(0)＝α²I中的取值(I为单位阵)， $\widehat{\mathrm{\θ}}\left(0\right)=\mathrm{\ϵ}$ (α为尽可能大的数，一般可取100以上，ε为充分小的实向量)

一般取 $\left\{\begin{array}{c}\widehat{\mathrm{\θ}}\left(0\right)=\stackrel{\→}{0}\\ p\left(0\right)={\mathrm{\α}}^{2}I\end{array}\right.$

trace[ψ(t)]是矩阵ψ(t)的迹。

从而可得ξ_f(j)的平滑估值。

${\widehat{\mathrm{\ξ}}}_f\left(j\right)={\mathrm{\Δy}}_f\left(j\right)-{\mathrm{\Φ}}_1^T\left(j\right){\widehat{\mathrm{\θ}}}_1(t+1),j=t-1,\·\·\·,t-n_c$

${\widehat{\mathrm{\σ}}}_{{\mathrm{\ξ}}_f}^2\left(t\right)={\widehat{\mathrm{\σ}}}_{{\mathrm{\ξ}}_f}^2(t-1)+\frac{1}{t}[{\widehat{\mathrm{\ξ}}}_f^2\left(t\right)-{\widehat{\mathrm{\σ}}}_{{\mathrm{\ξ}}_f}^2(t-1)]$

将平滑估值

代入Φ^T(t-1)中，得到用时变遗忘因子两段参数估计法给出的参数估值

时变遗忘因子两段参数估计解决了控制参数与噪声干扰紧耦合时参数估值变慢的问题，算法随着系统动态特征的变化调整遗忘因子，增大或减小记忆长度，提高了系统的稳定性。

5、根据权利要求1所述的基于广义预测控制的在线关键参数测量方法，其特征是：基于广义预测控制在线关键参数测量方法，就是利用带滤波器的广义预测控制自校正算法，计算、输出关键参数测量值，具体如下：

(1)读入参数N、M、λ、P(z^-1)、Q(z^-1)、T(z^-1)及参数估计算法中初始值α、η_σ(t)、m单位阵I。

(2)t＝0，P(t)＝α²I， $\widehat{\mathrm{\θ}}\left(t\right)=0;$

(3)用公式

${\widehat{\mathrm{\θ}}}_1\left(t\right)={\widehat{\mathrm{\θ}}}_1(t-1)+\frac{P(t-1){{\mathrm{\Φ}}_1}^T(t-1)[\mathrm{\Δ}{y}_f\left(t\right)-{\mathrm{\Φ}}_1^T(t-1){\widehat{\mathrm{\θ}}}_1(t-1)]}{\mathrm{\β}\left(t\right)+{\mathrm{\Φ}}_1(t-1)P(t-1){\mathrm{\Φ}}_1^T(t-1)}$

$\mathrm{\β}\left(t\right)=1-\frac{1}{N\left(t\right)}$

$N\left(t\right)=1+\frac{1+{{\mathrm{\Φ}}_1}^T(t-1)P(t-1){\mathrm{\Φ}}_1(t-1)}{{[{\mathrm{\Δy}}_f\left(t\right)-{\mathrm{\Φ}}_1^T(t-1){\widehat{\mathrm{\θ}}}_1(t-1)]}^2}\mathrm{\η}\left(t\right)$

记 $\mathrm{\ψ}\left(t\right)=P(t-1)-\frac{P(t-1){\mathrm{\Φ}}_1(t-1){\mathrm{\Φ}}_1^T(t-1)P(t-1)}{1+{\mathrm{\Φ}}_1^T(t-1)P(t-1){\mathrm{\Φ}}_1(t-1)}$

如果              trace[ψ(t)]/β(t)≤ρ

则                P(t)＝ψ(t)/β(t)

否则              P(t)＝ψ(t)

参数η(t)与ρ的选择同上，初值P(0)及

选法同上，从而可得ξ_f(j)的平滑估值。

${\widehat{\mathrm{\ξ}}}_f\left(j\right)={\mathrm{\Δy}}_f\left(j\right)-{\mathrm{\Φ}}_1^T\left(j\right){\widehat{\mathrm{\θ}}}_1(t+1)$

j＝t-1，…，t-n_c

给出

(4)将

代入公式

Φ^T(t-1)＝[-Δy_f(t-1)，...，-Δy_f(t-n_a)，Δu_f(t-1)，Δu_f(t-2)，...，

${\mathrm{\Δu}}_f(t-n_b-1),{\mathrm{\Δv}}_f(t-1),\·\·\·,{\mathrm{\Δv}}_f(t-n_d-1),{\widehat{\mathrm{\ξ}}}_f(t-1),\·\·\·,{\widehat{\mathrm{\ξ}}}_f(t-n_c)]$

的Φ^T(t-1)中。

利用公式

$N\left(t\right)=1+\frac{1+{\mathrm{\Φ}}^T(t-1)p(t-1)\mathrm{\Φ}(t-1)}{{[{\mathrm{\Δy}}_f\left(t\right)-{\mathrm{\Φ}}^T(t-1)\widehat{\mathrm{\θ}}(t-1)]}^2}\mathrm{\η}\left(t\right)$

$\mathrm{\β}\left(t\right)=1-\frac{1}{N\left(t\right)}$

$\widehat{\mathrm{\θ}}\left(t\right)=\widehat{\mathrm{\θ}}(t-1)+\frac{p(t-1)\mathrm{\Φ}(t-1)[{\mathrm{\Δy}}_f\left(t\right)-{\mathrm{\Φ}}^T(t-1)\widehat{\mathrm{\θ}}(t-1)]}{\mathrm{\β}\left(t\right)+{\mathrm{\Φ}}^T(t-1)p(t-1)\mathrm{\Φ}(t-1)}$

引进 $\mathrm{\ψ}\left(t\right)=p(t-1)-\frac{p(t-1)\mathrm{\Φ}(t-1){\mathrm{\Φ}}^T(t-1)p(t-1)}{1+{\mathrm{\Φ}}^T(t-1)p(t-1)\mathrm{\Φ}(t-1)}$

如果                  trace[ψ(t)]/β(t)≤ρ

则                    p(t)＝ψ(t)/β(t)

估计参数

(5)利用递推优化算法得到最优值

(6)利用公式

$\mathrm{\Δ}{u}_f\left(t\right)=[\mathrm{1,0},\·\·\·,0]\stackrel{\→}{u}$

则在t时刻的控制量为

                     u(t)＝T(Z^-1)[u_f(t-1)+Δu_f(t)]

计算t时刻的控制量(鼓风量)u(t)，把u(t)作为历史数据存储；

利用t时刻u(t)值计算被测量(烟气含氧量)y(t)的值；

输出烟气含氧量软测量值；

(7)置t＝t+1，返回步(2)。

Images