96年A题最优捕鱼计划

1、课程名称：数学建模班号 学号： 姓名:成绩:任课教师：郑勋烨日期：2011 年 5 月 29 日最优捕鱼摘要：为了保护人类赖以生存的自然环境，可再生资源（如渔业、林业资源等）的开发必须适度。在经济生活中，对于“经济人”而言都是追求利润最大化，因此，最优捕鱼策略不仅要求我们要考虑商家的经济效益，还要考虑自然环境的生态效益。一种合理、简化的策略是，在实现可持续收获 的前提下，追求最大产量或最佳效益。本文就渔业的可再生问题进行研究，建立模型，并利用相关 的 lingo 和 matlab 等数学软件对模型进行求解。关键词：最优捕捞策略生产能力成活率季节性集中产卵捕捞强度系数课程名称：数学建模班号 学号

2、： 姓名:成绩:任课教师：郑勋烨日期：2011 年 5 月 29 日一. 问题重述为了保护人类赖以生存的自然环境，可再生资源（如渔业、林业资源）的开发必须适度。一种 合理、简化的策略是，在实现可持续收获的前提下，追求最大产量或最佳效益。考虑对某种渔（鱼是鱼）的最优捕捞策略：假设这种鱼分 4 个年龄组，称 1 龄鱼，1 龄鱼。各年龄组每条鱼的平均重量分别为5.07,11.55,17.86,22.99（ 克），各年龄组鱼的自然死亡率均0.8（1/年），这种鱼为季节性集中产卵5繁殖，平均每条 龄鱼的产卵量为1.109 10（个），龄鱼的产卵量为这个数的一半，龄鱼和 龄鱼不产卵。产卵和孵化期为每年的最

3、后4 个月，卵孵化并成活为 1 龄鱼，成活率（龄鱼条数与产卵总量 n 之比）1.22 1011/（1.22 1011n）。渔业管理部门规定，每年只允许在产卵孵化期前的8 个月进行捕捞作业。如果每年投入的捕捞能力（如渔船数、下网次数等）固定不变，这时单位时间捕捞量将与各年龄组鱼群条数成正比，比例系数不妨称捕捞强度系数。通常使用 13mm 网眼的拉网，这种网只能捕捞龄鱼和龄鱼，其两个捕捞强度系数之比为 0.42:1。渔业上称这种方式为固定努力量捕捞。1）建立数学模型分析如何实现可持续捕获（即每年开始捕捞时渔场中各年龄组鱼群条数不变），并且在此前提下得到最高的年收获量（捕捞总重量）。2）某渔业公司承

4、包这种鱼的捕捞业务年，合同要求年后鱼群的生产能力不能受到太大破坏。已知承包时各年龄组鱼群的数量分别为：122,29.7,10.1,3.29（09条）。如果仍用固定努力量的捕捞方式，该公司应采取怎样的策略才能使总收获量最高。综上所述，原问题实质上是给出了各年龄组鱼群之间数量的变化规律，并给出了它们的自然死 亡率及捕捞和产卵的时间分布，并固定3、4 龄鱼捕捞能力的比值，要求选择一定的捕捞能力系数，使得各年龄组鱼的数量在各年开始的第一天条数不变（第一问），5 年后鱼群的生产能力不会有太大的破坏（第二问），并在此条件下，求到最大捕获量。二. 模型假设1.这种鱼在一年内的任何时间都会发生自然死亡，即死亡

5、是一个连续的过程；2.捕捞也是一个连续的过程，不是在某一时刻突然发生；3.鱼群的死亡率已考虑种群的相互竞争及环境等因素；4.3、4 龄鱼产卵集中在 9 月初期，到第二年初完成孵化；5.龄鱼到来年分别长一岁成 i+1 龄鱼，其中上一年存活下来的4 龄鱼仍是 4 龄鱼。三. 符号定义t-时间过程（年）T年份课程名称：数学建模班号 学号： 姓名:成绩:任课教师：郑勋烨日期：2011 年 5 月 29 日 =t-时间间隔M（t）i龄鱼在时间t的数量（个）r-自然死亡率（0.8/年）n-年产卵总量a-4 龄鱼的产卵量（1.109 105），3 龄鱼为a/2b-卵的成活率1.22 1011/（1.22 1

6、011n）k44龄鱼的捕捞强度（k3：3龄鱼为0.42k4）决策变量四. 模型的建立与求解4.1 .问题分析4.1.1. 题目简述及解释:1）题目简述鱼的基本情况：年龄组1 龄鱼2 龄鱼3 龄鱼4 龄鱼平均重量（克）5.0711.5517.8622.99产卵量00a/2a（个）a）各年龄组鱼的自然死亡率均为0.8（1 /年）b）产卵和孵化期为每年最后4 个月成活率=1 龄鱼条数一 1.22 1011成率产卵总量 n 1.22 1011n捕鱼的限定：a）渔业管理部门规定：只允许在产卵孵化期前的 8 个月进行捕捞作业b）固定努力量捕捞：投入的捕捞能力固定不变，这时单位时间捕捞量将与各年龄组鱼群条数

7、成正比，比例系数不妨称捕捞强度系数。 使用 13mm 网眼的拉网，只能捕捞 龄鱼和 龄鱼，其两个捕捞强度系数之比 为 0.42:1课程名称：数学建模班号 学号： 姓名:成绩:任课教师：郑勋烨日期：2011 年 5 月 29 日2）自然死亡率：由题目及假设给出的鱼的自然死亡率是指平均死亡率，即单位时间鱼群死亡数量与现有鱼群数量的比例系数。由死亡是一个数量的连续变化过程，同时1 , 2 龄鱼在每年数量的变化及 3, 4 龄鱼在后 4 个月的数量的变化只与死亡率有关。由此可知，各龄鱼的变化满足：课程名称：数学建模班号 学号： 姓名:成绩:dNi(t)N t,i十,3,4; dt3)捕捞强度系数：单位

8、时间内捕捞量与各年龄组鱼群条数的比例系数，由于只能对捕捞期内进行捕捞。所以，捕捞强度系数影响了 3, 4 龄鱼在前八个月内内的数量变化:设k4= k，则k3= 0.42k故有dN/t)dt4)卵的成活率：1，2 龄鱼不产卵，3，4 龄鱼在每年的后四个月产卵，假设了在9 月初一次产卵，因此可将每年的产卵量 n 表示为：心.109 1055N3(|)吨)11, 1.22 10b=-H又有成活率1.22 10，所以每年年初 1 龄鱼的数目为1.221011111.221011n4.1.2.问题一分析：对于问题一，商家的目标是得到最大的捕捞量，故由捕捞量得到目标函数。又要实现可持续捕 获，即每年开始捕

9、捞时渔场中各年龄组鱼群条数不变。可由每年年初鱼的数目不变得到约束条件。 进而求最优化问题。对于捕捞量，只能对 3，4 龄鱼在 1 到 8 月进行捕捞。只求前八个月的对3、4 龄鱼的捕捞量即为全年的捕捞量；对于鱼的数量变化，1，2 龄鱼不能被捕捞，仅受自然死亡率的影响；而3，4 龄鱼不仅受自然死亡率的影响，还受捕捞强度系数的影响。对于得到的最优化问题，可以用lingo 软件和 matlab 软件进行求解。4.1.3.问题二分析：问题二中，合同要求 5 年后鱼群的生产能力不能受到太大破坏，又要使总收益最高，这就 有可能发生满足了前者满足不了后者之类的情况。对该问题，我们可以考虑在 5 年的合同期内

10、，捕捞强度不变，先确定一个收益最高的捕捞策略， 再计算在这种策略下，5年后各个鱼群的数量，进而确定是否对生产能力产生了太大破坏。但从理 论分析可知，5 年后将在鱼群尽可能接近可持续鱼群的情况下来使捕捞量达到最大。(1)3，4 龄鱼进行捕捞，并且只在则有(4)课程名称：数学建模班号 学号： 姓名:成绩:对于破坏大小，由于各龄鱼群的数量在很大程度上受承包初1 龄鱼影响，我们以第六年初1 龄鱼群数量的变化量与承包时鱼群数量初值之比表示。不妨假定该比值小于5%寸，即可认为鱼群的生任课教师：郑勋烨日期：2011 年 5 月 29 日课程名称：数学建模班号 学号： 姓名： 成绩:日期：2011 年 5 月

11、 29 日产能力没有受到太大破坏。同理可按照问题一中的解法进行求解。4.2 .模型建立4.2.1 问题一模型1可知 1 , 2 龄鱼的生长只受自然死亡率的影响，由此可知1 , 2 龄鱼的生长的微分方程满足方程(1):咽一啊=1,2;dt可得：Ni(t)=Noe ,i ,2;每年年初龄鱼的数量。(5)T 年的 i 龄鱼在 T+1 年变为 i+1 龄鱼,则有，NMT+1)=eNi(T)2 而对于 3, 4 龄鱼，在前八个月生长不仅受自然生长率和捕捞强度系数的影响，而后四个月 仅受自然生长率的影响。以一年为一个时间单位，则这一时间单位可以分为前八个月和后四个月两个阶段。1)前八个月 3、4 龄鱼生长

12、的微分方程满足：dNi=_(r +ki)Nidt可以得到，汕二Nhk|)ti =3,4由于每年的捕捞只在 1 到 8 月进行,为kiNi(t)，则 i 龄鱼的年捕捞量为:,N0为年初龄并且只能捕到 3,鱼的数量(8)4 龄鱼，所以任意一个时刻的捕捞量03KNi(t)dt 二宀 N(1-eFki r)0kr2)后四个月 3、4 龄鱼生长的微分方程满足方程(1)：dNi-rNi,i =3,4dt可得：Ni(t) = Noe_rt,i=3,4;每年年初 龄鱼的数量22.(r .42k)(Hk)其中产卵量 n =0.5aN3e3- aN4e311又孵化存活率苗所以年初 1 龄鱼的总量汕厂1)n(T)3

13、 根据以上分析，我们可以建立非线性规划模型:2 2目标函数：喰十86。际3dt 22.99际4(9)(10)(11)(12)(13)(14)约束条件：任课教师：郑勋烨课程名称：数学建模班号 学号： 姓名： 成绩:日期：2011 年 5 月 29 日课程名称：数学建模班号 学号： 姓名:成绩:2 2一(0.42k卡)一(k卡)n = 0.5aN3e3+aN4e3N =b nN2=b neN3= bne_(2/3 0.42k 3r)e_ /处r)4.2.2 问题二模型针对渔业公司的 5 年捕捞计划，我们利用已得到的迭代方程在已知各个年龄组的鱼的初始值的 前提下，可迭代求出各龄鱼群第 i年的鱼量的分

14、布的函数。n(t 1) =a/2N3(t)e3 r)2/3aN4(t)e3 r)2/31.22011n(t +1)Ni(t 1)右1.2210 n(t 1)rN2(t 1)=eNdt)N3(t+1)=eN2(t)N4(t+1) = N3(t)e+ N4(t)e4& =0.42k4同时写出目标函数:4.3 .模型求解4.3.1 问题一求解4.3.1.1 由 4.2.1 中的 3,我们可将目标函数和约束条件转化为： 目标函数为：2f220.42k丄6(0.8七42k)k(丄282 4)/(吋8)二(0.8*)max =17.86bne (1-e3) 22.99bne228k4).(1-e3

15、)(1-e3)0.42k+0.80.8 + k约束条件:11,1.22汇10b =-11-N4=bn整个生存过程满足的关系式为5max=C N3(t)17.86t吕丄(1e 小r)2/3) C N4(t) 22.99k3rtdk4k4r(1-(k4-r)2/3e(15)课程名称：数学建模班号 学号： 姓名:成绩:1.2210+ n64_(Mk北.28k屉11-(0.28 k七；)e33n =1.22 1011(a (0.5 e32-) -1)-(-0.8)1-e34.3.1.2 然后我们利用 lingo 软件和 matlab 软件分别进行求解。1） lingo 程序：（附录 1）任课教师：郑勋

16、烨日期：2011 年 5 月 29 日课程名称：数学建模班号 学号： 姓名:成绩:直接运行，输出结果为:Local optimal soluti on found at iteratio n:92Objective value:0.3887076E+12Variable Value Reduced CostK 17.36292-1.034723N 0.6078067E+130.000000RowSlack or SurplusDual Price10.3887076E+121.00000020.0000000.1258406E-02k =17.36292,f = 0.3887076E+122）

17、 matlab 程序：（附录 2）首先，建立 fun 1.m、max1.m 1.m、picture.m 文件，运行 picture.m 像：由图像，我们看出kL川一，这与事实是相符的。然后，在命令窗口中输入1.m 文件中命令，运行出结果我们观察到k17,20，fL，因此，我们改变 k 的取值区间即步长，建立了文件 确的结果。从运行结果看：k =17.3000 或 k =17.4000，取最大值 f =3.8871e+011文件，画出 n 关于 k 的图k1,17）f2.m、3.m,从而求得更精课程名称：数学建模班号 学号： 姓名:成绩:最后，建立主程序 3.m17.3,17.4（程序见附录 2

18、）任课教师：郑勋烨日期：2011 年 5 月 29 日课程名称：数学建模班号 学号： 姓名： 成绩:任课教师：郑勋烨日期：2011 年 5 月 29 日从运行结果看：k17.3,17.4, f =3.8871e+011问题一所求结果为：匕17.3,17.4 时，f 取最大值 3.8871&0114.3.1.3 结果分析113 龄鱼的捕捞强度为 7.29/年；4 龄鱼的捕捞强度为 17.36/年：最优可持续捕捞量3.8871*e,可持续捕捞的鱼群大小（条数）：1龄1.1960 10112龄5.3740 10103龄2.4147 10104龄8.3967 107分析结果发现，4 龄鱼在年末

19、存活的数量占全部数量的比例相对很小。4.3.2 问题二求解4.3.2.1 将目标函数转化为:4.3.2.2 用 matlab 软件求解1、利用 matlab 软件，建立 fun.m 文件（见附录 3）；2、 首先，我们对k 1,20，建立 main.m 主程序（见附录 3）0.42k身F书.42kmax1 -e30.42k+rI d+忌扌*为.29+1%-f2r k i-f3r 0.42k j-r3r;：E 0.42k23.29e.3122e.329.7e.31.22 1011n 2_3r 3 01.22 1011n 2e122e11C21.22X10 Xn(2)(10.1 29.7e122e

20、右e1.22 10 n 21121.22 10 n 3牛*e+11e )1.22 10 n 2 J 1.22 10 n 3人3.29e弋挣丿+29.7号.42k10.1e弋3 3k10.1e-3r20.42k-3k-3.29-3r2k329.7e.32P.0.423.4r!：；2：0.42k!：k-f3n：；20.42k -2k“3片10.9 i3丿2(0.42kr)2(k r)3.29 m e340.28 k r)n(3) =29.7 a/2 e3这样 max 就变为关于 k 的函数,-(0.28 k-k5r)10.1 a e3 3_(4k Wr)3.29 a e3 3易于求解。课程名称：数

21、学建模班号 学号： 姓名： 成绩:任课教师：郑勋烨日期：2011 年 5 月 29 日根据运行结果，我们观察到其中，(1,18) fL18,20)fL课程名称：数学建模班号 学号： 姓名:成绩:=19, f19f19：怙：f18接着，改变 k 的取值区间及范围，建立主程序main1.m,main2.m，求解出更精确的结果。运行结果分析:k17.5,17.8, f =1.6056e+012问题二所求结果为：17.5,17.8 时，f 取最大值 1.605610124.3.2.3、验证 5 年后鱼群的生产能力有没有受到太大破坏 迭代求得第六年初各龄鱼群的数量为：11N1(6) =1.1956 10

22、10N2( 6)=5.3713 10N3(6) =2.4157 1010山(6) =8.2710 107第一年各龄鱼群的数量为第六年 1 龄鱼数量占第一年 1 龄鱼数量的比例为:4.324、结果分析捕捞强度在区间（17.5 , 17.8 ）内时（因为电脑精确度问题，暂时只能精确到这一区间）12捞量达到最大值1.6056 10。在这种捕捞强度下，5 年后 1 龄鱼数量占第一年 1 龄鱼数量的比例为，即可认为生产能力没有受到太大破坏。因此，求解出的结果即为最优捕鱼策略。五. 模型的评价本文采用了非线性规划的思想建立模型，通过求解有约束的非线性最大值问题，找到一组最优 解。问题一，在实现可持续捕获

23、（即每年开始捕捞时渔场中各年龄组鱼群条数不变）的前提下，用固定努力量的捕捞方式，确定捕捞策略以得到最大捕捞总质量。用微分、积分的方法来分析每年各龄鱼的数量，建立每年捕捞量的方程，用 软件分别求解，两个结果误差很小，肯定了结果的正确性。问题二，所求模型为五年鱼群生长模型的组合。由于所给的初始鱼群并不是可持续捕捞的鱼群，=17,f1712= 1.6051* e= 1.6053* e12= 1.6020* e1211N1(11.22 1010N2(1)=2.97 10N3(1) = 1.01 1010N4(1)=3.29 109N1(6)100% =98%Ni(1)，总捕lingo 软件与 matl

24、ab课程名称：数学建模班号 学号： 姓名:成绩:为了在五年内既得到最大的收益，又不破坏鱼群的生产能力，即五年后在达到产量最高的条件下使 得鱼群尽量接近可持续捕捞鱼群。 我们在五年内以同样的强度实现固定努力量的捕捞。 对于每年每 条龄鱼在每个时刻的条数， 我们可以用算法迭代求解出 n年的条数，从而比较第六年年初与初始时刻 任课教师： 郑勋烨日期：2011 年 5 月 29 日课程名称：数学建模班号 学号： 姓名： 成绩:任课教师：郑勋烨日期：2011 年 5 月 29 日条数的差值，得出生产能力的破坏度不显著。本模型采用连续模型的方法，成功地解决了可持续捕捞问题，得到了较为精确且合理的结果。六.

25、 模型的改进我们知道，原题中没有说明四龄以上的鱼如何处理。我们假设的是上一年存活下来的 4 龄鱼仍是4 龄鱼，而事实上还可以假设这种鱼只活到4 龄，以后它就死掉了。这对模型没有太大的差别，只是我们所做的假设的分析计算稍复杂，但计算结果也只是稍有差别，在我们模型的基础上，我们可以 假设鱼只能活到 4 龄，这样计算更简便一些。2 2max =17.86 ：k3N3(t)dt 22.99 jk4N4(t)dt约束条件:工（0.42 k十）上（k+）n = 0.5aN3e3+aN4e3汕=b n-N2=bneN3=b ner_（%0.42kd3r）N4二bne3我们可将上面的目标函数和约束条件转化为：

26、目标函数为：约束条件:111.22X0b1.22 1011n11-(0.28k戲)-(?k 0.28k吏)n -1.22 10 (a (0.5 e3e33)-1)用 lingo 软件进行求解，算法见附录4： lingo2 算法直接运行得：Local optimal soluti on found.Objective value: Exte nded solversteps: Total solver iteratio ns:目标函数:max =17.860.42 k0.42k 0.8bn e.6(1(0.8 0.42k)+22.99kk 0.820.42k 3r)bne32-2(0.&

27、：；k)(e30.3887076E+125209课程名称：数学建模班号 学号： 姓名： 成绩:任课教师：郑勋烨日期：2011 年 5 月 29 日K 17.362930.000000Row Slack or Surplus10.3887076E+1220.0000000.1258410E-02即k =17.36293, f = 0.3887076E12而以我们的假设算得出的结果为k=17.36292,f=0.3887076E+12两个结果相差甚小，但改进的模型计算非常简便。七.参考文献赵静，但琦，数学建模与数学实验（第 3 版）高等教育出版社姜启源、谢金星、叶俊，数学模型（第三版）高等教育出版

28、社，2003刘来福，最优捕鱼策略问题答案评述，数学的实践与认识,1997 年 1 月八.附录附录 1:lin go1 算法：max=17.86*0.42*k/(0.8+0.42*k)*1.22*10A11/(1.22*10A11+ n)* n*exp(-1.6)*(1-exp(-2/3*(0.8+0.42*k)+22.99*k/(0.8+k)*1.22*10A11/(1.22*10A11+ n)* n*exp(-0.28*k-2.4)/(1-exp(-2/3*k-0.8)*(1-exp(-2/3*(0.8+k);n=1.22*10A11*(1.109*10A5*(0.5*exp(-0.28*k

29、-6.4/3)+exp(-(0.28+2/3)*k-8.8/3)/(1-exp(-2/3*k-0.8)-1);附录 2:用 matlab 求解，算法如下：picture1.mk=li nspace(1,20,20);n=1.22*10A11*(1.109*10A5*(0.5*exp(-0.28*k-6.4/3)+exp(-(0.28+2/3)*k-8.8/3)/(1-exp(-2/3*k-0.8)-1);plot(k ,n)N 0.6078058E+130.000000VariableValue Reduced CostDual Price1.000000课程名称：数学建模班号 学号： 姓名：

30、 成绩:任课教师：郑勋烨日期：2011 年 5 月 29 日fun 1.m 文件：function n=fun 1(k)n=1.22*10A11*(1.109*10A5*(0.5*exp(-0.28*k-6.4/3)+exp(-(0.28+2/3)*k-8.8/3)/(1-exp(-2/3*k-0.8)-1)max1.m 文件:function y=max1( n,k)y=17.86*0.42*k/(0.8+0.42*k)*1.22*10A11/(1.22*10A11+ n)* n*exp(-1.6)*(1-exp(-2/3*(0.8+0.42*k)+22.99*k/(0.8+k)*1.22*

31、10A11/(1.22*10A11+ n)*n *exp(-0.28*k-2.4)/(1-exp(-2/3*k-0.8)*(1-exp(-2/3*(0.8+k)我们对k 1,20,主程序 1.m:for k=1:1:20n=fun 1(k);f=max1( n,k);kfend主程序 2.m :for k=17:0.1:18n=fun 1(k);f=max1( n,k);kfend主程序 3.m:for k=17.3:0.01:17.4n=fun 1(k);f=max1( n,k);kfend附录 3:fun.m 文件:课程名称：数学建模班号 学号： 姓名： 成绩:任课教师：郑勋烨日期：201

32、1 年 5 月 29 日fun ction y=fu n(k3,k,m,l1,l2,l3)m1= 10.1*10A9+29.7*10A9*exp(-0.8)+122*10A9*exp(-1.6)+l1*exp(-1.6)+l2*exp(-1.6);m21=3.29*10A9+10.1*10A9*exp(-(0.8+2/3*k3)+3.29*10A9*exp(-(0.8+2/3*k)+29.7*10A9*exp(-(1.6+2/3*k3)+l3*exp(-(0.8+2/3*k)+122*10A9*exp(-(2.4+2/3*k3)+l3*exp(-2*(0.8+2/3*k)+29.7*10A9*

33、exp(-(2.4+2/3*k3+2/3*k);m22=l1*exp(-(2.4+2/3*k3)+122*10A9*exp(-(3.2+2/3*k3+2/3*k)+29.7*10A9*exp(-(3.2+2/3*k3+4/3*k)+l3*exp(-3*(0.8+2/3*k);y=17.86*k3*(1-exp(-(0.8+k3)*2/3)*m1/(0.8+k3)+22.99*k*(1-exp(-(0.8+k)*2/3)*(m21+m22)/(0.8+k);main .mfor k=0:1:20m=1.109*10A5;k3=0.42*k;a=29.7*10A9*0.5*exp(-1/3*(4+

34、2*k3)+10.1*10A9*exp(-(4/3+2/3*k3+2/3*k)+3.29*10A9*exp(-1/3*(4+4/3*k);课程名称：数学建模班号 学号： 姓名： 成绩：l1=1.22*10A11*m*(10.1*10A9*0.5*exp(-2/3*(0.8+k3)+3.29*10A9*exp(-2/3*(0.8+k)(1.22*10M1+m*(10.1*10A9*0.5*exp(-2/3*(0.8+k3)+3.29*10A9*exp(-2/3*(0.8+k);l2=1.22*10A11*m*a/(1.22*10A11+m*a);l3=10.1*10A9*exp(-(0.8+2/3*k3)+3.29*10A9*exp(-(0.8+2/3*k);y=fu n( k3,k,m,l1,l2,l3);kyendmai n1.mfor k=17:0.1:18m=1.109*10A5;k3=0.42*k;a=29.7*10