基于拉格朗日插值拟合、灰色关联度探究影响催化剂效率及其组合
摘要
C4烯烃在生活以及药物生产占领着举足轻重的地位。本文通过拉格朗日插值拟合、偏性相关分析、灰色关联度以及设置对照等方法探究在怎样的催化剂组合和温度下其利用率较高。此问题的研究可以大大提高反应速率，节约资源。
针对问题一，本题利用拉格朗日插值拟合算法，分别计算
正在上传…重新上传取消(k=0，1，2，3)，使用前四组数据插值拟合，然后使用剩余数据检验及结合模型进行数据修正。计算到3次方，即a+b++d.第二问采用偏相关分析，C4烯烃选择性与时间的t检验显著性为0.581，说明接受原假设，拒绝备择假设，呈现明显的负相关，即随着时间的变动，乙醇转化率逐渐下降。
针对问题二，利用灰色关联度探究C4 烯烃的选择性和乙醇转化率的变化趋势是否具有一致性，从而求得不同催化剂组合和温度的影响程度大小。计算后发现对于不同的催化剂组合，A7组对乙醇转化率的影响最大，关联度为0.603，B3组影响最小，关联度为0.440；A12组对乙醇转化率的影响最大，关联度为0.63，而A2组影响最小，关联度为0.439；
针对问题三，本题通过建立乘积模型，通过对照实验，划定取值区间，使用控制变量法，分别控制乙醇进气量，质量比等因素，采用单因素方差分析，算得P值为P=0.0183，效果为显著，通过实验组与对照组，最终得到最适宜的催化剂组合和温度，即：400度环境下，Co、SiO2、HAP比例为1:100:100。
针对问题四，增设五组实验，根据前三问的结论，设置1wt，Co/Sio2:HAP=1:1=200mg,前四组实验使乙醇进气量为唯一变量，分别设置为：0.1ml/min，0.6ml/min，0.12ml/min，0.18ml/min.第五组实验使Co/Sio2及HAP的质量为100mg，乙醇进气量为0.6ml/min，与第二组形成对照。
关键词：拉格朗日插值拟合、灰色关联度、双因素方差分析、对照实验
问题背景：
随着科技的日益增长，各种能源的利用不断增大，无论是工业生产还是日常生活都离不开化学能源的支撑。C4烯烃是在石油化工领域具有非常重要地位的基础原料，正是由于其地位，C烯烃被广泛用于化工、医药等许多方面的，所以，对于C4烯烃的研究重大意义以及极高的价值。
问题的提出：
在研究C4 烯烃的过程中，催化剂是必不可少的东西。因此，通过对催化剂的组合设计研究，探索使用乙醇催化偶合制备 C4烯烃的方法逐渐提上日程，现从以下几方面研究不同催化剂对制备C4烯烃的影响：
- 通过对催化剂的不同组合分别研究乙醇转化率、C4烯烃的选择性与温度的有什么样的关系。
- 在不同催化剂组合的作用下，温度对乙醇转化率的影响以及 C4烯烃选择性大小的影响。
- 在选择催化剂组合和温度的不同值，使制备C4烯烃在相同实验条件下 C4烯烃收率得倒最大化。
- 增加五组额外的实验，使C4烯烃的收率尽可能提高
二、问题分析
问题的研究对象是制备C4 烯烃的化学反应，该问题描述了不同组合催化剂在不同的温度下对结果的影响，并提出了在怎样的外部和催化条件下才能将C4 烯烃的获得率尽可能地提高。
2.1问题1
第一问：题干给出表格中不同催化剂下不同温度下乙醇转化率、C4 烯烃的选择性的数值，本题思路采用拉格朗日插值拟合，因为温度数据较少，采用此方法可以分得多个值，生成的点拟合曲线，拉格朗日插值拟合对实践中的某个物理量进行观测，在若干个不同的地方得到相应的观测值，拉格朗日插值法可以找到一个多项式，其恰好在各个观测的点取到观测到的值。这样的多项式称为拉格朗日（插值）多项式。
第二问：双因素方差分析,题干给出了350度时给定的某种催化剂组合的测试数据，在不同时间长度下，比较乙醇转化率(%)和C4烯烃选择性，根据这两因素通过分析研究中不同来源的变异对总变异的贡献大小，从而确定可控因素对研究结果影响力的大小。比较不同时间和反应结果的关联度。
2.2问题2
本题要探究在不同的外部温度下，不同催化剂组合对于C4 烯烃的选择性和乙醇转化率，需要先选择一组对照组，通过与其它组的对比进行相关性分析，最后比较出对乙醇转化率和C4烯烃选择性影响最大、最小的组，采用灰色关联度方法。在系统发展过程中，若两个因素变化的趋势具有一致性，即同步变化程度较高，即可谓二者关联程度较高；反之，则较低。
2.3问题3
第一问：由题意，选择最优的催化剂组合与温度，在相同的实验条件下，使C4烯烃的收率最高，而C4烯烃的收率等于乙醇转化率乘上C4烯烃的选择性，故而此题应使乙醇的转化率及C4烯烃的选择性尽可能高，建立乘积模型，通过对照实验，划定取值区间，通过对照优化算法，寻求最优解。
第二问：整体思路承接该小题第一小问，但需要在温度这个自变量上加以约束，及T<350，而后通过乘积模型，利用对照实验所得数据，优化关系式，得到最优解。
2.4问题4
由问题二、问题三的结果可选出相对适合的催化剂（Co/Si02和HAP装料比）组合，但对于乙醇进气量这一影响因素的探究较少。为了尽可能使C4烯烃选择性尽可能高，设置5组在较合适的相同Co/Si02和HAP装料比下，对乙醇进气量采用黄金分割法增设对照试验。
- 假设压强、海拔、湿度和当日温度等对实验结果不造成额外影响。
- 假设反应所产生的温度变化对实验设定温度没有影响。
3.假设反应物的量在称量时没有实际数值差异
三、符号说明
| 符号 | 正在上传…重新上传取消 | i | j | t | r | ξ |
| 定义 | 拉格朗日插值多项式 | 21组不同的催化剂组合 | 每组催化剂组合的不同温度 | 统计量 | 相关系数 | 相关度 |
5.1.1拉格朗日插值拟合
本题研究乙醇转化率与温度的关系，C4烯烃选择性与温度的关系，因数据量较少，故采用拉格朗日插值算法，进行科学扩充，后再将数据拟合。
拉格朗日插值算法：
正在上传…重新上传取消例：对F0--->A1组的五条数据进行拉格朗日插值拟合
通过拉格朗日插值法拟合绘图：
温度与乙醇转换率的关系
温度与烯烃选择性的关系
其他组别图像请见附录
通过对比，发现所有组别乙醇转化率以及C4烯烃选择性在一定的温度区间内，都随温度升高而升高。
最后通过算法拟合得出结果可以得出如下模型：
乙醇转化率的模型：
C4选择性的模型：
5.1.2偏性相关分析
对附件2中350度时给定的催化剂组合在实验中不同时间的测试结果进行分析。观察附件2中数据不难发现，在350度时，随着反应的进行，乙醇转化率逐渐下降，乙烯选择性、乙醛选择性和其他缓慢上升，碳数为4-12脂肪醇缓慢下降，C4烯烃选择性、甲基苯甲醛在微小的范围内进行波动，影响较小。当实验进行了273分钟时，乙醇转化率趋于稳定。
从图中可清晰看到在350度的温度环境下，各物质随时间的变化过程，在240min-270min间乙醇转化率已趋于平稳，故而推测此时乙醇偶合支配C4烯烃的反应已达到平衡。
某个指标变化时，可能会对其他指标产生影响，例如：乙醇转化率变化时，其他指标也会随之变动，因此，为了消除其他指标的影响，本题采用偏相关分析得到两指标之间的相关性。
即为计算每两个指标之间的偏相关系数，控制其余指标，除去其余变量的影响，真正反应两个指标之间的相关系数。本题中选取时间、乙醇转化率和C4烯烃选择性作为研究的对象，公式为：
式中，是变量与的简单相关系数，是变量与的简单相关系数, 是变量与的简单相关系数
利用SPSS计算不同指标之间的偏相关系数。得到各指标之间的偏相关系数表如下表5-4所示：
表5-4
| 偏相关分析 | |||
| C4烯烃选择性 | 时间 | ||
| C4烯烃选择性 | 偏相关系数 | 1 | 0.287 |
| 时间 | 显著性(双尾) | - | 0.581 |
| 偏相关分析 | |||
| 乙醇转化率 | 时间 | ||
| 乙醇转化率 | 偏相关系数 | 1 | -0.967 |
| 显著性(双尾) | - | 0.002** | |
| 时间 | 偏相关系数 | -0.967 | 1 |
| 显著性(双尾) | 0.002** | - | |
对偏相关性系数进行显著性检验，将原假设为变量之间无相关性。设统计量t：
对该统计量进行假设检验，则当显著性大于0.05时，说明接受原假设，该变量与其他变量无相关性，当显著性小于0.05时，说明拒绝原假设，该变量与其他变量有相关性。查阅相关资料，当r=1，说明变量之间完全相关，此时变量之间具有较强的正相关关系或负相关关系，r>0.8，说明变量之间具有高度相关性，r<0.3,说明变量之间有低度相关性，其余的情况变量之间具有中度相关性。
其中，C4烯烃选择性与时间的t检验显著性为0.581，说明接受原假设，则C4烯烃选择性与时间不存在明显的关系，与附件2中的数据符合，即随着时间的变动，C4烯烃选择性在较小得范围内波动。乙醇转化率与时间的t检验显著性为0.002，说明拒绝原假设，则乙醇转化率与时间存在明显的关系，又偏相关系数为-0.967，说明两个变量之间呈明显的负相关关系，与附件2中的数据符合，即随着时间的变动，乙醇转化率逐渐下降。
5.2灰色关联度方法
要探究某一变量对实验结果的影响，就需要做对比实验。
我们先确定参考组与比较组，以第一组数据为参考数列，其他二十组的数据为比较数列。
其次，对变量进行均值化处理，其中k对应各组的温度，i对应不同的催化剂组合（即一个特征）
计算关联系数
计算关联度：因为关联系数是比较数列与参考数列在不同温度下的关联程度值，所以它的数不止一个，而信息过于分散不便于进行整体性比较。因此要将同一组不同温度的关联系数集中为一个值，即求其平均值，作为比较数列与参考数列间关联程度的数量表示，关联度公式如下：
乙醇转化率关联度排序
| A7 | A13 | A1 | A5 | A3 | A14 | A11 | A6 | B1 | A12 |
| 0.603 | 0.587 | 0.58131 | 0.58052 | 0.56884 | 0.55180 | 0.55089 | 0.54535 | 0.53619 | 0.52500 |
| B5 | A8 | A2 | B6 | A4 | B4 | A9 | B3 | A10 | B2 |
| 0.5218 | 0.515 | 0.507858 | 0.503192 | 0.491804 | 0.487457 | 0.480236 | 0.448453 | 0.445291 | 0.439953 |
C4烯烃选择性关联度排序
| A12 | A10 | A1 | B2 | B6 | A8 | A3 | A13 | A11 | A14 |
| 0.632452 | 0.625382 | 0.586025 | 0.565950 | 0.551257 | 0.548447 | 0.546642 | 0.529772 | 0.528793 | 0.524696 |
| B4 | B5 | B3 | B1 | A4 | A6 | A7 | A9 | A5 | A2 |
| 0.521020 | 0.511145 | 0.507523 | 0.495391 | 0.493922 | 0.489408 | 0.481340 | 0.465057 | 0.453571 | 0.438979 |
结果分析
计算得到对于不同的催化剂组合，A7组对乙醇转化率的影响最大，关联度为0.603，B3组影响最小，关联度为0.440；A12组对乙醇转化率的影响最大，关联度为0.63，而A2组影响最小，关联度为0.439；
5.3设置对照
选择最优的催化剂组合与温度，在相同的实验条件下，使C4烯烃的收率最高，而C4烯烃的收率等于乙醇转化率乘上C4烯烃的选择性，故而此题应使乙醇的转化率及C4烯烃的选择性尽可能高。
设置对照实验，控制单一变量，如下表：
选取A7、A8、A9、A12、B1、B5，这六组实验数据，控制催化剂组合相同、乙醇进气速率相同，在不同温度、不同的质量比速率，与乙醇转化率比较分析：
使用单因素分析，p值为0.0183，结果为显著，则具有价值。
继续使用控制变量法，如图：
对数据统计后，开始第二组对照实验，控制乙醇进气量的不同，进行分析比对，使用单因素方差分析，取得p值获得结论。
Matlib单因素方差分析数据：
对数据分析可知，最合适的组合应为400度环境下，1wt，催化剂的质量比为1：1，此时既可以保证催化剂的活性保持最大，同时也可使乙醇转化率及C4乙烯的选择性最大，及C4乙烯的收率最高。温度低于350度时，其他条件不变，温度则要尽可能高。
5.4增设实验设计
增设五组实验，根据前三问的结论，设置1wt，Co/Sio2:HAP=1:1=200mg,前四组实验使乙醇进气量为唯一变量，分别设置为：0.1ml/min，0.6ml/min，0.12ml/min，0.18ml/min.第五组实验使Co/Sio2及HAP的质量为100mg，乙醇进气量为0.6ml/min，与第二组形成对照。
通过C1-C4组实验对比，可以确定乙醇进气量的最佳区间，C2与C5的对照，则可以得出质量比对实验的影响。
该五组对比试验可以对以上三问的结果做实际数据的支撑，是模型更加真是可靠。同时保证了单因素变量，便于分析其他产物对C4烯烃收率的影响。
- 吕绍沛. 乙醇偶合制备丁醇及C_4烯烃[D].大连理工大学,2018.
- 李庆杨. 数值分析. 华中科技大学出版社. 第五版
7.1模型的优点
1.基于拉格朗日插值算法的关系模型，解决的数据量较少而无法准确的预测判断的问题，准确的得出究乙醇转化率、C4 烯烃的选择性与温度的关系，结果合理，数据可靠，直观实用。
2.在灰色关联度算法算法模型中，将各个因素关联，计算关联度，并加以排序，合理预测，基于数据分析，数据真实，富有新意。
3.在基于对照实验的乘积模型中，数据以前两问为基础构建，整体性强，架构清晰，逻辑顺畅。
7.2模型的缺点
1.对乙醇偶合制备C4烯烃的反应中，未考虑除主要产物外其他产物是否会对C4烯烃的影响。
2.在模型的优化中，对一些数据通过实际情况进行了估计，存在一定的改进空间。
7.3模型的改进
1.建立副产物对主产物的影响模型，使灰色关联度模型更加精确。
2.对拉格朗日插值模型进行更高次的分析，是模型更加精确。
7.3模型的推广
在其他一些化学反应过程中，可使用拉格朗日插值模型，对数据合理扩充分析，而后使用灰色关联度模型计算关联度的强弱。
附录
附录一
实现代码
问题一
import matplotlib.pyplot as plt
from pylab import mpl
mpl.rcParams['font.sans-serif'] = ['SimHei']
mpl.rcParams['axes.unicode_minus'] = False
Q1_T = []
Q1_Z = []
Q1_C = []
list1 = [5,10,17,23,29,34,39,44,49,54,59,64,69,74,79,84,90,96,102,108,115]
list2 = [5,5,7,6,6,5,5,5,5,5,5,5,5,5,5,5,6,6,6,6,6]
group1 = []
group2 = []
group_all = []
with open ('附件1.csv','r') as f1:
lines=f1.readlines()[1:]
for line in lines:
arr = line.split(',')
if(arr[0]!=''):
group1.append(arr[0])
group2.append(arr[2]+','+arr[3]+','+arr[4])
index1 = 0
index2 = 0
t = []
c = []
z = []
for g in group2:
g_=g.split(',')
if(index1 t.append(int(g_[0])) z.append(float(g_[1])) c.append(float(g_[2])) index1+=1 if(index2==20): if(index1==114): Q1_T.append(t) Q1_Z.append(z) Q1_C.append(c) else: index2+=1 Q1_T.append(t) Q1_Z.append(z) Q1_C.append(c) t = [] z = [] c = [] t.append(int(g_[0])) z.append(float(g_[1])) c.append(float(g_[2])) index1+=1 def Parameters(data_x,data_y,size): parameters=[] i=0; while i < size: j = 0; temp = 1; while j < size: if(i != j): temp*=data_x[i]-data_x[j] j+=1; parameters.append(data_y[i]/temp) i += 1; return parameters def Calculate(data_x,parameters,x): returnValue=0 i = 0; while i < len(parameters): temp = 1 j = 0; while j< len(parameters): if(i!=j): temp *=x-data_x[j] j+=1 returnValue += temp * parameters[i] i += 1 return returnValue def Draw(data_x,data_y,new_data_x,new_data_y,n): plt.plot(new_data_x, new_data_y, label="拟合曲线", color="black") plt.scatter(data_x,data_y, label="离散数据",color="red") plt.scatter(300,n, label="真实数据", color="green") mpl.rcParams['font.sans-serif'] = ['SimHei'] mpl.rcParams['axes.unicode_minus'] = False plt.title("拉格朗日插值拟合数据") plt.legend(loc="upper left") plt.show() for i in range(21): x = Q1_T[i] y = Q1_Z[i] parameters=Parameters(x,y,list2[i]) if(list2[i]==5): if(i<4): datax=[245,255,265,275,285,295,305,315,325,335,345,355] else: datax=[245,260,275,290,305,320,335,350,365,380,395,410] elif(list2[i]==6): datax=[245,260,275,290,305,320,335,350,365,380,395,410] else: datax=[245,270,295,320,345,370,395,420,445,470] datay=[] for temp in datax: datay.append(Calculate(x,parameters,temp)) x.append(300) y.append(Calculate(x,parameters,300)) Draw(x,y,datax,datay,y[2]) Q1_T = [] Q1_Z = [] Q1_C = [] list1 = [5,10,17,23,29,34,39,44,49,54,59,64,69,74,79,84,90,96,102,108,115] group1 = [] group2 = [] group_all = [] with open ('附件1.csv','r') as f1: lines=f1.readlines()[1:] for line in lines: arr = line.split(',') if(arr[0]!=''): group1.append(arr[0]) group2.append(arr[2]+','+arr[3]+','+arr[5]) index1 = 0 index2 = 0 t = [] c = [] z = [] for g in group2: g_=g.split(',') if(index1 t.append(int(g_[0])) z.append(float(g_[1])) c.append(float(g_[2])) index1+=1 if(index2==20): if(index1==114): Q1_T.append(t) Q1_Z.append(z) Q1_C.append(c) else: index2+=1 Q1_T.append(t) Q1_Z.append(z) Q1_C.append(c) t = [] z = [] c = [] t.append(int(g_[0])) z.append(float(g_[1])) c.append(float(g_[2])) index1+=1 l = [] for i in Q1_C: count = 0 for j in i: if(count<4): l.append(j) count+=1 a1 = [-1,-3,-3,1] b1 = [900,2625,2550,-825] c1 = [-269375,-761250,-718125,226250] d1 = [26812500,73125000,67031250,-20625000] overa1 = [] overb1 = [] overc1 = [] overd1 = [] a11=0 b11 = 0 c11=0 d11=0 count = 0 for i in l: if(count<4): a11 = a11+a1[count]*i b11 = b11+b1[count]*i c11 = c11+c1[count]*i d11 = d11+d1[count]*i count+=1 else: overa1.append(a11) overb1.append(b11) overc1.append(c11) overd1.append(d11) count=0 a11=0 b11=0 c11=0 d11=0 a11 = a11+a1[count]*i b11 = b11+b1[count]*i c11 = c11+c1[count]*i d11 = d11+d1[count]*i count+=1 overa1.append(a11) overb1.append(b11) overc1.append(c11) overd1.append(d11) sum1=0 sum2=0 sum3=0 sum4=0 for i in overa1: sum1+=i for i in overb1: sum2+=i for i in overc1: sum3+=i for i in overd1: sum4+=i ave1=sum1/21 ave2=sum2/21 ave3=sum3/21 ave4=sum4/21 print(ave1/93750) print(ave2/93750) print(ave3/93750) print(ave4/93750) 问题二 import pandas as pd Q1_T = [] Q1_Z = [] Q1_C = [] group1 = [] group2 = [] list1 = [5,10,17,23,29,34,39,44,49,54,59,64,69,74,79,84,90,96,102,108,115] with open ('附件1.csv','r') as f1: lines=f1.readlines()[1:] for line in lines: arr = line.split(',') if(arr[0]!=''): group1.append(arr[0]) group2.append(arr[2]+','+arr[3]+','+arr[4]) index1 = 0 index2 = 0 t = [] c = [] z = [] for g in group2: g_=g.split(',') if(index1 t.append(int(g_[0])) z.append(float(g_[1])) c.append(float(g_[2])) index1+=1 if(index2==20): if(index1==114): Q1_T.append(t) Q1_Z.append(z) Q1_C.append(c) else: index2+=1 Q1_T.append(t) Q1_Z.append(z) Q1_C.append(c) t = [] z = [] c = [] t.append(int(g_[0])) z.append(float(g_[1])) c.append(float(g_[2])) index1+=1 x=pd.Dataframe(Q1_C) x_mean=x.mean(axis=1) for i in range(x.index.size): x.iloc[i,:] = x.iloc[i,:]/x_mean[i] ck=x.iloc[0,:] cp=x.iloc[1:,:] t=pd.Dataframe() for j in range(cp.index.size): temp=pd.Series(cp.iloc[j,:]-ck) t=t.append(temp,ignore_index=True) mmax=t.abs().max().max() mmin=t.abs().min().min() print('mmax',mmax) print('mmin',mmin) rho=0.5 ksi=((mmin+rho*mmax)/(abs(t)+rho*mmax)) r=ksi.sum(axis=1)/ksi.columns.size result=r.sort_values(ascending=False) print(result) 问题三 X=[2.07 5.85 14.97 19.68 36.80 1.4 3.5 6.9 19.9 44.5 1.4 3.4 6.7 19.3 43.6 2.8 4.4 6.2 16.2 45.1 0.4 0.6 1.1 3.3 6.0 21.1 0.5 1.1 3.0 6.1 9.6 33.5 2.8 7.5 12.6 15.9 27.0 63.2 ]; group=[ones(1,5),2.*ones(1,10),3.*ones(1,5),4.*ones(1,6),5*ones(1,6),6*ones(1,6)]; [p,anovatab,stats]=anova1(X,group);%单因素方差分析 fa=finv(0.95,anovatab{2,3},anovatab{3,3});%计算fa F=anovatab{2,5};%F值 if p<=0.01 && F>fa disp('非常显著') elseif p<=0.05 && F>fa disp('显著') else disp('不显著') end
