1、 1数 学 模 型课程设计学号：姓名：2012 年 月 日2摘要本文研究了在渔业可持续捕捞的前提下，实现最优捕捞策略的问题。根据假设，鱼群的死亡和捕捞均视为连续的过程。建立微分方程，分析一年中不同时间点各龄鱼生长的变化，得到了各龄鱼的生长关于时间 t 的微分方程。 建立相关模型，得到目标函数和约束条件，并利用相关的 lingo 和 matlab 等数学软件对模型进行求解。针对问题一，根据题目中的要求，即每年开始捕捞时渔场中各年龄组鱼群不变为约束条件，以年最大捕捞量为目标函数，建立了非线性规划模型。通过 matlab 编程计算得出 ，此时对应的捕捞强度系数 k =17.3000 或 k =17.
2、4000。在模型的1f 3.870后续分析中，我得出了实现可持续捕捞的必要条件是 。 17.3,4针对问题二，可以考虑在 5 年的合同期内，捕捞强度不变，先确定一个收益最高的捕捞策略，再计算在这种策略下，5 年后各个鱼群的数量，进而确定是否对生产能力产生了太大破坏。但从理论分析可知，5 年后将在鱼群尽可能接近可持续鱼群的情况下来使捕捞量达到最大。对于破坏大小，由于各龄鱼群的数量在很大程度上受承包初 1龄鱼影响，我以第六年初 1 龄鱼群数量的变化量与承包时鱼群数量初值之比表示。不妨假定该比值小于 5%时，即可认为鱼群的生产能力没有受到太大破坏。由此求出捕捞强度在区间 内时，总捕捞量达到最大值 。
3、并且在模型的后续分17.58（ ， ） 12.605析中，证明了此捕捞强度下，鱼群生产能力并未遭到很大破坏，很容易就能恢复到最大稳定状态。关键词：最优捕捞策略 成活率 捕捞强度 3一问题重述为了保护人类赖以生存的自然环境，可再生资源（如渔业、林业资源）的开发必须适度。一种合理、简化的策略是，在实现可持续收获的前提下，追求最大产量或最佳效益。考虑对某种渔（鲳鱼）的最优捕捞策略：假设这种鱼分 4 个年龄组，称 1 龄鱼，4 龄鱼。各年龄组每条鱼的平均重量分别为 5.07,11.55,17.86,22.99(克)，各年龄组鱼的自然死亡率均 0.8(1/年)，这种鱼为季节性集中产卵繁殖，平均每条 4
4、龄鱼的产卵量为 (个)，3 龄鱼的产卵量为这51.09个数的一半，1 龄鱼和 2 龄鱼不产卵。产卵和孵化期为每年的最后 4 个月，卵孵化并成活为 1 龄鱼，成活率（龄鱼条数与产卵总量 n 之比） 。11.2/(.20n)渔业管理部门规定，每年只允许在产卵孵化期前的 8 个月进行捕捞作业。如果每年投入的捕捞能力（如渔船数、下网次数等）固定不变，这时单位时间捕捞量将与各年龄组鱼群条数成正比，比例系数不妨称捕捞强度系数。通常使用 13mm 网眼的拉网，这种网只能捕捞 3 龄鱼和 4 龄鱼，其两个捕捞强度系数之比为 0.42:1。渔业上称这种方式为固定努力量捕捞。1)建立数学模型分析如何实现可持续捕获
5、（即每年开始捕捞时渔场中各年龄组鱼群条数不变） ，并且在此前提下得到最高的年收获量（捕捞总重量） 。2)某渔业公司承包这种鱼的捕捞业务 年，合同要求 年后鱼群的生产能力不能受到太大破坏。已知承包时各年龄组鱼群的数量分别为： 。如果912,9.70,3.2(1)条仍用固定努力量的捕捞方式，该公司应采取怎样的策略才能使总收获量最高。二模型假设1、这种鱼在一年内的任何时间都会发生自然死亡，即死亡是一个连续的过程；2、捕捞也是一个连续的过程，不是在某一时刻突然发生；3、鱼群的死亡率已考虑种群的相互竞争及环境等因素；4、4 龄鱼产卵集中在 9 月初期,到第二年初完成孵化；45、 龄鱼到来年分别长一岁成
6、i+1 龄鱼，其中上一年存活下来的 4 龄鱼仍是 4 龄鱼。三符号定义-时间过程(年 )t-年份T-时间间隔t-i 龄鱼在时间 t 的数量（个）()iN-自然死亡率(0.8/年)r-年产卵总量n-4 龄鱼的产卵量( )，3 龄鱼为a 51.09/2a-卵的成活率b 12/(.2n)-4 龄鱼的捕捞强度( 龄鱼为 ) 决策变量 4k 3k： 4.k四问题分析4.1、题目简述及解释：1、题目简述鱼的基本情况:年龄组 1 龄鱼 2 龄鱼 3 龄鱼 4 龄鱼平均重量（克） 5.07 11.55 17.86 22.99产卵量（个） 0 0 a/2 aa）各年龄组鱼的自然死亡率均为 0.8(1年)b）产卵
7、和孵化期为每年最后 4 个月 c）捕鱼的限定：11.20=nn龄 鱼 条 数成 活 率 产 卵 总 量5a)渔业管理部门规定:只允许在产卵孵化期前的 8 个月进行捕捞作业 b)固定努力量捕捞： 投入的捕捞能力固定不变，这时单位时间捕捞量将与各年龄组鱼群条数成正比，比例系数不妨称捕捞强度系数。使用 13mm 网眼的拉网，只能捕捞 3 龄鱼和 4 龄鱼，其两个捕捞强度系数之比为 0.42:1 2、自然死亡率：由题目及假设给出的鱼的自然死亡率是指平均死亡率，即单位时间鱼群死亡数量与现有鱼群数量的比例系数。由死亡是一个数量的连续变化过程，同时1、2龄鱼在每年数量的变化及3、4龄鱼在后4个月的数量的变化
8、只与死亡率有关。由此可知，各龄鱼的变化满足：（1）()i1,234iidNtrt，3、捕捞强度系数：单位时间内捕捞量与各年龄组鱼群条数的比例系数，由于只能对 3、4 龄鱼进行捕捞，并且只在捕捞期内进行捕捞。所以，捕捞强度系数影响了 3、4 龄鱼在前八个月内内的数量变化：设 ，则4k30.42k故有（2）()(,34iidNtrxt4、卵的成活率：1、2 龄鱼不产卵，3、4 龄鱼在每年的后四个月产卵，假设了在 9 月初一次产卵，因此可将每年的产卵量 n 表示为：（3）5342.109.()N又有成活率 ，所以每年年初 1 龄鱼的数目为1.2bn=（4）11.200n4.2、问题一分析：6对于问题
9、一，商家的目标是得到最大的捕捞量，故由捕捞量得到目标函数。又要实现可持续捕获，即每年开始捕捞时渔场中各年龄组鱼群条数不变。可由每年年初鱼的数目不变得到约束条件。进而求最优化问题。对于捕捞量，只能对 3、4 龄鱼在 1 到 8 月进行捕捞。只求前八个月的对 3、4 龄鱼的捕捞量即为全年的捕捞量；对于鱼的数量变化，1、2 龄鱼不能被捕捞，仅受自然死亡率的影响；而 3、4 龄鱼不仅受自然死亡率的影响，还受捕捞强度系数的影响。 对于得到的最优化问题，可以用 lingo 软件和 matlab 软件进行求解。4.3、问题二分析：问题二中，合同要求 5 年后鱼群的生产能力不能受到太大破坏，又要使总收益最高，
10、这就有可能发生满足了前者满足不了后者之类的情况。对该问题，可以考虑在 5 年的合同期内，捕捞强度不变，先确定一个收益最高的捕捞策略，再计算在这种策略下，5 年后各个鱼群的数量，进而确定是否对生产能力产生了太大破坏。但从理论分析可知，5 年后将在鱼群尽可能接近可持续鱼群的情况下来使捕捞量达到最大。对于破坏大小，由于各龄鱼群的数量在很大程度上受承包初 1 龄鱼影响，我们以第六年初 1 龄鱼群数量的变化量与承包时鱼群数量初值之比表示。不妨假定该比值小于 5%时，即可认为鱼群的生产能力没有受到太大破坏。同理可按照问题一中的解法进行求解。五、模型的建立与求解5.1、问题一模型1可知 1、2 龄鱼的生长只
11、受自然死亡率的影响，由此可知 1、2 龄鱼的生长的微分方程满足方程（1）： ,12;iidNrt可得：（5）0(),;rtitei其中 每年年初龄鱼的数量。 0NT 年的 i 龄鱼在 T+1 年变为 i+1 龄鱼，则有7（6）1()()ri iNTe2而对于 3、4 龄鱼，在前八个月生长不仅受自然生长率和捕捞强度系数的影响，而后四个月仅受自然生长率的影响。以一年为一个时间单位，则这一时间单位可以分为前八个月和后四个月两个阶段。1）前八个月 3、4 龄鱼生长的微分方程满足：（7）()3,4iidNrkit可以得到（8）()0()irktite其中 为年初龄鱼的数量 0N由于每年的捕捞只在 1 到
12、 8 月进行，并且只能捕到 3、4 龄鱼，所以任意一个时刻的捕捞量为 ，则 i 龄鱼的年捕捞量为：()ikt（9）22 ()3300()1ikriikNtder2）后四个月 3、4 龄鱼生长的微分方程满足方程（1）： ,34iirdt可得：(10)0(),;rtiNtei 其中 每年年初龄鱼的数量 0其 中 产 卵 量(11) 22(0.4)()334.5rkrknaNeae又 孵 化 存 活 率(12)1.20bn所以年初 1 龄鱼的总量(13)()(NTb83根据以上分析，我们可以建立非线性规划模型：目标函数：(14)2233400max17.86().9()kNtdkNtd约束条件：22
13、(0.4)()3341223(/30.4)4.51krkrrkrnaeaeNbebn5.2、问题二模型针对渔业公司的 5 年捕捞计划，我们利用已得到的迭代方程在已知各个年龄组的鱼的初始值的前提下，可迭代求出各龄鱼群第 i 年的鱼量的分布的函数。整个生存过程满足的关系式为 3 334()2/()2/4112132(/)(2/3)434()/.0()0.krkrrkrkrntaNteatenttteNttNtek同时写出目标函数： 3 45 5()2/ ()2/33 43 41 1max()7.861()2.91(15)kr krt tNe ek 5.3、模型求解问题一求解我们可将目标函数和约束条
14、件转化为：目标函数为： 2 22(0.84) (0.8)(0.8)1.6 (0.28.4)3 330.42max17.86.91k kkkkbnebneee 约束条件：9128.(0.)6.43(0.28)1 (.)20(.511kkbnenae然后利用 lingo 软件和 matlab 软件分别进行求解。1）lingo 程序：（附录 1）直接运行,输出结果为：Local optimal solution found at iteration: 92Objective value: 0.3887076E+12Variable Value Reduced CostK 17.36292 -1.03
15、4723N 0.6078067E+13 0.000000Row Slack or Surplus Dual Price1 0.3887076E+12 1.0000002 0.000000 0.1258406E-02k17.3629,f0.876E+2）matlab 程序：（附录 2）首先,建立 fun1.m、max1.m、1.m、picture.m 文件, 运行 picture.m 文件,画出 n 关于 k 的图像：10由上面图像，可以看出 ，这与事实是相符的。,knA然后，在命令窗口中输入 1.m 文件中命令,运行出结果我们观察到 ， ,1,7kfA， ,因此，可以改变 k 的取值区间即步长
16、 ,建立了文件 2.m、3.m,从而求17,20kf得更精确的结果。从运行结果看：k =17.3000 或 k =17.4000，取最大值 1f 3.870最后，建立主程序 3.m （程序见附录 2）17.,4结果如下：，17.3,4k1f 3.80问题一所求结果为： 17,.4ff 3.870k时 ， 取 最 大 值3) 结果分析3龄鱼的捕捞强度为7.29/年；4龄鱼的捕捞强度为17.36/年：最优可持续捕捞量,可持续捕捞的鱼群大小（条数）：1f .8701 龄 2 龄 3 龄 4 龄数量（条） .960105.374102.4778.39610分析结果发现，4 龄鱼在年末存活的数量占全部数
17、量的比例相对很小。问题二求解将目标函数转化为： 1120.4 222322 20.4 0.43333. .030.4max1(1971 )(.9. 9.7rkrr r rrkrkrkrknnkeeee 2.4342230.4. .1.40.3 971901.2 rkrrkrr rne 24240. 30.4383 .1.) kkkrk ee11其中22(0.4) ()33(2)10./.9krkrnaeme5 545(.8)(0.8)()33971.9r rkraae这样 max 就变为关于 k 的函数,易于求解。用 matlab 软件求解1、利用 matlab 软件,建立 fun.m 文件(
18、见附录 3);2、首先，对 ，建立 main.m 主程序（见附录 3） 120k根据运行结果,观察得到 , , ,18kfA1820kfA其中 12178121978,.6053,.kff接着，改变 k 的取值区间及范围，建立主程序 main1.m,main2.m，求解出更精确的结果。运行结果分析：，17.5,8k12f .605问题二所求结果为： 127.,8ff .605k时 ， 取 最 大 值验证 5 年后鱼群的生产能力有没有受到太大破坏迭代求得第六年初各龄鱼群的数量为：11N(6).950102N(6)5.37324748第一年各龄鱼群的数量为11().0102().973N4N3第六
19、年 1 龄鱼数量占第一年 1 龄鱼数量的比例为： 1(6)0%98结果分析12捕捞强度在区间 内时，总捕捞量达到最大值 。在这种捕捞17.58（ ， ） 12.605强度下，5 年后 1 龄鱼数量占第一年 1 龄鱼数量的比例为 98%，即可认为生产能力没有受到太大破坏。因此，求解出的结果即为最优捕鱼策略。六模型的改进我们知道，原题中没有说明四龄以上的鱼如何处理。我们假设的是上一年存活下来的4龄鱼仍是4龄鱼，而事实上还可以假设这种鱼只活到4龄，以后它就死掉了。这对模型没有太大的差别，只是我们所做的假设的分析计算稍复杂，但计算结果也只是稍有差别，在我们模型的基础上，我们可以假设鱼只能活到4龄，这样
20、计算更简便一些。目标函数：2233400max17.86().9()kNtdkNtd约束条件： 22(0.4)()3341223(0.43)4.5krkrrkrnaeaebNebn我们可将上面的目标函数和约束条件转化为：目标函数为： 2 22(0.84) (0.43)(0.8)1.63 30.42max17.86.91.8k krkk kbnebnee约束条件： 16.428.(0.28)(0.)13320(.51kkbnnaee用lingo软件进行求解，算法见附录4：lingo2算法直接运行得：Local optimal solution found.Objective value: 0.3
21、887076E+12Extended solver steps: 5Total solver iterations: 20913Variable Value Reduced CostK 17.36293 0.000000N 0.6078058E+13 0.000000Row Slack or Surplus Dual Price1 0.3887076E+12 1.0000002 0.000000 0.1258410E-02即 117.3629,0.8760kf而以我们的假设算得出的结果为 12k7.392,f0.8760两个结果相差甚小，使改进后的模型计算非常简便。七模型的评价本文采用了非线性
22、规划的思想建立模型，通过求解有约束的非线性最大值问题，找到一组最优解。问题一，在实现可持续捕获(即每年开始捕捞时渔场中各年龄组鱼群条数不变)的前提下，用固定努力量的捕捞方式，确定捕捞策略以得到最大捕捞总质量。用微分、积分的方法来分析每年各龄鱼的数量，建立每年捕捞量的方程，用lingo软件与matlab软件分别求解，两个结果误差很小，肯定了结果的正确性。问题二，所求模型为五年鱼群生长模型的组合。由于所给的初始鱼群并不是可持续捕捞的鱼群，为了在五年内既得到最大的收益，又不破坏鱼群的生产能力，即五年14后在达到产量最高的条件下使得鱼群尽量接近可持续捕捞鱼群。我们在五年内以同样的强度实现固定努力量的捕
23、捞。对于每年每条龄鱼在每个时刻的条数，我们可以用算法迭代求解出n年的条数，从而比较第六年年初与初始时刻条数的差值，得出生产能力的破坏度不显著。本文有着一些优缺点，其优缺点如下：主要优点1. 推理过程严谨，理论性强。 2. 适用范围广，除鱼群最优捕捞问题外，同样适用于部分其他优化问题。 3. 稳定性好，模型本身不存在误差，误差主要来源软件计算误差。 主要缺点 1. 部分假设可能不符合实情。这些假设可能与实际情况不符，导致结果不准确。 2. 模型还有很多因素未纳入考虑范围，因此实用性方面可能有一定的缺陷。 3. 此问建立的模型并未考虑鱼场的经济效益问题（即未考虑最佳效益这一部分） 。而经济效益与很
24、多因素，包括生产成本、市场需求、市场价格等有关。因此仅仅考虑最优捕捞量是不足以得到最佳经济效益的。如果能搜集更多有关这方面的信息，就可以对模型进行修改，得出关于捕捞强度 k 与经济收益率的数学模型，以求出最佳经济效益。八参考文献【1】姜启源、谢金星、叶俊 数学模型 高等教育出版社 2003 年 8 月【2】刘来福 最优捕鱼策略问题答案评述 数学的实践与认识 1997 年 1 月【3】苏金明、阮沈用 mtlab 实用教程 电子工业出版社 2008 年 2 月15附录附录 1：lingo1 算法：max=17.86*0.42*k/(0.8+0.42*k)*1.22*1011/(1.22*1011+
25、n)*n*exp(-1.6)*(1-exp(-2/3*(0.8+0.42*k)+22.99*k/(0.8+k)*1.22*1011/(1.22*1011+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*1011*(1.109*105*(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);附录 2：用 matlab 求解，算法如下：picture1.mk=linspace(1,20,20);n=1.22*1011*(1.1
26、09*105*(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);16plot(k,n)fun1.m 文件:function n=fun1(k)n=1.22*1011*(1.109*105*(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*1011/(1.22*1011+n)*n*exp(-1
27、.6)*(1-exp(-2/3*(0.8+0.42*k)+22.99*k/(0.8+k)*1.22*1011/(1.22*1011+n)*n*exp(-0.28*k-2.4)/(1-exp(-2/3*k-0.8)*(1-exp(-2/3*(0.8+k)我们对 ，1,20k主程序 1.m：for k=1:1:20n=fun1(k);f=max1(n,k);kfend主程序 2.m：for k=17:0.1:18n=fun1(k);f=max1(n,k);k17fend主程序 3.m：for k=17.3:0.01:17.4n=fun1(k);f=max1(n,k);kfend附录 3：fun.m
28、 文件:function y=fun(k3,k,m,l1,l2,l3)m1=10.1*109+29.7*109*exp(-0.8)+122*109*exp(-1.6)+l1*exp(-1.6)+l2*exp(-1.6);m21=3.29*109+10.1*109*exp(-(0.8+2/3*k3)+3.29*109*exp(-(0.8+2/3*k)+29.7*109*exp(-(1.6+2/3*k3)+l3*exp(-(0.8+2/3*k)+122*109*exp(-(2.4+2/3*k3)+l3*exp(-2*(0.8+2/3*k)+29.7*109*exp(-(2.4+2/3*k3+2/3
29、*k);m22=l1*exp(-(2.4+2/3*k3)+122*109*exp(-(3.2+2/3*k3+2/3*k)+29.7*109*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*105;k3=0.42*k;a=29.7*109*0.5*exp(-181/3*(4+2*k3)+10.1*109*exp(-(4/
30、3+2/3*k3+2/3*k)+3.29*109*exp(-1/3*(4+4/3*k);l1=1.22*1011*m*(10.1*109*0.5*exp(-2/3*(0.8+k3)+3.29*109*exp(-2/3*(0.8+k)/(1.22*1011+m*(10.1*109*0.5*exp(-2/3*(0.8+k3)+3.29*109*exp(-2/3*(0.8+k);l2=1.22*1011*m*a/(1.22*1011+m*a);l3=10.1*109*exp(-(0.8+2/3*k3)+3.29*109*exp(-(0.8+2/3*k);y=fun(k3,k,m,l1,l2,l3);
31、kyendmain1.mfor k=17:0.1:18m=1.109*105;k3=0.42*k;a=29.7*109*0.5*exp(-1/3*(4+2*k3)+10.1*109*exp(-(4/3+2/3*k3+2/3*k)+3.29*109*exp(-1/3*(4+4/3*k);l1=1.22*1011*m*(10.1*109*0.5*exp(-2/3*(0.8+k3)+3.29*109*exp(-2/3*(0.8+k)/(1.22*1011+m*(10.1*109*0.5*exp(-2/3*(0.8+k3)+3.29*109*exp(-2/3*(0.8+k);l2=1.22*1011*
32、m*a/(1.22*1011+m*a);l3=10.1*109*exp(-(0.8+2/3*k3)+3.29*109*exp(-(0.8+2/3*k);y=fun(k3,k,m,l1,l2,l3);kyendmain2.m19for k=17.5:0.01:17.8m=1.109*105;k3=0.42*k;a=29.7*109*0.5*exp(-1/3*(4+2*k3)+10.1*109*exp(-(4/3+2/3*k3+2/3*k)+3.29*109*exp(-1/3*(4+4/3*k);l1=1.22*1011*m*(10.1*109*0.5*exp(-2/3*(0.8+k3)+3.29
33、*109*exp(-2/3*(0.8+k)/(1.22*1011+m*(10.1*109*0.5*exp(-2/3*(0.8+k3)+3.29*109*exp(-2/3*(0.8+k);l2=1.22*1011*m*a/(1.22*1011+m*a);l3=10.1*109*exp(-(0.8+2/3*k3)+3.29*109*exp(-(0.8+2/3*k);y=fun(k3,k,m,l1,l2,l3);kyend附录 4：lingo2 算法：max=17.86*0.42*k/(0.8+0.42*k)*1.22*1011/(1.22*1011+n)*n*exp(-1.6)*(1-exp(-2/3*(0.8+0.42*k)+22.99*k/(0.8+k)*1.22*1011/(1.22*1011+n)*n*exp(-0.28*k-2.4) *(1-exp(-2/3*(0.8+k);n=1.22*1011*(1.109*105*(0.5*exp(-0.28*k-6.4/3)+exp(-(0.28+2/3)*k-8.8/3)-1);
