CN104459089B - 一种高浓磨浆系统游离度的软测量方法 - Google Patents

Abstract

一种高浓磨浆系统游离度的软测量方法涉及化学机械制浆生产过程高浓磨浆系统游离度的测量技术领域，具体地说是涉及一种基于粒子群算法优化ε-支持向量回归机(PSO-ε-SVR)的高浓磨浆系统游离度的软测量方法。本发明提供一种高浓磨浆系统游离度的软测量方法，该方法能够实现对高浓磨浆机出浆游离度的实时测量，提高人工测量的效率和精度，为工业实时控制和实现节能降耗提供技术支持。本发明方法由硬件平台及测量软件组成，其中硬件平台核心包括高浓磨浆系统三盘磨浆机、压力传感器、震动传感器、加速度传感器、流量传感器、分料螺旋等。

Description

一种高浓磨浆系统游离度的软测量方法

技术领域

本发明涉及化学机械制浆生产过程高浓磨浆系统游离度的测量技术领域，具体地说是涉及一种基于粒子群算法优化ε-支持向量回归机(PSO-ε-SVR)的高浓磨浆系统游离度的软测量方法。

背景技术

造纸过程由制浆和造纸两大环节串联而成。制浆环节的主要功能是从植物纤维原材料生产出具有特定形态的纤维，而造纸环节的功能主要是以特定形态的纤维为原料生产出各种纸制产品。制浆造纸工业在我国是重点能耗行业之一，生产过程中消耗大量热能、电能和水能。特别是制浆过程的能耗及其所产生的纸浆的质量直接关系到后续造纸环节的能耗及产品质量。

按现有工艺，制浆分为机械浆、化学浆和化学机械浆三大类。其中，化学机械浆生产过程主要包括：木片预处理，高浓磨浆，粗浆漂白，粗浆洗涤，低浓磨浆，渣浆磨浆，纤维回收，浆料浓缩脱水等过程。高浓磨浆作为化学机械浆生产的关键步骤，其用电量占到整个工厂用电的35％以上，决定着纸浆质量和生产能耗。

然而，高浓磨浆过程具有多变量、强耦合和非线性的特征，使得高浓磨浆过程的机理分析、建模存在很大困难。目前，研究的假设性强，得到的高浓磨浆模型缺乏通用性。此外，研究主要集中在低浓磨浆过程、单盘磨浆机和磨盘改进上，对于高浓磨浆系统三盘磨浆机的建模和分析还没有公开报道。

纤维网滤水性(游离度)是检测纸浆通过铜网的滤水程度快慢的一个指标。作为评价纸浆质量的关键指标之一，该指标不仅关系着最终纸张的质量(透光性，吸水性等)，同时与整个生产过程的能耗密切相关。为此，使用游离度来评价纤维质量十分合理。

公开资料表明，目前还没有传感器对于高浓磨浆系统游离度直接进行测量，主要依靠离线间歇人工测量。该测量方法存在着较大的时滞性和误差，进而无法实时对制浆生产进行控制，阻碍实现节能降耗和纸浆质量控制。为此，提出一种基于粒子群算法优化ε-支持向量机(PSO-ε-SVR)的高浓磨浆系统游离度软测量的技术。本发明可通过常规在线测量方法提供的模型输入数据，给出输入数据对应的在线估计值，为制浆生产过程的优化操作和优化运行提供关键质量指标，进而实现质量控制和节能降耗。

发明内容

为了克服高浓磨浆过程中游离度无法直接测量，离线人工测量时效性差、精度低等缺点，本发明提供一种高浓磨浆系统游离度的软测量方法，该方法能够实现对高浓磨浆机出浆游离度的实时测量，提高人工测量的效率和精度，为工业实时控制和实现节能降耗提供技术支持。

本发明方法由硬件平台及测量软件组成，其中硬件平台核心包括高浓磨浆系统三盘磨浆机、压力传感器、震动传感器、加速度传感器、流量传感器、分料螺旋等。

一种最优参数的支持向量回归机高浓磨浆系统游离度软测量包括以下步骤：(1)辅助变量选择和模型输入变量确定，(2)PSO-ε-SVR软测量模型的训练和软测量实现。

(1)辅助变量选择和模型输入变量确定

需要软测量的高浓磨浆系统纸浆关键参数为游离度(ml)。根据生产工艺及机理以及可测可观和变量间相关性分析(因子分析、协方差分析)，选择i个软测量辅助变量(i＝9)，包括：

设定产量u₁(adt/d)；

磨浆机功率u₂(Mw)；

传动侧稀释水量u₃(l/min)；

非传动侧稀释水量u₄(l/min)；

非传动侧震动加速度u₅(％)；

非传动侧震动u₆(mm/s)；

非传动侧磨盘间隙u₇(mm)；

传动侧磨盘压力u₈(bar)；

非传动侧磨盘压力u₉(bar)。

确定上述9个软测量辅助变量为软测量模型的输入变量，游离度(ml)为软测量模型的输出变量。

(2)PSO-ε-SVR软测量模型的训练和软测量实现

(A)开始：变量初始化；

(B)选择进行软测量模型训练，转至(C)，读取软测量训练样本；选择高浓磨浆系统游离度软测量，转至(K)，读取软测量模型最优参数(惩罚因子C和核函数参数δ)；

(C)读取软测量训练样本：从数据库导入或输入训练集其中，软测量模型输入：x_i∈Rⁿ，i＝1,2,…,l，软测量模型输出：l为训练集样本数量；

(D)数据预处理：

采用尖峰滤波算法将高浓磨浆系统的尖峰异常数据进行剔除；由于选取的9个辅助变量存在着量纲差，对每个变量分别进行中心标准化处理，即每类样本数据减去对应样本的平均值，然后除以其样本方差，即：

计算均值：

$\stackrel{\‾}{x_i}=\frac{1}{l}\underset{j=1}{\overset{l}{\Σ}}{x}_{ij}---\left(1\right)$

计算方差：

${\mathrm{\σ}}_{x_i}^2=\frac{1}{l-1}\underset{j=1}{\overset{l}{\Σ}}(x_{ij}-{\stackrel{\‾}{x}}_i)---\left(2\right)$

标准化：

${\tilde{x}}_{ij}=\frac{x_{ij}-{\stackrel{\‾}{x}}_i}{{\mathrm{\σ}}_{x_i}}---\left(3\right)$

其中，x_ij为第i类变量的第j个训练样本，l为训练样本数，为第i类变量的均值，为中心标准化后的训练样本，表示第i类变量的方差，表示第i类变量的标准差；

(E)确定ε-SVR软测量模型的初始化参数，包括：惩罚因子C和误差容限值ε，核函数的类型及其核函数参数δ。惩罚因子C取20，误差容限值ε取0.01，核函数取高斯径向基函数K(x,x')＝exp(-||x-x'||²/δ²)，其中核函数参数δ＝1；

(F)确定ε-SVR待优化参数惩罚因子C的搜索范围为[0,100]，核函数参数δ的搜索范围为[0,50]，粒子群个体的数目40，粒子的最大迭代寻优次数600，粒子的惩罚因子C迭代速度的最大值为25，粒子的核函数参数δ迭代速度的最大值为25，粒子的初始位置、初始速度、个体极值以及全局极值随机生成；

(G)ε-SVR模型训练：

ε-SVR计算高浓磨浆系统游离度的回归超平面方程y＝(ω·x)+b，其中ω为其法方向，b为其截距。

设模型样本的第i个目标输出为y_i，则ε-SVR将拟合问题等价转化为如下所述的二次规划问题：

$\begin{array}{cc}\underset{w,b}{\min }& \frac{1}{2}\left|\right|\mathrm{\ω}|{|}^2+C\underset{i=1}{\overset{l}{\mathrm{\Σ}}}({\mathrm{\ξ}}_i+{\mathrm{\ξ}}_i^{*})\\ s.t.& \left(\right(\mathrm{\ω}\·x_i)+b)-y_i\≤\mathrm{\ϵ}+{\mathrm{\ξ}}_i,i=1,...,l;\\ & y_i-\left(\right(\mathrm{\ω}\·x_i)+b)\≤\mathrm{\ϵ}+{\mathrm{\ξ}}_i^{*},i=1,...,l;\\ & {\mathrm{\ξ}}_j^{*}\≥0,i=1,...,l.\end{array}---\left(4\right)$

其中，C为惩罚因子，ε为误差容限值，ξ_i和(简记为)为松弛变量。

求解式(4)中二次规划问题的解后，得到回归超平面方程：

y＝(ω·x)+b(5)

将Lagrange函数引入上述问题，导出其对偶问题：

$\begin{array}{cc}\underset{{\mathrm{\α}}^{(*)}\∈R^{2l}}{\min }& \frac{1}{2}\underset{i,j=1}{\overset{l}{\mathrm{\Σ}}}({\mathrm{\α}}_i^{*}-{\mathrm{\α}}_i)({\mathrm{\α}}_j^{*}-{\mathrm{\α}}_j)(x_i\·x_j)\\ & +\mathrm{\ϵ}\underset{i=1}{\overset{l}{\mathrm{\Σ}}}({\mathrm{\α}}_i^{*}+{\mathrm{\α}}_i)-\underset{i=1}{\overset{l}{\mathrm{\Σ}}}{y}_i({\mathrm{\α}}_i^{*}-{\mathrm{\α}}_i)\\ s.t.& \underset{i=1}{\overset{l}{\mathrm{\Σ}}}({\mathrm{\α}}_i^{*}-{\mathrm{\α}}_i)=0,\\ & 0\≤{\mathrm{\α}}_i^{(*)}\≤C,i=1,...,l.\end{array}---\left(6\right)$

解得Lagrange乘子向量 ${\mathrm{\α}}^{(*)}={({\mathrm{\α}}_1,{\mathrm{\α}}_1^{*},...,{\mathrm{\α}}_l,{\mathrm{\α}}_l^{*})}^T.$

计算b：选择位于开区间(0,C)中的α^(*)的分量α_j或若选到的是α_j，则有：

$b=y_j-\underset{i=1}{\overset{l}{\Σ}}({\mathrm{\α}}_i^{*}-{\mathrm{\α}}_i)(x_i\·x_j)+\mathrm{\ϵ}---\left(7\right)$

若选择的是则有：

$b=y_k-\underset{i=1}{\overset{l}{\Σ}}({\mathrm{\α}}_i^{*}-{\mathrm{\α}}_i)(x_i\·x_k)-\mathrm{\ϵ}---\left(8\right)$

构造回归方程：

$y=\underset{i=1}{\overset{l}{\Σ}}({\mathrm{\α}}_i^{*}-{\mathrm{\α}}_i)(x_i\·x_k)+b---\left(9\right)$

其中，y为训练样本输出。(x_i·x_k)是核函数，这里取高斯径向基函数K(x,x')＝exp(-||x-x'||²/δ²)，是对应的拉格朗日乘子的各个分量。

(H)计算适应度：

设定优化的目标函数，并计算适应度，对每组惩罚因子C和核函数参数δ的效果进行评价，通过相应的误差函数计算适应度。对于误差大的粒子，认为其适应度小，粒子的适应度f_p表示为：

f_p＝1/(e_p+1)(10)

其中，e_p是回归方程的误差函数，即：

$e_p=\frac{1}{l}\underset{i=1}{\overset{l}{\Σ}}{({\hat{y}}_i-y_i)}^2---\left(11\right)$

其中，是回归方程的预测输出，y_i是回归方程的目标输出；

更新每个粒子的速度和位置

$\begin{array}{c}{v}_{id}^k={\mathrm{wv}}_{id}^{k-1}+c_1{\mathrm{rand}}_1(p_{id}-x_{id}^{k-1})+c_2{\mathrm{rand}}_2(p_{gd}-x_{id}^{k-1})\\ x_{id}^k=x_{id}^{k-1}+v_{id}^{k-1}\end{array}---\left(12\right)$

其中，表示待更新粒子的速度，表示待更新粒子的位置，p_id表示待更新粒子的个体最优值，p_gd表示整个粒子群的全局最优值。w为惯性权重，c₁、c₂为加速度常数，rand₁和rand₂为[0,1]范围内的随机数；

对粒子p，如果更新后的适应度值大于原个体的最优值，则更新粒子的个体最优值p_id，即：

p_id＝f_p(13)

如果粒子p的个体最优值p_id大于原来粒子群的全局最优值p_gd，则更新原来的全局最优值p_gd：

p_gd＝p_id(14)

作为一种优选方案：本发明还包括以下步骤：

(I)建模效果评价：

判定建模效果是否满足要求。若满足，结束寻优，得到惩罚因子C和核函数参数δ的最优值，进入步骤(J)，保存ε-SVR惩罚因子C和核函数参数δ的最优参数；否则返回步骤(F)，继续迭代寻优，直到达到最大寻优代数，结束寻优；

(J)保存ε-SVR惩罚因子C和核函数参数δ的最优参数；

(K)读取ε-SVR惩罚因子C和核函数参数δ的最优参数；

(L)读取软测量模型的输入数据；

(M)将输入变量数据中心标准化处理后，使用ε-SVR惩罚因子C和核函数参数δ的最优参数进行高浓磨浆系统游离度在线软测量；

(N)游离度软测量结果显示：在软测量系统人机界面上显示游离度软测量的结果；

(O)游离度软测量结果保存：将本次软测量的输入输出数据保存到相应的数据库，供查询、修正使用；

(P)判断软测量是否结束，结束转至步骤(Q)，否则转至步骤(L)继续进行游离度软测量；

(Q)结束：完成游离度软测量。

作为另一种优选方案：本发明还包括以下步骤：

定期将游离度离线间歇式数据输入到训练集中，更新游离度ε-支持向量回归机软测量模型。

本发明的一种高浓磨浆系统的游离度软测量方法具有以下技术效果：

本发明针对目前高浓磨浆过程中无法直接测量游离度的现状，在没有设备测量游离度的前提下，选取对高浓磨浆系统游离度有直接影响的9个操作变量和检测变量作为软测量模型的输入数据，实现对高浓磨浆系统游离度的软测量，其测量精度高，为实现整个造纸过程的节能降耗提供理论依据和模型支持。

附图说明：

图1是本发明选择的研究对象：高浓磨浆系统三盘磨浆机的结构及13个主要操作变量和检测变量

图2是本发明的一种基于粒子群优化ε-支持向量回归机(PSO-ε-SVR)的高浓磨浆系统游离度软测量算法的流程图

图3是本发明选择的高浓磨浆系统制浆生产过程

高浓磨浆过程中，经预处理后的木片经过传动侧分料螺旋和非传动侧分料螺旋，传动侧负载感应螺旋和非传动侧负载感应螺旋进入高浓磨浆机，经过高浓磨机磨浆后的浆料进入旋风分离器从而实现气浆分离，便形成了粗浆。

图4是本发明寻得最优参数对游离度进行软测量的效果图

具体实施方式

如图1所示，高浓磨浆系统三盘磨浆机由中间一个具有双齿面结构的定盘与两侧两个定盘构成双磨室结构。三盘磨浆机在高转速时运行平稳，不会出现动盘偏斜等问题，在不提高转速及增大磨盘直径的情况下，磨盘面积增加一倍，具有提高产量、改进磨浆质量、方便热能回收等优点。

本发明提供一种高浓磨浆系统三盘磨浆机游离度的软测量方法。该方法能够实现对高浓磨浆机出浆游离度的实时测量，提高人工测量的效率和精度，为工业实时控制和实现节能降耗提供技术支持。

参照图1、图2，一种基于粒子群算法优化ε-支持向量回归机(PSO-ε-SVR)的高浓磨浆系统游离度软测量的方法，其特征在于有常规测量系统、数据采集器、PSO-ε-SVR软测量软件以及运行软测量软件的计算机系统构成。流量计、压力机、震动计、加速度计等常规测量仪表安装于高浓磨浆系统的各个相应位置。数据采集器连接常规测量系统，并通过通信总线连接运行软测量软件的计算机系统。常规测量系统主要包括如下常规测量仪表：

二个流量计(春辉LWJY-LWGY)，分别用于在线测量高浓磨浆系统传动侧稀释水量，非传动侧稀释水量；

一个振动加速度传感器(YMK151A100)，用于在线测量高浓磨浆系统非传动侧的震动；

一个加速度计(PCBSA51SC)，用于在线测量高浓磨浆系统非传动侧震动加速度的大小；

一个电涡流位移传感器(航振HZ891XL)，用于在线测量高浓磨浆系统非传动侧定盘与动盘之间的间隙；

二个压力计(GP50211)，用于在线测量高浓磨浆系统传动侧磨盘压力和非传动侧磨盘压力。

此外，常规测量系统还包括如下两个测量仪：

一个高浓磨浆系统设定产量测量仪：通过传动侧分料螺旋和非传动侧分料螺旋转速设定来设定高浓磨浆系统产量；

一个高浓磨浆系统磨浆机功率测量仪：通过各子系统能耗计量得到高浓磨浆系统磨浆机功率。

本发明实现方法包括：(1)辅助变量选择和模型输入变量确定；(2)PSO-ε-SVR软测量模型的训练和软测量实现。

(1)辅助变量选择和模型输入变量确定

需要软测量的高浓磨浆系统纸浆关键参数为游离度(ml)。根据生产工艺及机理以及可测可观和变量间相关性分析(因子分析、协方差分析)，选择i个(i＝9)软测量辅助变量包括：

设定产量u₁(adt/d)；

磨浆机功率u₂(Mw)；

传动侧稀释水量u₃(l/min)；

非传动侧稀释水量u₄(l/min)；

非传动侧震动加速度u₅(％)；

非传动侧震动u₆(mm/s)；

非传动侧磨盘间隙u₇(mm)；

传动侧磨盘压力u₈(bar)；

非传动侧磨盘压力u₉(bar)。

确定上述9个软测量辅助变量为为软测量模型的输入变量，游离度(ml)作为软测量模型的输出变量。

(2)PSO-ε-SVR软测量模型的训练和软测量实现

(A)开始：变量初始化；

(B)选择进行软测量模型训练，转至(C)，读取软测量训练样本；选择高浓磨浆系统游离度软测量，转至(K)，读取软测量模型最优参数(惩罚因子C和核函数参数δ)；

(C)读取软测量训练样本：从数据库导入或输入训练集其中，软测量模型输入：x_i∈Rⁿ，i＝1,2,…,l，软测量模型输出：l为训练集样本数量；

(D)数据预处理：

采用尖峰滤波算法将高浓磨浆系统的尖峰异常数据进行剔除；由于选取的9个辅助变量存在着量纲差，对每个变量分别进行中心标准化处理，即每类样本数据减去对应样本的平均值，然后除以其样本方差，即：

计算均值：

${\stackrel{\‾}{x}}_i=\frac{1}{l}\underset{j=1}{\overset{l}{\Σ}}{x}_{ij}---\left(15\right)$

计算方差：

${\mathrm{\σ}}_{x_i}^2=\frac{1}{l-1}\underset{j=1}{\overset{l}{\Σ}}(x_{ij}-{\stackrel{\‾}{x}}_i)---\left(16\right)$

标准化：

${\tilde{x}}_{ij}=\frac{x_{ij}-{\stackrel{\‾}{x}}_i}{{\mathrm{\σ}}_{x_i}}---\left(17\right)$

其中，x_ij为第i类变量的第j个训练样本，l为训练样本数，为第i类变量的均值，为中心标准化后的训练样本，表示第i类变量的方差，表示第i类变量的标准差；

(E)确定ε-SVR软测量模型的初始化参数，包括：惩罚因子C和误差容限值ε，核函数的类型及其核函数参数δ。惩罚因子C取20，误差容限值ε取0.01，核函数取高斯径向基函数K(x,x')＝exp(-||x-x'||²/δ²)，其中核函数参数δ＝1；

(F)确定ε-SVR待优化参数惩罚因子C的搜索范围为[0,100]，核函数参数δ的搜索范围为[0,50]，粒子群个体的数目40，粒子的最大迭代寻优次数600，粒子的惩罚因子C迭代速度的最大值为25，粒子的核函数参数δ迭代速度的最大值为25，粒子的初始位置、初始速度、个体极值以及全局极值随机生成；

(G)ε-SVR模型训练：

ε-SVR计算高浓磨浆系统游离度的回归超平面方程y＝(ω·x)+b，其中ω为其法方向，b为其截距。

设模型样本的第i个目标输出为y_i，则ε-SVR将拟合问题等价转化为如下所述的二次规划问题：

$\begin{array}{cc}\underset{w,b}{\min }& \frac{1}{2}\left|\right|\mathrm{\ω}|{|}^2+C\underset{i=1}{\overset{l}{\mathrm{\Σ}}}({\mathrm{\ξ}}_i+{\mathrm{\ξ}}_i^{*})\\ s.t.& \left(\right(\mathrm{\ω}\·x_i)+b)-y_i\≤\mathrm{\ϵ}+{\mathrm{\ξ}}_i,i=1,...,l;\\ & y_i-\left(\right(\mathrm{\ω}\·x_i)+b)\≤\mathrm{\ϵ}+{\mathrm{\ξ}}_i^{*},i=1,...,l;\\ & {\mathrm{\ξ}}_j^{*}\≥0,i=1,...,l.\end{array}---\left(18\right)$

其中，C为惩罚因子，ε为误差容限值，ξ_i和(简记为)为松弛变量。

求解上述二次规划问题的解后，得到回归方程：

y＝(ω·x)+b(19)

将Lagrange函数引入上述问题，导出其对偶问题：

$\begin{array}{cc}\underset{{\mathrm{\α}}^{(*)}\∈R^{2l}}{\min }& \frac{1}{2}\underset{i,j=1}{\overset{l}{\mathrm{\Σ}}}({\mathrm{\α}}_i^{*}-{\mathrm{\α}}_i)({\mathrm{\α}}_j^{*}-{\mathrm{\α}}_j)(x_i\·x_j)\\ & +\mathrm{\ϵ}\underset{i=1}{\overset{l}{\mathrm{\Σ}}}({\mathrm{\α}}_i^{*}+{\mathrm{\α}}_i)-\underset{i=1}{\overset{l}{\mathrm{\Σ}}}{y}_i({\mathrm{\α}}_i^{*}-{\mathrm{\α}}_i)\\ s.t.& \underset{i=1}{\overset{l}{\mathrm{\Σ}}}({\mathrm{\α}}_i^{*}-{\mathrm{\α}}_i)=0,\\ & 0\≤{\mathrm{\α}}_i^{(*)}\≤C,i=1,...,l.\end{array}---\left(20\right)$

解得Lagrange乘子向量

计算b：选择位于开区间(0,C)中的α^(*)的分量α_j或若选到的是α_j，则有：

$b=y_j-\underset{i=1}{\overset{l}{\Σ}}({\mathrm{\α}}_i^{*}-{\mathrm{\α}}_i)(x_i\·x_j)+\mathrm{\ϵ}---\left(21\right)$

若选择的是则有：

$b=y_k-\underset{i=1}{\overset{l}{\Σ}}({\mathrm{\α}}_i^{*}-{\mathrm{\α}}_i)(x_i\·x_k)-\mathrm{\ϵ}---\left(22\right)$

构造回归方程：

$y=\underset{i=1}{\overset{l}{\Σ}}({\mathrm{\α}}_i^{*}-{\mathrm{\α}}_i)(x_i\·x_k)+b---\left(23\right)$

其中，y为训练样本输出。(x_i·x_k)是核函数，这里取高斯径向基函数K(x,x')＝exp(-||x-x'||²/δ²)，是对应的拉格朗日乘子的各个分量。

(H)计算适应度：

设定优化的目标函数，并计算适应度，对每组惩罚因子C和核函数参数δ的效果进行评价，通过相应的误差函数计算适应度。对于误差大的粒子，认为其适应度小，粒子的适应度f_p表示为：

f_p＝1/(e_p+1)(24)

其中，e_p是回归方程的误差函数，即：

$e_p=\frac{1}{l}\underset{j=1}{\overset{l}{\Σ}}{({\hat{y}}_j-y_j)}^2---\left(25\right)$

其中，是回归方程的预测输出，y_i是回归方程的目标输出；

更新每个粒子的速度和位置

$\begin{array}{c}{v}_{id}^k={\mathrm{wv}}_{id}^{k-1}+c_1{\mathrm{rand}}_1(p_{id}-x_{id}^{k-1})+c_2{\mathrm{rand}}_2(p_{gd}-x_{id}^{k-1})\\ x_{id}^k=x_{id}^{k-1}+v_{id}^{k-1}\end{array}---\left(26\right)$

其中，表示待更新粒子的速度，表示待更新粒子的位置，p_id表示待更新粒子的个体最优值，p_gd表示整个粒子群的全局最优值。w为惯性权重，c₁为加速度常数，rand₁和rand₂为[0,1]范围内的随机数；

对粒子p，如果更新后的适应度值大于原个体的最优值，则更新粒子的个体最优值p_id，即：

p_id＝f_p(27)

如果粒子p的个体最优值p_id大于原来粒子群的全局最优值p_gd，则更新原来的全局最优值p_gd：

p_gd＝p_id(28)

(I)建模效果评价：

判定建模效果是否满足要求。若满足，结束寻优，得到惩罚因子C和核函数参数δ的最优值，进入步骤(J)，保存ε-SVR惩罚因子C和核函数参数δ的最优参数；否则返回步骤(F)，继续迭代寻优，直到达到最大寻优代数，结束寻优；

(J)保存ε-SVR惩罚因子C和核函数参数δ的最优参数；

(K)读取ε-SVR惩罚因子C和核函数参数δ的最优参数；

(L)读取软测量模型的输入数据；

(M)将输入变量数据中心标准化处理后，使用ε-SVR惩罚因子C和核函数参数δ的最优参数进行高浓磨浆系统游离度在线软测量；

(N)游离度软测量结果显示：在软测量系统人机界面上显示游离度软测量的结果；

(O)游离度软测量结果保存：将本次软测量的输入输出数据保存到相应的数据库，供查询、修正使用；

(P)判断软测量是否结束，结束转至步骤(Q)，否则转至步骤(L)继续进行游离度软测量；

(Q)结束：完成游离度软测量。

作为优选的一种方案：所述软测量方法还包括以下步骤：

定期将游离度离线间歇式数据输入到训练集中，更新游离度ε-支持向量回归机软测量模型。

图4为一段时间内，高浓磨浆系统游离度软测量的效果，可以看出游离度的预测值与其实际值基本一致，误差比较小，且变化趋势基本一致。此外，本发明方法速度快、精度高、泛化能力强且有严格的数学解释，相比于其他方法有较高的优越性。因此本发明是一种具有很高实用价值的、低成本的高浓磨浆系统游离度计量手段。

可以理解的是，以上关于本发明的具体描述，仅用于说明本发明而并非受限于本发明实施例所描述的技术方案，本领域的普通技术人员应当理解，仍然可以对本发明进行修改或等同替换，以达到相同的技术效果；只要满足使用需要，都在本发明的保护范围之内。

Claims (2)

1.一种高浓磨浆系统游离度的软测量方法，包括以下步骤：

(1)辅助变量选择和模型输入变量确定

需要软测量的高浓磨浆系统纸浆关键参数为游离度，ml：根据生产工艺及机理和变量间相关性分析，所述相关性分析为因子分析、协方差分析，选择9个软测量辅助变量，包括：设定产量u₁，adt/d；磨浆机功率u₂，Mw；传动侧稀释水量u₃，l/min；非传动侧稀释水量u₄，l/min；非传动侧震动加速度u₅，％；非传动侧震动u₆，mm/s；非传动侧磨盘间隙u₇，mm；传动侧磨盘压力u₈，bar；非传动侧磨盘压力u₉，bar；确定上述9个软测量辅助变量为软测量模型的输入变量，软测量模型的输出变量为游离度，ml；

(2)PSO-ε-SVR软测量模型的训练和软测量实现

(A)开始：变量初始化；

(B)选择进行软测量模型训练，转至(C)，读取软测量训练样本；选择高浓磨浆系统游离度软测量，转至(K)；

(C)读取软测量训练样本：从数据库导入或输入训练集其中，软测量模型输入：x_i∈Rⁿ，i＝1,2,…,l，软测量模型输出：i＝1,2,…,l，l为训练集样本数量；

(D)数据预处理：

采用尖峰滤波算法将高浓磨浆系统的尖峰异常数据进行剔除；9个辅助变量存在着量纲差，对每个变量分别进行中心标准化处理，即每类样本数据减去对应样本的平均值，然后除以其样本方差，即：

计算均值：

${\stackrel{\‾}{x}}_i=\frac{1}{l}\underset{j=1}{\overset{l}{\Σ}}{x}_{ij}---\left(1\right)$

计算方差：

${\mathrm{\σ}}_{x_i}^2=\frac{1}{l-1}\underset{j=1}{\overset{l}{\Σ}}(x_{ij}-{\stackrel{\‾}{x}}_i)---\left(2\right)$

标准化：

${\tilde{x}}_{ij}=\frac{x_{ij}-{\stackrel{\‾}{x}}_i}{{\mathrm{\σ}}_{x_i}}---\left(3\right)$

其中，x_ij为第i类变量的第j个训练样本，l为训练样本数，为第i类变量的均值，为中心标准化后的训练样本，表示第i类变量的方差，表示第i类变量的标准差；

(E)确定ε-SVR软测量模型的初始化参数，包括：惩罚因子C和误差容限值ε，核函数的类型及其核函数参数δ；惩罚因子C取20，误差容限值ε取0.01，核函数取高斯径向基函数K(x,x')＝exp(-||x-x'||²/δ²)，其中核函数参数δ＝1；

(F)确定ε-SVR待优化参数惩罚因子C的搜索范围为[0,100]，核函数参数δ的搜索范围为[0,50]，粒子群个体的数目为40，粒子的最大迭代寻优次数为600，粒子的惩罚因子C迭代速度的最大值为25，粒子的核函数参数δ迭代速度的最大值为25，粒子的初始位置、初始速度、个体极值以及全局极值随机生成；

(G)ε-SVR模型训练：

ε-SVR计算高浓磨浆系统游离度的回归超平面方程y＝(ω·x)+b，其中ω为其法方向，b为其截距；

设模型样本的第j个目标输出为y_j，则ε-SVR将拟合问题等价转化为如下的二次规划问题：

$\begin{array}{cc}\underset{\mathrm{\ω},b}{min}& \frac{1}{2}\left|\right|\mathrm{\ω}|{|}^2+C\underset{i=1}{\overset{l}{\Σ}}({\mathrm{\ξ}}_i+{\mathrm{\ξ}}_i^{*})\end{array}---\left(4\right)$

s.t.((ω·x_i)+b)-y_i≤ε+ξ_i,i＝1,…,l；

$y_j-((\mathrm{\ω}\·x_i)+b)\≤\mathrm{\ϵ}+{\mathrm{\ξ}}_i^{*},i=1,...,l;$

${\mathrm{\ξ}}_i^{*}\≥0,i=1,...,l.$

其中，C为惩罚因子，ε为误差容限值，ξ_i和为松弛变量，简记为

求解式(4)中二次规划问题的解(ω,b,ξ)后，得到回归超平面方程：

y＝(ω·x)+b(5)

将Lagrange函数引入上述问题，导出其对偶问题：

$\begin{array}{c}\underset{{\mathrm{\α}}^{(*)}\∈R^{2l}}{\min }\frac{1}{2}\underset{i,j=1}{\overset{l}{\Σ}}({\mathrm{\α}}_i^{*}-{\mathrm{\α}}_i)({\mathrm{\α}}_j^{*}-{\mathrm{\α}}_j)(x_i\·x_j)\\ +\mathrm{\ϵ}\underset{i=1}{\overset{l}{\Σ}}({\mathrm{\α}}_i^{*}+{\mathrm{\α}}_i)-\underset{i=1}{\overset{l}{\Σ}}{y}_i({\mathrm{\α}}_i^{*}-{\mathrm{\α}}_i)\end{array}---\left(6\right)$

$\begin{array}{cc}s.t.& \underset{i=1}{\overset{l}{\Σ}}({\mathrm{\α}}_i^{*}-{\mathrm{\α}}_i)=0,\end{array}$

$0\≤{\mathrm{\α}}_i^{(*)}\≤C,i=1,...,l.$

解得Lagrange乘子向量 ${\mathrm{\α}}^{(*)}={({\mathrm{\α}}_1,{\mathrm{\α}}_1^{*},...,{\mathrm{\α}}_l,{\mathrm{\α}}_l^{*})}^T$

计算b：选择α^(*)的分量α_j或若选到的是α_j，则有：

$b=y_j-\underset{i=1}{\overset{l}{\Σ}}({\mathrm{\α}}_i^{*}-{\mathrm{\α}}_i)(x_i\·x_j)+\mathrm{\ϵ}---\left(7\right)$

若选择的是则有：

$b=y_k-\underset{i=1}{\overset{l}{\Σ}}({\mathrm{\α}}_i^{*}-{\mathrm{\α}}_i)(x_i\·x_k)-\mathrm{\ϵ}---\left(8\right)$

构造回归方程：

$y=\underset{i=1}{\overset{l}{\Σ}}({\mathrm{\α}}_i^{*}-{\mathrm{\α}}_i)(x_i\·x_k)+b---\left(9\right)$

其中，y为训练样本i的输出；(x_i·x_k)是核函数，这里取高斯径向基函数K(x,x')＝exp(-||x-x'||²/δ²)，是对应的拉格朗日乘子的各个分量；

(H)计算适应度：

设定优化的目标函数，并计算适应度，对每组惩罚因子C和核函数参数δ的效果进行评价，通过相应的误差函数计算适应度；对于误差大的粒子，认为其适应度小，粒子的适应度f_p表示为：

f_p＝1/(e_p+1)(10)

其中，e_p是回归方程的误差函数，即：

$e_p=\frac{1}{l}\underset{i=1}{\overset{l}{\Σ}}{({\hat{y}}_i-y_i)}^2---\left(11\right)$

其中，是回归方程的预测输出，y_i是回归方程的目标输出；

更新每个粒子的速度和位置

$\begin{array}{c}{v}_{id}^k={\mathrm{wv}}_{id}^{k-1}+c_1{\mathrm{rand}}_1(p_{id}-x_{id}^{k-1})+c_2{\mathrm{rand}}_2(p_{gd}-x_{id}^{k-1})\\ x_{id}^k=x_{id}^{k-1}+v_{id}^{k-1}\end{array}---\left(12\right)$

其中，表示待更新粒子的速度，表示待更新粒子的位置，p_id表示待更新粒子的个体最优值，p_gd表示整个粒子群的全局最优值；w为惯性权重，c₁、c₂为加速度常数，rand₁和rand₂为[0,1]范围内的随机数；

对粒子p，如果更新后的适应度值大于原个体的最优值，则更新粒子的个体最优值p_id，即：

p_id＝f_p(13)

如果粒子p的个体最优值p_id大于原来粒子群的全局最优值p_gd，则更新原来的全局最优值p_gd：

p_gd＝p_id(14)

还包括以下步骤：

(I)建模效果评价：

判定建模效果是否满足要求；若满足，结束寻优，得到惩罚因子C和核函数参数δ的最优值，进入步骤(J)，保存ε-SVR惩罚因子C和核函数参数δ的最优参数；否则返回步骤(F)，继续迭代寻优，直到达到最大寻优代数，结束寻优；

(J)保存ε-SVR惩罚因子C和核函数参数δ的最优参数；

(K)读取ε-SVR惩罚因子C和核函数参数δ的最优参数；

(L)读取软测量模型的输入数据；

(M)将输入变量数据中心标准化处理后，使用ε-SVR惩罚因子C和核函数参数δ的最优参数进行高浓磨浆系统游离度在线软测量；

(N)游离度软测量结果显示：在软测量系统人机界面上显示游离度软测量的结果；

(O)游离度软测量结果保存：将本次软测量的输入输出数据保存到相应的数据库，供查询、修正使用；

(P)判断软测量是否结束，结束转至步骤(Q)，否则转至步骤(L)继续进行游离度软测量；

(Q)结束：完成游离度软测量。

2.根据权利要求1所述一种高浓磨浆系统游离度的软测量方法，其特征在于还包括以下步骤：

定期将游离度离线间歇式数据输入到训练集中，更新游离度ε-支持向量回归机软测量模型。