一类受污染的渔业资源管理模型的最优脉冲控制

Optimal Impulse Control for a Class of Polluted Fishery Resource Management Models

摘要: 本文建立了一类处于污染环境中的渔业捕捞脉冲控制模型。考虑环境受到污染后，鱼类会吸收环境中的污染物，发生持续死亡现象。同时在目标函数中加入处理受污染的鱼类的额外成本，利用最优脉冲控制原理进行理论分析，得出了该模型的最优脉冲捕获时刻和相应的种群水平，最后利用数值模拟，验证所得结论。

Abstract: The paper establishes a fishing pulse control model in a polluted environment, taking into account the absorption of pollutants by fish and their subsequent mortality. Additionally, the objective function incorporates the extra cost of dealing with contaminated fish, and theoretical analysis is conducted using the optimal impulse control principle to determine the optimal capture time and corresponding population level. Finally, numerical simulation is employed to validate the conclusions.

1. 引言

当今世界面临着人口日益增长和资源日渐匮乏的双重挑战，可持续发展理念得到了国际社会的普遍认同。我国政府也依据此制定了国家社会经济发展的长期战略目标。如何对可再生资源进行适当的管理是许多学者十分关心的问题，他们通过研究种群资源的脉冲模型 [1] [2] [3] 来实现这一目标。文献 [4] 给出了脉冲微分系统控制问题相关理论。Blaquiere于1997年将博弈论同脉冲微分方程相结合，得到了最优脉冲控制原理 [5] 。在渔业资源管理领域，赵 [6] 从生态学角度切入，构建了一个Logistic模型，运用最优脉冲控制原理，深入探讨了对于该模型的优化脉冲管理问题，并成功得出了模型的最优脉冲捕获时刻以及相应的捕获水平。与此同时，李 [7] 分别建立连续捕捞与脉冲捕捞相结合的渔业模型，从经济学角度利用最优脉冲控制原理得出了使管理者获得最大经济收益的最优脉冲策略。郑等 [8] [9] 对于麦蚜生态系统和人工湿地中的藻类种群进行了最优脉冲控制分析，为生态系统的污染、病害控制等方面提供了理论指导。作为人类赖以生存的重要资源之一，渔业资源在全球范围内发挥着重要的经济、社会和环境作用。渔业资源通常具有一定再生能力，如果对其合理开发使用，应当是取之不尽，用之不竭的。然而，由于过度捕捞、海洋污染、气候变化等多重因素的影响，全球渔业资源正面临严重威胁。特别是2023年7月日本将福岛核污水排入太平洋后，海洋污染日益严峻。在这一背景下，如何在污染环境中进行合理捕捞，已成为我国渔民越来越关注的问题。在其他资源管理系统 [10] - [16] 的优化控制中，大部分是分析种间竞争，脉冲投放等方式对系统的造成影响，很少考虑污染对于优化目标函数的影响，基于此，本文构建了一个在污染环境下，鱼类种群遵循Logistic增长的渔业资源管理模型，本文考虑污染会对于鱼类和管理者同时造成影响，污染会导致鱼类发生持续不断的死亡现象，环境污染会使管理者需要额外付出处理受污染鱼类的成本。通过运用最优脉冲控制原理，研究了使渔民收益最大的最优脉冲捕捞时刻以及相对应的最优捕捞水平，为渔民进行合理的渔业捕捞提供了理论指导。

2. 模型建立

假设鱼种群的个体具有相同的动态特性并且增长服从Logistic模型规律，建立模型为：

$\{\begin{array}{l}\frac{\text{d}N\left(t\right)}{\text{d}t}=rN\left(t\right)\left(1-\frac{N\left(t\right)}{K}\right),\\ N\left(0\right)=N_0.\end{array}$ (1)

其中， $N\left(t\right)$ 表示时间t时渔场中鱼类种群数量，N₀表示渔场中鱼类的初始数量， $N_0>0$ ，r表示渔场中鱼类的内禀增长率(出生率减去死亡率)，K表示鱼类种群的环境容纳量。在渔场受到污染后，鱼类会将污染物质吸收入体，导致发生持续不断死亡的现象。

此时被污染后的渔业模型为：

$\{\begin{array}{l}\frac{\text{d}N\left(t\right)}{\text{d}t}=rN\left(t\right)\left(1-\frac{N\left(t\right)}{K}\right)-\varphi N\left(t\right),\\ N\left(0\right)=N_0.\end{array}$ (2)

其中 $\varphi$ 表示鱼群的受污染致死率。在现实生活中，渔民会在某一时刻出海大量捕捞鱼类种群，经过脉冲捕获后的模型为：

$\{\begin{array}{l}\frac{\text{d}N\left(t\right)}{\text{d}t}=rN\left(t\right)\left(1-\frac{N\left(t\right)}{K}\right)-\varphi N\left(t\right),t\ne {\tau }_k,\\ \Delta N\left(t\right)=-EN\left(t\right),t={\tau }_k,\\ N\left(0\right)=N_0.\end{array}$ (3)

其中， ${\tau }_k$ 表示脉冲捕获时刻，E表示脉冲捕获系数。

渔民运营渔场的目的通常是为了捕捞贩卖鱼群，以此获取经济收入。考虑渔民管理渔场需要持续付出养殖成本，渔场环境受污染后，渔民处理受污染的鱼类需要付出额外的成本，则假设此时总经济利润(总收入 − 总支出)为：

$W\left(T\right)=\mu \left(IEN\left({\tau }_k\right)-CEN\left({\tau }_k\right)-C_0\right)-{\displaystyle {\int }_0^{T}RN\left(t\right)\text{d}t}$ ，

其中R表示管理鱼群的单位成本， $\mu \in \left[0,1\right]$ 为脉冲变量， $\mu =0$ 时表示无脉冲发生， $\mu =1$ 时表示有脉冲发生。I表示鱼的售价，C表示处理受污染鱼类额外成本，C₀表示进行捕捞的固定成本。设性能指标 $V\left(N_0,s,N(\cdot )\right)=-W\left(T\right)$ ，为了使渔民能获得最大经济收入，即使性能指标V取得最小值。

综上所述，此时可以得到一个环境受污染影响的渔业资源管理模型：

$\{\begin{array}{l}\frac{\text{d}N\left(t\right)}{\text{d}t}=rN\left(t\right)\left(1-\frac{N\left(t\right)}{K}\right)-\varphi N\left(t\right), t\ne {\tau }_k,\\ \Delta N\left(t\right)=-EN\left(t\right), t={\tau }_k,\\ N\left(0\right)=N_0,\\ V\left(N_0,s,N(\cdot )\right)={\displaystyle {\int }_0^{T}RN\left(t\right)\text{d}t}-\mu \left(IEN\left({\tau }_k\right)-CEN\left({\tau }_k\right)-C_0\right).\end{array}$ (4)

3. 模型分析

模型(1)在初值状态下的解为：

$N\left(t\right)=\frac{\left(r-\varphi \right)K}{r+c{\text{e}}^{\left(\varphi -r\right)t}}$ ， $c=\frac{\left(r-\varphi \right)K-N_0}{N_0}$ 。 (5)

根据最优脉冲控制原理 [5] ，此模型的最优脉冲控制问题可以分为两个部分：连续型控制部分，用 $f_i\left(t,N\left(t\right)\right)$ ， $i=1,2$ 和 $\mu \left(N\left(t\right)\right)$ 表示；脉冲型控制部分，用 $g_i\left(t,N\left(t\right)\right)$ ， $i=1,2$ 和 $\mu \left(N\left(t\right)\right)$ 表示。模型可表示为：

$\{\begin{array}{l}{f}_0\left(t,N\left(t\right),u\left(N\left(t\right)\right)\right)=RN\left(t\right),\\ g_0\left(t,N\left(t\right),\mu \left(N\left(t\right)\right)\right)=-\mu \left(IEN\left(t\right)-CEN\left(t\right)-C_0\right),\\ f_1\left(t,N\left(t\right),u\left(N\left(t\right)\right)\right)=rN\left(t\right)\left(1-\frac{N\left(t\right)}{K}\right)-\varphi N\left(t\right),\\ g_1\left(t,N\left(t\right),\mu \left(N\left(t\right)\right)\right)=-\mu EN\left(t\right).\end{array}$ (6)

$\lambda \left(t\right)$ 表示伴随变量， ${\lambda }_0\left(t\right)=1$ ，由于不考虑末端时间，所以 $\lambda \left(T\right)=0$ 。令连续哈密顿函数为 $H\left(\lambda ,N,u\right)$ ，脉冲哈密顿函数为 $H_c\left(\lambda ,N,\mu \right)$ ，由庞特里亚金极小值原理可得：

$H\left(\lambda ,N,\mu \right)={\lambda }_0\left(t\right)f_0\left(N,u\right)+{\lambda }_1\left(t\right)f_1\left(N,u\right)=RN+{\lambda }_1\left(t\right)\left(rN\left(1-\frac{N}{K}\right)-\varphi N\right)$ ，

$H_c\left(\lambda ,N,\mu \right)={\lambda }_0\left(t\right)g_0\left(N,\mu \right)+{\lambda }_1\left(t\right)g_1\left(N,\mu \right)=-\mu \left[\left(I-C\right)EN-C_0+{\lambda }_1\left({\tau }_c+0\right)EN\right]$ 。

4. 最优脉冲控制

接下来对模型进行最优脉冲控制分析。

由文献 [17] 中的庞特里亚金极小值原理有：

$\frac{\text{d}{\lambda }_1}{\text{d}t}=-\frac{\partial H}{\partial N}=-R-{\lambda }_1\left(t\right)\left(r-\frac{2rN\left(t\right)}{K}\right)+{\lambda }_1\left(t\right)\varphi =-R-{\lambda }_1\left(t\right)\left(r-\varphi -\frac{2rN\left(t\right)}{K}\right)$ ， (7)

${\lambda }_1\left({\tau }_c\right)={\lambda }_1\left({\tau }_c+0\right)+\frac{\partial H_c}{\partial N}=\left(1-\mu E\right){\lambda }_1\left({\tau }_c+0\right)+\mu \left(C-I\right)E$ 。 (8)

其中 $N\left(t\right)$ 为模型(1)的解。

将式(5)代入(7)中得：

$\frac{\text{d}{\lambda }_1}{\text{d}t}=-R-{\lambda }_1\left(t\right)\left(r-\varphi \right)\left(1-\frac{2r}{r+c{\text{e}}^{\left(\varphi -r\right)t}}\right)$ 。

解此微分方程可得：

${\lambda }_1\left(t\right)=\frac{R}{r\left(r-\varphi \right)}\left(\frac{1}{c+r{\text{e}}^{\left(r-\varphi \right)t}}-\frac{1}{c+r{\text{e}}^{\left(r-\varphi \right)T}}\right)\frac{{\left(c+r{\text{e}}^{\left(r-\varphi \right)t}\right)}^2}{{\text{e}}^{\left(r-\varphi \right)t}}$ ，

${\lambda }_1\left({\tau }_c+0\right)=\frac{R}{r\left(r-\varphi \right)}\left(\frac{1}{c+r{\text{e}}^{\left(r-\varphi \right){\tau }_c}}-\frac{1}{c+r{\text{e}}^{\left(r-\varphi \right)T}}\right)\frac{{\left(c+r{\text{e}}^{\left(r-\varphi \right){\tau }_c}\right)}^2}{{\text{e}}^{\left(r-\varphi \right){\tau }_c}}$ 。

因为模型中没有施加连续控制，即此时不必考虑连续型哈密顿函数H。现对脉冲型哈密顿函数进行分析。

由文献 [5] 中引理1的(ii)有：

${\min }_{\mu \in K_{\mu }}{H}_c\left(x^{*}\left(t\right),\mu \right)=H_c\left(x^{*}\left(t_c\right),{\pi }^{*}\left(x^{*}\left(t\right)\right)\right)$ ，

为了使这个条件成立，因为

$H_c\left(\lambda ,N,\mu \right)=-\mu \left[\left(I-C\right)EN-C_0+{\lambda }_1\left({\tau }_c+0\right)EN\right]$ ，

所以有以下三种情况：

1) 若 $\left(I-C\right)EN-C_0+{\lambda }_1\left({\tau }_c+0\right)EN>0$ ，则 $\mu =1$ ；

2) 若 $\left(I-C\right)EN-C_0+{\lambda }_1\left({\tau }_c+0\right)EN<0$ ，则 $\mu =0$ ；

3) 若 $\left(I-C\right)EN-C_0+{\lambda }_1\left({\tau }_c+0\right)EN=0$ ， $\mu$ 值不定。

此时定义边界

$\Omega \left({\tau }_c,N\right)=\left(I-C\right)EN-C_0+{\lambda }_1\left({\tau }_c+0\right)EN$ ，

设集合：

$\{\begin{array}{l}{V}_1=\left\{\left({\tau }_c,N\left({\tau }_c\right)\right)|\Omega \left({\tau }_c,N\left({\tau }_c\right)\right)>0,D=\left[0,T\right]\times {ℝ}_{+}\right\},\\ V_2=\left\{\left({\tau }_c,N\left({\tau }_c\right)\right)|\Omega \left({\tau }_c,N\left({\tau }_c\right)\right)<0,D=\left[0,T\right]\times {ℝ}_{+}\right\},\\ V_3=\left\{\left({\tau }_c,N\left({\tau }_c\right)\right)|\Omega \left({\tau }_c,N\left({\tau }_c\right)\right)=0,D=\left[0,T\right]\times {ℝ}_{+}\right\}.\end{array}$

V₁表示有脉冲捕捞发生的区域，V₂表示无脉冲捕捞发生的区域，V₃表示临界情况。

由模型(3)得：

$N\left({\tau }_c+0\right)=\left(1-E\right)N\left({\tau }_c\right)$ 。

此时发生脉冲捕获，即 $\mu ={\pi }^{*}\left(N\right)=1$ 。

根据文献 [5] 中引理的(iii)得：

$\begin{array}{l}{\min }_{u\in K_u}H\left(\lambda \left({\tau }_c+0\right),N^{*}\left({\tau }_c+0\right),u\right)-{\min }_{u\in K_u}H\left(\lambda \left({\tau }_c\right),N^{*}\left({\tau }_c\right),u\right)\\ ={\frac{\partial H_c\left(N,\mu \right)}{\partial N_n}|}_{N=N^{*}\left({\tau }_c\right),\mu ={\pi }^{*}\left(N\right)}\\ =R{N}^{*}\left({\tau }_c+0\right)+\left(r-\varphi \right){\lambda }_1\left({\tau }_c+0\right)N^{*}\left({\tau }_c+0\right)-\frac{r}{K}{\lambda }_1\left({\tau }_c+0\right)N^{*}{\left({\tau }_c+0\right)}^2\\ -R{N}^{*}\left({\tau }_c\right)-\left(r-\varphi \right){\lambda }_1\left({\tau }_c\right)N^{*}\left({\tau }_c\right)+\frac{r}{K}{\lambda }_1\left({\tau }_c\right)N^{*}{\left({\tau }_c\right)}^2\\ =0.\end{array}$

将

$N\left({\tau }_c+0\right)=\left(1-E\right)N\left({\tau }_c\right)$ ， ${\lambda }_1\left({\tau }_c\right)=\left(1-\mu E\right){\lambda }_1\left({\tau }_c+0\right)+\mu \left(C-I\right)E$ ，

代入上式中可得：

$\begin{array}{l}-ER{N}^{*}\left({\tau }_c\right)+\left(r-\right)N^{*}\left({\tau }_c\right)-\left(r-\varphi \right)\left[\left(1-E\right){\lambda }_1\left({\tau }_c+0\right)+\left(C-I\right)E\right]N^{*}\left({\tau }_c\right)\\ -\frac{r}{K}{\lambda }_1\left({\tau }_c+0\right){\left(1-E\right)}^2{N}^{*}{\left({\tau }_c\right)}^2+\frac{r}{K}\left[\left(1-E\right){\lambda }_1\left({\tau }_c+0\right)+\left(C-I\right)E\right]N^{*}{\left({\tau }_c\right)}^2=0\end{array}$

化简后可得：

$ER+\left(C-I\right)E\left(r-\varphi \right)=\frac{r}{K}\left[\left(E-E^2\right){\lambda }_1\left({\tau }_c+0\right)+\left(C-I\right)E\right]N^{*}\left({\tau }_c\right)$ 。

进而得到 $N^{*}\left({\tau }_c\right)=\frac{K\left[R+\left(C-I\right)\left(r-\varphi \right)\right]}{r\left[\left(1-E\right){\lambda }_1\left({\tau }_c+0\right)+\left(C-I\right)\right]}=\frac{K\left[R+\left(C-I\right)\left(r-\varphi \right)\right]}{r\left[{\lambda }_1\left({\tau }_c\right)-\left(C-I\right)\left(1-E\right)\right]}$ 。

令曲线MN为上式所表示的曲线，那么如果模型(3)的解曲线与曲线MN相交，且两曲线的交点在区域V₁内，则在此点所对应的 ${\tau }_c$ 时刻时，渔民对鱼群进行脉冲捕捞可以使渔民获得最高的经济收入。即此时 ${\tau }_c$ 为最优捕获时刻，其对应的种群水平为最优捕获水平。如果解曲线与曲线MN的交点落在区域V₂内，则表示此时不会发生脉冲捕捞。

接下来，求解最优捕获时刻，令 $N^{*}\left({\tau }_c\right)=N\left(t\right)$ 得到下式：

$r\left(1-E\right)R\frac{\left(c+r{\text{e}}^{\left(r-\varphi \right)t}\right)\left({\text{e}}^{\left(r-\varphi \right)T}-{\text{e}}^{\left(r-\varphi \right)t}\right)}{\left(c+r{\text{e}}^{\left(r-\varphi \right)T}\right){\text{e}}^{\left(r-\varphi \right)t}}+\left(C-I\right)\left(r-\varphi \right)=\frac{c+r{\text{e}}^{\left(r-\varphi \right)t}}{{\text{e}}^{\left(r-\varphi \right)t}}\left[R+\left(C-I\right)\left(r-\varphi \right)\right]$ 。

整理后得到：

$\begin{array}{l}\left(1-E\right)R{r}^2{\left({\text{e}}^{\left(r-\varphi \right)t}\right)}^2+\left(2Rrc-r^{2}RE{\text{e}}^{\left(r-\varphi \right)T}+crER\right){\text{e}}^{\left(r-\varphi \right)t}+C^2\left[R+\left(r-\varphi \right)\left(C-I\right)\right]\\ +Cr\left(r-\varphi \right)\left(C-I\right){\text{e}}^{\left(r-\varphi \right)T}+rcER{\text{e}}^{\left(r-\varphi \right)T}=0\end{array}$

所以令 $A=\left(1-E\right)R{r}^2$ ， $B=2Rrc-r^{2}RE{\text{e}}^{\left(r-\varphi \right)T}+crER$ ，

$D=C^2\left[R+\left(r-\varphi \right)\left(C-I\right)\right]+Cr\left(r-\varphi \right)\left(C-I\right){\text{e}}^{\left(r-\varphi \right)T}+rcER{\text{e}}^{\left(r-\varphi \right)T}$

由韦达定理可知，若 $\Delta =B^2-4AD>0$ ，且 $\frac{-B+\sqrt{\Delta }}{2A}>0$ 时，则由上式可解得模型的最优脉冲捕捞时刻 ${\tau }_c^{*}$ 为：

${\tau }_c^{*}=\frac{\ln \left(\frac{-B+\sqrt{B^2-4AD}}{2A}\right)}{r-\varphi }$ ，设 $\theta =\frac{-B+\sqrt{B^2-4AD}}{2A}$ ，

得到

${\lambda }_1\left({\tau }_c+0\right)=\frac{R}{r\left(r-\varphi \right)}\left(\frac{1}{c+r\theta }-\frac{1}{c+r{\text{e}}^{\left(r-\varphi \right)T}}\right)\frac{{\left(c+r\theta \right)}^2}{\theta }$ 。

将上式中的 ${\lambda }_1\left({\tau }_c+0\right)$ 定义为Z，则可以得出符合的最优脉冲捕捞水平：

$N^{*}\left({\tau }_c^{*}\right)=\frac{K\left[R+\left(C-I\right)\left(r-\varphi \right)\right]}{r\left[\left(1-E\right)Z+\left(C-I\right)\right]}$ 。

综上所述，可以归纳出以下定理：

定理1.1：对于模型(3)，如果使以下条件成立：

1) $\Delta =B^2-4AD>0$ ，且 $\theta =\frac{-B+\sqrt{B^2-4AD}}{2A}>0$ ，

2) $\left({\tau }_c^{*},N^{*}\left({\tau }_c^{*}\right)\right)\in V_1$ 。

则模型(3)存在最优脉冲捕获时刻 ${\tau }_c^{*}$ 和相应的最优捕获水平 $N^{*}\left({\tau }_c^{*}\right)$ ，

${\tau }_c^{*}=\frac{\ln \left(\frac{-B+\sqrt{B^2-4AD}}{2A}\right)}{r-\varphi }$ ， $N^{*}\left({\tau }_c^{*}\right)=\frac{K\left[R+\left(C-I\right)\left(r-\varphi \right)\right]}{r\left[\left(1-E\right)Z+\left(C-I\right)\right]}$ 。

其中

$A=\left(1-E\right)R{r}^2$ ， $B=2Rrc-r^{2}RE{\text{e}}^{\left(r-\varphi \right)T}+crER$ ， $D=C^2\left[R+\left(r-\varphi \right)\left(C-I\right)\right]+Cr\left(r-\varphi \right)\left(C-I\right){\text{e}}^{\left(r-\varphi \right)T}+rcER{\text{e}}^{\left(r-\varphi \right)T}$ ， $Z=\frac{R}{r\left(r-\varphi \right)}\left(\frac{1}{c+r\theta }-\frac{1}{c+r{\text{e}}^{\left(r-\varphi \right)T}}\right)\frac{{\left(c+r\theta \right)}^2}{\theta }$ ， $\theta =\frac{-B+\sqrt{B^2-4AD}}{2A}$ 。

5. 数值模拟

本节的目的是验证得出的理论结果。

对于主要结论，选取参数：

$r=0.2$ ， $\varphi =0.01$ ， $N_0=200$ ， $K=6000$ ， $T=50$ ，

$E=0.5$ ， $R=1$ ， $I=5$ ， $C=1$ ， $C_0=1000$ 。

在MATLAB中进行数值模拟后得到图1。

Figure 1. Optimal impulse control curve

图1. 最优脉冲控制曲线

由图1可知，模型(2)的初始解曲线为 $N\left(t\right)$ ，最优脉冲收获曲线为MN，两条曲线在区域V₁内相交于点A，即此时可以发生最优脉冲捕获，此时的时刻为 ${\tau }_c$ ，为渔民可获得最大经济收益的最优脉冲捕获时刻，此时的最优捕获水平为 $N\left({\tau }_c\right)$ 。

6. 主要结论

针对本文的单种群渔业资源管理模型，若能够求解出模型的解析解，那么我们可以通过构建适当的性能指标函数，并结合庞特里亚金极小值原理与最优脉冲控制原理，进一步解出哈密顿函数中的协态变量。通过协态变量，我们可以推导出最优脉冲控制曲线的解析式。最终，通过解析解曲线与最优脉冲控制曲线的交点，我们可以精确地确定模型的最优脉冲控制时刻以及对应的脉冲捕获种群水平。本文所得结论能为处于环境受污染情境下的渔业资源管理模型提供最优脉冲捕捞策略，从而为资源管理者进行合理捕捞活动提供坚实的理论支撑。

NOTES

^*通讯作者。