辛烷值损失预测模型的构建与优化

Construction and Optimization of Octane Loss Prediction Model

摘要: 汽油是小型车辆的主要燃料，汽油燃烧产生的尾气排放对大气环境有重要影响。针对我国汽油在精制过程中辛烷值损失的问题，首先使用Lasso回归对工业样本数据进行预处理，再运用多元线性回归建立辛烷值损失预测模型，最后基于数据包络分析法(DEA)对汽油精制过程中的主要操作变量进行优化。该优化模型具有实际应用意义，为汽油清洁化过程降低辛烷值损失的研究提供了理论依据。

Abstract: Gasoline is the major fuel for small-sized vehicles, and the gasoline-fueled vehicle emissions generated by gasoline combustion exert an important impact on the atmosphere. Regarding the loss of octane value amid the process of gasoline refining in our country, t Lasso regression is used to process the industrial specimen data beforehand, and then a prediction model for octane loss is established through multiple linear regression. Finally, the major operating variables used in gasoline refining are optimized based on Data Envelopment Analysis (DEA) model. The optimized model is of practical significance for application and proffers theoretical ground for the reduction of the loss of octane value during the refining of gasoline.

1. 引言

目前，小型车辆的主要燃料是汽油，汽油燃烧排放的产物一氧化碳、氮氧化合物、碳氢化合物和微粒等对大气环境造成了重大影响。为此，世界主要大国均制定了日益严格的汽油质量标准。汽油清洁化的重点是降低汽油中的硫和烯烃含量，同时尽量保持其辛烷值。辛烷值(RON)是衡量汽油燃烧性能的最重要指标，目前除考虑改进发动机的构造和操作条件外，主要依靠提高汽油辛烷值，借以提高发动机压缩比的办法来增加燃烧功率和节约汽油。然而，现有技术在对催化裂化汽油进行脱硫和降烯烃过程中，普遍降低了汽油辛烷值 [1]。辛烷值每降低1个单位，造成经济损失约150元/吨。以一个年产100万吨的催化裂化汽油精制装置为例，若能降低0.3个单位RON损失，其经济效益将达到4500万元。

为满足日益严格的清洁汽油标准而不断降低硫和烯烃含量的需求，科研人员们在汽油清洁化领域开展了大量的研究工作。如冯连坤指出汽油辛烷值损失的影响因素为预加氢反应器温度、加氢脱硫反应器入口温度和轻汽油抽出量 [2]；刘畅采用随机森林算法去寻找降低汽油精制过程中的辛烷值损失模型中的主要变量 [3]，而本文使用Lasso回归方法寻找主要变量。对于有不同取值的属性的数据，取值划分较多的属性会对随机森林产生更大的影响，所以随机森林在这种数据上产出的属性权值是不可信的。Lasso算法是一种带有惩罚因子的线性模型估计方法，该方法的本质是约束各个回归系数的绝对值之和小于某个常数的条件下，最小化回归方程的残差平方和，同时阈值的设定又可以收缩每个估计的参数值 [4] [5] [6]。Lasso方法可以有效地估计回归模型中的各个参数，同时也可以较好地解决变量间的多重线性问题。

本文首先通过Lasso回归找出精制过程中的主要操作变量，再采用多元线性回归方法建立辛烷值损失预测模型。多元线性回归不需要很复杂的计算，在数据量大的情况下依然运行速度很快，并且可以根据系数给出每个变量的理解。本文最后一部分运用了数据包络分析方法(DEA)对主要操作变量进行调整来达到降低辛烷值损失的目的。DEA是一个对于多投入多产出的多个决策单元的效率评价方法，在汽油精制过程中，将多个主要操作变量作为输入，辛烷值的损失作为输出，利用输出型CCR模型来计算并进行调整。

2. 变量降维

2.1. 数据预处理

本文具体数据来源于某石化企业的催化裂化汽油精制脱硫装置运行4年所积累的历史数据。由于工厂得到的原始数据存在一定数据缺失和数据失真的情况，需要排除对数据中的坏值或者短缺值，对失真的数据进行修正。

数据预处理的具体步骤如下：

Step 1：遍历每个样本，查看样本是否大多变量为空值，若是则剔除该样本；查看变量是否全列为空值，若是则删除该列，若部分是，用平均值补充空值。

Step 2：根据各变量取值范围来剔除不符样本。

Step 3：根据 $3\sigma$ 准则删除不符样本。

Step 4：对于Matlab，xlswrite (filename, data)为保存到excel的命令。

$3\sigma$ 准则：首先设对被测量变量进行等精度测量，得到 $x_1,x_2,\cdots ,x_n$，算出其算术平均值x及剩余误差 $v_i=x_i-x\left(i=1,2,\cdots ,n\right)$，并按贝塞尔公式算出标准误差 $\sigma$。若某个测量值 $x_b$ 的剩余误差 $v_b\left(1\le b\le n\right)$ 满足 $\left|v_b\right|=\left|x_b-x\right|>3\sigma$，则认为 $x_b$ 是含有粗大误差值的坏值，应予剔除。

2.2. 数据降维

在解决实际问题时，自变量过多对模型的构建弊大于利，因此，一般通过数据降维来简化模型。本文中原有367个操作变量，因为变量数量过多且大于样本数据数量，主成分分析法等常规方法不能很好地实行降维，因此本文利用构造惩罚函数的Lasso回归对多个变量进行筛选，选出具有代表性且独立性的19个建模主要变量。降维的好处有：1) 减少数据存储所需的空间，节约成本。2) 减少数据处理与建模的时间，提高效率。3) 提高该算法的性能，因为会有一些算法在此300维的数据上变现不佳。4) 能更直观地看出降维的结果，有助于数据可视化。

2.2.1. Lasso回归

Lasso (Least absolute shrinkage and selection operator)方法是以缩小变量集(降阶)为思想的压缩估计方法，在拟合广义线性模型的同时进行变量筛选和复杂调整。通过构造一个惩罚函数，可以将变量的系数进行压缩并使某些回归系数变为0，进而达到变量选择的目的。普通线性模型

$Y=X\beta +\epsilon$ (1)

其中，响应变量 $Y={\left(y_1,y_2,\cdots ,y_n\right)}^{\text{T}}$，协变量 $X=\left(X^{\left(1\right)},X^{\left(2\right)},\cdots ,X^{\left(d\right)}\right)$，对于每一个 $X^{\left(j\right)}$，有 $X^{\left(j\right)}={\left(x_1^j,x_2^j,\cdots ,x_n^j\right)}^{\text{T}}$。假设每个 $x_i^j\left(i=1,2,\cdots ,n;j=1,2,\cdots ,d\right)$ 均中心化和标准化，则随机误差项 ${\epsilon }_i~N\left(0,{\sigma }^2\right)$， $\left(i=1,2,\cdots ,n\right)$， $\epsilon ={\left({\epsilon }_1,{\epsilon }_2,\cdots ,{\epsilon }_n\right)}^{\text{T}}$，回归系数 $\beta ={\left({\beta }_1,{\beta }_2,\cdots ,{\beta }_d\right)}^{\text{T}}$。

当X为列满秩设计矩阵时，回归系数 $\beta$ 可由普通最小二乘估计方法求得

${\hat{\beta }}_{OLS}=\arg \underset{\beta \in R^d}{\min }{\Vert Y-X\beta \Vert }^2={\left(X^{\text{T}}X\right)}^{-1}{X}^{\text{T}}Y$

惩罚方法：当设计矩阵X不满足列满秩时，将不能采用普通最小二乘法来求解回归系数，此时应引入惩罚方法。该方法可以同时实现变量选择和参数估计，在参数估计时，通过将部分参数压缩为0来达到变量选择的目的。惩罚方法是取惩罚似然函数最小时的值作为回归系数的估计值，即

$\hat{\beta }=\arg \underset{\beta \in R^d}{\min }\left({\Vert Y-X\beta \Vert }^2+P_{\lambda }\left(\left|\beta \right|\right)\right)$，

其中惩罚项

$P_{\lambda }\left(\left|\beta \right|\right)=\lambda {\displaystyle \underset{j=1}{\overset{d}{\sum }}{\left|{\beta }_j\right|}^m},m\ge 0$，

$\lambda$ 为调节参数(也可以为向量)。当 $m=2$ 时，得 $L_2$ 惩罚项(Ridge惩罚)。

Lasso方法则是在普通线性模型中增加 $L_1$ 惩罚项。具体算法如下：

1) 令 ${\sigma }_i=sign\left({\hat{\beta }}^0\right)$，这里 ${\hat{\beta }}^0$ 为通常的最小二乘估计。

2) 计算 $\hat{\beta }=\arg \underset{\beta }{\min }\left\{{\displaystyle \underset{i=1}{\overset{N}{\sum }}{\left(y_i-{\displaystyle \underset{j}{\sum }{\beta }_j{x}_{ij}}\right)}^2}\right\}$， $\text{s}\text{.t}\text{.}G\beta \le t$，这里 $G={\sigma }_i^{\text{T}}$。

3) 验证是否满足 ${\displaystyle \underset{j}{\sum }\left|{\hat{\beta }}_j\right|}\le t$，若满足则停止， $\hat{\beta }$ 即为所求，否则继续。

4) 令 ${\sigma }_j=sign\left(\hat{\beta }\right)$ 加入到G中，即令 $G=\left(\begin{array}{l}{\sigma }_i^{\text{T}}\\ {\sigma }_j^{\text{T}}\end{array}\right)$，然后返回(2)。

LASSO回归复杂度调整的程度由参数λ来控制，λ越大对变量较多的线性模型的惩罚力度就越大，从而最终获得一个变量较少的模型。LASSO回归与Ridge回归同属于一个被称为Elastic Net的广义线性模型家族。这一家族的模型除了相同作用的参数λ之外，还有另一个参数α来控制应对高相关性数据时模型的性状。LASSO回归α = 1，Ridge回归α = 0，一般Elastic Net模型0 < α < 1。

2.2.2. 主要变量筛选结果

首先，从与辛烷损失值相关的367个变量中筛选出主要变量。设辛烷损失值为因变量，原料经过一系列工艺流程得到最终产品，此处研究的是产品的性质，不宜作为自变量，所以选取作为因变量使用。辛烷值损失可近似于原料辛烷值与产品辛烷值的差值，所以若把原料性质和产品性质一起进行筛选，那原料辛烷值和产品辛烷值对辛烷值损失值的影响是巨大的，会让操作变量所造成的影响几乎不显现，这与我们所求的实际目标不符。因此我们去除2个产品的性质变量，只对365个变量使用Lasso回归，根据其λ的变化利用Matlab编程实现。由于变量太多，从Matlab运行结果的参数系数图中不能看出具体的系数情况。使用交叉验证选出最佳的λ，其筛选的主要变量如下，见表1。

Table 1. Main variables

表1. 主要变量

3. 辛烷值(RON)损失预测

在解决实际问题时，实验结果往往受多个因素的影响，这时可以通过多元回归分析求出试验指标(因变量) y与多个试验因素(自变量) $x_j\left(j=1,2,\cdots ,m-1\right)$ 之间的近似函数关系 $y=f\left(x_1,x_2,\cdots ,x_{m-1}\right)$，亦称多重回归或复回归。当影响因素与因变量之间呈线性关系时，所进行的回归分析即为多元线性回归。多元线性回归是多元统计分析中的一种重要方法，被广泛应用于社会、经济、技术以及众多自然科学领域的研究中 [7] [8] [9] [10] [11]。其最主要的思想是多个影响因素与因变量之间是线性关系，设y是可观测的随机变量，它受 $m-\text{1}$ 个非随机因素 $x_1,x_2,\cdots ,x_k$ 和随机因素 $\epsilon$ 的影响。若y与 $x_1,x_2,\cdots ,x_{m-1}$ 有如下线性关系

$y={\beta }_0+{\beta }_1{x}_1+{\beta }_2{x}_2+\cdots +{\beta }_{m-1}{x}_{m-1}+\epsilon$

其中 ${\beta }_0,{\beta }_1,{\beta }_2,\cdots ,{\beta }_{m-1}$ 是未知参数； $\epsilon$ 是均值为0、方差为 ${\sigma }^2>0$ 的不可观测的随机变量，称为误差项，并通常假定 $\epsilon ~N\left(0,{\sigma }^2\right)$，称

$\hat{y}=b_0+b_1{x}_{i1}+b_2{x}_{i2}+\cdots +b_{m-1}{x}_{i,m-1},i=1,2,\cdots ,n$

为多元线性回归模型，且称y为因变量， $x_1,x_2,\cdots ,x_{m-1}$ 为自变量。

本文明确定义研究的自变量与因变量，将辛烷值的损失值设为因变量Y，自变量为问题二中筛选出的19个变量，用SPSS软件计算这些变量间的关联度，计算结果中变量的VIF均小于5。

VIF是方差膨胀因子，是容差的倒数，说明这些变量之间不存在多重共线性，关联度很低。因此，先选择多元线性回归分析方法进行建模。通过MATLAB计算，得：

从中可以看出：回归方程是

$\begin{array}{c}y=-8.898656+0.025845x_1+0.017971x_2-0.000324x_3\\ -0.000008x_4-0.000280x_5-0.000447x_6-0.001817x_7\\ +0.014382x_8-0.000011x_9-0.002475x_{10}-0.000566x_{11}\\ +0.000006x_{12}-0.088715x_{13}-0.010613x_{14}-0.022320x_{15}\\ +0.001933x_{16}-0.249645x_{17}+0.025012x_{18}+0.015629x_{19}\end{array}$

由方差分析表的显著性可以看出，变量 $x_1,x_2,\cdots ,x_{19}$ 对y的影响都是显著的，说明选取的变量具有一定代表性，此方程的剩余标准差S接近0，但可决系数R = 0.5682，说明回归效果一般，有较大改进的余地。

从325个样本中选取最后50个样本作为上述模型的误差分析，可得出该模型的预测值与RON损失值的真实数据平均误差为0.0981。

4. 方案优化

针对主要操作变量的优化问题，要求在保证产品硫含量不大于5 μg/g的前提下，使得样本的辛烷值(RON)损失降幅大于30%。本文采用数据包络分析方法(DEA)进行方案优化，DEA是一种针对多投入多产出问题的效率评价以及优化模型。

4.1. DEA模型

4.1.1. DEA模型概述

DEA模型，又称为数据包络分析方法，是1978年由美国运筹学家A. Charnes和W.W. Cooper提出的一种定量分析方法，常用的有CCR和BCC等。DEA以相对效率概念为基础，以凸分析和线形规划为工具，应用数学规划模型计算比较决策单元之间的相对效率，对评价对象做出评价 [12] [13] [14]。DEA能充分考虑对于决策单元本身最优的投入产出方案，因而能够更理想地反映评价对象自身的信息和特点，同时对于评价复杂系统的多投入多产出分析具有独到之处。

4.1.2. DEA原理

设有n个决策单元，每个决策单元都有m个输入和s个输出，分别记第j个决策单元的输入和输出为

$X_j={\left(x_{1j},x_{2j},\cdots ,x_{mj}\right)}^{\text{T}},j=1,2,\cdots ,n$

$Y_j={\left(y_{1j},y_{2j},\cdots ,y_{sj}\right)}^{\text{T}},j=1,2,\cdots ,n$

其中 $x_{ij}>0$， $y_{rj}>0$ ( $i=1,2,\cdots ,m$， $r=1,2,\cdots ,s$ )分别表示第j个决策单元 ${\text{DMU}}_j$ 的第i种类型的输入量和第r种类型的输出量。DEA的产出指向型C²R模型如下所示

$\begin{array}{l}\min E_{j0}=w^{\text{T}}{X}_{j0}\\ \text{s}\text{.t}\text{.}{w}^{\text{T}}{X}_j-{\mu }^{\text{T}}{Y}_j\ge 0,j=1,2,\cdots ,n\\ {\mu }^{\text{T}}{Y}_{j0}=1\\ w\ge 0,\mu \ge 0\end{array}$ (2)

利用对偶原理得到产出指向型的对偶模型

$\begin{array}{l}\max z\\ \text{s}\text{.t}\text{.}{\lambda }^{\text{T}}{X}_j\le X_{j0},j=1,2,\cdots ,n\\ {\lambda }^{\text{T}}{Y}_j\ge z{Y}_{j0},j=1,2,\cdots ,n\\ \lambda \ge 0,z\in E^n\end{array}$ (3)

判断决策单元的有效性，本质上是判断它是否位于生产可能集的生产前沿面上。如果决策单元不为DEA有效，通过计算、求解，可以对原有的投入向量和产出向量进行调整，使其成为DEA有效，我们称经过调整后的点为决策单元在生产前沿面上的“投影”，即为优化过程。求解最优解的模型如下：

设 ${\lambda }_0$， $S^{-0}$， $S^{+0}$， $z^0$ 为上述模型的最优解，令

${\hat{X}}_0=X_0-S^{-0}$

${\hat{Y}}_0=z{Y}_0+S^{+0}$

称 $\left({\hat{X}}_0,{\hat{Y}}_0\right)$ 为 ${\text{DMU}}_{j_0}$ 在生产可能集 $T_{C^{2}R}$ 的生产前沿面上的“投影”。

4.2. 优化结果

根据优化过程中原料、待生吸附剂和再生吸附剂的性质保持不变，此时所剩变量共19个。根据工艺要求，19个变量的取值范围应作为约束条件，再依据优化目标，得到其他约束条件如下：

$\{\begin{array}{l}y\le \hat{y}\ast 0.7\\ S\le 5,\end{array}$

其中 $\hat{y}$ 为原始数据中的辛烷(RON)损失值。

本文将19个主要变量作为输入指标，辛烷值和硫含量作为两个输出指标，总共325个样本，这满足数据包络分析法的使用要求。利用CCR产出模型计算后对各个样本的操作变量进行调整，优化操作方案使其满足约束条件，采用Matlab编程实现，部分计算结果见表2。

Table 2. Results of Matlab running

表2. Matlab运行结果

由于分析的样本数据过多，以上为截取部分数据的计算结果。为达优化目标，通过包络分析法分析出每组样本数据的调整情况，图中数据即为每组样本数据所进行优化操作变量后的值。

优化方案是每个数据样本中19个变量所对应原数据值与此图中对应变量数据值的调整区间。例：用图中编号2样本数据为例，为达到优化目标，x₂变量建议下调至35.00525，x₃变量建议下调至470.9836，其他变量操作相同。

5. 结束语

本文研究了一般石化企业的催化裂化汽油精制脱硫装置的模型构建。首先对原始数据进行预处理，基于处理后的数据运用Lasso回归方法通过Matlab软件对变量进行降维处理。然后，通过Matlab语言编程分析降维后的变量与RON损失的线性关系，建立RON损失预测模型，并通过预测结果与实际值的对比验证了RON损失预测模型的可靠性。最后基于数据包络分析法(DEA)对汽油精制过程中的主要操作变量进行优化。建立的模型对降低辛烷损失值具有重要的理论意义和实际应用价值。