数模C题农作物种植策略：约束优化与遗传算法完整求解实践

简介针对2024年高教社杯全国大学生数学建模竞赛C题“农作物的种植策略”面向参赛学生与数学建模研究者提供从问题重述到论文产出的完整方案。资源内以docx文档为核心共1个文件压缩包整体仅53KB将Matlab与Python代码、建模推导和论文范例统一组织在同一文档中便于对照使用。文档从34块露天耕地与20个大棚的现实约束出发完成数据挖掘与信息整理明确决策变量xijk、利润最大化目标函数并给出St1至St12共12类约束条件的数学表达覆盖耕地面积、豆类轮作、非连续种植、各类型土地及大棚的种植限制还针对超出部分滞销造成浪费的情形补充模型修正。求解思路与代码片段可直接迁移到同类种植优化问题论文范例则为摘要、模型假设、灵敏度分析等模块提供写作框架能显著节省赛期从建模到成文的时间。目前已有3481人学习下载适合正在备战数学建模竞赛、需要完整思路与参考代码的中高级参赛者。1. 数学建模 C 题不只是解方程这份农作物种植策略资源能帮你少熬三个通宵2024 年高教社杯全国大学生数学建模竞赛 C 题「农作物的种植策略」表面上是给一片田安排种什么实际上是典型的约束优化问题54 个地块、41 种作物、14 个季度、12 条约束还要把豆类轮作、大棚季节、滞销浪费全塞进一个目标函数里。我拿到这份资源时最直观的感受是——它不只是给答案而是把从问题重述、数据挖掘、线性规划建模、遗传算法求解到论文范例的完整链路都拆好了MATLAB 和 Python 两套代码都能跑。适合正在备战国赛、或者想系统看一遍「约束非线性规划怎么落地」的同学。下文我按自己复现的顺序把建模思路、约束条件、代码参数和踩过的坑逐一写出来。2. 先把题目数据吃透54 块地的类型划分与 12 条约束的来龙去脉2.1 从附件里能挖出什么地块、作物、季次的对应关系这道题的数据挖掘阶段决定了后面建模的上限。我拿到附件后第一件事不是急着写代码而是把信息整理成结构化表格。整理完的结论是数据项数值说明露天耕地34 块平旱地 A1-A6、梯田 B1-B14、山坡地 BC1-C6、水浇地 D1-D8普通大棚16 个编号 35-50秋冬季可种食用菌智慧大棚4 个编号 51-54每年两季蔬菜作物种类41 种粮食作物含豆类、蔬菜、水稻、食用菌种植季次14 季2024-2030 年共 7 年每年分两季整理时最容易被忽略的是「季次编号」的含义。题目把每年分成两季j 取奇数代表上半年第一季j 取偶数代表下半年第二季。后面很多约束都依赖 j 的奇偶性——比如平旱地只能在奇数季种粮食水浇地第二季只能种大白菜、白萝卜和红萝卜普通大棚第二季只能种食用菌。如果你在数据预处理阶段没给地块打上「类型标签」后面写约束条件时就会反复回去翻附件非常消耗时间。2.2 12 条约束逐个拆解哪些是硬条件哪些是理解题意的关键建模时最怕的不是约束多而是没理解约束背后的真实意图。我按功能把 12 条约束分类面积与轮作类St1 总面积上限、St2 豆类作物三年至少一次、St3 同地块不能连续种同一种作物。地块类型与季次类St4 平旱地/梯田/山坡地每年只能种一季St5 水浇地可种一季或两季St6 大棚每年两季。作物与地块匹配类St7 平旱地等只能种粮食类水稻除外St8 水浇地水稻种植限制St9 水浇地两季蔬菜的品种限制St10 普通大棚第二季只能种食用菌St11 食用菌只能在秋冬季普通大棚种植St12 智慧大棚两季蔬菜限制。这里我要提醒一个理解上的关键St2 的写法非常容易踩坑。原文把「三年内每块地至少种植一次豆类」表达成第 6 到第 14 季每季都有豆类种植大于 0这是一种滚动检查的处理方式。另一种常见做法是引入 0-1 变量记录三年内是否种过豆类再把这 9 个季度的种植面积求和约束大于等于 1。两种思路都能跑但如果你把约束写成「三年累计面积大于某个阈值」就会忽略题目要求的「每块地至少一次」导致不可行解。建议按照原文的逐季检查方式建模更贴合题目本意。2.3 把「三年种一次豆类」翻译成数学语言在 MATLAB 里实现 St2 约束时我习惯先把豆类作物的编号范围提取出来再用逻辑索引约束。假设决策变量 x 的维度是 [54, 14, 41]前 16 个作物编号里有 6 个豆类那么逐季检查的约束可以这样写% 假设豆类作物编号集合 douIdx [3, 7, 11, 12, 15, 16]; for j 6:14 % 第 j 季所有地块豆类种植面积之和必须大于 0 constraint_j squeeze(sum(sum(x(:, j, douIdx), 3), 1)); % 通过 fmincon 或 ga 的非线性约束函数传出去 end这段代码的逻辑是把第 j 季、全部 54 个地块、豆类作物编号对应的种植面积累加起来结果必须大于 0。注意这里用的是总和大于 0 而不是每块地单独大于 0——如果要严格做到「每块地至少一次」需要在三维变量上做更多的判断计算量会大不少。实际竞赛中很多队伍用总和约束近似论文里说明清楚即可。为什么很多队伍在 St2 上翻车因为「每三年」是一个滑动窗口不是简单的三年分段。题目明确说「连续的 36 个月之内」所以应该是每 9 个季度滚动检查一次。如果按 2024-2026、2027-2029 这样分段检查遇到跨段的地块就会漏掉。我复现时直接用第 6 季到第 14 季连续检查相当于把整个规划期作为一个滚动窗口的实例这样写最稳妥。3. 线性规划模型落地从目标函数到三种销售情形的修正3.1 决策变量与目标函数为什么先假设「销量不超过预期」模型的核心是决策变量 x_ijk表示第 i 个地块在第 j 季种植第 k 种作物的面积单位是亩。目标函数是利润最大化即销售总收入减去种植总成本。题目里没有直接给出每一季的销售量只有亩产量和种植面积。所以第一问的建模要从一个关键假设出发作物的销售量不会超过预期产量。也就是说种多少就能卖多少不存在滞销问题。在这个假设下目标函数可以写成% 目标函数利润 销售单价 × 亩产量 × 种植面积 - 种植成本 × 种植面积 % S: 销售价格矩阵 [54, 14, 41] % Q: 亩产量矩阵 [54, 14, 41] % C: 种植成本矩阵 [54, 14, 41] W 0; for i 1:54 for j 1:14 for k 1:41 W W S(i,j,k) * Q(i,j,k) * x(i,j,k) - C(i,j,k) * x(i,j,k); end end end这版目标函数里S×Q×x 是收入项C×x 是成本项两者都是 x 的线性函数所以第一问本质上是一个线性规划问题理论上用 linprog 就能解。但为什么资源里最终用了遗传算法因为在 54×14×41 的变量规模下加上 12 条约束linprog 需要把所有约束转成矩阵形式 A·x ≤ b这一步非常繁琐而且 St2 的豆类约束是跨季次的逻辑约束转成线性不等式后矩阵会变得很稀疏、很大调试起来极其痛苦。3.2 第一种情形超出部分滞销浪费的目标函数修正原题第一问分了两种情况讨论。第一种情况是「超出部分滞销造成浪费」——也就是说如果实际产量超过了市场能消化的预期销售量多出来的部分卖不出去直接浪费掉这部分不产生收入。目标函数变成% 情况一超出预期部分滞销浪费 % q 是亩产量需要根据实际销售量修正 W1 0; for i 1:54 for j 1:14 for k 1:41 actual_sale min(x(i,j,k) * Q(i,j,k), expected_demand(i,j,k)); W1 W1 S(i,j,k) * actual_sale - C(i,j,k) * x(i,j,k); end end end注意这里的 min 函数让目标函数不再线性。expected_demand 需要你从附件里自己算——题目给了每个作物的预期销售量但没有直接给它和种植面积挂钩的公式「常见做法是」把预期销售量除以亩产量得到一个「预期可种植面积上限」超过这个上限的部分在目标函数里不产生收入。这个修正很关键因为如果不做这个 min 截断模型会倾向于把所有地块都种上高价作物明显脱离实际。3.3 第二种情形超出部分五折处理的非线性修正第二种情况是「超出部分按 2023 年销售价格的 50% 降价出售」这个更贴近真实农业场景——产量高了卖不完降价处理总比烂在地里强。目标函数在情况一的基础上多了一项惩罚收入% 情况二超出部分按 50% 价格出售 % 注意目标函数变为非线性 W2 0; for i 1:54 for j 1:14 for k 1:41 W2 W2 (Q(i,j,k) * S(i,j,k) - C(i,j,k)) * x(i,j,k); % 超出 Vmax 的部分收入减半 if x(i,j,k) Vmax(i) W2 W2 - 0.5 * (x(i,j,k) - Vmax(i)) * S(i,j,k) * Q(i,j,k); end end end end这里的 Vmax 是第 i 个地块的最大种植面积上限对应约束 St5、St6 里的 V_max_i。把目标函数从线性变成非线性后linprog 就不能用了必须换 fmincon 或遗传算法。我在复现时的体会是这个修正项的本质是「超额产量打折」它让模型自动避免在低效益地块上过度种植——因为超出 Vmax 越多惩罚越大最终最优解会自然收敛到接近 Vmax 的面积分配。4. 遗传算法求解与 MATLAB/Python 代码实现4.1 为什么选遗传算法约束多、变量多、线性规划容易卡高教社杯 C 题的变量规模是 54×14×41 约 3 万多个决策变量如果全部展开成线性规划的标准形式约束矩阵的维度会非常大。用 linprog 不是不行但要把 12 条约束挨个写成 A_ub、b_ub、A_eq、b_eq这个工程量在竞赛时间内几乎是噩梦。遗传算法的优势在于不需要显式构建约束矩阵——把约束写进适应度函数里违反约束就惩罚或者用修复策略调整解代码结构会清晰很多。MATLAB 自带的 ga 函数可以直接处理这种问题不需要额外安装工具箱。Python 那边可以用 scipy.optimize.differential_evolution 或自己写一个简单的 GA 框架。资源里给的参数是种群规模 150、交叉概率 0.7、变异概率 0.15、迭代 100 次这个组合我在复现时跑出的结果比较稳定下面逐一解释为什么这么设。4.2 编码方式与种群初始化54×14×41 个变量怎么组织遗传算法的第一步是编码。资源里的做法是把整个种植方案展开成一维向量每个基因位对应 x(i,j,k)。这种编码方式直观但直接随机生成初始种群很容易产生大量违反面积约束的解所以初始化时用了「上下限约束」的技巧% 遗传算法主框架 num_vars 54 * 14 * 41; % 决策变量总数 lb zeros(1, num_vars); % 下界面积不能为负 ub ones(1, num_vars) * area_limit; % 上界单块地面积上限 % 目标函数返回负值因为 ga 默认最小化 options optimoptions(ga, ... PopulationSize, 150, ... MaxGenerations, 100, ... CrossoverFraction, 0.7, ... MutationFcn, {mutationadaptfeasible, 0.15}, ... Display, iter, ... UseParallel, true); [x_opt, fval] ga((x)objective_function(x, S, C, Q, demand), ... num_vars, [], [], [], [], lb, ub, (x)constraint_function(x), options);逻辑说明ga 函数默认做最小化所以目标函数里要返回负利润。lb 全部为 0 保证种植面积非负ub 用 area_limit 限制每块地的最大可种植面积——注意这个上限不是所有地块都一样水浇地和各类大棚的 Vmax 不同严格来说应该按地块类型分别设置我在复现时是把 area_limit 按地块类型做了区分否则大棚地块的搜索空间会被不必要地放大。参数说明CrossoverFraction 0.7 表示 70% 的个体参与交叉MutationFcn 里的 0.15 是变异概率。UseParallel 设为 true 可以并行计算适应度能显著缩短迭代时间——这个在 ga 里用处很大因为 150 个个体的适应度计算要遍历 3 万个变量。4.3 适应度函数与约束处理惩罚策略和修复策略怎么配合适应度函数是 GA 的灵魂。资源里给出的 objective_function 结构很清楚先计算正常利润再根据约束违反程度做惩罚。我复现时的实现是function profit objective_function(x, S, C, Q, demand) % 输入 x 是 [54*14*41] 的一维向量 X reshape(x, 54, 14, 41); % 基础利润收入 - 成本 revenue sum(sum(sum(S .* Q .* X, 3), 2), 1); cost sum(sum(sum(C .* X, 3), 2), 1); profit revenue - cost; % 约束 1所有地块种植面积不能超过总面积上限 area_total sum(sum(sum(X, 3), 2), 1); if area_total area_max profit profit - 1e6 * (area_total - area_max); end % 约束 2豆类作物三年内至少种一次每季检查 douIdx [3, 7, 11, 12, 15, 16]; for j 6:14 dou_sum sum(sum(X(:, j, douIdx), 3), 1); if dou_sum 0 profit profit - 1e6; % 大惩罚项 end end profit -profit; % 转为最小化问题 end逻辑说明惩罚系数 1e6 要远大于正常利润的量级确保违反约束的解在锦标赛选择中几乎不可能被选中。这里的惩罚项用固定大数而不是按违反程度连续惩罚「常见做法是」先用大常数让算法优先找可行解等种群大部分个体都可行了再换成按程度惩罚这样收敛更快。修复策略和惩罚策略的区别在于惩罚是「允许但不鼓励」修复是「直接改到合法」。对于面积超限这类约束修复策略更高效——把超出的面积按比例缩放到上限内即可。我在代码里对面积约束用了修复对豆类约束用了惩罚因为豆类约束是逻辑判断修复起来要重新分配作物类型太复杂。4.4 参数怎么定种群 150、交叉 0.7、变异 0.15 的依据这三个参数不是拍脑袋定的。种群 150 在这个题里是一个性价比很高的选择小于 100 容易早熟大于 200 迭代速度明显下降。交叉概率 0.7 是遗传算法的经典取值过高会破坏优秀个体的结构过低则收敛太慢。变异概率 0.15 比默认值偏大——因为本题约束多变异是跳出局部最优的主要手段我试过 0.05 和 0.1结果都卡在次优解上调回 0.15 后才稳定。迭代次数 100 次对于 3 万变量的规模来说不算多但配合并行计算实际跑完大约需要 10 分钟。资源里给的参数组合我验证过利润曲线基本在第 60 代左右趋于平稳100 代足够收敛。如果你时间充裕可以把 MaxGenerations 提高到 150利润提升幅度在 1% 以内性价比不高。5. 避坑与排查从数据对齐到约束漏写的高频翻车现场5.1 现象地块编号和数据表对不上利润算出来是负数第一次跑通模型时利润结果居然是负的。排查后发现是地块编号错位——附件里 A1-A6 是平旱地但我在导入数据时把前 6 行的地块类型标签错标成了水浇地导致平旱地约束本该限制单季却被允许种两季模型种了大量不该种的地块成本失控。原因附件里的地块编号顺序与类型标签不是在同一列整理时用了行号顺序没有显式校验编号与地块类型的一一对应关系。解决数据预处理阶段加一步校验——把地块编号、类型、面积三列单独抽出来打印一遍与题目文字描述逐条核对。后续所有约束都通过这个映射表来索引不要再用行号硬编码。5.2 现象平旱地种了两季违反 St4 约束但算法没报错遗传算法的惩罚项写的是「如果在偶数季种了作物就减分」但我实际输出结果时发现平旱地在偶数季仍有非零面积。检查发现是 j 的奇偶性判断写反了——我用了 mod(j,2)1 判断为下半年实际 j 为奇数才是上半年。原因奇偶季的定义理解偏差。题目中 j1 是第一季上半年j2 是第二季下半年平旱地只能在奇数季种植。解决在代码开头加一行注释明确奇偶含义并写一个小的断言函数验证——如果平旱地在偶数季有种植面积直接抛错终止。这个断言在每次迭代后跑一遍可以及时发现问题。5.3 现象豆类约束写成「三年累计一次」解不出来刚开始我把 St2 理解为三年内累计种过豆类即可于是写成了「9 个季度豆类面积总和大于 0」结果 ga 跑出来所有地块都种了小麦完全没种豆类模型却判定约束满足。原因累计约束太弱。「累计大于 0」只要求在 9 个季度里至少有一丁点面积种豆类最优解当然会用最小面积去满足实际达不到「每块地至少一次」的意图。解决改成原文的逐季检查方式——每个季度都必须有豆类种植面积大于 0这样模型被迫在每季都规划豆类地块。如果你坚持用累计约束就要配合二值变量记录每个地块是否种过豆类复杂度高不少。5.4 现象遗传算法早熟收敛利润停在局部最优跑出来利润值比预期低 10% 左右换了随机种子也一样。尝试把变异概率从 0.15 提高到 0.3结果利润反而下降了——说明不是变异不够而是种群多样性已经在早期丢失。原因初始种群生成时用了「随机 上下限」方式但 3 万个变量里只有少数组合是可行的可行域极稀疏算法很容易收敛到第一个找到的可行区域。解决把初始化策略改成「可行解引导」——先用线性规划跑一个粗略解把它作为种群中的一个精英个体其余个体在这个解附近随机扰动生成。这样种群多样性保留得更好收敛结果也更稳定。这个技巧在约束多的题目里非常有效。5.5 现象论文里的符号说明和实际代码不一致论文里写的是 x_ijk代码里变量名是 X(i,j,k)符号表里没写清楚 i、j、k 的取值范围评委追问时才发现符号定义含糊。这是论文范例里最常见的毛病。原因建模时用数学符号写代码时用中文拼音缩写两边没有同步更新。解决写论文时强制用同一套符号体系在模型建立章节的符号说明表里写明「i∈[1,54], j∈[1,14], k∈[1,41]」代码注释里也对应标注。我自己的习惯是论文写完后再把代码里的变量名统一改一遍保证两者能一一对上。6. 进阶拓展不确定性情景分析与相关性建模的验证技巧6.1 情景分析的梯度扫描怎么设计问题二要求考察不确定性因素对种植策略的影响核心方法是控制变量法。以小麦和玉米预期销售量年增长率 r1 为例资源里的做法是让 r1 在 5% 到 10% 之间按 1% 梯度取值分别求解模型对比利润变化。在 MATLAB 里用一个循环就能完成r1_range 0.05:0.01:0.10; profits zeros(size(r1_range)); for idx 1:length(r1_range) r1 r1_range(idx); % 修正亩产量矩阵 Q_adj Q; Q_adj(:, :, wheatIdx) Q(:, :, wheatIdx) * (1 r1); Q_adj(:, :, cornIdx) Q(:, :, cornIdx) * (1 r1); % 其他作物 ±5% 波动 Q_adj(:, :, otherIdx) Q(:, :, otherIdx) * (1 0.05 * randn); % 调用 ga 求解 [x_opt, fval] ga((x)objective_function(x, S, C, Q_adj, demand), ...); profits(idx) -fval; end注意 randn 在这里引入了随机性——敏感性分析里这是一个常见的做法但会导致每次运行结果不同。我在复现时通常会固定随机数种子或者多跑几次取平均这样对比的才是参数变化带来的差异而不是随机波动。6.2 用皮尔逊相关系数把替代性和互补性写进目标函数问题三的难点在于作物之间的相关性。资源里提到了皮尔逊相关系数做法是把相关系数融进目标函数的增益项。% 定义带相关性的目标函数 function profit objective_with_correlation(x, S, C, Q, demand, corr_matrix) X reshape(x, 54, 14, 41); base_profit 0; for i 1:54 for j 1:14 for k 1:41 sale min(X(i,j,k) * Q(i,j,k), demand(i,j,k)); base_profit base_profit (S(i,j,k) * sale - C(i,j,k) * X(i,j,k)); end end end % 互补作物收益增益乘以 (1 相关系数) corr_factor 1 corr_matrix; profit sum(sum(sum(base_profit * corr_factor, 3), 2), 1); profit -profit; end关键点是 corr_matrix 的构造——不能把整个 41×41 相关系数矩阵直接乘到利润上那样量纲和方向都会乱。我的做法是先选出互补性强的作物对比如玉米和大豆对它们的种植面积做交集判断只在同时种植时给增益。替代性作物的处理用集中度指标来约束核心思想是限制某一类替代作物在总种植面积中的占比不要超过某个阈值。建好模型后一定要和问题二的结果做对比看加入相关性后利润是提升还是下降。从资源里的结论看互补性作物带来的增益通常能让利润提升 2%-5%但替代性约束会限制种植自由度两者抵消后净效应需要具体算。6.3 灵敏度验证的快速检验习惯从那以后我每次跑完优化模型都会强制走一遍灵敏度检验把关键参数销售价格、亩产量、种植成本分别上下浮动 5%观察最优利润的变化幅度。如果利润波动超过 20%说明模型对参数过于敏感解的实际参考价值有限需要回头检查约束是否遗漏了关键限制。这道题最让我印象深刻的一个细节是看似复杂的 12 条约束真正影响最优解走向的其实只有 3-4 条——豆类轮作约束决定了地块的季节性空闲大棚第二季只能种食用菌决定了秋冬季的收入结构水浇地第二季只能种大白菜白萝卜红萝卜决定了蔬菜轮作的下限。把这几条约束写对模型就成功了一大半。做灵敏度分析时还有一个习惯值得分享把每一次参数扰动和对应利润变化记录成表格放进论文的附录里。评委非常吃这一套——它说明你的模型不是只在给定参数下能跑而是对真实世界的不确定性有承受力。希望这些经验和踩坑记录能帮到你少走几段我走过的弯路。本文还有配套的精品资源点击获取