Bifurcation analysis and optimal control of a delayed single-species fishery economic model

1.   Introduction

The fishery economic model is a kind of model based on biological population- and economic-related equations. It investigates the development of regulations of the biological population and economic income under the influence of the economic value of fishing behavior. By using reasonable parameters to build a model, the properties of the model are studied by using mathematical theory and methods. Therefore, we can reasonably explain the phenomenon in the actual fishing, predict the future trend and put forward guiding opinions on the reasonable maximization of fishing benefits. Some authors ^[1,2] established basic harvesting models to solve optimal harvesting problems. Jerry et al. ^[3] established a harvesting economic model and chose the fishing effort variation rate as the controllor to solve a nonlinear problem of optimal control. Conrad et al. ^[4] proposed and analyzed a mathematical model to study the dynamics of a fishery resource system in an aquatic environment that consists of two zones. Ami et al. ^[5] considered the optimal management problem by defining the stock density of the resource and the sum of discounted benefits as the biological indicator and economic indicator, respectively Bairagi et al. ^[6] performed a qualitative study of the bioeconomic management of a fishery in the presence of some infection.

In ^[7], the author established a single-species fishery economic model based on a logistic model with a constant harvest. They found that the stock change of fish is related to its own growth and the catch rate. Additionally, they found that the price change is affected by the market positive linear demand function and the catch rate. It has a positive correlation with the market positive linear demand function and a negative correlation with the catch rate. So, the fishery economic model is as follows;

where the variables $x(t)$ and $p(t)$ denote the fish stock and the unit price of the stock at time $t$, respectively. $D(p)$ is the positive linear demand function such that $D(p) = a-p(t)\ge0$ ^[8]. The parameter $r$ is the intrinsic growth rate of the biomass. $k$ is the carrying capacity of the environment. $Y$ is the catch rate. $a$ is the market capacity. $s$ is the price speed adjustment. $r$, $k$, $Y$, $a$ and $s$ are positive constants.

Indeed, the catch rate is affected by the total revenue and the total cost. Assume that the catch rate is positively correlated with the total revenue and negative correlated with the total cost. If the cost is high, the economic benefits of engaging in fishery production will be reduced accordingly, the attraction of engaging in relevant industries will be reduced, the human capital will flow to other fields, and the catch rate will be reduced accordingly. So, only considering the catch rate as a constant simply is not very reasonable. Therefore the catch rate should be a variable at time $t$, and it is represented by $y(t)$. So the following fishery economic model is obtained:

where $\beta$ is the response coefficient, $c$ is the cost of catching. $\beta$ and $c$ are positive constants.

Considering that the development of the system depends on both the current state and the past state, it is necessary to consider time delays in the model. Chakraborty et al. ^[9] introduced a single discrete gestation delay in a differential-algebraic bioeconomic system and investigated Hopf bifurcation in the neighborhood of the coexisting equilibrium point. Song et al. ^[10] determined the direction of the Hopf bifurcation and the stability of the bifurcating periodic solutions of a predator-prey system with two delays. Liu et al. ^[11] proposed a delayed Gause predator-prey model with Michaelis-Menten type harvesting and derived the conditions of local stability and Hopf bifurcation. Zhang et al. ^[12] determined the direction of Hopf bifurcation and the stability of bifurcated periodic solutions of a bioeconomic predator-prey model.

In real life, there is a period of time from fishermen's capture to sale. That means that the fishermen's capture rate is also affected by the price at a certain time in the past, which is recorded as $\tau_{1}$ in the following model. Considering that there is a certain time delay between the actual market information and the market information acquired by buyers, a time delay in the process of the price affecting market demand exists, which is recorded as $\tau_{2}$. So we get the following delayed single-species fishery economic model:

where $\tau_{1}\geq0$ and $\tau_{2}\geq0$ are the time delays, and the other parameters are similar to ones of System (1.2).

The initial conditions of the delayed single-species fishery economic system given by System (1.3) are

with $\tau = max{\{\tau_{1}, \tau_{2}\}}$.

The rest of the paper is organized as follows: In Section 2, the conditions of local stability and Hopf bifurcation are discussed as functions of the time delays in different intervals. In Section 3, we investigate the direction of Hopf bifurcation and the stability of bifurcated periodic solutions. In Section 4, using cost control, the optimal harvesting of fish stocks is considered to maximize the net profit and minimize the waste caused by hoarding for speculation, while ensuring the sustainable survival of fish stocks. Besides, a control system with time delays is established, which is about guiding interventions and aims to reduce the waste. A realistic Penaeus vannamei's cultivation model simulation is demonstrated to prove the validity of the theoretical analysis in Section 5. Finally, the conclusions are presented in Section 6.

2.   Local stability analysis and bifurcation

In this section, one concentrates on the local stability and Hopf bifurcation phenomenon around the equilibria of System (1.3), as the time delays take different values.

System (1.3) has non-negative equilibria $S_0 = (0, 0, a)$, $\hat{S}_0 = (k, 0, a)$ and $S = (x, y, p)$, where $p = c$ and $y = a-c$, and $x$ satisfies the following equation:

After a simple calculation, it can be obtained that $x_{1} = \frac{k}{2}+ \frac{k}{2r}\sqrt{r^{2}-\frac{4r}{k}(a-c)}$, $x_{2} = \frac{k}{2}- \frac{k}{2r}\sqrt{r^{2}-\frac{4r}{k}(a-c)}$. There are two positive interior equilibria $S_{1, 2} = (x_{1, 2}, a-c, c)$ when $a-\frac{rk}{4} < c < a$. The two internal equilibria of the system are merged into one, which is called the degenerate equilibrium $(x^*_{0}, a-c, c)$ when $c = a-\frac{rk}{4}$. And there is no interior equilibrium when $c < a-\frac{rk}{4}$.

The characteristic equation of System (1.3) at the equilibrium $S_0$ can be expressed as follows:

and it always has two characteristic values $\lambda_1 = r > 0$ and $\lambda_2 = \beta(a-c) > 0$. So $S_0$ is always unstable. The characteristic equation of System (1.3) at the equilibrium $\hat{S}_0$ can be expressed as follows:

and it always has two characteristic values $\lambda_1 = -r < 0$ and $\lambda_2 = \beta(a-c) > 0$. So $\hat{S}_0$ is always unstable.

And the characteristic equation of System (1.3) at the equilibria $S_{1, 2}$ can be expressed as follows:

The characteristic Eq (2.1) always has a characteristic value $\lambda^* = r(1-\frac{2x}{k})$, and $\lambda^*\mid_{S_1} < 0$, $\lambda^*\mid_{S_2} > 0$. So $S_2$ is always unstable. In order to study the stability of the equilibrium $S_1$, we need to further study other characteristic values.

The other characteristic values are satisfied:

In the following section, we will investigate the dynamic behavior of System (1.3) around $S_1$ (we denote $S_{*}$) according to the time delays $\tau_{1}$ and $\tau_{2}$ with different values, respectively.

Case 1. $\tau_{1} = \tau_{2} = 0$.

Equation (2.2) becomes

The two eigenvalues have negative real parts, so we can get the following theorem:

Theorem 1. System (1.3)is asymptotically stablearound the interior equilibrium $S_{*}$.

Remark 1. It means that in the absence of a time delay, the system will remain stable, that is, the fish stock will remain at a stable quantity without extinction; the capture rate and the unit price of the stock will also remain stable, and they are all stable at the value of the internal equilibrium.

Case 2. $\tau_{1} = \tau_{2} = \tau > 0$.

Equation (2.2) becomes

It is assumed that for some values of $\tau > 0$, there exists a real number $\omega > 0$ such that $\lambda = \pm i\omega$ are two purely imaginary roots of Eq (2.3). Substituting $\lambda = i\omega$ into Eq (2.3) and separating out the real and imaginary parts, it follows that

From Eq (2.4), one obtains that

There is a unique real and positive solution $\omega^{*2}$ of the above equation.

So the characteristic Eq (2.3) has a pair of purely imaginary roots $\pm i\omega^{*}$ for $\tau_{k}^{*}$. The values of $\tau_{k}^{*}$ are calculated as follows:

From ^[13], by using a method of analyzing the exponential polynomial zero distribution, we get System (1.3) is locally asymptotically stable around $S_{*}$ for $\tau\in(0, \tau_{0}^{*})$.

From Eq (2.3), we can verify the following transversality conditions ^[14]

Theorem 2. System (1.3) is asymptotically stable around the equilibrium $S_{*}$for $\tau\in(0, \tau_{0}^{*})$ and unstable for $\tau > \tau_{0}^{*}$.The system undergoes a bifurcation at $\tau = \tau_{0}^{*}$.

Case 3. $\tau_{1} = 0$ and $\tau_{2} > 0$.

Equation (2.2) becomes

Substituting $\lambda = i\omega$ into Eq (2.6) and separating out the real and imaginary parts, it follows that

Similar to Case 2, we get

The above equation has two real and positive roots

So, the characteristic equation only has two pairs of imaginary roots $\pm i\omega_{\pm}^{*}$. The values of $\tau_{2k\pm}^{*}$ are given by Eq (2.8), which is calculated from Eq (2.7).

By calculating from Eq (2.6), it is obtained that

Subsequently, we can verify the following transversality conditions

By using the Butlers lemma ^[15], we can obtain the following theorem:

Theorem 3. System (1.3) is asymptotically stable around the equilibrium $S_{*}$for $\tau_{2}\in(0, \tau_{20+}^{*})\bigcup(\bigcup\limits_{k = 0}^{j-1}(\tau_{2k-}^*, \tau_{2k+1+}^*))$, andunstable for $\tau_{2}\in(\bigcup\limits_{k = 0}^{j-1}(\tau_{2k+}^*, \tau_{2k-}^*))\bigcup(\tau_{2j+}^{*}, +\infty), j > 0$.The system undergoes a Hopf bifurcation at $\tau_{2} = \tau_{20\pm}^{*}$.

Case 4. $\tau_{1} > 0$ and $\tau_{2} = 0$.

The calculations are similar to that for Case 1, so we will only list the theorem.

where $\omega^{*}$ is satisfied with

Theorem 4. System (1.3) is asymptotically stable around the equilibrium $S_{*}$for $\tau_{1}\in(0, \tau_{10}^{*})$ andunstable for $\tau_{1} > \tau_{10}^{*}$.The system undergoes a Hopf bifurcation at $\tau_{1} = \tau_{10}^{*}$.

Case 5. $\tau_{1} > 0$, $\tau_{2}\in(0, \tau_{20+}^{*})$, $\tau_{1}\neq\tau_{2}$.

In this subsection, $\tau_{1}$ is considered to be a bifurcation parameter and $\tau_{2}$ is confined to a range of $(0, \tau_{20+}^{*})$, where $\tau_{20+}^{*}$ is determined by Eq (2.8). Substituting $\lambda = i\omega$ into Eq (2.2) and separating out the real and imaginary parts, it follows that

From Eq (2.10), ithe following is obtained:

where

Denote $F'(\omega) = 4\omega^{3}+3A_1\omega^{2}+2A_2\omega$. Equation (2.11) has a unique real and positive root $\omega^{*}$ if Inequality (2.12) is satisfied. Otherwise, Eq (2.11) has one positive and real root at least.

where $\omega_* = \frac{-3A_1+\sqrt{9A_1^2-24A_2}}{12}$. So, Eq (2.2) has imaginary roots if Inequality (2.12) is satisfied, and the roots of Eq (2.2) have a negative real part when $\tau_{1}\in(0, \tilde{\tau}_{0})$. $\tilde{\tau}_{0}$ is given by the following equation:

Then, differentiating Eq (2.2) at $\tau_{1} = \tilde{\tau}_{0}$ and separating out the real and imaginary parts, it follows that

where

Thus, we get

if

If Inequality (2.12) is not satisfied, it can only hold that System (1.3) is locally asymptotically stable for $\tau_{1}\in(0, \hat{\tilde{\tau}}_{0})$, where $\hat{\tilde{\tau}}_{0}$ stands for the smallest $\tilde{\tau}_{0}$ corresponding to all positive roots of Eq (2.11).

Theorem 5. If Eqs (2.12) and (2.14) are satisfied, System (1.3) is asymptotically stable for$\tau_{1}\in(0, \tilde{\tau}_{0})$ andunstable for $\tau_{1} > \tilde{\tau}_{0}$, ($\tilde{\tau}_{0}$ is defined in Eq (2.13)).The system undergoes a Hopf bifurcation at $\tau_{1} = \tilde{\tau}_{0}$.

Case 6. $\tau_{2} > 0$, $\tau_{1}\in(0, \tau_{10}^{*})$, $\tau_{1}\neq\tau_{2}$.

$\tau_{10}^{*}$ is determined by Eq (2.9). The calculations are similar to those for Case 5, so we will only state the theorem.

Theorem 6. If Eqs (2.15) and (2.17) are satisfied, the system isasymptotically stable around the equilibrium $S_{*}$ for $\tau_{2}\in(0, \tilde{\tau}_{20})$($\tilde{\tau}_{20}$ is defined in Eq (2.16))and undergoes a Hopf bifurcation at $\tau_{2} = \tilde{\tau}_{20}$.

and $\omega^{*}$ is satisfied with

where

Remark 2. If the time delays are nonzero, corresponding to each case, when the time delays are small, the system will remain stable; when they are slightly greater than the bifurcation values, the system will undergo Hopf bifurcation. At this time, the fish stock, the capture rate and the unit price of the stock will fluctuate within a certain range and cannot reach stability. Because of these fluctuations, the fish stock cannot be stably supplied to the market, and the price fluctuation of fishery resources will also be relatively large, which is not conducive to the stability of the fishery economy. When the time delays are particularly large, the system will be unstable. At this time, the fluctuation ranges of the fish stock, the capture rate and the unit price will be larger and larger, which is not conducive to the stability of the fishery economy.

3.   Stability of bifurcated periodic solutions

In this part, we investigate the direction of Hopf bifurcation and the stability of bifurcated periodic solutions. For the delay differential equations, when the conditions of the Hopf bifurcation theorem are satisfied, the calculation formulas for determining the direction of Hopf bifurcation and the stability of bifurcated periodic solutions can be given by applying the central manifold theory and gauge type method.

We choose $\tau_{1} = \tilde{\tau}_{0}+\mu$ as the bifurcation parameter, $\mu = 0$ is the Hopf bifurcation value of System (1.3) and $\tilde{\tau}_{0}$ is defined in Eq (2.13). Without the loss of generality, we assume that $\tau_{2} < \tilde{\tau}_{0}$, such that $-1 < -\frac{\tau_{2}}{\tilde{\tau}_{0}} < 0$. Let

and $t = t\tilde{\tau}_{0}$, where $(x_{*}, y_{*}, p_{*})$ is the interior equilibrium $S_1$. System (1.3) is expressed in the phase space $\mathbb{C}: = \mathbb{C(}[-1, 0], \mathbb{R}^{3})$ which has the following vector form:

where $u_{t} = (u_{1}(t), u_{2}(t), u_{3}(t))^{T}\in \mathbb{R}^{3}$, $L_{\mu}:\mathbb{C}\rightarrow \mathbb{R}^{3}$ and $F: \mathbb{R}\times \mathbb{C}\rightarrow \mathbb{R}^{3}$ are the linear part and nonlinear part respectively,

where $\Phi = (\Phi_{1}, \Phi_{2}, \Phi_{3})^{T}\in \mathbb{C}$. By the Riesz representation theorem, there exists a $3\times3$ matrix $\eta(\theta, \mu):[-1, 0]\rightarrow \mathbb{R}^{3\times3}$ of bounded variation functions, such that

and $\eta(\theta, \mu)$ in our system can be selected as follows:

where $F_{1} = -\frac{r}{k}\Phi_{1}^{2}(0)$, $F_{2} = \frac{\beta}{2}\Phi_{2}(0)\Phi_{3}(-1)$ and $F_{3} = 0$. The infinitely small generator $A(\mu)$ corresponding to the linearization part of System (1.3) is

Define $R\Phi(\theta)$ as

Then, Eq (3.1) is equivalent to the abstract ordinary differential equation

Defining the formal adjoint matrix of $A(\mu_{0})$ as $A^{*}$

For $\Psi\in \mathbb{C}([0, 1], \mathbb{R}^{3})$ and $\Phi\in \mathbb{C}([-1, 0], \mathbb{R}^{3})$, we define a bilinear form suitable for the complex vector

Since $A^{*}$ is the formal adjoint matrix of $A(\mu_{0})$, we can obtain that

Since $\pm i\omega^{*}\tilde{\tau}_{0}$ are eigenvalues of $A(0)$, they are also eigenvalues of $A^{*}$. We suppose that $q(\theta) = (q_{1}, q_{2}, 1)^{T} e^{i\omega^{*}\tilde{\tau}_{0}\theta}$ is the eigenvector of $A(0)$ corresponding to $i\omega^{*}\tilde{\tau}_{0}$, and $q^{*}(s) = D(q_{1}^{*}, q_{2}^{*}, 1)e^{i\omega^{*}\tilde{\tau}_{0}s}$ is the eigenvector of $A^{*}(0)$ corresponding to $-i\omega^{*}\tilde{\tau}_{0}$, where $D$ is a nonzero coefficient. And they satisfy $\langle q^{*}, q\rangle = 1$. It can be obtained that

from the formal adjoint matrix $A^{*}$. Through a simple calculation, we get

$q_{1} = \frac{-\frac{i\omega^{*}}{s}-e^{-i\omega^{*}\tilde{\tau}_{0}}}{r(1-\frac{2x_{*}}{k})-i\omega^{*}}$, $q_{2} = -\frac{i\omega^{*}}{s}-e^{-i\omega^{*}\tilde{\tau}_{0}}$. Similarly, we can obtain $q_{1}^{*} = \frac{s\omega^{*}i e^{i\omega^{*}\tau_{2}}+\omega^{*2}}{\beta(a-c)e^{i\omega^{*}\tilde{\tau}_{0}}}-s$, $q_{2}^{*} = \frac{se^{i\omega^{*}\tau_{2}}-i\omega^{*}}{\beta(a-c)e^{i\omega^{*}\tilde{\tau}_{0}}}$. Because $\langle q^{*}, q\rangle = 1$, we get

Next, we will realize spectral decomposition. By the center manifold reduction ^[16], we define

where $u_{t}$ is the solution of Eq (3.1) for $\mu = 0$. We denote

so, $W(t, \theta)\in Q_{\pm i\omega^{*}\tilde{\tau}_{0}}$. $z$, $\bar{z}$ are local coordinates of the center manifold $C_{0}$ in the direction of $q^{*}$ and $\bar{q^{*}}$, where $z$ and $\bar{z}$ are conjugate complex numbers. On the central manifold $C_{0}$, $W(t, \theta) = W(z, \bar{z}, \theta)$, and we only consider the real solutions here, where $W(z, \bar{z}, \theta)$ can be written in the form of a power series for $z$ and $\bar{z}$:

The solution $u_{t}\in C_{0}$ of Eq (3.1) can be calculated from

where

The above equation can be rewritten as

where $g(z, \bar{z})$ is denoted by

So, we can obtain that

From Eq (3.7), we can obtain that

where $u_{t}(\theta) = (u_{1t}(\theta), u_{2t}(\theta), u_{3t}(\theta))^{T}$, and $W(t, \theta) = (W^{(1)}(t, \theta), W^{(2)}(t, \theta), W^{(3)}(t, \theta))^{T}$. We can calculate the parameter expression from Eq (3.10), which can determine the direction of Hopf bifurcation, the stability of periodic solutions and the increase or decrease of the period of bifurcating periodic solutions. It can be further obtained that

From above, it can be found that $g_{21}$ depends on the coefficients $W_{20}(\theta)$ and $W_{11}(\theta)$ of $W(z, \bar{z}, \theta)$, while $g_{20}$, $g_{11}$, $g_{02}$ do not depend on $W(z, \bar{z}, \theta)$. Next, we calculate $W_{20}(\theta)$ and $W_{11}(\theta)$.

For Eq (3.7), we get $\dot{W} = \dot{u}_{t}-\dot{z}q-\dot{\bar{z}}\bar{q}$. Combining with the definitions of the infinitely small generator $A(\mu)$, $R\Phi(\theta)$ and that given by (3.9), we can get

In addition, on the central manifold $C_{\mu_{0}}$, we can obtain that

Through comparing the coefficients of the items $\frac{z^{2}}{2}$ and $z\bar{z}$ between Eq (3.12) and Eq (3.11), we can obtain that

and

where

For $\theta\in[-1, 0)$, from Eq (3.13), we can obtain that

Substituting $q(\theta) = q(0)e^{i\omega^{*}\tilde{\tau}_{0}\theta}$ into Eq (3.15), we can obtain that

For $\theta = 0$, from Eq (3.13), we can obtain that

Substitute Eq (3.16) into the above equation to obtain

which gives

For $\theta\in[-1, 0)$, from Eq (3.14), we can obtain that

From further calculation, we can obtain that

Similarly, we can obtain that

Then, we can judge the properties of Hopf bifurcation by the parameters $\mu_{2}$, $\beta_{2}$ and $T_{2}$. We denote

then,

So we can obtain the conclusion as follows.

Theorem 7. The parameters $\mu_{2}$, $\beta_{2}$ and $T_{2}$ determine theproperties of Hopf bifurcation. $\mu_{2}$ determines the direction of the Hopfbifurcation: if $\mu_{2} > 0$ ($\mu_{2} < 0$), the Hopf bifurcation is supercritical(subcritical); $\beta_{2}$ determines the stability of the bifurcating periodicsolution: the bifurcating periodic solution is stable (unstable) if $\beta_{2} < 0$($\beta_{2} > 0$); $T_{2}$ determines the period of the bifurcating periodicsolution: the period will increase (decrease) if $T_{2} > 0$ ($T_{2} < 0$).

4.   Optimal control problems

In this section, we investigate the optimal control strategy based on System (1.3). The optimal control problem is to seek the basic principle of control and the limitation of the available principle ^[17,18]. It is widely used in aerospace, electronic information-related fields, bioengineering, economic management and other fields ^[19,20,21]. In this section, two optimal control problems with different control variables are considered.

4.1.   Optimal harvesting and cost control

We first consider System (1.3) with time delays $\tau_{1} = \tau_{2} = 0$. To control the cost, we consider the optimal harvesting of fish stocks to maximize the net economic income and minimize the waste caused by the stored amount of fishery resources. Assume that the probability that the harvesting is greater than the market demand is $p_1$($0\leq p_1\leq1$). When the market demand is greater than the harvesting, then the sales is equal to harvesting and the total economic income is $y(t)(1-p_1)p(t)$; otherwise, the sales is equal to market demand and the total economic income is $(a-p(t))p_1p(t)$. And the total cost is $y(t)c(t)$. Then the net economic income can be expressed as follows:

In the process of marine fish breeding, the management of the fishing ground is bound to produce corresponding resource management costs. When the monthly fish stock is greater than the catch rate, the excess fish stock is hoarded for speculation. The larger the value of this part, the more management costs will be incurred in theory. We do not want the waste caused by the stock of fish resources. The waste caused by monthly hoarding for speculation is represented by$x(t)-y(t)$, and $A$ is the corresponding positive weight parameter. For the fishing process, the fishery administration department can adjust the costs incurred during the fishing process to a certain extent by adjusting the relevant tariff preferences, ocean fishing vessel construction subsidies, fuel subsidies and other support policies. The optimal control problem involves determining the control variable $c$ to maximize the objective function

where $\delta$ is the instantaneous annual discount rate and $p(x)$, $x(t)$ and $y(t)$ are state variables of System (1.3). Using the Pontryagin maximization principle, we construct the Hamiltonian function as follows

where $\lambda_{1}$, $\lambda_{2}$ and $\lambda_{3}$ are adjoint variables and $c$ is the control variable, which can be varied within the range $a-\frac{rk}{4} < c(t) < a$. The condition for a singular control to be optimal is

From the above equation, we can obtain that $\lambda_{2}(t) = -\frac{1}{\beta}e^{-\delta t}$. For the adjoint variables $\lambda_{1}$, $\lambda_{2}$ and $\lambda_{3}$, the following equations can be satisfied:

By substituting $\lambda_{2}(t) = -\frac{1}{\beta}e^{-\delta t}$ into Eq (4.3), it can be obtained that

Then adding the integrating factor $e^{-st}$ ^[22], we can obtain that

By substituting $\lambda_{2}$ and $\lambda_{3}$ into Eq (4.2), we can obtain that

By substituting Eq (4.5) into Eq (4.1), we can obtain that

Equation (4.6) is the optimal control path $S_{\delta}(x_{\delta}, y_{\delta}, p_{\delta})$ of System (1.2) under the control of cost. And the optimal cost and the optimal harvesting are as follows:

Equation (4.6) provides an equation for the singular path and gives the optimal equilibrium levels of $x = x_{\delta}, y = y_{\delta}$ and $p = p_{\delta}$. $\lambda_{i}e^{\delta t}, i = 1, 2, 3$ represents the shadow prices along the singular path. From Eqs (4.4) and (4.5) and $\lambda_{2}(t) = -\frac{1}{\beta}e^{-\delta t}$, it may be concluded that these shadow prices may remain constant over the time interval in an optimal equilibrium when they satisfy the strict transversality condition at $\infty$. Further, they remain bounded when $t\rightarrow \infty$. There exist an optimal control $c_{\delta}$ and corresponding solutions $x_{\delta}, y_{\delta}$ and $p_{\delta}$ that maximize the objective function. The optimal control path can involve any combination of the three variables $x_{\delta}, y_{\delta}$ and $p_{\delta}$ that satisfy the constraints.

4.2.   Optimal control with guiding interventions

The fisheries and fishery administrations play a crucial role in the process of cultivation, harvesting and the market circulation of aquatic products. They put forward guiding interventions to assist cultivation enterprises and individual farmers in finding a reasonable scale for farming scientifically, and they publicize marketing strategies for aquatic products reasonably and moderately. By doing so, companies and individuals will reduce the waste caused by hoarding for speculation or income loss caused by too few products, and maximize their interests, boosting the economic growth of the country.

According to the above idea, System (1.3) is generalized by incorporating two controls. The controls $u_1(t)$ and $u_2(t)$ stand for the influence of guiding interventions on finding reasonable scale for farming and marketing strategies, respectively. Regarding finding a reasonable scale for farming, the government adjusts the breeding scale by increasing subsidies to fish fry suppliers and feed suppliers, thereby reducing the loan interest of enterprises of this product, and opening directional loans. Regarding marketing strategies, the regulatory authorities adjust the market price of the commodity by increasing import tariffs, limiting the import quantities of similar agricultural products, implementing the export tax rebate policy and refunding the domestic tax when the commodity is declared for export. At the same time, it also affects the capture rate for fishermen and farms. $C_1$ and $C_2$ are positive weight parameters. So, we obtain the following system:

The state vector of System (4.7) is given by

Considering the biological significance, we assume that the initial conditions of System (4.7) are given as follows:

where $\tau = max(\tau_1, \tau_2)$. The optimal control problem involves determining the control variables $u_1(t)$ and $u_2(t)$ to minimize the objective function

where $x(t)-y(t)$ is the waste caused by hoarding for speculation and $A$ is positive weight parameter. The squares of the control variables reflect the severity of the side effects of guiding interventions on finding a reasonable scale for farming and marketing strategies; $B_1$ and $B_2$ are positive weight parameters that are associated with them. The variables $u_1(t), u_2(t)\in U$ and $U$ represent the control set defined by

where $u_i^{max}$ represents the maximum attainable values of $u_{i}$. To solve the optimal control problem, we must prove the existence of the optimal control first. It will be proved if the following conditions are satisfied:

$(1)$ The set of controls and state variables is nonempty.

$(2)$ The control space is closed and convex.

$(3)$ The right side of System (4.7) is bounded by a linear function with the state and control.

$(4)$ The integrand in the objective function is convex with respect to the input controls $u_1$ and $u_2$.

$(5)$ There exists a constant $D_1 > 1$ and positive numbers $D_2$ and $D_3$ such that the integrand of the objective functional satisfies

By neglecting the negative terms in System (4.7), we obtain

From the third equation of the above system, the solution of $p(t)$ in finite time is bounded and exists according to the comparison theorem ^[22]. System (4.8) can be rewritten in the vector form as follows:

where $p^{max}$ is the maximum of $p$. This system is linear in finite time with bounded coefficients. Then the solutions of this linear system are uniformly bounded. Therefore, the solutions of the nonlinear system given by System (4.8) are bounded and exist. Hence, Condition $1$ is satisfied. Condition $2$ is satisfied by the definition of $U$. System (4.7) can be rewritten as

where $X(t) = (x(t), y(t), p(t))^T, X_{\tau_1}(t) = (x(t-\tau_1), y(t-\tau_1), p(t-\tau_1))^T$, $X_{\tau_2} = (x(t-\tau_2), y(t-\tau_2), p(t-\tau_2))^T$ are the vectors of the state variables. According to the Hölder inequality ^[23], $X_1, X_{\tau_11}, X_{\tau_21}, X_2, X_{\tau_12}$ and $X_{\tau_22}$ exist and satisfy

where

$Q = max\{\frac{2r}{k}x^{max}, \beta p^{max}\} < \infty$. This means that $G(X)$ is uniformly Lipschitz continuous, Condition $3$ is satisfied. Condition $4$ is verified by the definition. Finally,

where $D_2 = max\{\frac{B_1}{2}, \frac{B_2}{2}\}$, $D_1 = 2$, $D_3 > 0$. Condition $5$ is satisfied. So the following theorem is obtained.

Theorem 8. There exist optimal control functions $u_{1\delta}$, $u_{2\delta}$ and a setof corresponding solutions $x_{\delta}(t)$, $y_{\delta}(t)$, $p_{\delta}(t)$ sothat $J(u_{1\delta}, u_{2\delta}) = minJ(u_1, u_2)$, $u_1, u_2\in U$.

Then, using Pontryagin's maximization principle, we construct the Hamiltonian function as follows

Let $\chi_{[0, t_f-\tau]}(t)$ be the indicator function of the interval $[0, t_f-\tau]$

By using optimality, we can obtain that

So, we can find that

with $x(0) = x_0$, $y(0) = y_0$, $p(0) = p_0$. The adjoint equations and transversality conditions are obtained using

It is obtained that

where $(x_\delta, y_\delta, p_\delta)$ is the optimal solution for the optimal controls $u_{1\delta}(t)$, $u_{2\delta}(t)$, and they satisfy System (4.7).

5.   Numerical simulations

Penaeus vannamei is an important commercial shrimp and there are many experiments and studies on its cultivation ^[24,25]. There are also some reports on the breeding, fishing and marketing of Penaeus vannamei ^[26]. In this section, the above theoretical results will be applied to the culture of Penaeus vannamei.

For System (1.3), we take $r = 0.86$, $k = 700$, $\beta = 0.4$, $c = 50.07$, $a = 88.92$, $s = 0.05$ and the set initial value $X_0(x_0, y_0, p_0) = (650, 38, 50)$. We take the month as the time unit. The value of $c$ is the unit cost (ten thousand yuan) and the value of $a$ is the monthly market demand (tons); their values are derived from the preliminary treatment of the fishing situation investigation results of Penaeus vannamei in ^[26], where $c$ and $a$, respectively, are the average cost and total sales of Penaeus vannamei in China. The value $r$ is converted and averaged according to the survival rate of the seedlings in several environments in a culture experiment report ^[24] to obtain an estimated value. The above-mentioned parameter values are in line with the definition of the parameters in our model. By calculation, it can be obtained that $S_{1} = (651.4596, 38.8500, 50.0700)$ and $S_2 = (48.5404, 38.8500, 50.0700)$. When $\tau_{1} = \tau_{2} = 0$, there is only one equilibrium $S_{1}$, according to Theorem 1, the internal equilibrium point $S_{1}$ is globally asymptotically stable, which can be seen in Figure 1. The points near $S_2$ always reach a steady-state $S_1$, which can be seen in Figure 2. From Case 1 and Eq (2.5), it can be obtained that $\omega^* = 2.8100$ and $\tau_{0}^{*}\approx0.0637$. So, the interior equilibrium $S_{1}$ remains stable for $\tau < \tau_{0}^{*}$, which can be seen in Figure 3. As $\tau$ increases through $\tau_{0}^{*}$, Hopf bifurcation occurs. The bifurcating periodic solution exists for $\tau$ slightly larger than $\tau_{0}^{*}$ which can be seen in Figure 4. The phase space trajectory shows a limit cycle corresponding to the periodic solution in the solution curves, which forms around the fixed point. For $\tau > \tau_{0}^{*}$, the interior equilibrium becomes unstable, which can be seen in Figure 5. The solution of System (1.3) varies as the bifurcation parameter $\tau$ for $\tau_{1} = \tau_{2}$, which is shown in Figure 6. And from this figure, we can observe that the solution changes from stable to unstable as the bifurcation parameter $\tau$ increases through the bifurcation value $\tau_{0}^{*}$. From Eq (2.8), it can be obtained that $\tau_{20+}^*\approx0.1767, \tau_{20-}^*\approx0.2276$. From Eq (2.9), it can be obtained $\tau_{10}^*\approx0.2054$. The graphs that show dynamical variation of System (1.3) when $\tau_1\neq\tau_2\neq0$ are similar to $\tau_{1} = \tau_{2}$.

For case 5, we fix $\tau_2 = 0.0670 < \tau_{20+}^*$; by calculation, it can be obtained that $\omega^*\approx2.8640, \tilde{\tau}_0\approx0.0598$. Based on our analysis in Section 3, we can compute the crucial values with the help of Matlab to get $\mu_2 = -2.7672\times10^{-4}$, $\beta_2 = 7.4396\times10^{-4}$, $T_2 = 2.514610^{-5}$. Therefore, according to Theorem 7, when $\tau_1 = \tau_{0}^{*} = 0.0598$, it can be concluded that the direction of the local Hopf bifurcation is subcritical additionally, the nontrivial periodic solutions bifurcating from the interior equilibrium $S_{1} = (651.4596, 38.8500, 50.0700)$ are unstable and increase on the center manifold, which can be seen in Figure 7. The solution of System (1.3) varies according to the bifurcation parameter $\tau_1$, for $\tau_{2} = 0.0670$, which is shown in Figure 8.

The optimal control problem of System (4.7) cannot be solved analytically; ,consequently, reliable numerical methods are essentially required. It is transformed into a nonlinear programming problem. The rough steps are as follows. Let $A = 0.001$, $B_1 = 1000$, $B_2 = 10$, $C_1 = C_2 = 10$. Under the assumption that there exist a step size $h > 0$ and integers $(n, g_1, g_2)\in\mathbb{N}^3$ with $t_f = nh$, $\tau_1 = g_1h$ and $\tau_2 = g_2h$, let $\tau_1 = 0.05$, $\tau_2 = 0.5$, $h = 0.01$, $t_f = 300$, so that $n = 3000, g_1 = 1, g_2 = 10$ and other parameters are the same as above. We set $X(\theta) = (x^{\theta/h}, y^{\theta/h}, p^{\theta/h}) = X_0, \theta\in[-\tau, 0]$ and $\lambda_j(\theta) = \lambda_j^{\theta/h} = 0$, where, $\theta\in[t_f, t_f+\tau], j = 1, 2, 3$, where $\tau = max\{\tau_1, \tau_2\}$. By using backward difference approximation, adjoint functions with transversality conditions can be obtained as

By utilizing combinations of the forward and backward difference approximations, it can be derived that

The optimal harvest controls are updated by values of the state and adjoint variables

where $u_1^{max} = u_2^{max} = 1$. Repeat the above steps for $i = 0, ..., n-1$. From these steps, we obtain Figure 9. From this figure we can get that under the reasonable controls, the waste caused by hoarding for speculation is reduced. We also get that the interior equilibrium of System (1.3) is stable when the cost is $c = 56.8426$ and the time delay $\tau_{1} = \tau_{2} < \tau_{0}^{*}$, which is shown in Figure 9.

6.   Conclusions

In this paper, a delayed single-species fishery economic model was established. It is considered that the catch rate was affected by the total revenue and the total cost. Two time delays were added to the system; they were found to be related to the effects of price on the catch rate and the market requirements respectively. For the delayed system, the local stability behaviors around the interior equilibrium point were discussed. It could be seen that the time delays had a great influence on the stability of the system. When the time delays are too large, the system will change from stable to unstable. This is reflected in the sharp fluctuations of fishery resources and prices in fisheries, which we do not want to see. If we can reduce these two time delays as much as possible, we can increase the stability of the fishery industry. For example, we would be able to feedback market information to fishermen and farms in real time, improve price transparency and reduce the time delay of price in the fishing process. For another example, the pricing of the current year is affected by the previous pricing experience, resulting in a time delay of the price. The supervision department can lead the scientific research institutes to make scientific predictions on the market trend of that year, guide the market price and reduce the time delay. Using the central manifold theory, we obtained the direction of Hopf bifurcation, stability and existence of the bifurcated periodic orbit.

Using Pontryagin's maximum principle, we obtained the optimal cost strategy. According to the optimal cost strategy, guidances could be given to regulatory agencies' subsidies and tax adjustments to fisheries. This will only ensure the maximization of economic benefits and minimize the waste, but also realize the sustainable development of the population. We also got the stable dynamical variation of the delayed system under the control of guiding interventions. This showed that there would be less waste caused by hoarding for speculation if fisheries and fishery administrations could assist cultivation enterprises and individual farmers in finding a reasonable scale for farming scientifically, and publicize marketing strategies for aquatic products reasonably and moderately. And, from the simulation results, compared with that before adding the optimal control with guiding interventions, the waste caused by hoarding is indeed reduced on the premise of ensuring the stability of the system.

There are some topics remaining to refine our modeling and further the analysis. These include the influences of other species on the target species in polycultural mode and the influences of environmental complexity caused by algae, fungi and other microorganisms in water on population growth, as well as the global stability analysis and other bifurcation in the sense of mathematical theory.

Acknowledgments

This work was supported by the National Natural Science Foundation of China (61703083 and 61673100). The authors gratefully acknowledge the reviewers for their comments and suggestions, which have greatly improved the presentation of this work.

Conflict of interest

All authors declare that there is no conflicts of interest for this study.