2021国赛B题：乙醇偶合制备C4烯烃的BP神经网络与粒子群寻优

简介本资源为2021年高教社杯全国大学生数学建模竞赛B题「乙醇偶合制备C4烯烃」的二等奖完整论文面向备战数模国赛的本科生与建模爱好者尤其适合需要参考优秀获奖论文结构与建模思路的参赛者。压缩包内仅含1个PDF文件大小约3.93MB即论文全文涵盖摘要、问题重述、模型建立与求解及附录代码便于直接阅读与对照学习。目前已有4262人学习下载热度较高。论文围绕乙醇高效制备C4烯烃的工艺条件展开针对四个子问题分别采用Newton插值刻画温度与乙醇转化率、C4烯烃选择性的关系建立多元线性回归模型分析催化剂成分与温度的影响并引入BP神经网络结合粒子群算法PSO优化C4烯烃收率最终给出最优催化剂组合与温度区间并设计补充实验方案。读者可从中获取完整的建模框架、算法选型依据、检验方法与论文写作范式是数模国赛备赛的实用参考资料。1. 从 2021 国赛 B 题说起乙醇偶合制备 C4 烯烃到底在算什么2021 年全国大学生数学建模竞赛 B 题给了一批乙醇偶合制备 C4 烯烃的实验数据核心诉求就一句话在不同催化剂组合和温度下C4 烯烃收率怎么变怎么找到最优工艺条件。做过这道题的人都知道它表面是化工题骨子里是数据建模题——温度、乙醇浓度、催化剂配比这些自变量和乙烯、C4 烯烃收率这些因变量之间既不是简单线性关系也不是纯黑箱能糊弄过去的。这道题当年拿二等奖的队伍通常不是靠堆模型复杂度赢的而是把数据清洗、插值补点、回归拟合、非线性寻优这条链路走扎实了。我后来复盘过很多次发现真正拉开差距的是三件事温度区间内缺失点的处理方式、回归模型对非线性段的拟合能力、以及寻优时参数边界设得合不合理。这篇笔记就按这条链路拆开讲从数据预处理到 BP 神经网络拟合再到粒子群算法找最优条件每一步都给可复现的代码和参数说明。适合正在做化工数据建模、或者想拿这道题练手的同学也适合想搞清楚 BP 神经网络和粒子群算法怎么在真实数据上落地的人。2. 数据预处理与 Newton 插值把散点补成能用的曲线2.1 先搞清楚原始数据长什么样题目给的附件数据一般是按催化剂组合分组的每组里有温度、乙醇转化率、C4 烯烃选择性、C4 烯烃收率这几列。收率等于转化率乘选择性这个关系要先验证一遍如果对不上说明数据里有需要处理的异常值。我一般会先做三件事检查缺失值、看每个催化剂组合下的温度覆盖范围、画散点图看趋势。import pandas as pd import numpy as np import matplotlib.pyplot as plt # 读取附件数据假设文件名为 data.xlsx df pd.read_excel(data.xlsx) # 检查缺失值和基本统计 print(df.isnull().sum()) print(df.describe()) # 验证收率 转化率 * 选择性 df[check] df[乙醇转化率] * df[C4烯烃选择性] / 100 print((df[check] - df[C4烯烃收率]).abs().max()) # 按催化剂组合分组画散点 for name, group in df.groupby(催化剂组合): plt.scatter(group[温度], group[C4烯烃收率], labelname) plt.xlabel(温度) plt.ylabel(C4烯烃收率) plt.legend() plt.show()这段代码的关键在第三步如果check和实际收率差得超过 1 个百分点说明数据录入或计算口径有问题得回去核对。分组画散点是为了看每个组合的温度点是不是均匀分布很多队伍在这一步就发现某些组合只有三四个温度点直接拟合会严重过拟合。2.2 Newton 插值补点的适用边界温度点稀疏的时候常见做法是用插值把曲线补密。Newton 插值适合等距节点计算量小但有个坑高次插值会出现龙格现象两端震荡得厉害。我一般把插值阶数控制在 3 到 5 之间超过 5 就改用分段三次样条。下面是一个 Newton 插值的实现输入是温度数组和对应的收率数组输出是加密后的温度-收率对。def newton_interpolation(x, y, x_new): x, y: 原始节点 x_new: 待插值点 n len(x) # 计算差商表 diff np.zeros((n, n)) diff[:, 0] y for j in range(1, n): for i in range(n - j): diff[i][j] (diff[i1][j-1] - diff[i][j-1]) / (x[ij] - x[i]) # 计算插值结果 result np.zeros_like(x_new, dtypefloat) for k, xv in enumerate(x_new): val diff[0][0] term 1.0 for i in range(1, n): term * (xv - x[i-1]) val diff[0][i] * term result[k] val return result # 对某一组数据插值 group df[df[催化剂组合] A1].sort_values(温度) x group[温度].values y group[C4烯烃收率].values x_new np.linspace(x.min(), x.max(), 50) y_new newton_interpolation(x, y, x_new)差商表是 Newton 插值的核心diff[0][i]就是第 i 阶差商。参数上要注意x必须严格递增否则差商计算会除零x_new的范围不能超出原始温度区间外推段没有物理意义。插值完一定要画图对比原始点和插值曲线如果曲线在两端翘得离谱就说明阶数太高了降到 3 阶再试。2.3 插值之后还要做的一件事插值补出来的点不能直接拿去训练模型因为插值点之间是强相关的会让回归模型高估自己的拟合能力。我的习惯是把插值点只用于画趋势图和确定温度区间真正训练 BP 神经网络时还是用原始实验点或者用插值点做交叉验证的补充。另外如果某个催化剂组合的温度范围和其他组合差太多建议单独建模不要混在一起训练否则温度这个变量会被稀释掉。3. 多元线性回归打底先知道线性部分能解释多少3.1 回归模型怎么设在上一章把数据补密、趋势看清楚之后下一步不是直接上神经网络而是先用多元线性回归探底。原因很简单如果线性模型就能解释 80% 以上的方差那非线性模型的提升空间有限没必要把问题搞复杂。我一般会把温度、乙醇浓度、催化剂配比作为自变量C4 烯烃收率作为因变量先跑一个基准回归。import statsmodels.api as sm # 构造自变量矩阵温度做中心化处理 X df[[温度, 乙醇浓度, 催化剂配比]].copy() X[温度] X[温度] - X[温度].mean() X sm.add_constant(X) y df[C4烯烃收率] model sm.OLS(y, X).fit() print(model.summary())中心化处理是为了让截距项有物理意义不然截距就是温度为零时的收率没有参考价值。summary()里重点看三个数R-squared、各变量的 p 值、以及残差的正态性检验。如果某个变量 p 值大于 0.05说明它对收率的影响不显著可以考虑去掉或者换成非线性项。3.2 什么时候该加交互项和平方项温度对收率的影响通常不是线性的低温段收率随温度上升快高温段可能反而下降。这时候要在回归里加温度的平方项甚至温度和其他变量的交互项。我一般会先画收率对温度的散点如果明显是个倒 U 型就加温度^2如果不同催化剂组合的曲线斜率不一样就加温度 * 催化剂配比交互项。# 加入平方项和交互项 X2 df[[温度, 乙醇浓度, 催化剂配比]].copy() X2[温度] X2[温度] - X2[温度].mean() X2[温度平方] X2[温度] ** 2 X2[温度_催化剂] X2[温度] * X2[催化剂配比] X2 sm.add_constant(X2) model2 sm.OLS(y, X2).fit() print(model2.summary())加完平方项后 R-squared 一般会涨但如果涨得太多而样本量又小就要警惕过拟合。判断标准是调整后的 R-squared 有没有同步提升如果调整 R-squared 反而降了说明加的项不值得。3.3 回归残差告诉你的信息回归跑完不要只看 R-squared残差图才是最有信息量的。把预测值做横轴、残差做纵轴如果残差呈现喇叭口形状说明存在异方差这时候要么对因变量做变换要么改用加权最小二乘。如果残差在某个温度区间系统性偏正或偏负说明线性模型在这个区间失效了这正是后面 BP 神经网络要补的地方。我一般会把残差绝对值大于两倍标准差的点标出来回去核对原始数据是不是记录错了。4. BP 神经网络拟合结构、参数和训练技巧4.1 网络结构怎么定BP 神经网络做函数拟合结构不用太深。这道题的自变量一般不超过 5 个我通常用一层隐藏层神经元个数在 8 到 15 之间试。隐藏层激活函数用 tanh 或 relu输出层用线性激活因为收率是连续值。下面是一个用 PyTorch 搭的最小可用网络。import torch import torch.nn as nn import torch.optim as optim class BPNet(nn.Module): def __init__(self, input_dim, hidden_dim): super(BPNet, self).__init__() self.fc1 nn.Linear(input_dim, hidden_dim) self.relu nn.ReLU() self.fc2 nn.Linear(hidden_dim, 1) def forward(self, x): x self.relu(self.fc1(x)) x self.fc2(x) return x # 假设输入是温度、乙醇浓度、催化剂配比三个特征 input_dim 3 hidden_dim 12 net BPNet(input_dim, hidden_dim) criterion nn.MSELoss() optimizer optim.Adam(net.parameters(), lr0.01)隐藏层神经元个数不是越多越好12 个是我在这类化工数据上比较常用的起点。如果训练集 loss 降不下去先加神经元如果训练集 loss 很低但验证集 loss 高说明过拟合要减神经元或者加正则化。4.2 训练循环和早停训练的时候要把数据分成训练集和验证集比例大概 8:2。每轮记录训练 loss 和验证 loss验证 loss 连续 20 轮不下降就早停防止过拟合。from sklearn.model_selection import train_test_split from sklearn.preprocessing import StandardScaler # 准备数据 X df[[温度, 乙醇浓度, 催化剂配比]].values y df[C4烯烃收率].values.reshape(-1, 1) scaler StandardScaler() X scaler.fit_transform(X) X_train, X_val, y_train, y_val train_test_split(X, y, test_size0.2, random_state42) X_train torch.FloatTensor(X_train) y_train torch.FloatTensor(y_train) X_val torch.FloatTensor(X_val) y_val torch.FloatTensor(y_val) best_val_loss float(inf) patience 20 counter 0 for epoch in range(1000): net.train() optimizer.zero_grad() output net(X_train) loss criterion(output, y_train) loss.backward() optimizer.step() net.eval() with torch.no_grad(): val_output net(X_val) val_loss criterion(val_output, y_val) if val_loss best_val_loss: best_val_loss val_loss counter 0 torch.save(net.state_dict(), best_model.pth) else: counter 1 if counter patience: print(fEarly stop at epoch {epoch}) break标准化是必须的温度数值在几百乙醇浓度在几十不标准化的话梯度下降会非常慢。早停的 patience 设 20 是我试出来的经验值太小容易停在局部最优太大浪费时间。4.3 学习率调度和批量大小的选择学习率一开始设 0.01如果 loss 震荡就降到 0.001。批量大小在这类小数据集上直接用全批量就行没必要搞 mini-batch。如果数据量超过一千条可以用 32 或 64 的批量。另外Adam 优化器对学习率不敏感但如果你换成 SGD学习率要调到 0.001 以下不然很容易发散。5. 粒子群算法寻优在拟合面上找最优工艺条件5.1 粒子群算法的参数怎么设BP 神经网络训练好之后它就是一个可调用的函数输入温度、乙醇浓度、催化剂配比输出预测收率。粒子群算法的任务是在这些自变量的取值范围内找到让预测收率最大的组合。粒子群的核心参数有三个粒子数、惯性权重、学习因子。我一般设粒子数 30惯性权重从 0.9 线性降到 0.4学习因子都设 2.0。import numpy as np def pso_optimize(model, scaler, bounds, num_particles30, max_iter100): model: 训练好的 BP 网络 scaler: 标准化器 bounds: 每个自变量的取值范围列表形式 [(min, max), ...] dim len(bounds) # 初始化粒子位置和速度 positions np.random.uniform( low[b[0] for b in bounds], high[b[1] for b in bounds], size(num_particles, dim) ) velocities np.random.uniform(-1, 1, size(num_particles, dim)) # 个体最优和全局最优 pbest positions.copy() pbest_score np.full(num_particles, -np.inf) gbest positions[0].copy() gbest_score -np.inf for iteration in range(max_iter): # 惯性权重线性递减 w 0.9 - 0.5 * iteration / max_iter for i in range(num_particles): # 预测收率 x_scaled scaler.transform(positions[i].reshape(1, -1)) x_tensor torch.FloatTensor(x_scaled) model.eval() with torch.no_grad(): score model(x_tensor).item() if score pbest_score[i]: pbest_score[i] score pbest[i] positions[i].copy() if score gbest_score: gbest_score score gbest positions[i].copy() # 更新速度和位置 r1 np.random.rand(num_particles, dim) r2 np.random.rand(num_particles, dim) velocities (w * velocities 2.0 * r1 * (pbest - positions) 2.0 * r2 * (gbest - positions)) positions positions velocities # 边界处理 for d in range(dim): positions[:, d] np.clip(positions[:, d], bounds[d][0], bounds[d][1]) return gbest, gbest_score惯性权重从 0.9 降到 0.4 是为了前期探索、后期收敛。学习因子 2.0 是经典取值调大容易早熟调小收敛慢。边界处理用np.clip直接截断比反射边界简单效果也够用。5.2 寻优结果怎么验证粒子群找到的最优组合不能直接信要做两件事验证。第一把这个组合代回原始实验数据附近看有没有实际实验点支持如果最优温度落在两个实验点中间要说明这是插值预测的结果。第二换一组随机种子重新跑粒子群如果两次结果差很多说明拟合面不平滑需要回去检查 BP 网络的训练质量。# 假设 bounds 是 [(200, 400), (0.5, 2.0), (1, 5)] bounds [(200, 400), (0.5, 2.0), (1, 5)] best_pos, best_score pso_optimize(net, scaler, bounds) print(f最优条件: 温度{best_pos[0]:.1f}, 乙醇浓度{best_pos[1]:.2f}, 催化剂配比{best_pos[2]:.2f}) print(f预测收率: {best_score:.2f})如果最优温度贴着边界比如正好是 400说明真实最优可能在边界外需要扩大搜索范围重新跑。如果最优收率比所有实验点的收率都高很多要警惕过拟合导致的虚高。6. 避坑与排查这道题最容易翻车的五个地方6.1 插值点混入训练集导致 R-squared 虚高现象回归模型 R-squared 跑到 0.99但预测新数据时误差很大。原因把 Newton 插值补出来的点当成真实实验点放进训练集插值点之间强相关模型相当于在背答案。解决训练集只用原始实验点插值点只用于画图和确定温度区间验证集从原始点里划。6.2 BP 网络不标准化直接训练现象loss 一直不降或者降得很慢。原因温度在 200 到 400 之间乙醇浓度在 0.5 到 2 之间量纲差两个数量级梯度下降被大数值特征主导。解决训练前用StandardScaler对自变量做标准化预测时记得用同一个 scaler 做逆变换。6.3 粒子群早熟收敛到局部最优现象粒子群跑了几十轮就不动了找到的最优解明显不是全局最优。原因惯性权重降得太快或者粒子数太少。解决惯性权重从 0.9 降到 0.4粒子数加到 30 以上如果还不行就加变异操作每轮随机重置 10% 的粒子位置。6.4 温度边界设得太窄现象最优温度贴着搜索边界。原因初始 bounds 是根据实验数据的最小最大值设的但真实最优可能在实验范围之外。解决把温度边界往外扩 10% 到 20%重新跑粒子群如果最优还在边界上继续扩直到最优落在区间内部。6.5 忽略催化剂组合的类别差异现象所有催化剂组合混在一起建模预测误差大。原因不同催化剂组合的反应机理不一样温度-收率曲线形状不同混在一起模型学不到统一规律。解决按催化剂组合分组建模或者把催化剂组合做 one-hot 编码作为额外输入特征。7. 进阶技巧用交叉验证选隐藏层神经元个数BP 神经网络的隐藏层神经元个数是这道题里最玄学的参数。我一开始靠试后来改成用 5 折交叉验证来选。具体做法是对每个候选神经元个数8、10、12、15、20跑 5 折交叉验证取验证集 MSE 的平均值选最小的那个。下面是一个可复用的代码框架。from sklearn.model_selection import KFold def cv_select_hidden(X, y, hidden_list, epochs500): kf KFold(n_splits5, shuffleTrue, random_state42) results {} for hidden_dim in hidden_list: mse_list [] for train_idx, val_idx in kf.split(X): X_tr, X_val X[train_idx], X[val_idx] y_tr, y_val y[train_idx], y[val_idx] net BPNet(X.shape[1], hidden_dim) optimizer optim.Adam(net.parameters(), lr0.01) criterion nn.MSELoss() X_tr_t torch.FloatTensor(X_tr) y_tr_t torch.FloatTensor(y_tr) X_val_t torch.FloatTensor(X_val) y_val_t torch.FloatTensor(y_val) for epoch in range(epochs): net.train() optimizer.zero_grad() loss criterion(net(X_tr_t), y_tr_t) loss.backward() optimizer.step() net.eval() with torch.no_grad(): val_pred net(X_val_t) mse criterion(val_pred, y_val_t).item() mse_list.append(mse) results[hidden_dim] np.mean(mse_list) return results # 使用 hidden_list [8, 10, 12, 15, 20] results cv_select_hidden(X, y, hidden_list) best_hidden min(results, keyresults.get) print(f最佳隐藏层神经元个数: {best_hidden})这个框架的关键是每折都重新初始化网络不能用同一个网络跑五折否则验证集的信息会泄漏到训练里。另外epochs设 500 是折中值如果数据量大可以降到 200数据量小可以加到 1000。跑完交叉验证后用最佳神经元个数在全量数据上重新训练一次作为最终模型。我自己的习惯是交叉验证选出来的神经元个数再手动加 2 到 3 个作为最终值因为交叉验证是在子集上跑的全量数据上稍微大一点更稳。这个技巧帮我在这道题上把验证集 MSE 降了大概 15%比盲目试参数靠谱得多。希望帮到你。本文还有配套的精品资源点击获取