CN102629104B - 一种水泥回转窑煅烧预测控制系统及方法 - Google Patents

Abstract

一种水泥回转窑煅烧预测控制方法及系统，其方法是：(1)采集水泥回转窑煅烧过程的现场数据并对其分类；(2)对其数据分别进行模型辨识，并将其有机结合，建立预测模型；(3)用预测模型对煅烧过程的历史和未来的数据信息预测出煅烧过程的输出，并利用模型输出误差反馈校正得到闭环的煅烧过程的预测输出；(4)依据闭环预测输出和参考输出轨迹，构造出非线性目标函数，用序列二次规划法对目标函数求最优解，得到煅烧控制量的预测值。系统包括与水泥回转窑连接的智能检测仪表和执行器、数据存储装置及上位机，其中上位机嵌入了煅烧预测控制算法。本发明能够适应水泥回转窑煅烧过程的动态性，多变量之间耦合性、非线性、时滞性，得到良好的控制效果。

Description

一种水泥回转窑煅烧预测控制系统及方法

技术领域

本发明涉及水泥生产过程的先进控制领域，特别是涉及一种水泥回转窑煅烧过程预测控制系统和方法。

背景技术

水泥回转窑煅烧是水泥生产过程中一个很典型很重要的环节，对它的操作质量直接影响着水泥的品质、产能和生产成本。水泥回转窑煅烧操作的好坏很大程度上又决定于对它的控制方法，所以水泥回转窑煅烧的先进控制方法也就越来越受到国内外的关注和重视。

水泥回转窑煅烧之所以难以控制，主要在于其具有复杂的变量耦合性、非线性和时滞性，加之其最重要的被控变量——窑内温度，无法准确选取其检测点也是难点之一，因为传统的经验调节法、PID调节法、模糊控制法等控制方法都很难得到较好且稳定的控制效果。由于国家对经济结构进行调整，对节能减排提出了较高的要求，这就使得水泥回转窑煅烧过程控制的研究显得越来越重要。

许多国内外的工艺、自控专家对此做了大量的研究工作，提出了许多先进的控制理论，也取得了一些令人鼓舞的进展。但是鉴于水泥回转窑煅烧过程的复杂性和先进控制理论对控制对象的局限性，使得一些先进的控制算法和控制方案，未能在水泥回转窑煅烧过程控制中取得实质的突破性进展。

发明内容

为了克服已有的水泥回转窑控制方案的不能适应煅烧过程变量之间强耦合性、强烈的非线性、变量之间的时滞性和不能得到良好控制效果的不足，本发明提供一种能够解决煅烧过程的变量耦合性、大时滞性、系统非线性问题，并能得到良好控制效果的基于组合模型的水泥回转窑煅烧过程预测控制方法及系统。

本发明解决其技术问题所采用的技术方案是：

一种基于BP-ARX时滞模型的水泥回转窑煅烧过程预测控制方法，该预测控制方法包括以下步骤：

步骤一：确定煅烧过程的主要输入输出量并获得分类辨识数据，影响煅烧过程的主要输入量窑头喷煤量和高温风机挡板开度，输出量窑内NO_x含量和氧气含量，采集相应变量的现场运行数据存放于数据存储装置，并将其分类为稳态数据和动态数据；

步骤二：建立预测模型，利用BP神经网络对分类出的稳态数据进行辨识得到水泥回转窑煅烧过程的非线性BP稳态模型，利用最小二乘法对分类出的动态数据进行辨识得到煅烧过程的线性时滞ARX动态模型，并通过将稳态模型和动态模型进行有机并联组合，得到水泥回转窑煅烧过程的MIMO时滞系统BP-ARX预测模型；

步骤三：预测未来输出状态，在当前时刻k，利用水泥回转窑煅烧过程的历史和未来数据的输入输出信息，通过基于BP-ARX模型的预测控制模型，并将其进行迭代对水泥回转窑煅烧过程未来一段时间的输出状态窑内NO_x含量y_1n(k+j)、氧气含量y_2n(k+j)进行预测；

步骤四：误差反馈校正，将从步骤三中得到的回转窑煅烧过程未来的输出状态y_in(k+j)和k时刻的输出误差e_i(k)相加，得到基于k时刻回转窑煅烧过程未来的闭环预测输出窑内NO_x含量y_1c(k+j)、氧气含量y_2c(k+j)；

步骤五：设定输出量参考轨迹，为使煅烧过程未来输出量能够沿着预先设定轨迹平稳地达到设定值，引入参考输出轨迹窑内NO_x含量y_1r(k+j)、氧气含量y_2r(k+j)；

步骤六：滚动优化求解，将得到的预测闭环输出y_ic(k+j)与设定的参考轨迹y_ir(k+j)进行比较，构建带约束的二次型目标函数，并通过序列二次规划法对其进行滚动优化求解，计算出当前时刻应加于水泥回转窑的喷煤量u₁(k)、高温风机挡板开度u₂(k)；

步骤七：根据步骤六计算出的当前时刻输入量控制现场智能执行器，实现现场水泥回转窑煅烧过程的自动控制。

上述步骤一中所述确定影响水泥回转窑煅烧过程的主要输入输出量，输入量为窑头喷煤量和高温风机挡板开度，输出量为窑内NO_x含量和窑内氧气含量，选择的依据为：水泥回转窑煅烧的质量主要取决于对回转窑窑内温度的准确控制和窑内煤的燃烧状况，由于对窑内温度不易直接测量，但窑内温度的高低和窑内NO_x的含量存在一种固定关系，温度越高NO_x含量越高，并且气体在窑内的流动速度很快，用窑尾NO_x的含量可以快速地反映出窑内温度，同时窑尾O₂含量可以反映出窑内煤的燃烧情况，既保证了煅烧过程温度的平稳，也兼顾了单产能耗。

上述步骤一中所述的数据分类是将采集的现场数据划分为稳态数据和动态数据，划分依据是：当回转窑煅烧过程的输入量喷煤量和高温风机挡板开度中任意一个发生阶跃变化时，输出量NO_x含量和O₂含量从阶跃时刻到最终都达到稳定状态时的数据为动态数据；输出量从稳定状态到下一个阶跃发生前的数据为稳态数据。

上述步骤二所述的水泥回转窑煅烧过程的MIMO时滞系统预测模型是稳态模型和含有时滞的动态模型并联而成，稳态模型是通过模糊神经网络根据回转窑煅烧过程的稳态数据辨识获得，非线性时滞动态模型是通过最小二乘法对煅烧过程的动态数据辨识获得。

上述步骤三所述的预测模型迭代求解水泥回转窑煅烧过程未来一段时间内的输出预测值，采用矩阵分离变换的方法来完成对MIMO时滞系统未来预测输出量状态y_in(k+j)的求解。

上述步骤六所述预测控制器在滚动时域下采用序列二次规划法对带有约束的二次型目标函数进行求解，保证求出的煅烧过程预测输入解是在约束条件下的全局最优解。

本发明还公开了基于BP-ARX时滞模型的水泥回转窑煅烧过程预测控制系统，该系统包括现场智能仪表、智能执行器、数据存储装置及上位机，智能仪表和智能执行器直接与现场相连，智能仪表和智能执行器又通过现场总线与DCS系统、上位机依次相连。

智能仪表用于采集水泥回转窑煅烧过程的输出量，即窑内NO_x含量和氧气含量并将其传送到数据存储装置；

数据存储装置为DCS系统的数据存储装置，所述的DCS系统包括数据接口、控制站和数据存储装置，完成对智能仪表上传信息的采集和存储，并将上位机的命令下达至智能执行器；

上位机用于运行回转窑煅烧过程预测控制算法，根据水泥回转窑煅烧过程的输出量计算出当前时刻应加入回转窑内煤和风的控制量，并通过DCS系统对智能执行器进行调节来实现现场回转窑煅烧过程的自动控制。所述的回转窑煅烧过程预测控制算法，是指上述方法中所说的步骤一至六。

本发明的有益效果主要表现在：

1．选择窑内NO_x含量和O₂含量作为水泥回转窑煅烧过程的控制输出量，既保证了煅烧过程温度的平稳，也兼顾了单产能耗；

2．采用稳态和动态分离辨识，再有机并联组合的模型辨识方法，并采用矩阵分离变换的方法将时滞添加到MIMO系统的预测模型中，能够适应水泥回转窑煅烧过程的变量耦合性、大时滞性、系统非线性，实现了对水泥回转窑煅烧过程的平稳控制；

3．充分利用预测控制技术的优点，引入参考轨迹，反馈校正和滚动优化技术，获取更多的水泥回转窑煅烧过程的运行信息，实现了回转窑煅烧过程的自动控制，取得良好控制效果；

4．利用序列二次规划法对带有约束的二次型目标函数进行滚动优化求解，易得的全局最优解；操作简单，适应性强。

附图说明

图1为基于BP-ARX时滞模型的水泥回转窑煅烧过程预测控制方法流程图。

图2为本发明的BP-ARX并联时滞模型框图。

图3为本发明提出的水泥回转窑煅烧过程预测控制系统的方框图。

图4为本发明提出的水泥回转窑煅烧过程预测控制系统的现场接线图。.

具体实施例

下面结合附图对本发明做进一步的详细描述。

如图1所示的预测控制方法的流程图和图3所示控制系统方框图：

1．建立基于BP-ARX模型的预测模型

BP-ARX预测模型是一个由非线性BP稳态模型和线性ARX动态模型并联而成的组合模型，详细的模型结构关系如图2所示。建立该BP-ARX模型大体分两步：训练数据获取和模型辨识。

（1）．训练数据获取

本实施例对水泥厂的回转窑煅烧过程中NO_x含量、O₂含量、喷煤量以及高温风机挡板开度四个变量通过图3中DCS系统的数据存储装置（11）对现场数据进行采集，采样周期为60s，共采集20000组数据。将所采集数据分成动态、静态两种数据，划分依据是：当输入量喷煤量u₁和高温风机挡板开度u₂中任意一个发生阶跃变化时，输出量NO_x含量y₁和O₂含量y₂从阶跃时刻到最终都达到稳定状态时的数据选为动态数据；输出量从稳定状态到下一个阶跃发生前的数据选为稳态数据。经分类后得到动态数据7341组，稳态数据12659组。

（2）．模型辨识

参考图2中模型结构，分为三部分：BP稳态模型（6）、动态增益K（7）和ARX动态模型（8），用BP稳态模型就出当前时刻系统的动态增益K，再利用动态增益K实时调整ARX动态模型参数，完成稳态模型和动态模型的有机并联组合，建立代表水泥回转窑煅烧过程综合特性的BP-ARX模型。

非线性BP稳态模型辨识

BP神经网络选取三层网络结构，其中，输入层变量的维数为n_u=2，隐含层的输出变量维数为p=10，即隐含层有10个神经元，输出层的输出变量维数为n_y=2,W_1,1表示输入层到隐含层的连接权值，W_2,1表示隐含层到输出层的连接权值，θ₁、θ₂分别表示隐含层和输出层神经元的阈值。BP算法的具体训练步骤如下：

1）随机初始化权值矩阵W_1,1，W_2,1以及神经元阈值向量θ₁，θ₂；

2）由训练数据中第l个样本数据输入值u_l计算出对应的隐含层的输出值，如式（1）所示。

$x_l\left(i\right)=f_1(\underset{k}{\overset{\mathrm{nu}}{\mathrm{\Σ}}}{W}_{\mathrm{1,1}}(i,k)u_l\left(k\right)+{\mathrm{\θ}}_1\left(i\right))---\left(1\right)$

式中：x_l(i)表示隐含层中第i个单元的输出，其中i=1,2,…,p，k=1,2,…,n_u，W_1,1(i,k)表示权值矩阵W_1,1中位于i行k列的数值，即从第k个输入神经元到第i个输出单元的权值，f₁为隐含层的S型正切传递函数，函数表达式

其中 $\hat{x}=\underset{k}{\overset{\mathrm{nu}}{\mathrm{\Σ}}}{W}_{\mathrm{1,1}}(i,k)u_l\left(k\right)+{\mathrm{\θ}}_1\left(i\right);$

3）由2）步中得到隐含层输出值x_l计算对应的输出层的输出值:

$y_l=f_2(\underset{i=1}{\overset{p}{\mathrm{\Σ}}}{W}_{\mathrm{2,1}}(j,i)x_l\left(i\right)+{\mathrm{\θ}}_2\left(j\right))---\left(2\right)$

式中：j=1,2,…,n_y，输出层的传递函数f₂为一个线性函数，即：

其中C为常数， $\hat{x}=\underset{i=1}{\overset{p}{\mathrm{\Σ}}}{W}_{\mathrm{2,1}}(j,i)x_l\left(i\right)+{\mathrm{\θ}}_2\left(j\right);$

4）计算输出层误差：

$e_l\left(j\right)=y_l\left(j\right)[1-y_l\left(j\right)][\hat{y}\left(j\right)-y\left(j\right)]---\left(3\right)$

其中

为输出层中第j个单元的期望输出值。

5）计算隐含层误差：

$e_{x_l}\left(i\right)=x_l\left(i\right)[1-x_l\left(i\right)]\underset{j=1}{\overset{n_y}{\mathrm{\Σ}}}{W}_{\mathrm{2,1}}(j,i)e\left(j\right)---\left(4\right)$

6）循环调整隐含层到输出层的权值矩阵W_2,1中的每个单权值W_2,1(j,i)：

W_2,1(j,i)=W_2,1(j,i)+αx_l(i)e_l(j)  （5）

式中：α为学习率，0<α<1；

7）修正输出层的阈值θ₂(j)：

θ₂(j)=θ₂(j)+αe_l(j)  （6）

8）循环调整输入层到隐含层的权值矩阵W_1,1中每个单权值W_1,1(i,k)：

$W_{\mathrm{1,1}}(i,k)=W_{\mathrm{1,1}}(i,k)+{\mathrm{\&beta;u}}_l\left(j\right)e_{x_l}\left(i\right)---\left(7\right)$

式中：β为学习率，0<β<1；

9）修正隐含层的阈值θ₁(i)：

${\mathrm{\&theta;}}_1\left(i\right)={\mathrm{\&theta;}}_1\left(i\right)+{\mathrm{\&beta;e}}_{x_l}---\left(8\right)$

10）重复循环计算2）～9）步，直至对于所有的输入输出样本数据j=1,2,…,n_y，l=1,2,…,N的训练误差为足够小或者为零时停止训练，其中N为样本稳态数据总数，本实施例中N=12659。

综上所述，可以利用上述的BP算法对给定的输入、输出样本数据进行训练，最终训练出如下反映输入到输出关系的神经网络稳态模型：

y_s=f₂[W_2,1·f₁(W_1,1·u_s+θ₁)+θ₂]  （9）

式中： $u_s={\left[\begin{array}{cccc}{u}_{1s}& u_{2s}& \&CenterDot;\&CenterDot;\&CenterDot;& u_{n_{u}s}\end{array}\right]}^T$ 是系统的稳态输入， $y_s={\left[\begin{array}{cccc}{y}_{1s}& y_{2s}& \&CenterDot;\&CenterDot;\&CenterDot;& y_{n_{y}s}\end{array}\right]}^T$ 是系统的稳态输出，f₁和f₂分别是隐含层和输出层的传递函数，θ₁和θ₂分别是隐含层和输出层神经元的阈值，W_1,1和W_2,1分别是输入层到隐含层和隐含层到输出层的权值。

上述的BP稳态模型是通过离线训练获得，在工程实际应用中可以根据系统的实际运行数据每隔一段时间对稳态输入输出样本数据进行更新，重新训练该BP模型，以使该稳态模型尽可能逼近现场实际系统的稳态模型。

线性ARX时滞动态模型辨识

线性动态模型辨识是利用最小二乘法的辨识方法。首先对筛选出来动态数据进行处理，分别用当前时刻的输入输出值减去该时刻的输入输出稳态值获得动态变化增量数据，再对变化增量数据进行模型辨识。辨识建立水泥回转窑煅烧过程的线性动态ARX模型。如式（10）所示：

Δy_λ(k)=A₁ ^(λ)Δy(k-1)+A₂ ^(λ)Δy(k-2)+B₁ ^(λ)Δu(k-τ^(λ)-1)+B₂ ^(λ)Δu(k-τ^(λ)-2)  （10）

式中，Δy(k)=y(k)-y_s(k)，Δu(k)=u(k)-u_s(k)，Δy(k)=[Δy₁(k)  Δy₂(k)  …  Δy_ny(k)]^T,Δu(k)=[Δu₁(k)  Δu₂(k)  …  Δu_nu(k)]^T，其中u(k-τ^(λ))=[u₁(k-τ_λ1)  u₂(k-τ_λ2)  …  u_nu(k-τ_λnu)]^Ty(k)=[y₁(k)  y₂(k)  …  y_ny(k)]^T、分别是系统的输入向量和输出向量，u_s(k)、y_s(k)为系统的稳态输入输出值，结构与系统输入输出向量相同；其中，τ^(λ)表示时滞矩阵τ中的第λ行，τ_λβ表示时滞矩阵τ中的第λ行第β列的数值，

Δy_λ(k)分别表示取矩阵A_i、B_i和Δy(k)的第λ行，i=1,2；

ny为被控量个数ny=2，nu为控制量个数nu=2。

BP-ARX预测模型建立

BP-ARX模型是先将系统的稳态和动态特性进行单独辨识，最后再通过增益将BP稳态模型与动态ARX模型进行有机地结合，如图2所示，形成一个并联的组合模型。最终用BP-ARX模型来完全体现出系统的特性。稳态模型和动态模型具体结合过程如下，

在式（10）ARX模型中令

A₀=I-A₁-A₂  （11）

式中：I为n_y×n_y维的单位矩阵。则ARX模型的稳态增益矩阵K

为：

$K=\frac{B_1+B_2}{I-A_1-A_2}={(I-A_1-A_2)}^{-1}\&CenterDot;(B_1+B_2)=A_0^{-1}\&CenterDot;(B_1+B_2)---\left(12\right)$

式中：A₀按式(11)进行计算，A₁、A₂、B₁、B₂是ARX模型(10)中的参数。

因ARX模型式（10）描述的是系统输入增量与输出增量之间的变化关系，所以其稳态增益式（12）既是描述系统输入与输出之间变化关系的BP稳态模型式（12）的动态增益。

保持模型式（10）中的时滞τ和A₁、A₂不变，令

$B_i^{\&prime;}=A_0\&CenterDot;\frac{B_i}{B_1+B_2}=(I-A_1-A_2)\&CenterDot;[B_i./(B_1+B_2)]---\left(13\right)$

式中：B_i./(B₁+B₂)表示矩阵B_i与矩阵(B₁+B₂)中对应位置的元素相除，i=1,2。

将式（13）代入ARX模型式（10），则可得到下述ARX模型，该模型的稳态增益矩阵中各个元素的值都将恒等于1。

$\mathrm{\&Delta;}{y}_{\mathrm{\&lambda;}}\left(k\right)=A_1^{\left(\mathrm{\&lambda;}\right)}\mathrm{\&Delta;y}(k-1)+A_2^{\left(\mathrm{\&lambda;}\right)}\mathrm{\&Delta;y}(k-2)$   （14）

$+B_1^{\&prime;\left(\mathrm{\&lambda;}\right)}\mathrm{\&Delta;u}(k-{\mathrm{\&tau;}}^{\left(\mathrm{\&lambda;}\right)}-1)+B_2^{\&prime;\left(\mathrm{\&lambda;}\right)}\mathrm{\&Delta;u}(k-{\mathrm{\&tau;}}^{\left(\mathrm{\&lambda;}\right)}-2)$

式中：B′₁和B′₂按式（13）所得取值，其它符号的取值与模型式（10）相同。

再将动态ARX模型变换为含有稳态增益K参数的BP-ARX模型，如式（15）所示：

$\mathrm{\&Delta;}{y}_{\mathrm{\&lambda;}}\left(k\right)=A_1^{\left(\mathrm{\&lambda;}\right)}\mathrm{\&Delta;y}(k-1)+A_2^{\left(\mathrm{\&lambda;}\right)}\mathrm{\&Delta;y}(k-2)+B_{11}^{\left(\mathrm{\&lambda;}\right)}\mathrm{\&Delta;u}(k-{\mathrm{\&tau;}}^{\left(\mathrm{\&lambda;}\right)}-1)+B_{12}^{\left(\mathrm{\&lambda;}\right)}\mathrm{\&Delta;u}(k-{\mathrm{\&tau;}}^{\left(\mathrm{\&lambda;}\right)}-2)$   （15）

$+B_{21}^{\left(\mathrm{\&lambda;}\right)}\mathrm{\&Delta;u}(k-{\mathrm{\&tau;}}^{\left(\mathrm{\&lambda;}\right)}-1)+B_{22}^{\left(\mathrm{\&lambda;}\right)}\mathrm{\&Delta;u}(k-{\mathrm{\&tau;}}^{\left(\mathrm{\&lambda;}\right)}-2)$

式中，

B₁₁=A₀·[B′₁./(B′₁+B′₂).*K^c]，B₁₂=A₀·[B′₂./(B′₁+B′₂).*K^c]

B₂₁=A₀·[B′₁./(B′₁+B′₂).*K′]，B₂₂=A₀·[B′₂./(B′₁+B′₂).*K′]

K′=(Kⁿ-K^c)./(u_s(k+1)-u_s(k))

${\mathrm{\&Delta;u}}^2(k-{\mathrm{\&tau;}}^{\left(\mathrm{\&lambda;}\right)})={\left[\begin{array}{cccc}{\mathrm{\&Delta;u}}_1^2(k-{\mathrm{\&tau;}}_{\mathrm{\&lambda;}1})& {\mathrm{\&Delta;u}}_2^2(k-{\mathrm{\&tau;}}_{\mathrm{\&lambda;}2})& \&CenterDot;\&CenterDot;\&CenterDot;& {\mathrm{\&Delta;u}}_{\mathrm{nu}}^2(k-{\mathrm{\&tau;}}_{\mathrm{\&lambda;nu}})\end{array}\right]}^T,$ B′₁和B′₂是由式（13）中系数变换而得，K^c、Kⁿ分别是稳态输入向量u_s(k)、u_s(k+1)由BP稳态模型求出的对应动态增益矩阵，其它符号的取值与式(2)的取值相同，B₁*B₂表示B₁和B₂中对应位置的元素相乘，K./u_s表示K中每一行的元素和向量u_s中对应位置的元素相除。

上式（15）是基于k-1时刻以前的煅烧过程中输入输出量状态预测输出k时刻回转窑煅烧过程的输出增量状态，如何得到在预测时域P和控制时域M内的煅烧过程输出状态是BP-ARX预测模型的关键。本发明采用矩阵分离变换的方法来获得预测时域内的煅烧过程预测输出状态，首先构造出一组递推式。

S(n)=A₁·S(n-1)+A₂S(n-2)

T(n)=S(n)·B₁₁+S(n-1)·B₁₂

T_g(n)=S(n)·B₂₁+S(n-1)·B₂₂

式中：n=3,4,…,P+1，S为n_y×n_y维矩阵，T、T_g为n_y×n_u维矩阵，亦即

S(1)=I，S(2)=A₁，T(1)=B₁₁，T_g(1)=B₂₁，I为n_y×n_y维的单位矩阵，其它参数与式（15）中相同。

然后，依据所有输入量对第一个输出量的时延，即喷煤量u₁和高温风机挡板开度u₂对窑内NO_x含量y₁的延时，得到y₁在预测时域P内的预测输出值

再依据喷煤量u₁和高温风机挡板开度u₂对窑内O₂含量y₂的延时，得到y₂在预测时域P内的预测输出值

具体表达式如下。

当j<M+d_λ(j=1,2,…,P)时：

$\mathrm{\&Delta;}{\hat{y}}_{\mathrm{\&lambda;}}(k+j|k)={\left[S(j+1)\right]}^{\left(\mathrm{\&lambda;}\right)}\mathrm{\&Delta;y}\left(k\right)+{[S\left(j\right)\&CenterDot;A_2]}^{\left(\mathrm{\&lambda;}\right)}\mathrm{\&Delta;y}(k-1)$

$+{[S\left(j\right)\&CenterDot;B_{12}]}^{\left(\mathrm{\&lambda;}\right)}\mathrm{\&Delta;u}(k-{\mathrm{\&tau;}}^{\left(\mathrm{\&lambda;}\right)}-1)+\underset{i=1}{\overset{j}{\mathrm{\&Sigma;}}}{\left[T\left(i\right)\right]}^{\left(\mathrm{\&lambda;}\right)}\mathrm{\&Delta;u}(k+j-i-{\mathrm{\&tau;}}^{\left(\mathrm{\&lambda;}\right)})$

$+{[S\left(j\right)\&CenterDot;B_{22}]}^{\left(\mathrm{\&lambda;}\right)}{\mathrm{\&Delta;u}}^2(k-{\mathrm{\&tau;}}^{\left(\mathrm{\&lambda;}\right)}-1)+\underset{i=1}{\overset{j}{\mathrm{\&Sigma;}}}{\left[T_g\left(i\right)\right]}^{\left(\mathrm{\&lambda;}\right)}{\mathrm{\&Delta;u}}^2(k+j-i-{\mathrm{\&tau;}}^{\left(\mathrm{\&lambda;}\right)})$

当j≥M+d_λ(j=1,2,…,P)时：

$\mathrm{\&Delta;}{\hat{y}}_{\mathrm{\&lambda;}}(k+j|k)={\left[S(j+1)\right]}^{\left(\mathrm{\&lambda;}\right)}\mathrm{\&Delta;y}\left(k\right)+{[S\left(j\right)\&CenterDot;A_2]}^{\left(\mathrm{\&lambda;}\right)}\mathrm{\&Delta;y}(k-1)$

$+{[S\left(j\right)\&CenterDot;B_{12}]}^{\left(\mathrm{\&lambda;}\right)}\mathrm{\&Delta;u}(k-{\mathrm{\&tau;}}^{\left(\mathrm{\&lambda;}\right)}-1)+\underset{i=j-d_{\mathrm{\&lambda;}}-M+1}{\overset{j}{\mathrm{\&Sigma;}}}{\left[T\left(i\right)\right]}^{\left(\mathrm{\&lambda;}\right)}\mathrm{\&Delta;u}(k+j-i-{\mathrm{\&tau;}}^{\left(\mathrm{\&lambda;}\right)})$

$+{[S\left(j\right)\&CenterDot;B_{22}]}^{\left(\mathrm{\&lambda;}\right)}{\mathrm{\&Delta;u}}^2(k-{\mathrm{\&tau;}}^{\left(\mathrm{\&lambda;}\right)}-1)+\underset{i=j-d_{\mathrm{\&lambda;}}-M+1}{\overset{j}{\mathrm{\&Sigma;}}}{\left[T_g\left(i\right)\right]}^{\left(\mathrm{\&lambda;}\right)}{\mathrm{\&Delta;u}}^2(k+j-i-{\mathrm{\&tau;}}^{\left(\mathrm{\&lambda;}\right)})$

其中 $d_{\mathrm{\&lambda;}}=\max \{{\mathrm{\&tau;}}_{\mathrm{\&lambda;}1},{\mathrm{\&tau;}}_{\mathrm{\&lambda;}2},\&CenterDot;\&CenterDot;\&CenterDot;,{\mathrm{\&tau;}}_{{\mathrm{\&lambda;n}}_u},\}(\mathrm{\&lambda;}=\mathrm{1,2},\&CenterDot;\&CenterDot;\&CenterDot;,n_y),$ 即d_λ代表时滞矩阵中第λ行所有元素中的最大值，其它符合含义与式（15）中相同表示第λ个输出量的预测值，[S(j+1)]^(λ)表示矩阵[S(j+1)]中的第λ行λ=1,2，其它符号依此类推。

最后将上式中的系统矩阵[S(j+1)]、[S(j)·A₂]、[S(j)·B₁₂]、[T(i)]、[S(j)·B₂₂]、[T_g(i)]进行分行交叉变换组合得到，含有时滞和稳态增益K的预测模型向量表达式，如式（16）所示：

$\mathrm{\&Delta;}\hat{Y}\left(k\right)=G_{11}\mathrm{\&Delta;U}\left(k\right)+G_{12}\mathrm{\&Delta;}{U}^2\left(k\right)+G_{21}\mathrm{\&Delta;}{U}_2\left(k\right)+G_{22}{\mathrm{\&Delta;U}}_2^2\left(k\right)+F_1\mathrm{\&Delta;y}\left(k\right)+F_2\mathrm{\&Delta;y}(k-1)---16)$

式中，

G₁₁,G₁₂∈R^{(ny·P)×(nu·M)}；ΔU(k),ΔU²(k)∈R^(nu·M)×1，F₁,F₂∈R^(ny·P)×ny； $G_{21},G_{22}\&Element;R^{(\mathrm{ny}\&CenterDot;P)\&times;[\mathrm{nu}\&CenterDot;(\mathrm{ny}+\underset{i=1}{\overset{\mathrm{ny}}{\mathrm{\&Sigma;}}}{d}_i)]};$ $\mathrm{\&Delta;}{U}_2\left(k\right),{\mathrm{\&Delta;U}}_2^2\left(k\right)\&Element;R^{[\mathrm{nu}\&CenterDot;(\mathrm{ny}+\underset{i=1}{\overset{\mathrm{ny}}{\mathrm{\&Sigma;}}}{d}_i)]\&times;1};$ Δy(k),Δy(k-1)∈R^ny×1；d_λ=max{τ_λ1,τ_λ2,…,τ_λnu}(λ=1,2,…,ny)。其中P为预测时域，M为控制时域，τ₁₂为u₂对y₁的时滞步数，为被控量预测值的增量矩阵，Δy(k),Δy(k-1)被控量当前步k和上一步k-1时刻的增量矩阵，ΔU(k),ΔU²(k)为控制量预测值的增量矩阵，ΔU₂(k),

为控制量过去时滞时间段内的增量矩阵，G₁₁,G₁₂,G₂₁,G₂₂,F₁,F₂是包含有稳态增益K和动态ARX模型参数A₁,A₂,B₁,B₂的参数矩阵。

2．预测未来输出状态

在当前时刻k利用由1中所辨识出的式（16）BP-ARX预测模型对未来一段时间内水泥回转窑煅烧过程的输出量状态进行预测，得到输出量的预测增量值

再由预测增量值经位移矩阵转换，得到输出量的预测值为：

${\hat{Y}}_n\left(k\right)=V\&CenterDot;\mathrm{\&Delta;}\hat{Y}\left(k\right)+W\&CenterDot;y\left(k\right)---\left(17\right)$

式中 ${\hat{Y}}_n\left(k\right)={\left[\begin{array}{cccc}{y}_n(k+1|k)& y_n(k+2|k)& \&CenterDot;\&CenterDot;\&CenterDot;& y_n(k+P|k)\end{array}\right]}^T,$ 其中y_n(k+i|k)表示基于时刻k的未来第i时刻的回转窑煅烧过程预测输出量状态，y_1n(k+i|k)和y_2n(k+i|k)的行向量，i=1,2,…,P，其中P为预测时域，

$W={\left[\begin{array}{cccc}I& I& \&CenterDot;\&CenterDot;\&CenterDot;& I\end{array}\right]}_{(n_y\&times;P)\&times;n_y}^T,$ I为n_y×n_y的单位矩阵，y(k)为k时刻的煅烧过程的输出向量值。

3．误差反馈闭环输出

由图1中BP-ARX模型（1）得到的回转窑煅烧过程预测输出值和真实的水泥回转窑（5）的输出值之间总是存在一定的误差，为克服模型输出的误差以及在运行状态下的各种干扰，引入输出量的反馈校正是必要的，计算k时刻的输出误差e(k)=y(k)-y_n(k|k-1)，将误差补偿到模型的预测输出值上，就可以得到闭环预测输出值：

${\hat{Y}}_c\left(k\right)={\hat{Y}}_n\left(k\right)=\mathrm{He}\left(k\right)---\left(18\right)$

其中Y_c(k)=[y_c(k+1|k)  y_c(k+2|k)  …  y_c(k+P|k)]^T，其中y_c(k+i|k)表示基于时刻k的未来第i时刻的回转窑煅烧过程闭环预测输出量状态，y_c(k+i|k)是一个包含了y_1c(k+i|k)和y_2c(k+i|k)的行向量，H为校正矩阵，取校正矩阵中的参数矢量为 $H={\left[\begin{array}{cccc}{h}_{j1}& h_{j2}& \&CenterDot;\&CenterDot;\&CenterDot;& h_{{\mathrm{jn}}_y}\end{array}\right]}_{n_y\&times;n_y},h_j={\mathrm{diag}\left\{\begin{array}{cccc}{h}_{j1}& h_{j2}& \&CenterDot;\&CenterDot;\&CenterDot;& h_{{\mathrm{jn}}_y}\end{array}\right\}}_{n_y\&times;n_y}(j=\mathrm{1,2},\&CenterDot;\&CenterDot;\&CenterDot;,P),$ 取h_ji为1。

4．设定参考轨迹

引入参考轨迹的目的是使系统输出能沿着其设定轨迹平稳地达到设定值，参考轨迹的形式如下：

$Y_r={\left[\begin{array}{cccc}{y}_r^T(k+1)& y_r^T(k+2)& \&CenterDot;\&CenterDot;\&CenterDot;& y_r^T(k+P)\end{array}\right]}^T---\left(19\right)$

式中：y_r(k+j)=Cy_r(k+j-1)+(I-C)y_sp，(j=1,2,…,P)，y_r(k)=y(k)，其中y(k)是系统在k时刻的实际测量输出值， $C=\mathrm{diag}{\left\{\begin{array}{cccc}{c}_1& c_2& \&CenterDot;\&CenterDot;\&CenterDot;& \end{array}{c}_{n_y}\right\}}_{n_y\&times;n_y}$ 是柔化系数矩阵，单个输出量的柔化系数0≤c_i≤1，(i=1,2,…,ny)，柔化系统对闭环系统的动态特性和鲁棒性起很关键的作用，y_sp=[y_1sp  y_2sp  …  y_nysp]^T是系统输出量的设定值向量。

5．非线性预测控制器输出

如图1中预测控制器（4）的运算实质是对目标函数求最优解的过程，依据输入输出量的预测值和历史值构建带有约束的二次型目标函数，如式（20）所示：

$\min J\left(k\right)={\left|\right|{\hat{Y}}_c\left(k\right)-Y_r\left(k\right)\left|\right|}_Q^2+{\left|\right|\mathrm{\&Delta;U}\left(k\right)\left|\right|}_R^2$

s.t.Δu_min≤Δu(k+i)≤Δu_max  （20）

u_min≤u(k+i)≤u_max,(i=0,1,…,M-1)

y_min≤y(k+j)≤y_max,(j=1,2,…,P)

式中，Q=qI、R=rI分别是输出误差和控制增量的加权矩阵，I是单位矩阵，

是被控量修正后的预测值矩阵，Y_r(k)是被控量的参考轨迹矩阵，Δu是控制量增量矩阵，u和y为控制量矩阵和被控量矩阵，Δu_min,Δu_max,u_min,u_max,y_min,y_max为对应变量的上下限约束值。

对带有约束的优化问题目标函数式（20）的求优，采用序列二次规划（SQP）算法进行求解。求得式（20）中最有控制输入增量ΔU^*(k)，取其第一组最优控制增量Δu(k)，可获得当前时刻的最优控制输入值：

u(k)=u(k-1)+Δu(k)  （21）

式中， $\mathrm{\&Delta;u}\left(k\right)={\left[\begin{array}{cccc}\mathrm{\&Delta;}{u}_1\left(k\right)& \mathrm{\&Delta;}{u}_2\left(k\right)& \&CenterDot;\&CenterDot;\&CenterDot;& \mathrm{\&Delta;}{u}_{n_u}\left(k\right)\end{array}\right]}^T.$ 将求得的当前时刻的控制输入u(k)作用到系统中去，再进行下一步的优化求解，如此循环执行下去即完成了预测控制的在线滚动优化求解。

综上所述，整个控制方法流程如下：

1．通过DCS系统中数据存储装置采集现场水泥回转窑煅烧过程变量数据，并将其筛选分类；

2．选取BP稳态网络模型和动态ARX模型的初始参数，然后分别通过BP算法和最小二乘法对其进行离线模型训练，再将训练得到的BP和ARX模型通过增益进行并联结合，最后将BP-ARX预测模型投入在线运行；

3．选择设定预测控制参数P，M，C，Q，R，以及输入输出变量的约束值Δu_min,Δu_max,u_min,u_max,y_min,y_max。

4．在当前时刻k，依据式（19）计算出水泥回转窑煅烧过程中参考输出量NO_x含量y_1r(k+j)和O₂含量y_2r(k+j)；

5．计算输出量误差e(k)=y(k)-y_n(k|k-1)；

6．由BP-ARX预测模型输出NO_x含量y_1n(k+j)和O₂含量y_2n(k+j)，经反馈校正后生成闭环预测输出NO_x含量y_1c(k+j)和O₂含量y_2c(k+j)

7．利用序列二次规划法求解二次型目标函数J，获得最优解Δu(k+j-1)，选取当前时刻k作为输入量信号u(k)，作为水泥回转窑煅烧过程煤和高温风机挡板开度的输入量，然后转至第4步，直到完成整个控制。

实现上述基于BP-ARX模型的预测控制方法的预测控制系统，如图3所示，该预测控制系统由智能仪表（9）、上位机（13）、智能执行器（14）和DCS系统组成，DCS系统又由数据接口（10）、数据存储装置（11）和控制站（12）组成。智能仪表（9）和智能执行器（14）直接和水泥回转窑（5）相连，DCS系统由通过现场总线分别和智能仪表和智能执行器相连，上位机又通过总线与DCS系统相连。

如图4所示，现场智能仪表（9）直接采集水泥回转窑（5）煅烧过程输出量NO_x含量和O₂含量的状态，并将实时信息通过DCS系统中数据接口传送至数据存储装置并将其存储；上位机（13）中嵌入了上述的预测控制算法，并通过总线从DCS系统中的数据存储装置中采集所需的模型训练数据，又通过自己的预测控制算法计算出水泥回转窑煅烧过程的预测输入量；将预测出的回转窑煅烧过程的输入量通过DCS系统调节智能执行器（14），通过智能执行器（14）将输入量的预测值作用到水泥回转窑（5）上，此时水泥回转窑的煅烧过程的输出状态再由智能仪表（9）往上反映，如此循环下去，实现水泥回转窑煅烧过程在预测控制方法的调节下，在预测控制系统的支撑下，完成自动控制。

Claims (7)

1.一种基于BP-ARX时滞模型的水泥回转窑煅烧的预测控制方法，其特征在于该预测控制方法包括以下步骤：

步骤一：确定煅烧过程的主要输入输出量并获得分类辨识数据，影响煅烧过程的主要输入量窑头喷煤量和高温风机挡板开度，输出量窑内NO_x含量和氧气含量，采集相应变量的现场运行数据存放于数据存储装置，并将其分类为稳态数据和动态数据；

步骤二：建立预测模型，利用BP神经网络对分类出的稳态数据进行辨识得到水泥回转窑煅烧过程的非线性BP稳态模型，利用最小二乘法对分类出的动态数据进行辨识得到煅烧过程的线性时滞ARX动态模型，并通过将稳态模型和动态模型进行有机并联组合，得到水泥回转窑煅烧过程的MIMO时滞系统BP-ARX预测模型；

步骤三：预测未来输出状态，在当前时刻k，利用水泥回转窑煅烧过程的历史和未来数据的输入输出信息，通过基于BP-ARX模型的预测控制模型，并将其进行迭代对水泥回转窑煅烧过程未来一段时间的输出状态窑内NO_x含量y_1n(k+j)、氧气含量y_2n(k+j)进行预测；

步骤四：误差反馈校正，将从步骤三中得到的回转窑煅烧过程未来的输出状态y_in(k+j)和k时刻的输出误差e_i(k)相加，得到基于k时刻回转窑煅烧过程未来的闭环预测输出窑内NO_x含量y_1c(k+j)、氧气含量y_2c(k+j)；

步骤五：设定输出量参考轨迹，为使煅烧过程未来输出量能够沿着预先设定轨迹平稳地达到设定值，引入参考输出轨迹窑内NO_x含量y_1r(k+j)、氧气含量y_2r(k+j)；

步骤六：滚动优化求解，将得到的预测闭环输出y_ic(k+j)与设定的参考轨迹y_ir(k+j)进行比较，构建带约束的二次型目标函数，并通过序列二次规划法对其进行滚动优化求解，计算出当前时刻应加于水泥回转窑的喷煤量u₁(k)、高温风机挡板开度u₂(k)；

步骤七：根据步骤六计算出的当前时刻输入量控制现场智能执行器，实现现场水泥回转窑煅烧过程的自动控制。

2.根据权利要求1所述的一种基于BP-ARX时滞模型的水泥回转窑煅烧预测控制方法，其特征在于步骤一所述的影响水泥回转窑煅烧过程的主要变量，是采用窑内NO_x含量间接反映窑内烧成带温带，用窑内氧气含量反映窑内煤的燃烧状况。

3.根据权利要求1所述的一种基于BP-ARX时滞模型的水泥回转窑煅烧预测控制方法，其特征在于步骤二所述的水泥回转窑煅烧过程的MIMO时滞系统预测模型是稳态模型和含有时滞的动态模型并联而成，非线性稳态模型是通过模糊神经网络根据回转窑煅烧过程的稳态数据辨识获得，线性时滞动态模型是通过最小二乘法对煅烧过程的动态数据辨识获得。

4.根据权利要求1所述的一种基于BP-ARX时滞模型的水泥回转窑煅烧预测控制方法，其特征在于步骤三所述的预测模型迭代求解水泥回转窑煅烧过程未来一段时间内的输出预测值，采用矩阵分离变换的方法来完成对MIMO时滞系统未来预测输出量状态y_in(k+j)的求解。

5.根据权利要求1所述的一种基于BP-ARX时滞模型的水泥回转窑煅烧预测控制方法，其特征在于步骤六所述预测控制器在滚动时域下采用序列二次规划法对带有约束的二次型目标函数进行求解，保证求出的煅烧过程预测输入解是在约束条件下的全局最优解。

6.一种实现权利要求1所述的一种基于BP-ARX时滞模型的水泥回转窑煅烧过程预测控制方法的控制系统，其特征在于：所述的控制系统包括现场智能仪表、智能执行器、数据存储装置及上位机；智能仪表用于采集水泥回转窑煅烧过程的输出量，即窑内NO_x含量和氧气含量并将其传送到数据存储装置；所述的数据存储装置为DCS系统的数据存储装置，所述的DCS系统包括数据接口、控制站和数据存储装置；上位机用于运行回转窑煅烧过程预测控制算法，根据水泥回转窑煅烧过程的输出量计算出当前时刻应加入回转窑内煤和风的控制量，并通过对智能执行器的调节来实现现场回转窑煅烧过程的自动控制。

7.如权利要求6所述的基于BP-ARX时滞模型的水泥回转窑煅烧过程预测控制系统，其特征在于：智能仪表和智能执行器直接与现场相连，智能仪表和智能执行器又通过现场总线与DCS系统、上位机依次相连。

Images