CN100507460C - 基于脉冲响应模板和参数优化的动态软测量建模方法 - Google Patents

Abstract

本发明属于生产过程软测量技术领域。其特征是，基于脉冲响应原理建立卷积模型，脉冲响应参数作为主要的辨识对象，在生产过程机理分析或实际经验的基础上产生一组脉冲响应模板以及模板参数约束范围，使用优化方法辨识脉冲响应模板参数，完成动态软测量模型参数的优化求解。本发明的优点：1)适应动态变化生产过程下的软测量建模和应用的需求。2)适用于难测主导变量往往是人工取样化验情况下的动态软测量建模，对易测辅助变量也仅需要相对于难测主导变量采样时及前一个过渡过程内的相应采样数据，解决了现有的以神经网络模型为代表的的基于实验模型的动态软测量方法只能用于连续等周期采样的问题。3)使得脉冲响应参数辨识仅限定在指定的可以根据机理分析和经验确定的模板类和模板参数约束范围内，极大减少了学习参数以节省学习时间。

Description

基于脉冲响应模板和参数优化的动态软测量建模方法

技术领域

本发明涉及生产过程中一种基于脉冲响应模板和优化方法的动态软测量方法。该方法适合于生产过程的通用软测量建模。

背景技术

随着生产技术的发展和生产过程的日益复杂，为确保生产装置安全、保证产品质量，推动产品质量的直接闭环控制、质量约束和安全约束控制的广泛应用，对产品质量指标等目前还不可视的生产装置重要过程变量提出了对某些重要难测过程变量的实时测量或估计要求。由于技术或经济上的原因，很难通过在线过程仪表进行测量，因此，软测量及其应用技术是解决这些难测过程变量实时测量的关键。

软测量或软仪表方法，是对一些难测的重要变量，选择另外一些容易测量的变量，通过构成某种数学关系用计算机软件来推断和估计，以用来代替在线过程仪表功能。按建模方法的不同，软测量方法可以分为两类：一类称为基于工艺机理建模的方法。通常在对工艺机理较为清晰的应用场合，由于工程背景明确，与一般工艺设计和计算关系密切，相应的软测量模型较为简单。由于不同生产过程工艺机理不同，很难形成通用的软测量技术，建模代价较高。另一类是基于实验建模方法。该类方法是形成通用软测量技术的重要途径，获得了更为广泛的研究和应用。其中，在生产过程获得广泛应用的主要有两种方法：一种是基于回归分析方法。例如多线性回归(MLR)、基于主元分析的主元回归法((PCA-PCR)或偏最小二乘法(PLS)等；另一种是基于神经网络的方法。例如多层前向网络(MFN)、径向基网络(RBF)、延迟动态神经网络、动态递归神经网络等。

目前广泛应用的基于实验建模的软测量方法基本上是静态软测量方法。这是由于工业传感器技术或成本等方面的原因，输出变量常常是离线采样数据(通常，在一个过渡过程很可能只能得到有限几个数据，甚至一个)，采样频度很难满足软测量模型的动态建模。由于静态软侧量模型学习和在线计算，都要求输人、输出变量平稳，这对于输人变化对产品质量这样一个存在大延迟和容量滞后的过程，要求在采样前很长的过渡过程周期内都保持平稳是不现实的，因此静态软测量方法得到的估计精度很难满足生产的要求。

目前使用的动态软测量模型主要采用动态神经网络，通过在网络两端加上时间延迟单元引入输入与输出的动态过程信息，形式上转化为静态神经网络。同样由于稀疏离线采样，网络模型中所需的输出时延信号是无法获取的，即使不考虑输出时延，面对一个多变量系统，大量的输入时延将使得构造的网络极其复杂，训练模型的时间会很长，同时也难以保证训练精度及其泛化性能。因此上述动态神经网络模型很难在线使用。

文献[马勇，黄德先，金以慧.动态软测量建模方法初探.化工学报，2005，(56)8：1516-1519]对基于卷积模型的动态软测量方法进行了尝试性的研究。由于该方法在仿模预处理单元中所给出估计的脉冲响应曲线过于简单，使得其与实际情况偏差甚远，后期又通过大量的学习网络进行补偿校正，因此该方法不仅训练非常复杂，而且在不同的生产过程中很难形成一种满足通用的应用方法。

发明内容

本发明的目的：提出一种通用的基于脉冲响应模板和参数优化的动态软测量建模方法，解决对生产过程中重要难测变量的在线实时测量，克服了现有软测量方法在上述背景下无法或很难推广使用的问题。

本发明的特征在于，该方法是在计算机中依次按以下步骤进行：

步骤(1)、初始化

在生产过程工艺机理分析或实际经验的基础上构建一个难测主导变量o和n个易测辅助变量u₁，u₂，...，u_n构成的n个独立子系统，n≥1，该各个子系统的结构类型为下述两个结构类型之中的任何一种：

一阶惯性环节加纯延迟： $G\left(s\right)=\frac{{\mathrm{Ke}}^{-\mathrm{\τs}}}{\mathrm{Ts}+1}$

二阶惯性环节加纯延迟； $G\left(s\right)=\frac{K_p{e}^{-{\mathrm{\τ}}_{1}s}}{(T_{1}s+1)(T_{2}s+1)}$

其中，K，K_p为比例系数，

T，T₁，T₂为时间常数，

τ，τ₁为纯时延常数，

K，K_p，T，T₁，T₂，τ，τ₁均为待估的模型参数，每个模型参数的估计区间在生产过程工艺机理分析或实际经验的基础上给定；

步骤(2)、从在线采集的数据库以及离线化验数据库中分析搜集相关的易测辅助变量u₁，u₂，...，u_n与难测主导变量o构成M组数据样本用于动态软测量模型的训练，并按采样时间递增的顺序形成下述相应的时间序列：

$[u_1\left(Q_{1j}\right),...,u_1\left(P_j\right);u_2\left(Q_{2j}\right),...,u_2\left(P_j\right);...;u_n\left(Q_{\mathrm{nj}}\right),...,u_n\left(P_j\right)]\⇔o\left(P_j\right),$ j＝1，...，M

其中，：左边表示系统输入，右边表示系统输出，

M为总的采样次数，代表样本数，

第j组数据样本： $[u_1\left(Q_{1j}\right),...,u_1\left(P_j\right);u_2\left(Q_{2j}\right),...,u_2\left(P_j\right);...;u_n\left(Q_{\mathrm{nj}}\right),...,u_n\left(P_j\right)]\⇔o\left(P_j\right),$

表示易测辅助变量u₁的第j个样本的采样时刻，Q_1j与P_j之间形成等差数列，u₁(Q_1j)，..，u₁(P_j)包含L₁个采样数据并作为第1个子系统的输入；Q_2j，...，P_j表示易测辅助变量u₂的第j个样本的采样时刻，Q_2j与P_j之间形成等差数列，u₂(Q_2j)，...，u₂(P_j)包含L₂个采样数据并作为第2个子系统的输入；依此类推，Q_nj，...，P_j表示易测辅助变量u_n的第j个样本的采样时刻，Q_nj与P_j之间形成等差数列，u_n(Q_nj)，...，u_n(P_j)包含L_n个采样数据并作为第n个子系统的输入；P_j表示难测主导变量o的第j个样本的采样时刻，o(P_j)为P_j时刻的采样值并作为总系统的输出，

用u_s(Q_sj)，...，u_s(P_j)表示第s个子系统的输入，s∈[1，n]，该子系统的L_s个采样数据的采样总时间不短于易测辅助变量u_s与难测主导变量o之间的一个动态过渡过程时间；

步骤(3)、对步骤(2)得到的样本数据进行坏点剔除以及标准化处理，其步骤如下：

步骤(3.1)、用相邻点绝对误差阈值法剔除坏点：从样本的第二个点开始，依次计算各点与相邻点的绝对差，若该差值超过某个设定的阈值，则剔除该点，并用被剔除点位置前后两个样本点的平均值代替；

步骤(3.2)、按下式对步骤(3.1)得到的样本数据进行标准化处理：

x^o＝(x-m_x)/σ_x，

其中，x为相应样本原始数据，x^o为相应样本对应的标准化数据，m_x、σ_x分别为该组数据样本的算术平均值和标准差；

步骤(4)、对步骤(3)得到的各组数据样本，基于脉冲响应原理对每个子系统的M组数据样本建立卷积模型，再把各个子系统的卷积模型求和，得到一个全系统的脉冲响应模型D1：

$D1:o\left(P_j\right)=\underset{s=1}{\overset{n}{\mathrm{\Σ}}}\left\{\underset{i=0}{\overset{L_s}{\mathrm{\Σ}}}{h}_s\left(i\right)u_s(P_j-i)\right\},$ j＝1，.，M，

其中，n为子系统的总数，

L_s为第s个子系统的脉冲响应的长度参数，

h_s(i)为待估的第s个子系统度的脉冲响应参数；

步骤(5)、根据步骤(1)给出的各子系统的结构类型，用Ψ(Ω_s)产生一组对应的脉冲响应D2及其中的参数范围D3：

D2：h_s(i)＝Ψ(Ω_s)，s＝1，...，n，i＝1，2，...，L_s

其中，Ω_s为代表子系统结构类型的模板参数：K，或K_p，或T，或T₁，或T₂，或τ，或τ₁，Ψ(·)为脉冲响应函数，

D3：Ω_s-Π^low≤Ω_s≤Ω_s+Π^upp，

其中，Π^low为最小阈值，Π^low>0，设定值，

Π^upp为最大阈值，Π^upp>0，设定值，

L_s由香农采样确定，其取值范围D4为：

D4：L_s-a≤L_s≤L_s+b，s＝1，...，n；20≤a≤50；30≤b≤60；

步骤(6)、利用差分进化算法与改进差分进化算法搜索最佳模板参数，实现动态软测量模型优化求解，其步骤如下：

步骤(6.1)、优化目标函数D5表示为：

$D5:\max :J=1/\underset{j=1}{\overset{M}{\mathrm{\Σ}}}{(y\left(P_j\right)-o\left(P_j\right))}^2,$

其中， $y\left(P_j\right)=\underset{s=1}{\overset{n}{\mathrm{\Σ}}}\left\{\underset{i=0}{\overset{L_s}{\mathrm{\Σ}}}\stackrel{\‾}{h_s}\left(i\right)u_s(P_j-i)\right\},$ j＝1，.，M，表示优化算法计算的输出，

$\stackrel{\‾}{h_s}\left(i\right)=\mathrm{\Ψ}\left(\stackrel{\‾}{{\mathrm{\Ω}}_s}\right),$

为优化算法给出的子系统结构类型的模板参数，模板参数指的是用D3优化时用的模型参数，

o(P_j)为对应的真实模型输出，由步骤(2)得到；

步骤(6.2)、按以下步骤进行优化计算：

步骤(6.2.1)、初始化

设定：用NP表示种群规模，30≤NP≤60，maxiter表示最大迭代步数，300≤maxiter≤1500，ε表示优化算法精度值，0.01≤ε≤0.1，D3表示优化参数矢量约束范围，F表示差分进化算法中缩放比例，F＝1，C表示交叉率，C＝0.5，FF_min、FF_max表示改进差分进化算法中最小缩放比例以及最大缩放比例，FF_min＝0.1，FF_max＝1.8，CC_min、CC_max表示改进差分进化算法中最小交叉率以及最大交叉率，CC_min＝0.1，CC_max＝0.8，根据步骤(5)设定L_s初值为L_s-a，L_s的取值上限

为L_s+b，QQ为L_s的递增步长，为5到10之间的一个整数；

步骤(6.2.2)、群体即种群规模初始化：

n个子系统中的模板参数的总数构成一个v维的个体矢量，初始种群规模由NP个v维的个体矢量组成，其中第k个个体矢量X_k是记为

$\{x_{k1},x_{k2},...,x_{\mathrm{kv}}\}=$

的数据，

其中，v为一个个体矢量的维数，即模板参数的总数，K，T，τ，K_p，T₁，T₂，τ₁的上标表示子系统序号，若子系统为一阶系统时，子系统的模型参数为K，T，τ，若子系统为二阶系统时，子系统的模型参数为K_p，T₁，T₂，τ₁，

迭代步数iter＝0，

根据下式生成初始群体：

$x_{\mathrm{kr}}^{\mathrm{iter}}=x_{\mathrm{kr}}^{\min }+\mathrm{rand}\left(\mathrm{0,1}\right)\×(x_{\mathrm{kr}}^{\max }-x_{\mathrm{kr}}^{\min });,$

其中，

代表第k个个体矢量的第r个分量所对应模板参数的最小值，

为其最大值，

表示第iter步时第k个个体矢量的第r个分量的计算值，rand(0，1)表示0与1之间的随机数；

步骤(6.2.3)、按步骤(6.1)中的D5式计算每个个体的适应度值，其中第k个个体X_k的适应度值

确定种群中最大的适应度值 $J_{\max }^{\mathrm{iter}}=\max \left(J_k^{\mathrm{iter}}\right),$ k＝1，...，NP，并记录最大的适应度值个体为gBest；

步骤(6.2.4)、判断是否满足精度要求，或者是否达到最大迭代步数，若条件

或iter>maxiter满足，则输出个体gBest，并执行步骤(6.2.10)；否则，执行步骤(6.2.5)；步骤(6.2.5)、根据下式对种群中的每个个体执行变异操作：

$X_k^{*}=X_{a1}^{\mathrm{iter}}+F\×(X_{a2}^{\mathrm{iter}}-X_{a3}^{\mathrm{iter}}),$

其中，

表示第k个个体的变异操作结果，a1，a2，a3是与k不相同的自然数，1≤a1，a2，a3≤NP，差分进化算法中F＝1，在改进差分进化算法中F＝FF_max-[(FF_max-FF_min)×iter]/maxiter；

步骤(6.2.6)、根据下式对种群中的每个个体的每维分量执行交叉操作，其中第k个个体的第r维分量的交叉操作为：

其中，

表示步骤(6.2.5)中第k个个体的变异操作结果的第r维分量，表示第k个个体的第r维分量的交叉操作结果，差分进化算法中C＝0.5，在改进差分进化算法中，C＝CC_min+[(CC_max-CC_min)×iter]/maxiter；

步骤(6.2.7)、根据下式对步骤(6.2.6)中的每个新生个体的模板参数的每个分量进行边界处理，其中第k个个体的第r维分量的边界处理为：

若 $x_{\mathrm{kr}}^{\mathrm{\Δ}}>x_{\mathrm{kr}}^{\max }$

则 $x_{\mathrm{kr}}^{\mathrm{\Δ}}=x_{\mathrm{kr}}^{\min }+2.0\×\frac{p}{q}\×(x_{\mathrm{kr}}^{\max }-x_{\mathrm{kr}}^{\min }),$

若 $x_{\mathrm{kr}}^{\mathrm{\&Delta;}}<x_{\mathrm{kr}}^{\min }$

则 $x_{\mathrm{kr}}^{\mathrm{\&Delta;}}=x_{\mathrm{kr}}^{\min }+2.0\&times;\frac{q}{p}\&times;(x_{\mathrm{kr}}^{\max }-x_{\mathrm{kr}}^{\min })$

其中， $p=x_{\mathrm{kr}}^{\mathrm{\&Delta;}}-x_{\mathrm{kr}}^{\max },$ $q=x_{\mathrm{kr}}^{\mathrm{\&Delta;}}-x_{\mathrm{kr}}^{\min };$

步骤(6.2.8)、根据下式选择下一次迭代中要保留的个体，其中第k个个体为：

其中，表示步骤(6.2.7)中得到的第k个个体，

表示第iter次迭代的第k个个体，

表示第k个个体

的适应度值；

步骤(6.2.9)、iter＝iter+1，迭代步数加1，返回到步骤(6.2.3)，继续执行下一次迭代操作；

步骤(6.2.10)、如果优化结果不满足精度要求，则增加L_s的值L_s＝L_s+QQ，如果L_s大于

执行步骤(6.2.11)，否则回到步骤(6.2.2)；

步骤(6.2.11)：根据最终的模板参数优化结果，利用D2式生成动态软测量模型中每个子系统的脉冲响应参数；

步骤(7)、根据步骤(6)中得到的脉冲响应参数，建立如下动态软测量模型进行在线应用估计，第t时刻的输出估计值为：

$D6:y\left(t\right)=\underset{s=1}{\overset{n}{\mathrm{\&Sigma;}}}\left\{\underset{i=0}{\overset{L_r}{\mathrm{\&Sigma;}}}{h}_s\left(i\right)u_s(t-i)\right\}.$

子系统结构类型为一阶系统时，D2与D3分别表示为：

$D2:\left\{\begin{array}{c}{h}_s\left(i\right)=\frac{K^s}{T^s}\exp \left(\frac{-1}{T^s}(i-{\mathrm{\&tau;}}^s)\right),i={\mathrm{\&tau;}}^s,...,L_s;s\&Element;[1,n]\\ h_s\left(i\right)=0,i<{\mathrm{\&tau;}}^s\end{array}\right.,$

D3： ${\mathrm{\&zeta;}}_{s1}^{\min }\&le;T^s\&le;{\mathrm{\&zeta;}}_{s1}^{\max },$ ${\mathrm{\&zeta;}}_{s2}^{\min }\&le;K^s\&le;{\mathrm{\&zeta;}}_{s2}^{\max },$ ${\mathrm{\&zeta;}}_{s3}^{\min }\&le;{\mathrm{\&tau;}}^s\&le;{\mathrm{\&zeta;}}_{s3}^{\max },$ s∈[1，n]，

其中，T^s，K^s，τ^s为第s个子系统结构类型的模板参数，分别表示时间常数、比例系数、纯时延常数。模板参数范围通常为：

${\mathrm{\&zeta;}}_{s1}^{\min }=(0.01~0.1)T_1^s;$ ${\mathrm{\&zeta;}}_{s1}^{\max }=(5~8)T_1^s;$ ${\mathrm{\&zeta;}}_{s2}^{\min }=(0.05~0.15)K^s,$ ${\mathrm{\&zeta;}}_{s2}^{\max }=(6~10)K^s,$

${\mathrm{\&zeta;}}_{s3}^{\min }=(0.01~0.1){\mathrm{\&tau;}}^s,$ ${\mathrm{\&zeta;}}_{s3}^{\max }=(2~5){\mathrm{\&tau;}}^s,$

其中，

分别表示第s个子系统的时间常数的最小以及最大值，

分别表示第s个子系统的比例系数的最小以及最大值，

分别表示第s个子系统的纯时延常数的最小以及最大值。

子系统结构类型为二阶系统时，D2与D3分别表示为：

$D2:\left\{\begin{array}{c}{h}_s\left(i\right)=\frac{K_p^s}{T_1^s-T_2^s}\{\exp \left(\frac{-1}{T_1^s}(i-{\mathrm{\&tau;}}_1^s)\right)-\exp \left(\frac{-1}{T_2^s}\right(i-{\mathrm{\&tau;}}_1^s\left)\right)\},i={\mathrm{\&tau;}}_1^s+1,...,L_s;s\&Element;[1,n]\\ h_s\left(i\right)=0,i\&le;{\mathrm{\&tau;}}_1^s\end{array}\right.,$

D3： ${\mathrm{\&xi;}}_{s1}^{\min }\&le;T^s\&le;{\mathrm{\&xi;}}_{s1}^{\max },$ ${\mathrm{\&xi;}}_{s2}^{\min }\&le;T_2^s\&le;{\mathrm{\&xi;}}_{s2}^{\max },$ ${\mathrm{\&xi;}}_{s3}^{\min }\&le;K_p^s\&le;{\mathrm{\&xi;}}_{s3}^{\max },$ ${\mathrm{\&xi;}}_{s4}^{\min }\&le;{\mathrm{\&tau;}}_1^s\&le;{\mathrm{\&xi;}}_{s4}^{\max },$ s∈[1，n]

其中，

为第s个子系统结构类型的模板参数，分别表示时间常数1、时间常数2、比例系数、纯时延常数。其中模板参数范围通常为：

${\mathrm{\&xi;}}_{s1}^{\min }=(0.01~0.1)T_1^s,$ ${\mathrm{\&xi;}}_{s1}^{\max }=(5~8)T_1^s;$ ${\mathrm{\&xi;}}_{s2}^{\min }=(0.01~0.1)T_2^s,$ ${\mathrm{\&xi;}}_{s2}^{\max }=(5~8)T_2^s;$

${\mathrm{\&xi;}}_{s3}^{\min }=(0.05~0.15)K_p^s,$ ${\mathrm{\&xi;}}_{s3}^{\max }=(6~10)K_p^s;$ ${\mathrm{\&xi;}}_{s4}^{\min }=(0.01~0.1){\mathrm{\&tau;}}_1^s,$ ${\mathrm{\&xi;}}_{s4}^{\max }=(2~5){\mathrm{\&tau;}}_1^s,$

其中，

分别表示第s个子系统的时间常数1的最小以及最大值，

分别表示第s个子系统的的时间常数2的最小以及最大值，

分别表示第s个子系统的的比例系数的最小以及最大值，

分别表示第s个子系统的纯时延常数的最小以及最大值。

本发明涉及基于脉冲响应模板和参数优化的动态软测量建模方法，其优点在于：

1)该方法是一种动态软测量方法，适应动态变化生产过程下的软测量建模和应用的需求，可以比静态软测量方法应用在动态生产过程的估计精度有明显改进。

2)该方法对难测主导变量的采样无特殊限制(例如无需采样周期很短、等周期采样或大量的采样样本等)，这样才能适用于难测主导变量往往是人工取样化验情况下的动态软测量建模，对易测辅助变量也仅需要相对于难测主导变量采样时及前一个过渡过程内的相应采样数据，解决了现有的以神经网络模型为代表的的基于实验模型的动态软测量方法只能用于连续等周期采样的问题。

3)该方法由于利用了脉冲响应模板和优化方法，使得脉冲响应参数辨识限定在指定的可以根据机理分析和经验确定的模板类和模板参数约束范围内，有效改进了具有噪声情况下模型的过学习导致的泛化性及鲁棒性差的问题和降低了模型学习时间，尤其对于多变量系统，大量脉冲响应参数的模型参数学习转化为对有限个模板参数学习，极大减少了学习参数以节省学习时间。

附图说明

图1.生产过程系统的软测量采样原理图，其中易测辅助变量在线采样，采样周期应满足分钟级，通常为1分钟；难测主导变量离线采样，采样周期可较长，通常为小时级，例如6～8小时。

图2.基于脉冲响应模板和参数优化的动态软测量建模方法原理图。

图3.基于脉冲响应模板和参数优化的动态软测量建模方法主程序图。

图4.仿真实例测试结果比较。上图为脉冲响应参数估计与真实值的比较；下图为难测变量估计与真实值的比较。

具体实施方式

下面结合一个仿真实例对本发明做进一步说明。

为接近实际生产过程，选定2输入—单输出生成包含两个子系统的过程系统，子系统为为二阶惯性加滞后环节，各子系统的模型参数为：

${\mathrm{\&Omega;}}_1=\{K_p^1,T_1^1,T_2^1,{\mathrm{\&tau;}}_1^1\}=\left\{\mathrm{120,15,12,2}\right\},$ ${\mathrm{\&Omega;}}_2=\{K_p^2,T_1^2,T_2^2,{\mathrm{\&tau;}}_1^2\}=\left\{\mathrm{30,5,10,5}\right\};$

用Matlab仿真软件模拟训练数据样本：两个易测辅助变量的输入类型均为随机均匀分布，模拟采样周期为一分钟，分别产生两组48000个点的序列(模拟400小时的采样数据)，基于上述给定的子系统的模型参数，利用两组输入产生两组对应的过程子系统输出序列，并将对应序列点进行求和同时外加20％的随机噪声，得到48000个点的系统总输出序列；对总输出序列进行稀疏采样，采样周期为P＝480分钟(8小时)，产生M＝100个点作为模型输出o(P_j)，用于模拟离线采样数据；简单起见，规定两个子系统的脉冲响应长度初值L₁＝L₂＝100，对应的上限值为 $L_1^{\max },L_2^{\max }=200;$ 根据说明书中步骤(2)产生所需的训练样本数据。

脉冲响应模板转化函数为：

$\mathrm{\&Xi;}\left\{\begin{array}{c}{h}_s\left(i\right)=\frac{K_p^s}{T_1^s-T_2^s}\{\exp \left(\frac{-1}{T_1^s}(i-{\mathrm{\&tau;}}_1^s)\right)-\exp \left(\frac{-1}{T_2^s}\right(i-{\mathrm{\&tau;}}_1^s\left)\right)\},i={\mathrm{\&tau;}}_1^s+1,...,L_s;s\&Element;[1,2]\\ h_s\left(i\right)=0,i\&le;{\mathrm{\&tau;}}_1^s\end{array}\right.,$

其中，模板参数为：

$x_{s1}^{\min }\&le;T_1^s\&le;x_{s1}^{\max },$ $x_{s2}^{\min }\&le;T_2^s\&le;x_{s2}^{\max },$ $x_{s3}^{\min }\&le;K_p^s\&le;x_{s3}^{\max },$ $x_{s4}^{\min }\&le;{\mathrm{\&tau;}}_1^s\&le;x_{s4}^{\max },$ s∈[1，2]

模板参数范围为：

$\mathrm{\&Theta;}:\left\{\begin{array}{c}{x}_{s1}^{\min }=0.01T_1^s,x_{s1}^{\max }=8T_1^s;x_{s2}^{\min }=0.01T_2^s,x_{s2}^{\max }=8T_2^s;\\ x_{s3}^{\min }=0.05K_p^s,x_{s3}^{\max }=10K_p^s;x_{s4}^{\min }=0.01{\mathrm{\&tau;}}_{s4}^{\min }=0.01{\mathrm{\&tau;}}_1^s,x_{s4}^{\max }=5{\mathrm{\&tau;}}_1^s\end{array}\right.$

差分进化算法与改进差分进化算法的参数设置：种群规模NP＝30；最大迭代步数maxiter＝500；目标函数值ε＝0.05；优化参数矢量约束范围为Θ。差分进化算法中缩放比例F＝1.0，交叉率C＝0.5；改进差分进化算法中缩放比例FF_min＝0.1，FF_max＝1.8，交叉率CC_min＝0.1，CC_max＝0.8。

训练步骤如下：

步骤A1、初始化

设定：差分进化算法与改进差分进化算法的参数设置：种群规模NP＝30；最大迭代步数maxiter＝500；优化算法精度值ε＝0.05；优化参数矢量约束范围为Θ。差分进化算法中缩放比例F＝1.0，交叉率C＝0.5；改进差分进化算法中缩放比例FF_min＝0.1，FF_max＝1.8，交叉率CC_min＝0.1，CC_max＝0.8；初值L1＝L2＝100，对应的上限值为QQ为L₁，L₂的递增步长，QQ＝10；

步骤A2、群体即种群规模初始化：

2个子系统中的模板参数的总数构成一个8维的个体矢量，初始种群规模由30个8维的个体矢量组成，其中第k个个体矢量X_k是记为

迭代步数iter＝0：

根据Θ生成初始群体：

$x_{\mathrm{kr}}^{\mathrm{iter}}=x_{\mathrm{kr}}^{\min }+\mathrm{rand}\left(\mathrm{0,1}\right)\&times;(x_{\mathrm{kr}}^{\max }-x_{\mathrm{kr}}^{\min });,$

其中，

代表第k个个体矢量的第r个分量所对应模板参数的最小值，

为其最大值，表示第iter步时第k个个体矢量的第r个模板参数的计算值；

步骤A3、按说明书步骤(6.1)中D5式计算每个个体的适应度值，其中第k个个体X_k的适应度值

确定种群中最大的适应度值 $J_{\max }^{\mathrm{iter}}=\max \left(J_k^{\mathrm{iter}}\right),$ k＝1，...，30，并记录最大的适应度值个体为gBest；

步骤A4、判断

是否满足精度要求，或者是否达到最大迭代步数，若条件 $J_{\max }^{\mathrm{iter}}>\mathrm{\&epsiv;}$ 或iter>maxiter满足，则输出个体gBest，并执行步骤A10；否则，执行步骤A5；

步骤A5、根据下式对种群中的每个个体执行变异操作：

$X_k^{*}=X_{a1}^{\mathrm{iter}}+F\&times;(X_{a2}^{\mathrm{iter}}-X_{a3}^{\mathrm{iter}}),$

其中，

表示k第个个体的变异操作结果，a1，a2，a3是与k不相同的自然数，1≤a1，a2，a3≤30，在改进差分进化算法中F＝1.8-(1.7×iter)/500；

步骤A6、根据下式对种群中的每个个体的每维分量执行交叉操作，其中第k个个体的第r维分量的交叉操作为：

其中，表示步骤A5中第k个个体的变异操作结果的第r维分量，

表示第k个个体的第r维分量的交叉操作结果，在改进差分进化算法中，C＝0.1+(0.7×iter)/500；

步骤A7、根据下式对步骤A6中的每个新生个体的模板参数的每个分量进行边界处理，其中第k个个体的第r维分量的边界处理为：：

若 $x_{\mathrm{kr}}^{\mathrm{\&Delta;}}>x_{\mathrm{kr}}^{\max }$

则 $x_{\mathrm{kr}}^{\mathrm{\&Delta;}}=x_{\mathrm{kr}}^{\min }+2.0\&times;\frac{p}{q}\&times;(x_{\mathrm{kr}}^{\max }-x_{\mathrm{kr}}^{\min })$

若 $x_{\mathrm{kr}}^{\mathrm{\&Delta;}}<x_{\mathrm{kr}}^{\min }$

则 $x_{\mathrm{kr}}^{\mathrm{\&Delta;}}=x_{\mathrm{kr}}^{\min }+2.0\&times;\frac{q}{p}\&times;(x_{\mathrm{kr}}^{\max }-x_{\mathrm{kr}}^{\min })$

其中， $p=x_{\mathrm{kr}}^{\mathrm{\&Delta;}}-x_{\mathrm{kr}}^{\max };$ $q=x_{\mathrm{kr}}^{\mathrm{\&Delta;}}-x_{\mathrm{kr}}^{\min };$

步骤A8、根据下式选择下一次迭代中要保留的个体，其中第k个个体为：

其中，表示步骤A7中得到的第k个个体，

表示第iter次迭代的第k个个体，表示第k个个体

的适应度值；

步骤A9、iter＝iter+1，迭代步数加1，返回到步骤A3，继续执行下一次迭代操作；

步骤A10、如果优化结果不满足精度要求，则增加L₁，L₂的值L₁＝L₁+10，L₂＝L₂+10，如果L₁，或L₂大于200，执行步骤A11，否则回到步骤A2；

步骤A11：根据最终的模板参数优化结果，利用

式生成动态软测量模型中每个子系统的脉冲响应参数；

在本例中当L₁＝L₂＝150时，模型训练精度最高。重新随机生成200组数据样本，进行估计精度与泛化性能测试。

测试动态软测量模型为：

$y\left(k\right)=\underset{s=1}{\overset{2}{\mathrm{\&Sigma;}}}\left\{\underset{i=0}{\overset{150}{\mathrm{\&Sigma;}}}{h}_s\left(i\right)u_s(k-i)\right\}$

误差衡量指标为均方误差衡量：

$\mathrm{MSE}=\sqrt{\frac{1}{200}\underset{i=1}{\overset{200}{\mathrm{\&Sigma;}}}{[o\left(i\right)-y\left(i\right)]}^2}$

其中y(i)，o(i)分别表示模型估计值以及真实输出值。

当噪声为20％时，仿真实例测试结果如图4所示。上图为脉冲响应参数估计与真实值的比较。下图为难测变量估计与真实值的比较，其中，基于脉冲响应模板与改进差分进化优化方法的动态软测量方法(MDE-IRP-DSSM)的测试MSE为1.4725；基于脉冲响应模板与差分进化优化方法的动态软测量方法(DE-IRP-DSSM)的测试MSE为1.6100。从测试结果来看，本发明方法在大噪声环境下估计精度较高，可以满足生产过程的要求。

Claims (3)

1.基于脉冲响应模板和参数优化的动态软测量建模方法，其特征在于，该方法是在计算机中依次按以下步骤进行：

步骤(1)、初始化

在生产过程工艺机理分析或实际经验的基础上构建一个难测主导变量o和n个易测辅助变量u₁，u₂，...，u_n构成的n个独立子系统，n≥1，该各个子系统的结构类型为下述两个结构类型之中的任何一种：

一阶惯性环节加纯延迟： $G\left(s\right)=\frac{{\mathrm{Ke}}^{-\mathrm{\&tau;s}}}{\mathrm{Ts}+1}$

二阶惯性环节加纯延迟： $G\left(s\right)=\frac{K_p{e}^{-{\mathrm{\&tau;}}_{1}s}}{(T_{1}s+1)(T_{2}s+1)}$

其中，K，K_p为比例系数，

T，T₁，T₂为时间常数，

τ，τ₁为纯时延常数，

K，K_p，T，T₁，T₂，τ，τ₁均为待估的模型参数，每个模型参数的估计区间在生产过程工艺机理分析或实际经验的基础上给定；

步骤(2)、从在线采集的数据库以及离线化验数据库中分析搜集相关的易测辅助变量u₁，u₂，..，u_n与难测主导变量o构成M组数据样本用于动态软测量模型的训练，并按采样时间递增的顺序形成下述相应的时间序列：

$[u_1\left(Q_{1j}\right),...,u_1\left(P_j\right);u_2\left(Q_{2j}\right),...,u_2\left(P_j\right);...;u_n\left(Q_{\mathrm{nj}}\right),...,u_n\left(P_j\right)]\&DoubleLeftRightArrow;o\left(P_j\right),j=1,...,M$

其中，左边表示系统输入，右边表示系统输出，

M为总的采样次数，代表样本数，

第j组数据样本： $[u_1\left(Q_{1j}\right),...,u_1\left(P_j\right);u_2\left(Q_{2j}\right),...,u_2\left(P_j\right);...;u_n\left(Q_{\mathrm{nj}}\right),...,u_n\left(P_j\right)]\&DoubleLeftRightArrow;o\left(P_j\right),$

Q_1j，..，P_j表示易测辅助变量u₁的第j个样本的采样时刻，Q_1j与P_j之间形成等差数列，u₁(Q_1j)，...，u₁(P_j)包含L₁个采样数据并作为第1个子系统的输入；Q_2j，..，P_j表示易测辅助变量u₂的第j个样本的采样时刻，Q_2j与P_j之间形成等差数列，u₂(Q_2j)，...，u₂(P_j)包含L₂个采样数据并作为第2个子系统的输入；依此类推，Q_nj，...，P_j表示易测辅助变量u_n的第j个样本的采样时刻，Q_nj与P_j之间形成等差数列，u_n(Q_nj)，...，u_n(P_j)包含L_n个采样数据并作为第n个子系统的输入；P_j表示难测主导变量o的第j个样本的采样时刻，o(P_j)为P_j时刻的采样值并作为总系统的输出，

用u_s(Q_sj)，...，u_s(P_j)表示第s个子系统的输入，s∈[1，n]，该子系统的L_s个采样数据的采样总时间不短于易测辅助变量u_s与难测主导变量o之间的一个动态过渡过程时间；

步骤(3)、对步骤(2)得到的样本数据进行坏点剔除以及标准化处理，其步骤如下：

步骤(3.1)、用相邻点绝对误差阈值法剔除坏点：从样本的第二个点开始，依次计算各点与相邻点的绝对差，若该差值超过某个设定的阈值，则剔除该点，并用被剔除点位置前后两个样本点的平均值代替；

步骤(3.2)、按下式对步骤(3.1)得到的样本数据进行标准化处理：

x^o＝(x-m_x)/σ_x，

其中，x为相应样本原始数据，x^o为相应样本对应的标准化数据，m_x、σ_x分别为该组数据样本的算术平均

值和标准差；

步骤(4)、对步骤(3)得到的各组数据样本，基于脉冲响应原理对每个子系统的M组数据样本建立卷积模型，再把各个子系统的卷积模型求和，得到一个全系统的脉冲响应模型D1：

$D1:o\left(P_j\right)=\underset{s=1}{\overset{n}{\mathrm{\&Sigma;}}}\left\{\underset{i=0}{\overset{L_s}{\mathrm{\&Sigma;}}}{h}_s\left(i\right)u_s(P_j-i)\right\},j=1,...,M,$

其中，n为子系统的总数，

L_s为第s个子系统的脉冲响应的长度参数，

h_s(i)为待估的第s个子系统度的脉冲响应参数；

步骤(5)、根据步骤(1)给出的各子系统的结构类型，用Ψ(Ω_s)产生一组对应的脉冲响应D2及其中的参数范围D3：

D2：h_s(i)＝Ψ(Ω_s)，s＝1，...，n，i＝1，2，...，L_s

其中，Ω_s为代表子系统结构类型的模板参数：K，或K_p，或T，或T₁，或T₂，或τ，或τ₁，

Ψ(·)为脉冲响应函数，

D3：Ω_s-Π^low≤Ω_s≤Ω_s+Π^upp，

其中，Π^low为最小阈值，Π^low>0，设定值，

Π^upp为最大阈值，Π^upp>0，设定值，

L_s由香农采样确定，其取值范围D4为：

D4：L_s-a≤L_s≤L_s+b，s＝1，...，n；20≤a≤50；30≤b≤60；

步骤(6)、利用差分进化算法与改进差分进化算法搜索最佳模板参数，实现动态软测量模型优化求解，其步骤如下：

步骤(6.1)、优化目标函数D5表示为：

$D5:\max :J=1/\underset{j=1}{\overset{M}{\mathrm{\&Sigma;}}}{(y\left(P_j\right)-o\left(P_j\right))}^2,$

其中， $y\left(P_j\right)=\underset{s=1}{\overset{n}{\mathrm{\&Sigma;}}}\left\{\underset{i=0}{\overset{L_s}{\mathrm{\&Sigma;}}}\stackrel{\&OverBar;}{h_s}\left(i\right)u_s(P_j-i)\right\},$ j＝1，...，M，表示优化算法计算的输出，

$\stackrel{\&OverBar;}{h_s}\left(i\right)=\mathrm{\&Psi;}\left(\stackrel{\&OverBar;}{{\mathrm{\&Omega;}}_s}\right),\stackrel{\&OverBar;}{{\mathrm{\&Omega;}}_s}$ 为优化算法给出的子系统结构类型的模板参数，模板参数指的是用D3优化时用的模型参数，

o(P_j)为对应的真实模型输出，由步骤(2)得到；

步骤(6.2)、按以下步骤进行优化计算：

步骤(6.2.1)、初始化

设定：用NP表示种群规模，30≤NP≤60，maxiter表示最大迭代步数，300≤maxiter≤1500，ε表示优化算法精度值，0.01≤ε≤0.1，D3表示优化参数矢量约束范围，F表示差分进化算法中缩放比例，F＝1，C表示交叉率，C＝0.5，FF_min、FF_max表示改进差分进化算法中最小缩放比例以及最大缩放比例，FF_min＝0.1，FF_max＝1.8，CC_min、CC_max表示改进差分进化算法中最小交叉率以及最大交叉率，CC_min＝0.1，CC_max＝0.8，根据步骤(5)设定L_s初值为L_s-a，L_s的取值上限

为L_s+b，QQ为L_s的递增步长，为5到10之间的一个整数；

步骤(6.2.2)、群体即种群规模初始化：

n个子系统中的模板参数的总数构成一个v维的个体矢量，初始种群规模由NP个v维的个体矢量组成，其中第k个个体矢量X_k是记为

$\{x_{k1},x_{k2},...,x_{\mathrm{kv}}\}=$

的数据，

其中，v为一个个体矢量的维数，即模板参数的总数，K，T，τ，K_p，T₁，T₂，τ₁的上标表示子系统序号，若子系统为一阶系统时，子系统的模型参数为K，T，τ，若子系统为二阶系统时，子系统的模型参数为K_p，T₁，T₂，τ₁，

迭代步数iter＝0，

根据下式生成初始群体：

$x_{\mathrm{kr}}^{\mathrm{iter}}=x_{\mathrm{kr}}^{\min }+\mathrm{rand}\left(\mathrm{0,1}\right)\&times;(x_{\mathrm{kr}}^{\max }-x_{\mathrm{kr}}^{\min });,$

其中，

代表第k个个体矢量的第r个分量所对应模板参数的最小值，

为其最大值，表示第iter步时第k个个体矢量的第r个分量的计算值，rand(0，1)表示0与1之间的随机数；

步骤(6.2.3)、按步骤(6.1)中的D5式计算每个个体的适应度值，其中第k个个体X_k的适应度值

确定种群中最大的适应度值 $J_{\max }^{\mathrm{iter}}=\max \left(J_k^{\mathrm{iter}}\right),$ k＝1，...，NP，并记录最大的适应度值个体为gBest；

步骤(6.2.4)、判断

是否满足精度要求，或者是否达到最大迭代步数，若条件 $J_{\max }^{\mathrm{iter}}>\mathrm{\&epsiv;}$ 或iter>maxiter满足，则输出个体gBest，并执行步骤(6.2.10)；否则，执行步骤(6.2.5)；

步骤(6.2.5)、根据下式对种群中的每个个体执行变异操作：

$X_k^{*}=X_{a1}^{\mathrm{iter}}+F\&times;(X_{a2}^{\mathrm{iter}}-X_{a3}^{\mathrm{iter}}),$

其中，

表示第k个个体的变异操作结果，a1，a2，a3是与k不相同的自然数，1≤a1，a2，a3≤NP，差分进化算法中F＝1，在改进差分进化算法中F＝FF_max-[(FF_max-FF_min)×iter]/maxiter；

步骤(6.2.6)、根据下式对种群中的每个个体的每维分量执行交叉操作，其中第k个个体的第r维分量的交叉操作为：

其中，

表示步骤(6.2.5)中第k个个体的变异操作结果的第r维分量，

表示第k个个体的第r维分量的交叉操作结果，差分进化算法中C＝0.5，在改进差分进化算法中，C＝CC_min+[(CC_max-CC_min)×iter]/maxiter；

步骤(6.2.7)、根据下式对步骤(6.2.6)中的每个新生个体的模板参数的每个分量进行边界处理，其中第k个个体的第r维分量的边界处理为：

若 $x_{\mathrm{kr}}^{\mathrm{\&Delta;}}>x_{\mathrm{kr}}^{\max }$

则 $x_{\mathrm{kr}}^{\mathrm{\&Delta;}}=x_{\mathrm{kr}}^{\min }+2.0\&times;\frac{p}{q}\&times;(x_{\mathrm{kr}}^{\max }-x_{\mathrm{kr}}^{\min })$

若 $x_{\mathrm{kr}}^{\mathrm{\&Delta;}}<x_{\mathrm{kr}}^{\min }$

则 $x_{\mathrm{kr}}^{\mathrm{\&Delta;}}=x_{\mathrm{kr}}^{\min }+2.0\&times;\frac{q}{p}\&times;(x_{\mathrm{kr}}^{\max }-x_{\mathrm{kr}}^{\min })$

其中， $p=x_{\mathrm{kr}}^{\mathrm{\&Delta;}}-x_{\mathrm{kr}}^{\max },$ $q=x_{\mathrm{kr}}^{\mathrm{\&Delta;}}-x_{\mathrm{kr}}^{\min };$

步骤(6.2.8)、根据下式选择下一次迭代中要保留的个体，其中第k个个体为：

其中，

表示步骤(6.2.7)中得到的第k个个体，

表示第iter次迭代的第k个个体，

表示第k个个体

的适应度值；

步骤(6.2.9)、iter＝iter+1，迭代步数加1，返回到步骤(6.2.3)，继续执行下一次迭代操作；

步骤(6.2.10)、如果优化结果不满足精度要求，则增加L_s的值L_s＝L_s+QQ，如果L_s大于

执行步骤(6.2.11)，否则回到步骤(6.2.2)；

步骤(6.2.11)：根据最终的模板参数优化结果，利用D2式生成动态软测量模型中每个子系统的脉冲响应参数；

步骤(7)、根据步骤(6)中得到的脉冲响应参数，建立如下动态软测量模型进行在线应用估计，第t时刻的输出估计值为：

$D6:y\left(t\right)=\underset{s=1}{\overset{n}{\mathrm{\&Sigma;}}}\left\{\underset{i=0}{\overset{L_s}{\mathrm{\&Sigma;}}}{h}_s\left(i\right)u_s(t-i)\right\}.$

2.根据权利要求1所述的基于脉冲响应模板和参数优化的动态软测量建模方法，其特征在于，

子系统结构类型为一阶系统时，D2与D3分别表示为：

$D2:\left\{\begin{array}{c}{h}_s\left(i\right)=\frac{K^s}{T^s}\exp \left(\frac{-1}{T^s}(i-{\mathrm{\&tau;}}^s)\right),i={\mathrm{\&tau;}}^s,...,L_s;s\&Element;[1,n]\\ h_s\left(i\right)=0,i<{\mathrm{\&tau;}}^s\end{array}\right.,$

$D3:{\mathrm{\&zeta;}}_{s1}^{\min }\&le;T^s\&le;{\mathrm{\&zeta;}}_{s1}^{\max },$ ${\mathrm{\&zeta;}}_{s2}^{\min }\&le;K^s\&le;{\mathrm{\&zeta;}}_{s2}^{\max },$ ${\mathrm{\&zeta;}}_{s3}^{\min }\&le;{\mathrm{\&tau;}}^s\&le;{\mathrm{\&zeta;}}_{s3}^{\max },$ s∈[1，n]，

其中，T^s，K^s，τ^s为第s个子系统结构类型的模板参数，分别表示时间常数、比例系数、纯时延常数；模板参数范围通常为：

${\mathrm{\&zeta;}}_{s1}^{\min }=(0.01~0.1)T_1^s;$ ${\mathrm{\&zeta;}}_{s1}^{\max }=(5~8)T_1^s;$ ${\mathrm{\&zeta;}}_{s2}^{\min }=(0.05~0.15)K^s,$ ${\mathrm{\&zeta;}}_{s2}^{\max }=(6~10)K^s,$

${\mathrm{\&zeta;}}_{s3}^{\min }=(0.01~0.1){\mathrm{\&tau;}}^s,$ ${\mathrm{\&zeta;}}_{s3}^{\max }=(2~5){\mathrm{\&tau;}}^s,$

其中，分别表示第s个子系统的时间常数的最小以及最大值，分别表示第s个子系统的比例系数的最小以及最大值，

分别表示第s个子系统的纯时延常数的最小以及最大值；

子系统结构类型为二阶系统时，D2与D3分别表示为：

$D2:\left\{\begin{array}{c}{h}_s\left(i\right)=\frac{K_p^s}{T_1^s-T_2^s}\{\exp \left(\frac{-1}{T_1^s}(i-{\mathrm{\&tau;}}_1^s)\right)-\exp \left(\frac{-1}{T_2^s}(i-{\mathrm{\&tau;}}_1^s)\right)\},i={\mathrm{\&tau;}}_1^s+1,...,L_s;s\&Element;[1,n]\\ h_s\left(i\right)=0,i<{\mathrm{\&tau;}}_1^s\end{array}\right.,$

$D3:{\mathrm{\&xi;}}_{s1}^{\min }\&le;T_1^s\&le;{\mathrm{\&xi;}}_{s1}^{\max },$ ${\mathrm{\&xi;}}_{s2}^{\min }\&le;T_2^s\&le;{\mathrm{\&xi;}}_{s2}^{\max },$ ${\mathrm{\&xi;}}_{s3}^{\min }\&le;K_p^s\&le;{\mathrm{\&xi;}}_{s3}^{\max },$ ${\mathrm{\&xi;}}_{s4}^{\min }\&le;{\mathrm{\&tau;}}_1^s\&le;{\mathrm{\&xi;}}_{s4}^{\max },$ s∈[1，n]

其中，为第s个子系统结构类型的模板参数，分别表示时间常数1、时间常数2、比例系数、纯时延常数；其中模板参数范围通常为：

${\mathrm{\&xi;}}_{s1}^{\min }=(0.01~0.1)T_1^s,$ ${\mathrm{\&xi;}}_{s1}^{\max }=(5~8)T_1^s;$ ${\mathrm{\&xi;}}_{s2}^{\min }=(0.01~0.1)T_2^s,$ ${\mathrm{\&xi;}}_{s2}^{\max }=(5~8)T_2^s;$

${\mathrm{\&xi;}}_{s3}^{\min }=(0.05~0.15)K_p^s,$ ${\mathrm{\&xi;}}_{s3}^{\max }=(6~10)K_p^s;$ ${\mathrm{\&xi;}}_{s4}^{\min }=(0.01~0.1){\mathrm{\&tau;}}_1^s,$ ${\mathrm{\&xi;}}_{s4}^{\max }=(2~5){\mathrm{\&tau;}}_1^s,$

其中，

分别表示第s个子系统的时间常数1的最小以及最大值，分别表示第s个子系统的的时间常数2的最小以及最大值，分别表示第s个子系统的的比例系数的最小以及最大值，

分别表示第s个子系统的纯时延常数的最小以及最大值。

3.根据权利要求1所述的基于脉冲响应模板和参数优化的动态软测量建模方法，其特征在于，步骤(6)中的差分进化算法或改进差分进化算法还采用粒子群算法PSO、遗传算法GA优化方法代替。

Images