汽油辛烷值优化建模
汽油燃烧产生的尾气排放对大气环境有重要影响。汽油清洁化重点是降低汽油中的
硫、烯烃含量，同时尽量保持其辛烷值。然而汽油精制工艺是一个多变量的强非线性、强
耦联性的复杂系统。为控制精制过程中产品的硫含量、建立汽油辛烷值预测模型和优化操作变量降低辛烷值(RON)损失，本文充分利用线性、非线性相关分析、互信息理论、神经网络、智能寻优算法等数据挖掘技术研究了如何降低汽油精制过程中辛烷值损失的问题。
针对问题一，依据附件二给出的“样本确定方法”对附件三：285 号和 313 号样本原始数据进行数据整定。经整定，两个样本原始数据中不存在残缺数据较多、全部为空值或部分为空值的位点（操作变量），不需作数据处理；依据变量的最大最小限幅，筛选出 8个变量，剔除4个涉及样本数超过42.15%的变量；将整定好的样本替换附件一中的对应样本，同时利用拉依达准则剔除4个存在非操作变量异常值的样本。最终形成321个样本、364个变量的数据文件。
针对问题二，汽油精制工艺是一个多变量的强非线性、强耦联性的复杂系统，首先模型一利用简单线性相关分析和偏相关分析，得到产品辛烷值(RON)和原料辛烷值(RON)的强相关性，相关系数为 0.973。然后建立基于互信息的降维模型，将剩余变量依据与产品辛烷值(RON)互信息值的大小排序，剔除线性相关较强的变量后，选取前14个变量和原料辛烷值(RON)共计 15 个主要变量。模型二基于遗传算法—BP 神经网络的降维模型利用了神经网络对非线性问题的良好适应性，采用二进制编码形成自变量组合，通过 BP 神经网络计算得到产品辛烷值(RON)预测值的均方误差作为该自变量组合的适应度，经遗传算法寻优得到最优的 14 个主要变量的组合。对比两个模型，发现两组主要变量有 8 个完全相同，但模型一的剩余主要变量之间存在更弱的相关性，符合主要变量的独立性标准。
针对问题三，基于问题二模型一所得到的15个主要变量，建立模型三BP神经网络模型预测产品辛烷值(RON)。单隐含层BP神经网络预测误差在（-0.36%，0.55%）之间，双隐含层 BP 神经网络预测误差在（-0.4%，0.35%）之间。模型四小波神经网络预测模型的预测误差在（-0.15%，0.3%）。对比得知，小波神经网络的预测精度最高，但对比产品辛烷值预测值与实测值的整体趋势，发现双隐含层 BP 神经网络的预测误差分布更均匀，精度与小波神经网络预测模型接近。因此选定双隐含层BP神经网络作为预测模型。
针对问题四，基于问题三的产品辛烷值(RON)预测模型，同样建立产品硫含量的 BP神经网络预测模型。建立模型五基于 BP 神经网络预测模型的粒子群寻优算法，优化主要变量操作方案。优化过程中，除原料辛烷值(RON)的速度为 0 以外，粒子中的其他主要变辛烷值(RON)测量前两个小时内，以1次/3min的调整频次所允许的最大调整次数40次。假定合理的产品辛烷值目标区间为（0.5,0.7），最终优化得到全部 325 个样本的最优操作条件。
针对问题五，利用问题四建立的模型五，根据优化过程中存储的粒子位置矩阵和速度
矩阵，倒推得到操作变量最优解的变化轨迹，再利用问题三的 BP 神经网络预测模型得到过程中每个操作变量组合的产品辛烷值(RON)和硫含量的变化轨迹，并给出具体的操作变量轨迹。结果表明：133号样本的产品辛烷值(RON)由初始的88.09个单位提高到88.71个
关键词：互信息理论、非线性相关分析、辛烷值预测、BP神经网络、粒子群寻优算法
目录
目录 ......................................................................................................................................... 1
一、 问题重述 .................................................................................................................... 2
二、 问题分析 .................................................................................................................... 3
三、 模型假设 .................................................................................................................... 4
四、 模型建立与求解 ........................................................................................................ 5
4.1 问题一：样本数据处理 .......................................................................................... 5
4.2 问题二：寻找建模主要变量 .................................................................................. 6
4.3 问题三：产品辛烷值预测 .................................................................................... 16
4.4 问题四：主要变量操作方案的优化 .................................................................... 23
4.5 问题五：模型的可视化展示 ................................................................................ 26
五、 模型检验与评价 ...................................................................................................... 28
5.1 模型的检验 ............................................................................................................ 28
5.2 模型的优点 ............................................................................................................ 28
5.3 模型的缺点 ............................................................................................................ 28
六、 模型推广 .................................................................................................................. 29
七、 结论 .......................................................................................................................... 30
八、 参考文献 .................................................................................................................. 31
附录 ............................................................................................................................... 32
一、 问题重述
汽油是小型车辆的主要燃料，汽油燃烧产生的尾气排放对大气环境有重要影响。汽油清洁化重点是降低汽油中的硫、烯烃含量，同时尽量保持其辛烷值。为进一步提高汽油产品的品质，依据采集的数据样本通过数据挖掘技术建立辛烷值(RON)损失预测模型，给出每个样本的优化操作条件，在保证汽油产品脱硫效果的前提下，使得汽油辛烷值(RON)损失降幅在30%以上。
问题一：参考近 4 年的工业数据的预处理结果，依“样本确定方法”对 285 号和 313号数据样本进行预处理并将处理后的数据分别加入到附件一相应的样本号中，供后续研究使用。
问题二：由于催化裂化汽油精制过程是连续的，虽然操作变量每3分钟就采样一次，
但辛烷值（因变量）的测量比较麻烦，一周仅2次无法对应。但根据实际情况认为辛烷值
的测量值是测量时刻前两小时内操作变量的综合效果，因此预处理中取操作变量两小时内
的平均值与辛烷值的测量值对应，产生了325个样本。
建立降低辛烷值(RON)损失模型涉及包括7个原料性质、2个待生吸附剂性质、2个再
生吸附剂性质、2个产品性质等变量以及另外354个操作变量（共计367个变量），采用
先降维后建模的方法，有利于忽略次要因素，发现并分析影响模型的主要变量与因素。因
此，根据提供的325个样本数据，通过降维的方法从367个操作变量中合理筛选出建模主
要变量，使之尽可能具有代表性、独立性。
问题三：建立辛烷值(RON)损失预测模型。采用上述样本和建模主要变量，通过数据挖掘技术建立辛烷值(RON)损失预测模型，并进行模型验证。
| g/ | g | 的前提下， |
利用上述模型获得325个数据样本中，辛烷值(RON)损失降幅大于30%的样本对应的主要变量优化后的操作条件（优化过程中原料、待生吸附剂、再生吸附剂的性质保持不变，以它们在样本中的数据为准）。
问题五：模型的可视化展示：工业装置为了平稳生产，优化后的主要操作变量（即问题二中的主要变量）往往只能逐步调整到位，对133号样本（原料性质、待生吸附剂和再生吸附剂的性质数据保持不变，以样本中的数据为准），以图形展示其主要操作变量优化调整过程中对应的汽油辛烷值(RON)和硫含量的变化轨迹。
二、 问题分析
问题一要求对附件三的数据按照“样本确定方法”进行处理，使得处理后的数据符合
附件一的形式。本问将根据“样本确定方法”中的数据整定方法对285号和313号数据样
本进行整定，并放入附件一对应的样本号中。同时对附件一中极个别异常值所在的样本应用拉依达准则予以剔除，形成便于数据挖掘的样本数据。
问题二要求对模型涉及的 364 个变量进行降维，找到与降低产品辛烷值(RON)损失最相关的多个变量，并说明选择变量的过程及合理性。问题二中的样本变量包括原料性质、
待生吸附剂性质、再生吸附剂性质、产品性质以及操作变量，考虑模型的目标结果应为产品性质的两个变量，即产品辛烷值(RON)和硫含量。如果将产品性质与其他变量当作一个函数映射关系，那么将两者作为模型因变量，其余364个变量作为自变量，并对因变量和自变量作相关性分析，选择相关系数较大的变量作为主要变量。由于炼油工艺过程的复杂
性等原因，其操作变量与因变量是高度非线性的关系，而原料性质等其余自变量与因变量可能是线性关系。可以采用相关性分析，选择与产品辛烷值(RON)和硫含量线性相关度较高的变量；采用基于互信息的非线性相关分析，选择与产品辛烷值(RON)和硫含量非线性相关度较高的变量，并采用遗传算法进行优化选择，从而得到建模的主要变量。
问题三要求建立辛烷值(RON)损失预测模型。采用整定后的 321 个样本和建模主要变量，通过数据挖掘技术建立辛烷值(RON)损失预测模型，并进行模型验证。本问建立BP神经网络模型预测产品辛烷值(RON)。考虑单隐含层BP神经网络、双隐含层BP神经网络以及小波神经网络分别进行预测并予以比较。
| 问题四要求对主要变量的操作方案进行优化。在保证产品硫含量不大于 5.00 | g/ | g | 的 |
前提下，优化得到辛烷值(RON)损失降幅大于 30%的样本对应的主要变量优化后的操作条件。利用产品辛烷值(RON)预测模型和产品硫含量的BP神经网络预测模型，作为基于BP神经网络预测模型的粒子群寻优算法的适应度函数，优化主要变量操作方案。
问题五要求对优化模型进行可视化展示，以图形展示其主要操作变量优化调整过程中
对应的汽油辛烷值(RON)和硫含量的变化轨迹。利用问题四建立的模型，根据优化过程中
存储的粒子位置矩阵和速度矩阵，倒推得到操作变量最优解的变化轨迹，再利用问题三的BP 神经网络预测模型得到过程中每个操作变量组合的产品辛烷值(RON)和硫含量的变化轨迹，并给出具体的操作变量轨迹。
三、 模型假设
| 1、 2、 3、 4、 5、 | 假设测量数据都是真实数据，存在各种可能的误差。 假设操作变量的取值范围完全参照附件四。 | ||
| 假设任意操作变量 | ix 一次调整的幅度只能是 | 。 ix | |
| 假设任意操作变量 | ix 每次正向调整或负向调整方向的概率是50%。 | ||
| 假设操作变量按照与采集频率相同的频率进行调整，即1次/3min。 | |||
四、 模型建立与求解
4.1 问题一：样本数据处理
4.1.1 思路概述
原始数据是从实际操作过程中直接测量得到，由于装置和测量方法等存在误差、部分
位点存在问题，因此需要依据题目要求对样本数据进行处理。问题一要求参考附件一的预处理结果，以附件二所提到的“样本确定方法”中所提及的5点要求，对附件三中的2个原始样本数据进行筛选处理，并替代附件一中相应样本。
4.1.2 数据处理
首先对于附件三中的数据按照附件二“样本确定方法”进行处理，将附件三中285号的操作变量以及313号的操作变量读入MATLAB中，生成40行354列的矩阵，然后对矩阵进行操作。
第一个要求是对于残缺数据较多的位点进行删除处理。经程序搜索测量得到的数据，
并无残缺数据，虽然有部分零值，但这些零值在操作范围之内，因此是合理的。
第二个要求是删除325个样本中数据全部为空值的位点。将附件一中的数据全部读入MATLAB中生成325行368列的数据，经程序搜索，不存在空值的情况。因此，不需要对附件一中325个样本数据作删除位点的处理。
第三个要求是对于部分数据为空值的位点，空值处用其前后两个小时数据的平均值代替。实际上，由于285号以及313号样本数据无空值，因此不需要删除或使用前后两个小时数据的平均值进行替代。
第四个要求是根据工艺要求与操作经验，总结出原始数据变量的操作范围，然后采用
最大最小的限幅方法剔除一部分不在此范围的操作变量。首先依据附件四的354个操作变量取值范围得到每个操作变量的最大、最小值，对于超过限幅的数据所对应的操作变量进行剔除。将第二个要求建立的 325×368 的矩阵编入 MATLAB 代码中，按其顺序从 1~368依次进行列编号，设定最大、最小值作为筛选的范围，将超出筛选范围的样本数据挑选出
来，挑选得到的异常数据见表4.1。
表4.1 最大、最小值检验异常的样本
| 异常变量代号 | S-ZORB.AT_5201.P V | S-ZORB.SIS_LT_1 001.PV | S-ZORB.AI_2903.P V | S-ZORB.FT_1202.T OTAL |
| 异常变量标号 | 36 | 63 | 98 | 120 |
| 异常样本数 | 147 | 325 | 315 | 10 |
| 异常变量代号 | S-ZORB.FT_1204.T OTAL | S-ZORB.TE_2001. DACA | S-ZORB.AT-0012.D ACA.PV | S-ZORB.CAL.LEV EL.PV |
| 异常变量标号 | 125 | 252 | 337 | 357 |
| 异常样本数 | 137 | 1 | 9 | 1 |
从中可以看出，第36列（S-ZORB.AT_5201.PV）、第63列（S-ZORB.SIS_LT_1001.PV）、第98列（S-ZORB.AI_2903.PV）和第125列（S-ZORB.FT_1204.TOTAL）的异常数据所对应的样本数量较多，应当予以删除；而第120列（S-ZORB.FT_1202.TOTAL）、第252列（S-ZORB.TE_2001.DACA）、第 337 列（S-ZORB.AT-0012.DACA.PV）和第 357 列（S-ZORB.CAL.LEVEL.PV）的异常数据所对应的样本数量均较少，对整体样本数影响不大，可能属于部分测量误差，因此可以保留。在将有问题的操作变量剔除后，得到处理后的数组为 364 列。其具体组成是 7 个原料性质、2 个待生吸附剂性质、2 个再生吸附剂性质、3 个产品性质和 350 个操作变量，一共组成 364 列数组。值得注意的是，附件一的数据原本在285、313行已经存在数据，用整定完成的附件三的样本数据替代。
分析，计算与产品中含硫量和辛烷值的相关系数，其部分代表性结果如下所示：
表4.5 产品性质的简单相关分析
| 产品硫含量皮尔逊相关性 | 产品辛烷值皮尔逊相关性 | |
| 原料硫含量 | 0.233 | 0.495 |
| 原料辛烷值 | 0.158 | 0.973 |
| 饱和烃 | -0.327 | -0.465 |
| 烯烃 | 0.370 | 0.395 |
| 芳烃 | -0.180 | 0.093 |
| 溴值 | 0.081 | -0.129 |
| 密度 | -0.004 | 0.261 |
| … | … | … |
| 8.0MPa氢气至反吹氢压缩机出口 | -0.205 | -0.180 |
| D101原料缓冲罐压力 | -0.399 | -0.214 |
可明显看出除了产品辛烷值(RON)与原料辛烷值(RON)存在明显线性相关性，其余变量与2个因变量之间无明显线性相关。
考虑到变量间可能存在线性相关，如果不作偏相关分析，可能会造成信息的重复表达、
主要变量过多等问题。因此对原料性质、待生吸附剂和再生吸附剂性质等 11 个变量做偏
相关分析，其结果见表4.6：
表4.6 产品性质偏相关分析
| 产品硫含量 | 产品辛烷值(RON) | |||
| 偏相关性 | 显著性 | 偏相关性 | 显著性 | |
| 原料硫含量 | 0.105 | 0.064 | 0.112 | 0.048 |
| 原料辛烷值(RON) | 0.007 | 0.899 | 0.965 | 0.000 |
| 饱和烃 | 0.029 | 0.613 | -0.070 | 0.216 |
| 烯烃 | 0.031 | 0.585 | -0.069 | 0.225 |
| 芳烃 | 0.027 | 0.636 | -0.068 | 0.227 |
| 溴值 | 0.075 | 0.182 | -0.046 | 0.420 |
| 密度 | 0.038 | 0.498 | 0.044 | 0.442 |
| 待生吸附剂焦炭 | 0.091 | 0.109 | 0.151 | 0.007 |
| 待生吸附剂S | -0.035 | 0.534 | -0.102 | 0.072 |
| 再生吸附剂焦炭 | -0.065 | 0.253 | -0.220 | 0.000 |
| 再生吸附剂S | 0.058 | 0.307 | 0.176 | 0.002 |
偏相关分析计算结果表明产品辛烷值(RON)与原料辛烷值(RON)存在明显线性相关性。偏相关分析剔除了自变量之间的线性相关性，使得筛选出的变量独立性更高，更具代
表性。在简单相关性分析中，产品辛烷值与原料辛烷值的线性相关系数为 0.973，偏相关分析中，线性相关系数为0.965，且显著性系数为0，满足显著性检验，代表其具有显著代表性。因此可以将原料辛烷值(RON)作为主要变量之一。
将321个样本中的产品辛烷值(RON)与原料辛烷值(RON)绘制简单散点图（见图4.1），可发现两者之间呈现明显线性关系，因此上述计算分析是准确的，而其余自变量与因变量的相关系数较低，以线性相关度第二高的硫含量为例（见图 4.2），其线性相关程度不明显，不能通过线性分析确定其是否属于主要变量。
图4.1 产品辛烷值与原料辛烷值线性分布图
图4.2 产品辛烷值与原料硫含量线性分布图
（2）基于互信息[1]的相关分析
汽油精制过程中的变量具有高度非线性和相互强耦联性，传统的线性降维方法包括主成分分析法等无法计算变量间的非线性关系。而互信息方法是用来描述变量间相互关联大小，既可以表现变量间的线性关系，又可以反映出变量间的非线性关系。当变量间的互信息值越大时，即代表其相互耦联性越强。
式中：x是变量X 的可能取值。当x的分布高度偏向某个特定值时，即x的取值不确定性小，熵值较低；若其分布均匀则代表其不确定性较大，熵值较大。
代码
% 本程序主要用于问题一的数据处理
clc,clear;
data1=textread('11.txt');
d31=xlsread('3-1.xlsx');
d32=xlsread('3-2.xlsx');
d11=xlsread('1-1.xlsx');
dz=xlsread('zdata.xlsx');
d3leng=size(d31,2);
d1leng=size(d11,1);
d3leng2=size(d31,1);
dele=[];
…………
% 将285数据以及313数据存储至附件1数据中
wu=zeros(2,d3leng);
for i=1:d3leng
wu(1,i)=d11(285,i)-chuli(i);
wu(2,i)=d11(313,i)-chuli2(i);
end
for i=1:d3leng
if ~isnan(chuli(i))
d11(285,i)=chuli(i);
end
if ~isnan(chuli2(i))
d11(313,i)=chuli2(i);
end
end
[m1 n1]=size(dz);
% 将非操作变量进行3sigema处理
cun1=[];
sig=[];
xba1=[];
des=zeros(m1,14);
for i=1:14
sum=0;
num=0;
cun=[];
sum2=0;
znum=0;
zba=0;
sum3=0;
for j=1:m1
sum=sum+dz(j,i);
num=num+1;
cun=[cun dz(j,i)];
end
xba=sum/num;
xba1=[xba1 xba];
cun=cun-xba;
for k=1:num
sum2=sum2+cun(k)^2;
end
sigema=sqrt(sum2/(num-1));
sig=[sig sigema];
for j=1:m1
if abs(cun(j))<=3*sigema
des(j,i)=1;
end
end
end
% 将关于原料辛烷值、产品辛烷值、和辛烷损失值超过3sigema的样本删除 del1(:,1)=des(:,2);
del1(:,2)=des(:,9);
del1(:,3)=des(:,10);
del2(:)=del1(:,1)+del1(:,2)+del1(:,3);
j=0;
cc1=[];
for i=1:d1leng
if del2(i)==3
j=j+1;
aa(j,:)=dz(i,:);
else
cc1=[cc1 i];
end
end
% 将数据异常的22、49、84、111列删除
bb=aa;
bb(:,[36 63 98 125])=[];
% 输出整定后的数据
xlswrite('321.xlsx',bb);
问题二：模型一主要程序
% 本程序用于挑选计算所用的变量
clc,clear
ppz=xlsread('321.xlsx');
pp2=xlsread('bzh.xlsx');
% 将数据进行从新处理，将产品辛烷值、硫、损伤值作为1,2,3,存储于pp1中 m1=size(ppz,1);
m2=size(ppz,2);
pp1(:,1)=ppz(:,9);
pp1(:,2)=ppz(:,8);
pp1(:,3)=ppz(:,10);
pp1(:,4:10)=ppz(:,1:7);
pp1(:,11:m2)=ppz(:,11:m2);
% 选择三个线性先关系数强的变量
cc=corrcoef(pp1);
dd=abs(cc(1,:));
nn1=100;
maxx=[];
for i=1:nn1
[m,p1]=max(dd);
maxx=[maxx; m p1];
dd(p1)=0;
end
% 保证选择的变量之间非强相关
% 数据标准化处理
pp=zeros(m1,m2-3);
for i=1:m2
pp(:,i)=zscore(pp1(:,i));
end
p=pp(:,4:m2);
t1=pp1(:,1);
t2=pp1(:,2);
t3=pp1(:,3);
t1=(t1-min(t1))/(max(t1)-min(t1));%zscore(t1);
t2=zscore(t2);
t3=zscore(t3);
u1 = t1;
for k=4:m2
%p(:,k)=zscore(p(:,k));
u2 = pp1(:,k);
wind_size = size(u1,1);
mi = calmi(u1, u2, wind_size);
MIC(k)=mi;
end
mic2=MIC;
% 寻找最大值
max1=[];
nu1=m2-3;
for i=1:nu1
[m,p1]=max(mic2);
max1=[max1; m p1];
mic2(p1)=0;
end
% 按照互相关系数大小输出数据
xlswrite('max1.xlsx',max1);
% 互相关数据的关联性
for i=1:m2-3
gl(:,i)=pp1(:,max1(i,2));
end
gl1=corrcoef(gl);
xlswrite('guanlian.xlsx',gl1);
nu=90; %数据个数
a1=zeros(m1,nu);
q1=1;
q2=2;
max2=max1(1,:);
aa=pp1(:,max2(1,2));
a1=[];
for i=1:nu
j=0;
pan1=1;
while pan1>=0.8
j=j+1;
q2=q1+j;
s1=pp(:,max1(q1,2));
bb=aa;
s2=pp(:,max1(q2,2));
bb=[bb s2];
length=size(bb,2);
des=corrcoef(bb);
% 使得变量之间的相关系数尽可能大
pan = corrcoef(s1,s2);
% pan1=abs(pan(1,2));
pan1=0;
for k=1:length-1
if des(length,k)>0.5
pan1=1+pan1;
end
end
end
q1=q2;
max2=[max2;max1(q2,1) max1(q2,2)];
aa=[aa pp1(:,max1(q2,2))];
end
xlswrite('max22.xlsx',max2);
% 将样本数据输出
for i=1:nu
a1(:,i)=pp1(:,max2(i,2));
扫一扫
1260
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
