简介面向全国大学生数学建模竞赛B题参赛者这份资源是一篇获得二等奖的完整论文聚焦乙醇偶合制备C4烯烃的工艺优化问题。论文以牛顿插值、多元线性回归、BP神经网络与粒子群算法为建模主线完整覆盖四个问题先用牛顿插值刻画温度与乙醇转化率及C4烯烃选择性的关系再用多元线性回归分析催化剂组分与温度的影响并对回归结果进行F检验与T检验继而以BP神经网络建立C4烯烃收率预测模型借助粒子群算法搜索最优催化剂组合与温度最后给出允许追加五次实验时的补充设计方案附录代码齐全便于复现关键模型与结果。压缩包内含一个PDF文件大小约三点九三兆字节适合竞赛备赛、建模方法研习或相关课题参考该资源已有四千二百六十四人学习浏览二等奖方案在算法组合、论文结构与图表呈现上均有示范性尤其适合希望学习插值拟合与智能优化结合应用的读者。1. 2021 年数模国赛B题「乙醇偶合制备 C4 烯烃」这道题到底在考什么2021 年数模国赛 B 题「乙醇偶合制备 C4 烯烃问题」表面挂的是化学工程的名头实际是一道标准的实验数据回归与寻优题题目给出不同催化剂组合和温度下的乙醇转化率、C4 烯烃选择性数据要求你回答三件事——转化率随催化剂组合怎么变、选择性随温度怎么变、最优工艺条件落在哪个区间。拿二等奖的关键不是催化机理背得多熟而是能不能把「催化剂组合」这种实验编号拆成可计算的物理变量再保证交叉验证不虚高。适合有 Python 基础、想稳定拿奖的队伍即使没学过催化把数据切分、特征设计和寻优边界做扎实就能超过大多数参赛组。这类题最常翻车的点是数据切分方式错误导致验证结果好看、实际外推失灵后面几章我会把这一条单独拿出来讲透。2. 破题先从数据表开始变量拆解、评估口径和建模选型2.1 附件数据里真正能用的变量把催化剂编号还原成物理量打开附件 1 之后表格结构其实非常规整每一行是一次反应结果包含催化剂组合编号、温度、乙醇转化率、C4 烯烃选择性。真正麻烦的是「催化剂组合」这一列它不是一个数而是一组实验条件的代称。翻附录说明就能看到每个编号背后至少对应四样东西钴负载量、HAP 负载量、Co/SiO2 与 HAP 的装填质量比以及装填方式是混合还是分层。C 组实验还会额外改变乙醇浓度。原始字段拆解后的建模字段类型作用催化剂组合Co_load钴负载量数值影响脱氢/偶合活性位数量催化剂组合HAP_loadHAP 负载量数值影响酸碱性与分散度催化剂组合mass_ratioCo/SiO2 与 HAP 质量比数值核心可调工艺变量催化剂组合mode_mix混合/分层0/1影响床层接触方式温度T数值最主要的反应条件乙醇浓度ethanol_conc数值C 组中变化A/B 组固定我一般拿到数据后的第一步不是建模而是写一个函数把「催化剂组合」列解析成上面这张表里的六列并和原来的温度、转化率、选择性拼成一张新表。这里有个原则能做数值化的不要用类别编码能拆成物理量的不要保留编号。因为后面寻优是要外推到数据里没有出现过的催化剂组合的——如果模型里只有编号外推这一步就完全没法做。2.2 评估口径先定住GroupKFold 按催化剂组合切分这道题数据规模不大但并不是普通回归里的独立同分布样本。同一个催化剂组合的温度点是一组实验做出来的它们共享这个组合的特性如果不做处理直接随机切分训练集和验证集同一个组合的温度点会被同时分到两边模型在验证时相当于已经见过这个组合的其它温度点R² 会虚高到不真实。from sklearn.model_selection import GroupKFold # 关键groups 传催化剂组合编号而不是行索引 groups df[catalyst_id] gkf GroupKFold(n_splits5) for train_idx, val_idx in gkf.split(X, y, groupsgroups): X_train, X_val X.iloc[train_idx], X.iloc[val_idx] y_train, y_val y.iloc[train_idx], y.iloc[val_idx] # 在这里训练模型并记录验证集分数这段代码的含义是把每一个催化剂组合当作一个整体要么全部进训练集要么全部进验证集。这样验证集里的数据点是模型从未见过的「新催化剂组合」评估出的分数才代表真实外推能力。注意gkf.split里第一个参数是特征矩阵第二个是目标值第三个groups必须是按行顺序排列的组编号。提示GroupKFold 是这道题交叉验证的底线。如果你用了普通 KFold后面的模型调参、寻优全部建立在虚高分数上评审一换数据就露馅。2.3 建模选型对比为什么首选多项式回归而不是深度学习选模型之前先明确三个约束数据量只有一两百行样本远小于特征可能数需要外推到温度区间内部但没做过实验的点论文里要能解释每个自变量对结果的影响方向。把这三个约束摆出来深度神经网络第一个排除——一两百条数据训练神经网络基本是过拟合到噪声随机森林和梯度提升树拟合能力强但对实验范围外的外推无能为力因为它们本质上是分段常数预测SVR 可以做非线性拟合但核函数参数调起来很玄学系数解释更是黑匣子。方案优势风险定位多元线性回归稳定、可解释欠拟合严重基线模型多项式回归 Ridge可解释、能表达弯曲与交互阶数过高会过拟合主模型SVR非线性拟合强核参数难调、外推不稳对比模型随机森林 / 梯度提升非线性强、不用做特征工程无法外推、曲面不平滑对比模型动力学机理模型物理意义明确参数多、数据不足难辨识不推荐神经网络表达力强样本量太小、不可解释不推荐所以整篇建模的主线是温度项用多项式保留物理趋势催化剂变量做成编码和显式交互项最后用带 Ridge 正则的线性回归做拟合。模型形式大致是 y β₀ β₁T β₂T² Σᵢ(γᵢxᵢ δᵢTxᵢ) ε其中 xᵢ 表示拆解后的催化剂变量。这样写出来的每个系数都有明确物理含义评委看论文时能顺着你的解释走而不是面对一个黑匣子。3. 乙醇转化率模型从单因素温度曲线到响应曲面3.1 先看图转化率随温度为什么是单调饱和曲线建模前先画散点图横轴温度纵轴乙醇转化率按催化剂组合着色。几乎每个组合都呈现出同一个规律低温区转化率低且上升缓慢中温区快速拉升高温区进入平台甚至略有回落。这个形状的物理解释很容易写进论文低温段反应速率受动力学控制温度升高带来的收益有限中温段主反应速率显著加快转化率近乎线性上升高温段乙醇接近完全转化转化率已经触及上限。单看某一个催化剂组合用二次多项式甚至 S 型函数都能拟合得不错。但全局模型要同时吸收温度和催化剂变量如果每个组合单独建一个 S 型函数S 型函数的参数平台值、拐点温度又变得难以统一描述。更实用的是直接做响应曲面把温度的平方项、催化剂变量的一次项、温度与催化剂的交互项全部放进同一个模型里让数据自己决定各项系数。3.2 用多项式回归加 Ridge 把转化率模型跑起来实际训练时我用的是 sklearn 的 Pipeline顺序是PolynomialFeatures 展开特征、StandardScaler 标准化、Ridge 回归。关键代码和参数如下。import pandas as pd from sklearn.preprocessing import StandardScaler, PolynomialFeatures from sklearn.pipeline import Pipeline from sklearn.linear_model import Ridge from sklearn.model_selection import GroupKFold, GridSearchCV # 字段说明 # T 反应温度Co_load 钴负载量HAP_load HAP负载量 # mass_ratio Co/SiO2 与 HAP 的装填比mode_mix 混合/分层 features [T, Co_load, HAP_load, mass_ratio, mode_mix] X df[features] y_conv df[ethanol_conv] pipe Pipeline([ # 先展开平方与交互项再统一量纲最后做带正则的线性回归 (poly, PolynomialFeatures(degree2, include_biasFalse)), (scaler, StandardScaler()), (ridge, Ridge()) ]) # 用 GroupKFold 并做网格搜索 cv GroupKFold(n_splits5) gs GridSearchCV(pipe, { poly__degree: [2, 3], ridge__alpha: [0.1, 1.0, 10.0] }, cvcv, scoringr2) # 注意groups 必须传进 fit否则 GroupKFold 不生效 gs.fit(X, y_conv, groupsdf[catalyst_id]) print(最佳参数:, gs.best_params_) print(交叉验证R2:, gs.best_score_.round(3))Pipeline 里的顺序有讲究。PolynomialFeatures 先做会生成 T²、T×mass_ratio 这类多项式列StandardScaler 再把这些列归一化避免 T² 动辄上十万的数值把 Ridge 的正则惩罚带偏最后 Ridge 回归负责收缩系数。alpha是正则强度alpha 越大系数越向零收缩过拟合风险越低但太大也会把有用的交互项压没。degree从 2 开始试如果 GroupKFold 的 R² 不够再试 3数据量只有一两百行时degree4 基本必过拟合直接不要碰。3.3 参数调整和失败排查按现象找原因如果发现训练集 R² 很高、交叉验证 R² 很低优先怀疑 degree 太高或 alpha 太小把 degree 降回 2、alpha 调到 10 再跑一次。如果换了新催化剂组合预测明显失真先查是不是切分方式写成了普通 KFold把cv换成 GroupKFold 后重新评估虚高的分数会立刻掉下来这才是真实水平。还有一种情况是某个催化剂组合的残差全部朝同一个方向偏比如预测值系统性偏低这通常不是模型的问题而是特征拆解时漏了一项——比如「装填方式」其实有三种状态但你只编码了混合/分层两种。回附录核对组合定义补一列哑变量就好。高温区残差特别大的话说明三次项被少数高温点带跑了建议单独对温度做样条变换而不是全局加三次项。4. C4 烯烃选择性模型峰值区间、交互项与残差分析4.1 选择性为什么不是单调的温度与装填比在抢同一个反应网络转化率是单调饱和曲线C4 烯烃选择性却完全是另一副面孔随着温度升高先是上升到某个温度区间达到峰值然后回落。这个峰的存在直接决定了寻优不能只看高温区。工程上可以这样解释低温段乙醇转化率太低乙醇脱水生成的中间产物量少偶合反应没有充分发生C4 选择性自然不高中温段主反应和偶合反应速率同时加快选择性被拉起来高温段深度脱氢、裂解等副反应加速把已经生成的 C4 烯烃继续消耗掉选择性反而下降。这里还有一个隐藏的交互作用不同装填比下选择性峰值所在的温度位置不一样。原因是 Co/SiO2 提供脱氢和偶合活性HAP 提供酸碱性与分散作用两者比例不同最优温度窗口就不同。这个现象意味着模型里必须显式构造「温度 × 装填比」的交互项否则峰值位置拟合不准。4.2 手动构造交互项让模型系数能解释而不是黑匣子碰运气很多队伍会把所有变量丢进 PolynomialFeatures 让它自动生成几十列交互项再一股脑喂给线性回归。这样做出来的系数正负号互相打架几乎没法解释。我一般会手动构造选择性模型的特征矩阵只留真正有物理意义的项。这样每一列对应一句人话评委问起来能答得上来。import pandas as pd from sklearn.preprocessing import StandardScaler from sklearn.linear_model import RidgeCV from sklearn.model_selection import GroupKFold, cross_val_score # 手动构造选择性模型的特征矩阵 X_sel pd.DataFrame({ T: df[T], T2: df[T] ** 2, # 处理峰值没有二次项就没有峰 T3: df[T] ** 3, # 处理峰两侧不对称 ratio: df[mass_ratio], ratio2: df[mass_ratio] ** 2, Co: df[Co_load], HAP: df[HAP_load], mode: df[mode_mix], T_ratio: df[T] * df[mass_ratio], # 温度×装填比交互 T2_ratio: df[T] ** 2 * df[mass_ratio], # 峰位置随装填比移动 Co_HAP: df[Co_load] * df[HAP_load], # 负载量间的配合 }) y_sel df[C4_select] pipe_sel Pipeline([ (scaler, StandardScaler()), (ridge, RidgeCV(alphas[0.1, 1.0, 10.0])) ]) cv GroupKFold(n_splits5) scores cross_val_score(pipe_sel, X_sel, y_sel, cvcv, groupsdf[catalyst_id], scoringr2) print(C4选择性CV R2:, scores.mean().round(3))这段代码的重点是特征设计。T3单独放进来是因为选择性峰两侧不对称如果只有 T²峰的形状被强制对称拟合残差会集中在峰附近。T_ratio和T2_ratio是这道题的关键交互项含义是「温度对选择性的影响随装填比连续变化」。模型拟合后这两个系数的正负号直接告诉你在温度升高时该加大还是减小 Co/SiO2 与 HAP 的装填比这个结论可以原封不动写进论文的催化剂设计建议。RidgeCV的好处是自带 alpha 选择交叉验证时它会自动挑出合适的正则强度省去手动网格搜索。4.3 残差分析按温度和装填比分色画图不要只盯着 RMSE模型跑完之后我一般不会急着看寻优结果先花十分钟做残差分析。做法是拿训练好的模型回代所有实验点把真实值减预测值然后按温度和装填比分别画散点图。import matplotlib.pyplot as plt pipe_sel.fit(X_sel, y_sel) pred pipe_sel.predict(X_sel) resid y_sel - pred # 左图残差随温度分布颜色映射装填比 # 右图残差随装填比分布颜色映射温度 fig, ax plt.subplots(1, 2, figsize(12, 4)) ax[0].scatter(df[T], resid, cdf[mass_ratio], cmapviridis) ax[0].axhline(0, colorgray, ls--) ax[1].scatter(df[mass_ratio], resid, cdf[T], cmapplasma) ax[1].axhline(0, colorgray, ls--)先看残差是否围绕 0 均匀分布、有没有漏斗形扩散。如果残差在峰温附近全部同号比如峰处的选择性预测系统性偏低说明 T³ 项或交互项还不够可以加一列T × ratio²再试。如果残差随装填比明显分层说明某个主效应漏进了误差里先核对特征矩阵有没有拼错列。这一步做到位论文里的模型诊断部分就有实质性内容可写而不是一句「模型拟合良好」带过。血泪经验是实验数据不是随机化采集的系统误差会藏在某个因素里只看整体 RMSE 永远发现不了。5. B题最容易翻车的5个坑现象、原因、解决确定好主模型之后真正的分水岭不在算法复杂度而在细节处理。以下 5 个坑按「踩的人头数」排序几乎每个都对应一类实际扣分点。避开它们二等奖的底线就守住了。5.1 把建模题做成了化学机理题现象论文开篇花两页写乙醇脱水、乙烯齐聚的反应网络画了一堆分子结构式正文的模型却只是对温度做个二次拟合中间没有任何推导关联。原因题目给的产物数据只有转化率和 C4 烯烃选择性没有中间产物浓度机理网络无法定量标定写得再详细也验证不了。解决机理只用来解释系数的符号。比如高温段 C4 选择性下降你就写「T²项系数为负这与高温裂解副反应加剧的趋势一致」把机理收敛成一句话的佐证而不是把反应网络当成建模主体。5.2 催化剂编号被直接数字化现象打开数据表看到「催化剂组合」列是 A1、A2 这样的编号顺手用 LabelEncoder 变成 1、2、3 就丢进回归模型还跑出了漂亮的 R²。原因类别型标签被当成有序连续量模型会学出「编号 2 介于编号 1 和 3 之间」这种完全无物理意义的规律。解决回到附录说明把每个编号拆成 Co 负载量、HAP 负载量、装填比、装填方式四列。如果个别组合确实拆不出完整信息就用 one-hot 编码兜底但注意 one-hot 变量只能代表已出现的组合不能用来外推新组合。5.3 随机切分导致验证分数虚高现象KFold 交叉验证 R² 高达 0.98但拿到新催化剂组合预测时明显失真寻优结果也落不到合理区间。原因同一组合的多个温度点被切到训练集和验证集两侧模型在验证时已经见过同一个催化剂组合的相邻温度点属于数据泄漏。解决统一改用 GroupKFold按催化剂组合整体切分。写论文时把这个切分方式明确写出来不要只贴一个漂亮的 R² 数字评审看到 GroupKFold 会比看到 0.98 更信任你的结果。5.4 目标函数只优化 C4 烯烃选择性现象寻优结果给出一个温度很高的点C4 选择性确实最高但转化率掉到三成收率反而不如中温区。原因题目关心的是 C4 烯烃收率收率 乙醇转化率 × C4 烯烃选择性只优化选择性会把搜索往高温带偏高温区虽然选择性高但转化率太低没有实际意义。解决把转化率模型和选择性模型分别建好后目标函数统一写成yield conv_model.predict(X) * sel_model.predict(X) / 100在这个复合目标上寻优。画收率等值线图时峰的位置才会落在合理区间。5.5 支撑材料程序跑不通现象评审拿到支撑材料运行Python 直接报错——文件路径不存在、缺库、随机种子没固定复现不了论文表里的数字。原因程序在你自己电脑上能跑是因为目录结构、已安装包和评审环境不一致随机抽样过程每次运行结果不同论文数字和程序输出对不上。解决程序里不写绝对路径把用到的数据列直接以 CSV 内嵌进代码目录或用相对路径读取所有随机过程固定random_state42程序跑完后打印论文里引用的 RMSE、R² 和最优工艺条件让评审一眼能对得上号。这道题支撑材料占分不低代码可复现性比模型技巧更容易拿分。6. 把二等奖方案往前再推一步插值校验与贝叶斯寻优主模型稳了之后想从二等奖往一等奖够一够可以在验证和寻优上再加两道保险。第一道是用实验点插值给模型做边界校验。选定温度和装填比作为两个主变量用 scipy 的 griddata 在实验散点上做二维插值得到一组不依赖回归假设的曲面再和你的模型预测曲面叠对比。如果两个曲面的峰值位置明显不一致说明局部数据太稀模型在那个区域的预测不可信寻优结果应该降级为「建议实验验证」而不是板上钉钉的结论如果两个曲面峰值一致就可以放心写进结论。from scipy.interpolate import griddata import numpy as np points df[[T, mass_ratio]].values values df[C4_yield].values grid_T, grid_ratio np.meshgrid( np.linspace(df[T].min(), df[T].max(), 50), np.linspace(df[mass_ratio].min(), df[mass_ratio].max(), 50)) interp griddata(points, values, (grid_T, grid_ratio), methodcubic)第二道是用贝叶斯优化替代网格搜索。二维网格搜索步长取小了计算量大取大了容易跳过峰。用 Optuna 定义目标函数返回负的 C4 烯烃收率搜索温度和装填比的范围严格限定在实验数据覆盖区间内。超出覆盖范围的预测只是外推只能给趋势性建议不能写成确定结论。def objective(trial): T trial.suggest_float(T, 300, 450) ratio trial.suggest_float(ratio, 0.1, 1.0) # 其余特征取训练集中位数predict 时拼成完整特征向量 return -(conv_model.predict(xx) * sel_model.predict(xx) / 100)第三道是从模型系数反推催化剂设计方向。比如T_ratio交互项系数为正说明温度升高时增大装填比更有利于 C4 选择性为负则相反。把这一条写进论文结尾的「催化剂设计启示」比单纯列一个最优温区更有工程说服力。这也是从「算出结果」到「讲清规律」的关键一步。我当年做这道题时第一版用的就是普通 KFold交叉验证分数漂亮得不真实换成按催化剂组合切分后瞬间掉了几个点。但也正是那个「难看」的分数逼着我重新检查特征拆解和交互项构造最后拿到的模型才对得起评审的验证。建模题最怕的不是模型不够高级而是验证方式自欺欺人。希望帮到你。本文还有配套的精品资源点击获取