2024年国赛高教杯数学建模C题农作物的种植策略解题全过程文档及程序

2024年国赛高教杯数学建模

C题农作物的种植策略

原题再现

根据乡村的实际情况，充分利用有限的耕地资源，因地制宜，发展有机种植产业，对乡村经济的可持续发展具有重要的现实意义。选择适宜的农作物，优化种植策略，有利于方便田间管理，提高生产效益，减少各种不确定因素可能造成的种植风险。
某乡村地处华北山区，常年温度偏低，大多数耕地每年只能种植一季农作物。该乡村现有露天耕地1201亩，分散为 34 个大小不同的地块，包括平旱地、梯田、山坡地和水浇地4种类型。平旱地、梯田和山坡地适宜每年种植一季粮食类作物；水浇地适宜每年种植一季水稻或两季蔬菜。该乡村另有16个普通大棚和4个智慧大棚，每个大棚耕地面积为0.6亩。普通大棚适宜每年种植一季蔬菜和一季食用菌，智慧大棚适宜每年种植两季蔬菜。同一地块（含大棚）每季可以合种不同的作物。详见附件 1。
根据农作物的生长规律，每种作物在同一地块（含大棚）都不能连续重茬种植，否则会减产；因含有豆类作物根菌的土壤有利于其他作物生长，从 2023 年开始要求每个地块（含大棚）的所有土地三年内至少种植一次豆类作物。同时，种植方案应考虑到方便耕种作业和田间管理，譬如：每种作物每季的种植地不能太分散，每种作物在单个地块（含大棚）种植的面积不宜太小，等等。2023 年的农作物种植和相关统计数据见附件 2。

请建立数学模型，研究下列问题：
问题 1 假定各种农作物未来的预期销售量、种植成本、亩产量和销售价格相对于 2023 年保持稳定，每季种植的农作物在当季销售。如果某种作物每季的总产量超过相应的预期销售量，超过部分不能正常销售。请针对以下两种情况，分别给出该乡村 2024~2030 年农作物的最优种植方案，将结果分别填入 result1_1.xlsx 和 result1_2.xlsx 中（模板文件见附件 3）。
(1) 超过部分滞销，造成浪费；
(2) 超过部分按 2023 年销售价格的 50%降价出售。
问题 2 根据经验，小麦和玉米未来的预期销售量有增长的趋势，平均年增长率介于5%~10% 之间，其他农作物未来每年的预期销售量相对于 2023 年大约有±5%的变化。农作物的亩产量往往会受气候等因素的影响，每年会有±10%的变化。因受市场条件影响，农作物的种植成本平均每年增长 5%左右。粮食类作物的销售价格基本稳定；蔬菜类作物的销售价格有增长的趋势，平均每年增长5% 左右。食用菌的销售价格稳中有降，大约每年可下降1%~5%，特别是羊肚菌的销售价格每年下降幅度为5%。
请综合考虑各种农作物的预期销售量、亩产量、种植成本和销售价格的不确定性以及潜在的种植风险，给出该乡村 2024~2030 年农作物的最优种植方案，将结果填入 result2.xlsx 中（模板文件见附件 3）。
问题 3 在现实生活中，各种农作物之间可能存在一定的可替代性和互补性，预期销售量与销售价格、种植成本之间也存在一定的相关性。请在问题 2 的基础上综合考虑相关因素，给出该乡村 2024~2030 年农作物的最优种植策略，通过模拟数据进行求解，并与问题 2 的结果作比较分析。
附件 1 乡村现有耕地和农作物的基本情况
附件 2 2023 年乡村农作物种植和相关统计数据
附件 3 须提交结果的模板文件（result1_1.xlsx，result1_2.xlsx，result2.xlsx）

整体求解过程概述(摘要)

本文研究最大化利用土地资源，建立栽种策略优化模型，利用贪心算法、随机扰动、蒙特卡洛、灵敏度检验等方法求解科学土地管理、超额出售、多因素时间波动、农作物替代性、互补性以及相关性等问题。
针对问题一，定义第 t 年的第 i 季度时，在第 j 块地种植第 k 种作物种植面积为决策变量，构建了以种植经济效益最大化为目标函数，可耕种地面积、实际可售量、连作方式、地块及作物栽种等限制为约束的种植策略线性规划模型。为科学管理土地，满足种植地不宜太分散的目标，传入参数 p、 q，分别约束每块土地最多种植作物数量，每种作物最多可种植块数量。综合考虑 p、q 尽可能小与收益尽可能大。对题目给出的数据进行预处理，统一数据格式便于读取，并统计数据生成成本与产量的三维数据表。最终，使用求解器求解：超出部分滞销结果为40244799.20元，超出部分折价结果为56325297.78元；使用贪心策略求解：超出部分滞销结果为36843378.8元，超出部分折价结果为46724871.84元。这两种求解方法各有优缺点，求解器求解结果更优，但求解慢，而贪心算法则相反，根据具体需求选择方法，两种销售情况会导致结果产生较大差异，原因归结于收益大的作物在降价后仍可保持高收益，会被高频率大面积种植。
针对问题二，考虑作物亩产量、预计销量和销售价格的波动因素，为增强风险应对能力，以最大化种植经济效益期望为目标函数，设定决策变量为波动场景下的可行策略对应的种植经济效益，增加其余参数的时间维度，在延续上一问的约束条件基础上，构建了随机规划模型。对于超过部分滞销情况，使用蒙特卡洛算法生成100组符合正态分布的随机参数序列，使用求解器在随机序列下求出100组种植策略。将分布函数离散化，可得到每组随机序列对应的概率，接着对每种规划策略进行扰动。最终得到100组随机扰动策略下，种植经济效益均值最高的种植规划策略，其中，抗波动性最强方案的种植经济效益均值为：51667432.02。
针对问题三，基于问题二模型，目标函数与决策变量保持不变，分析作物相关性及销量-价格-成本相关性对变量的影响。对相关性强且作物类型相同的作物进行替代，比较目标函数对其替换程序的灵敏度检测其替代性，最终选择用小麦替代谷子，青椒替代辣椒，对互补性强的植物进行软约束，使其尽可能协同耕种，提升效率，如豆类轮作可提升总体产量。接着，根据销量-价格-成本关系，根据销量推算合理的成本与价格。综合考虑上述因素后，在模拟数据下求解最优种植方案的种植经济效益均值为：53023389.86，相较于第二问结果更优，符合优化的目标。

模型假设：

1. 假设当季种植的农作物在当季销售，无库存。
2. 假设问题一中每种农作物的未来预期销量、种植成本、亩产量和售价相较于2023年保持稳定。
3. 假设问题二相关销量、售价等变量波动符合正态分布。
4. 假设问题三中各种农作物预期销售量与销售价格、种植成本之间存在一定相关性。

问题分析：

问题一的分析
问题一中，假定后续每一年的策略规划中，其预期销售量、种植成本、亩产量都与2023年相同。基于此分别考虑超过部分滞销和按照50%折价出售的种植策略情况，可以构建线性规划模型，对连种约束和豆类种植等约束进行限制后，使用求解器进行求解或构建贪心策略，依次遍历每块土地，选取对当前土地性价比最高的作物优先进行种植。同时考虑豆类种植约束，要求任何连续三年内豆类作物种植面积大于当前土地面积，则可以保证三年内每亩地都种过一次豆类作物。
问题二的分析
问题二中，在问题一构建的优化模型基础上，增加了作物预期销售量、种植成本、销售价格的变化条件。为输出最优的种植方案，需要在各种不确定因素和种植风险中选择一个能够在不同的预期销售量、种植成本、销售价格的变化条件中相对都有较好表现的种植策略。求解过程可以考虑使用蒙特卡洛算法生成多个随机序列，模拟不同的现实情况。
问题三的分析
问题三中，基于问题二构建的增加变量扰动的优化模型，进一步考虑农作物之间的可替代性和互补性，并综合考虑预期销售量与销售价格、种植成本之间的相关性，从而使构建的种植策略优化模型更加接近现实情况。首先，可依据农作物对总收入的灵敏度进行分析替代，从而减少农作物种类，优化种植结构。同时，对预期销售量与销售价格、种植成本进行相关性约束，从而构建增加了变量扰动和作物关联因素的种植优化模型。

模型的建立与求解整体论文缩略图

在这里插入图片描述

全部论文及程序请见下方“ 只会建模 QQ名片” 点击QQ名片即可

部分程序代码：(代码和文档not free)

import pandas as pd
import gurobipy as gp
from gurobipy import GRB
import openpyxl

# 读取Excel文件中的地块面积数据
file1 = '附件1.xlsx'
data1 = pd.read_excel(file1, sheet_name='乡村的现有耕地')
data2 = pd.read_excel(file1, sheet_name='乡村种植的农作物')

file2 = '附件2.xlsx'
data3 = pd.read_excel(file2, sheet_name='2023年的农作物种植情况')
data4 = pd.read_excel(file2, sheet_name='2023年统计的相关数据')

# 已知数据（需根据实际情况初始化）
T = 7  # 年数
I = 2  # 季节数
J = 54  # 地块数
K = 41  # 作物种类数
p = 4

M = 100000
S = data1['地块面积/亩'].tolist()
Lk = data2['Ik'].tolist()
Price = [[3.25, 7.5, 8.25, 7, 6.75,
          3.5, 3, 6.75, 6, 7.5, 40, 1.5,
          3.25, 5.5, 3.5, 7, 8, 6.75, 6.5,
          3.75, 6.25, 5.5, 5.75, 5.25, 5.5,
          6.5, 5, 5.75, 7, 5.25, 7.25, 4.5,
          4.5, 4, 0, 0, 0, 0, 0, 0],
         [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
          0, 0, 0, 0, 9.6, 8.1, 7.8, 4.5, 7.5,
          6.6, 6.9, 6.8, 6.6, 7.8, 6, 6.9,
          8.4, 6.3, 8.7, 5.4, 5.4, 4.8, 2.5,
          2.5, 3.25, 57.5, 19, 16, 100]]
Request = [[15700, 21850, 22400, 33040, 6975,
            170840, 132750, 71400, 30000, 12500,
            1500, 35100, 36000, 14000, 10000, 21000,
            36480, 26880, 6480, 30000, 35400, 43200,
            0, 1800, 3600, 4050, 4500, 34400, 9000, 1500,
            1200, 3600, 1800, 0, 0, 0, 0, 0, 0],
           [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
            0, 0, 0, 0, 0, 810, 2160, 900, 810, 0, 0,
            0, 0, 0, 0, 810, 2160, 900, 810, 0, 0]]

df1 = pd.read_excel('cost.xlsx', sheet_name='第一季')
df2 = pd.read_excel('cost.xlsx', sheet_name='第二季')
Cost1 = df1.values.transpose()
Cost2 = df2.values.transpose()
Cost = [Cost1, Cost2]

df3 = pd.read_excel('Produce.xlsx', sheet_name='第一季')
df4 = pd.read_excel('Produce.xlsx', sheet_name='第二季')
Produce1 = df3.values.transpose()
Produce2 = df4.values.transpose()
Produce = [Produce1, Produce2]

model = gp.Model("Crop_Planting")

# 决策变量
X = model.addVars(T, I, J, K, vtype=GRB.CONTINUOUS, name="X")
Y = model.addVars(T, I, J, K, vtype=GRB.BINARY, name="Y")
Z = model.addVars(T, I, K, vtype=GRB.CONTINUOUS, name="Z")
Z_rice = model.addVars(T, range(27, 35), vtype=GRB.BINARY, name="Z_Rice")

# 定义目标函数
model.setObjective(
    gp.quicksum(Price[i][k] * Z[t, 1, k] - gp.quicksum(Cost[i][j][k] * X[t, 1, j, k] for j in range(J))
               for t in range(T) for i in range(I) for k in range(K)),
    GRB.MAXIMIZE
)

# 约束1：销量不超过作物总产量
model.addConstrs((Z[t, 1, k] <= gp.quicksum(Produce[i][j][k] * X[t, i, j, k] for j in range(J))
                 for t in range(T) for i in range(I) for k in range(K)), name="Production_Limit")

# 约束2：销量不超过市场需求
model.addConstrs((Z[t, 1, k] <= Request[i][k] for t in range(T) for i in range(I) for k in range(K)), name="Demand_Limit")

# 约束3：是否种植该作物
model.addConstrs((X[t, 1, j, k] <= M * Y[t, 1, j, k]
                 for t in range(T) for i in range(I) for j in range(J) for k in range(K)),
                name="X_UpperBound_Y")

model.addConstrs((X[t, 1, j, k] >= 0.01 * Y[t, 1, j, k]
                 for t in range(T) for i in range(I) for j in range(J) for k in range(K)),
                name="X_LowerBound_Y")

# 约束4：每块地每季度种植面积总和不能超过地块总面积
for t in range(T):
    for i in range(I):
        for j in range(J):
            model.addConstr(gp.quicksum(X[t, i, j, k] for k in range(K)) <= S[j], name=f"Area_{t}-{i}-{j}")

# 约束5：三年内必须至少种植一次豆类作物
model.addConstrs((gp.quicksum(X[t, i, j, k] * Lk[k] for t in range(2) for i in range(I) for k in range(K)) >= S[j]
                 for j in range(J)),
                name="Legume_First_Two_Years")

for j in range(J):
    for t in range(T - 2):  # 以3年为单位进行检查
        model.addConstr(gp.quicksum(X[tt, 1, j, k] * Lk[k] for tt in range(t, t + 3) for i in range(I) for k in range(K)) >= S[j], name=f"Legume_{j}_{t}")

# 约束6：同一种作物在同一片土地上不能连续两个季度种植
model.addConstrs((X[t, 1, j, k] * X[t, i + 1, j, k] <= S[j]
                 for t in range(T) for j in range(J) for k in range(K) for i in range(I - 1)),
                name="No_Consecutive_Planting")
model.addConstrs((X[t, i + 1, j, k] * X[t + 1, i, j, k] <= S[j]
                 for t in range(T - 1) for j in range(J) for k in range(K) for i in range(I - 1)),
                name="No_Consecutive_Planting")

# 约束7：最多种植p种作物
model.addConstrs((gp.quicksum(Y[t, i, j, k] for k in range(K)) <= p
                 for t in range(T) for i in range(I) for j in range(J)),
                name="Max_Three_Crops")

# 添加约束：每种作物最多种在q块地上
model.addConstrs((gp.quicksum(Y[t, i, j, k] for j in range(J)) <= q
                 for t in range(T) for i in range(I) for k in range(K)),
                name="Max_Five_Plots_Per_Crop")

# 约束 8：确保粮食作物在连续年份的第一季不能连种
model.addConstrs((X[t, 0, j, k] + X[t + 1, 0, j, k] <= S[j]
                 for t in range(T - 1) for j in range(J) for k in range(1, 16)),
                name="No_Consecutive_Years_For_Grain")

# 约束 9：编号为 1-26 的土地在第二季不种植任何作物
model.addConstrs((X[t, 1, j, k] == 0
                 for t in range(T) for j in range(26) for k in range(K)),
                name="No_Planting_Second_Season_For_Lands_1_26")
# 约束：编号为 1-26 的土地上只能种植编号为 1-15 的作物
model.addConstrs((X[t, i, j, k] == 0
                 for t in range(T) for i in range(I) for j in range(26) for k in range(15, 41)),
                name="No_Planting_Crops_16_41_On_Lands_1_26")

# 约束：编号为 1-15 的作物只能种植在编号为 1-26 的土地上
model.addConstrs((X[t, i, j, k] == 0
                 for t in range(T) for i in range(I) for j in range(26, J) for k in range(15)),
                name="No_Planting_Crops_1_15_On_Lands_27_54")

# 约束：编号为 27-34 的土地种植水稻
model.addConstrs((gp.quicksum(X[t, i, j, k] for i in range(I) for k in range(K) if k == 15) <= M * Z_rice[t, j]
                 for t in range(T) for j in range(27, 35)),
                name="Rice_Planting_Only_Once")

# 确保水稻只能种植在单季
model.addConstrs((gp.quicksum(X[t, i, j, 15] for i in range(I)) <= S[j]
                 for t in range(T) for j in range(27, 35)),
                name="Single_Season_Rice")

# 添加约束 1：如果种植了水稻，则第二季不种植任何作物
model.addConstrs((gp.quicksum(X[t, i, j, k] for k in range(K)) <= M * (1 - Z_rice[t, j])
                 for t in range(T) for j in range(27, 35)),
                name="No_Second_Season_If_Rice")

# 添加约束 2：第一季只能种植 17-34 号作物
model.addConstrs((gp.quicksum(X[t, 0, j, k] for k in range(16, 35)) == gp.quicksum(X[t, 0, j, k] for k in range(16, 35))
                 for t in range(T) for j in range(27, 35)),
                name="First_Season_Crops_17_34")

# 添加约束 3：第二季只能种植 35-37 号作物
model.addConstrs((gp.quicksum(X[t, 1, j, k] for k in range(34, 38)) == gp.quicksum(X[t, 1, j, k] for k in range(34, 38))
                 for t in range(T) for j in range(27, 35)),
                name="Second_Season_Crops_35_37")

# 添加约束：编号为 35-37 的作物只能种植在编号为 27-34 的土地上
model.addConstrs((X[t, i, j, k] == 0
                 for t in range(T) for i in range(I) for j in range(26) for k in range(34, 37 + 1)),
                name="No_Planting_Crops_35_37_On_Lands_1_26")

# 添加约束 4：编号为 38-41 的作物只能在 35-50 号地的第二季种植
model.addConstrs((X[t, i, j, k] == 0
                 for t in range(T) for j in range(35) for k in range(37, 41)),
                name="No_Planting_Crops_38_41_On_Lands_1_34")

# 添加约束 5：编号为 38-41 的作物只能种植在第二季
model.addConstrs((X[t, 0, j, k] == 0
                 for t in range(T) for j in range(35, 51) for k in range(37, 41)),
                name="No_Planting_Crops_38_41_First_Season")
# 设置相对Gap
model.setParam('MIPGap', 0.01)

# 优化模型
model.optimize()

# 输出结果
if model.status == GRB.OPTIMAL:
    print(f"Optimal solution found with objective value: {model.objVal}")
    for i in range(I):
        for j in range(J):
            for k in range(K):
                if X[t, i, j, k].x > 0:
                    print(f"Year {t + 1}, Season {i + 1}, Land {j + 1}, Crop {k + 1}: {X[t, i, j, k].x} acres planted")
    print(f"Optimal solution found with objective value: {model.objVal}(元)")