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 = 4M = 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)) <= pfor 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)) <= qfor 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] == 0for 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] == 0for 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] == 0for 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] == 0for 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] == 0for 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] == 0for 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}(元)")