【实战练习】汽油辛烷值优化建模（二）（问题一的详细讲解）

1.问题回顾

问题1为数据处理，需参考近4年的工业数据(见附件一“325个数据样本数据.xlsx”)的预处理结果，依“样本确定方法”（附件二）对285号和313号数据样本进行预处理（原始数据见附件三“285号和313号样本原始数据.xlsx”）并将处理后的数据分别加入到附件一中相应的样本号中，供下面研究使用。

2.问题1的分析

（1）准备工作

首先明确附件二的数据处理方法，将附件三中285号样本数据和313号样本数据汇总整理，并读入Python中形成矩阵，然后在Python中对矩阵进行数据处理的步骤。

附件二中给出的数据处理步骤分为五步：

第一步是对于只含有部分时间点的位点，如果其残缺数据较多，无法补充，将此类位点删除
第二步是删除样本中数据全部为空值的位点。
第三步是对于部分数据为空值的位点，空值处用其前后两个小时数据的平均值代替。
第四步是要求根据工艺要求与操作经验，总结出原始数据变量的操作范围，然后采用最大最小的限幅方法剔除一部分不在此范围的样本。
第五步是根据3σ3\sigma3σ准则去除异常值。

题目给出了285 号和313 号样本所采集的原始数据，包括其各自的7 个原料性质、2个待生吸附剂性质、2 个再生吸附剂性质、2 个产品性质等变量以及另外354 个操作变量的原始数据。

我们的算法操作步骤展示如下：

285 号和313 号样本主要原始数据汇总结果如下：

数据类别	285 号	313 号
采样日期	2017.07.17	2017.05.15
操作变量采样时段	6:00:00-8:00:00	6:00:00-8:00:00
辛烷值采样时间点	8:00:00	8:00:00
操作变量数据采集量	40x354	40x354

（2）异常数据的剔除

1）时间不连续数据：

因 285 号和 313 号的采样日期处于2017.04-2019.09 时间范围内，故其数据采集频次应为 3min/次。遍历两样本原始数据，检测第i组和第i+1组数据采集时间T的差值，若存在：
∣Tt+1−Tt∣>3min|T_{t+1}-T_{t}|>3\quad min ∣Tt+1−Tt∣>3min
则认为该时间点数据缺失。之后判断时间段缺失程度，若其残缺数据较多，无法补充，则将此类位点删除。若只有小段数据缺失，可将此部分数据标为 NaN，之后按附件二的要求进行填充。

我们看一下285号样本的数据：

import numpy as np
import pandas as pd
data=pd.read_excel('附件三：285号和313号样本原始数据.xlsx',sheet_name='操作变量285')
data=pd.read_excel('附件三：285号和313号样本原始数据.xlsx',sheet_name='操作变量285')
data

我们得到的样本数据为：

times=data.iloc[:,0]
times

from datetime import datetime, date
Seconds=[]
for i in range(1,40):time_i=data.iloc[i-1,0]time_i_1=data.iloc[i,0]time_i_struct = datetime.strptime(time_i.strip(), "%Y-%m-%d %H:%M:%S")time_i_1_struct = datetime.strptime(time_i_1.strip(), "%Y-%m-%d %H:%M:%S")seconds = (time_2_struct - time_1_struct).secondsSeconds.append(seconds/60)
Seconds

我们得到的区间为：

第一眼看去几乎全为3！

for i in Seconds:counts=0if i>3.0:counts+=1else:counts=0
counts

经过上面的检验，我们发现285号样本数据没有时间位点超过3分钟的数据。

同理应用于313号样本，同样没有时间位点超过3分钟的数据点。

综上所述，附件三中不存在因时间采集时间不连续导致的缺失数据。

2）全为空值的样本

经过查找，不存在空值。

3）替换空值

不存在空值的数据，故这一步不需要处理。

4）超范围数据

此部分数据处理工作需根据工艺要求与操作经验总结出原始数据变量的操作范围，然后采用最大最小的限幅方法剔除一部分不在此范围的样本。

首先依据附件四“354 个操作变量信息”确定第 j 个数据变量对应操作范围的最小值和最大值。

我们首先看一下附件四中的数据区间：

使用Excel等软件是很难把区间分开的，因此我们需要使用程序的思想。

我们导入附件四的内容：

data_range=pd.read_excel('附件四：354个操作变量信息.xlsx')
data_range.head()

import re
def get_min_range_value(data):try:if data[0]=='-':return -float(data.split('-')[1])else:return float(data.split('-')[0])except:print(data.split('-'))def get_max_range_value(data):if ('(' in data) or (')' in data):try:temp=re.search('\((.*?)\)',data).group(1)except:temp=re.search("（(.*?)）",data).group(1)return float(temp)try:return float(data.split('-')[-1])except:print(data)

data_range['min_region']=data_range.apply(lambda x:get_min_range_value(x['取值范围']),axis=1)
data_range['max_region']=data_range.apply(lambda x:get_max_range_value(x['取值范围']),axis=1)
data_range

我们展示结果为：

如上图所示，我们成功将区间分成两列！

我们检查285号样本数据不在范围内的数据点：

data=pd.read_excel('./源数据/附件三：285号和313号样本原始数据.xlsx',sheet_name='285操作变量')
data

data=data.iloc[:,1:]
data

def check_data(data,min_values,max_values):if (data > max_values) or (data < min_values):return np.nanelse:return data
for i in range(data_range.shape[0]):names=data_range.iloc[i,1]data_min=data_range.iloc[i,6]data_max=data_range.iloc[i,7]data[names]=data[names].apply(lambda x:check_data(x,data_min,data_max))
data

我们已经把不在范围内的点替换成空值了。

data.isnull().sum()

因为数据很长，我们无法准确看到究竟那一列数据空值点，我们做一步查找：

data.isnull().sum()[data.isnull().sum()!=0]

我们发现这三列数据存在严重问题！全部不在区间范围内。

同样的方法，我们检查313样本点。

data_1=pd.read_excel('./源数据/附件三：285号和313号样本原始数据.xlsx',sheet_name='313操作变量')
data_1=data_1.iloc[:,1:]
data_1

def check_data(data_1,min_values,max_values):if (data_1 > max_values) or (data_1 < min_values):return np.nanelse:return data_1
for j in range(data_range.shape[0]):names=data_range.iloc[j,1]data_min=data_range.iloc[j,6]data_max=data_range.iloc[j,7]data_1[names]=data_1[names].apply(lambda x:check_data(x,data_min,data_max))
data_1

data_1.isnull().sum()[data_1.isnull().sum()>0]

结果如下：

S-ZORB.AT_5201.PV,39
S-ZORB.PT_9403.PV,3
S-ZORB.FT_9402.PV,1
S-ZORB.PDC_2502.PV,20
S-ZORB.FC_2501.PV,4
S-ZORB.SIS_LT_1001.PV,40
S-ZORB.PC_6001.PV,2
S-ZORB.AI_2903.PV,40
S-ZORB.PDC_2607.PV,2
S-ZORB.FT_1204.PV,2
S-ZORB.FT_1204.TOTAL,40
S-ZORB.PC_3101.DACA,1
S-ZORB.PT_2501.DACA,8
S-ZORB.PT_2502.DACA,10
S-ZORB.FC_2432.DACA,6
S-ZORB.FT_2431.DACA,8
S-ZORB.PDI_2801.DACA,3
S-ZORB.PDI_2301.DACA,3
S-ZORB.BS_LT_2401.PV,14
S-ZORB.PC_2401.DACA,4
S-ZORB.PC_2401B.DACA,3
S-ZORB.PC_2401B.PIDA.SP,3
S-ZORB.PC_2401B.PIDA.OP,4
S-ZORB.PC_2401.PIDA.OP,7
S-ZORB.PC_2401.PIDA.SP,6
S-ZORB.PDT_2409.DACA,2

超范围异常数据剔除结果如下表所示：

我们来看附件1的数据：

data_325_all=pd.read_excel('./源数据/325个样本数据.xlsx')
data_325_all_cao_zuo=data_325_all.iloc[:,16:]
data_325_all_cao_zuo

def check_data(data_325_all_cao_zuo,min_values,max_values):if (data_325_all_cao_zuo > max_values) or (data_325_all_cao_zuo < min_values):return np.nanelse:return data_325_all_cao_zuo
for j in range(data_range.shape[0]):names=data_range.iloc[j,1]data_min=data_range.iloc[j,6]data_max=data_range.iloc[j,7]data_325_all_cao_zuo[names]=data_325_all_cao_zuo[names].apply(lambda x:check_data(x,data_min,data_max))
data_325_all_cao_zuo

data_325_all_cao_zuo.isnull().sum()[data_325_all_cao_zuo.isnull().sum()>0]

我们得到不在范围内的数据为：

5）拉以达准则检验

设对被测量变量进行等精度测量，得到x1,x2,⋯,xnx_1\;,\;x_2\;,\;\cdots\;,\;x_nx1,x2,⋯,xn，算出其算术平均值xxx及剩余误差vi=xi−xv_i=x_i-xvi=xi−x，i=1,2,⋯,ni=1\;,\;2\;,\;\cdots\;,ni=1,2,⋯,n

按照贝塞尔公式:

σ=[1n−1∑i=1nvi2]12={[∑i=1nxi2−(∑i=1nxi)2/n]/(n−1)}12\sigma=[\frac{1}{n-1}\sum_{i=1}^nv_i^2]^{\frac{1}{2}}=\{[\sum_{i=1}^nx_i^2-(\sum_{i=1}^nx_i)^2/n]/(n-1)\}^{\frac{1}{2}} σ=[n−11i=1∑nvi2]21={[i=1∑nxi2−(i=1∑nxi)2/n]/(n−1)}21

算出标准误差σ\sigmaσ，若某个测量值xbx_bxb的剩余误差vbv_bvb且b∈[1,b]b\in[1,b]b∈[1,b]

∣vb∣=∣xb−x∣>3σ|v_b|=|x_b-x|>3\sigma ∣vb∣=∣xb−x∣>3σ

则认为xbx_bxb是含有粗大误差值的坏值，应予剔除。

因为285和313样本的非操作变量只有一条结果，所以满足3σ\sigmaσ准则，接下俩检验操作变量是否满足：

def three_sigma(df_col):"""df_col:DataFrame数据的某一列"""rule = (df_col.mean() - 3 * df_col.std() > df_col) | (df_col.mean() + 3 * df_col.std() < df_col)index = np.arange(df_col.shape[0])[rule]out_range_index=[pd.DataFrame(df_col.iloc[index]).columns,pd.DataFrame(df_col.iloc[index]).shape[0]]return out_range_index

我们首先看一下285号样本：

out_range_285_idx=[]
for i in range(data.shape[1]):df_col=data.iloc[:,i]out_range_285=three_sigma(df_col)out_range_285_idx.append(out_range_285)
out_range_285_idx

输出的结果为：

数量比较混乱，我们优化一下：

counts=0
for m in range(len(out_range_285_idx)):if out_range_285_idx[m][1]==0:counts+=1else:counts+=0
counts

说明285号样本全部都是符合3σ\sigmaσ准则的！

接下来我们看一下313号样本：

out_range_313_idx=[]
for i in range(data_1.shape[1]):df_col_1=data_1.iloc[:,i]out_range_313=three_sigma(df_col_1)out_range_313_idx.append(out_range_313)
out_range_313_idx

counts=0
for n in range(len(out_range_313_idx)):if out_range_313_idx[n][1]==0:counts+=1else:counts+=0
counts

我们的结果为：

说明我们有354-312=42条数据异常！查找出异常点的索引值和异常个数：

#%%
index_313=[]
for k in range(354):if out_range_313_idx[k][1]!=0:index_313.append((out_range_313_idx[k][0],out_range_313_idx[k][1]))
index_313

我们的输出结果为：

0,"Index(['S-ZORB.FC_2801.PV'], dtype='object')",1
1,"Index(['S-ZORB.PC_2105.PV'], dtype='object')",2
2,"Index(['S-ZORB.PT_9402.PV'], dtype='object')",1
3,"Index(['S-ZORB.FT_9401.PV'], dtype='object')",1
4,"Index(['S-ZORB.AC_6001.PV'], dtype='object')",1
5,"Index(['S-ZORB.PC_1301.PV'], dtype='object')",1
6,"Index(['S-ZORB.PT_1201.PV'], dtype='object')",2
7,"Index(['S-ZORB.FC_1202.PV'], dtype='object')",1
8,"Index(['S-ZORB.PC_1202.PV'], dtype='object')",2
9,"Index(['S-ZORB.FC_3101.PV'], dtype='object')",1
10,"Index(['S-ZORB.PDC_2607.PV'], dtype='object')",2
11,"Index(['S-ZORB.FT_1204.PV'], dtype='object')",1
12,"Index(['S-ZORB.TE_1101.DACA'], dtype='object')",1
13,"Index(['S-ZORB.TE_1107.DACA'], dtype='object')",1
14,"Index(['S-ZORB.TE_1106.DACA'], dtype='object')",1
15,"Index(['S-ZORB.LT_1501.DACA'], dtype='object')",3
16,"Index(['S-ZORB.PT_2502.DACA'], dtype='object')",2
17,"Index(['S-ZORB.PDT_2503.DACA'], dtype='object')",1
18,"Index(['S-ZORB.LT_1002.DACA'], dtype='object')",1
19,"Index(['S-ZORB.PT_1101.DACA'], dtype='object')",1
20,"Index(['S-ZORB.PT_6003.DACA'], dtype='object')",1
21,"Index(['S-ZORB.FT_3702.DACA'], dtype='object')",1
22,"Index(['S-ZORB.PT_2607.DACA'], dtype='object')",1
23,"Index(['S-ZORB.PDT_2606.DACA'], dtype='object')",2
24,"Index(['S-ZORB.ZT_2634.DACA'], dtype='object')",2
25,"Index(['S-ZORB.PC_2401B.PIDA.OP'], dtype='object')",2
26,"Index(['S-ZORB.PDT_2409.DACA'], dtype='object')",1
27,"Index(['S-ZORB.PDT_3502.DACA'], dtype='object')",1
28,"Index(['S-ZORB.PT_2106.DACA'], dtype='object')",2
29,"Index(['S-ZORB.PT_7510B.DACA'], dtype='object')",1
30,"Index(['S-ZORB.PT_7505B.DACA'], dtype='object')",1
31,"Index(['S-ZORB.PT_7107B.DACA'], dtype='object')",1
32,"Index(['S-ZORB.PT_1601.DACA'], dtype='object')",1
33,"Index(['S-ZORB.TE_5008.DACA'], dtype='object')",1
34,"Index(['S-ZORB.AT-0010.DACA.PV'], dtype='object')",1
35,"Index(['S-ZORB.AT-0011.DACA.PV'], dtype='object')",1
36,"Index(['S-ZORB.AT-0013.DACA.PV'], dtype='object')",2
37,"Index(['S-ZORB.PT_2106.DACA.PV'], dtype='object')",2
38,"Index(['S-ZORB.FT_1204.DACA.PV'], dtype='object')",2
39,"Index(['S-ZORB.TE_1106.DACA.PV'], dtype='object')",1
40,"Index(['S-ZORB.TE_1107.DACA.PV'], dtype='object')",1
41,"Index(['S-ZORB.TE_1101.DACA.PV'], dtype='object')",1

下图是拉以达准则检验异常数据结果：

3. 缺失数据的处理

在此步中，我们进一步需检查数据集中缺失值的分布情况。一般情况下，当某变量的缺失值在数据集中的占比过高时，往往说明其包含的信息过少，需对其进行删除操作。若缺失值占比较低，则可通过拉格朗日插值对数据进行补充。定义数据缺失值比率α\alphaα为：

α=位点的缺失数据个数位点的全部数据数\alpha=\frac{位点的缺失数据个数}{位点的全部数据数} α=位点的全部数据数位点的缺失数据个数

首先计算各位点数据的缺失值比率。将计算值与缺失值比率的阈值（20%）相比，按照其是否超过阈值将缺失数据分为两类：

（1）缺失值比率低的数据；

（2）数据缺失值比率高的数据。

接下来对两类数据的处理方法分别进行讨论。

1）数据缺失值比率低

对于缺失值比率较低的数据缺失情况，可认为各位点变量在较短的时间间隔内是按照线性关系变化的，故可利用数据缺失前后的数据进行两点拉格朗日插值，实现数据的补充。

设某缺失段有nnn个缺失数据，缺失段前临数据为xax_axa，后临数据为xbx_bxb，则第iii个填充数据xix_ixi（i=1,,2,,⋯,ni=1\;,\;,2\;,\;,\cdots\;,\;ni=1,,2,,⋯,n）的计算公式为：

xi=xa+xb−xan+1ix_i=x_a+\frac{x_b-x_a}{n+1}i xi=xa+n+1xb−xai

313处理的数据如下：

0,"Index(['S-ZORB.FC_2801.PV'], dtype='object')",S-ZORB.FC_2801.PV
37           892.1379
1,"Index(['S-ZORB.PC_2105.PV'], dtype='object')",S-ZORB.PC_2105.PV
20           5.701880
21           5.665528
2,"Index(['S-ZORB.PT_9402.PV'], dtype='object')",S-ZORB.PT_9402.PV
7            0.59776
3,"Index(['S-ZORB.FT_9401.PV'], dtype='object')",S-ZORB.FT_9401.PV
0           155.0762
4,"Index(['S-ZORB.AC_6001.PV'], dtype='object')",S-ZORB.AC_6001.PV
23           2.028183
5,"Index(['S-ZORB.PC_1301.PV'], dtype='object')",S-ZORB.PC_1301.PV
37           3.131695
6,"Index(['S-ZORB.PT_1201.PV'], dtype='object')",S-ZORB.PT_1201.PV
11           2.270846
38           2.248786
7,"Index(['S-ZORB.FC_1202.PV'], dtype='object')",S-ZORB.FC_1202.PV
4           45.07979
8,"Index(['S-ZORB.PC_1202.PV'], dtype='object')",S-ZORB.PC_1202.PV
11           2.240562
38           2.217983
9,"Index(['S-ZORB.FC_3101.PV'], dtype='object')",S-ZORB.FC_3101.PV
12           379.3788
10,"Index(['S-ZORB.PDC_2607.PV'], dtype='object')",S-ZORB.PDC_2607.PV
26            29.55237
32            25.89298
11,"Index(['S-ZORB.FT_1204.PV'], dtype='object')",S-ZORB.FT_1204.PV
23           6.597014
12,"Index(['S-ZORB.TE_1101.DACA'], dtype='object')",S-ZORB.TE_1101.DACA
3              132.862
13,"Index(['S-ZORB.TE_1107.DACA'], dtype='object')",S-ZORB.TE_1107.DACA
3              120.916
14,"Index(['S-ZORB.TE_1106.DACA'], dtype='object')",S-ZORB.TE_1106.DACA
3             147.8135
15,"Index(['S-ZORB.LT_1501.DACA'], dtype='object')",S-ZORB.LT_1501.DACA
5             -1.236727
17            -1.284955
25            -1.234872
16,"Index(['S-ZORB.PT_2502.DACA'], dtype='object')",S-ZORB.PT_2502.DACA
0              0.195299
33             0.199662
17,"Index(['S-ZORB.PDT_2503.DACA'], dtype='object')",S-ZORB.PDT_2503.DACA
24             -0.144057
18,"Index(['S-ZORB.LT_1002.DACA'], dtype='object')",S-ZORB.LT_1002.DACA
37             -4.99151
19,"Index(['S-ZORB.PT_1101.DACA'], dtype='object')",S-ZORB.PT_1101.DACA
37             3.124788
20,"Index(['S-ZORB.PT_6003.DACA'], dtype='object')",S-ZORB.PT_6003.DACA
26            -0.066202
21,"Index(['S-ZORB.FT_3702.DACA'], dtype='object')",S-ZORB.FT_3702.DACA
13             34.85133
22,"Index(['S-ZORB.PT_2607.DACA'], dtype='object')",S-ZORB.PT_2607.DACA
32             0.144194
23,"Index(['S-ZORB.PDT_2606.DACA'], dtype='object')",S-ZORB.PDT_2606.DACA
32              14.55408
34              16.26591
24,"Index(['S-ZORB.ZT_2634.DACA'], dtype='object')",S-ZORB.ZT_2634.DACA
26             49.16373
32             48.34832
25,"Index(['S-ZORB.PC_2401B.PIDA.OP'], dtype='object')",S-ZORB.PC_2401B.PIDA.OP
19                 33.33333
37                 33.33333
26,"Index(['S-ZORB.PDT_2409.DACA'], dtype='object')",S-ZORB.PDT_2409.DACA
4              20.64958
27,"Index(['S-ZORB.PDT_3502.DACA'], dtype='object')",S-ZORB.PDT_3502.DACA
36              16.03315
28,"Index(['S-ZORB.PT_2106.DACA'], dtype='object')",S-ZORB.PT_2106.DACA
20             5.708737
21             5.671792
29,"Index(['S-ZORB.PT_7510B.DACA'], dtype='object')",S-ZORB.PT_7510B.DACA
15              0.009181
30,"Index(['S-ZORB.PT_7505B.DACA'], dtype='object')",S-ZORB.PT_7505B.DACA
13              0.052802
31,"Index(['S-ZORB.PT_7107B.DACA'], dtype='object')",S-ZORB.PT_7107B.DACA
37              3.142406
32,"Index(['S-ZORB.PT_1601.DACA'], dtype='object')",S-ZORB.PT_1601.DACA
12             2.412431
33,"Index(['S-ZORB.TE_5008.DACA'], dtype='object')",S-ZORB.TE_5008.DACA
31             63.01546
34,"Index(['S-ZORB.AT-0010.DACA.PV'], dtype='object')",S-ZORB.AT-0010.DACA.PV
9                0.669095
35,"Index(['S-ZORB.AT-0011.DACA.PV'], dtype='object')",S-ZORB.AT-0011.DACA.PV
32                0.777408
36,"Index(['S-ZORB.AT-0013.DACA.PV'], dtype='object')",S-ZORB.AT-0013.DACA.PV
3                 0.816696
34                0.822257
37,"Index(['S-ZORB.PT_2106.DACA.PV'], dtype='object')",S-ZORB.PT_2106.DACA.PV
20                5.708737
21                5.671792
38,"Index(['S-ZORB.FT_1204.DACA.PV'], dtype='object')",S-ZORB.FT_1204.DACA.PV
24                     0.0
25                     0.0
39,"Index(['S-ZORB.TE_1106.DACA.PV'], dtype='object')",S-ZORB.TE_1106.DACA.PV
3                147.8135
40,"Index(['S-ZORB.TE_1107.DACA.PV'], dtype='object')",S-ZORB.TE_1107.DACA.PV
3                 120.916
41,"Index(['S-ZORB.TE_1101.DACA.PV'], dtype='object')",S-ZORB.TE_1101.DACA.PV
3                 132.862

2）数据缺失值比率高

如果缺失值比例高于阈值，说明其缺失的数据量信息较多，仅通过拉格朗日插值、取均值等操作都不能较好地还原正常的位点采集数据。对于这种当数据的缺失量过大而对整个数据分析造成影响的情况，应该剔除数据缺失值比例过高的位点，以减少该段数据后续的数据挖掘工作带来干扰。

缺失值比例结果如下所示：

4.样本确定

本题目标为降低 S-Zorb 装置产品辛烷值损失，故确定样本的主要依据为样品的辛烷值数据。由于辛烷值的测定数据相对于操作变量数据而言相对较少，而且辛烷值的测定往往滞后，故采取以辛烷值数据测定的时间点为基准时间，取其前 2 个小时的操作变量数据的平均值作为对应辛烷值的操作变量数据作为确定某个样本的方法。

我们综上所述，删除位点S-ZORB.AT_5201.PV、S-ZORB.SIS_LT_1001.PV、S-ZORB.AI_2903.PV以及S-ZORB.FT_1204.TOTAL

285 号和 313 号样本的辛烷值数据测定的时间点均为 8:00:00，故需取前两小时（6:00:00-8:00:00）各位点的数据的平均值作为 285 号和 313 号样本值。即本步需对进行了上述预处理操作后得到数据集中每个位点的数据取平均值并填入附件一。

第一问完工！

处理后的285和313数据：
链接：https://pan.baidu.com/s/112pGRNLEJKeRSXnb2olf2A
提取码：6666