数学模型课程设计.doc

1、 摘要本文研究了在渔业可持续捕捞的前提下，实现最优捕捞策略的问题。根据假设，鱼群的死亡和捕捞均视为连续的过程。建立微分方程，分析一年中不同时间点各龄鱼生长的变化，得到了各龄鱼的生长关于时间t的微分方程。 建立相关模型，得到目标函数和约束条件，并利用相关的lingo和matlab等数学软件对模型进行求解。针对问题一，根据题目中的要求，即每年开始捕捞时渔场中各年龄组鱼群不变为约束条件，以年最大捕捞量为目标函数，建立了非线性规划模型。通过matlab编程计算得出，此时对应的捕捞强度系数k =17.3000或k =17.4000。在模型的后续分析中，我得出了实现可持续捕捞的必要条件是。 针对问题二，可

2、以考虑在5年的合同期内，捕捞强度不变，先确定一个收益最高的捕捞策略，再计算在这种策略下，5年后各个鱼群的数量，进而确定是否对生产能力产生了太大破坏。但从理论分析可知，5年后将在鱼群尽可能接近可持续鱼群的情况下来使捕捞量达到最大。对于破坏大小，由于各龄鱼群的数量在很大程度上受承包初1龄鱼影响，我以第六年初1龄鱼群数量的变化量与承包时鱼群数量初值之比表示。不妨假定该比值小于5%时，即可认为鱼群的生产能力没有受到太大破坏。由此求出捕捞强度在区间内时，总捕捞量达到最大值。并且在模型的后续分析中，证明了此捕捞强度下，鱼群生产能力并未遭到很大破坏，很容易就能恢复到最大稳定状态。关键词：最优捕捞策略 成活率

3、 捕捞强度 一问题重述为了保护人类赖以生存的自然环境，可再生资源（如渔业、林业资源）的开发必须适度。一种合理、简化的策略是，在实现可持续收获的前提下，追求最大产量或最佳效益。考虑对某种渔（鲳鱼）的最优捕捞策略：假设这种鱼分4个年龄组，称1龄鱼，4龄鱼。各年龄组每条鱼的平均重量分别为5.07,11.55,17.86,22.99(克)，各年龄组鱼的自然死亡率均0.8(1/年)，这种鱼为季节性集中产卵繁殖，平均每条4龄鱼的产卵量为(个)，3龄鱼的产卵量为这个数的一半，1龄鱼和2龄鱼不产卵。产卵和孵化期为每年的最后4个月，卵孵化并成活为1龄鱼，成活率（龄鱼条数与产卵总量n之比）。渔业管理部门规定，每年

4、只允许在产卵孵化期前的8个月进行捕捞作业。如果每年投入的捕捞能力（如渔船数、下网次数等）固定不变，这时单位时间捕捞量将与各年龄组鱼群条数成正比，比例系数不妨称捕捞强度系数。通常使用13mm网眼的拉网，这种网只能捕捞3龄鱼和4龄鱼，其两个捕捞强度系数之比为0.42:1。渔业上称这种方式为固定努力量捕捞。1)建立数学模型分析如何实现可持续捕获（即每年开始捕捞时渔场中各年龄组鱼群条数不变），并且在此前提下得到最高的年收获量（捕捞总重量）。2)某渔业公司承包这种鱼的捕捞业务年，合同要求年后鱼群的生产能力不能受到太大破坏。已知承包时各年龄组鱼群的数量分别为： 。如果仍用固定努力量的捕捞方式，该公司应采取

5、怎样的策略才能使总收获量最高。二模型假设1、这种鱼在一年内的任何时间都会发生自然死亡，即死亡是一个连续的过程；2、捕捞也是一个连续的过程，不是在某一时刻突然发生；3、鱼群的死亡率已考虑种群的相互竞争及环境等因素；4、4龄鱼产卵集中在9月初期,到第二年初完成孵化；5、龄鱼到来年分别长一岁成i+1龄鱼，其中上一年存活下来的4龄鱼仍是4龄鱼。三符号定义 -时间过程(年) -年份 -时间间隔 -i龄鱼在时间t的数量（个） -自然死亡率(0.8/年) -年产卵总量 -4龄鱼的产卵量()，3龄鱼为 -卵的成活率 -4龄鱼的捕捞强度(龄鱼为) 决策变量 四问题分析4.1、题目简述及解释：1、题目简述鱼的基本

6、情况:年龄组1龄鱼2龄鱼3龄鱼4龄鱼 平均重量（克）5.0711.5517.8622.99产卵量（个）00a/2aa）各年龄组鱼的自然死亡率均为0.8(1年)b）产卵和孵化期为每年最后4个月 c）捕鱼的限定：a)渔业管理部门规定: 只允许在产卵孵化期前的8个月进行捕捞作业 b)固定努力量捕捞： 投入的捕捞能力固定不变，这时单位时间捕捞量将与各年龄组鱼群条数成正比，比例系数不妨称捕捞强度系数。使用13mm网眼的拉网，只能捕捞3龄鱼和4龄鱼，其两个捕捞强度系数之比为0.42:1 2、 自然死亡率： 由题目及假设给出的鱼的自然死亡率是指平均死亡率，即单位时间鱼群死亡数量与现有鱼群数量的比例系数。由死

7、亡是一个数量的连续变化过程，同时1、2龄鱼在每年数量的变化及3、4龄鱼在后4个月的数量的变化只与死亡率有关。由此可知，各龄鱼的变化满足： （1）3、 捕捞强度系数：单位时间内捕捞量与各年龄组鱼群条数的比例系数，由于只能对3、4龄鱼进行捕捞，并且只在捕捞期内进行捕捞。所以，捕捞强度系数影响了3、4龄鱼在前八个月内内的数量变化：设，则故有 （2）4、 卵的成活率：1、2龄鱼不产卵，3、4龄鱼在每年的后四个月产卵，假设了在9月初一次产卵，因此可将每年的产卵量n表示为： （3）又有成活率，所以每年年初1龄鱼的数目为 （4）4.2、问题一分析：对于问题一，商家的目标是得到最大的捕捞量，故由捕捞量得到目标

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

9、捞强度不变，先确定一个收益最高的捕捞策略，再计算在这种策略下，5年后各个鱼群的数量，进而确定是否对生产能力产生了太大破坏。但从理论分析可知，5年后将在鱼群尽可能接近可持续鱼群的情况下来使捕捞量达到最大。对于破坏大小，由于各龄鱼群的数量在很大程度上受承包初1龄鱼影响，我们以第六年初1龄鱼群数量的变化量与承包时鱼群数量初值之比表示。不妨假定该比值小于5%时，即可认为鱼群的生产能力没有受到太大破坏。同理可按照问题一中的解法进行求解。五、模型的建立与求解5.1、问题一模型1可知1、2龄鱼的生长只受自然死亡率的影响，由此可知1、2龄鱼的生长的微分方程满足方程（1）：可得： （5）其中每年年初龄鱼的数量。

10、 T年的i龄鱼在T+1年变为i+1龄鱼，则有 （6）2而对于3、4龄鱼，在前八个月生长不仅受自然生长率和捕捞强度系数的影响，而后四个月仅受自然生长率的影响。以一年为一个时间单位，则这一时间单位可以分为前八个月和后四个月两个阶段。1） 前八个月3、4龄鱼生长的微分方程满足： （7）可以得到 （8）其中为年初龄鱼的数量 由于每年的捕捞只在1到8月进行，并且只能捕到3、4龄鱼，所以任意一个时刻的捕捞量为，则i龄鱼的年捕捞量为： （9）2） 后四个月3、4龄鱼生长的微分方程满足方程（1）：可得： (10)其中每年年初龄鱼的数量 (11) (12)所以年初1龄鱼的总量 (13)3根据以上分析，我们可以建

11、立非线性规划模型：目标函数： (14)约束条件： 5.2、问题二模型针对渔业公司的5年捕捞计划，我们利用已得到的迭代方程在已知各个年龄组的鱼的初始值的前提下，可迭代求出各龄鱼群第i年的鱼量的分布的函数。整个生存过程满足的关系式为同时写出目标函数：5.3、模型求解问题一求解我们可将目标函数和约束条件转化为：目标函数为：约束条件：然后利用lingo软件和matlab软件分别进行求解。1）lingo 程序：（附录1）直接运行,输出结果为：Local optimal solution found at iteration: 92 Objective value: 0.3887076E+12 Varia

12、ble Value Reduced Cost K 17.36292 -1.034723 N 0.6078067E+13 0.000000 Row Slack or Surplus Dual Price 1 0.3887076E+12 1.000000 2 0.000000 0.1258406E-022）matlab程序：（附录2）首先,建立fun1.m、max1.m、1.m、picture.m文件,运行picture.m文件,画出n关于k的图像：由上面图像，可以看出 ，这与事实是相符的。然后，在命令窗口中输入1.m文件中命令,运行出结果我们观察到，,，,因此，可以改变k的取值区间即步长,建立了

13、文件2.m、3.m,从而求得更精确的结果。从运行结果看：k =17.3000或k =17.4000，取最大值最后，建立主程序3.m （程序见附录2）结果如下：，问题一所求结果为：3) 结果分析3龄鱼的捕捞强度为7.29/年；4龄鱼的捕捞强度为17.36/年：最优可持续捕捞量,可持续捕捞的鱼群大小（条数）：1龄2龄3龄4龄数量（条）分析结果发现，4龄鱼在年末存活的数量占全部数量的比例相对很小。问题二求解将目标函数转化为：其中这样max就变为关于k的函数,易于求解。用matlab软件求解1、利用matlab软件,建立fun.m文件(见附录3);2、首先，对，建立main.m主程序（见附录3） 根据

14、运行结果,观察得到,其中接着，改变k的取值区间及范围，建立主程序main1.m,main2.m，求解出更精确的结果。运行结果分析：，问题二所求结果为：验证5年后鱼群的生产能力有没有受到太大破坏迭代求得第六年初各龄鱼群的数量为： 第一年各龄鱼群的数量为 第六年1龄鱼数量占第一年1龄鱼数量的比例为：结果分析捕捞强度在区间内时，总捕捞量达到最大值。在这种捕捞强度下，5年后1龄鱼数量占第一年1龄鱼数量的比例为98%，即可认为生产能力没有受到太大破坏。因此，求解出的结果即为最优捕鱼策略。六模型的改进我们知道，原题中没有说明四龄以上的鱼如何处理。我们假设的是上一年存活下来的4龄鱼仍是4龄鱼，而事实上还可以

15、假设这种鱼只活到4龄，以后它就死掉了。这对模型没有太大的差别，只是我们所做的假设的分析计算稍复杂，但计算结果也只是稍有差别，在我们模型的基础上，我们可以假设鱼只能活到4龄，这样计算更简便一些。目标函数：约束条件：我们可将上面的目标函数和约束条件转化为：目标函数为：约束条件：用lingo软件进行求解，算法见附录4：lingo2算法直接运行得：Local optimal solution found. Objective value: 0.3887076E+12 Extended solver steps: 5 Total solver iterations: 209 Variable Value

16、 Reduced Cost K 17.36293 0.000000 N 0.6078058E+13 0.000000 Row Slack or Surplus Dual Price 1 0.3887076E+12 1.000000 2 0.000000 0.1258410E-02即而以我们的假设算得出的结果为两个结果相差甚小，使改进后的模型计算非常简便。七模型的评价本文采用了非线性规划的思想建立模型，通过求解有约束的非线性最大值问题，找到一组最优解。问题一，在实现可持续捕获(即每年开始捕捞时渔场中各年龄组鱼群条数不变)的前提下，用固定努力量的捕捞方式，确定捕捞策略以得到最大捕捞总质量。用微分、

17、积分的方法来分析每年各龄鱼的数量，建立每年捕捞量的方程，用lingo软件与matlab软件分别求解，两个结果误差很小，肯定了结果的正确性。问题二，所求模型为五年鱼群生长模型的组合。由于所给的初始鱼群并不是可持续捕捞的鱼群，为了在五年内既得到最大的收益，又不破坏鱼群的生产能力，即五年后在达到产量最高的条件下使得鱼群尽量接近可持续捕捞鱼群。我们在五年内以同样的强度实现固定努力量的捕捞。对于每年每条龄鱼在每个时刻的条数，我们可以用算法迭代求解出n年的条数，从而比较第六年年初与初始时刻条数的差值，得出生产能力的破坏度不显著。本文有着一些优缺点，其优缺点如下：主要优点1. 推理过程严谨，理论性强。 2.

18、 适用范围广，除鱼群最优捕捞问题外，同样适用于部分其他优化问题。 3. 稳定性好，模型本身不存在误差，误差主要来源软件计算误差。 主要缺点 1. 部分假设可能不符合实情。这些假设可能与实际情况不符，导致结果不准确。 2. 模型还有很多因素未纳入考虑范围，因此实用性方面可能有一定的缺陷。 3. 此问建立的模型并未考虑鱼场的经济效益问题（即未考虑最佳效益这一部分）。而经济效益与很多因素，包括生产成本、市场需求、市场价格等有关。因此仅仅考虑最优捕捞量是不足以得到最佳经济效益的。如果能搜集更多有关这方面的信息，就可以对模型进行修改，得出关于捕捞强度k与经济收益率的数学模型，以求出最佳经济效益。八参考文

19、献【1】姜启源、谢金星、叶俊 数学模型 高等教育出版社 2003年8月【2】刘来福 最优捕鱼策略问题答案评述 数学的实践与认识 1997年1月【3】苏金明、阮沈用 mtlab实用教程 电子工业出版社 2008年2月附录附录1：lingo1算法：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*k-0.8)*(1-exp(-2

20、/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.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);plot(k,n)fun1.m文件:function n=fun1(k)n=1.22*1011*(1.1

21、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)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.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

22、.m：for k=1:1:20 n=fun1(k); f=max1(n,k); k fend主程序2.m：for k=17:0.1:18 n=fun1(k); f=max1(n,k); k fend主程序3.m：for k=17.3:0.01:17.4 n=fun1(k); f=max1(n,k); k fend附录3：fun.m文件: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

23、(-(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*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(-

24、(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:20 m=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

25、.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); k yendmain1.mfor k=17:0.1:18 m=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*e

26、xp(-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); k yendmain2.mfor k=17

27、.5:0.01:17.8 m=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*m*a/(1.22*

28、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); k yend附录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); .18