Optimal control problem arises from illegal poaching of southern white rhino mathematical model

In this paper, a novel dynamical population model of a southern white rhino with legal and illegal poaching activity is introduced. The model constructed is based on a predator–prey model with southern white rhinos as prey and humans (hunters) as predators. We divide the southern white rhino population into three classes based on their horn condition. We investigate the existence and the stability of the equilibrium points, which depend on some threshold functions. From an analytical result, it is trivial that arresting as many hunters as possible helps conserve white rhinos, but it comes at a high cost. Therefore, an optimal strategy is needed. The optimal control is then constructed using Pontryagin’s minimum principle and solved numerically with an iterative forward–backward method. Optimal control simulations are given to provide additional insight into the dynamics of the model. Analysis of the cost function effectiveness is conducted using the ACER (Average Cost–Effectiveness Ratio) and ICER (Incremental Cost–Effectiveness Ratio) indicator method. The results show that the hunter population can be more easily controlled with a time-dependent hunter arrest rate rather than by treating it as a constant.

In order to increase the population of the rhino, a translocation strategy has been implemented [4]. A study about the impact of human intervention on different translocation strategies is discussed in [4,5]. Other efforts can also be applied to preserve white rhinos. One of the initiatives is to overcome poaching. It involves microchipping rhino horns so the government can track the horn and arrest the hunters [6]. An overpoaching of white rhinos could lead to an environmental change. Climate and environmental changes can also influence the dynamics of rhino population [7][8][9]. Because of overpoaching, there has not been any significant increase in the southern white rhino population since 2011 in Kruger National Park, South Africa [10].
According to [11], one way to learn about real-world problems is by constructing a mathematical model. The issue of rhino poaching, as described, can be compared to predatorprey interactions with white rhinos as prey and humans (hunters) as predators. Alfred J. Lotka in 1925 [12] and Vito Volterra in 1926 [13] proposed a simple model, known as the Lotka-Volterra model, to understand the interaction involving two species, namely one predator and one prey. Several previous studies have been constructed a more complex predator-prey mathematical model based on the Lotka-Volterra model in [14][15][16][17].
Optimal control theory is usually applied to control the spread of a disease in an epidemiology model, as proposed in [18,19]. It can also be used in a predator-prey model with infections among prey or predator populations, known as an eco-epidemiology model. For example, optimal control is used in [20,21] to minimize infected prey by harvesting or treating the infected population as the control variables.
The authors of [22] discuss the idea of optimal frequency to dehorning the rhinos. Their result suggests that dehorning should be done annually under some circumstances. The authors of [5] propose a mathematical model to describe a strategy to harvest rhinos such that the population can still survive. According to the above explanation and to the best of our knowledge, there are few mathematical models concerning the dynamics of southern white rhinos, especially since it was affected by illegal and legal poaching. Hence, unlike previously mentioned references, the aim of this article is to give an alternative perspective of the poaching problem of the southern white rhino, which, in this case, is from a mathematical point of view, especially with a dynamical system approach. This paper discusses how the illegal poaching could affect the dynamics of the southern white rhino. Intervention from the government, to imprison hunters and increase legal hunting, is considered in this model to counter overpoaching. Mathematical analysis on the model is conducted rigorously regarding the existence and bifurcation type of the model. An optimal control problem was analyzed to give the optional strategy against the hunters.
The rest of this paper is organized as follows. In Sect. 2, we present a four-dimensional system for interaction between white rhinos and hunters based on the Lotka-Volterra model. Intervention to overcome poaching by arresting hunters is also included in the model. In Sect. 3, we mainly discuss the existence and stability of the equilibrium points. An optimal control problem to design hunter arrest rates is formulated using Pontryagin's minimum principle in Sect. 4. The goal is to minimize the population of hunters and, at the same time, minimize the cost of implementing such means of intervention. In Sect. 5, we provide some numerical simulations based on the analytical results and the optimization algorithm. The paper ends with some conclusions and potential future directions in Sect. 6.

The construction of the mathematical model
A white rhino poaching model with intervention to control the number of hunters will be constructed carefully in this section. To begin our model construction, we separate the white rhino population into three classes based on their horn status. The first class is the juvenile rhino (denoted by N 1 (t)), where a rhino in this class has a small horn. We assume that hunters are not interested in hunting this kind of rhino due to its relatively small horn, which yields minimal profits. The second class is the adult rhino, which has a horn that is ready to be hunted (denoted by N 2 (t)). In our model, this class is the only one to be hunted for their horns. In some cases, hunters (mainly from legal hunting) do not kill the rhino. Instead, they only take the horn and let the rhino live. Therefore, it is important to introduce the last class of the rhino, which comprises dehorned rhinos. The adult rhinos in this class have been hunted for their horns but were allowed to live (denoted by N 3 (t)). Therefore, we have a total rhino population given by . The other population that will interact in our model is the hunter (human) population, which will act as the predator of rhinos, denoted by M(t).
Several assumptions that have been used in this article in the purpose of the model construction are given as follows: (A 1 ) Southern white rhinos live in groups. Each group may contain up to 14 rhinos. Therefore, we assume that the rhino population grows logistically to its carrying capacity K with a growth rate of r in the absence of hunters.
(A 2 ) Only adult rhinos can produce newborns, which in our case are the N 2 and N 3 populations.
(A 3 ) A transition from N 1 to N 2 population is due to the growth of the horn. Therefore, let α * be the transition rate from juvenile rhino N 1 to adult rhino N 2 .
(A 4 ) Interaction between hunter and adult rhinos is a form of predation. If β * is the success rate of hunters hunting rhinos per hunter, then β * N 2 M presents the predation of hunters on adult rhinos.
(A 5 ) As mentioned before, not all hunted rhinos will die. We assume that a fraction p of them will survive (pβ * N 2 M) after being hunted and transferred to the dehorned population N 3 , while the rest ((1p)β * N 2 M) will die.
(A 6 ) It takes approximately 3 years for adult rhinos N 3 to regrow their horns. Therefore, there is a transition (denoted by δ * ) from N 3 to N 2 due to the regrowth of the rhino horn.
(A 7 ) Let μ * be the natural death rate of rhinos in each class.
(A 8 ) The recruitment of new hunters comes from the profit of selling the horn from a dead rhino due to hunting activity, at a rate of ((1p)β * N 2 M). We assume that the recruitment rate of new hunters is proportional to the economic gain of this form of hunting. Hence, assuming γ * as the conversion parameter from dying rhinos due to hunting, recruitment of new hunters per unit time is given by γ * ((1p)β * N 2 M).
(A 9 ) Hunter population can decrease due to its natural drop-out rate due to increased awareness, causing the hunters to voluntarily stop hunting (denoted by ξ * ), and also due to being arrested by the government (denoted by h * ).
Based on the above assumptions and the interaction diagram in Fig. 1, the predatorprey model interaction between southern white rhinos and hunters is given by a four- Figure 1 Interaction diagram between white rhinos and hunters. A white rhino poaching model with intervention to control the number of hunters. The "in-coming" arrows indicate an increase in population, whereas "out-going" arrows indicate a decrease in population dimensional nonlinear ordinary differential equation as follows: which is supplemented by a nonnegative initial condition. Please note that all parameters in system (1) are assumed to be nonnegative.
A non-dimensionalization process is carried out using the transformations of variables and parameters in Table 1. The non-dimensional form of system (1) is used to simplify the analytical calculation for the existence and local stability analysis. Therefore, this form is only used in some subsections in Sect. 3.
Before we analyzed the model, we non-dimensionalized the parameters first, by substituting the transformations of parameters and variables in Table 1 to system (1). With this approach, we can reduce the number of parameters from ten to only seven. We obtain the non-dimensional form of system (1), where all variables and parameters involved have no units, which written as follows: Thus, it is trivial to show that the solution for each variable in system (2) will always be nonnegative in D for all time. So, system (2) is both mathematically and biologically wellposed. Please see [23,24] for further examples of the positively invariant region proof.

Model analysis
In this section, we will analyze the existence and local stability criteria of equilibrium points of system (2).

Equilibrium points
The equilibrium points of system (2) are obtained by making the right-hand side of system (2) equal to zero and solving the system with respect to each variable. We get three equilibrium points as follows: (E 1 ) The extinction equilibrium is given by This equilibrium describes the extinction of white rhino and hunter populations. (E 2 ) The hunter-free equilibrium is given by .
This equilibrium describes a condition when hunters disappear over time and only the white rhino population is left in the ecosystem. (E 3 ) The coexistence equilibrium is given by and y † is a positive root of the second degree polynomial and This equilibrium describes a condition when both populations exist in the ecosystem. Note that all of these equilibrium points will have a biological meaning if E 1 , E 2 , E 3 ≥ 0. It is easy to determine the conditions for the existence of E 1 and E 2 , while the existence of E 3 depends on the positive root of P(y). Therefore, based on the sum and product signs of the roots of P(y), E 3 has the following possibilities regarding its existence: 1. Unique coexistence equilibrium E 3 whenever R 2 > 1; 2. Two coexistence equilibrium E 3 whenever R 2 < 1, a 1 < 0 and a 2 1 -4a 2 a 0 > 0; 3. No coexistence equilibrium otherwise. If we consider the second case in the existence conditions of E 3 , where we have R 2 < 1, then which is positive since all the parameters used in the model are nonnegative. From (4), it follows that coefficient a 1 can be written in the form of a min < a 1 < a max , where Since a min and a max are positive, a 1 is always positive whenever R 2 < 1. Based on the above analysis, we can conclude that the second case will not happen if we use nonnegative parameters. Furthermore, from the discussion in this subsection, we have the following results.
Theorem 3.1 System (2) has three types of equilibrium. The extinction equilibrium E 1 will always exist, the hunter-free equilibrium E 2 can exist only when R 1 > 1, and the unique coexistence equilibrium E 3 exists whenever R 2 > 1.
From the relation between R 1 and R 2 in equation (3), and Theorem 3.1, we have the following corollary.

Stability of the equilibrium points
Stability analysis of the equilibrium points is conducted using the linearization technique with the Jacobian matrix approach and then evaluating the eigenvalues of the Jacobian matrix, which are observed in the equilibrium points. The equilibrium is stable when the real parts of all the eigenvalues are negative [25]. The Jacobian matrix of system (2) is given by

Theorem 3.3
The extinction equilibrium E 1 of system (2) is locally asymptotically stable if R 1 < 1.
Note that from Theorem 3.1 we have that E 2 exists only when R 1 > 1 which contradicts the stability of E 1 . Therefore, we can conclude that the existence of E 2 causes E 1 to become unstable, and vice versa. This result is stated in the following corollary.

Stability of the hunter-free equilibrium E 2
Next, we analyze the local stability of E 2 . The Jacobian matrix around E 2 takes the form where The other two eigenvalues are given by the roots of the second degree polynomial: Polynomial P 2 (λ) will have 2 negative roots if R 1 > 1. However, based on Theorem 3.1, E 2 exists whenever R 1 > 1. Therefore, λ 3 and λ 4 are always negative. From this analysis, the result is stated in the following theorem.

Theorem 3.5
The hunter-free equilibrium E 2 of system (2) is locally asymptotically stable if R 2 < 1.
Similar to the relation between E 1 and E 2 , from Theorem 3.1, the existence of E 3 causes E 2 to become unstable, and vice versa. This result is stated in the following corollary.

Corollary 3.6
Whenever E 1 is stable, then E 2 and E 3 do not exist. Furthermore, the hunterfree equilibrium E 2 is stable if and only if the coexistence equilibrium E 3 does not exist.

Stability of the coexistence equilibrium E 3
Different from E 1 and E 2 , the stability analysis for E 3 becomes very complex and is not tractable directly using the previous approach. Therefore, the center manifold theory in [26] will be used to obtain the local asymptotical behavior of E 3 near R 2 = 1.
To conduct the local stability of E 3 , we will show the existence of a single zero eigenvalue of the linearization of system (2) around E 2 . First, let h be the bifurcation parameter, x 1 = y 1 , x 2 = y 2 , x 3 = y 3 , y = y 4 , and denote Y = (y 1 , y 2 , y 3 , y 4 ) T . System (2) can be written in the form of dY dt = (f 1 , f 2 , f 3 , f 4 ) T as follows: The linearization around E 2 (see (5)) evaluated at h =h is The characteristic polynomial of matrix A is given by It is clear that zero is a simple eigenvalue of A. Hence, we can continue to analyze the local stability of E 3 of system (6) for h nearh.
To compute a right eigenvector w, we consider the system Aw = 0. We obtained that a right eigenvector associated with zero eigenvalue is w = (w 1 , w 2 , w 3 , w 4 ) T , where To simplify the expressions for the components, let w 3 = αμp(α + μ) 2 (R 1 -1). Based on the existence condition of E 2 , which is R 1 > 1, we have that w 3 > 0. So, a right eigenvector (7) becomes It can be shown that w 2 < 0 and w 4 > 0. The component w 2 is negative but this is acceptable, since it corresponds to the positive second entry of E 2 . Furthermore, we compute the left eigenvector v by solving vA = 0. The left eigenvector which is associated with zero eigenvalue is given by Next, we will compute a and b in Theorem 4.1 of [26].
Since v 1 = v 2 = v 3 = 0, we only need to compute the partial derivative of f 4 . For system (6), the associated nonzero partial derivatives of f 4 at E 2 are as follows: , The other second derivatives appearing in the formula for a and b are all zero. From the right eigenvector (8), left eigenvector, and result (9), algebraic calculations show that: Therefore, based on item (ii) of Theorem 4.1 in [26], we can conclude that we have a change of stability between E 2 and E 3 when R 2 = 1. To be precise, E 3 is stable if R 2 > 1, while E 2 becomes unstable. This result is stated in Theorem 3.7.
Theorem 3.7 The coexistence equilibrium E 3 of system (2) is locally asymptotically stable whenever it exists, that is, when R 2 > 1. (3), if R 1 < 1, then R 2 < 1. In other words, R 1 < 1 < R 2 is impossible. It implies that the extinction equilibrium E 1 (stable if R 1 < 1) and the coexistence equilibrium E 3 (stable if R 2 > 1) cannot be stable at the same time.

Numerical example on bifurcation diagram of system (2)
Here in this subsection, we show the existence of a forward bifurcation based on previous results in Theorems 3.5 and 3.7 as shown in Fig. 2, through a numerical example by creating a bifurcation diagram around R 2 = 1. To draw the bifurcation diagram (the graph of y † as a function of p), we substitute K = 1000, r = 1 32 , β * = 0.8 1000 , δ * = 1 3 , γ * = 0.3, μ * = 1 40 , ξ * = 1 65 , h * = 0, α * = 1 3 , into the non-dimensional version of the system. We choose p as the bifurcation parameter.
As it can be seen in Fig. 2(a), using γ * = 0.3, we have that R 2 = 1 when p = 0.507. Starting from p = 0, when p is increasing, y in E 3 exists and is monotonically decreasing until it reaches p = 0.507. When y exists, at the same time, it is also stable, but y = 0 in E 2 is unstable. When p = 0.507, a change of stability occurs. We found that y disappeared when p leaves p = 0.507 and increases; at the same time, y = 0 in E 2 becomes stable. On the other hand, in Fig. 2(b), we use γ * = 0.9 to present a larger conversion of poaching profits into the new recruitment of new hunters. When γ * = 0.9, we found that R 2 = 1 when p = 0.84. It means that a larger γ * can enlarge the interval of the existence of y in E 3 . It can also be seen that when γ * = 0.9, y becomes non-monotonically decreasing with respect to p. As a reminder, p presents the proportion of dehorned surviving adult rhinos. This result tells us that enlarging γ * (highly related to the high price of horns in the market, low economic level of humans which increase the desire for hunting, etc.) will change the interpretation of p to the coexistence size of y in E 3 . From Fig. 2(b), starting from p = 0, when p increases, the size of y in E 3 increases until it reaches a specific point of p, before decreasing to zero as p approaches 0.84.
Remark When the desire for hunting increases, which may be caused by high prices for horns, a greater survival rate for rhinos dehorned by legal hunting will increase the number  Figure (a) is presented using γ * = 0.3, while (b) is for γ * = 0.9. Red and blue curves present E 3 and E 2 , respectively. The solid curve presents the stable equilibrium, while the dotted curve is for the unstable equilibrium of hunters in the equilibrium level in p ∈ [0, p critical ), and decrease when p ∈ [p critical , 1], where p critical is p such that ∂y ∂p = 0. Therefore, when the price of horns is high enough, increasing legal hunting can be counterproductive on the number of hunters in the field.

Summary of the existence and local stability for system (1)
Until this stage, we have obtained the results of existence and local stability of all equilibrium points for system (2). Based on Theorem 3.1, the existence of equilibrium points depend on R 1 and R 2 . From Theorem 3.3, E 1 , which describes the extinction of white rhino and hunter populations, is stable if R 1 < 1. From Theorem 3.5, E 2 , which describes the extinction of hunter populations, will exist if R 1 > 1 and stable if R 2 < 1. From Theorem 3.7, E 3 , which describes both populations existing in ecosystem, will exist and be stable if R 2 > 1.
Note that the parameters of systems (1) and (2) have the same description. However, system (1) is more easily interpreted than system (2). Therefore, in this subsection, we will consider system (1). The dimensional forms of R 1 and R 2 are denoted by R * 1 and R * 2 , respectively, and are given by: Each threshold (R 1 , R 2 ) for system (2) is being substituted with transformation in Table 1 to get the dimensional form of these thresholds, which are related to the original model (1). The summary of the existence and local stability conditions of system (1) is given in Table 2 and will be interpreted in this subsection.
Based on Table 2, the extinction equilibrium E 1 condition tells us that if R * 1 is less than 1, then the system cannot sustain the southern white rhino population, which also has an impact on the extinction of the hunter population. From (10a), we can see that the value of R * 1 depends on the transition rate from juvenile to adult rhino (α * ), intrinsic growth rate (r), and the natural death rate (μ * ). Furthermore, since ∂R * 1 ∂α * > 0, ∂R * 1 ∂r > 0, and ∂R * 1 ∂μ * < 0, we can conclude that the slower the juvenile white rhinos mature into adults (α * ), the slower the growth rate of white rhinos (r), and the faster the white rhinos die (μ * ), the easier it becomes for both populations to become extinct (E 1 is locally stable).
Ideally, we want a hunter-free equilibrium E 2 in the real world. We need E 2 to exist and be locally stable. Based on Table 2, to achieve this, the condition R * 2 < 1 < R * 1 must be satisfied, where hunters disappear over time and only white rhinos remain in the ecosystem.
Based on Table 2, the coexistence equilibrium E 3 condition tells us that if R 2 is greater than 1, then hunters remain and rhino poaching still occurs. From (10b), the value of R * 2 depends on almost all parameters (see Table 1) except for the transition rate due to the Table 2 Summary of the existence and local stability conditions of each equilibrium points of system (1). The results are based on those of system (2) as stated in Theorems 3.1, 3.3, 3.5, and 3.7 Extinction equilibrium E 1 = (0, 0, 0, 0) Hunter-free equilibrium E 2 = (N 1 , N 2 , 0, 0) Table 3 Sensitivity of R * 2 . A positive sign (+) means that the curve of the parameter toward R * 2 is monotonically increasing, so a larger parameter value will lead to greater R * 2 . Meanwhile, a negative sign (-) means the exact opposite regrowth of rhino horns (δ * ). For example, from ∂R * 2 ∂α * > 0, we know that whenever α * increases, R * 2 also increases. On the other hand, since ∂R * 2 ∂h * < 0, we know that as h * increases, R * 2 decreases. So, the faster the juveniles mature (α * ) and the fewer hunters are arrested by the government (h * ), the easier it is for both populations to exist. Table 3 gives the rest of the sensitivity results of R * 2 . From the sensitivity results of R * 1 and R * 2 , a similar interpretation can be applied to determine the conditions when R * 2 < 1 < R * 1 is fulfilled, so that the hunter-free equilibrium E 2 of system (1) is locally stable. If R * 2 < 1 < R * 1 is not fulfilled, then either the extinction equilibrium E 1 or the coexistence equilibrium E 3 will be locally stable. For an additional insight of these results, autonomous simulations are conducted in Sect. 5.

Optimal control problem
One of the obstacles in arresting hunters by microchipping rhino horns or increasing the number of forest rangers is the high cost of such means of intervention. Therefore, to overcome the problem, system (1) can be developed into an optimal control problem. The intervention parameter h * , which was previously set as a constant, now changes into a control variable that depends on time, denoted by h * (t). The purpose is to minimize the population of hunters (y(t)) and also minimize the cost incurred from implementing such intervention (h * (t)).
First, assume that the incurred cost from hunters (ω M ) is directly proportional to the numbers of hunters M(t), so the cost can be represented as a linear function. However, since a wider area of intervention leads to a nonlinear cost increase, the cost (ω h ) of intervention h * (t) is represented by a quadratic function. Therefore, we consider the objective function of system (1) with arresting hunters as follows: Cost due to the number of hunters Cost for arresting the hunters dt, where ω M > 0, ω h > 0 are weighted constants, and t f is the final time of simulation.
Next, we construct the Hamiltonian function H = g + z · f based on Pontryagin's minimum principle [27]. The Hamiltonian consists of the sum of the integrand function (g) of the objective function (11) and the inner products of the state system (f) in system (1) with the adjoint variables z k (t), k = 1, 2, 3, 4. The Hamiltonian function (H) can be stated as follows: Differentiating the Hamiltonian function H with respect to the state variables yields the adjoint system given by: with the transversality condition for adjoint variables z k (t f ) = 0, k = 1, 2, 3, 4.
To obtain the optimal conditions, we differentiate the Hamiltonian (12) with respect to the control variable and set the equation to zero: Solving (14) with respect to the control variable, we obtain h * (t) = z 4 (t)M(t) 2ω h , which must satisfy h min ≤ h * (t) ≤ h max for all t ∈ [t 0 , t f ]. Since the admissible control variable should be bounded by h min and h max , the optimal control variableĥ(t) can be written as:ĥ where h min and h max are the lower and upper bounds of allowed intervention, which are related to the limitations of budget, resources, etc. The computations of the optimal control and state values are performed using a Runge-Kutta method. The algorithm used is the iterative forward-backward method, which is summarized as follows: An initial estimate for the control h * should be determined. The next steps require the solution of a system of ordinary differential equations by any suitable numerical scheme. The state variables are solved forward in time using system (1), integrating from t 0 to t f . The results obtained for the state variables are plugged into the adjoint equations (13). Solving these adjoint equations requires backward integration from t f down to t 0 , which is supplemented with the final time condition of the adjoint. Both the state and adjoint values are then used to update the control (15). By applying this updated control value to replace the initial guess of control, the method afterwards is repeated until the current value of the control variable converges [27].

Numerical simulations
In this section, several scenarios are implemented in numerical simulations to confirm and illustrate the results from the previous sections.

Autonomous simulations
To understand the role of the model parameters to R * 1 , R * 2 , and the dynamics of system (1), we perform autonomous simulations to illustrate the results given in Sect. 3. Note that in this section, we only focus on the long-term effects to see the stability tendency of system (1).
The first analysis is conducted to determine the effects of the transition rate from juvenile to adult rhinos (α * ) to the dynamics of system (1). We use three different variations of α * , namely 1 11 , 1 7 , and 1 3 . A set of the other parameters are given as follows: From (10a) and (10b), R * 1 and R * 2 depend on parameter α * . From Table 2, the value of R * 1 determines the stability of E 1 and the existence of E 2 , while E 2 determines the stability of E 2 and E 3 . For the set of parameters in these simulations, the result on the stable equilibrium points of system (1) is given in Table 4. It can be seen in Table 4 that for the smallest value of α * , the extinction equilibrium E 1 will be stable. For the biggest value of α * , the coexistence equilibrium E 3 is the stable one. This confirms the analysis in Sect. 3.4 of the sensitivity results of R * 1 and R * 2 towards α * . To be precise, the condition of α * to guarantee the existence and stability of hunter-free equilibrium E 2 is given by α min < α * < α max , where In our numerical experiment, we get that α min = 0.1 and α max = 0.1835. It means that, for the set of parameters in this simulation, the transition from a juvenile white rhino to an adult white rhino needs to take at least 5.45 years and no more than 10 years for system (1) for the hunter-free equilibrium E 2 to be stable.

Figure 3
The effects of the changes in parameter α * . The blue curve (stable to E 1 ) represents that it takes 11 years for white rhino to become adult, the black curve (stable to E 2 ) represents that it takes 7 years, while the red curve (stable to E 3 ) represents that it takes 3 years The illustrations for these autonomous simulations are given in Figs. 3 and 4. Figure 3 shows the results for each compartment with the initial values (N 1 , N 2 , N 3 , M) = (30, 70, 0, 5) and shows that for each value of α * , system (1) will approach different equilibrium points. From Fig. 3(a), we can see that in the beginning, the slower rate of the transition from juvenile to adult is causing a higher number of juvenile white rhinos. However, the reproduction of white rhinos depends only on adult rhinos. Thus, for a longer time, the slowest rate (represented by the blue curve) prevents white rhinos from reproducing. This is why the blue curve is stable to E 1 .
Further, Fig. 4 shows the dynamics of adult white rhino (N 2 ), dehorned white rhinos (N 3 ), and hunters (M) populations using two different initial values, (N 2 , N 3 , M) = (70, 0, 5) and (140, 1, 10). This figure tells us that if all parameters remain the same but have different initial values, system (1) does not lead to different equilibrium points.
Next, we conduct a simulation to determine how the numbers of hunted white rhinos that survived (the changes in parameter p) and hunters arrested by the government (the changes in parameter h * ) affect the dynamics of system (1). From (10a), parameters p and h * do not change the value of R * 1 . So, we will see how R * 2 can be determined by relying on p and h * , qualitatively. A set of the other parameters are given as follows: such that R * 1 = 2.1276 > 1. Based on Table 2, for any value of p ∈ [0, 1] and h * with this set of parameters, system (1) will never reach a stable extinction equilibrium E 1 . When we substitute all parameter values into R * 2 = 1, we have which is illustrated in Fig. 5. From Fig. 5, when p > 0.8579, the population always achieves a hunter-free situation (R * 2 < 1), and thus it is unnecessary to arrest hunters. When p < 0.8579, we must arrest hunters to achieve the condition R * 2 < 1, with the minimum rate of h * is given by h min = 0.0929 -0.1083p. This condition indicates that there is a possibility that arresting hunters is not needed to eliminate illegal hunting. It is done by stocking up on horns. It means that we need to cut the horns and let the rhino live (higher values of p). Then, the number of dehorned rhinos (N 3 ) will increase, which reduces the chance of illegal hunters hunting and killing adult white rhino with horns (N 2 ). Note that, when the price of horns is high enough, increasing legal hunting can be counterproductive on the number of hunters in the field as explained in Sect. 3.3.
Finally, we compare the effects of different values of p and h * to each compartment in system (1) with the values from Fig. 5 taken into consideration. The result of this autonomous simulation is given in Fig. 6 with the initial values (N 1 , N 2 , N 3 , M) = (100, 600, 0, 10).
From Fig. 6, the red curve with p = 0.2, h * = 0 yields R * 2 = 1.7725 > 1 and, according to Table 2, system (1) will reach a stable coexistence equilibrium E 3 . For the blue curve, we increase the value of p, which becomes p = 0.9, h * = 0. The blue curve is a stable hunterfree equilibrium E 2 with R * 2 = 0.8175 < 1. Meanwhile, for the black curve, we increase the value of h * , which becomes p = 0.2, h * = 0.4, yielding R * 2 = 0.3319 < 1 and a stable E 2 . To make it clearer, we compare the red and blue curves (the changes in parameter p) in Fig. 7, while the red and black curves (the changes in parameter h * ) are compared in Fig. 8.  Table 3 about the sensitivity of R * 2 towards p and h * , which is that as p and h * increase, R * 2 decreases.  Table 3 that larger p and h * will reduce R 2 . The dashed red line in (b) is the critical line when R * 2 = 1 based on (16). The red area defines the coexistence condition, while the blue area defines an ecosystem with no hunters Further, the condition of p to guarantee the existence and stability of hunter-free equilibrium E 2 is given by In our numerical experiment (when we set h * = 0), we get that p > 0.8579. It means that, for the set of parameters in this simulation, the hunter needs to let as much as 85.79% of the total of adult white rhinos that were successfully hunted live, for system (1) to tend to a hunter-free equilibrium E 2 . Meanwhile, the condition of h * to guarantee the existence and stability of hunter-free E 2 is given by

Figure 6
The effects of the changes in parameters p and h * . The red curve (stable to E 3 ) represents that the proportion of hunted rhino being left alive is p = 0.2 with no intervention or no hunters being arrested (h = 0), the blue curve (stable to E 2 ) represents the proportion of p = 0.9 and h = 0, and the black curve (stable to E 2 ) indicates that the proportion is p = 0.2 and the intervention is carried out with constant rate h * = 0.4

Figure 7
Autonomous simulation of p. The dynamic of N 2 , N 3 , and M of system (1) with two different initial values, (N 2 , N 3 , M) = (600, 0, 10) and (400, 10, 5). For both initial values, the blue curve (with p = 0.9) tends to hunter-free equilibrium E 2 = (451.06, 0, 0), and the red curve (with p = 0.2) tends to coexistence equilibrium E 3 = (80.13, 1. 27, 35.59) In our numerical experiment (when we set p = 0.2), we get that h * > 0.0712. It means that, for the set of parameters in this simulation, we need to use intervention to arrest the hunters with minimum constant rate h * = 0.0712, for system (1) to tend to a hunter-free equilibrium E 2 .

Optimal control simulations
Here, the optimal control simulations are conducted to determine an optimal hunter arrest rate as a time-dependent variable which is based on the results and algorithm in Sect. 4. We choose the interval value of intervention (h * (t)) to be 0 ≤ h * (t) ≤ 1 with the time interval of 0 ≤ t ≤ 10 years in these simulations. To find a balance between each component in the objective function (11), we choose the weighted parameter values used on the objective function to be ω M = 10 and ω h = 1.  First, for the purpose of comparing the reduction of hunters and killed rhinos, the simulations are conducted without control (h * = 0), with constant control (h * ), and with time-dependent control (h * (t)). The initial values of each compartment are N 1 (0) = 100, N 2 (0) = 600, N 3 (0) = 0, M(0) = 10, with a set of parameters given as follows: Several reduction results are given in Table 5 with the illustrations in Fig. 9. We have chosen the constant control h * = 0.6908 which is equivalent to the average of time-dependent control h * (t) for these simulations.
From Table 5, a large number of hunters without intervention causes the most cost, that is, 1668.2890. The result of the intervention carried out as a time-dependent variable successfully reduced the average total hunter population by 93.29%. Meanwhile, the constant control reduced the average total hunter population by 89.90%. In addition, the cost of time-dependent control (117.8970) is less than that of constant control (173.1947).
The dynamics of the time-dependent (h * (t)) control is shown in the blue curve of Fig. 9(f ). Some form of intervention is required if hunters are present. In this simulation, for about 5 years, this intervention is carried out in the upper bound rate which is higher than a constant rate (red curve of Fig. 9(f )). However, when few hunters remain, the in- Figure 9 Control optimal simulation. The dynamic of each compartments, killed white rhino, and the intervention. The black curve represents no intervention, h * = 0, the red curve represents the intervention is carried out with constant rate, h * = 0.6908, and the blue curve represents the time-dependent intervention, h * (t) tervention gradually decreasing over time to the lower bound at the end of the simulation period.
In both constant and time-dependent control cases, the decrease in the number of hunters ( Fig. 9(e)) reduces the risk of white rhinos being hunted and killed ( Fig. 9(d)), thereby leading to fewer dehorned white rhinos (Fig. 9(c)) and eventually resulting in a white rhino population which is dominated by adult rhinos that still have horns ( Fig. 9(b)). From both Table 5 and Fig. 9, we can conclude that time-dependent control is more effective at reducing costs and hunter population than constant control. So, starting from this, we only consider the time-dependent control (h * (t)).
Next, we investigate the effect of varying ω M and ω h . These coefficients are the balancing cost factors or weights on the costs caused by hunters M and the costs of implementing the controls h * (t). Weighted parameter ω h represents the cost of buying weapons or providing security personnel, while ω M is the cost due to the existence of hunters such as media campaigns, efforts to strengthen the law, etc. We consider the following combination of weighted parameters making up three control strategies: 1. Strategy A: the baseline, ω M = 10 and ω h = 1. The results of the optimal control simulations using three strategies are shown in Fig. 10. We only show the effect of the given controls (Fig. 10(c)) on dead white rhinos ( Fig. 10(a)) and the illegal hunter populations (Fig. 10(b)). The solutions for other compartments in system (1) are similar to the ones portrayed in the blue curves of Fig. 9.
Agusto and Leite in [23] used cost-effectiveness analysis to determine the vaccination and personal protection control strategies and the benefits gained from implementing these controls in the meningitis model. Here, we will using the same approach to determine the most cost-effective strategy to reduce the number of illegal hunters and dead rhinos.
The cost-effectiveness analysis implemented in [23] uses three approaches, namely, the infection averted ratio (IAR), the average cost-effectiveness ratio (ACER), and the incre-  mental cost-effectiveness ratio (ICER). However, since this model in system (1) does not have a recovered compartment, we will not take IAR into consideration. The ACER is calculated against the non-intervention scenarios. The total cost produced by the intervention is estimated using the objective function given in (11). Instead of total number of infections averted, we use the total number of rhinos that were prevented from being killed. It means the difference between the total of white rhinos killed in the absence of control compared to that with time-dependent control. With some adjustments, ACER formula in Sect. 4 of [23] becomes ACER = Total cost produced by the intervention Total number of rhinos that have been prevented from being killed .
The ACER for each strategy is determined and the results are given in Table 6. According to Table 6, the cheapest way will be to employ Strategy C, while the best policy in terms of preventing the loss of white rhinos is Strategy B. The most cost-effective strategy, with the smallest ACER, will be Strategy C, followed by Strategy B, and the least cost-effective strategy is Strategy A.
For more clarity, we can also calculate the ICER of all strategies. The ICER is the additional cost per additional health outcome. In this model, health outcome is defined by the number of rhinos that avoided being killed. To compare two or more competing intervention strategies incrementally, one intervention is compared with the next-least-effective alternative [23]. Thus, the ICER formula is given by ICER = Difference costs in strategies i and j Difference in total number of rhino deaths averted in strategies i and j .
The ICER calculation start from the strategy which given the least number of averted killed white rhinos. Using these simulation results in Table 6, we have that Strategy C prevented the fewest rhino deaths, followed by Strategy A, and Strategy B which prevented the most rhino deaths. The ICER is computed as follows: The comparison (18b) shows a cost saving of 1058 for Strategy B over Strategy C. Using the same reasoning, Strategy B will be more costly and less effective to implement compared to Strategy C. Following this analysis, we conclude that Strategy C (cheaper costs caused by hunters) is the most cost-effective strategy.
Repeating the entire process, we can determine the next most cost-effective strategy between Strategies A and B. The ICER in (17c) is obtained by comparing Strategies A and B. From the negative ICER value, Strategy B strongly dominates Strategy A. This analysis implies that Strategy B (cheaper costs of implementing the control) is the next most costeffective strategy after Strategy C.

Conclusions
This paper proposed a white rhino poaching model with intervention to control the number of hunters in the system (1) and obtained a non-dimensional model in system (2). Firstly, the existence and stability of the equilibrium were discussed. System (1) has extinction (E 1 ), hunter-free (E 2 ), and coexistence (E 3 ) equilibria, along with their local stability conditions. We do not want any illegal hunting to happen, so we need system (1) to have a stable hunter-free equilibrium E 2 . The results show that to guarantee the existence and local stability of E 2 , the condition R * 2 < 1 < R * 1 must be achieved. The analytical result regarding the dependency of the final size of the hunter population in the coexistence equilibrium with proportion parameter (p), could trigger a counterproductive effect on controlling the hunter population. We find that it is possible for any increase in this proportion to increase the number of hunters, whenever the lure of being a hunter is also high. This trend is presented by γ * (the conversion parameter), which can be related to the high price of horns, lack of security, high demand of horns, etc.
To reduce the cost caused by the hunter population and intervention, system (1) was developed into an optimal control problem. The hunter arrest rate parameter, which was previously set as constant, is changed into control variables that depended on time. An optimal control problem to investigate hunter arrest rates has been formulated in the objective function (11). We have derived a Hamiltonian function to solve the optimization problem numerically. Our optimal control simulation states that the time-dependent control is much more effective rather than constant control. To give a better understanding on the effect of weight parameters in the objective function, we performed three simulation scenarios regarding the composition of the weight parameters and analyzed them using the concept of ACER and ICER. We found that reducing the weighted cost which presents a cost that is related to the high number of hunters will give a better result in controlling the number of hunters in the field. This cost is probably related to media campaign on the dangers of illegal hunting, effort to strengthen the law, etc.
To sum up our results, intervention by the government must consider conditions in the field such as intervention costs, sale prices for rhino horns, hunter populations, etc., so that the expenses incurred are optimal. Controlling the number of hunters as a time-dependent variable gave lower costs and succeeded in reducing the total number of hunters better than the constant control or without any control at all. Based on numerical simulations, when hunters are present, some form of intervention is required. However, when the number of hunters begins to decrease, the rate of intervention must be reduced. All of our optimal control simulations suggest exerting maximum effort to imprison the hunters since the beginning of simulation, and starting to decrease it when the number of hunters is already small enough.
Despite all contributions of our article to the problem of illegal poaching in the southern white rhino population, some limitations need to be stated. One limitation is that the continuous time-dependent control in our article is hard to be implemented directly in the field. To implement the intervention directly in the field, it is easier to do it as a constant effort for a specific time interval. Several authors have been implementing this in epidemiological models [19,28]. Another limitation is that we do not consider the government's possibility to open up the privatization of rhino hunting. In some cases, this private conservation may permit rhinos to be hunted if they meet certain conditions, but only to cut rhino horns, not to kill the rhinos. This consideration can be accommodated by adding variables in the model. Besides, the influence of media campaigns, the availability of horns on the market, and fluctuations in the price of horns are very influential on the intensity of hunting. Therefore, awareness-based models can be considered as modification of this article.