★ FEATURED ARTICLE 简介2024年数学建模国赛C题完整代码与数据包面向参赛学生、数据分析学习者及建模爱好者围绕“作物种植优化”场景提供从数据清洗到模型求解的整套Python实现方案。压缩包共51个文件约2.4MB其中17个py脚本构成核心代码覆盖异常值处理、相关性分析、聚类分析、遗传算法、线性规划及结果可视化等完整流程13张png图片展示了商品单价横向柱状图、相关性热图、QQ图、箱线图、散点图矩阵等诊断与输出结果9个xlsx和3个csv文件用于存储中间数据、优化结果及附件表格另有README、.idea配置等辅助文件目录结构清晰。内容按问题一、问题二、问题三分模块组织每个问题都配有独立脚本与对应输出清晰还原国赛C题从数据探索到模型建优的求解思路尤其适合准备数模竞赛的团队对照实操、快速迁移方法。目前已有5961人浏览学习可作为备赛冲刺和算法实战的优质参考。1. 2024数学建模国赛C题代码和数据先理清这场比赛到底考什么2024数学建模国赛C题自带一个天然矛盾题目名字叫“农作物的种植策略”听起来是农业题但真正比拼的是数据整理和优化建模能力。赛后的搜索热词里“代码”和“数据”排在最前说明多数队伍的实际困难不是数学公式而是怎么把一堆Excel表格变成能喂给求解器的输入。这道题的考点可以概括成三句话从冗余数据里识别出决策变量、把文字要求改写成约束条件、在不确定条件下给出可解释的最优方案。本文按这个链路完整拆一遍你不需要有农业背景只要会Python和数据框架操作就能跟着复现。2. 拆题与数据梳理地块、作物、价格、成本建模前先画出四张表2.1 题目附件的常见结构四类数据表如何对位2024年国赛C题的数据量在近几年国赛里不算最大但胜在表多、字段杂。通常你拿到的附件包括地块信息表、作物信息表、未来几年的预测价格表以及题目正文里以文字形式给出的种植规则。地块表里一般有地块编号、地块类型、地块面积作物表里有作物编号、作物类别、单产、种植成本、补贴等字段价格表则按年份给出一组预测价格区间。我的习惯是先把所有表读进内存然后手动画一张数据地图哪张表提供哪些字段、哪些字段之间有外键关系、哪些字段是建模用不到的冗余列。这一步看起来琐碎但它决定了后面的代码是清爽还是灾难。import pandas as pd # 读取题目附件的三个核心表请替换成实际文件名 df_land pd.read_excel(land_info.xlsx) # 地块信息 df_crop pd.read_excel(crop_info.xlsx) # 作物信息 df_price pd.read_excel(price_info.xlsx) # 预测价格 # 快速体检查看字段名和前两行数据 print(地块表字段, df_land.columns.tolist()) print(df_land.head()) print(作物表字段, df_crop.columns.tolist()) print(df_crop.head()) print(价格表字段, df_price.columns.tolist()) print(df_price.head())逻辑说明先跑一遍head()而不是直接建模是为了尽早发现列名不匹配、数据单位异常等问题。国赛附件里的列名经常带中文和特殊符号比如“单产kg/亩”、“成本元/亩”这些括号和单位字符在后续代码里容易踩坑所以我通常会在这里顺手把列名重命名成纯英文字段。参数说明land_info.xlsx、crop_info.xlsx、price_info.xlsx是我按常见命名习惯起的占位文件名你实际操作时一定要替换成附件里的真实文件名。head()默认显示前5行如果你发现列名有不可见字符可以用df.columns [c.strip() for c in df.columns]做一次清洗。2.2 决策变量定义连续面积变量是首选别一上来就搞0-1这道题的核心决策是“种什么、种多少、哪年种”所以决策变量最自然的定义是连续面积变量[ x_{ijt} \text{第 } i \text{ 块地上第 } t \text{ 年种植作物 } j \text{ 的面积亩} ]为什么不用0-1变量面积种植本身是可分的种2.5亩小麦在现实中完全合理没必要引入整数约束。用了0-1变量后模型从线性规划LP升成混合整数规划MILP求解难度直接上一个台阶。除非题目明确要求“每种作物只能整块地种”或者“要么不种、要么种满”否则用连续变量是更稳妥的选择。另一个设计点是时间维度 (t)。如果只按单年建模轮作约束无法表达如果把多年放在同一个模型里变量数会成倍增加。我的建议是保留时间维因为C题的核心约束之一就是年际轮作删掉时间维等于把最有区分度的考点砍了。变量数量需要心里有数假设有20块地、15种作物、5年变量总数是 (20 \times 15 \times 5 1500) 个这个规模对线性规划求解器来说是轻松任务。真正让模型变复杂的是约束条件不是变量数。2.3 约束条件逐项拆解面积、轮作、占比三类底线约束从C题的题目背景看至少有三类约束是必须写进模型的。第一类是面积约束每块地所有作物的种植面积之和不能超过地块面积。第二类是轮作约束同一地块相邻两年不能种植同一种作物或者同一科作物这是农业实践里的重茬禁忌。第三类是占比约束题目通常规定豆类作物种植面积不得低于总种植面积的一定比例这与土壤肥力修复有关。轮作约束的不等式写法如下[ x_{ijt} x_{ij(t1)} \leq A_i, \quad \forall i, j, t ]这里 (A_i) 是第 (i) 块地的面积。因为 (x_{ijt}) 本身不会超过 (A_i)这个约束的实际含义是同一块地、同一种作物两年中最多只有一年能种。它没有用二元变量但通过把两年面积和限制在总面积以内有效阻止了连作方案。这是一个很典型的“用线性约束表达逻辑约束”的技巧在论文里值得写清楚。占比约束的写法类似[ \sum_{\text{豆类作物 } j} x_{ijt} \geq \theta \cdot \sum_{\text{所有作物 } j} x_{ijt}, \quad \forall i, t ]这里的 (\theta) 是题目给出的最低比例。注意这个约束两边都含变量需要整理成标准形式把左边减到右边得到非正不等式再手动乘以负号变成 (\leq) 形式才能放进线性规划求解器的 (A_{ub}) 矩阵。3. 确定性建模与求解用 scipy 线性规划跑出第一个最优方案3.1 把目标函数和约束组装成标准形式第一问通常是在不考虑产量和价格波动的前提下制定最优种植方案本质是一个线性规划问题。目标函数是总收益最大化收益由三部分组成售价收入减去种植成本如果题目给了补贴政策再把补贴加进去。用数学形式表达就是[ \max Z \sum_{i,j,t}(p_{jt} \cdot y_{ijt} - c_{ijt}) \cdot x_{ijt} ]其中 (p_{jt}) 是第 (t) 年作物 (j) 的预测售价(y_{ijt}) 是单产(c_{ijt}) 是单位面积种植成本。在scipy.optimize.linprog里所有问题默认转成最小化所以目标系数取负号。组装顺序是整个建模的关键先确定变量索引顺序再按同一顺序生成目标系数、约束矩阵、边界。我用的是“地块-年份-作物”三维打平原因是这样生成面积约束行时最直观——每行只需要在一个连续的区间内填1。3.2 核心代码从 DataFrame 到 linprog 求解的一条龙流程下面给出一个可以直接运行的框架你只需把文件路径和字段名替换成自己题目的附件配置。import numpy as np import pandas as pd from scipy.optimize import linprog # 1. 数据读取与字段预处理 df_land pd.read_excel(land_info.xlsx) df_crop pd.read_excel(crop_info.xlsx) df_price pd.read_excel(price_info.xlsx) # 统一列名建议手动改成简洁英文 df_land.columns [land_id, area, land_type] df_crop.columns [crop_id, crop_type, yield, cost] df_price.columns [crop_id, year, price] land_list df_land[land_id].tolist() crop_list df_crop[crop_id].tolist() year_list [2024, 2025, 2026, 2027, 2028] # 按题目实际年份修改 n_land len(land_list) n_crop len(crop_list) n_year len(year_list) n_vars n_land * n_year * n_crop # 2. 构建收益系数矩阵 # 变量顺序: x[i][t][j] - 展开后索引 (i * n_year t) * n_crop j revenue np.zeros((n_land, n_year, n_crop)) price_lookup df_price.set_index([crop_id, year])[price].to_dict() yield_lookup df_crop.set_index(crop_id)[yield].to_dict() cost_lookup df_crop.set_index(crop_id)[cost].to_dict() for i in range(n_land): for t, yr in enumerate(year_list): for j, crop in enumerate(crop_list): p price_lookup.get((crop, yr), 0) y yield_lookup.get(crop, 0) c cost_lookup.get(crop, 0) revenue[i, t, j] p * y - c # 3. 目标系数linprog 默认最小化取负 c -revenue.reshape(-1) # 4. 面积约束: 每个地块每年的种植面积 地块面积 A_ub [] b_ub [] for i in range(n_land): for t in range(n_year): row np.zeros(n_vars) start (i * n_year t) * n_crop row[start : start n_crop] 1 # 该地块该年份的所有作物变量求和 A_ub.append(row) b_ub.append(df_land[area].iloc[i]) # 5. 轮作约束: 同一地块同一种作物相邻两年不同时种植 for i in range(n_land): for t in range(n_year - 1): for j in range(n_crop): row np.zeros(n_vars) idx1 (i * n_year t) * n_crop j idx2 (i * n_year (t 1)) * n_crop j row[idx1] 1 row[idx2] 1 A_ub.append(row) b_ub.append(df_land[area].iloc[i]) # 6. 变量边界: 面积不能为负 bounds [(0, None)] * n_vars # 7. 求解 res linprog(c, A_ubnp.array(A_ub), b_ubnp.array(b_ub), boundsbounds, methodhighs) if res.success: x res.x.reshape(n_land, n_year, n_crop) print(最优总收益, round(-res.fun, 2)) # 输出一个地块的第一年种植结果用于人工验证 print(第 1 块地 2024 年各作物种植面积, x[0, 0]) else: print(求解失败, res.message)逻辑说明代码分七个步骤前两步是数据准备和收益系数计算中间四步是组装约束矩阵最后一步是求解。price_lookup用set_index加to_dict建了一个以作物年份为键的字典这样查价格的时间复杂度是 O(1)避免在循环里反复用布尔筛选数据量大时能明显提速。参数说明methodhighs是 scipy 1.6 以上版本默认推荐的求解器底层是 HiGHS 开源实现比旧的simplex方法数值稳定性好、速度更快。bounds全设为(0, None)表示变量下界为0、上界不限。如果你的题目有“某种作物种植面积不能超过某个上限”的要求在这里改成(0, max_area)即可。3.3 结果校验求解器说最优不代表方案能用求解器返回successTrue只是第一步我还需要做三项人工校验否则论文里可能出现低级错误。第一把解还原成表格后随机抽几个“地块-年份”组合手算面积总和确认没有超过地块面积。第二检查是否有某个作物的种植面积为0——如果某些作物在所有地块上都是0说明它们在当前价格和产量下没有竞争力这个结果在逻辑上没问题但要在论文里解释为什么淘汰了这些作物不能假装它们不存在。第三把总收益跟题干的参考数值对比如果数量级不对大概率是单位换算出错。这三步建议写成脚本自动跑而不是用眼睛看# 重建解表并按地块年份汇总面积 result_df pd.DataFrame([ { land: land_list[i], year: year_list[t], crop: crop_list[j], area: x[i, t, j], } for i in range(n_land) for t in range(n_year) for j in range(n_crop) ]) summary result_df.groupby([land, year])[area].sum().reset_index() merged summary.merge(df_land[[land_id, area]].rename(columns{land_id: land}), onland) merged[violation] merged[area] merged[area_y] 1e-6 # 面积校验 invalid merged[merged[violation]] if len(invalid) 0: print(存在面积超限的记录) print(invalid) else: print(所有地块-年份组合均满足面积约束)这段代码先按“地块年份”汇总种植面积再和地块总面积做对比1e-6是浮点误差容忍度。如果你的方案里有面积超限记录不要慌先检查约束矩阵生成部分——大概率是索引错位导致某些变量没有被约束覆盖。4. 不确定性场景产量和价格波动下怎么保证方案不是纸上谈兵4.1 确定性最优解在波动面前有多脆弱第一问的确定性方案在理想参数下表现完美但第二问一旦引入产量和价格的不确定性问题性质就变了。设想两个场景第一种某作物价格下跌20%第二种某作物遭遇产量歉收30%。你会发现第一问的最优方案在这两个场景下的实际收益远低于预期甚至可能亏本。原因在于确定性模型把参数当作已知常数求解器会毫不犹豫地把资源集中在“看起来最赚钱”的作物上。而真实世界里的种植决策需要的是稳健性——即使某个作物价格下跌整体方案也不会崩盘。国赛第二问考察的正是你如何把这种不确定性的考量融入优化模型。常见做法有以下三种蒙特卡洛模拟加期望值优化、随机规划两阶段或多阶段、鲁棒优化。国赛题目通常不要求你用太高级的理论但你要能把方法讲清楚、结果解释得合理。4.2 蒙特卡洛模拟生成场景集分布假设与随机种子控制蒙特卡洛模拟的核心是给单产和价格设定合理的随机分布然后生成大量场景。单产我一般用三角分布下界和上界取基准值的±20%30%众数取基准值本身。价格也用三角分布但如果题目给了预测区间直接用区间均匀分布也可以。三角分布优于均匀分布的原因它把“最可能值是基准值”的先验信息融入采样不容易产生极端离谱的场景。下面这段代码生成500个场景供后续的期望值计算使用。import numpy as np # 固定随机种子保证论文结果可复现 rng np.random.default_rng(42) n_scen 500 base_yield df_crop[yield].values # 基准单产长度 n_crop base_price df_price[df_price[year] 2024].set_index(crop_id).reindex(crop_list)[price].values # 波动区间上下浮动 20% low_y base_yield * 0.8 high_y base_yield * 1.2 low_p base_price * 0.8 high_p base_price * 1.2 scen_yield np.zeros((n_scen, n_crop)) scen_price np.zeros((n_scen, n_crop)) for s in range(n_scen): scen_yield[s] rng.triangular(low_y, base_yield, high_y) scen_price[s] rng.triangular(low_p, base_price, high_p) # 打印前3个场景的第一种作物检查采样合理性 print(前3个场景的作物单产前5个品种) print(scen_yield[:3, :5]) print(前3个场景的作物价格前5个品种) print(scen_price[:3, :5])逻辑说明rng.triangular(low, mode, high)是 numpy 内置的三角分布采样函数三个参数分别是最小值、众数、最大值返回数组形状和输入一致。这里low_y、base_yield、high_y本身是数组所以一次调用就能生成一整行场景数据不需要在作物维度循环。参数说明n_scen500是场景数量。500个场景足以让期望值收敛到小数点后两位再多会拖慢后面的优化计算。default_rng(42)是新版 numpy 推荐的随机数接口不要用旧版np.random.seed()前者是局部生成器不会污染全局状态。固定种子是论文复现的基础务必保留。4.3 从上到下用场景平均系数替代确定性系数跑一版稳健解有了场景集之后最常见的第一步是计算每个作物的期望单产和期望价格然后用期望值替代确定性模型里的原始参数。这个方法在学术上叫 Sample Average ApproximationSAA思路很直接既然无法预知哪个场景会发生就让方案对所有场景的平均表现最优。实际上SAA的数学形式比直观理解要多走一步。如果第二阶段的决策不能根据场景调整那么SAA退化为直接用期望系数求解的确定性等价模型很容易实现# 计算期望单产和期望价格按场景取平均 avg_yield scen_yield.mean(axis0) # 长度 n_crop avg_price scen_price.mean(axis0) # 长度 n_crop # 重新组装收益系数其余约束完全复用第3章的代码 revenue_stoc np.zeros((n_land, n_year, n_crop)) for i in range(n_land): for t in range(n_year): for j in range(n_crop): revenue_stoc[i, t, j] avg_yield[j] * avg_price[j] - cost_lookup.get(crop_list[j], 0) # 求解复用A_ub, b_ub, bounds c_stoc -revenue_stoc.reshape(-1) res_stoc linprog(c_stoc, A_ubA_ub, b_ubb_ub, boundsbounds, methodhighs) if res_stoc.success: x_stoc res_stoc.x.reshape(n_land, n_year, n_crop) # 把稳健方案放回每个具体场景里计算实际收益分布 actual_rev np.zeros(n_scen) for s in range(n_scen): scen_coef np.zeros((n_land, n_year, n_crop)) for i in range(n_land): for t in range(n_year): for j in range(n_crop): scen_coef[i, t, j] scen_yield[s, j] * scen_price[s, j] - cost_lookup.get(crop_list[j], 0) actual_rev[s] np.sum(x_stoc * scen_coef) print(500个场景下的收益均值, round(actual_rev.mean(), 2)) print(5%分位数最差运气下的收益, round(np.percentile(actual_rev, 5), 2))逻辑说明这段代码先求出期望系数并求解然后把解固定下来回放到500个具体场景中计算实际收益分布。actual_rev的分布特征是第二问论的关键数据如果收益的5%分位数不是负数说明方案在最差情况下仍然保本这个结论比单纯报一个期望值更有说服力。参数说明percentile(actual_rev, 5)计算的是5%分位数代表最差的5%场景下的收益水平。论文里建议同时报告均值、标准差、5%分位数三个指标分别对应期望表现、波动程度、风险底线。如果你的方案在5%分位数为负说明稳健性不够需要考虑引入更保守的约束比如强制要求某些保收作物必须种植一定面积。5. 2024 国赛 C 题避坑清单数据处理和建模求解中 5 个常见的翻车现场5.1 数据单位不统一收益算出来差一个量级现象总收益算出来比题目参考值大了十倍或小了十倍检查代码逻辑却没有发现明显的错误。原因面积单位“亩”和“公顷”混用价格单位“元/斤”“元/公斤”“元/吨”没有统一换算。国赛附件里不同表格的单位习惯不一致这是最隐蔽的坑。解决在数据处理开头就做统一换算面积一律换算成“亩”价格一律换算成“元/公斤”单产一律换算成“公斤/亩”。换算完以后打印一份关键字段的汇总统计人工检查量级合理性——每亩收益通常在几百元到几千元之间如果某个作物的每亩收益超过一万或者为负数说明换算或字段对应有问题。5.2 漏掉轮作约束方案被评委一句话问住现象第一问求出的最优方案里同一种作物在同一块地上连续多年大面积种植答辩时被质疑不符合农业常识。原因题目文字里对轮作的约束分散在各段描述里建模时只顾着面积约束把轮作要求漏掉了。解决把轮作约束的合法性质检写成独立函数在求解前和求解后各跑一次。更重要的在生成A_ub矩阵后随机抽几行打印非零位置用肉眼确认轮作约束确实在矩阵里且索引正确。这个“抽行检查”的习惯比我预想的更能避免低级错误。5.3 随机数没固定种子论文结果复现不出来现象同样一份代码两次运行得到不同的“最优收益”被队友质疑结果可靠性。原因蒙特卡洛模拟环节没有固定随机种子或者混用了np.random.seed()和default_rng()两种接口导致抽样结果不一致。解决全程只使用np.random.default_rng(固定整数)并在论文的参数设置表里明确写出场景数量和随机种子。这是一个成本极低但边际收益极高的习惯——评审看到“固定种子场景数”的写法会默认你具备基本的可复现意识。5.4 约束矩阵行数爆炸求解时间从秒级变成分钟级现象把多年数据和占比约束全部加入后linprog的求解时间从几秒涨到几分钟甚至提示内存不足。原因约束矩阵行数接近数万甚至数十万其中大量行是线性相关的冗余约束。虽然 HiGHS 求解器有预求解机制但过大的矩阵依然拖慢迭代。解决先做约束合并。同一个地块的多个年份面积约束如果形式完全一致可以合并成更紧凑的区块矩阵。另外可以按年份把模型拆成多个独立子问题分别求解再汇总C题的多年规划模型如果没有年份间约束除轮作外分年求解完全可行速度会快一个量级。5.5 求解器返回 infeasible却找不到是哪条约束出错现象res.status等于 2不可行约束矩阵检查了好几遍都看不出问题。原因约束方向写反或边界设置错误。典型的例子是把“种植面积不超过地块面积”写成了“不小于地块面积”或者把b_ub传成了负数。解决用一个“全变量赋0.5”的探测向量做快速验证。把所有变量都设成0.5代入每一条约束看哪条约束被违反。违反的那条就是问题所在检查它的系数和常数项即可。这个方法比肉眼扫矩阵快得多十分钟能解决的问题没必要耗一小时。6. 进阶验证敏感性分析、求解器横向对比与方案呈现的三个技巧敏感性分析是让方案从“能用”变成“有说服力”的分水岭。做法是把产量和价格分别按±10%、±20%扰动重新求解记录最优收益的变化幅度。如果产量波动10%导致收益波动超过30%说明方案对产量极度敏感论文里就应该建议决策者重点关注抗风险品种如果敏感性很低方案就是稳健的可以大胆声称“在合理波动范围内仍然保持最优性”。# 敏感性分析封装函数复用第3章的求解逻辑 def solve_with_factor(yield_f, price_f): rev np.zeros((n_land, n_year, n_crop)) for i in range(n_land): for t in range(n_year): for j in range(n_crop): y base_yield[j] * yield_f p base_price[j] * price_f c cost_lookup.get(crop_list[j], 0) rev[i, t, j] p * y - c res linprog(-rev.reshape(-1), A_ubA_ub, b_ubb_ub, boundsbounds, methodhighs) return -res.fun if res.success else np.nan for fy in [0.8, 0.9, 1.0, 1.1, 1.2]: val solve_with_factor(fy, 1.0) print(f产量因子 {fy:.1f}: 最优收益 {val:.2f})求解器选择上scipy.optimize.linprog的highs方法能覆盖国赛绝大部分线性规划场景。如果你想把求解换成开源求解器建议用PuLP加 HiGHS 后端或者直接用scipy.optimize.milp处理混合整数扩展。商业求解器在国赛这种规模下不是必需除非你遇到了特别复杂的整数规划模型。最后一个技巧是结果呈现。论文里不要只贴“总收益最大”的最终结果至少给三张表最优种植方案表、各年份收益明细表、敏感性分析结果表。方案表按地块分行、按年份分列每格填作物名称和面积。收益明细表按年份汇总总收益和主要作物贡献占比。敏感性分析表用行列交叉的热力图为佳颜色深浅一眼看出哪些参数是命脉。我个人的习惯是每一版求解结果都导出到results/目录文件名带上时间戳和关键参数说明。这个习惯救过我一次——有一次改了价格数据后发现最优收益反而更高靠对比旧结果文件才发现是成本字段被误改了一个单位。数据和代码分开存放结果文件只读不覆盖是一个朴素但可靠的工程习惯。希望这个习惯也能帮到你。本文还有配套的精品资源点击获取
阅读完成 · 觉得有帮助？
