基于蒙特卡洛模拟的农作物动态种植策略
Dynamic Crop Planting Strategy Based on Monte Carlo Simulation
摘要: 在农村发展和环境保护的时代背景之下,农村耕地资源合理利用和可持续化发展至关重要。本文主要研究农作物动态种植策略,首先以利润最大化为目标,结合某区域气候条件、耕地资源、种植成本和销售价格等因素,建立了线性规划模型,制定最优的农作物种植策略。接着,采用蒙特卡洛模拟处理不确定因素,建立农作物种植的动态规划模型,并利用经济学上的交叉弹性公式来确定农作物之间的互补性和可替代性,讨论了可替代性的种植最优方案。基于上述建立的三种模型利用给定的数据制定了某区域未来五年的最优种植策略,在一定程度上为农村农作物的种植提供了可靠的理论依据和实践指导。
Abstract: Under the backdrop of rural development and environmental protection, the rational utilization and sustainable development of rural arable land resources are of paramount importance. This paper primarily focuses on the study of dynamic crop planting strategies. First, with the goal of profit maximization, a linear programming model is established by integrating factors such as regional climate conditions, arable land resources, planting costs, and sales prices to formulate the optimal crop planting strategy. Next, Monte Carlo simulation is employed to address uncertainties, leading to the development of a dynamic programming model for crop planting. Additionally, the cross-elasticity formula from economics is utilized to determine the complementarity and substitutability between crops, and the optimal planting scheme based on substitutability is discussed. Using the three models established above and the given data, the optimal planting strategy for a specific region over the next five years is formulated, providing, to a certain extent, a reliable theoretical foundation and practical guidance for rural crop cultivation.
文章引用:常芳欣, 冯爱芬, 黄璐锴, 轩晗. 基于蒙特卡洛模拟的农作物动态种植策略[J]. 应用数学进展, 2025, 14(6): 361-374. https://doi.org/10.12677/aam.2025.146326

1. 引言

在农村发展和环境保护的时代背景之下,农村耕地资源合理利用和可持续化发展至关重要[1]。本文主要以华北山区的某乡村为例进行农作物种植的优化研究,数据参考2024年全国大学数学建模竞赛C题数据[2]。已知该地区温度常年来偏低且可以作为耕地的土地极其有限,现在有露天的耕地和非露天耕地,露天耕地分为平旱地、梯田、山坡地和水浇地四种类型共1201亩并被分为了34个不同地块,非露天的地分为普通大棚和智慧大棚,分别有16个和4个且每种大棚地有0.6亩。对于不同的地块可以种植的农作物不同,例如水稻地只能种水稻或者蔬菜,水稻不能种植在干旱地,梯田和山坡地。那么对于不同的地块我们需要种植可以适应该地块条件并且可以存活的植物。那么对于同类植物,由于需要从土地吸取的养分相同,那么在种植的时候就要避免连续种植相同的植物,防止由于植物的种植使得土地失去肥力。同时我们了解到种豆类有利于补充土地的养分,那么我们就要满足在三年内所有的土地豆类至少要种植一次。同时,考虑到方便耕种作业和田间管理,每种作物的种植地不能太分散,单个地块种植面积不宜太小。依据2023年种植数据,需要为该乡村2024~2030年制定最优种植方案。

首先,在假定各种农作物未来预期销售量、种植成本、亩产量和销售价格相对于2023年保持稳定的情况下,分别考虑超过部分滞销造成浪费和超过部分按2023年销售价格的50%降价出售两种情况下的最优种植方案。其次,考虑各种农作物预期销售量、亩产量、种植成本和销售价格的不确定性以及潜在种植风险来制定最优种植方案。然后,综合考虑农作物之间的可替代性和互补性以及预期销售量与销售价格、种植成本之间的相关性,给出最优种植策略并与前面结果作比较分析。

2. 基本假设

基本假设如下:

(1) 同一地块(含大棚)每季可以合种不同的作物。

(2) 每种作物在同一地块(含大棚)都不能连续重茬种植,否则会减产。

(3) 从2020年开始每个地块(含大棚)的所有土地三年内至少种植一次豆类作物。

(4) 每种作物每季种植地块不能太分散并且单个作物种植面积不宜太小。

(5) 气候变化具体情况以及市场需求突发变化暂不考虑。

(6) 同一地块可以合种不同的作物,互为互补的作物对资源的竞争是互不干扰的。为了简化计算,假设每种作物的亩产量在混合种植且互为互补品的情况下增长10%,反之亩产量不变。

3. 基于线性规划的最优种植方案

在农作物相关指标稳定的情况下,研究作物总产量超过预期销售量的两种处理方式对种植方案的影响:

(1) 滞销浪费情况:超出的部分直接扔掉。

(2) 降价销售情况:超出部分按2023年销售价格的50%降价出售。

3.1. 数据预处理

Shapiro-Wilk正态性检验法是一种用概率回归直线的斜率来建立统计量进行该假设检验,计算简单,功效高[3]

Shapiro-Wilk正态性检验法[4]的具体步骤如下:统计假设H0:子样来自正态总体。

(1) 按照子样数据的大小,从小到大排列成 X ( 1 ) , X ( 2 ) ,, X ( n1 ) , X ( n )

(2) 检查Shapiro-Wilk系数 a in 表,查出对应于n的各 a in 的值。

(3) 计算统计量W的值。

(4) 选定检验水平 α ,按n α W分布表,查对应的 W( n,α ) 的值。

(5) 下结论:当 WW( n,α ) 时,否定H0,即总体不是正态分布的。反之则认为总体是正态分布的。

利用Shapiro-Wilk法对亩产量、种植成本和销售单价的数据进行正态性检验,根据拟合程度判断均不符合正态分布。利用箱型图判断种植成本和销售单价存在异常值,由于数据来源于真实统计,异常值可能是由于农作物稀有或难种植导致。如羊肚菌售价较高且被市场认可,因此这些异常值不作处理。

3.2. 决策变量

在解决问题的过程中,我们将第i种类型的土地种植第j类农作物在第k季的产品种植面积即 x ijkl 作为决策变量。其中k表示季节,单季为0,第一季和第二季分别为1和2。l为年份下标,范围为0, 1, ..., 7。j是根据附件1中的乡村种植农作物的编号来赋值的,j的范围为1, 2, ..., 41。i是根据附件1中给出的乡村现有耕地,我们对其进行了编号i的范围为1, 2, ..., 54,如下表1所示:

Table 1. The plots and numbers for crop cultivation

1. 作物种植的地块及编号

编号

地块名称

1

A1

2

A2

……

……

3

A54

由于在约束条件时我们需要对各区的面积进行约束,我们可以根据附件1的乡村现有耕地对各区所占面积进行计算,各区总面积如图1所示。

3.3. 目标函数

第一种情况超出的部分不盈利 δ 取为0,第二种情况超出部分半价出售 δ 取为0.5,以最大化利润为目标函数(利润 = 总收益 − 成本),建立线性规划模型[5]

maxZ= i=1 54 j=1 41 k=0 2 l=1 7 x ijkl b ijk ( P ijk Q ijk ) +δ P ijk max{ i=1 54 j=1 41 k=0 2 l=1 7 x ijkl B ijk Y ijk ,0 } (1)

其中,Z其中,Z表示总收益, Y ijk 表示第i种类型的土地种植第j类农作物在第k季的农作物预期销售量, P ijk 表示第i种类型的土地种植第j类农作物在第k季的销售单价, Q ijk 表示第i种类型的土地种植第j类农作物在第k季的农作物成本, δ 为链接两种情况的变量,若为第一种则 δ 为0,若为第二种情况则 δ 为0.5。 B ijk 为第i种类型的土地种植第j类农作物在第k季的产量, b ijkl 为第i种类型的土地种植第j类农作物在第k季的每亩预期产量。

Figure1. Distribution area by region

1. 各地区分布面积图

3.4. 约束条件

3.4.1. 种植作物约束

不同的土地适宜生长的植物不同,也存在同一种植物无法在多个类型的土地上生存,那么就会出现在第i块地上无法种植第j类型的农作物。

(1) 不可种植在水浇地和大棚的农作物

x ijk =0,i=27,28,,54;j=1,2,,15

公式含义为这些农作物j在此范围内不可生长在水浇地和大棚地。

(2) 除了水浇地外均不可种植水稻

x ijk =0,i=1,2,,26;i=35,36,,54;j=16

(3) 水浇地和普通大棚第二季不可种植的农作物

x ijk =0,i=27,28,,50;j=17,18,,34;k=2

(4) 仅水稻第二季可种农作物

x ijk =0,i27,28,,34;j35,36,37;k2

(5) 仅有智慧大棚第二季种该类农作物

x ijk =0,i51,52,53,54;j38,39,40,41;k2

3.4.2. 土地限制

(1) 露天田地总亩数

根据题目要求,露天总亩数不可超过1201亩。

i=1 34 x ijk 1201

(2) 每个大棚面积限制

x ijk 0.6,i=35,36,,54

(3) 根据附件每个区不可超过区总面积

(a) 对A区而言(A区总面积为365)

i=1 6 x ijk 365

(b) 对B区而言(B区总面积为619)

i=7 20 x ijk 619

(c) 对C区而言(C区总面积为108)

i=21 26 x ijk 108

(d) 对D区而言(D区总面积为109)

i=27 34 x ijk 109

3.4.3. 种植要求

(1) 两季或两年连种

(a) 只针对干旱地、梯田、山坡地的单季作物

x ij0l x ij0( l+1 ) ,i=1,2,,27

(b) 智慧大棚两季种的农作物不相同

x ij1l x ij2l ,i=51,52,,54;j=17,18,,34

(c) 对于水浇地来说,只能种植两季蔬菜类或者单季水稻

x ij0l 0 x ij0( l+1 ) =0, x ij1( l+1 ) x ij2( l+1 ) ,l=0,1,,7

第l年种植单季水稻,那么l + 1年需要不种植单季水稻种植两季蔬菜,且蔬菜两季不相同。

(2) 三年至少种植一次豆类,且所有土地都至少种植过一次豆类

(a) 针对干旱地,梯田,山坡地

x ijkl + x ijk( l+1 ) + x ijk( l+2 ) >0,i=1,2,,26;j=1,2,,5;l=0,1,,5

(b) 针对水浇地和普通大棚

由附件1可以了解到水浇地和普通大棚只有在第一季可以种植豆类植物,第二季不可以,且可种植豆类植物有限只能分为豇豆,刀豆,芸豆这三类。

x ijkl + x ijk( l+1 ) + x ijk( l+2 ) >0,i=27,28,,50;j=17,18,19;k=1;l=0,1,,5

(c) 针对智慧大棚

由附加1可知智慧大棚第一季第二季均有可能种植豆类植物,豆类植物也仅在豇豆,刀豆,芸豆这三类中选择。

x ijkl + x ijk( l+1 ) + x ijk( l+2 ) >0,i=51,52,53,54;j=17,18,19;k=1,2;l=0,1,,5

(3) 种植密度

为了方便对田地的管理且有利于作物生长,每种作物的种植面积不能过小。 x min 为最小的种植面积。

x ijkl x min ,i,j,k,l

3.4.4. 销量与产量

销售量一定要满足销量小于或等于产量的约束条件。

Y ijkl x ijkl × B ijkl

3.5. 结果分析

基于两种不同的情况,并且在满足预期销售产量大于产量的条件下,针对性地定制相应的最优种植方案,以使每种情况下都能达到利润最大化。在满足相关约束条件的情况下,我们使用python编程软件对此问题进行求解并将结果输出到Excel文件中,使用SPSSPRO网站对输出的Excel文件进行绘图展示。

图2图3可以看出在目标函数改变的情况下,种植作物种类的选择以及作物的种植量都会随之发生相应的改变。其中玉米和小麦的波动变化情况较大,但是种植量也最多。上述情况反映,小麦和玉米在露天耕地一年交替一次种植的效果最好,此种情况下所获得的利润也相对较大。

4. 基于蒙特卡洛模拟[6]的动态规划最优种植方案

根据经验,农作物的预期销售量、亩产量、种植成本和销售价格存在不确定性,同时还需要考虑种植风险和相关种植要求。为了实现乡村经济的可持续发展,需要充分利用耕地资源,为此结合动态规划模型和蒙特卡洛模拟,对自变量变化情况进行模拟,从而选择适宜的农作物并实现种植策略的优化。

Figure 2. Results for Case 1

2. 情况一结果展示图

Figure 3. Results for Case 2

3. 情况二结果展示图

4.1. 蒙特卡洛模拟

4.1.1. 阶段划分

由于上一年的销售量、亩产量、种植成本和销售价格会对下一年的数据造成影响,因此在这里把2024至2030年划分为7个阶段,每个阶段为1年。

4.1.2. 状态定义

每年都需要每个阶段各种农作物的种植面积、预期销售量、亩产量、种植成本以及销售价格的状态。

4.2. 动态规划模型的建立

4.2.1. 随机模拟

使用蒙特卡洛模拟来处理农作物的预期销售量、亩产量、种植成本和销售价格的不确定因素。

假设小麦和玉米的预期销售量的平均年增长率介于5%~10%之间,其他农作物的预期销售量年增长率变化区间为±5%,并对其各自变化的区间利用蒙特卡洛模拟进行随机模拟。

假设天气对农作物亩产量的影响在±10%之间,对每个阶段的亩产量利用蒙特卡洛进行随机模拟。此外,假设农作物的种植成本平均每年增长5%,蔬菜类作物的销售价格平均每年增长5%,食用菌的销售价格每年下降1%~5%,特别是羊肚菌的销售价格每年下降幅度为5%,并采用同样的方法对此数据进行随机模拟。

4.2.2. 目标函数和约束条件

目标函数为最大化期望收益,同时考虑潜在的种植风险。收益利用某类农作物某季销售单价和种植亩产量及种植面积的乘积除去种植成本,然后对所有产品收益进行一个求和。以最大化利润为目标函数(利润 = 总收益 − 成本),建立线性规划模型。

maxZ= i=1 54 j=1 41 k=0 2 l=1 7 x ijkl b ijk ( P ijk Q ijk ) +δ P ijk max{ i=1 54 j=1 41 k=0 2 l=1 7 x ijkl B ijk Y ijk ,0 } (2)

其中 δ 为0.5。

根据各变量的随机变化情况得到动态规划模型的约束条件,如下:

(1) 预期销售产量约束

针对玉米小麦:

Y jl = Y j( l+1 ) ×( 1+m ),j=6,7;l=0,1,,7 (3)

其中 Y jl 表示j类农作物在l年的预期销售量,m为利用蒙特卡洛模拟在小麦玉米预期产量销售的年变化率范围内模拟出来的变化值。

针对其他农作物:

Y jl = Y j0 ×( 1+n ),j=1,,5,8,,41;l=1,,7 (4)

n为利用蒙特卡洛模拟在其他农作物预期产量销售的年变化率范围内模拟出来的变化值。

(2) 亩产量

由于农作物亩产量受天气影响,亩产量变化为:

B j1 = B j( l1 ) ×( 1+p ) (5)

p为利用蒙特卡洛模拟在农作物亩产量的年变化率范围内模拟出来的变化值。

(3) 种植成本

每年种植成本逐步以5%增长,公式如下:

Q jl = Q j( l1 ) ×( 1+q ),j,l=1,2,,7 (6)

(4) 销售价格

对于粮食类处于相对稳定状态则不发生变化:

P jl = P j( l1 ) ,j=1,2,,17;l=1,2,,7 (7)

对于蔬菜类每年价格变化公式:

P jl = P j( l1 ) ×( 1+s ),j=17,18,,37;l=1,2,,7 (8)

对于羊肚菌每年价格变化公式:

P jl = P j( l1 ) ×( 1s ),j=41;l=1,2,,6 (9)

对于除了羊肚菌以外的食用菌每年价格的变化公式:

P jl = P j( l1 ) ×( 1t ),j=38,39,40;l=1,2,,6 (10)

t为利用蒙特卡罗模拟在除了羊肚菌以外的食用菌销售价格的年变化率范围内模拟出来的变化值。

4.3. 结果分析

在综合考虑各种农作物的预期销售量、亩产量、种植成本和销售价格的不确定性以及潜在的种植风险的情况下,可以得出在波动因素影响下的作物分布图,如图4所示,小麦和玉米的种植量大幅度减少,由图5可以看出其他作物的变化,如黑豆的波动幅度较小,植被覆盖率最低的地块是A1地块的白灵菇,最高的是F2地块的红豆。整体看来,各个地块间不同种类的豆类或菌类的覆盖率差异不大。

Figure 4. Crop distribution map under fluctuating factors

4. 在波动因素影响下的作物分布图

由于上述图片不易观测出来种植量小的作物种植变化情况,在此对种植量小的作物利用SPSSPRO软件进行绘图展示。

Figure 5. Spatial distribution of leguminous crops

5. 豆类作物分布图

对实验结果进行验证,计算得到粮食类作物利润的变异系数为13.86%,其它类作物的利润变异系数为24.83%。由于传统粮食类作物的变异系数在10%~15%,高附加值经济作物(如蔬菜、菌类)的变异系数在20%~30%,这表明模型是较为稳定的。

5. 基于蒙特卡洛模拟的动态规划最优种植方案

5.1. 作物之间的替代性和互补性

可以利用交叉弹性来判断互为替代品和互补品的两种作物的关系。交叉弹性[7]的计算公式:

AB= % B % A (11)

其中AB为交叉弹性,若AB为正,则表示两个农作物互为替代品,若AB为负,则表示两种农作物互为互补品。%BB作物销售量需求代数,%AA作物价格代数。

下面给出 % B % A 的计算公式:

%B= Δ b ijkl Y ijkl %A= Δ P ijkl P (12)

P 为随机销售价格。销售量和价格的变化是动态规划模型的蒙特卡洛生成的变化率计算得到新的预期销售量和价格减去实际销售量和价格。

5.2. 互补品与替代品判断

根据对上述的计算,我们得到了41类农作物互为互补品和替代品的结果,部分结果如表2表3所示。

Table 2. Crops that are complements to each other

2. 互为互补品的作物

作物

黄豆

红豆

绿豆

南瓜

香菇

西红柿

红薯

玉米

生菜

菠菜

高粱

Table 3. Crops that are substitutes for each other

3. 互为替代品的作物

种植编号

白萝卜

红萝卜

榆黄菇

香菇

白灵菇

黄豆

黑豆

红豆

绿豆

爬豆

表2表3可知,互为互补品的作物中只有部分作物可以满足题中种植地块对种植作物的要求,如玉米和黄豆、玉米和红豆、玉米和绿豆以及高粱和绿豆等,因此其他作物虽然互为互补物也不再纳入考虑范围。互为替代品的作物中仅粮食类作物之间存在可替代性,其它类作物均不能满足种植地块对作物的硬性要求,所以不纳入考虑范围。

假设当互为互补品的作物种植在同一地块时,亩产量增加10%。基于这个假设运用Python代码得到的结果显示为利润不存在明显波动。

5.3. 销售量、价格与成本的相关性

通过SPSSPRO对处理后的总预期销售量数据,即销售单价数据,种植成本数据进行相关性分析,得到的相关系数及相关性强度如图6所示:

Figure 6. Correlation coefficient heatmap

6. 相关性系数热力图

图6可知,总预期销售量与随机销售单价的相关性和总预期销售量与种植成本的相关性最强。因此利用这两个相数引入相关性矩阵。

这里我们利用经济基础上的公式来表示三个变量之间的关系:

P ijkl = P ijkl ×( 1+α× Y ijkl ) Q ijkl = Q ijkl ×( 1+α× Y ijkl ) (13)

其中 P ijkl Q ijkl 为新生成的随机销售单价和成本。 P ijk 表示第i种类型的土地种植第j类农作物在第k季的销售单价, Y ijk 表示第i种类型的土地种植第j类农作物在第k季的农作物预期销售量, Q ijk 表示第i种类型的土地种植第j类农作物在第k季的农作物成本。

5.4. 优化模型建立

  • 目标函数仍然为利润最大化:

maxZ= i=1 54 j=1 41 k=0 2 l=1 7 x ijkl × b ijk ×( P ijk Q ijk ) +δ× P ijk ×max{ i=1 54 j=1 41 k=0 2 l=1 7 x ijkl × B ijk Y ijk ,0 } (14)

其中 δ 为0.5。

  • 约束条件

在动态规划模型的基础上增加相关性的约束条件:

P ijkl = P ijkl ×( 1+α× Y ijkl ) Q ijkl = Q ijkl ×( 1+α× Y ijkl ) (15)

5.5. 结果分析

利用lingo软件进行求解,得到基于作物之间的替代性与互补性的最优化种植方案部分结果展示如表4

Table 4. Crop allocation table based on substitution and complementary effects

4. 经替代性与互补性考虑下的作物种植分布表

种植地块

作物编号

种植季次

年份

种植量

1

6

0

7

80

2

6

0

7

55

3

6

0

7

35

4

6

0

7

72

5

6

0

7

68

6

6

0

7

55

利用SPSSPRO软件绘制出不同地块类型下的作物种植分布。

基于作物之间的替代性与互补性的最优化种植方案得出的结果和动态规划模型结果进行对比分析,可以得出在考虑其他因素情况下的作物的种植情况大体趋势相同,波动变化也大致相同。但是在某些地块上,种植的作物和作物种植量存在较大的差异。并且基于作物之间的替代性与互补性的最优化种植方案得出的大麦和玉米的种植量有所提高。效果展示如图7图8

由图可以看出,结果变化微小,所以建立的线性规划模型比较稳健。

6. 模型检验

利用敏感性分析的方法,选取这种作物的亩产量这一数据,对其值进行增加或者减少,来观测亩产量变化情况对结果的影响,具体结果如下表5所示。

由于改变亩产量的值,利润变化微小。由此得出所建立的线性规划模型较为稳健。

Figure 7. Comprehensive distribution map of crop planting quantities

7. 综合考虑下作物种植量分布图

Figure 8. Legume and root crop spatial distribution

8. 豆类和薯类分布图

6.1. 模型优点

  • 针对不同的情况建立针对性模型,在外部资源和条件不变的情况下采用单阶段的线性规划模型,在外部资源和条件随着时间和空间发生变化采用动态规划模型并结合蒙特卡洛模拟,在考虑作物之间的相关性和替补性时引入弹性系数和相关性来讨论建立优化模型。

Table 5. Sensitivity analysis results

5. 敏感性分析结果

亩产量浮动

Δ利润

Change_−10%

106949.1601

Change_−20%

100909.8701

Change_−30%

94869.58006

Change_10%

119195.3401

Change_20%

125484.2301

Change_30%

131777.1201

  • 对原始数据进行了全面的处理,包括正态性检验、描述性统计和异常值处理。

  • 充分考虑实际情况中各种因素的制约,如地块类型限制、不能重茬种植、豆类作物种植限制、种植面积限制等。

对模型结果分析,该模型能够为乡村制定阶段性最优种植策略,具有一定的实践意义和参考。

6.2. 模型缺点

  • 采用蒙特卡洛模拟对不确定因素进行处理具有一定的局限性。

  • 仅采用亩产量这一个因素变化对模型进行敏感性分析。在实际情况下,我们应该考虑多种因素对模型的影响。

基金项目

河南省“专创融合”特色示范课程项目(2023),河南科技大学教育教学改革研究与实践项目(2024BK089),河南科技大学大学生创新创业计划项目(2024238)。

NOTES

*通讯作者。

参考文献

[1] 仝芮宁. 基于多目标线性规划的濮阳市土地利用优化分析[J]. 国土与自然资源研究, 2023(3): 13-17.
[2] 全国大学生数学建模竞赛组委会. 全国大学生数学建模竞赛赛题, C题[EB/OL].
https://www.mcm.edu.cn/html_cn/node/a0c1fb5c31d43551f08cd8ad16870444.html, 2024-09-05.
[3] 方涛, 姚应生. 夏皮罗-威尔克方法对不等精度观测列正态性的检验[J]. 矿山测量, 1990(2): 17-19.
[4] 张纪泉. 总体分布的正态性检验——介绍夏皮罗-威尔克的W检验法[J]. 中国纤检, 1982(5): 34-40.
[5] 司守奎, 孙兆亮. 数学建模算法与应用[M]. 第2版. 北京: 国防工业出版社, 2011.
[6] 胡雯, 靳芙蓉, 张珍明. 基于蒙特卡洛模拟的农田土壤微塑料阈值及污染风险评估[J/OL]. 中国无机分析化学, 2025: 1-14.
http://kns.cnki.net/kcms/detail/11.6005.O6.20250422.2058.005.html, 2025-05-30.
[7] 肖成浩. 考虑需求交叉弹性和制造商风险偏好的供应链联盟定价决策研究[D]: [硕士学位论文]. 北京: 北京物资学院, 2024.