2024数学建模国赛C题：Python+Gurobi实现整数规划与轮作约束

简介面向2024年国赛C题整理而成的代码与数据包内容覆盖作物种植优化问题的完整分析流程适合备战全国大学生数学建模竞赛的学生、指导教师以及想了解数据分析和优化建模的初学者。压缩包共五十一个文件以十七个Python脚本为主体配合十三张可视化图表、九个Excel结果表、三个CSV文件另有Markdown说明与项目配置整体约2.4MB脚本涵盖数据清洗、异常值检验、相关性分析、聚类分析、线性规划与遗传算法求解等关键环节结构清晰可按题目顺序逐一复现。自发布以来已有五千九百六十人浏览学习。通过这份资源可以获得一套可运行的赛题解法既能查看从原始附件到结果表、诊断图的完整分析脉络也能对照图表和表格理解每一步建模依据是实战演练和赛前复习中很实用的参考资料。1. 2024数学建模国赛C题代码和数据的正确打开方式2024数学建模国赛C题常被叫成“一眼看懂、一写就跪”几十个地块、几十种作物、七年的种植安排表面是数据题实际上是一套带地块限制、轮作约束和阶梯销售价格的整数规划模型。官方发下来的是数据和题目原文你要交的是每块地每年种什么、怎么收益最大。下面按“代码和数据”这条线拆哪些列是约束、哪些列是系数模型怎么建代码怎么写才能跑通坑在哪。适合手上有附件但卡在建模阶段的本科生也适合想快速复现一份可改方案的研究生。这道题的难点在把题目原文转成约束地块面积种满、同一地块不连种同一作物、超出订单量部分降价卖。三条规则看着简单写成模型时各有各的坑。下面从读数据开始一句一句立模型。2. 先读懂C题的数据和隐藏约束地块、作物、订单和轮作怎么写进模型2.1 把附件1到附件4读成四张约束表再谈建模常见做法是官方给四个Excel附件地块信息、上年种植情况、销售量与价格、作物类型。它们不是用来查数的是四类模型输入。我一般会先把每张表当什么用写在一张纸上再打开Excel附件关键字段在模型里当什么附件1地块信息地块编号、地块面积、地块类型/土质面积类约束和地块适宜性约束附件2上年种植情况地块、作物名称、面积初始状态决定2024年不能种什么附件3销售量与价格作物名称、亩产、销售单价、种植成本、销售量目标函数系数和分段销售拐点附件4作物类型作物名称、类型是否豆类豆类轮作规则的判断依据注意附件3里的“销售量”决定分段点在哪个位置这个字段的作用域后面要单独确认它是全区订单量还是每亩订单量直接影响超产折扣算多少。附件2则特别容易被人忽略很多人把它当成历史数据展示实际上它是第一年轮作约束的右手边。读表顺序我建议固定成“附件1→附件2→附件3→附件4”。附件1先把地块个数和面积讲清楚后面所有约束都以地块为最小单位附件2给你2023年每块地种了什么这是2024年计划的地基附件3提供目标函数的系数附件4只是补充豆类判断放最后读不影响理解。2.2 轮作约束定义0-1变量把“不能连种”写成不等式轮作约束的正确写法依赖一组0-1变量。设 y_{l,c,t}1 表示地块 l 在第 t 年种植作物 c那么每年每块地只能种一种作物就是Σ_c y_{l,c,t} 1“同一地块不能连续种植同一种作物”写成y_{l,c,t} y_{l,c,t1} ≤ 1同一个地块和同一个作物组合上任意相邻两年最多出现一次。注意是“地块作物”组合不是“全区域某作物面积不能超过去年”。我见过有队伍把轮作约束写成“今年玉米总面积不能大于去年玉米总面积”这是总量约束不是轮作评阅一眼就能看出来。豆类规则看起来特殊实际不冲突“种过豆类的地块次年可种任何作物”指的是豆类在轮作周期里可以插入任意地块不需要额外排斥豆类作物自身仍然遵守间隔约束。附件4的作用就是标记哪些作物属于豆类让代码可以对豆类单独打标。如果题目允许“豆类十年内只能出现一次”这种更强限制再用额外约束处理C题一般不要求。2.3 决策变量是连续还是整数面积用实数是否种植用0-1种植亩数 x_{l,c,t} 我建议用非负连续变量是否种植用0-1变量 y_{l,c,t}。原因是土地可以任意切分不需要把亩数做成整数段连续变量配合0-1变量已经能把“种什么、种多少”表达清楚。若把面积也做成整数变量数量会扩大几个量级求解时间从分钟级变成小时级而且结果并不会更实用。数量级不用怕34个地块、7年、20种作物0-1变量大约 34×7×204760 个连续变量同一量级Gurobi这类求解器几秒到几十秒就能出可行解。如果题目要求上报整数亩可以在连续解基础上四舍五入再人工调整而不是一上来就让求解器算整数面积。3. 用PythonGurobi跑通最小求解从Excel到可提交方案的完整代码3.1 先读数据再打印列名不要急着建模import pandas as pd # 文件都放在 data/ 目录下文件名以官网下载为准 land pd.read_excel(data/附件1.xlsx, header0) plant pd.read_excel(data/附件2.xlsx, header0) sell pd.read_excel(data/附件3.xlsx, header0) crop_type pd.read_excel(data/附件4.xlsx, header0) for df, name in zip([land, plant, sell, crop_type], [land, plant, sell, crop_type]): df.columns [str(c).strip().replace(\n, ) for c in df.columns] print(name, df.columns.tolist()) print(df.head())这一步的目的是把表头统一成干净字符串。Excel表头经常带前后空格、换行或不可见字符直接引用列名会报KeyError而报错位置又往往不在读取这行会让你误以为是模型写错。header0表示第一行当列名strip去掉首尾空格replace(\n,)把换行清掉。打印前几行是让自己确认每个字段的单位和含义后面建模全靠这份打印结果。3.2 把Excel字段变成Python字典面积、价格、成本和订单量L list(land[地块编号]) C list(crop_type[作物名称]) Y range(2024, 2031) # 规划年份2024到2030 area dict(zip(land[地块编号], land[地块面积])) price dict(zip(sell[作物名称], sell[销售单价])) # 元/吨 cost dict(zip(sell[作物名称], sell[种植成本])) # 元/亩 order dict(zip(sell[作物名称], sell[销售量])) # 全区域订单量 yield_ dict(zip(sell[作物名称], sell[亩产])) # 吨/亩 bean dict(zip(crop_type[作物名称], crop_type[作物类型].astype(str).str.contains(豆类)))用字典而不是反复切片DataFrame是为了建模时按作物名O(1)取数。注意C要以附件4的作物名称为准然后用断言检查它和附件3的作物集合是否一致不一致说明两个附件里的叫法有差异先统一再往下走。yield_的单位要和order单位匹配如果附件里亩产是公斤、订单量是吨这里就要先做除以1000的换算。单位问题会在第4章重点说但建模前先核对能避开大半。3.3 创建变量、面积约束和轮作约束三行代码扛住核心难点import gurobipy as gp from gurobipy import GRB m gp.Model(C2024_crop) x m.addVars(L, C, Y, vtypeGRB.CONTINUOUS, namex) # 种植亩数 y m.addVars(L, C, Y, vtypeGRB.BINARY, namey) # 是否种植 for l in L: for t in Y: # 每年每块地面积和必须等于地块面积 m.addConstr(gp.quicksum(x[l, c, t] for c in C) area[l], namefarea_{l}_{t}) # 每年每块地只能种一种作物 m.addConstr(gp.quicksum(y[l, c, t] for c in C) 1, namefone_crop_{l}_{t}) for c in C: # 面积超过0时y必须为1x上限是地块面积 m.addConstr(x[l, c, t] area[l] * y[l, c, t], nameflink_{l}_{c}_{t}) # 轮作同一地块同一作物不能出现在相邻两年 for l in L: for c in C: for t in range(2024, 2030): m.addConstr(y[l, c, t] y[l, c, t 1] 1, namefrotation_{l}_{c}_{t}) # 第一年轮作要看2023年历史种植 hist dict(zip(plant[地块编号], plant[作物名称])) for l in L: if l in hist and hist[l] in C: m.addConstr(y[l, hist[l], 2024] 0, namefhist_{l})link约束是把连续变量x和0-1变量y绑定只要某种作物在这块地的种植面积超过0对应的y就必须等于1。由于one_crop约束每年只能选一种作物link约束不会把多块地的面积搞混。area约束用等式是因为题目要求每块地每年都必须种满如果题目允许休耕且休耕有收益/成本这里再改成不等式加休耕变量。轮作约束的range(2024,2030)覆盖2024到2029t1正好压到2030不会漏最后一年。历史种植约束是很多队伍第一次跑出infeasible的根源。如果2023年种的是某种菜类不在附件4的作物列表里hist判断会把它跳过这时要回到题目确认这种作物是否允许原地换种再决定要不要人工补一条“2024年不种该菜类”的约束。3.4 目标函数里的分段销售线性化订单量以内原价超出部分五折u m.addVars(C, Y, nameu) # 按订单价售出的产量 v m.addVars(C, Y, namev) # 超产折价售出的产量 obj gp.LinExpr() for c in C: for t in Y: total gp.quicksum(x[l, c, t] * yield_[c] for l in L) # 全区总产量 m.addConstr(u[c, t] v[c, t] total, namefsplit_{c}_{t}) m.addConstr(u[c, t] order[c], nameforder_cap_{c}_{t}) # 订单量以内按原价超出部分按50% obj price[c] * u[c, t] 0.5 * price[c] * v[c, t] # 种植成本按亩数扣除与最终销量无关 obj - gp.quicksum(cost[c] * x[l, c, t] for l in L) m.setObjective(obj, GRB.MAXIMIZE) m.optimize() if m.status GRB.OPTIMAL: out pd.DataFrame([ {地块: l, 年份: t, 作物: c, 面积: round(x[l, c, t].X, 2)} for l in L for t in Y for c in C if x[l, c, t].X 1e-6 ]) out.to_excel(C2024_plan.xlsx, indexFalse) print(最优收益, round(m.objVal, 2))这里没有用min/max直接写非线性函数而是拆成u和v两个变量uv等于总产量u不超过订单量。因为目标函数是最大化v只会取到真正超出订单量且必须折价的那部分产量不会虚高。0.5是当年C题的超产折扣系数其他年份和赛题要看原文不要照抄。order[c]的作用域在前面提过这里再强调一次它是全区域订单量不是每亩订单量否则超产部分会被严重高估或低估。求解完后导出方案表至少要包含地块、年份、作物、面积四列后续检查轮作和订单量都靠这张表。round两位小数是为了防止浮点误差把面积写成0.0000001这种值筛选 x1e-6 同理。提示模型一直不可行时先写一句 m.write(model.lp) 把LP文件导出用文本编辑器搜约束名比盲改代码快得多。4. C题避坑指南五个让代码跑不通和结果不对的常见问题4.1 模型不可行轮作约束撞上历史种植数据现象只加面积约束时模型秒出结果一加轮作约束就报infeasible且没有任何提示是哪条约束出错。原因最常见是2023年历史种植没有接进模型导致2024年计划里出现和2023年完全相同的作物另一个常见原因是地块适宜作物约束漏了某些地块只允许种一两种作物而你在全集上要求它七年内不能连种几种约束叠加后无解。解决先把历史种植约束按3.3的方式加进去再把附件1里“适宜作物/地块类型”列读出来对不允许的作物直接限制x0或者只在允许组合上建变量。跑之前打印 len(land)、len(C) 核对数据规模。4.2 收益数字怪订单量作用域和单位是两大黑匣子现象结果里某种作物总产量上千吨但收益只有几万元或者总收益夸张到几十亿明显不符合常识。原因如果order是每亩订单量而你把它当成全区订单量超产部分全被折价处理收益被低估如果亩产单位是公斤、订单量单位是吨目标函数里产量数量级差了1000倍数字自然离谱。我见过最惨的情况是价格单位是元/斤、亩产单位是公斤/亩换算系数混在一起模型还“正常”算完最后收益比真实值大了近10倍。解决建模前打印量级自查。比如玉米订单量如果是几千吨说明它是全区订单量如果是几百公斤说明量纲不对或作用域理解错了。单位统一的原则是所有进目标函数的量都换成“吨”和“元”。亩产公斤就除以1000价格元/斤就乘以2换成元/公斤再乘以1000换成元/吨。宁可多写一行换算也不要写在注释里骗自己。4.3 方案一眼假地块限制没落全休耕变量不要乱加现象方案里同一块地几年内反复种同一种作物或者水稻出现在明显不应该种水稻的地块上还有些队为了让收益更高让部分地块“自愿撂荒”结果总成本下降、收益上升看起来很聪明实际上已经违背题目。原因附件1的地块类型/土质限制没有转成约束或者面积约束用了“≤地块面积”而不是“地块面积”给模型留了休耕的缝。很多队把“≤”当成给求解器更多自由但题目明确要求每块地每年都种植这个自由度不该给。解决面积约束用等式地块限制在建变量阶段就处理。允许种植集合可以用附件1的字段生成allow {1号地: [玉米, 大豆], ...}然后只在允许组合上addVar。如果题目真的允许休耕官方会在正文里给补贴标准那时的处理是引入休耕变量并把它当成一种“种植方案”参与轮作而不是直接删面积约束。4.4 轮和“不能重茬”被写成总量限制现象代码能跑结果检查时发现某块地连续三年种玉米但全区域玉米总面积每年不同于是以为没问题。原因轮作约束是“地块作物”级别很多队伍图省事写成“年度总量约束”比如今年玉米总面积不得超过去年。这个约束不限制单块地连作评阅时只要抽查连续两年种植表就能发现违规。解决直接保留 y[l,c,t] y[l,c,t1] 1 的逐地块写法不要为了减少约束数量把轮作升级成总量限制。总量限制可以作为合理性检查指标不能替代轮作。4.5 最优解出现“同一地块一年种两种作物”的假象现象输出表某地块同一年出现两行不同作物面积加起来正好等于地块面积甚至还有交叉面积。原因one_crop约束和link约束都加了但如果x变量的下标用了“地块、年份”而不是“地块、作物、年份”或者导出时把不同作物的x累计在了一起就会出现这种假象。另外浮点误差也可能让本应为0的x残留小值。解决导出表里先按“地块年份”分组检查每个组内作物种类数是否都为1筛选阈值1e-6能滤掉数值噪声。如果确实出现两种作物面积都大于阈值回头检查one_crop约束是否使用了正确的变量维度。5. 让结果经得起追问影子价格、灵敏度分析与一页附件表模型跑出最优解只是第一步论文还差两块为什么某些作物不种、方案对价格波动是否可靠。前者用影子价格解释后者用灵敏度分析补证据。影子价格可以看订单量约束的价值。求解完对每种作物打印 order_cap 这条约束的Pi值for c in C: constr m.getConstrByName(forder_cap_{c}_{2024}) print(c, constr.Pi)Pi接近0说明订单量不是瓶颈这种作物多种少种对收益影响不大Pi明显大于0说明订单量每增加一吨总收益就能增加一截论文里可以写“若未来订单量提升优先扩种该作物”。注意整数规划下的对偶值并不是严格定义稳妥做法是先固定整数变量再重解LP把LP影子价格当边际贡献讨论这个边界在论文里最好说明。灵敏度分析我习惯跑三组场景价格×0.9、价格×1.1、订单量×0.9各重解一次记录总收益和主粮面积。如果种植结构基本不动结论是“方案主要由轮作和面积约束驱动价格扰动影响有限”如果玉米面积跟着价格大幅起伏就在论文里承认收益对价格敏感并给出应对策略。这张表放附录评阅印象分通常不错。我现在拿到这类题第一件事不是写求解器代码而是先把附件1到附件4的字段和单位打印出来逐条对着题目原文过一遍因为所有坑都藏在数据列的表述里。代码是最后二十分钟的事前面两小时都在确认约束和量纲。希望帮到你。本文还有配套的精品资源点击获取