CN104503236A - 一种基于回归模型的分解炉出口温度滑模控制方法 - Google Patents

Abstract

本发明公开了一种基于回归模型的分解炉出口温度滑模控制方法，首先根据水泥预分解工艺流程及现场操作人员经验，选取喂煤量作为模型的输入变量，根据现场实际情况，选定分解炉出口温度840℃～860℃为典型工况；然后根据历史数据，建立基于回归分析的分解炉出口温度数学模型；最后采用自适应趋近率求取最优控制量，建立自适应滑模控制器，具有较强的鲁棒性和不变性。本发明可准确实现分解炉出口温度的控制，为实现工业现场分解炉的优化控制提供了新思路。建立自适应滑模控制器，具有较强的鲁棒性和不变性，可准确实现分解炉出口温度的控制，为实现工业现场分解炉的优化控制提供了新思路。

Description

一种基于回归模型的分解炉出口温度滑模控制方法

技术领域

本发明属于工业自动化技术领域，尤其涉及一种基于回归模型的分解炉出口温度滑模控制方法。

背景技术

我国水泥总产量位居世界首位，年产量高达20多亿吨，占全球总量的50％以上，因而利用自动化技术实现水泥行业节能降耗成为当前水泥生产的研究重点与热点。水泥生产的预分解过程是水泥生产的核心环节之一，其煤耗大，约占整个水泥生产过程煤耗总量的60％。因此，对该环节实施优化控制对水泥企业实现节能降耗具有重大的意义。

作为预分解技术的核心设备，分解炉担负着预分解系统中繁重的燃烧、热传递和物料分解的任务。由于生料预分解过程的工况条件变化频繁并且测控点少，这使得在实际生产中经常出现分解炉温度大幅波动的现象。温度过高容易引起预热器结皮，影响窑系统正常运行；温度过低，则造成入窑分解率过低，增加窑系统负担，不能充分发挥分解炉的作用。因此，分解炉出口温度的控制，既对水泥企业实现节能降耗具有重要的意义，也影响着水泥生产的正常进行。

为了实现分解炉出口温度的控制，建立合适的分解炉出口温度数学模型是十分重要的。文献(费德诺.水泥预分解煅烧过程分解炉的数值建模.粉末技术，2007，117(1)：81-85.)从反应动力学出发，建立煤粉燃烧和碳酸盐分解的动力学数学模型。并没有考察影响分解炉温度时变的主要因素与分解炉温度之间的关系。文献(姚维孟，颜丈俊，诸静.水泥回转窑分解炉温度的模糊控制[J].自动化仪表，2001，22(2)：35-37.)提出了一种分解炉温度的模糊控制方法。丈献(汪海峰，诸静.模糊预测控制在水泥生产中的应用[J].工业仪表与自动化装置，2003，02：10-13.)给出了一种改进的分解炉温度的模糊预测控制算法，并给出其与传统模糊控制的运行效果比较分析，证明其具有更好的控制精度。然而，上述文献中，模糊控制的论域是通过手动在线调整预测误差及预测误差变化量的量化因子而获得的，因此，控制量在论域中的精确值很难及时取得。

发明内容

本发明的目的在于提供一种基于回归模型的分解炉出口温度滑模控制方法，旨在解决现有的分解炉出口温度滑模控制方法存在的模糊控制论域是通过手动在线调整预测误差及预测误差变化量的量化因子而获得，控制量在论域中的精确值很难及时取得，造成方法精确度较低的问题。

本发明是这样实现的，一种基于回归模型的分解炉出口温度滑模控制方法，该基于回归模型的分解炉出口温度滑模控制方法包括以下步骤：

步骤一，根据分解炉出口温度与喂煤量之间的关系，采用回归分析算法，建立分解炉在典型工况下的数学模型；

分解炉数学模型为：

$y(n+1)={\widehat{\mathrm{\α}}}_1+{\widehat{\mathrm{\α}}}_{2}y\left(n\right)+{\widehat{\mathrm{\α}}}_{3}y(n-1)+{\widehat{\mathrm{\α}}}_{4}u\left(n\right)$

式中，y(n+1)为n+1分解炉出口温度，y(n)为n分解炉出口温度，y(n-1)为n-1分解炉出口温度，α₁，α₂，α₃，α₄为待辨识的模型参数，u(n)为n时刻喂煤量；

步骤二，根据分解炉数学模型，采用自适应趋近率求取最优控制量，建立分解炉的滑模控制器，从而实现分解炉出口温度控制；

自适应趋近率为：

$s(k+1)-s\left(k\right)=-\mathrm{qTs}\left(k\right)-\frac{s\left(k\right)}{2}\mathrm{Tsgn}\left(s\left(k\right)\right)$

式中，q＝0.1，T＝0.01，s(k)，s(k+1)为k及k+1时刻的切换函数：

s(k)＝C_eE＝C_e(R(k)-x(k))，

其中，C_e＝[10，1]，R(k)＝[r(k)，dr(k)]，r(k)为位置指令，dr(k)为r(k)的微分，x(k)为分解炉出口温度初值；

控制率为：

$u\left(k\right)={\left(C_{e}B\right)}^{-1}(C_{e}R(k+1)-C_e\mathrm{Ax}\left(k\right)-C_e\mathrm{\θ}-s\left(k\right)+\mathrm{qTs}\left(k\right)+\frac{s\left(k\right)}{2}\mathrm{Tsgn}\left(s\left(k\right)\right))$

式中， $A=\left[\begin{array}{cc}-0.0045& -0.0002\\ 1& 0\end{array}\right],$ $B=\left[\begin{array}{c}9.6468\\ 0\end{array}\right],$ $\mathrm{\θ}=\left[\begin{array}{c}718.9233\\ 0\end{array}\right],$ R(k+1)为k+1时刻的位置指令，sgn(·)为符号函数。

进一步，在步骤一之前需要得出分解炉出口温度与喂煤量之间的关系，并选取分解炉出口温度840℃～860℃作为分解炉的工况，喂煤量变化范围为13.06t～15.1t。

进一步，分解炉出口温度与喂煤量之间的关系：当分解炉处于正常工作状态时，分解炉出口温度随喂煤量的增加而升高。

进一步，在步骤一中分解炉出口温度数学模型的输入变量选取的工艺流程如下：

步骤一，生料由提升机喂入C1～C2级旋风筒的的连接管道；

步骤二，生料由来自C2级旋风筒的热风带入C1级旋风筒进行气固热交换，再由C1级旋风筒底部的锁风阀排出，进入C2～C3级旋风筒的连接管道，又被气流带入C2级旋风筒内继续气固热交换，如此反复；

步骤三，预热后的物料经由C4级旋风筒锥部进入分解炉，煤粉则从分解炉中部给煤口进入分解炉；

步骤四，从分解炉出来的物料经分解炉上部的鹅颈管进入C5级旋风筒；最后，由C5级旋风筒的锥部进入回转窑进行煅烧。

进一步，分解炉数学模型获取的方法包括：

根据分解炉出口温度工作于840℃-860℃为典型工况，即分解炉出口温度在840℃至860℃时，分解炉工作效率最高，采用回归分析算法对分解炉出口温度在840℃-860℃区间内建模；

分解炉出口温度与喂煤量采样值的内在联系是线性的，即由控制系统的测量装置对分解炉出口气体温度及喂煤量进行采样，得到n组测量数据，结构形式如下：

$\left\{\begin{array}{c}y\left(2\right)={\mathrm{\α}}_1+{\mathrm{\α}}_{2}y\left(1\right)+{\mathrm{\α}}_{3}y\left(0\right)+{\mathrm{\α}}_{4}u\left(1\right)\\ y\left(3\right)={\mathrm{\α}}_1+{\mathrm{\α}}_{2}y\left(2\right)+{\mathrm{\α}}_{3}y\left(1\right)+{\mathrm{\α}}_{4}u\left(2\right)\\ ......\\ y\left(n\right)={\mathrm{\α}}_1+{\mathrm{\α}}_{2}y(n-1)+{\mathrm{\α}}_{3}y(n-2)+{\mathrm{\α}}_{4}u(n-1)\\ y(n+1)={\mathrm{\α}}_1+{\mathrm{\α}}_{2}y\left(n\right)+{\mathrm{\α}}_{3}y(n-1)+{\mathrm{\α}}_{4}u\left(n\right)\end{array}\right.---\left(1\right)$

式中，y₀，y₂，…，y_n+1为t₀至t_n+1时刻的分解炉出口温度，u₁，u₂，…u_n为t₁至t_n时刻的喂煤量，α₁，α₂，α₃，α₄为待辨识的模型参数：

将(1)式改写成矩阵形式，得到：

y＝Aα    (2)

式中：

y＝[y(2)，y(3)，…，y(n+1)]^T

α＝[α₁，α₂，α₃，α₄]^T

$A={\left[\begin{array}{cccc}1& y\left(1\right)& y\left(0\right)& u\left(1\right)\\ 1& y\left(1\right)& y\left(1\right)& u\left(2\right)\\ .& .& .& .\\ .& .& .& .\\ .& .& .& .\\ 1& y\left(n\right)& y(n-1)& u\left(n\right)\end{array}\right]}_{n\×4}$

根据最小二乘估计，得到回归方程为：

$\hat{y}(n+1)={\mathrm{\α}}_1+{\mathrm{\α}}_{2}y\left(n\right)+{\mathrm{\α}}_{3}y(n-1)+{\mathrm{\α}}_{4}u\left(n\right)---\left(3\right)$

回归系数为：

α＝(A^TA)^-1A^Ty    (4)

由此得，基于最小二乘学习算法的分解炉出口温度数学模型为：

$y(n+1)={\widehat{\mathrm{\α}}}_1+{\widehat{\mathrm{\α}}}_{2}y\left(n\right)+{\widehat{\mathrm{\α}}}_{3}y(n-1)+{\widehat{\mathrm{\α}}}_{4}u\left(n\right).$

进一步，自适应趋近率和控制率公式的获取方法：

根据所建立的数学模型，由式得如下状态方程：

A(z^-1)y(k+1)＝B(z^-1)u(k)+θ    (7)

其中：

A(z^-1)＝1-0.0045z^-1-0.0002z^-2

B(z^-1)＝9.6468

θ＝718.9233

将式(7)改写成如下形式：

x(k+1)＝Ax(k)+Bu(k)+θ    (8)

y(k)＝Cx(k)

其中：

$A=\left[\begin{array}{cc}-0.0045& -0.0002\\ 1& 0\end{array}\right];$ $B=\left[\begin{array}{c}9.6468\\ 0\end{array}\right];$ C＝[1 0]； $\mathrm{\θ}=\left[\begin{array}{c}718.9233\\ 0\end{array}\right]$

切换函数选取为：

s(k)＝C_eE＝C_e(R(k)-x(k))    (9)

其中，C_e＝[c，1]；R＝[r(k)，dr(k)]，r(k)为位置指令，式(9)被改写成如下形式：

s(k+1)＝C_e(R(k+1)-x(k+1))

＝C_e(R(k+1)-Ax(k)-Bu(k)-θ)    (10)

＝C_eR(k+1)-C_eAx(k)-C_eBu(k)-C_eθ

控制率为：

u(k)＝(C_eB)^-1(C_eR(k+1)-C_eAx(k)-C_eθ-s(k+1))    (11)

基于指数的离散趋近率为：

s(k+1)＝s(k)+T(-ε sgn(s(k))-qs(k))    (12)

其中，ε＞0；q＞0；1-qT＞0；T为采样时间，T<<1.0，根据(12)，式(11)改写成：

u(k)＝(C_eB)^-1(C_eR(k+1)-C_eAx(k)-C_eθ-s(k)-ds(k))    (13)

其中，ds(k)＝-εT sgn(s(k))-qTs(k)；

根据式(12)，得到：

$\begin{array}{c}s(k+1)=(1-\mathrm{qT})s\left(k\right)-\mathrm{\&epsiv;T}\frac{s\left(k\right)}{\left|s\left(k\right)\right|}\\ =(1-\mathrm{qT}-\frac{\mathrm{\&epsiv;T}}{\left|s\left(k\right)\right|})s\left(k\right)=\mathrm{ps}\left(k\right)\end{array}---\left(14\right)$

由式(14)有如下三种情况：

当 $\left|s\left(k\right)\right|>\frac{\mathrm{\&epsiv;T}}{2-\mathrm{qT}},$ 有：

$p>1-\mathrm{qT}-\frac{\mathrm{\&epsiv;T}(2-\mathrm{qT})}{\mathrm{\&epsiv;T}}=1-\mathrm{qT}-(2-\mathrm{qT})=-1$

因此，|p|＜1；|s(k+1)|＜|s(k)|；|s(k)|是递减的；

当 $\left|s\left(k\right)\right|<\frac{\mathrm{\&epsiv;T}}{2-\mathrm{qT}},$ 有：

$p<1-\mathrm{qT}-\frac{\mathrm{\&epsiv;T}(2-\mathrm{qT})}{\mathrm{\&epsiv;T}}=1-\mathrm{qT}-(2-\mathrm{qT})=-1$

因此，|p|＞1；|s(k+1)|＞|s(k)|；|s(k)|是递增的；

当 $\left|s\left(k\right)\right|=\frac{\mathrm{\&epsiv;T}}{2-\mathrm{qT}},$ 有：

$p=1-\mathrm{qT}-\frac{\mathrm{\&epsiv;T}(2-\mathrm{qT})}{\mathrm{\&epsiv;T}}=1-\mathrm{qT}-(2-\mathrm{qT})=-1$

因此，|p|＝1；|s(k+1)|＝|s(k)|；|s(k)|是震荡的；

显然，|s(k)|需要满足：

$\left|s\left(k\right)\right|>\frac{\mathrm{\&epsiv;T}}{2-\mathrm{qT}}---\left(15\right)$

式(15)改写成：

$\mathrm{\&epsiv;}<\frac{1}{T}(2-\mathrm{qT})\left|s\left(k\right)\right|---\left(16\right)$

取ε＝|s(k)|/2，式(12)改写成：

自适应趋近率： $s(k+1)-s\left(k\right)=-\mathrm{qTs}\left(k\right)-\frac{s\left(k\right)}{2}\mathrm{Tsgn}\left(s\left(k\right)\right)$

控制率为：

$\begin{array}{c}u\left(k\right)={\left(C_{e}B\right)}^{-1}(C_{e}R(k+1)-C_e\mathrm{Ax}\left(k\right)-C_e\mathrm{\&theta;}\\ -s\left(k\right)+\mathrm{qTs}\left(k\right)+\frac{s\left(k\right)}{2}\mathrm{Tsgn}\left(s\left(k\right))\right)\end{array}.$

本发明提供的基于回归模型的分解炉出口温度滑模控制方法，建立自适应滑模控制器，具有较强的鲁棒性和不变性，可准确实现分解炉出口温度的控制，具体优势如下：由图4-5可知，所建模型能够准确反映分解炉出口温度变化状态，误差在±1℃以内，满足工业生产要求；由图6-8可知，所设计的滑模控制器能够准确调节喂煤量以实现分解炉的出口温度的控制，调节时间在1分钟以内，稳态误差回零，控制准确；由图9可知，自适应率ε在控制初期较大，随着时间的增加，其值逐渐减小回零，既避免了系统的抖振，又提高了系统到达切换面的速度。因此，本发明为实现工业现场分解炉的优化控制提供了新思路；

与现有技术相比，本发明的有益效果为：

1、本发明中建立了水泥分解炉出口温度在840℃-860℃时的数学模型，该温度范围为分解炉的典型工况，所建模型能够充分反映分解炉在典型工况下的温度变化情况，具有一定的针对性。

2、本发明中考虑了样本数据间量纲差异，采用了样本数据归一化，能够更好的挖掘出数据间的关系。

3、本发明中控制器的构建采用了滑模变结构控制方法，响应速度快，物理实现简单。

4、本发明中采用自适应指数趋近率求取最优控制量，系统状态不发生震荡，渐进趋向平衡点零，控制输入信号平滑无抖振。

附图说明

图1是本发明实施例提供的基于回归模型的分解炉出口温度滑模控制方法流程图；

图2是本发明实施例提供的分解炉温度控制原理图；

图3是本发明实施例提供的水泥预分解工艺流程图；

图4是本发明实施例提供的分解炉出口温度数学模型拟合曲线图；

图5是本发明实施例提供的分解炉出口温度数学模型误差曲线图；

图6是本发明实施例提供的分解炉出口温度控制曲线图；

图7是本发明实施例提供的分解炉出口温度跟踪误差曲线图；

图8是本发明实施例提供的控制量(喂煤)变化曲线图；

图9为是本发明实施例提供的自适应率ε变化曲线图。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。

下面结合附图及具体实施例对本发明的应用原理作进一步描述。

水泥分解炉温度控制原理图如图2所示；

如图1所示，本发明实施例的基于回归模型的分解炉出口温度滑模控制方法包括以下步骤：

S101：根据水泥预分解工艺流程及历史数据，得出分解炉出口温度与喂煤量之间的关系，并选取分解炉出口温度840℃～860℃作为分解炉的典型工况；

S102：根据步骤S101中所述的关系，采用回归分析算法，建立分解炉在典型工况下的数学模型；

S103：根据步骤S102中所述的数学模型，采用自适应趋近率求取最优控制量，建立分解炉的滑模控制器，从而实现分解炉出口温度控制。

下面对每个步骤作进一步详细说明：

步骤S101中：水泥预分解工艺流程及分解炉出口温度与喂煤量之间的关系；

本发明中分解炉出口温度数学模型的输入变量是根据水泥预分解工艺流程选取的，其工艺流程图如图3所示，具体可分为如下步骤：

步骤一，生料由提升机喂入C1～C2级旋风筒的的连接管道；

步骤二，生料由来自C2级旋风筒的热风带入C1级旋风筒进行气固热交换，再由C1级旋风筒底部的锁风阀排出，进入C2～C3级旋风筒的连接管道，又被气流带入C2级旋风筒内继续气固热交换，如此反复；

步骤三，预热后的物料经由C4级旋风筒锥部进入分解炉，煤粉则从分解炉中部给煤口进入分解炉；由于煤粉颗粒很小并且充分同物料混合，因此分解炉中的煤粉以无焰状态燃烧；煤粉燃烧所释放的热量被碳酸盐吸收，导致碳酸盐吸热而发生分解反应；

步骤四，从分解炉出来的物料经分解炉上部的鹅颈管进入C5级旋风筒；最后，由C5级旋风筒的锥部进入回转窑进行煅烧；

分析该工艺过程可知，在分解炉的预分解过程中，因为存在着很多反应，导致影响分解炉温度变化的变量繁多，依据分解炉结构、生产工艺及内部燃烧机理的分析可以发现在气候环境、原料的硬度、颗粒大小等因素一定的情况下，分解炉出口温度受生料下料量、三次风量及喂煤量的影响；再根据现场操作人员的经验，在分解炉正常工作时，三次风阀门开度通常保持恒定；再者，由于生料下料口距分解炉较远，依靠调节喂料改变分解炉温度具有较大延迟，因此，现场操作人员通常依靠调节喂煤量来改变分解炉温度；当喂煤量增加时，分解炉出口温度也相应升高，因此选取喂煤量作为模型的输入变量，建立分解炉出口温度单入单出数学模型；

步骤S102中：分解炉处于典型工况下的数学模型：

根据步骤S101中的数据分析，可以得到分解炉出口温度与喂煤量间的关系如表1所示：

表1 分解炉出口温度与喂煤量及生料喂料量之间的关系

根据现场人员工作经验，分解炉出口温度工作于840℃-860℃为其典型工况，即分解炉出口温度在840℃至860℃时，分解炉工作效率最高；因此，采用回归分析算法对分解炉出口温度在840℃-860℃区间内建模；

由于分解炉出口温度在一天内变化范围很大(820℃至880℃)，然而，所选取的建模区间温差在20℃以内，因此可以假设分解炉出口温度与喂煤量采样值的内在联系是线性的；即由控制系统的测量装置对分解炉出口气体温度及喂煤量进行采样，得到n组测量数据，其结构形式如下：

$\left\{\begin{array}{c}y\left(2\right)={\mathrm{\&alpha;}}_1+{\mathrm{\&alpha;}}_{2}y\left(1\right)+{\mathrm{\&alpha;}}_{3}y\left(0\right)+{\mathrm{\&alpha;}}_{4}u\left(1\right)\\ y\left(3\right)={\mathrm{\&alpha;}}_1+{\mathrm{\&alpha;}}_{2}y\left(2\right)+{\mathrm{\&alpha;}}_{3}y\left(1\right)+{\mathrm{\&alpha;}}_{4}u\left(2\right)\\ ......\\ y\left(n\right)={\mathrm{\&alpha;}}_1+{\mathrm{\&alpha;}}_{2}y(n-1)+{\mathrm{\&alpha;}}_{3}y(n-2)+{\mathrm{\&alpha;}}_{4}u(n-1)\\ y(n+1)={\mathrm{\&alpha;}}_1+{\mathrm{\&alpha;}}_{2}y\left(n\right)+{\mathrm{\&alpha;}}_{3}y(n-1)+{\mathrm{\&alpha;}}_{4}u\left(n\right)\end{array}\right.---\left(1\right)$

式中，y₀，y₂，…，y_n+1为t₀至t_n+1时刻的分解炉出口温度，u₁，u₂，…u_n为t₁至t_n时刻的喂煤量，α₁，α₂，α₃，α₄为待辨识的模型参数；

将(1)式改写成矩阵形式，得到

y＝Aα    (2)

式中：

y＝[y(2)，y(3)，…，y(n+1)]^T

α＝[α₁，α₂，α₃，α₄]^T

$A={\left[\begin{array}{cccc}1& y\left(1\right)& y\left(0\right)& u\left(1\right)\\ 1& y\left(1\right)& y\left(1\right)& u\left(2\right)\\ .& .& .& .\\ .& .& .& .\\ .& .& .& .\\ 1& y\left(n\right)& y(n-1)& u\left(n\right)\end{array}\right]}_{n\&times;4}$

根据最小二乘估计，可得到回归方程为：

$\hat{y}(n+1)={\mathrm{\&alpha;}}_1+{\mathrm{\&alpha;}}_{2}y\left(n\right)+{\mathrm{\&alpha;}}_{3}y(n-1)+{\mathrm{\&alpha;}}_{4}u\left(n\right)---\left(3\right)$

回归系数为：

α＝(A^TA)^-1A^Ty    (4)

由此可得，基于最小二乘学习算法的分解炉出口温度数学模型为：

$y(n+1)={\widehat{\mathrm{\&alpha;}}}_1+{\widehat{\mathrm{\&alpha;}}}_{2}y\left(n\right)+{\widehat{\mathrm{\&alpha;}}}_{3}y(n-1)+{\widehat{\mathrm{\&alpha;}}}_{4}u\left(n\right)---\left(5\right)$

核心程序编写如下：

clear all

clc                    ％清屏

％读取文件名为840-860的数据表中数据

caiyangzhi＝xlsread(′C：\Users\asus\Desktop\840-860.xls′)；

ceshishuju＝xlsread(′C：\Users\asus\Desktop\840-860.xls′)；

a(1：200，1)＝caiyangzhi(1：200，1)；

b(1：200，1)＝caiyangzhi(1：200，2)；

c(1：200，1)＝caiyangzhi(1：200，3)；

d(1：200，1)＝caiyangzhi(1：200，4)；

e(1：200，1)＝caiyangzhi(1：200，5)；

％求取各变量的最大值及最小值

minb＝min(b)；maxb＝max(b)；

minc＝min(c)；maxc＝max(c)；

mind＝min(d)；maxd＝max(d)；

mine＝min(e)；maxe＝max(e)；

％数据归一化处理

for i＝1：1：200

b(i，1)＝(b(i，1)-minb)/(maxb-minb)；

c(i，1)＝(c(i，1)-minc)/(maxc-minc)；

d(i，1)＝(d(i，1)-mind)/(maxd-mind)；

e(i，1)＝(e(i，1)-mine)/(maxe-mine)；

end

A1(1：200，1)＝a；

A1(1：200，2)＝b；

A1(1：200，3)＝c；

A1(1：200，4)＝d；

Y1(1：200，1)＝e；

[b1，bint1，r1，rint1，stats1]＝regress(Y1，A1，0.05)；

％取采样数据表中前100行数据作为测试数据

Q2＝caiyangzhi(1：100，5)；

Q1＝A1(1：100，1：4)*b1；

％数据反归一化

for k＝1：1：100

Q1(k，1)＝Q 1(k，1)*(maxe-mine)+mine；

End

％绘图

figure(1)

x＝1：100；

plot(x，Q 1，′r-o′，x，Q2，′b--+′)

figure(2)

x＝1：100；

plot(Q1-Q2，′r′)；

原始建模数据如下表2所示：

表2 原始建模数据

综上所述，分解炉出口温度的数学模型为：

$y(n+1)={\widehat{\mathrm{\&alpha;}}}_1+{\widehat{\mathrm{\&alpha;}}}_{2}y\left(n\right)+{\widehat{\mathrm{\&alpha;}}}_{3}y(n-1)+{\widehat{\mathrm{\&alpha;}}}_{4}u\left(n\right)---\left(6\right)$

式中，

$\begin{array}{cc}{\widehat{\mathrm{\&alpha;}}}_1=718.9233;& {\widehat{\mathrm{\&alpha;}}}_2=-0.0045;\\ {\widehat{\mathrm{\&alpha;}}}_3=-0.0002;& {\widehat{\mathrm{\&alpha;}}}_4=9.6468;\end{array}$

步骤S103中：分解炉出口温度的滑模控制器：

根据步骤S102中所建立的数学模型，由式(6)可得如下状态方程：

A(z^-1)y(k+1)＝B(z^-1)u(k)+θ    (7)

其中，

A(z^-1)＝1-0.0045z^-1-0.0002z^-2

B(z^-1)＝9.6468

θ＝718.9233

将式(7)改写成如下形式：

x(k+1)＝Ax(k)+Bu(k)+θ    (8)

y(k)＝Cx(k)

其中：

$A=\left[\begin{array}{cc}-0.0045& -0.0002\\ 1& 0\end{array}\right];$ $B=\left[\begin{array}{c}9.6468\\ 0\end{array}\right];$ C＝[1 0]； $\mathrm{\&theta;}=\left[\begin{array}{c}718.9233\\ 0\end{array}\right]$

切换函数选取为：

s(k)＝C_eE＝C_e(R(k)-x(k))    (9)

其中，C_e＝[c，1]；R＝[r(k)，dr(k)]，r(k)为位置指令；式(9)可被改写成如下形式：

s(k+1)＝C_e(R(k+1)-x(k+1))

＝C_e(R(k+1)-Ax(k)-Bu(k)-θ)    (10)

＝C_eR(k+1)-C_eAx(k)-C_eBu(k)-C_eθ

控制率为：

u(k)＝(C_eB)^-1(C_eR(k+1)-C_eAxk)-C_eθ-s(k+1))    (11)

基于指数的离散趋近率为：

s(k+1)＝s(k)+T(-εsgn(s(k))-qs(k))    (12)

其中，ε＞0；q＞0；1-qT＞0；T为采样时间，T<<1.0；根据(12)，(11)可改写成：

u(k)＝(C_eB)^-1(C_eR(k+1)-C_eAx(k)-C_eθ-s(k)-ds(k))    (13)

其中，ds(k)＝-εT sgn(s(k))-qTs(k)；

根据式(12)，可以得到：

$\begin{array}{c}s(k+1)=(1-\mathrm{qT})s\left(k\right)-\mathrm{\&epsiv;T}\frac{s\left(k\right)}{\left|s\left(k\right)\right|}\\ =(1-\mathrm{qT}-\frac{\mathrm{\&epsiv;T}}{\left|s\left(k\right)\right|})s\left(k\right)=\mathrm{ps}\left(k\right)\end{array}---\left(14\right)$

由式(14)有如下三种情况：

1、当 $\left|s\left(k\right)\right|>\frac{\mathrm{\&epsiv;T}}{2-\mathrm{qT}},$ 有：

$p>1-\mathrm{qT}-\frac{\mathrm{\&epsiv;T}(2-\mathrm{qT})}{\mathrm{\&epsiv;T}}=1-\mathrm{qT}-(2-\mathrm{qT})=-1$

因此，|p|＜1；|s(k+1)|＜|s(k)|；|s(k)|是递减的；

2、当 $\left|s\left(k\right)\right|<\frac{\mathrm{\&epsiv;T}}{2-\mathrm{qT}},$ 有：

$p<1-\mathrm{qT}-\frac{\mathrm{\&epsiv;T}(2-\mathrm{qT})}{\mathrm{\&epsiv;T}}=1-\mathrm{qT}-(2-\mathrm{qT})=-1$

因此，|p|＞1；|s(k+1)|＞|s(k)|；|s(k)|是递增的；

3、当 $\left|s\left(k\right)\right|=\frac{\mathrm{\&epsiv;T}}{2-\mathrm{qT}},$ 有：

$p=1-\mathrm{qT}-\frac{\mathrm{\&epsiv;T}(2-\mathrm{qT})}{\mathrm{\&epsiv;T}}=1-\mathrm{qT}-(2-\mathrm{qT})=-1$

因此，|p|＝1；|s(k+1)|＝|s(k)|；|s(k)|是震荡的；

显然，|s(k)|需要满足：

$\left|s\left(k\right)\right|>\frac{\mathrm{\&epsiv;T}}{2-\mathrm{qT}}---\left(15\right)$

式(15)可以改写成：

$\mathrm{\&epsiv;}<\frac{1}{T}(2-\mathrm{qT})\left|s\left(k\right)\right|---\left(16\right)$

取ε＝|s(k)|/2，式(12)可以改写成：

$s(k+1)-s\left(k\right)=-\mathrm{qTs}\left(k\right)-\frac{s\left(k\right)}{2}\mathrm{Tsgn}\left(s\left(k\right)\right)---\left(17\right)$

控制率为：

$\begin{array}{c}u\left(k\right)={\left(C_{e}B\right)}^{-1}(C_{e}R(k+1)-C_e\mathrm{Ax}\left(k\right)-C_e\mathrm{\&theta;}\\ -s\left(k\right)+\mathrm{qTs}\left(k\right)+\frac{s\left(k\right)}{2}\mathrm{Tsgn}\left(s\left(k\right))\right)\end{array}---\left(18\right)$

稳定性分析：

由式(14)可得：

$\begin{array}{c}(s(k+1)-s\left(k\right))\mathrm{sgn}\left(s\left(k\right)\right)\\ =(-\mathrm{qTs}\left(k\right)-\frac{\left|s\left(k\right)\right|}{2}\mathrm{Tsgn}\left(s\left(k\right)\right))\mathrm{sgn}\left(s\left(k\right)\right)\\ =-(q+0.5)T\left|s\left(k\right)\right|<0\end{array}$

$\begin{array}{c}(s(k+1)+s\left(k\right))\mathrm{sgn}\left(s\left(k\right)\right)\\ =((2-\mathrm{qT})s\left(k\right)-\frac{\left|s\left(k\right)\right|}{2}\mathrm{Tsgn}\left(s\left(k\right)\right))\mathrm{sgn}\left(s\left(k\right)\right)\\ =(2-0.5T-\mathrm{qT})\left|s\left(k\right)\right|>0\end{array}$

因此，滑模趋近式(17)满足滑动模态的存在性和到达条件，所设计的控制器是稳定的；

核心程序编写如下：

％清屏

clear all；

close all；

％设置采样时间及参数

ts＝0.01；

A＝[-0.0002，-0.0045；1，0]；

B＝[9.6468；0]；

C＝[1，0]；

D＝[718.9233；0]；

％设置初值

x＝[840；840]；

r_1＝855；r_2＝850；

c＝10；q＝0.1；

Ce＝[c，0]；

for k＝1：1：8000

time(k)＝k*ts；

％设置期望输出温度

r(k)＝855.0；

％外推法

dr(k)＝(r(k)-r_1)/ts；

dr_1＝(r_1-r_2)/ts；

r1(k)＝2*r(k)-r_1；

dr1(k)＝2*dr(k)-dr_1；

R＝[r(k)；dr(k)]；

R1＝[r1(k)；dr1(k)]；

E＝R-x；

e(k)＝E(1)；

de(k)＝E(2)；

s(k)＝Ce*E；

％自适应ε设定值

eq(k)＝abs(s(k))/2；

ds(k)＝-eq(k)*ts*sign(s(k))-q*ts*s(k)；

u(k)＝inv(Ce*B)*(Ce*R1-Ce*A*x-Ce*D-s(k)-ds(k))；

xx＝x(1)；

x＝A*x+B*(u(k))+D；

y(k)＝x(1)；

x(2)＝xx；

％参数更新

r_2＝r_1；

r_1＝r(k)；

end

％绘图

％控制曲线

figure(1)

plot(time，r，′r′，time，y，′b′，′linewidth′，2)；

xlabel(′Time(s)′)；ylabel(′Outlet temperature(℃)′)；

％控制误差

figure(2)

plot(time，r-y)；

xlabel(′Time(s)′)；ylabel(′Tracking error(℃)′)；

％控制量(喂煤量)

figure(3)

plot(time，u，′r′，′linewidth′，2)；

xlabel(′Time(s)′)；ylabel(′u′)；

％自适应率变化

figure(4)；

plot(time，eq，′r′，′linewidth′，2)；

xlabel(′time(s)′)；ylabel(′adaptive\epsilon′)；

下面介绍本发明的具体实施例；

以山东某水泥厂生产线数据为依据，建立基于回归分析的分解炉出口温度数学模型，另取100组数据进行测试，验证所建模型的可靠性；在此基础之上，建立相应的滑模控制器，设置控制器参数C_e＝[10，1]；q＝0.1；T＝0.01；

图4显示了分解炉出口温度实际值与所建模型输出的拟合值对比结果；图5为相应的建模误差；图6为分解炉出口温度控制曲线图；图7为分解炉出口温度跟踪误差曲线图；图8为控制量(喂煤)变化曲线图；图9为自适应率ε变化曲线图。

由图4-5可知，所建模型能够准确反映分解炉出口温度变化状态，误差在±1℃以内，满足工业生产要求；由图6-8可知，所设计的滑模控制器能够准确调节喂煤量以实现分解炉的出口温度的控制，调节时间在1分钟以内，稳态误差回零，控制准确；由图9可知，自适应率ε在控制初期较大，随着时间的增加，其值逐渐减小回零，既避免了系统的抖振，又提高了系统到达切换面的速度。。

以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。

Claims (6)

1.一种基于回归模型的分解炉出口温度滑模控制方法，其特征在于，该基于回归模型的分解炉出口温度滑模控制方法包括以下步骤：

步骤一，根据分解炉出口温度与喂煤量之间的关系的关系，采用回归分析算法，建立分解炉在典型工况下的数学模型；

分解炉数学模型为：

$y(n+1)={\widehat{\mathrm{\&alpha;}}}_1+{\widehat{\mathrm{\&alpha;}}}_{2}y\left(n\right)+{\widehat{\mathrm{\&alpha;}}}_{3}y(n-1)+{\widehat{\mathrm{\&alpha;}}}_{4}u\left(n\right)$

式中，y(n+1)为n+1分解炉出口温度，y(n)为n分解炉出口温度，y(n-1)为n-1分解炉出口温度，α₁，α₂，α₃，α₄为待辨识的模型参数，u(n)为n时刻喂煤量；

步骤二，根据分解炉数学模型，采用自适应趋近率求取最优控制量，建立分解炉的滑模控制器，从而实现分解炉出口温度控制；

自适应趋近率为：

$s(k+1)-s\left(k\right)=-\mathrm{qTs}\left(k\right)-\frac{s\left(k\right)}{2}\mathrm{Tsgn}\left(s\left(k\right)\right)$

式中，q＝0.1，T＝0.01，s(k)，s(k+1)为k及k+1时刻的切换函数：

s(k)＝C_eE＝C_e(R(k)-x(k))

其中，C_e＝[10，1]，R(k)＝[r(k)，dr(k)]，r(k)为位置指令，dr(k)为r(k)的微分，x(k)为分解炉出口温度初值；

控制率为：

$u\left(k\right)={\left(C_{e}B\right)}^{-1}(C_{e}R(k+1)-C_e\mathrm{Ax}\left(k\right)-C_e\mathrm{\&theta;}-s\left(k\right)+\mathrm{qTs}\left(k\right)+\frac{s\left(k\right)}{2}\mathrm{Tsgn}\left(s\left(k\right)\right))$

式中， $A=\left[\begin{array}{cc}-0.0045& -0.0002\\ 1& 0\end{array}\right],B=\left[\begin{array}{c}9.6468\\ 0\end{array}\right],\mathrm{\&theta;}=\left[\begin{array}{c}718.9233\\ 0\end{array}\right],$ R(k+1)为k+1时刻的位置指令，sgn(·)为符号函数。

2.如权利要求1所述的基于回归模型的分解炉出口温度滑模控制方法，其特征在于，在步骤一之前需要得出分解炉出口温度与喂煤量之间的关系，并选取分解炉出口温度840℃～860℃作为分解炉的工况，喂煤量变化范围为13.06t～15.1t。

3.如权利要求1所述的基于回归模型的分解炉出口温度滑模控制方法，其特征在于，分解炉出口温度与喂煤量之间的关系：当分解炉处于正常工作状态时，分解炉出口温度随喂煤量的增加而升高。

4.如权利要求1所述的基于回归模型的分解炉出口温度滑模控制方法，其特征在于，在步骤一中分解炉出口温度数学模型的输入变量选取的工艺流程如下：

步骤一，生料由提升机喂入C1～C2级旋风筒的的连接管道；

步骤二，生料由来自C2级旋风筒的热风带入C1级旋风筒进行气固热交换，再由C1级旋风筒底部的锁风阀排出，进入C2～C3级旋风筒的连接管道，又被气流带入C2级旋风筒内继续气固热交换，如此反复；

步骤三，预热后的物料经由C4级旋风筒锥部进入分解炉，煤粉则从分解炉中部给煤口进入分解炉；

步骤四，从分解炉出来的物料经分解炉上部的鹅颈管进入C5级旋风筒；最后，由C5级旋风筒的锥部进入回转窑进行煅烧。

5.如权利要求1所述的基于回归模型的分解炉出口温度滑模控制方法，其特征在于，分解炉数学模型获取的方法包括：

根据分解炉出口温度工作于840℃-860℃为典型工况，即分解炉出口温度在840℃至860℃时，分解炉工作效率最高，采用回归分析算法对分解炉出口温度在840℃-860℃区间内建模；

分解炉出口温度与喂煤量采样值的内在联系是线性的，即由控制系统的测量装置对分解炉出口气体温度及喂煤量进行采样，得到n组测量数据，结构形式如下：

$\left\{\begin{array}{c}y\left(2\right)={\mathrm{\&alpha;}}_1+{\mathrm{\&alpha;}}_{2}y\left(1\right)+{\mathrm{\&alpha;}}_{3}y\left(0\right)+{\mathrm{\&alpha;}}_{4}u\left(1\right)\\ y\left(3\right)={\mathrm{\&alpha;}}_1+{\mathrm{\&alpha;}}_{2}y\left(2\right)+{\mathrm{\&alpha;}}_{3}y\left(1\right)+{\mathrm{\&alpha;}}_{4}u\left(2\right)\\ ......\\ y\left(n\right)={\mathrm{\&alpha;}}_1+{\mathrm{\&alpha;}}_{2}y(n-1)+{\mathrm{\&alpha;}}_{3}y(n-2)+{\mathrm{\&alpha;}}_{4}u(n-1)\\ y(n+1)={\mathrm{\&alpha;}}_1+{\mathrm{\&alpha;}}_{2}y\left(n\right)+{\mathrm{\&alpha;}}_{3}y(n-1)+{\mathrm{\&alpha;}}_{4}u\left(n\right)\end{array}\right.---\left(1\right)$

式中，y₀，y₂，…，y_n+1为t₀至t_n+1时刻的分解炉出口温度，u₁，u₂，…u_n为t₁至t_n时刻的喂煤量，α₁，α₂，α₃，α₄为待辨识的模型参数：

将(1)式改写成矩阵形式，得到：

y＝Aα   (2)

式中：

y＝[y(2)，y(3)，…，y(n+1)]^T

α＝[α₁，α₂，α₃，α₄]^T

$A={\left[\begin{array}{cccc}1& y\left(1\right)& y\left(0\right)& u\left(1\right)\\ 1& y\left(1\right)& y\left(1\right)& u\left(2\right)\\ .& .& .& .\\ .& .& .& .\\ .& .& .& .\\ 1& y\left(n\right)& y(n-1)& u\left(n\right)\end{array}\right]}_{n\&times;4}$

根据最小二乘估计，得到回归方程为：

$\hat{y}(n+1)={\mathrm{\&alpha;}}_1+{\mathrm{\&alpha;}}_{2}y\left(n\right)+{\mathrm{\&alpha;}}_{3}y(n-1)+{\mathrm{\&alpha;}}_{4}u\left(n\right)---\left(3\right)$

回归系数为：

α＝(A^TA)^-1A^Ty   (4)

由此得，基于最小二乘学习算法的分解炉出口温度数学模型为：

$y(n+1)={\widehat{\mathrm{\&alpha;}}}_1+{\widehat{\mathrm{\&alpha;}}}_{2}y\left(n\right)={\widehat{\mathrm{\&alpha;}}}_{3}y(n-1)+{\widehat{\mathrm{\&alpha;}}}_{4}u\left(n\right).$

6.如权利要求1所述的基于回归模型的分解炉出口温度滑模控制方法，其特征在于，自适应趋近率和控制率公式的获取方法：

根据所建立的数学模型，由式得如下状态方程：

A(z^-1)y(k+1)＝B(z^-1)u(k)+θ   (7)

其中：

A(z^-1)＝1-0.0045z^-1-0.0002z^-2

B(z^-1)＝9.6468

θ＝718.9233

将式(7)改写成如下形式：

x(k+1)＝Ax(k)+Bu(k)+θ   (8)

y(k)＝Cx(k)

其中：

$A=\left[\begin{array}{cc}-0.0045& -0.0002\\ 1& 0\end{array}\right];B=\left[\begin{array}{c}9.6468\\ 0\end{array}\right];C=\left[\begin{array}{cc}1& 0\end{array}\right];\mathrm{\&theta;}=\left[\begin{array}{c}718.9233\\ 0\end{array}\right]$

切换函数选取为：

s(k)＝C_eE＝C_e(R(k)-x(k))   (9)

其中，C_e＝[c，1]；R＝[r(k)，dr(k)]，r(k)为位置指令，式(9)被改写成如下形式：

s(k+1)＝C_e(R(k+1)-x(k+1))

＝C_e(R(k+1)-Ax(k)-Bu(k)-θ)   (10)

＝C_eR(k+1)-C_eAx(k)-C_eBu(k)-C_eθ

控制率为：

u(k)＝(C_eB)^-1(C_eR(k+1)-C_eAx(k)-C_eθ-s(k+1))   (11)

基于指数的离散趋近率为：

s(k+1)＝s(k)+T(-εsgn(s(k))-qs(k))   (12)

其中，ε＞0；q＞0；1-qT＞0；T为采样时间，T＜＜1.0，根据(12)，式(11)改写成：

u(k)＝(C_eB)^-1(C_eR(k+1)-C_eAx(k)-C_eθ-s(k)-ds(k))   (13)

其中，ds(k)＝-εTsgn(s(k))-qTs(k)；

根据式(12)，得到：

$\begin{array}{c}s(k+1)=(1-\mathrm{qT})s\left(k\right)-\mathrm{\&epsiv;T}\frac{s\left(k\right)}{\left|s\left(k\right)\right|}\\ =(1-\mathrm{qT}-\frac{\mathrm{\&epsiv;T}}{\left|s\left(k\right)\right|})s\left(k\right)=\mathrm{ps}\left(k\right)\end{array}---\left(14\right)$ 由式(14)有如下三种情况：

当 $\left|s\left(k\right)\right|>\frac{\mathrm{\&epsiv;T}}{2-\mathrm{qT}},$ 有：

$p>1-\mathrm{qT}-\frac{\mathrm{\&epsiv;T}(2-\mathrm{qT})}{\mathrm{\&epsiv;T}}=1-\mathrm{qT}-(2-\mathrm{qT})=-1$

因此，|p|＜1；|s(k+1)|＜|s(k)|；|s(k)|是递减的；

当 $\left|s\left(k\right)\right|<\frac{\mathrm{\&epsiv;T}}{2-\mathrm{qT}},$ 有：

$p<1-\mathrm{qT}-\frac{\mathrm{\&epsiv;T}(2-\mathrm{qT})}{\mathrm{\&epsiv;T}}=1-\mathrm{qT}-(2-\mathrm{qT})=-1$

因此，|p|＞1；|s(k+1)|＞|s(k)|；|s(k)|是递增的；

当 $\left|s\left(k\right)\right|=\frac{\mathrm{\&epsiv;T}}{2-\mathrm{qT}},$ 有：

$p=1-\mathrm{qT}-\frac{\mathrm{\&epsiv;T}(2-\mathrm{qT})}{\mathrm{\&epsiv;T}}=1-\mathrm{qT}-(2-\mathrm{qT})=-1$

因此，|p|＝1；|s(k+1)|＝|s(k)|；|s(k)|是震荡的；

显然，|s(k)|需要满足：

$\left|s\left(k\right)\right|>\frac{\mathrm{\&epsiv;T}}{2-\mathrm{qT}}---\left(15\right)$

式(15)改写成：

$\mathrm{\&epsiv;}<\frac{1}{T}(2-\mathrm{qT})\left|s\left(k\right)\right|---\left(16\right)$

取ε＝|s(k)|/2，式(12)改写成：

自适应趋近率： $s(k+1)-s\left(k\right)=-\mathrm{qTs}\left(k\right)-\frac{s\left(k\right)}{2}\mathrm{Tsgn}\left(s\left(k\right)\right)$

控制率为：

$\begin{array}{c}u\left(k\right)={\left(C_{e}B\right)}^{-1}(C_{e}R(k+1)-C_e\mathrm{Ax}\left(k\right)-C_e\mathrm{\&theta;}\\ -s\left(k\right)+\mathrm{qTs}\left(k\right)+\frac{s\left(k\right)}{2}\mathrm{Tsgn}\left(s\left(k\right)\right))\end{array}.$