首页 学术研究∶土地资源利用与保护

模拟退火法优化约束条件下造林规划设计的研究


吴承祯 洪伟

(福建林学院 资源与环境系,福建 南平 353001)

提 要 模拟退火法是一种新的优化算法,文中首次提出用该法确定在资金、计划、苗木品种和数量等约束条件下小班造林规划设计最优方案,并阐述其数学模型和具体方法。笔者用 97 个造林小班资料,应用模拟退火法对在约束条件下的小班造林规划设计进行了方案优化,结果表明该法优化约束条件下造林规划设计效果很理想。

关键词 模拟退火法;造林规划设计

中图分类号 S757.4

文献标识码 A

文章编号 1000-3037(2000)01-0086-05

  造林规划设计是林业建设的先行和基础工作,是科学造林和最优经营管理的重要组成部分。造林规划设计中方案的优化不仅要遵循适地适树的自然法则,而且要根据当地、当时社会经济条件、林业生产计划和经营目的,来选择最优的造林规划设计方案,做到在满足经济、计划、物质等方面的条件下,合理利用林地,最大限度地挖掘林地自然生产力,充分发挥人工林的生产效益,实现森林的持续经营与林业的可持续发展。这是当前造林规划设计中的一个重要研究课题,具有重大现实意义和应用价值。为此,前人从不同的技术角度,对造林规划设计的方案优化做了大量研究,取得了许多重要成果[1-4],但它们均具有一定的局限性。因为在造林规划设计的实践中,小班造林树种的选择除了要考虑某种立地条件下某一树种的适宜程度外,还常常要考虑投资金额的限制、指令性计划的约束以及造林苗木品种和数量的限制。鉴于此,笔者曾首次给出了 0-1 规划方法确定约束条件下小班造林规划设计最优化的具体方法,并成功地解决了 10 个小班的造林规划问题[4]。但当造林小班数量大于 50 甚至更多时,0-1 规划方法求解就需搜索很长时间方有可能得到最优方案。那么,在解决某个大林场或某个县(市)整个大地域范围内的约束条件下小班造林规划设计及其最优设计问题时,就存在最优化方法上的不足。因此,有必要提出优化的方法。

  造林规划设计问题本质上是一个优化问题。模拟退火法是人们从自然界固体退火过程中得到启发并从中抽象出来的一种随机优化算法,是解决优化问题的一种新途径。因此,本文在介绍模拟退火法机理的基础上,首次提出用模拟退火法优化大地域作业区在经济、计划、物质等方面的约束条件下为获得最佳经济效益的小班造林树种选择的最优方案,具有重要理论价值和现实意义,属首次报道。

1 在约束条件下造林规划设计的数学模型

  在实际的造林设计工作中,当用科学方法做出设计之后,具体实施过程又常常受到各种条件的约束,迫使修改原规划设计方案。在满足这些约束条件下,为达到获得最大经济收益目标,编制整个林场或县、市作业区内造林树种选择的最佳方案,我们把它称为“在约束条件下造林规划设计问题”[4]。在约束条件下造林规划设计考虑资金、苗木品质和数量以及计划等约束条件,其具体数学模型见文献 4。

2 模拟退火法

  Kirkpatrick 等[5]于 1983 年首先提出了模拟退火法,它是人们从自然界固体退火过程中得到启发并从中抽象出来的一种随机优化算法。模拟退火法用于求解优化问题的出发点是基于物理中固体物质的退火过程与一般优化问题间的相似性。在对固体物质进行退火处理时,常先将它加温使其粒子可自由运动,以后随着温度的逐渐下降,粒子逐渐形成低能态晶格。若在凝结点附近的温度下降速率足够慢,则固体物质定会形成最低能量的基态,优化问题也存在类似过程。解空间中每一点代表一个解,不同的解有不同的目标函数。该算法最为显著的特征是以一定的概率接受使目标函数值增大的移动,所以能够从局部最优解的“陷阱”中爬出来而不会简单地终止于一局部最优解上,即具有全局收敛性。并且在理论上已经证明了只要系统过程满足一定的要求(系统温度无限趋于零度且在每一温度下模拟充分),则算法将以概率 1 渐近收敛于全局最优解[6]。图 1 为典型的模拟退火法计算框图,可以看出,该算法包括以下要素:

nr20000115-1.gif

图 1 模拟退火法计算框图
Fig.1 Calculation chart of simulated annealing

  (1) 相邻状态的产生 相邻状态是指从当前状态经过一次移动(对任一优化变量进行一次扩大、缩小或对流程结构进行一次调整)所能达到的状态。这是模拟退火法中一个非常重要的概念,因为相邻状态产生的有效性直接关系到算法的有效性。一般的相邻状态产生方法是借助于一个随机数发生器,随机地选取一优化变量,然后对其进行放大或缩小。但是此法存在一些不足之处,本文采用这样一种改进的相邻状态产生方法[6]:首先,将当前状态的相邻结构集(相邻状态的集合)进行排列,然后按照这一随机排列顺序逐一执行各步移动。这样保证了当前状态的每一相邻状态均能被实现一次,也使得低温下每一状态仅被实现一次便可找到能够使目标函数值下降的移动,节省了计算时间。

  (2) 初始温度 T0 温度 T 在模拟退火法中具有决定性作用(称为模拟退火法的控制参数),它直接控制着退火的走向(即系统的优化方向)。由随机移动的接受准则可知,T 很大时新状态的接受概率很高。但初始温度 T0 的选取不能过高也不能过低,过高则以后的过程会有大量的时间浪费在因初始温度过高而接受的使目标函数值上升的移动上;T0 过低又会使算法的“爬山”能力减弱而可能终止于局部最优解。一般的 T0 确定方法是使初始温度 T0 下随机移动的接受比率落在某一给定的范围内(例如:0.8≤X0≤0.9)。确定 T0 的经验法则是:选定一个大值作为 T0 的当前值,并进行若干次变换,若接受比率 x 小于预定的初始接受比率 X0(可取 X0=0.8),则将当前 T0 值加倍。以 T0 新的当前值重复上述过程,直至得到X>X0的 T0 值。

  (3)每一温度 T 下随机移动的次数 每一温度下随机移动的次数取为相邻结构集尺寸的大小。

  (4) 降温进程 Tk+1=f(Tk) 选用 Aarts 和 Van Larrhoven[6]提出一种具有多项式收敛的降温进程:

Tk+1=Tk[1+Tk1n(1+δ)/3σ(Tk)]-1

其中σ为目标函数标准偏差,δ为给定数值,δ越小则降温越慢但模拟误差也越小,其取值不同对模拟退火算法的实验性能有显著影响。一般通过计算机调试而确定,本文经调试,δ取 0.10。

  (5)随机移动接受准则 采用 Metropolis 准则[7、8],随机移动接受概率

nr20000115-4.gif

其中△c 为目标函数值的变化。可以得知,系统温度 T 决定着随机移动的接受概率。温度越高则算法接受使目标函数值上升移动的能力越强,具有较强的“爬山”能力;温度很低则使目标函数值上升移动的接受概率很低。

  (6)终止准则 采用一个简单的终止准则:当 T=Tf 时终止计算(取 Tf=0.1)。

3 模拟退火法优化约束条件下造林规划设计实例

  用福建省尤溪县林业委员会用材林速生丰产林基地中的造林小班的具体资料进行在约束条件下的小班造林树种选择方案优化。

  福建省尤溪县位于东经 118.8°~119.6°,北纬 25.8°~26.4°。气候属中亚热带大陆性兼海洋性季风气候,年降水量 1599.6mm,年蒸发量 1323.4mm,相对湿度 83%,年平均气温 18.9℃,历年最大日降水量 131.7mm,3~6 月为多雨季节,4 个月降雨量占全年降雨量的 65%。地貌系闽中火山岩系中山地貌。用材林速生丰产林基地造林地土壤为发育在侏罗纪下流陆相盆地沉积砂岩上的红壤,平均坡度为 32°。林地植被以唐竹、芒萁为主,盖度为 63%。造林小班共 97 个,总面积为
618.20hm2,依据各小班坡位、坡向、坡度、表土层厚度等因素,应用文献 3 中的方法可以打印出各树种在基准年时累计材积和合理造林密度。由于用材林速生丰产基地造林树种主要是杉木(Cunninghamia lanceolata)和马尾松(Pinus massoniana),因此只考虑这两个树种,这些基础数据限于篇幅不列出。

  假设杉木每株苗木单价 (H1) 为 0.1 元,马尾松每株苗木单价 (H2) 为 0.05 元。杉木每公顷造林投资 (b1) 为 2700 元,马尾松每公顷造林投资 (b2) 为 1800 元(包括整地、造林及抚育等费用)。各树种的总材积(包括间伐材)平均每立方米杉木价格 (C1) 计为 480 元,马尾松 (C2) 为 360 元。假设对该作业区总投资为 160 万元 (B),对杉木苗木来说,最多可供应 180 万株 (D),但至少必须使用 120 万株 (E),马尾松苗木敞开供应。上级部门要求 20 年后该批造林小班必须收获杉木 4.5万m3 木材。以立地条件分析:第 5、59、97 小班极适宜种杉木;而第 10、51、96 小班极适宜种马尾松。这些条件在小班造林规划设计时都予以考虑。这样,目标函数为:

nr20000115-5.gif

约束条件:

nr20000115-6.gif
nr20000115-7.gif
nr20000115-9.gif
nr20000115-10.gif
nr20000115-11.gif

其中 Z 为基准年时创造价值(元);Vij 为第 i 小班种植第 j 树种在基准年(20 年)时累计材积(按文献 3 计算);Si 为第 i 小班面积;Cj 为第 j 树种(包括间伐和主伐收获的木材)每立方米平均计算价格;Xij 为第 i 小班种植第 j 树种的决策变量;Nij 为第 i 小班种植 j 树种的合理密度(需苗木株数)(按文献 3 计算);n 为小班数;m 为造林树种数。

  笔者根据模拟退火法的基本思想,应用 BASIC 语言编制了优化约束条件下造林规划设计的模拟退火法的计算机程序,其中 n=97,t0=100,δ=0.10,马尔可夫链的长度取为 (2n)2。程序在 AST-586 计算机上运行,收敛后,目标函数值趋于稳定,最优造林规划设计方案列于表 1。在采用模拟退火法优化约束条件下造林规划设计求解的尤溪县林业委员会 97 个小班的造林规划设计最优组合方案下,20 年累计杉木材积收获可达 4.72万m3,杉木苗木需要 1204010 株;马尾松材积收获为
3.17万m3,需要苗木 1150560 株;总投资为 1597349 元。在最优造林规划设计方案下,可望在 20 年时获得 3409.2 万元经济收益。

表1 造林规划设计表
Table 1 Tabulated culture program design



小班
面积
(hm2)
造林
树种
合理
密度
(株/hm2
需要苗
木数量
(株)
20 年木
材收获
(m3)


小班
面积
(hm2)
造林
树种
合理
密度
(株/hm2
需要苗
木数量
(株)
20 年木
材收获
(m3)
1 3.80 杉木 4200 15960 576.56 50 6.40 杉木 3300 21120 916.13
2 8.06 杉木 3600 29040 1316.24 51 3.53 马尾松 4800 16960 474.99
3 1.40 杉木 3600 5040 242.40 52 5.00 马尾松 4800 24000 817.20
4 5.00 杉木 4200 21000 747.60 53 11.00 杉木 3000 33000 1060.79
5 8.80 杉木 3600 31680 1478.91 54 18.40 马尾松 3600 66240 1405.94
6 4.19 杉木 3600 15120 471.64 55 2.53 杉木 4350 11020 496.89
7 3.00 杉木 3300 9900 429.35 56 2.46 马尾松 4800 11840 358.49
8 2.00 杉木 3300 6600 286.23 57 6.93 杉木 3600 24960 821.18
9 4.73 杉木 3600 17040 743.51 58 5.60 马尾松 4200 23520 654.61
10 6.46 马尾松 4800 31040 989.30 59 3.06 杉木 3900 11960 566.58
11 5.86 杉木 4200 24640 882.99 60 6.46 马尾松 4800 31040 902.68
12 9.13 杉木 3600 32880 1557.83 61 6.00 马尾松 4350 26100 703.89
13 9.39 杉木 3600 33840 1522.38 62 2.53 马尾松 4800 12160 316.20
14 5.73 杉木 3300 18920 806.51 63 7.26 马尾松 4350 31610 854.78
15 3.26 马尾松 4500 14700 398.08 64 8.33 杉木 3000 25000 917.13
16 2.33 马尾松 4200 9800 295.09 65 9.46 马尾松 3900 36920 996.98
17 8.80 马尾松 4050 35640 986.96 66 3.80 杉木 3300 12540 502.51
18 4.46 杉木 3000 13400 566.62 67 5.13 杉木 4200 21560 947.56
19 4.93 杉木 4200 20720 729.71 68 4.33 马尾松 4800 20800 635.25
20 3.60 马尾松 4200 15120 443.72 69 9.80 马尾松 4200 41160 1214.81
21 12.86 杉木 3000 38600 1494.98 70 11.26 马尾松 3900 43940 1217.81
22 6.26 马尾松 4200 26320 839.51 71 2.80 杉木 4200 11760 436.59
23 5.93 杉木 3000 19580 814.97 72 5.80 杉木 3900 22620 860.78
24 1.80 马尾松 4800 8640 270.97 73 7.66 杉木 3600 27600 942.54
25 7.73 马尾松 4050 31320 869.07 74 1.60 杉木 4350 6960 317.26
26 5.80 马尾松 3900 22620 612.13 75 4.93 杉木 3900 19240 806.26
27 6.06 马尾松 4500 27300 731.28 76 3.86 马尾松 4200 16240 491.61
28 10.40 杉木 3300 34320 1480.91 77 4.60 马尾松 4500 20700 595.40
29 3.93 马尾松 4200 16520 476.48 78 2.93 杉木 3600 10560 439.38
30 9.26 杉木 4200 38920 1381.38 79 8.26 杉木 3300 27280 1016.30
31 5.06 杉木 4200 21280 833.87 80 2.73 马尾松 4800 13120 423.08
32 3.13 马尾松 4800 15040 488.61 81 3.80 马尾松 4800 18240 545.32
33 11.20 杉木 3000 33600 1242.69 82 5.46 杉木 4200 22960 864.53
34 4.26 杉木 2700 11520 424.26 83 6.40 杉木 4050 25920 938.30
35 4.86 马尾松 4050 19710 599.77 84 4.86 马尾松 4500 21900 684.67
36 2.86 杉木 3000 8600 351.27 85 8.60 马尾松 4200 36120 1015.75
37 11.53 杉木 3600 41520 1582.26 86 9.53 杉木 3000 28600 1124.27
38 13.53 杉木 3000 40600 1561.27 87 2.46 杉木 4350 10730 484.40
39 12.33 杉木 3000 37000 1351.61 88 5.20 杉木 4200 21840 843.88
40 2.40 杉木 3300 7920 259.93 89 8.93 杉木 4200 37520 1314.81
41 9.80 马尾松 3600 41160 1102.21 90 7.60 杉木 3900 29640 953.84
42 4.19 杉木 3000 12600 754.49 91 11.46 马尾松 3900 44720 1038.36
43 5.40 马尾松 4200 22680 746.50 92 4.33 杉木 3900 16900 531.64
44 8.46 马尾松 3900 33020 877.06 93 7.13 马尾松 3900 27820 655.59
45 5.20 马尾松 4350 22620 659.99 94 3.20 马尾松 4800 15360 499.01
46 17.06 马尾松 3900 66560 1591.55 95 9.26 马尾松 3600 33360 801.75
47 5.06 杉木 4200 21280 818.44 96 15.80 马尾松 3600 56880 1447.36
48 7.06 杉木 3000 22800 957.14 97 3.40 杉木 3900 13260 419.99
49 7.13 杉木 3300 23540 906.18                  

4 讨论

  在约束条件下造林规划设计方案的优化是在反映立地条件的环境因子等估测出最佳造林树种及合理造林密度的计算机辅助造林设计的基础上,应用其定量数据,用模拟退火法优化 0-1 规划求解最优设计方案,达到在大区域范围内实现在约束条件下的造林规划设计,因此,研究结果对实现合理配置林地资源具有一定的现实指导意义。

  在约束条件下造林规划设计中采用模拟退火法优化求解,结果表明,模拟退火法可用来求解优化问题,但在实现过程中应注意 3 个问题:①怎样按某种概率过程产生新的搜索状态。产生新的搜索状态的范围越大、越广,则搜索的解空间越大,从而能够保证在全局中寻求最优,因此,为保证在全局中寻求最优,本研究产生新的搜索状态的概率过程采用随机生成方式;②根据当前温度及新状态与原状态的相应位置,如何确立新状态接受标准。本研究中采用 Metropolis 准则来确定新状态接受概率,即同时考虑目标函数值的变化Δc 和系统温度 T;③怎样选择初始温度 T0 及怎样更新温度、确定温度的下降过程。一般初始温度 T0 可经过计算机调试得到,本研究中经计算机调试后以 T0 取 100 为最理想;而对于温度下降过程,由于降温过程的数学模式较多,不同的研究对象可因具体情况选择不同的降温过程,本研究选择 Aarts 和 Van Larrhoven提出的一种具有多项式收敛的降温进程。以上 3 点影响模拟退火法的收敛性及收敛速度,且影响退火结束后以多大的概率使状态稳定在全局最小点。作为一种新的优化方法,其将在林业生产及科学研究中具有广阔的应用前景。

参考文献

1 徐智.造林规划中适地适树的经济评价方法探讨[J].林业经济,1984,(5):20~28.
2黑川泰亨.林业经营计画にょ为目标计画の利用にっ∽こ[J].日林志,1981,63(4):144~149.
3俞新妥,林思祖.计算机辅助造林设计系统[J].林业科学,1986,22(4):337~346.
4洪伟.在约束条件下造林设计的研究[J].福建林学院学报,1987,7(3):57~64.
5Krikpatrick s Gelett C, Veechi M. Optimization by simulated annealing [J]. Scie nce, 1983,220:671~680.
6Van Larrhoven P J M, Aarts E H L. Simulated Annealing: Theory and Application [M]. Netherland: Klu wer Academic Publisher,Dordrecht, 1987.
7Metropolis N, Rosenbluth A W, Rosenbluth M N. Equation of state calculations by fast computing mach ines [J]. J Chem phys, 1953,21:1087~1092.
8Emile Aarts, Jan Korst. Simulated Annealing and Boltzman Machines [M]. Chichester: John Wiley&Sons, 1990.6~17.

第一作者简介

吴承祯 (1970-),男,江西吉安人,副教授,在职博士,主要从事数量生态和林业系统工程方面的研究。

基金项目:福建省自然科学基金资助项目 (F991)。


本篇文章还收录在以下主题中:科学与技术
生态系统∶植被∶植树造林
关闭窗口