-
摘要: 针对一类在有限时间区间上重复作业的不确定非线性系统，本文提出一种重复学习控制方法，用于解决非参数不确定性问题.该方法采用死区修正学习律对期望控制输入与界函数进行估计，以避免参数的正向累加导致系统发散，并使该控制算法较少地依赖于系统信息，更方便于控制器的实现.基于Lyapunov方法所设计的控制器，保证了闭环系统所有信号的有界性，并使得跟踪误差收敛于死区界定的邻域.通过仿真算例及电机实验结果验证所提学习控制算法的有效性.
-
关键词:
- 重复学习控制 /
- 非参数不确定性 /
- 死区修正 /
- Lyapunov方法
Abstract: This paper presents a repetitive learning control method to handle the nonparametric uncertain problem for a class of uncertain nonlinear systems performing a given repetitive task over a finite time interval. The learning laws with dead-zone modification are adopted to estimate the desired control input and bound functions, which avoids the divergency of estimates due to the ceaseless positive accumulation and facilitates the implementation of the controller with less knowledge about the system dynamics. The repetitive learning controller is designed in terms of Lyapunov synthesis, so as to guarantee the boundedness of all closed-loop signals while ensuring the tracking error to converge to the pre-specified neighbourhood. Numerical results for an inverted pendulum system and the AC motor experiment are conducted to testify the effectiveness of the proposed learning control scheme. -
钢铁是国民经济的命脉,高炉炼铁则是钢铁生产中最为重要的环节之一. 如图 1所示,整个高炉炼铁系统分为高炉本体、给料系统、热风系统、煤粉喷吹系统、高炉煤气处理系统以及出铁系统等几个子系统,其中高炉本体又可分为炉喉、炉身、炉腰、炉腹和炉缸等几个部分. 高炉炼铁时,铁矿石、焦炭、溶剂按一定比例及布料制度逐层从高炉顶部装载到炉喉中,在其下部将预热的空气、氧气和煤粉鼓入炉缸中,空气、氧气、煤粉和焦炭在高温作用下发生一系列复杂物理化学反应,生成大量的高温还原性气体,这些还原性气体不断向上运动将铁从铁矿石中还原出来,而上行的气体最终变为高炉煤气从炉顶回收. 炉料则随着炉缸中焦炭的不断燃烧和铁水的不断滴落逐渐向下运动,在下降过程中,炉料经过加热、还原、熔化等一系列复杂的物理化学变化,最终生成液态的生铁和炉渣从出铁口排出 [1-2].
高炉炼铁的目标是实现冶炼过程的高产量与低能耗,为了实现这一目标,就应对高炉内部状态进行实时监测与有效控制. 然而高炉内部冶炼环境极其严酷,反应最剧烈的区域温度高达2 000多度,压强高达标准大气压的4倍左右,且伴随着固、液、气多相共存的状态,使高炉内部状态的实时监测难以实现,从而难以对高炉进行优化控制. 目前,被广泛用来间接反映高炉内部状态的指标为铁水质量参数,综合性的铁水质量指标通常采用硅含量([Si])、磷含量([P])、硫含量([S])和铁水温度(Molten iron temperature,MIT)来衡量. 硅含量是评定高炉炉况稳定性和生铁质量的重要指标,也是表征高炉热状态及其变化的标志之一,硅含量的变化可由多种高炉内部状态变化引起,比如,熔化区域的垂直位移、原料的变更和热量流动等 [3]. 硫和磷在钢材中均是有害元素,在高炉冶炼过程中应严格控制铁水中的硫含量和磷含量. 硫含量过高说明渣铁之间的脱硫反应不够彻底. 铁水温度是物理热的表征,是影响高炉稳定顺行及节能降耗的重要指标. 采用这四种铁水质量参数作为高炉内部状态的评判指标可以较全面地了解高炉内部的运行状态,为高炉的控制操作提供指导. 因此,要实现面向高炉铁水质量和稳定顺行的控制与优化,首先要建立铁水质量的动态模型.
高炉炼铁系统是一个物理化学反应复杂、多相多场耦合的非线性、大滞后、动态时变系统,因此,铁水质量的预测模型应该是一个非线性的动态模型. 目前,针对铁水质量建模的研究已有很多,包括基于机理分析的数学模型 [4]、基于规则的推理模型 [5]以及基于数据的统计模型 [6]. 近些年来,随着计算机的快速发展,数据驱动的建模方法已被广泛应用于解决铁水质量的建模问题. 数据驱动是一种黑箱建模方法,它的主要思想是利用一些数学工具和智能算法基于过程数据直接近似得到过程的输入输出关系而不需要任何先验知识 [7-9]. 已有的铁水质量数据驱动建模方法有偏最小二乘法 [10]、人工神经网络(Artificial neural network,ANN) [11-12]、支持向量回归机(Support vector regression,SVR) [13-14]等. 但是,这些铁水质量预测模型大多只是对单一铁水质量参数硅含量或铁水 温度的建模,并不能全面地反映高炉内部复杂的状态,且部分模型并未考虑高炉系统的非线性与动态特性. 另外,这些以铁水质量软测量为目标的输入输出智能模型用于解决铁水质量尤其是多元铁水质量控制系统设计时,总是存在诸多局限性和不足. 虽然,Zeng等 [15]采用子空间辨识方法建立了铁水硅含量的线性输入输出预测模型,Marutiram等 [16]采用线性回归的方法建立了铁水硅含量的线性ARMA (Auto regressive and moving average)预测模型,并且二者均针对建立的线性模型进行了预测控制,取得了一定的效果. 但是,高炉本身是一个极其复杂的非线性动态系统,用线性模型去近似非线性系统,并不能完全准确地表征铁水质量,因而,基于线性模型的线性控制方法也难于对多元铁水质量指标进行有效控制.
针对上述问题,为了建立一种既适合于控制器设计,又能充分体现高炉非线性动态特性,同时又无需任何先验知识的高炉炼铁多元铁水质量非线性预测模型,本文首先采用典型相关性分析与相关性分析相结合的方法选取与铁水质量相关性最强的高炉本体可控变量作为模型的输入变量,然后采用基于LS-SVM (Least square support vector machine)的非线性子空间辨识方法 [17]建立数据驱动的多元铁水质量非线性状态空间动态模型. 同时,将核函数表示的非线性特性用多项式函数拟合,以降低模型的复杂度. 最后基于实际工业现场数据进行工业试验和比较分析. 状态空间模型是现代控制理论的经典模型,针对此模型的控制方法研究已有很多,且理论也很成熟. 另外,Hammerstein模型由于其特殊的线性与非线性模块串联结构,相较于SVR和ANN等非线性智能模型,更易于设计非线性控制器. 所以本文建立的多元铁水质量模型可以直接应用于多元铁水质量的控制研究中,是一种面向控制的模型.
1. 建模策略
从提高产品质量、节约能源的角度而言,高炉系统控制与优化的主要对象是铁水质量参数,尤其是铁水[Si]和铁水温度,它们也是衡量高炉内热状态和稳定顺行的重要参数,其过高和过低对于燃料消耗和成本有较大的影响. 为此,对铁水质量参数进行建模意义重大,这也是实现铁水质量控制与优化的关键. 针对前述已有铁水质量建模存在的诸多问题,提出图 2所示的多元铁水质量非线性动态建模策略:
图 2 多元铁水质量参数非线性子空间建模策略Fig. 2 Strategy diagram of multivariable nonlinear dynamic modeling for molten iron quality1) 高炉本体参数较多,且变量之间存在着较强的相关性,若所有的变量均被用来建立铁水质量模型,势必会增大建模过程的计算复杂度,影响模型的准确性和有效性. 已有的铁水质量建模研究中,辅助变量的选取方法有机理分析、主成分分析(Principle component analysis,PCA)和相关性分析等 [11-12, 18-19],但是这些选取方法大多是针对单一铁水质量,并未考虑对多个铁水质量指标的综合影响,且已有研究中选取的变量多是针对建模问题的,在控制问题中有些变量并不适用. 所以,采用典型相关性分析方法选择与铁水质量指标相关性最大的几个可控变量,并用相关性分析方法对这些可控变量做进一步筛选,去除冗余变量.
2) 子空间辨识建模是一种由输入输出数据辨识系统状态空间模型的方法,它有三种基本算法,分别为N4SID (Numerical algorithms for subspace state space system identification) [20]、MOESP (Multi-variable output error state space) [21]、CVA (Canonical variate analysis) [22]. 子空间辨识可以在系统状态完全未知的情况下,仅由输入输出数据将系统的状态估计出来,从而构造系统的状态空间模型. 本文采用基于LS-SVM的非线性Hammerstein系统子空间建模方法 [17]建立多元铁水质量参数的状态空间模型. 该方法在传统线性子空间建模算法N4SID基础上,引入LS-SVM去解决无需任何先验知识的静态非线性环节的辨识问题,最终可建立一个能够较为准确刻画炼铁过程动态特性和非线性特性的状态空间模型. 该模型不但能准确地对多元铁水质量进行在线动态估计,并且可用于解决面向多元铁水质量的控制和优化问题.
3) 由于基于LS-SVM非线性Hammerstein子空间模型的非线性特性是用核函数表示,计算效率较低. 因此,采用多项式拟合方法 [23]拟合模型的非线性特性,即将核函数表示的非线性部分替换为计算较为简单的多项式表达式,以降低模型的计算复杂度.
2. 建模算法
2.1 基于典型相关分析和相关性分析的建模辅助变量选择
典型相关分析的基本思想与PCA类似,首先在两组变量${ X} = (x_1,x_2,\cdots,x_p)$,${ Y} = (y_1,y_2,\cdots ,y_q)$中分别找出变量的一个线性组合,即
$\left\{ \begin{array}{l} U_i ={ X}{ \alpha} ^{(i)} \ V_i = { Y}{ \beta} ^{(i)} \ \end{array} \right. $
(1) 式中,${ \alpha} ^{(i)},{ \beta} ^{(i)}$分别为线性组合的系数向量,使得两组线性组合之间具有最大的相关系数. 然后选取相关系数仅次于第一对线性组合并且与第一对线性组合不相关的第二对线性组合,如此继续下去,直到两组变量之间的相关性被提取完毕为止. $U_i ,V_i$的配对称为典型变量,它们之间的相关系数称为典型相关系数.
在进行辅助变量选择时,首先要对各对典型变量的典型相关系数进行显著性检验,若某一对典型变量的相关程度不显著,说明这对变量不具有代表性,舍弃这一对典型变量. 因为Ui中系数绝对值较大的n个变量对Ui 起决定性作用,所以若$U_i,V_i$的相关程度显著且又具有很大的相关系数,则这n个变量便与Vi,即${ Y}$中q个变量的线性组合具有很大的相关性. 因此,选取辅助变量时,先找出相关程度显著且典型相关系数较高的几对典型变量,然后选取这几个Ui中线性组合系数绝对值较大的几个变量作为候选辅助变量.
在高炉本体变量中,各个变量均具有一定的相关性,若两个变量相关性达到80 %以上,则可近似认为二者呈线性关系,不必将二者同时选为辅助变量. 因此,在确定辅助变量时,要用相关性分析的方法对候选辅助变量进一步缩减,缩减的原则是保留可控变量,舍弃不可控变量.
2.2 基于LS-SVM的非线性Hammerstein系统子空间建模
Hammerstein模型是一类非线性系统,由静态非线性环节和动态线性环节串联而成. 这种模型能较好地反映过程特征,可以描述一大类非线性系统,如无线发射器 [24]、聚丙烯牌号切换系统 [25]、燃料电池 [26]、生物系统 [27]等. 在状态空间模型输入端加入静态非线性函数 即构成Hammerstein模型,如式(2)所示,
${{\left\{ \begin{align} & {{x}_{t+1}}=A{{x}_{t}}+Bf({{u}_{t}})+{{\nu }_{t}} \\ & {{y}_{t}}=C{{x}_{t}}+Df({{u}_{t}})+v \\ \end{align} \right.}_{t}}$
(2) 式中,$t=1,\cdots ,N-1,{{u}_{t}}\in {{R}^{m}},{{y}_{t}}\in {{R}^{l}},{{x}_{t}}\in {{R}^{n}}$分别为系统的输入、输出和状态,$f(u)={{[{{f}_{1}}(u(1))\text{ }{{f}_{2}}(u(2))\cdots {{f}_{m}}(u(m))]}^{\text{T}}},{{\nu }_{t}}\in {{R}^{n}},{{v}_{t}}\in {{R}^{l}}$为零均值高斯白噪声序列向量,其协方差矩阵为
$\text{E}\left[ \begin{matrix} {{\nu }_{p}} \\ {{v}_{p}} \\ \end{matrix} \right]\left[ {{\nu }_{q}}^{\text{T}}{{v}_{q}}^{\text{T}} \right]=\left[ \begin{matrix} Q & S \\ {{S}^{\text{T}}} & R \\ \end{matrix} \right]{{\delta }_{pq}}$
(3) 在文献[17]中,Goethals等以预测为目的提出了一种针对如式(2)所示的多输入多输出Hammerstein系统的系统辨识方法. 通过将N4SID算法中的斜投影计算问题改写为一系列分量形式的LS-SVMs回归问题,从而对N4SID算法进行非线性扩展. Goethals等指出Hammerstein系统中的线性模型和式(2)中的静态非线性特性可以通过对LS-SVMs回归问题得到的矩阵进行低秩逼近得到. 具体建模过程描述如下.
2.2.1 计算斜投影
假设输入输出数据对$\left\{ {\left( {{ u_t},{ y_t} } \right)} \right\}_{t = 0}^{N - 1} $可以测量得到,定义$\Phi$为针对Hankel矩阵 [15]的矩阵算子,定义非线性函数$\rho ( \cdot )$为针对Z中每一个块矩阵$Z_{i}$ 的非线性算子. 根据子空间辨识建模原理,斜投影$O_i = {Y_{f/U_f} }W_p,W_p = [Y_p ;U_p]$可以通过$L_u,L_y,f$计算求得,式中$U_f,Y_f,U_p,Y_p$ 分别为未来输入输出数据和过去输入输出数据对应的Hankel矩阵 [15],而$L_u ,L_y,f$则可以通过一个最小二乘问题估计得出,即最小化式(4)中的残差E.
${{Y}_{f}}=\left[ \begin{matrix} {{L}_{u}} & {{L}_{y}} \\ \end{matrix} \right]\left[ \begin{matrix} {{\Phi }_{f}}({{U}_{0\left| 2i-1 \right.}}) \\ {{Y}_{p}} \\ \end{matrix} \right]+E$
(4) 式(4)亦可写为
$\begin{align} & {{Y}_{f}}(s,t)={{L}_{y}}(s,:){{Y}_{p}}(:,t)+ \\ & \sum\limits_{h=1}^{2i}{\sum\limits_{k=1}^{m}{{{L}_{u}}(s,(h-1)m+k){{f}_{k}}({{u}_{h+t-2}}(k))}}+ \\ & E(s,t) \\ \end{align}$
(5) 式中,$s = 1,\cdots,li,t = 1,\cdots,j$. 一旦式(4)、(5)中的$\hat L_u,\hat L_y,\hat f$的估计值通过最小二乘法确定下来,斜投影可计算如下:
$\begin{align} & {{O}_{i}}(s,t)={{{\hat{L}}}_{y}}(s,:){{Y}_{p}}(:,t)+ \\ & \sum\limits_{h=1}^{i}{\sum\limits_{k=1}^{m}{{{{\hat{L}}}_{u}}(s,(h-1)m+k){{{\hat{f}}}_{k}}({{u}_{h+t-2}}(k))}} \\ \end{align}$
(6) 然而式(5)中包含参数矩阵Lu与非线性函数f的相乘项,使得优化问题非凸,为此引入一组函数\${{g}_{h,s,k}}:R\to R$,
$\begin{align} & {{g}_{h,s,k}}{{f}_{k}}={{c}_{h,s,k}}{{f}_{k}} \\ & \begin{matrix} \text{s}.\text{t}. & {{c}_{h,s,k}}={{L}_{u}}(s,(h-1)m+k) \\ \end{matrix} \\ \end{align}$
(7) 式中,$h = 1,\cdots,2i$,$s = 1,\cdots,li$,$k = 1,\cdots,m$. 引入这组函数后,式(5)、(6)可写为
$\begin{align} & {{Y}_{f}}(s,t)={{L}_{y}}(s,:){{Y}_{p}}(:,t)+ \\ & \sum\limits_{h=1}^{2i}{\sum\limits_{k=1}^{m}{{{g}_{h,s,k}}({{u}_{h+t-2}}(k))+E(s,t)}} \\ \end{align}$
(8) $\begin{align} & {{O}_{i}}(s,t)={{{\hat{L}}}_{y}}(s,:){{Y}_{p}}(:,t)+ \\ & \sum\limits_{h=1}^{i}{\sum\limits_{k=1}^{m}{{{{\hat{g}}}_{h,s,k}}({{u}_{h+t-2}}(k))}} \\ \end{align}$
(9) 定义核函数${{k}^{k}}:R\times R\to R$,对所有$p,q = 0,\cdots,N - 1$和$k = 1,\cdots,m$,均有$k^k ({ u_p (k)},{ u_q (k)}) = \phi _k ({ u_p (k)})^{\rm T} \phi _k ({ u_q (k)})$,核函数矩阵${{K}^{k}}\in {{R}^{N\times N}}$,对所有$i,j = 1,\cdots,N$,均有$K^k (i,j) = k^k ({ u_{i - 1} (k)},{ u_{j - 1} (k)})$. 用$\omega _{h,s,k}^{\rm T} \phi _k $替换式(8)中的$g_{h,s,k} $得
$\begin{align} & {{Y}_{f}}(s,t)={{L}_{y}}(s,:){{Y}_{p}}(:,t)+ \\ & \sum\limits_{h=1}^{2i}{\sum\limits_{k=1}^{m}{\omega _{h,s,k}^{\text{T}}{{\phi }_{k}}({{u}_{h+t-2}}(k))+E(s,t),}} \\ & \forall s=1,\cdots ,li,t=1,\cdots ,j \\ \end{align}$
(10) 式(10)的形式可以方便地引入核函数及分量形式(Componentwise) LS-SVM的方法去估计参数.
考虑到将非线性函数展开成一组非线性函数的和并不唯一,引入中心约束条件$\sum_{t = 0}^{N - 1} {\bar f({ u_t} )} = 0$. 对式(2)通过引入如下状态变换$\left\{ \begin{align} & {{\xi }_{t}}=\psi ({{x}_{t}})={{x}_{t}}-{{(I-A)}^{-1}}B{{\delta }_{u}} \\ & {{\delta }_{y}}=(C{{(I-A)}^{-1}}B+D){{\delta }_{u}} \\ \end{align} \right.$,得到
$\left\{ \begin{align} & {{\xi }_{t+1}}=A{{\xi }_{t}}+B\bar{f}({{u}_{t}})+{{\nu }_{t}} \\ & {{y}_{t}}-{{\delta }_{y}}=C{{\xi }_{t}}+D\bar{f}({{u}_{t}})+{{v}_{t}} \\ \end{align} \right.$
(11) 式中,$\psi :{{R}^{n}}\to {{R}^{n}}$,${ \delta _u} $,${ \delta _y} $均为常数向量. 由于系统模型加入了一个新的参数${ \delta _y} $,因此式(10)变换为
$\begin{align} & {{Y}_{f}}(s,t)\text{+}[{{\text{1}}_{i}}\otimes {{\delta }_{y}}](s)={{L}_{y}}(s,:)({{Y}_{p}}(:,t)+{{1}_{i}}\otimes {{\delta }_{y}})+ \\ & \sum\limits_{h=1}^{2i}{\sum\limits_{k=1}^{m}{\omega _{h,s,k}^{\text{T}}{{\phi }_{k}}({{u}_{h+t-2}}(k))+E(s,t),}} \\ & \forall s=1,\cdots ,li,t=1,\cdots ,j \\ \end{align}$
式中$ \otimes $表示 Kronecker 积. 对于所有$h = 1,\cdots,2i$,$s = 1,\cdots,li$,均有$\omega _{h,s,k}^{\rm T} \phi _k = g_{h,s,k} = c_{h,s,k}^{} f_k $,则原LS-SVM回归问题可表示为一个有约束的优化问题
$\begin{align} & \underset{{{\omega }_{h,s,k}},{{L}_{y}},E,{{\delta }_{y}}}{\mathop{\min }}\,J({{\omega }_{h,s,k}},{{L}_{y}},E,{{\delta }_{y}})= \\ & \frac{1}{2}\sum\limits_{s=1}^{il}{\sum\limits_{h=1}^{2i}{\sum\limits_{k=1}^{m}{\omega _{h,s,k}^{\text{T}}{{\omega }_{h,s,k}}}}}\text{ + }\frac{\gamma }{2}\sum\limits_{s=1}^{il}{\sum\limits_{t=1}^{j}{E{{(s,t)}^{2}}}} \\ & \text{s}.\text{t}.{{Y}_{f}}(s,t)\text{+}[{{\text{1}}_{i}}\otimes {{\delta }_{y}}](s)= \\ & {{L}_{y}}(s,:)({{Y}_{p}}(:,t)+{{1}_{i}}\otimes {{\delta }_{y}})+ \\ & \sum\limits_{h=1}^{2i}{\sum\limits_{k=1}^{m}{\omega _{h,s,k}^{\text{T}}{{\phi }_{k}}({{u}_{h+t-2}}(k))}}+ \\ & E(s,t),\begin{matrix} & \forall s=1,\cdots ,li,t=1,\cdots ,j \\ \end{matrix} \\ & \sum\limits_{t=0}^{N-1}{\omega _{h,s,k}^{\text{T}}{{\phi }_{k}}({{u}_{t}}(k))=0,} \\ & \begin{matrix} & \forall h\text{ = }1,\cdots ,2i,s=1,\cdots ,li, \\ \end{matrix} \\ & \begin{matrix} & k=1,\cdots ,m \\ \end{matrix} \\ \end{align}$
(12) 通过构造式(12)的拉格朗日函数(Lagrangian),$L_y $和${ \delta _y} $的估计值可以从如下对偶系统求解
$\left[ {\begin{array}{*{20}{c}} 0&0&{{1^{\rm{T}}}}&0\\ 0&0&{{Y_p}}&0\\ 1&{Y_p^{\rm{T}}}&{{\kappa _p} + {\kappa _f} + {\gamma ^{ - 1}}}&S\\ 0&0&{{S^{\rm{T}}}}&T \end{array}} \right]\left[ {\begin{array}{*{20}{c}} \begin{array}{l} {\bar d}\\ L_y^{\rm{T}}\\ A\\ B \end{array} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} \begin{array}{l} 0\\ 0\\ Y_f^{\rm{T}}\\ 0 \end{array} \end{array}} \right]$
(13) 式中,$\bar d = ({ 1_i} \otimes I_l - L_y ({ 1_i} \otimes I_l )){ \delta _y} $,${ 1_j} $是长度为j且元素全为1的列向量,$A\in {{R}^{j\times li}},B\in {{R}^{2im\times li}}$分别为包含一组Lagrange乘子$\alpha _{s,t},\beta _{h,s}$的矩阵,$T,S$定义如下:
$\begin{align} & T:={{I}_{2i}}\otimes \left[ \begin{matrix} {{1}_{N}}^{\text{T}}{{K}^{1}}{{1}_{N}} & \cdots & 0 \\ \vdots & \ddots & \vdots \\ 0 & \cdots & {{1}_{N}}^{\text{T}}{{K}^{m}}{{1}_{N}} \\ \end{matrix} \right] \\ & S:=\left[ \begin{matrix} {{S}_{1}} & {{S}_{2}} & \cdots & {{S}_{2i}} \\ {{S}_{2}} & {{S}_{3}} & \cdots & {{S}_{2i+1}} \\ \vdots & \vdots & \ddots & \vdots \\ {{S}_{j}} & {{S}_{j+1}} & \cdots & {{S}_{N}} \\ \end{matrix} \right] \\ & {{S}_{q}}:=\sum\limits_{t=1}^{N}{\left[ \begin{matrix} {{K}^{1}}(t,q) & {{K}^{2}}(t,q) & \cdots & {{K}^{m}}(t,q) \\ \end{matrix} \right]} \\ \end{align}$
矩阵${{\kappa }_{p}}\in {{R}^{j\times j}},{{\kappa }_{f}}\in {{R}^{j\times j}}$,其元素分别为
$\kappa _p (p,q) = \sum\limits_{h = 1}^i {\sum\limits_{k = 1}^m {k^k ({ u_{h + p - 2} (k)},{ u_{h + q - 2} (k)})} } $
$\begin{align} & {{\kappa }_{f}}(p,q)=\sum\limits_{h=i+1}^{2i}{\sum\limits_{k=1}^{m}{{{k}^{k}}({{u}_{h+p-2}}(k),{{u}_{h+q-2}}(k)),}} \\ & p,q=1,\cdots ,j \\ \end{align}$
$g_{h,s,k} $的估计值以及斜投影分别如下
$\begin{align} & {{{\hat{g}}}_{h,s,k}}:R\to R:{{u}^{*}}(k)\to \\ & \sum\limits_{t=1}^{j}{{{\alpha }_{s,t}}{{k}^{k}}({{u}_{h+t-2}}(k),{{u}^{*}}(k))}+ \\ & {{\beta }_{h,s}}(k)\sum\limits_{t=0}^{N-1}{{{k}^{k}}({{u}_{t}}(k),{{u}^{*}}(k)),}\forall h,s \\ \end{align}$
(14) $\begin{align} & {{O}_{i}}={{L}_{u}}(:,1:li){{\Phi }_{{\underset{\raise0.3em\hbox{$\smash{\scriptscriptstyle-}$}}{f}}}}({{U}_{0\left| i-1 \right.}})+{{{\hat{L}}}_{y}}{{Y}_{p}}= \\ & {{A}^{\text{T}}}{{\kappa }_{p}}+B_{p}^{\text{T}}S_{p}^{\text{T}}+{{{\hat{L}}}_{y}}({{Y}_{p}}-({{1}_{i}}{{1}_{j}}^{\text{T}})\otimes {{\delta }_{y}}) \\ \end{align}$
(15) 式中,$B_p = B(1:mi,:)$,$S_p {\text{ = }}S(:,1:mi)$. $O_{i + 1} $的计算方法和Oi类似,即
$\begin{align} & {{O}_{i+1}}={{({{A}^{-}})}^{\text{T}}}{{(\kappa _{p}^{+})}^{\text{T}}}+{{(B_{p}^{-})}^{\text{T}}}{{(S_{p}^{-})}^{\text{T}}}+ \\ & L_{y}^{-}(Y_{p}^{+}-{{1}_{(i+1)}}{{1}_{j}}^{\text{T}}\otimes {{\delta }_{y}}) \\ \end{align}$
(16) 式中,{ $\kappa _p^ + (p,q) = \sum\nolimits_{h = 1}^{i + 1} {\sum\nolimits_{k = 1}^m {k^k ({ u_{h + p - 2} (k)},{ u_{h + q - 2 }}}}$ ${{{ (k)})} } $,$p,q = 1,\cdots,j$,$B_p^ - = B^ - (1:(i + 1)m,:)$,$S_p^ - = S^ - (:,1:(i + 1)m)$. $A^- ,B^- ,L_y^- $}可通过对偶系统求解.
2.2.2 系统状态估计
系统的状态序列$\tilde X_i,\tilde X_{i + 1} $的估计方法和线性子空间辨识的方法一样,可以通过对斜投影$O_i,O_{i +1} $进行SVD (Singular value decomposition)分解得到. 这些状态序列将在算法的后续步骤中用来估计系统矩阵和提取非线性特性$\bar f_k$.
2.2.3 提取系统矩阵和静态非线性函数
系统矩阵和静态非线性可由式(17)估计得出
$\begin{array}{l} (\hat A,\hat B,\hat C,\hat D,\hat f) = \arg \mathop {\min }\limits_{A,B,C,D,f} {\mkern 1mu} \left\| {\left[ {\begin{array}{*{20}{c}} {{{\tilde X}_{i + 1}}{Y_{i\left| i \right.}}} \end{array}} \right] - } \right.\\ \left[ {\begin{array}{*{20}{c}} A&B\\ C&D \end{array}} \right]\left. {\left[ {\begin{array}{*{20}{c}} \begin{array}{l} {{\tilde X}_i}\\ {\Phi _f}({U_{i\left| i \right.}}) \end{array} \end{array}} \right]} \right\|_F^2 \end{array}$
(17) 如同上节,这个最小二乘问题依然可以变为一个LS-SVM回归问题,定义
$\begin{array}{*{20}{l}} {{\chi _{i + 1}} = \left[ {\begin{array}{*{20}{c}} {\begin{array}{*{20}{l}} {{{\tilde X}_{i + 1}}}\\ {{Y_{i\left| i \right.}} - {\delta _y}} \end{array}} \end{array}} \right],{\Theta _{AC}} = \left[ {\begin{array}{*{20}{l}} A\\ C \end{array}} \right]}\\ {{\Theta _{BD}} = \left[ {\begin{array}{*{20}{c}} {\begin{array}{*{20}{l}} B\\ D \end{array}} \end{array}} \right]} \end{array}$
(18) 将$\Theta _{BD} (s,k) \bar f_k$用$\omega _{s,k}^{\rm T} \phi _k $替换,则
$\begin{align} & {{\chi }_{i+1}}(s,t)={{\Theta }_{AC}}(s,:){{{\tilde{X}}}_{i}}(:,t)+ \\ & \sum\limits_{k=1}^{m}{\omega _{s,k}^{\text{T}}{{\phi }_{k}}({{u}_{i+t-1}}(k))}+E \\ \end{align}$
式中,E为式(17)的残差. 则LS-SVM回归的原始问题可以写为
$\begin{align} & \underset{{{\omega }_{s,k}},E,{{\Theta }_{AC}}}{\mathop{\min }}\,J(\omega ,E)= \\ & \frac{1}{2}\sum\limits_{s=1}^{n+l}{\sum\limits_{k=1}^{m}{\omega _{s,k}^{\text{T}}}}{{\omega }_{s,k}}+\frac{{{\gamma }_{BD}}}{2}\sum\limits_{s=1}^{n+l}{\sum\limits_{t=1}^{j}{E{{(s,t)}^{2}}}} \\ & \text{s}.\text{t}.{{\chi }_{i+1}}(s,t)={{\Theta }_{AC}}(s,:){{{\tilde{X}}}_{i}}(:,t)+ \\ & \sum\limits_{k=1}^{m}{\omega _{s,k}^{\text{T}}{{\phi }_{k}}({{u}_{i+t-1}}(k)),} \\ & \begin{matrix} & \forall s=1,\cdots ,li,t=1,\cdots ,j \\ \end{matrix} \\ & \sum\limits_{t=0}^{N-1}{\omega _{s,k}^{\text{T}}{{\phi }_{k}}({{u}_{t}}(k))\text{ = }0,} \\ & \begin{matrix} & \forall s=1,\cdots ,li,k=1,\cdots ,m \\ \end{matrix} \\ \end{align}$
式中,$\gamma _{BD} $为一个正则化常数,与第2.2.1节中的$\gamma $相区别.
式(18)所示$\Theta _{AC} $中的系统矩阵$A,C$可由如下所示的对偶问题求出:
$\left[ {\begin{array}{*{20}{c}} 0&{{{\tilde X}_i}}&0\\ {\tilde X_i^{\rm{T}}}&{{\kappa _{BD}} + \gamma _{BD}^{ - 1}I}&{{S_{BD}}}\\ 0&{S_{BD}^{\rm{T}}}&{{T_{BD}}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} \begin{array}{l} \Theta _{AC}^{\rm{T}}\\ {A_{BD}}\\ {B_{BD}} \end{array} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} \begin{array}{l} 0\\ \chi _{i + 1}^{\rm{T}}\\ 0 \end{array} \end{array}} \right]$
(19) 式中,$\forall p,q = 1,\cdots,j$,$K_{BD}^{k}(p,q)={{k}^{k}}({{u}_{i+p-1}}(k),\{{{u}_{i+q-1}}(k))$,$\kappa _{BD} (p,q) = \sum\nolimits_{k = 1}^m {K_{BD}^k (p,q)} $,${{A}_{BD}}\in {{R}^{j\times (n+l)}},{{B}_{BD}}\in {{R}^{m\times (n+l)}}$,$T_{BD} {\text{ = diag\{ }}{ 1}^\text{T}_{N} K^k \times\%[-3.5mm] { 1_N}\},k = 1,\cdots,m$,$S_{BD} = \left[{S_{i + 1} \cdots S_{i + j} } \right]^{\text{T}} $. 系统矩阵$A,C$的估计值可以通过分解$\Theta _{AC} $得到,将式(19)与式(17)、式(18)相结合,得到下式
$\begin{align} & {{\Theta }_{BD}}(:,k)[{{{\bar{f}}}_{k}}({{u}_{0}}(k)),{{{\bar{f}}}_{k}}({{u}_{1}}(k)),\cdots ,{{{\bar{f}}}_{k}}({{u}_{N-1}}(k))]= \\ & A_{BD}^{\text{T}}{{K}^{k}}(i+1:i+j,:)+ \\ & B_{BD}^{\text{T}}(k,:)\sum\limits_{t=1}^{N}{K_{BD}^{k}(t,:)} \\ \end{align}$
(20) 对$\Theta _{BD} $中的系统矩阵$B,D$的第k列和非线性特性$\bar f_k $的估计可通过奇异值分解(SVD)得到.
2.3 多项式拟合模型静态非线性特性
由奇异值分解(SVD)获得非线性特性$\bar f_k $后,可以得到N组数据,即$\left\{ {\left( {{ u_t (k)},\bar f_k ({ u_t (k)})} \right)} \right\}_{t = 0}^{N - 1} $,将${ u_t (k)},\bar f_k ({ u_t (k)})$分别作为横、纵坐标画出非线性特性$\bar f_k $对应的散点图,将这些散点用一条多项式曲线拟合,即可获得非线性特性$\bar f_k $的具体函数表达式,具体算法如下:
首先,根据非线性特性$\bar f_k $的散点图判断出多项式的阶次m,以此决定$\bar f_k $的近似多项式表达式为:
${{\bar{f}}_{k}}(x)={{a}_{0}}+{{a}_{1}}x+{{a}_{2}}{{x}^{2}}+\cdots +{{a}_{m}}{{x}^{m}}$
然后,将多项式曲线拟合问题转化为如下最小二乘问题求解
$\underset{{{a}_{0}},{{a}_{1}},{{a}_{2}},\cdots ,{{a}_{m}}}{\mathop{\min }}\,\sum\limits_{t=0}^{N-1}{{{\left\| \sum\limits_{i=0}^{m}{{{a}_{i}}{{({{u}_{t}}(k))}^{i}}-{{{\bar{f}}}_{k}}({{u}_{t}}(k))} \right\|}^{2}}}$
(21) 利用多元函数求极值方法,分别对$a_i $求偏导,可得$m + 1$个方程,写成矩阵形式如式(22)所示.
$\begin{array}{l} \left[ {\begin{array}{*{20}{c}} N&{\sum\limits_{t = 0}^{N - 1} {{u_t}(k)} }& \cdots \\ {\sum\limits_{t = 0}^{N - 1} {{u_t}(k)} }&{\sum\limits_{t = 0}^{N - 1} {{{({u_t}(k))}^2}} }& \cdots \\ \vdots & \vdots & \ddots \\ {\sum\limits_{t = 0}^{N - 1} {{{({u_t}(k))}^m}} }&{\sum\limits_{t = 0}^{N - 1} {{{({u_t}(k))}^{m + 1}}} }& \cdots \end{array}} \right.\\ \left. {\begin{array}{*{20}{c}} {\sum\limits_{t = 0}^{N - 1} {{{({u_t}(k))}^m}} }\\ {\sum\limits_{t = 0}^{N - 1} {{{({u_t}(k))}^{m + 1}}} }\\ \vdots \\ {\sum\limits_{t = 0}^{N - 1} {{{({u_t}(k))}^{2m}}} } \end{array}} \right] \times \left[ {\begin{array}{*{20}{c}} \begin{array}{l} {a_0}\\ {a_1}\\ \vdots \\ {a_m} \end{array} \end{array}} \right] = \\ {\rm{ }}\left[ {\begin{array}{*{20}{c}} \begin{array}{l} \sum\limits_{t = 0}^{N - 1} {{{\bar f}_k}({u_t}(k))} \\ \sum\limits_{t = 0}^{N - 1} {{u_t}(k){{\bar f}_k}({u_t}(k))} \\ \vdots \\ \sum\limits_{t = 0}^{N - 1} {{{({u_t}(k))}^m}{{\bar f}_k}({u_t}(k))} \end{array} \end{array}} \right] \end{array}$
(22) 求解如式(22)所示的关于$a_0,a_1,\cdots,a_m$的线性方程组,即可获得非线性特性$\bar f_k $的多项式函数表达式. 从而,Hammerstein模型中用核函数表示的静态非线性特性可由多项式函数取代,降低了模型的复杂度.
2.4 建模算法实现步骤
基于前述讨论,最终的输入变量选择及非线性子空间辨识建模过程可以总结如下:
1) 根据式(1)针对输入输出数据求典型相关性系数.
2) 进行显著性检验和相关性分析,以此确定建模输入变量.
3) 根据式(15)和式(16)计算斜投影$O_i $和$O_{i + 1} $.
4) 通过对斜投影$O_i,O_{i + 1} $进行SVD分解估计出系统的状态.
5) 根据式(19)获得系统矩阵A,C及$A_{BD}$,$B_{BD} $的估计值.
6) 通过对式(20)进行秩为1的近似,获得系统矩阵$B,D$的第k列及非线性特性$\bar f_k$.
7) 采用多项式拟合的方法将核函数表示的非线性特性$\bar f_k$确定为具体的多项式函数表达式.
3. 工业实验
3.1 模型建立
文基于柳钢2\#高炉的实际生产数据进行工业试验,建立多元铁水质量模型. 由生产现场数据库可得18个高炉主体参数分别为: 送风比、热风压力、顶压、压差、顶压风量比、透气性、阻力系数、热风温度、富氧流量、富氧率、设定喷煤量、鼓风湿度、理论燃烧温度、标准风速、实际风速、鼓风动能、炉腹煤气量、炉腹煤气指数. 在建立高炉多元铁水质量模型之前,首先要从这18个高炉主体参数中确定出对多元铁水质量影响最大的几个辅助变量. 应用第2.1节基于典型相关分析和相关性分析的方法进行建模输入变量选择. 首先,采用典型相关分析的方法针对柳钢2\#高炉2013年10月10日到18日的200组(采样间隔为1 h)高炉本体数据与铁水质量数据进行分析,去除相关性不显著的1对典型变量,得到相关系数分别为0.715、0.571、0.531的3对典型变量,典型相关分析结果如表 1所示.
表 1 典型相关分析结果Table 1 The results of canonical correlation analysis0.563[Si]-0.453[P]- 0.638[Si]+0.899[P]- -0.857[Si]-0.368[P]- 0.447[S]+0.287MIT 0.986[S]-0.518MIT 0.801[S]-0.398MIT 送风比 -1.103 -3.254 0.716 热风压力 0.587 -1.109 -1.493 顶压0 .974 -0.911 -0.320 压差 0.495 -7.118 2.506 顶压风量 -0.948 -0.794 3.185 透气性- 2.687 -4.104 -2.782 阻力系数 -2.702 6.501 -5.392 热风温度 -5.375 14.185 4.230 富氧流量 4.914 -4.014 -0.054 富氧率 -7.284 8.123 -0.377 设定喷煤量 0.998 -7.335 -3.443 鼓风湿度 0.096 0.208 -0.071 理论燃烧温度 3.789 -12.762 -5.639 标准风速 -0.439 -3.953 -0.717 实际风速 1.816 2.624 3.444 鼓风动能 -0.438 -4.662 -2.728 炉腹煤气量 -1.306 14.932 4.029 炉腹煤气指数 -0.176 -1.503 -0.335 选出表 1中对输出典型变量影响较大的几个高炉主体参数作为候选辅助变量,即压差、阻力系数、热风温度、富氧流量、富氧率、设定喷煤量、理论燃烧温度和炉腹煤气量. 然后通过相关性分析的方法选出候选辅助变量中相关性较大的几组中的可调控变量,如富氧流量与富氧率的相关系数为0.999,选取富氧率;理论燃烧温度和富氧率的相关系数为0.78,舍弃理论燃烧温度. 压差和阻力系数均为模型间接计算值,考虑到压差和阻力系数均通过热风压力计算得出,且热风压力可控,所以选取热风压力替代压差和阻力系数. 炉腹煤气量由富氧流量、设定喷煤量等可调控量计算得来,为间接可调控量,所以选为辅助变量. 最终确定用于建模的可调可控辅助变量为热风压力、热风温度、富氧率、设定喷煤量和炉腹煤气量. 另外,为了反映高炉的非线性动态特性以及输入输出变量的时序关系,将上一采样时刻辅助变量测量值以及上一采样时刻多元铁水质量值,连同当前采样时刻输入辅助变量检测值共同作为动态模型的综合输入.
高炉是一个强噪声系统,因此在建模前首先要对工业现场采集到的实际生产数据进行数据预处理. 针对由于高炉炉况不稳定和检测仪器不精确造成的跳变数据,采用异常值检测算法剔除高炉生产过程中的噪声尖峰跳变数据; 然后采用移动平均滤波算法减弱训练数据中的高斯噪声干扰. 之后基于预处理后的数据,采用本文所提方法进行多元铁水质量参数的非线性状态空间模型建模,并采用遗传算法优化模型可调参数$\sigma,\gamma,\gamma _{BD} $. 当可调参数确定为$\sigma = 1,\gamma = 0.67,\gamma _{BD} = 1$时,得到的模型如式(23)所示:
$\left\{ {\begin{array}{*{20}l} {{ x_{t + 1}} = A{ x_t} + B\bar f({ u_t} )} \ {{ y_t} - { \delta _y} = C{ x_t} + D\bar f({ u_t} )} \ \end{array} } \right. $
(23) 式中:$A = \left[ {\begin{array}{*{20}{c}} {0.9556}&{ - 0.0508}&{0.1664}&{0.1121}\\ {0.0367}&{ - 0.8629}&{ - 0.1471}&{ - 0.0606}\\ { - 0.1443}&{ - 0.1605}&{ - 0.9077}&{ - 0.1780}\\ { - 0.1199}&{0.0295}&{0.0987}&{0.9188} \end{array}} \right]$,
$B = \left[ {\begin{array}{*{20}{c}} {0.1612}&{0.0862}&{ - 0.2201}&{0.0680}&{0.3208}\\ {0.0354}&{0.3049}&{ - 0.3723}&{0.0601}&{0.3720}\\ {0.2419}&{0.0999}&{ - 0.0251}&{0.0020}&{0.0946}\\ {0.0539}&{ - 0.0590}&{0.0919}&{0.0821}&{ - 0.1328} \end{array}} \right]$,
$C = \left[ {\begin{array}{*{20}{c}} { - 0.1882}&{ - 1.2661}&{0.2607}&{0.3425}\\ {0.0767}&{ - 1.5404}&{ - 0.4250}&{0.5169}\\ {1.4999}&{0.0204}&{ - 0.1522}&{ - 0.5670}\\ { - 0.6268}&{ - 0.5679}&{1.1473}&{1.0014} \end{array}} \right]$,
$D = \left[ {\begin{array}{*{20}{c}} { - 0.0481}&{0.6082}&{0.1316}&{ - 0.7082}&{ - 0.8085}\\ {0.6109}&{0.1919}&{ - 0.8541}&{ - 0.5200}&{ - 0.1834}\\ { - 0.4919}&{ - 0.6901}&{0.2145}&{ - 0.1354}&{ - 0.0373}\\ {0.5421}&{0.0569}&{0.1050}&{ - 0.4414}&{ - 0.2087} \end{array}} \right]$,
${\delta _y} = \left[ \begin{array}{l} - 0.3875\\ 0.4098\\ - 0.2016\\ - 0.1499 \end{array} \right]$,辨识得到的非线性特性散点分布如图 3所示,由图中散点分布形状选取多项式阶次为4,采用第2.2.4节的多项式拟合方法得到非线性函数如式(24)所示,相应拟合曲线如图 3所示.
3.2 建模及其估计效果
图 4为建立的非线性状态空间模型在多项式拟合非线性特性前后对多元铁水质量参数(即[Si],[P],[S]和MIT)的建模效果比较,图 5为建立的非线性状态空间模型在多项式拟合非线性特性前后对实际数据的多元铁水质量参数估计效果比较. 可以看出多项式拟合非线性特性的式(23)所示模型及用核函数表示非线性特性的模型得到的多元铁水质量预测值的大小和变化趋势几乎完全一致,两种表示形式的非线性子空间辨识建立的多元铁水质量模型均具有非常高的建模和估计精度,建模误差以及估计误差均非常小,且变化趋势与其实际值非常一致,即模型估计输出与其实际值已基本吻合. 另外,由于建立的模型既能对多元铁水质量输出进行准确估计,并且还能输出系统的状态估计曲线,因而可更好地用于各种控制算法的研究和设计.
图 4 基于非线性子空间辨识的多元铁水质量模型在多项式拟合非线性特性前后的建模结果Fig. 4 Modeling results of molten iron quality prediction with and without polynomial ¯tting
图 5 模型在多项式拟合非线性特性前后对多元铁水质量参数的实际估计效果Fig. 5 Estimation results of molten iron quality prediction with and without polynomial ¯tting$\bar{f}({{u}_{t}})=\left[ \begin{matrix} \begin{align} & -0.1074{{({{u}_{t}}(1))}^{4}}+0.1063{{({{u}_{t}}(1))}^{3}}+0.3132{{({{u}_{t}}(1))}^{2}}-0.2080({{u}_{t}}(1))-0.0089 \\ & -0.0701{{({{u}_{t}}(2))}^{4}}-0.0739{{({{u}_{t}}(2))}^{3}}+0.2320{{({{u}_{t}}(2))}^{2}}+0.3586({{u}_{t}}(2))-0.1099 \\ & 0.0454{{({{u}_{t}}(3))}^{4}}-0.1261{{({{u}_{t}}(3))}^{3}}-0.1029{{({{u}_{t}}(3))}^{2}}+0.2870({{u}_{t}}(3))-0.0177 \\ & 0.0499{{({{u}_{t}}(4))}^{4}}-0.1090{{({{u}_{t}}(4))}^{3}}-0.1598{{({{u}_{t}}(4))}^{2}}+0.4549({{u}_{t}}(4))\text{ + }0.0545 \\ & -0.0646{{({{u}_{t}}(5))}^{4}}-0.2026{{({{u}_{t}}(5))}^{3}}+0.2720{{({{u}_{t}}(5))}^{2}}+0.5439({{u}_{t}}(5))-0.2821 \\ \end{align} \\ \end{matrix} \right]$
(24) 为了进一步说明多项式拟合非线性特性的模型相比于核函数表示非线性特性模型的优越性,计算二者在实际数据上预测的均方根误差(Root mean square error,RMSE)及计算时间进行定量分析. 在实际工业测试数据计算得到的两种模型的估计RMSE和计算时间如表 2所示. 可以看出两种模型的预测精度相差非常小,核函数表示非线性特性的模型略好于多项式拟合非线性特性的模型,但是在计算时间上多项式拟合非线性特性的模型要明显优于核函数表示非线性特性的模型,计算效率提高了十几倍. 因此,综合以上的分析和评价,可以看出多项式拟合非线性特性的多元铁水质量非线性Hammerstein模型可以在损失非常小建模精度的前提下,大大降低模型复杂度,提高实际应用过程中的计算效率.
表 2 两种模型对各个铁水质量指标估计的均方根误差和计算时间比较Table 2 Comparison between two models in RMSE for molten iron quality prediction and computation timeRMSE Time (s) [Si] (%) [P] (%) [S] (%) MIT (±C) Kernel 0.0633 0.0032 0.0021 6.2388 0.01123 Polynomial 0.0664 0.0031 0.0022 6.6907 0.00061 为了进一步验证本文方法的优越性,将本文建模方法即多项式拟合非线性特性的模型(Nonlinear subspace model,N-SM)与文献[20]提出的线性子空间辨识(Linear subspace model,L-SM)方法、文献[28]提出的多输出LS-SVR (Multi-output least square vector regression,M-LS-SVR)方法以及常见的RLS-ARMA (Recursive least square-auto regressive and moving average)方法进行比较研究. 这几种方法根据实际工业数据对多元铁水质量的估计效果如图 6所示. 可以看出,本文模型得到的各个铁水质量参数的估计值相比其他方法能够更好地拟合真实值的变化趋势,且能保证较高的精度,具有更好的估计性能. 为了进一步说明这个问题,引入误差自相关函数来评价不同方法的建模性能. 众所周知,一个好的模型,它的估计误差自相关函数曲线应该近似为一个白噪声序列. 所以,画出不同建模方法得到的模型的估计误差自相关函数曲线,如图 7所示. 可以看出,相对于L-SM、M-LS-SVR以及RLS-ARMA这几种算法,本文模型4个铁水质量参数估计误差自相关曲线综合来说更加接近于白噪声序列.
图 6 不同建模方法对多元铁水质量的实际估计效果对比Fig. 6 Estimation results of molten iron quality prediction with di®erent models
图 7 不同模型估计误差自相关函数Fig. 7 Autocorrelation function of estimating error of di®erent models另外,为了进一步对比不同模型的估计精度,将铁水质量的实际值与模型估计值分别作为横、纵坐标,画出上述各个方法对各个指标的估计值与实际值散点图,并求出各个散点图回归线的斜率.由图 8可以看出,相较于L-SM、M-LS-SVR以及RLS-ARMA这几种建模方法,本文建立的模型得到的各个铁水质量参数散点更加集中地分布在 对角线附近,且回归线的斜率更加接近于1,这进一步表明本文模型得到的多元铁水质量估计值更加接近于各个指标参数的实际值. 评价一个模型的精度通常还会选用一些标准的统计指标,如均方误差(MSE)、均方根误差(RMSE)等. 本文选用RMSE作为评价指标对4种不同模型的估计精度进行定量评估. 在实际工业测试数据计算得到的各个模型的估计RMSE如表 3所示. 可以看出本文模型对各个铁水质量指标的RMSE数值均为最小,这进一步验证了建立的模型具有最好的估计精度.
图 8 不同建模方法的估计值与实际值散点图Fig. 8 Scatter diagram of estimated and actual value by di®erent models表 3 多元铁水质量估计值均方根误差比较Table 3 RMSE for molten iron quality predictionRMSE [Si] (%) [P] (%) [S] (%) MIT (±C) N-SM 0.0664 0.0031 0.0022 6.6907 L-SM 0.1584 0.0107 0.0073 8.4801 M-LS-SVR 0.0811 0.0058 0.0051 8.0037 RLS-ARMA 0.0721 0.0060 0.0044 7.2992 3.3 应用于多元铁水质量非线性预测控制
为了进一步验证本文方法建立的多元铁水质量模型的有效性,将建立的Hammerstein状态空间模型作为预测模型应用于多元铁水质量的非线性预测控制. 预测控制的性能指标如式(25)所示.
$\begin{align} & J={{\sum\limits_{k=1}^{{{N}_{p}}}{\left( {{r}_{t+k}}-{{{\tilde{y}}}_{t+k|t}} \right)}}^{\text{T}}}\left( {{r}_{t+k}}-{{{\tilde{y}}}_{t+k|t}} \right)\text{+} \\ & \sum\limits_{k=1}^{{{N}_{c}}}{\Delta u_{t+k-1}^{\text{T}}R\Delta {{u}_{t+k-1}}} \\ \end{align}$
(25) 选取预测时域$N_p $为1,控制时域$N_c $为1,控制项权重R为0.001. 多元铁水质量的初始设定值为[Si] = 0.45 %,[P] = 0.115 %,[S] = 0.025 %,MIT = 1 500 $^{\circ}$C,在20、40、60、80时刻分别更改设定值,将设定值分别更改为MIT = 1 510 $^{\circ}$C,[S] = 0.015 %,[P] = 0.125 %,[Si] = 0.50 %,在100、\120、140、160、180时刻分别对热风压力($u_1$)、热风温度($u_2$)、富氧率($u_3$)、设定喷煤量($u_4$)、炉腹煤气量($u_5$)加入较大的脉冲干扰,得到的多元铁水质量响应曲线如图 9所示,控制量曲线如图 10所示. 可以看出4个铁水质量参数均能够快速地跟随设定值的阶跃变化,且达到稳态时与设定值的偏差较小,在对输入变量加入较大干扰时,铁水质量参数能够从较大的波动中迅速地回到设定值,有效地抑制干扰.实验结果验证了基于本文模型的非线性预测控制在多元铁水质量控制上的有效性.
图 10 多元铁水质量预测控制控制量曲线Fig. 10 Control input curves of predictive control of molten iron quality parameters综合,可以看出本文建立的模型应用于多元铁水质量在线估计与预测时,具有较好的估计性能. 另外,本文方法是一种面向控制的非线性系统状态空间建模方法,建立的多元铁水质量非线性状态空间模型不仅可以用于多元铁水质量参数的在线估计,为工业现场的操作人员提供指导,并且还可进一步应用于铁水质量的控制与优化.
4. 结论
高炉炼铁的最终任务是生产出符合要求的铁水.因此,面向[Si]、[P]、[S]和铁水温度等多元铁水质量参数的建模、控制与优化历来是高炉炼铁生产及其控制工程非常重要的一环.本文针对现有铁水质量参数建模的诸多不足,提出一种面向控制的数据驱动高炉炼铁多元铁水质量参数非线性子空间建模方法. 首先,为了提高建模效率和降低计算复杂度,采用典型相关性分析与相关性分析相结合的方法从众多关联变量中提取与铁水质量相关性最强的5个关键可控变量作为建模输入变量;然后基于实际工业现场数据,采用基于LS-SVM的非线性Hammerstein系统子空间辨识方法,建立数据驱动的多元铁水质量非线性状态空间动态模型,并将核函数表示的非线性特性用多项式函数拟合,以降低模型的复杂度. 工业试验结果表明:建立的模型不仅能够对多元铁水质量参数进行准确在线估计,而且可很好地用于面向铁水质量参数的控制研究.
-
图 10 界函数估计$\hat l_{f, \, k}$
Fig. 10 Estimate of the bound function $\hat l_{f, \, k}$
图 11 界函数估计$\hat l_{g, \, k}$
Fig. 11 Estimate of the bound function $\hat l_{g, \, k}$
图 12 误差性能指标$J_k$ (其中三条虚线为控制器(31)的实验结果, 实线为控制器(11)的实验结果
Fig. 12 Error performance index $J_k$ (the three dotted lines are the result by controller (31), the solid line is the result by controller (11))
图 13 控制输入$u_k$ (其中三条虚线为控制器(31)的实验结果, 实线为控制器(11)的实验结果)
Fig. 13 Control input $u_k$ (the three dotted lines are the result by controller (31), the solid line is the result by controller (11))
图 14 参考输入的估计$\hat u_k$ (其中三条虚线为控制器(31)的实验结果, 实线为控制器(11)的实验结果)
Fig. 14 Control input $\hat u_k$ (the three dotted lines are the result by controller (31), the solid line is the result by controller (11))
图 15 误差性能指标$J_k$ (其中三条虚线为控制器(31)的实验结果, 实线为控制器(11)的实验结果)
Fig. 15 Error performance index $J_k$ (the three dotted lines are the result by controller $ (31)$, the solid line is the result by controller $ (11)$)
图 16 控制输入$u_k$ (其中三条虚线为控制器(31)的实验结果, 实线为控制器(11)的实验结果)
Fig. 16 Control input $u_k$ (the three dotted lines are the result by controller (31), the solid line is the result by controller (11))
-
[1] Arimoto S, Kawamura S, Miyazaki F. Bettering operation of robots by learning. Journal of Robotic Systems, 1984, 1(2):123-140 doi: 10.1002/(ISSN)1097-4563 [2] Mezghani M, Roux G, Cabassud M, Le Lann M V, Dahhou B, Casamatta G. Application of iterative learning control to an exothermic semibatch chemical reactor. IEEE Transactions on Control Systems Technology, 2002, 10(6):822-834 doi: 10.1109/TCST.2002.804117 [3] Tayebi A, Chien C J. A Unified adaptive iterative learning control framework for uncertain nonlinear systems. IEEE Transactions on Automatic Control, 2007, 52(10):1907-1913 doi: 10.1109/TAC.2007.906215 [4] Chien C J, Yao C Y. Iterative learning of model reference adaptive controller for uncertain nonlinear systems with only output measurement. Automatica, 2004, 40(5):855-864 doi: 10.1016/j.automatica.2003.12.009 [5] Xu J X. A survey on iterative learning control for nonlinear systems. International Journal of Control, 2011, 84(7):1275-1294 doi: 10.1080/00207179.2011.574236 [6] Norrlof M. An adaptive iterative learning control algorithm with experiments on an industrial robot. IEEE Transactions on Robotics and Automation, 2002, 18(2):245-251 doi: 10.1109/TRA.2002.999653 [7] Kuc T Y, Kwanghee N, Lee J S. An iterative learning control of robot manipulators. IEEE Transactions on Robotics and Automation, 1991, 7(6):835-842 doi: 10.1109/70.105392 [8] Choi J Y, Lee J S. Adaptive iterative learning control of uncertain robotic systems. IEEE Proceedings-Control Theory and Applications, 2000, 147(2):217-223 doi: 10.1049/ip-cta:20000138 [9] Sun M X, Ge S S, Mareels I M Y. Adaptive repetitive learning control of robotic manipulators without the requirement for initial repositioning. IEEE Transactions on Robotics, 2006, 22(3):563-568 doi: 10.1109/TRO.2006.870650 [10] French M, Rogers E. Non-linear iterative learning by an adaptive Lyapunov technique. International Journal of Control, 2000, 73(10):840-850 doi: 10.1080/002071700405824 [11] Tayebi A. Adaptive iterative learning control for robot manipulators. Automatica, 2004, 40(7):1195-1203 doi: 10.1016/j.automatica.2004.01.026 [12] Yin C K, Xu J X, Hou Z S. A high-order internal model based iterative learning control scheme for nonlinear systems with time-iteration-varying parameters. IEEE Transactions on Automatic Control, 2010, 55(11):2665-2670 doi: 10.1109/TAC.2010.2069372 [13] Sun M, Ge S S. Adaptive repetitive control for a class of nonlinearly parametrized systems. IEEE Transactions on Automatic Control, 2006, 51(10):1684-1688 doi: 10.1109/TAC.2006.883028 [14] Dong W J, Kuhnert K D. Robust adaptive control of nonholonomic mobile robot with parameter and nonparameter uncertainties. IEEE Transactions on Robotics, 2005, 21(2):261-266 doi: 10.1109/TRO.2004.837236 [15] Chen W S, Jiao L C. Adaptive tracking for periodically time-varying and nonlinearly parameterized systems using multilayer neural network. IEEE Transactions on Neural Networks, 2010, 21(2):345-351 doi: 10.1109/TNN.2009.2038999 [16] 庆吕庆.抑制初态误差影响的自适应迭代学习控制.自动化学报, 2015, 41(7):1365-1372 http://www.aas.net.cn/CN/abstract/abstract18710.shtmlLv Qing. Adaptive iterative learning control for inhibition effect of initial state random error. Acta Automatica Sinica, 2015, 41(7):1365-1372 http://www.aas.net.cn/CN/abstract/abstract18710.shtml [17] Li D, Li J M. Adaptive iterative learning control for nonlinearly parameterized systems with unknown time-varying delay and unknown control direction. International Journal of Automation and Computing, 2012, 9(6):578-586 doi: 10.1007/s11633-012-0682-9 [18] 陶洪峰, 霰学会, 杨慧中.输入饱和非线性系统的周期自适应补偿学习控制.自动化学报, 2014, 40(9):1998-2004 http://www.aas.net.cn/CN/abstract/abstract18471.shtmlTao Hong-Feng, Xian Xue-Hui, Yang Hui-Zhong. Periodic adaptive compensating learning control of nonlinear systems with saturated input. Acta Automatica Sinica, 2014, 40(9):1998-2004 http://www.aas.net.cn/CN/abstract/abstract18471.shtml [19] Zhang C L, Li J M. Adaptive iterative learning control of non-uniform trajectory tracking for strict feedback nonlinear time-varying systems. International Journal of Automation and Computing, 2014, 11(6):621-626 doi: 10.1007/s11633-014-0819-0 [20] Jin X, Huang D Q, Xu J X. Iterative learning control for systems with nonparametric uncertainties under alignment condition. In: Proceedings of the 51st Conference on Decision and Control (CDC). Maui, HI, USA: IEEE, 2012. 3942-3947 [21] Yu M, Huang D Q, He W. Robust adaptive iterative learning control for discrete-time nonlinear systems with both parametric and nonparametric uncertainties. International Journal of Adaptive Control and Signal Processing, 2016, 30(7):972-985 doi: 10.1002/acs.v30.7 [22] Zhang C L, Li J M. Adaptive iterative learning control for nonlinear pure-feedback systems with initial state error based on fuzzy approximation. Journal of the Franklin Institute, 2014, 351(3):1483-1500 doi: 10.1016/j.jfranklin.2013.11.018 [23] Xu J X, Jin X, Huang D Q. Composite energy function-based iterative learning control for systems with nonparametric uncertainties. International Journal of Adaptive Control and Signal Processing, 2014, 28(1):1-13 doi: 10.1002/acs.2380/full [24] Li X F, Huang D Q, Chu B, Xu J X. Robust iterative learning control for systems with norm-bounded uncertainties. International Journal of Robust and Nonlinear Control, 2016, 26(4):697-718 doi: 10.1002/rnc.v26.4 [25] Xu J X, Yan R. On initial conditions in iterative learning control. IEEE Transactions on Automatic Control, 2005, 50(9):1349-1354 doi: 10.1109/TAC.2005.854613 期刊类型引用(5)
1. 徐东波，武志涛. 基于延迟补偿的永磁直线同步电机伺服系统双闭环鲁棒重复控制. 电机与控制应用. 2024(11): 64-74 . 
百度学术2. 施卉辉，陈强. 一类不确定系统的自适应滑模迭代学习控制. 控制理论与应用. 2023(07): 1162-1171 . 百度学术3. 苏杰，曾喆昭. 非线性时变系统的自耦PID控制方法. 控制理论与应用. 2022(02): 299-306 . 百度学术4. 付晓东，陈力. 全柔性空间机器人运动振动一体化输入受限重复学习控制. 力学学报. 2020(01): 171-183 . 百度学术5. 陈强，余歆祺. 一类非参数不确定系统的自适应重复学习控制. 控制理论与应用. 2020(06): 1349-1357 . 百度学术其他类型引用(4)
-
下载:
下载:
计量
- 文章访问数: 2001
- HTML全文浏览量: 262
- PDF下载量: 748
- 被引次数: 9
