简介本资源为2021年高教社杯全国大学生数学建模竞赛B题「乙醇偶合制备C4烯烃」的二等奖完整论文面向备战数模国赛的本科生与指导教师尤其适合需要参考优秀获奖论文结构与建模思路的参赛者。压缩包内仅含1个PDF文件约3.93MB即整篇获奖论文全文涵盖摘要、问题重述、模型建立与求解、结果分析等完整章节。论文针对乙醇高效制备C4烯烃的工艺优化问题综合运用Newton插值刻画温度与乙醇转化率、C4烯烃选择性的关系以多元线性回归分析催化剂成分与温度的影响并借助BP神经网络预测C4烯烃收率再通过粒子群算法搜索最优催化剂组合与温度条件最后给出追加实验设计方案。目前已有4262人学习下载读者可从中获取完整的建模框架、公式推导与求解流程适合作为赛前研读与写作模仿的参考范本。1. 从 2021 国赛 B 题说起乙醇偶合制备 C4 烯烃到底在算什么2021 年全国大学生数学建模竞赛 B 题给了一批乙醇偶合制备 C4 烯烃的催化剂实验数据要求参赛队回答两个核心问题不同催化剂组合与温度如何影响乙醇转化率和 C4 烯烃选择性以及给定目标下如何搜索最优工艺条件。这道题当年让不少队伍在数据清洗和模型选择上翻了车因为数据里既有类别型变量催化剂组合又有连续型变量温度、装料量还夹杂着缺失值和量纲差异。如果你正在复现这道题或者手头有类似的化工工艺优化问题这篇笔记会按我实际做一遍的顺序把多元线性回归、BP 神经网络、粒子群算法和 Newton 插值这几块拼起来告诉你每一步为什么这么选、参数怎么定、哪里容易踩坑。适合已经学过 Python 基础、想拿真实赛题练手建模全流程的人。2. 数据预处理与多元线性回归基线先把能解释的部分榨干2.1 读数据与类别变量编码赛题附件通常是 Excel 或 CSV列名包含温度、催化剂组合编号、乙醇转化率、C4 烯烃选择性等。第一步不是急着上神经网络而是把数据读进来看看分布。我一般用 pandas 做快速体检重点看缺失比例和类别变量的取值个数。import pandas as pd import numpy as np # 读取附件数据sheet_name 按实际文件调整 df pd.read_excel(2021B.xlsx, sheet_name附件1) print(df.shape) print(df.isnull().sum()) # 看每列缺失数量 print(df[催化剂组合].nunique()) # 类别变量有多少种组合 print(df.describe()) # 连续变量的量纲和离群值初判逻辑说明isnull().sum()帮你决定缺失值是删还是插补nunique()决定类别变量用独热编码还是目标编码。参数上如果某个催化剂组合只有一两条样本独热编码会产生大量稀疏列这时候更适合按组合的物理成分拆成几个数值特征而不是硬编码成 21 个 0/1 列。2.2 多元线性回归作为可解释基线多元线性回归在这题里的价值不是拿高分而是给你一个可解释的下限。如果回归的 R² 只有 0.6说明温度和装料量对转化率的线性部分只能解释六成剩下的非线性交给后面的 BP 网络。做法上把类别变量做独热编码后拼进特征矩阵用 statsmodels 看每个系数的显著性。import statsmodels.api as sm # 对类别变量做独热编码drop_first 避免共线性 X pd.get_dummies(df[[温度, 装料量, 催化剂组合]], drop_firstTrue) y df[乙醇转化率] X sm.add_constant(X) # 加截距项 model sm.OLS(y, X).fit() print(model.summary()) # 重点看 P|t| 和 R-squared逻辑说明add_constant必须加否则截距被强行过原点系数会失真。drop_firstTrue是防止独热编码的虚拟变量陷阱。看结果时重点关注温度系数的符号和量级——如果温度系数为正且显著说明升温整体促进转化率这与化学直觉一致可以作为后续神经网络特征重要性的对照。参数上如果 VIF 大于 10说明类别变量之间或与温度存在共线性需要合并稀有类别。2.3 残差诊断决定要不要上非线性模型回归跑完不能只看 R²要画残差图。如果残差随温度呈现明显的 U 型或波浪形说明线性假设不成立这正是引入 BP 神经网络的信号。我一般用 matplotlib 画预测值对残差再用 Newton 插值把残差随温度的走势平滑出来看趋势。import matplotlib.pyplot as plt pred model.predict(X) resid y - pred plt.scatter(pred, resid, s8) plt.axhline(0, colorred, linewidth1) plt.xlabel(预测转化率) plt.ylabel(残差) plt.show()逻辑说明残差若在 0 附近随机分布线性模型够用若出现系统性弯曲说明温度与转化率之间存在非线性关系BP 网络的多层结构才有发挥空间。这一步是选型的分水岭不要跳过直接上神经网络否则你无法判断网络带来的提升是真实的还是过拟合。3. BP 神经网络建模结构、训练与验证集划分3.1 网络结构怎么定输入层、隐层和输出层BP 神经网络的核心是用反向传播调整权重让预测误差最小。这题输入维度取决于你编码后的特征数输出可以是乙醇转化率和 C4 烯烃选择性两个节点也可以分开建两个网络。我一般先建一个双输出网络因为两个目标共享同一组工艺条件共享隐层能捕捉共同的非线性模式。from sklearn.neural_network import MLPRegressor from sklearn.preprocessing import StandardScaler from sklearn.model_selection import train_test_split features pd.get_dummies(df[[温度, 装料量, 催化剂组合]], drop_firstTrue) targets df[[乙醇转化率, C4烯烃选择性]] scaler StandardScaler() X_scaled scaler.fit_transform(features) X_train, X_test, y_train, y_test train_test_split( X_scaled, targets, test_size0.2, random_state42 ) mlp MLPRegressor( hidden_layer_sizes(32, 16), # 两层隐层节点数递减 activationrelu, solveradam, learning_rate_init0.001, max_iter2000, early_stoppingTrue, validation_fraction0.15, random_state42 ) mlp.fit(X_train, y_train) print(mlp.score(X_test, y_test))逻辑说明hidden_layer_sizes(32, 16)是经验起点输入维度大约 20 出头时第一层 32 个神经元足够捕捉交互第二层 16 个做压缩。early_stoppingTrue配合validation_fraction0.15是防止过拟合的关键它会在验证损失连续不下降时提前停止。learning_rate_init0.001是 adam 的常用值如果损失震荡就降到 0.0005如果下降太慢就升到 0.005。标准化必须做因为温度是几百的量级装料量是个位数不缩放会导致梯度更新被大尺度特征主导。3.2 训练过程监控与超参数调整训练时不能只看最终得分要把 loss 曲线画出来。如果训练 loss 一直降但验证 loss 先降后升说明过拟合需要减层或加正则。MLPRegressor的loss_curve_属性可以直接取。import matplotlib.pyplot as plt plt.plot(mlp.loss_curve_, label训练损失) plt.xlabel(迭代轮次) plt.ylabel(损失) plt.legend() plt.show() # 如果验证损失需要单独看用 validation_fraction 的拆分手动训练逻辑说明loss_curve_只记录训练损失验证损失在 early_stopping 开启时不直接暴露所以更稳妥的做法是自己用train_test_split再切一个验证集手动循环训练并记录。参数上alpha控制 L2 正则强度默认 0.0001过拟合时调到 0.001 或 0.01batch_size默认 auto小数据集下设为 16 或 32 能让梯度更新更频繁。3.3 用交叉验证判断网络是否真的比回归好单次划分的测试得分波动很大必须做 K 折交叉验证才能下结论。我一般用 5 折比较 MLP 和线性回归的交叉验证 R²如果 MLP 只高 0.02 以内考虑到可解释性和训练成本我会优先保留回归模型。from sklearn.model_selection import cross_val_score from sklearn.linear_model import LinearRegression lr LinearRegression() mlp_cv cross_val_score(mlp, X_scaled, targets, cv5, scoringr2) lr_cv cross_val_score(lr, X_scaled, targets, cv5, scoringr2) print(MLP 交叉验证 R2:, mlp_cv.mean()) print(线性回归交叉验证 R2:, lr_cv.mean())逻辑说明cross_val_score对多输出回归默认返回每个输出的平均得分这里 targets 是两列得分是两列的平均 R²。如果 MLP 的均值高出回归 0.05 以上且标准差可控才值得在论文里作为主模型。否则把神经网络当作对照主模型仍用回归加插值这在评审眼里反而更稳。4. 粒子群算法寻优在什么空间里搜、约束怎么加4.1 粒子群算法原理与参数映射粒子群算法把每个候选工艺条件看作一个粒子粒子有位置和速度位置对应温度、装料量等决策变量速度决定下一步往哪飞。每个粒子记住自己历史最优位置群体记住全局最优位置迭代向这两个方向靠拢。在这题里决策变量是温度和装料量催化剂组合是离散的需要单独处理。import numpy as np # 决策变量温度 [200, 450]装料量 [0.5, 5.0] lb np.array([200, 0.5]) ub np.array([450, 5.0]) dim 2 n_particles 30 max_iter 100 # 初始化位置和速度 np.random.seed(42) pos lb np.random.rand(n_particles, dim) * (ub - lb) vel np.random.randn(n_particles, dim) * 0.1 pbest pos.copy() pbest_score np.full(n_particles, np.inf) gbest pos[0].copy() gbest_score np.inf逻辑说明n_particles30是中小规模搜索的常用值维度只有 2 时 20 到 40 都合理。vel初始化用标准正态乘 0.1避免初始速度过大导致粒子飞出边界。pbest_score初始化为无穷大保证第一次评估后一定更新。4.2 适应度函数把模型预测接进来适应度函数是粒子群和 BP 网络的接口。对每个粒子把位置解码成温度和装料量再拼上催化剂组合的编码送进训练好的 MLP 预测转化率和选择性按目标加权求和。如果目标是最大化 C4 烯烃收率适应度就是转化率乘选择性。def fitness(position, catalyst_code, mlp, scaler): temp, load position # 构造与训练时一致的特征向量 feat np.zeros((1, scaler.n_features_in_)) feat[0, 0] temp feat[0, 1] load # catalyst_code 对应的独热位按实际列顺序填入 feat[0, 2:] catalyst_code feat_scaled scaler.transform(feat) pred mlp.predict(feat_scaled)[0] conversion, selectivity pred[0], pred[1] return conversion * selectivity # 收率作为适应度逻辑说明特征向量的列顺序必须和训练时完全一致否则 scaler 的均值和方差会对错列。catalyst_code是预先编码好的独热向量寻优时对每种催化剂组合分别跑一次粒子群最后比较各组合的最优值。参数上如果希望兼顾转化率和选择性可以把适应度改成0.5*conversion 0.5*selectivity权重根据题目要求调整。4.3 迭代更新与边界处理粒子群的核心更新公式包括惯性项、个体认知项和社会认知项。惯性权重从 0.9 线性降到 0.4 是常见做法前期鼓励探索后期鼓励收敛。w_max, w_min 0.9, 0.4 c1, c2 1.5, 1.5 for it in range(max_iter): w w_max - (w_max - w_min) * it / max_iter for i in range(n_particles): score fitness(pos[i], catalyst_code, mlp, scaler) if score pbest_score[i]: pbest_score[i] score pbest[i] pos[i].copy() if score gbest_score: gbest_score score gbest pos[i].copy() r1, r2 np.random.rand(n_particles, dim), np.random.rand(n_particles, dim) vel (w * vel c1 * r1 * (pbest - pos) c2 * r2 * (gbest - pos)) pos pos vel # 边界截断 pos np.clip(pos, lb, ub)逻辑说明c1和c2都取 1.5 是让个体经验和群体经验权重相当如果发现早熟收敛就增大 c1 减小 c2。np.clip是最简单的边界处理比反弹法稳定但会让粒子贴在边界上如果最优解恰在边界附近需要检查是否真的到了物理极限。迭代 100 次对二维问题足够维度升高时要加到 200 以上。5. 避坑与排查这题最容易翻车的五个地方5.1 现象神经网络测试 R² 很高但寻优结果离谱原因训练集和测试集划分时没有按催化剂组合分层导致某些组合只在训练集出现网络对没见过的组合外推能力极差寻优时粒子飞到这些组合上得到虚高预测。解决用train_test_split的stratify参数按催化剂组合分层或者对每种组合单独留出测试样本。5.2 现象粒子群迭代几次就全部聚到同一点原因惯性权重下降太快或 c2 过大群体多样性丧失陷入局部最优。解决把w_min从 0.4 提到 0.5或者在速度更新后对 10% 的粒子加随机扰动强制跳出。也可以增大粒子数到 50。5.3 现象多元线性回归系数符号与化学直觉相反原因类别变量独热编码后存在完全共线性或者温度与装料量高度相关导致系数估计不稳定。解决检查 VIF删掉一个冗余的独热列或者改用岭回归Ridge(alpha1.0)让系数收缩。5.4 现象Newton 插值在数据稀疏区间剧烈震荡原因插值节点间距过大或分布不均高阶多项式在边缘产生龙格现象。解决只用 Newton 插值做趋势可视化不要用它做预测。如果必须插值改用分段三次 Hermite 插值或样条插值节点选在温度采样密集区。5.5 现象BP 网络训练损失不下降原因学习率过大导致震荡或者输入特征没有标准化或者激活函数选错。解决先确认StandardScaler已应用再把learning_rate_init降到 0.0001激活函数从 relu 换成 tanh 试一次。如果仍不降检查标签是否有 NaN。6. 用 Newton 插值做趋势验证与结果呈现的进阶技巧Newton 插值在这题里不是主力预测工具而是帮你把离散温度点上的转化率趋势平滑出来用于论文插图或验证神经网络预测是否合理。它的优势是新增节点时只需计算差商表的新一行不用重新拟合整个多项式。我一般先按温度排序取转化率列构造差商表再在密集温度网格上求值。def newton_interp(x, y, x_new): n len(x) # 计算差商表 coef np.zeros([n, n]) coef[:, 0] y for j in range(1, n): for i in range(n - j): coef[i][j] (coef[i1][j-1] - coef[i][j-1]) / (x[ij] - x[i]) # 逐项求值 result np.zeros_like(x_new, dtypefloat) for k in range(len(x_new)): term coef[0, 0] product 1.0 for j in range(1, n): product * (x_new[k] - x[j-1]) term coef[0, j] * product result[k] term return result # 按温度排序后取一组催化剂组合的数据 subset df[df[催化剂组合] 1].sort_values(温度) x_nodes subset[温度].values y_nodes subset[乙醇转化率].values x_dense np.linspace(x_nodes.min(), x_nodes.max(), 200) y_dense newton_interp(x_nodes, y_nodes, x_dense)逻辑说明差商表coef[0, j]就是 Newton 插值多项式的系数。x_new是你要评估的密集网格用于画平滑曲线。注意节点必须互异如果有重复温度点要先取平均。参数上节点数超过 8 个时高阶多项式容易震荡这时候只取温度范围中间的一段做插值边缘交给神经网络预测。验证方法是把 Newton 插值曲线和 BP 网络预测曲线画在同一张图上如果两条曲线在采样密集区重合度高说明网络学到了真实趋势如果在稀疏区分叉严重说明网络在那些区域不可信寻优结果要谨慎对待。我自己的习惯是任何寻优得到的最优温度都要回到原始数据里找最近的三个实验点看转化率是否支持这个结论。如果最优温度落在没有实验数据的空白区我会在论文里标注为模型外推并建议补充实验验证。这个习惯帮我避开了好几次把模型幻觉当成真实最优的翻车。希望帮到你。本文还有配套的精品资源点击获取