文章目录
例题例图代码展示例题
乙醇偶合制备 C4 烯烃
C4 烯烃广泛应用于化工产品及医药的生产，乙醇是生产制备 C4 烯烃的原料。
在制备过程中，催化剂组合（即：Co 负载量、Co/SiO2 和 HAP 装料比、乙醇浓度
的组合）与温度对 C4 烯烃的选择性和 C4 烯烃收率将产生影响（名词解释见附录）。
因此通过对催化剂组合设计，探索乙醇催化偶合制备 C4 烯烃的工艺条件具有非常
重要的意义和价值。
某化工实验室针对不同催化剂在不同温度下做了一系列实验，结果如附件 1 和
附件 2 所示。请通过数学建模完成下列问题：
(1) 对附件 1 中每种催化剂组合，分别研究乙醇转化率、C4 烯烃的选择性与温
度的关系，并对附件 2 中 350 度时给定的催化剂组合在一次实验不同时间的测试结
果进行分析。
(2) 探讨不同催化剂组合及温度对乙醇转化率以及 C4 烯烃选择性大小的影响。
(3) 如何选择催化剂组合与温度，使得在相同实验条件下 C4 烯烃收率尽可能
高。若使温度低于 350 度，又如何选择催化剂组合与温度，使得 C4 烯烃收率尽可
能高。
(4) 如果允许再增加 5 次实验，应如何设计，并给出详细理由。
例图
代码展示
# -*- codeing = utf-8 -*-# @Time : /9/11 20:45# @Author : slience_me# @File : LS_9_11_1.py# @Software : 数学建模# - * - coding: utf - 8 -*-import mathimport matplotlib.pyplot as pltimport numpy as npimport pylab as plfrom pylab import mplfrom scipy.optimize import leastsqmpl.rcParams['font.sans-serif'] = ['FangSong']mpl.rcParams['axes.unicode_minus'] = Falsedef func_Ls(X_list, Y_list, Y_list_, label_list, y_label, title_list, times, row):def get_ls(Xi, Yi):def regression_assessment(y_actual_, y_predicted_):y_actual = y_actual_y_predicted = y_predicted_# mae = mean_absolute_error(y_actual, y_predicted)# mape = np.mean(np.abs((y_predicted - y_actual) / y_actual)) * 100# mse = mean_squared_error(y_actual, y_predicted)sse = np.sum((y_actual - y_predicted) ** 2)ssr = np.sum((y_predicted - np.mean(y_actual)) ** 2)sst = np.sum((y_actual - np.mean(y_actual)) ** 2)r2 = 1 - sse / sst # r2_score(y_actual, y_predicted, multioutput='raw_values')# rmse = np.sqrt(mean_squared_error(y_actual, y_predicted))# count = np.size(y_predicted)# predictionMean = np.mean(y_predicted)# yMean = np.mean(y_actual)try:r = math.sqrt(r2)except ValueError:r = np.nan# print('mae：', mae)# print('mape：', mape)# print('mse：', mse)# print('r：', r)print('r2：', r2)# print('rmse：', rmse)print('sse：', sse)print('ssr：', ssr)print('sst：', sst)# print('count：', count)# print('predictionMean：', predictionMean)# print('yMean：', yMean)return r2def func1(p, x):k, b = preturn k * x + bdef func2(p, x):a, b, c = preturn a * x ** 2 + b * x + cdef func3(p, x):k, b, c = preturn k * np.exp(-b * x) + cdef error1(p, x, y, s):print(s)return func1(p, x) - y # x、y都是列表，故返回值也是个列表def error2(p, x, y, s):print(s)return func2(p, x) - y # x、y都是列表，故返回值也是个列表def error3(p, x, y, s):print(s)return func3(p, x) - y # x、y都是列表，故返回值也是个列表# TESTp01 = [100, 2]s1 = "Test the number of iteration" # 试验最小二乘法函数leastsq得调用几次error函数才能找到使得均方误差之和最小的k、bPara1 = leastsq(error1, p01, args=(Xi, Yi, s1)) # 把error函数中除了p以外的参数打包到args中k1, b1 = Para1[0]print("k=", k1, '\n', "b=", b1)p02 = [5, 2, 10]s2 = "Test the number of iteration" # 试验最小二乘法函数leastsq得调用几次error函数才能找到使得均方误差之和最小的a~cPara2 = leastsq(error2, p02, args=(Xi, Yi, s2)) # 把error函数中除了p以外的参数打包到args中a2, b2, c2 = Para2[0]print("a=", a2, '\n', "b=", b2, "c=", c2)p03 = [0, 0, 0]s3 = "Test the number of iteration" # 试验最小二乘法函数leastsq得调用几次error函数才能找到使得均方误差之和最小的a~cPara3 = leastsq(error3, p03, args=(Xi, Yi, s3)) # 把error函数中除了p以外的参数打包到args中k3, b3, c3 = Para3[0]print("k3=", k3, '\n', "b3=", b3)y1_predicted_ = []y2_predicted_ = []y3_predicted_ = []for i in Xi:print(i)y_temp_1 = k1 * i + b1y1_predicted_.append(y_temp_1)y_temp_2 = a2 * i ** 2 + b2 * i + c2y2_predicted_.append(y_temp_2)y_temp_3 = k3 * np.exp(-b3 * i) + c3y3_predicted_.append(y_temp_3)y1_predicted_ = np.array(y1_predicted_)y2_predicted_ = np.array(y2_predicted_)y3_predicted_ = np.array(y3_predicted_)r2_1 = regression_assessment(Yi, y1_predicted_)r2_2 = regression_assessment(Yi, y2_predicted_)r2_3 = regression_assessment(Yi, y3_predicted_)data1 = (k1, b1, a2, b2, c2, k3, b3, c3)data2 = (r2_1, r2_2, r2_3)return data1, data2def draw():plt.figure(figsize=(20, 30), dpi=80)num = 0for i in range(1, times):print(num)if i % 2 == 0:fig1 = pl.subplot(row, 2, i)pl.scatter(X_list[num], Y_list[num], color="red", label="Point", linewidth=3) # 画样本点data1, data2 = get_ls(X_list[num], Y_list[num])k1, b1, a2, b2, c2, k3, b3, c3 = data1[0], data1[1], data1[2], data1[3], data1[4], data1[5], data1[6], data1[7]r2_1, r2_2, r2_3 = data2[0], data2[1], data2[2],x = np.linspace(200, 500, 1000)y1 = k1 * x + b1y2 = a2 * x ** 2 + b2 * x + c2y3 = k3 * np.exp(-b3 * x) + c3plt.plot(x, y1, label="{} 线性, 可决系数为 {}".format(label_list[num], r2_1), linewidth=2) # 画拟合直线plt.plot(x, y2, label="{} 二次, 可决系数为 {}".format(label_list[num], r2_2), linewidth=2) # 画拟合直线plt.plot(x, y3, label="{} 指数 , 可决系数为 {}".format(label_list[num], r2_3), linewidth=2) # 画拟合直线plt.xlabel("温度变化")plt.ylabel("{}".format(y_label[1]))plt.title("{}".format(label_list[num] + " " + title_list[num]))plt.legend()num += 1print("绘制图{}".format(i))else:fig2 = pl.subplot(row, 2, i)pl.scatter(X_list[num], Y_list_[num], color="red", label="Point", linewidth=3) # 画样本点data1, data2 = get_ls(X_list[num], Y_list[num])k1, b1, a2, b2, c2, k3, b3, c3 = data1[0], data1[1], data1[2], data1[3], data1[4], data1[5], data1[6], \data1[7]r2_1, r2_2, r2_3 = data2[0], data2[1], data2[2],x = np.linspace(200, 500, 1000)y1 = k1 * x + b1y2 = a2 * x ** 2 + b2 * x + c2y3 = k3 * np.exp(-b3 * x) + c3plt.plot(x, y1, label="{} 线性, 可决系数为 {}".format(label_list[num], r2_1), linewidth=2) # 画拟合直线plt.plot(x, y2, label="{} 二次, 可决系数为 {}".format(label_list[num], r2_2), linewidth=2) # 画拟合直线plt.plot(x, y3, label="{} 指数 , 可决系数为 {}".format(label_list[num], r2_3), linewidth=2) # 画拟合直线plt.xlabel("温度变化")plt.ylabel("{}".format(y_label[0]))plt.title("{}".format(label_list[num] + " " + title_list[num]))plt.legend()print("绘制图{}".format(i))pl.tight_layout() # 布局方法pl.show() # 显示方法draw()x_A1 = np.array([250, 275, 300, 325, 350])y_A1_1 = np.array([2.07, 5.85, 14.97, 19.68, 36.80])y_A1_2 = np.array([34.05, 37.43, 46.94, 49.7, 47.21])x_A2 = np.array([250, 275, 300, 325, 350])y_A2_1 = np.array([4.60, 17.20, 38.92, 56.38, 67.88])y_A2_2 = np.array([18.07, 17.28, 19.6, 30.62, 39.1])x_A3 = np.array([250, 275, 300, 325, 350, 400, 450])y_A3_1 = np.array([9.7, 19.2, 29.3, 37.6, 48.9, 83.7, 86.4])y_A3_2 = np.array([5.5, 8.04, 17.01, 28.72, 36.85, 53.43, 49.9])x_A4 = np.array([250, 275, 300, 325, 350, 400])y_A4_1 = np.array([4.0, 12.1, 29.5, 43.3, 60.5, 88.4])y_A4_2 = np.array([9.62, 8.62, 10.72, 18.89, 27.25, 41.02])x_A5 = np.array([250, 275, 300, 325, 350, 400])y_A5_1 = np.array([14.8, 12.4, 20.8, 28.3, 36.8, 76.0])y_A5_2 = np.array([1.96, 6.65, 10.12, 13.86, 18.75, 38.23])x_A6 = np.array([250, 275, 300, 350, 400])y_A6_1 = np.array([13.4, 12.8, 25.5, 55.8, 83.3])y_A6_2 = np.array([3.3, 7.1, 7.18, 10.65, 37.33])x_A7 = np.array([250, 275, 300, 350, 400])y_A7_1 = np.array([19.7, 29.0, 40.0, 58.6, 76.0])y_A7_2 = np.array([5.75, 6.56, 8.84, 18.64, 33.25])x_A8 = np.array([250, 275, 300, 350, 400])y_A8_1 = np.array([6.3, 8.8, 13.2, 31.7, 56.1])y_A8_2 = np.array([5.63, 8.52, 13.82, 25.89, 41.42])x_A9 = np.array([250, 275, 300, 350, 400])y_A9_1 = np.array([2.1, 3.0, 4.7, 13.4, 40.8])y_A9_2 = np.array([5.4, 9.68, 16.1, 31.04, 42.04])x_A10 = np.array([250, 275, 300, 350, 400])y_A10_1 = np.array([0.3, 1.0, 1.7, 9.0, 28.6])y_A10_2 = np.array([2.19, 1.65, 2.17, 3.3, 10.29])x_A11 = np.array([250, 275, 300, 350, 400])y_A11_1 = np.array([0.2, 0.5, 1.6, 8.2, 32.6])y_A11_2 = np.array([0.1, 1, 1.82, 4.35, 7.93])x_A12 = np.array([250, 275, 300, 350, 400])y_A12_1 = np.array([1.4, 3.5, 6.9, 19.9, 44.5])y_A12_2 = np.array([6.17, 8.11, 11.22, 22.26, 36.3])x_A13 = np.array([250, 275, 300, 350, 400])y_A13_1 = np.array([1.3, 2.3, 4.1, 14.6, 40.0])y_A13_2 = np.array([5.19, 7.62, 12.74, 23.46, 27.91])x_A14 = np.array([250, 275, 300, 350, 400])y_A14_1 = np.array([2.5, 5.3, 10.2, 24.0, 53.6])y_A14_2 = np.array([1.89, 2.55, 3.61, 10.83, 22.3])x_B1 = np.array([250, 275, 300, 350, 400])y_B1_1 = np.array([1.4, 3.4, 6.7, 19.3, 43.6])y_B1_2 = np.array([6.32, 8.25, 12.28, 25.97, 41.08])x_B2 = np.array([250, 275, 300, 350, 400])y_B2_1 = np.array([2.8, 4.4, 6.2, 16.2, 45.1])y_B2_2 = np.array([3.26, 4.97, 9.32, 22.88, 38.7])x_B3 = np.array([250, 275, 300, 325, 350, 400])y_B3_1 = np.array([0.4, 0.6, 1.1, 3.3, 6.0, 21.1])y_B3_2 = np.array([2.85, 5.35, 7.61, 7.74, 13.81, 21.21])x_B4 = np.array([250, 275, 300, 325, 350, 400])y_B4_1 = np.array([0.5, 1.1, 3.0, 6.1, 9.6, 33.5])y_B4_2 = np.array([6.62, 6.62, 5.05, 8.33, 13.1, 21.45])x_B5 = np.array([250, 275, 300, 325, 350, 400])y_B5_1 = np.array([2.1, 3.8, 5.8, 9.8, 15.9, 45.0])y_B5_2 = np.array([4.3, 5.06, 7.92, 11.69, 15.34, 25.83])x_B6 = np.array([250, 275, 300, 325, 350, 400])y_B6_1 = np.array([2.8, 7.5, 12.6, 15.9, 27.0, 63.2])y_B6_2 = np.array([4.5, 4.79, 8.77, 16.06, 22.41, 30.48])x_B7 = np.array([250, 275, 300, 325, 350, 400])y_B7_1 = np.array([4.4, 7.9, 11.7, 17.8, 30.2, 69.4])y_B7_2 = np.array([4.08, 6.62, 12.86, 18.45, 25.05, 38.17])X_list = [x_A1, x_A2, x_A3, x_A4, x_A5, x_A6, x_A7, x_A8, x_A9, x_A10, x_A11, x_A12, x_A13, x_A14]X_list_ = [x_B1, x_B2, x_B3, x_B4, x_B5, x_B6, x_B7]Y_list = [y_A1_1, y_A2_1, y_A3_1, y_A4_1, y_A5_1, y_A6_1, y_A7_1, y_A8_1, y_A9_1, y_A10_1, y_A11_1, y_A12_1, y_A13_1,y_A14_1]Y_list_ = [y_A1_2, y_A2_2, y_A3_2, y_A4_2, y_A5_2, y_A6_2, y_A7_2, y_A8_2, y_A9_2, y_A10_2, y_A11_2, y_A12_2, y_A13_2,y_A14_2]Y_list__ = [y_B1_1, y_B2_1, y_B3_1, y_B4_1, y_B5_1, y_B6_1, y_B7_1]Y_list___ = [y_B1_2, y_B2_2, y_B3_2, y_B4_2, y_B5_2, y_B6_2, y_B7_2]label_list = ["A1", "A2", "A3", "A4", "A5", "A6", "A7", "A8", "A9", "A10", "A11", "A12", "A13", "A14"]label_list_ = ["B1", "B2", "B3", "B4", "B5", "B6", "B7"]title_list = ["200mg 1wt%Co/SiO2- 200mg HAP-乙醇浓度1.68ml/min","200mg 2wt%Co/SiO2- 200mg HAP-乙醇浓度1.68ml/min","200mg 1wt%Co/SiO2- 200mg HAP-乙醇浓度0.9ml/min","200mg 0.5wt%Co/SiO2- 200mg HAP-乙醇浓度1.68ml/min","200mg 2wt%Co/SiO2- 200mg HAP-乙醇浓度0.3ml/min","200mg 5wt%Co/SiO2- 200mg HAP-乙醇浓度1.68ml/min","50mg 1wt%Co/SiO2- 50mg HAP-乙醇浓度0.3ml/min","50mg 1wt%Co/SiO2- 50mg HAP-乙醇浓度0.9ml/min","50mg 1wt%Co/SiO2- 50mg HAP-乙醇浓度2.1ml/min","50mg 5wt%Co/SiO2- 50mg HAP-乙醇浓度2.1ml/min","50mg 1wt%Co/SiO2+ 90mg石英砂-乙醇浓度1.68ml/min，无HAP","50mg 1wt%Co/SiO2- 50mg HAP-乙醇浓度1.68ml/min","67mg 1wt%Co/SiO2- 33mg HAP-乙醇浓度1.68ml/min","33mg 1wt%Co/SiO2- 67mg HAP-乙醇浓度1.68ml/min"]title_list_ = ["50mg 1wt%Co/SiO2- 50mg HAP-乙醇浓度1.68ml/min","100mg 1wt%Co/SiO2- 100mg HAP-乙醇浓度1.68ml/min","10mg 1wt%Co/SiO2- 10mg HAP-乙醇浓度1.68ml/min","25mg 1wt%Co/SiO2- 25mg HAP-乙醇浓度1.68ml/min","50mg 1wt%Co/SiO2- 50mg HAP-乙醇浓度2.1ml/min","75mg 1wt%Co/SiO2- 75mg HAP-乙醇浓度1.68ml/min","100mg 1wt%Co/SiO2- 100mg HAP-乙醇浓度0.9ml/min"]y_label = ["乙醇转化率变化", "C4烯烃选择性变化"]func_Ls(X_list, Y_list, Y_list_, label_list, y_label, title_list, int(29), int(14))func_Ls(X_list_, Y_list__, Y_list___, label_list_, y_label, title_list_, int(15), int(7))
