最优捕鱼策略.docx
《最优捕鱼策略.docx》由会员分享,可在线阅读,更多相关《最优捕鱼策略.docx(15页珍藏版)》请在冰豆网上搜索。
最优捕鱼策略
最优捕鱼策略
摘要
为了保护人类赖以生存的自然环境,可再生资源(如渔业、林业资源等)的开发必须适度。
而在社会经济生活中,我们要使商业活动在一段时期内达到最大收益,因此我们要合理的开发资源,这时,我们不仅要考虑商业活动的当前经济效益,还要考虑生态效益及由此产生的对整体经济效益的影响。
本文就是对渔业这类可再生资源的开发问题进行研究,利用相关的数学软件进行求解。
对于问题一,我们考虑渔场生产过程中的各年龄组鱼群数量的制约因素,将其分为两大类,第1,2龄鱼群为一类,该鱼群数量变化在一年内只受自然死亡率制约,写出鱼群数量满足的微分方程;第3,4龄鱼群为一类,其数量变化在前8个月受捕捞强度和自然死亡率影响,后4个月只受自然死亡率的制约,分阶段写出写出鱼群数量满足的微分方程;根据微分方程,求出在某时刻各鱼群的数量表达式(类似于人口增长模型)。
因为捕捞是连续的,所以任意一个时刻的捕捞量为捕捞强度乘以鱼群的数量,又捕捞只在前8个月进行,则年捕捞量为前8个月各时刻鱼群数量的积分。
最后建立年总捕捞量的函数与生产过程中满足的关系式,转化为非线性规划模型,利用lingo和matlab软件分别求解。
对于问题二,题中已给出各年龄组鱼群的初始值,我们利用问题一中所得到的迭代方程,可迭代地求出第i年初各年龄组鱼群的数量;再根据问题一中的捕捞量表达式,可写出5年的捕捞总量表达式,以5年捕捞总量最大为前提,利用matlab软件求解出此时的捕捞强度,然后再验证在此捕捞强度下会不会使5年后鱼群的生产能力有太大的破坏。
最后,我们得出以下结论:
可持续捕获条件下,捕捞强度为17.36292时,达到最大捕捞总质量
g;5年后鱼群的生产能力不会有太大的破坏条件下,捕捞强度为
,达到最大最大捕捞总质量
关键词:
渔业最大收益捕捞策略生产能力生长率lingomatlab
一.问题重述
生态学表明,对可再生资源的开发策略应在事先可持续收获的前提下追求最大经济效益,考虑具有4个年龄组:
1龄鱼,……,4龄鱼的某种鱼。
该鱼类在每年后4个月季节性集中产卵繁殖。
而按规定,捕捞作业只允许在前8个月进行,每年投入的捕捞能力固定不变,单位时间捕捞量与各年龄组鱼群条数的比例称为捕捞强度系数。
使用只能捕捞3,4龄鱼的13mm网眼的拉网,其两个捕捞强度系数比为0.42:
1。
渔业上称这种方式为固定努力量捕捞。
该鱼群本身有如下数据:
各年龄组鱼的自然死亡率为0.8(1/年),其平均质量分别为5.07,11.55,17.86,22.99(单位:
g);1,2龄鱼不产卵,平均每条4龄鱼产卵量为
(个),3龄鱼为其一半;卵孵化的成活率为
(n为产卵总量);
有如下问题需要解决:
1.1.问题一就是在实现可持续捕获(即每年开始捕捞时渔场中各年龄组鱼群条数不变)的前提下,用固定努力量的捕捞方式,确定捕捞策略以得到最大捕捞总质量。
1.2.问题二就是给出了承包时各年龄组鱼群的数量,要求5年后鱼群的生产能力不会有太大的破坏,在用固定努力量的捕捞方式的前提下,确定捕捞策略,求出最大捕捞总质量。
综上所述,原问题实质上是给出了各年龄组鱼群之间数量的变化规律,并给出了它们的自然死亡率及捕捞和产卵的时间分布,并固定3、4龄鱼捕捞能力的比值,要求选择一定的捕捞能力系数,使得各年龄组鱼的数量在各年开始的第一天条数不变(第一问),5年后鱼群的生产能力不会有太大的破坏(第二问),并在此条件下,求到最大捕获量。
二.符号说明
三.模型假设
1.这种鱼在一年内的任何时间都会发生自然死亡,即死亡是一个连续的过程。
2.捕捞也是一个连续的过程,不是在某一时刻突然发生。
3.1、2龄鱼体形太小,不能被捕。
4.3、4龄鱼在一年中的后4个月的第一天集中一次产卵
5.i龄鱼到来年分别长一岁成为i+1龄鱼,i=1,2,3,其中上一年存活下来的4龄鱼仍是4龄鱼
四.模型的建立与求解
4.1.问题分析
4.1.1.对题中一些术语的解释:
对自然死亡率的理解:
本题中给出的鱼的自然死亡率是指平均死亡率,即单位时间鱼群死亡数量与现有鱼群数量的比例系数,它与环境等其它外在因素无关;这是一个有量纲的量,它既不是简单的百分率又不是简单的变化速率,实际上它是百分比率的变化率。
它应该理解为以每年死亡80%的速率减少,并不是在一年内恰好死亡80%。
另一方面,鱼群的数量是连续变化的,且1,2龄鱼在全年及3,4龄鱼在后4个月的数量只与死亡率有关。
由此可知,各龄鱼的变化满足:
对捕捞强度系数的理解:
捕捞强度系数是单位时间内捕捞量与各年龄组鱼群条数的比例系数,单位时间4龄鱼捕捞量与4龄鱼群总数成正比,捕捞强度系数是一定的,且只在捕捞期内(即每年的前8个月)捕捞3,4龄鱼。
所以,捕捞强度系数k影响了3,4龄鱼在捕捞期内的数量变化:
则有
对卵的成活率的理解:
1,2龄鱼不产卵,3,4龄鱼在每年的后四个月产卵,我们假设了在9月初一次产卵,因此可将每年的产卵量n表示为:
4.1.2.问题一分析:
对于问题一,要实现可持续捕获,即每年开始捕捞时渔场中各年龄组鱼群条数不变,因此我们要算出每年初各龄鱼组的数量。
4.1.1中已对自然死亡率,捕捞强度系数和卵成活率作出了解释,即1,2龄鱼仅受自然死亡率的影响;而3,4龄鱼不仅受自然死亡率的影响,还受捕捞强度系数的影响;因为该种鱼的最高寿命为4,所以在后四月中4龄鱼都不存活;,而对于1龄鱼的数量,是3,4龄鱼在前年的后4年产卵所存活下来的数量;对于捕捞量,题中规定只在1到8月才能捕捞,而且1,2龄鱼不被捕捞,所以主要来源于对3,4龄鱼的捕捞。
根据这些关系可列出一系列的方程,其中捕捞量作为目标函数,其他的作为约束条件,建立一个非线性规划模型,再然后用lingo软件和matlab软件进行求解。
4.1.3.问题二分析:
对于问题二,合同要求5年后鱼群的生产能力不能受到太大破坏,又要使总收益最高,这就有可能发生满足了前者满足不了后者之类的情况。
我们处理方法是先确定一个策略使其收益最高,再检验此捕鱼策略是否能保证5年后鱼群的生产能力不受到太大的破坏,若它让鱼群的生产能力受到了严重破坏,我们再求另外一种策略。
但从理论分析可知,5年后将在鱼群尽可能接近可持续鱼群的情况下来使捕捞量达到最大。
对于破坏大小,我们采用1龄鱼群数量变化率来衡量,即以第六年初1龄鱼群数量的变化量与承包时鱼群数量初值之比表示,因为2,3,4龄鱼群的数量在很大程度上受承包初1龄鱼影响,根据关系,可以知道5年后2,3,4龄鱼群的数量肯定会有较大变化。
只要该比值小于5%,我们就认为鱼群的生产能力没有受到太大破坏。
题中已经给了我们各年龄组的初始值,而问题一中也已得出一组迭代方程,我们利用这些迭代方程,求出各年的鱼量分布;同样可以根据问题一中捕捞量的表达式求出5年的总捕捞量,以此来确定我们的最优捕捞策略。
再然后我们通过验证来确定其5年后鱼群的生产能力有没有受到太大破坏。
4.2.模型建立
4.2.1问题一模型
1.由4.1.1中对自然死亡率的理解中的
(1)式,可知1,2龄鱼的生长只受自然死亡率的影响,由此可知1,2龄鱼的生长的微分方程满足方程
(1):
T年的i龄鱼在T+1年变为i+1龄鱼,
2.而对于3,4龄鱼的生长,在前八个月,他们的生长不仅受自然生长率的影响,还受捕捞强度系数的影响,而后四个月仅受自然生长率的影响。
我们以一年为一个时间单位,则这一时间单位可以分为两个阶段,见图
(1):
02/31
I:
捕捞期II:
产卵期
图
(1)
因此,
.前八个月3、4龄鱼生长的微分方程满足:
由于每年的捕捞只在1到8月进行,并且只能捕到3,4龄鱼,所以任意一个时刻的捕捞量为
,则年捕捞量为:
.后四个月3、4龄鱼生长的微分方程满足方程
(1):
所以年初1龄鱼的总量
3.根据以上分析,我们可以建立非线性规划模型:
目标函数:
约束条件:
4.2.2问题二模型
针对渔业公司的5年捕捞计划,我们利用已得到的迭代方程在已知各个年龄组的鱼的初始值的前提下,可迭代求出各龄鱼群第i年的鱼量的分布的函数。
整个生存过程满足的关系式为
同时写出目标函数:
4.3.模型求解
4.3.1问题一求解
4.3.1.1由4.2.1中的3,我们可将目标函数和约束条件转化为:
目标函数为:
约束条件:
4.3.1.2然后我们利用lingo软件和matlab软件分别进行求解。
1)lingo程序:
(附录1)
直接运行,输出结果为:
Localoptimalsolutionfoundatiteration:
92
Objectivevalue:
0.3887076E+12
VariableValueReducedCost
K17.36292-1.034723
N0.6078067E+130.000000
RowSlackorSurplusDualPrice
10.3887076E+121.000000
20.0000000.1258406E-02
2)matlab程序:
(附录2)
首先,建立fun1.m、max1.m、1.m、picture.m文件,运行picture.m文件,画出n关于k的图像:
由图像,,我们看出
,这与事实是相符的。
然后,在命令窗口中输入1.m文件中命令,运行出结果我们观察到
,
,
,
,因此,我们改变k的取值区间即步长,建立了文件2.m、3.m,从而求得更精确的结果。
从运行结果看:
k=17.3000或k=17.4000,取最大值f=3.8871e+011
最后,建立主程序3.m
(程序见附录2)
从运行结果看:
,f=3.8871e+011
问题一所求结果为:
4.3.1.3结果分析
3龄鱼的捕捞强度为7.29/年;4龄鱼的捕捞强度为17.36/年:
最优可持续捕捞量
可持续捕捞的鱼群大小(条数):
1龄
2龄
3龄
4龄
分析结果发现,4龄鱼在年末存活的数量占全部数量的比例相对很小。
4.3.2问题二求解
4.3.2.1将目标函数转化为:
其中
这样max就变为关于k的函数,易于求解。
4.3.2.2用matlab软件求解
1、利用matlab软件,建立fun.m文件(见附录3);
2、首先,我们对
,建立main.m主程序(见附录3)
根据运行结果,我们观察到
,
,
,
其中,
接着,改变k的取值区间及范围,建立主程序main1.m,main2.m,求解出更精确的结果。
运行结果分析:
,f=1.6056e+012
问题二所求结果为:
4.3.2.3、验证5年后鱼群的生产能力有没有受到太大破坏
迭代求得第六年初各龄鱼群的数量为:
第一年各龄鱼群的数量为
第六年1龄鱼数量占第一年1龄鱼数量的比例为:
4.3.2.4、结果分析
捕捞强度在区间(17.5,17.8)内时(因为电脑精确度问题,暂时只能精确到这一区间),总捕捞量达到最大值
。
在这种捕捞强度下,5年后1龄鱼数量占第一年1龄鱼数量的比例为98%,即可认为生产能力没有受到太大破坏。
因此,求解出的结果即为最优捕鱼策略。
五.模型的评价
我们采用了非线性规划的思想建立模型,通过求解有约束的非线性最大值问题,找到一组最优解。
问题一,在实现可持续捕获(即每年开始捕捞时渔场中各年龄组鱼群条数不变)的前提下,用固定努力量的捕捞方式,确定捕捞策略以得到最大捕捞总质量。
我们结合人口增长,地中海鲨鱼模型,用微分、积分的方法来分析每年各龄鱼的数量,建立每年捕捞量的方程,用lingo软件与matlab软件分别求解,两个结果误差很小,肯定了结果的正确性。
问题二,所求模型为五年(五组类似问题一的模型)鱼群生长模型的组合。
由于所给的初始鱼群并不是可持续捕捞的鱼群,为了在五年内既得到最大的收益,又不破坏鱼群的生产能力,即五年后在达到产量最高的条件下使得鱼群尽量接近可持续捕捞鱼群。
我们在五年内以同样的强度实现固定努力量的捕捞。
对于每年每条龄鱼在每个时刻的条数,我们可以用算法迭代求解出n年的条数,从而比较第六年年初与初始时刻条数的差值,得出生产能力的破坏度不显着。
本模型采用连续模型的方法,成功地解决了可持续捕捞问题,得到了较为精确且合理的结果。
六.模型的改进
我们知道,原题中没有说明四龄以上的鱼如何处理。
我们假设的是上一年存活下来的4龄鱼仍是4龄鱼,而事实上还可以假设这种鱼只活到4龄,以后它就死掉了。
这对模型没有太大的差别,只是我们所做的假设的分析计算稍复杂,但计算结果也只是稍有差别,在我们模型的基础上,我们可以假设鱼只能活到4龄,这样计算更简便一些。
目标函数:
约束条件:
我们可将上面的目标函数和约束条件转化为:
目标函数为:
约束条件:
用lingo软件进行求解,算法见附录4:
lingo2算法
直接运行得:
Localoptimalsolutionfound.
Objectivevalue:
0.3887076E+12
Extendedsolversteps:
5
Totalsolveriterations:
209
VariableValueReducedCost
K17.362930.000000
N0.6078058E+130.000000
RowSlackorSurplusDualPrice
10.3887076E+121.000000
20.0000000.1258410E-02
即
而以我们的假设算得出的结果为
两个结果相差甚小,但改进的模型计算非常简便。
七.参考文献
赵静,但琦,《数学建模与数学实验(第3版)》高等教育出版社
姜启源、谢金星、叶俊,《数学模型(第三版)》高等教育出版社,2003
刘来福,《最优捕鱼策略问题答案评述》,数学的实践与认识,1997年1月
八.附录
附录1:
lingo1算法:
max=17.86*0.42*k/(0.8+0.42*k)*1.22*10^11/(1.22*10^11+n)*n*@exp(-1.6)*(1-@exp(-2/3*(0.8+0.42*k)))+22.99*k/(0.8+k)*1.22*10^11/(1.22*10^11+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*10^11*(1.109*10^5*(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.m
k=linspace(1,20,20);
n=1.22*10^11*(1.109*10^5*(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文件:
functionn=fun1(k)
n=1.22*10^11*(1.109*10^5*(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文件:
functiony=max1(n,k)
y=17.86*0.42*k/(0.8+0.42*k)*1.22*10^11/(1.22*10^11+n)*n*exp(-1.6)*(1-exp(-2/3*(0.8+0.42*k)))+22.99*k/(0.8+k)*1.22*10^11/(1.22*10^11+n)*n*exp(-0.28*k-2.4)/(1-exp(-2/3*k-0.8))*(1-exp(-2/3*(0.8+k)))我们对
,
主程序1.m:
fork=1:
1:
20
n=fun1(k);
f=max1(n,k);
k
f
end
主程序2.m:
fork=17:
0.1:
18
n=fun1(k);
f=max1(n,k);
k
f
end
主程序3.m:
fork=17.3:
0.01:
17.4
n=fun1(k);
f=max1(n,k);
k
f
end
附录3:
fun.m文件:
functiony=fun(k3,k,m,l1,l2,l3)
m1=10.1*10^9+29.7*10^9*exp(-0.8)+122*10^9*exp(-1.6)+l1*exp(-1.6)+l2*exp(-1.6);
m21=3.29*10^9+10.1*10^9*exp(-(0.8+2/3*k3))+3.29*10^9*exp(-(0.8+2/3*k))+29.7*10^9*exp(-(1.6+2/3*k3))+l3*exp(-(0.8+2/3*k))+122*10^9*exp(-(2.4+2/3*k3))+l3*exp(-2*(0.8+2/3*k))+29.7*10^9*exp(-(2.4+2/3*k3+2/3*k));
m22=l1*exp(-(2.4+2/3*k3))+122*10^9*exp(-(3.2+2/3*k3+2/3*k))+29.7*10^9*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.m
fork=0:
1:
20
m=1.109*10^5;
k3=0.42*k;
a=29.7*10^9*0.5*exp(-1/3*(4+2*k3))+10.1*10^9*exp(-(4/3+2/3*k3+2/3*k))+3.29*10^9*exp(-1/3*(4+4/3*k));
l1=1.22*10^11*m*(10.1*10^9*0.5*exp(-2/3*(0.8+k3))+3.29*10^9*exp(-2/3*(0.8+k)))/(1.22*10^11+m*(10.1*10^9*0.5*exp(-2/3*(0.8+k3))+3.29*10^9*exp(-2/3*(0.8+k))));
l2=1.22*10^11*m*a/(1.22*10^11+m*a);
l3=10.1*10^9*exp(-(0.8+2/3*k3))+3.29*10^9*exp(-(0.8+2/3*k));
y=fun(k3,k,m,l1,l2,l3);
k
y
end
main1.m
fork=17:
0.1:
18
m=1.109*10^5;
k3=0.42*k;
a=29.7*10^9*0.5*exp(-1/3*(4+2*k3))+10.1*10^9*exp(-(4/3+2/3*k3+2/3*k))+3.29*10^9*exp(-1/3*(4+4/3*k));
l1=1.22*10^11*m*(10.1*10^9*0.5*exp(-2/3*(0.8+k3))+3.29*10^9*exp(-2/3*(0.8+k)))/(1.22*10^11+m*(10.1*10^9*0.5*exp(-2/3*(0.8+k3))+3.29*10^9*exp(-2/3*(0.8+k))));
l2=1.22*10^11*m*a/(1.22*10^11+m*a);
l3=10.1*10^9*exp(-(0.8+2/3*k3))+3.29*10^9*exp(-(0.8+2/3*k));
y=fun(k3,k,m,l1,l2,l3);
k
y
end
main2.m
fork=17.5:
0.01:
17.8
m=1.109*10^5;
k3=0.42*k;
a=29.7*10^9*0.5*exp(-1/3*(4+2*k3))+10.1*10^9*exp(-(4/3+2/3*k3+2/3*k))+3.29*10^9*exp(-1/3*(4+4/3*k));
l1=1.22*10^11*m*(10.1*10^9*0.5*exp(-2/3*(0.8+k3))+3.29*10^9*exp(-2/3*(0.8+k)))/(1.22*10^11+m*(10.1*10^9*0.5*exp(-2/3*(0.8+k3))+3.29*10^9*exp(-2/3*(0.8+k))));
l2=1.22*10^11*m*a/(1.22*10^11+m*a);
l3=10.1*10^9*exp(-(0.8+2/3*k3))+3.29*10^9*exp(-(0.8+2/3*k));
y=fun(k3,k,m,l1,l2,l3);
k
y
end
附录4:
lingo2算法:
max=17.86*0.42*k/(0.8+0.42*k)*1.22*10^11/(1.22*10^11+n)*n*@exp(-1.6)*(1-@exp(-2/3*(0.8+0.42*k)))+22.99*k/(0.8+k)*1.22*10^11/(1.22*10^11+n)*n*@exp(-0.28*k-2.4)*(1-@exp(-2/3*(0.8+k)));
n=1.22*10^11*(1.109*10^5*(0.5*@exp(-0.28*k-6.4/3)+@exp(-(0.28+2/3)*k-8.8/3))-1);