1. INTRODUCTION

The Bali Starling (Leucopsar rothschildi) is a rare and endemic species of Java and Bali. Its last natural habitat is in the West Bali National Park (WBNP) [3]. One of the main breeding centers is located next to the WBNP office, called the Tegal Bunder Wildlife Sanctuary Unit (US) at Jembrana, Bali. This breeding and release area is considered one of the original habitats of the Bali Starling. The success in increasing the wild population of the Bali Starling is largely due to the synergistic management of ex-situ (outside the habitat) and in-situ (inside the habitat) conservation efforts. Community enthusiasm for ex-situ breeding raises awareness and fosters a love for the conservation of the Bali Starling. The Bali Starling, an endangered species, faces numerous threats including habitat loss, capture for the pet trade, and forest fires ([31], [37], [30]). Effective population management is crucial for its survival, with studies highlighting the need for proper reserve management understanding of its habitat [37], and vulnerability to forest fires [30]. The species' captive breeding and conservation program has been a key focus, with assessments of its extinction risk and genetic management [12]. These studies collectively underscore the importance of a comprehensive approach to population management, including habitat protection, captive breeding, and conservation programs. We conducted this study to help determine the best carrying capacity for the Bali Starling breeding system so that the population can be sustained in the long term. Previous research, such as by Noerdjito [28] and Pramatana [31], shows that both habitat and population management play a big role in the success of these efforts. Noerdjito highlights the importance of restoring natural habitats, while Pramatana focuses on how artificial treatments can affect the birds' movement and distribution. Widyaputra [39] examine the production process management and appearance of Bali Starlings in captivity.

The Bali Starling (Leucopsar rothschildi) is critically endangered and endemic to Java and Bali. Its survival is threatened by habitat loss and illegal poaching. Immediate and effective conservation strategies are crucial to

^∗Corresponding Author

Received February 25th, 2025, Revised May 14th, 2025 (first), July 28th, 2025 (second), Accepted for publication September 10th, 2025. Copyright ©2025 Published by Indonesian Biomathematical Society, e-ISSN: 2549-2896, DOI:10.5614/cbms.2025.8.2.3

prevent extinction. While traditional enforcement approaches have been ineffective, initiatives embracing local context and involving diverse stakeholders have shown promise [19]. This research addresses the urgent need to enhance and optimize breeding programs, ensuring a sustainable increase in the Bali Starling population both in captivity and in the wild. By developing and applying an optimal control model, the research aims to improve the effectiveness of current conservation efforts, making them more efficient and impactful. Furthermore, Efficient allocation of limited resources is vital for the success of any conservation program. The breeding centers and habituation processes for the Bali Starling require substantial financial, human, and material resources [37]. This research is important as it aims to determine the optimal carrying capacity and cost-effective management strategies for the breeding system at US Tegal Bunder, WBNP. By minimizing costs and maximizing the population growth of the Bali Starling, the research ensures that conservation resources are used sustainably, allowing for long-term success and the possibility of replicating these strategies for other endangered species.

Mathematical models have played a critical role in understanding ecological and biological population dynamics, including the management of endangered species [23]. Numerous studies have developed deterministic and stochastic models to explore species survival, reproduction, and mortality in both natural and controlled environments [5], [27]. Optimal control theory, especially Pontryagin's Maximum Principle, has been widely applied to manage harvesting strategies, epidemic control, and resource allocation in ecological systems [24]. However, its application in the context of captive breeding programs for critically endangered species, such as the Bali Starling (Leucopsar rothschildi), remains limited.

Recent work by [40] and [4] have demonstrated how cost-sensitive control interventions can balance ecological goals and economic limitations, but most of these focus on general wildlife or prey–predator dynamics. Few studies explicitly integrate the captive-exposure-release phases of population management into a unified control framework. Moreover, the influence of stochastic environmental variations on captive population control is rarely incorporated. This study advances the current state of the art by developing a nonlinear dynamic system tailored to the Bali Starling breeding program and embedding a cost-optimized control policy that accounts for practical conservation constraints. The inclusion of stochastic perturbations further enhances the model's applicability to real-world ecological uncertainties.

This article presents several original contributions to the mathematical modeling of conservation systemsit introduces a novel compartmental model for the Bali Starling (Leucopsar rothschildi) that captures the key biological processes across wild, captive, and juvenile exposure phases. Furthermore, realistic simulations informed by field data from US Tegal Bunder, WBNP, are used to evaluate policy implications, offering a practical tool for decision-makers. These contributions provide a new bridge between mathematical theory and real-world conservation strategies for endangered species.

There are two main objectives of the paper, firstly is to present a dynamic mathematical model to describe the population dynamics of the Bali Starling (Leucopsar rothschildi) in the breeding center at US Tegal Bunder, WBNP. The model incorporated the growth, transfer, and habituation processes of the population, and were used to design an optimal control strategy. This strategy is to ensure the effective management of the population by balancing the breeding efforts, minimizing costs, and maintaining a sustainable population level. Secondly is to determine the actual carrying capacity and optimal feeding strategies within the Bali Starling breeding system at US Tegal Bunder, WBNP. By conducting simulations and stability analyses, the study will identify the maximum manageable population size and the minimum cost required for maintaining a healthy and stable population. This evaluation will provide valuable insights into the efficient allocation of resources and the potential for scaling up the breeding efforts for the conservation of the Bali Starling.

A range of mathematical models have been developed for bird population and management. Sunday [38] presents a model for bird harvesting in intensive poultry systems, while Yaegashi [42] focus on population management of fish-eating birds and vertebrate pest species, respectively. These models consider factors such as cost, human interventions, and the impact of management actions on bird populations. Ram´ırez [25] propose age-structured population growth models and simple population models, respectively, for specific bird species. Optimal control theory has been applied to various aspects of bird population management. Yaegashi [42] explore the use of stochastic and differential equations models, respectively, for sustainable population management. Guiver [16] and Guiver [17] introduce a novel management methodology for restocking declining populations, using integral control. These studies collectively highlight the potential of optimal control theory in addressing the complex challenges of bird population management.

Previous studies on the Bali Starling have highlighted the critical status of the species and the need for conservation efforts [3]; [36]; [18]). These efforts have included the use of customary village regulations (awig-awig) in Nusa Penida Islands [36] and the study of the bird's behavior and food sources [35]. The success of these conservation measures is evident in the increasing population of the Bali Starling in Nusa Penida Islands [36].

The conservation of endangered species like the Bali Starling can benefit from optimal control strategies and population modeling. Mathematical models have been developed to analyze population dynamics and breeding programs [14]. These approaches can enhance population growth and achieve higher population sizes compared to uncontrolled scenarios. Local community involvement and cultural practices have also played crucial roles in conservation efforts [18]. These studies collectively demonstrate the potential benefits of implementing control measures and considering diverse factors in the conservation of endangered species like the Bali Starling.

MODEL FORMULATION 2.

2.1. Dynamical System of Bali Starling Population

To model the dynamics of Bali Starling (Leucopsar rothschildi) population, we consider a compartmental system consisting of three subpopulations these are M(t) is the total population in the wild, P(t) is the population in captivity, and E(t) is the juvenile population in the habituation phase (exposure) as seen in Figure 1.

The diagram visualizes the change in the wild population M(t) over time. Term \(r\left(1-\frac{M}{K}\right)M\) represents the natural reproduction of the Bali Starling in the wild. It follows a logistic growth model, where r is the intrinsic growth rate, K is the environmental carrying capacity, the maximum number of individuals the habitat can support. As M approaches K, the growth rate slows due to limited resources, simulating densitydependent regulation. The natural mortality process \(-\delta M\) is the natural death rate or loss of individuals from the wild population due to aging, disease, predation, or other ecological pressures. The parameter \(\delta\)represents the mortality rate. The juvenile reintroduction process, \(+\gamma E\) captures (positive term) the return flow of individuals from the juvenile or exposure class E(t) back into the wild. Parameter \(\gamma\) quantifies the rate at which juveniles are successfully released or reintegrated into the wild population after a habituation period.

The captive population dynamics P, the arrow from \(M \to P\), indicates transfers from the wild to captivity, influenced by \(dP/dt = \alpha M - (\beta + \mu)P\), where \(\alpha\) is the percentage of wild birds transferred into captivity, \(\beta\)is the maturation rate of chicks from captivity into juvenile class, and \(\mu\) is the natural death rate. The juvenile (exposure) population (E(t)) is produced in captivity and move back to the wild \(dE/dt = \beta P - (\gamma + \mu)E\), where \(\beta P\) is the chicks maturing in captivity, \(\gamma\) is the transition from juvenile to wild, and \(\mu\) is the juvenile mortality.

The three primary biological processes are described by the following nonlinear differential equations

\[\frac{dM}{dt} = rM\left(1 - \frac{M}{K}\right) + \gamma E - \delta M, \tag{1}\] \[\frac{dP}{dt} = \alpha M - (\beta + \mu)P, \tag{2}\] \[\frac{dE}{dt} = \beta P - (\gamma + \mu)E. \tag{3}\]

\[\frac{dP}{dt} = \alpha M - (\beta + \mu)P,\tag{2}\]

\[\frac{dE}{dt} = \beta P - (\gamma + \mu)E. \tag{3}\]

All parameters and variables in the dynamical system, including their biological meaning and units are presented in Table 1. This system captures the dynamics of Bali Starling conservation by incorporating natural growth, transfer between compartments, and the impact of control interventions on the captive population.

2.2. Optimal Control Problem

The objective of this study is to determine an optimal control c(t) that minimizes the cost of managing the Bali Starling population in captivity, while ensuring the sustainability of the population. The control variable c(t) represents additional feeding and management efforts applied to the captive population P(t), as seen in Equation (6). It is explicitly shown as an arrow in Figure 1, it influences the transition process (transfer efficiency) and is penalized in the cost functional.

Figure 1: Flow diagram of the Bali Starling population dynamics is the representation of Equation (5)-(7). Arrows indicate transitions between compartments: M (wild), P (captive), and E (juvenile/exposure). The natural growth, transfer, maturation, and mortality processes are represented alongside the influence of the control variable c(t)(additional feeding).

Table 1: Description of variables and parameters in the dynamical system.

The cost functional to be minimized is given by

\[J(c) = \int_0^T \left( bP(t) + qc^2(t) \right) dt. \tag{4}\]

where b is the weight associated with the maintenance cost of the captive population, q is the weight representing cost of implementing the control, c(t) is the control function representing additional efforts (in our case is the feeding). The control variable c(t) is assumed to lie within the admissible set \(\mathcal{U} = \{c \in \{c\}\}\)\(L^{\infty}[0,T]:0\leq c(t)\leq 1\), reflecting practical constraints on the intensity of management intervention. Here, c=0 corresponds to a baseline level of feeding effort, while c=1 represents the maximum achievable intervention within the operational limits of the breeding program.

The first term of Equation (4) penalizes the size of the captive population, while the second term penalizes the cost of additional interventions. Now, the optimal control minimizes (J(c)), subject to the state equations

\[\frac{dM}{dt} = rM\left(1 - \frac{M}{K}\right) + \gamma E - \delta M, \tag{5}\] \[\frac{dP}{dt} = \alpha M - (c + \beta + \mu)P, \tag{6}\] \[\frac{dE}{dt} = (c + \beta)P - (\gamma + \mu)E. \tag{7}\]

\[\frac{dP}{dt} = \alpha M - (c + \beta + \mu)P,\tag{6}\]

\[\frac{dE}{dt} = (c+\beta)P - (\gamma + \mu)E. \tag{7}\] and initial conditions

\[M(0) = M_0, \quad P(0) = P_0, \quad E(0) = E_0.\] (8)

The term -c(t)P(t) represents the controlled transition of individuals from the captive population to the juvenile exposure class. Here, \(c(t) \in [0,1]\) models the intensity of additional management effort (e.g., feeding, training, habituation), which facilitates the progression of birds from captivity to exposure. A higher value of c(t) increases the outflow rate from P, accelerating preparation for eventual release into the wild.

To solve this optimal control problem, we apply Pontryagin's Maximum Principle to derive the associated Hamiltonian system and characterize the optimal control policy. The Hamiltonian H is formulated by combining the objective function and the state equations, weighted by the co-state variables \((\lambda_1, \lambda_2, \lambda_3)\) [21]:

\[H = bP + qc^{2} + \lambda_{1} \left( rM \left( 1 - \frac{M}{K} \right) + \gamma E - \delta M \right) + \lambda_{2} \left( \alpha M - (c + \beta + \mu)P \right) + \lambda_{3} \left( (\beta + c)P - (\gamma + \mu)E \right).\] \[(9)\]

The Pontryagin Maximum Principle provides necessary conditions for optimality in control problems. It states that for an optimal control \(c^*\), there exists a set of costate variables \(\lambda_1, \lambda_2, \lambda_3\) such that the Hamiltonian H is minimized with respect to the control variables at each point in time [1]. The principle involves the following steps: First consider the differential equations governing the system dynamics as in (5)-(7). Secondly, we define the differential equations for the costate variables by taking the negative partial derivative of the Hamiltonian with respect to each state variable. This gives

\[\frac{d\lambda_1}{dt} = -\frac{\partial H}{\partial M} = -\lambda_1 \left( r - \frac{2rM}{K} - \delta \right) - \lambda_2 \alpha. \tag{10}\]

\[\frac{d\lambda_2}{dt} = -\frac{\partial H}{\partial P} = -b + \lambda_2(\beta + \mu + c) - \lambda_3(\beta + c). \tag{11}\]

\[\frac{d\lambda_3}{dt} = -\frac{\partial H}{\partial E} = -\lambda_1 \gamma + \lambda_3 (\gamma + \mu). \tag{12}\]

To derive the optimal control law, we take the partial derivative of the Hamiltonian with respect to the control variable c and set it to zero

\[\frac{\partial \mathcal{H}}{\partial c} = 2qc - \lambda_2 P + \lambda_3 P = 0. \tag{13}\]

Solving (13) for c, we get

\[c^*(t) = \frac{(\lambda_2 - \lambda_3)P}{2a}. (14)\] with bounds

\[c^*(t) = \min\left(1, \max\left(0, \frac{(\lambda_2 - \lambda_3)P}{2q}\right)\right). \tag{15}\]

This is the optimal control law for the problem.

3. MODEL ANALYSIS

3.1. Equilibrium Points

Equilibrium points in population dynamics reveals a rich spectrum of behaviors in nonlinear systems. This includes stable equilibrium, cyclic oscillations, and even chaos. Our discussion in section follows the concept of equilibrium that lead to phenomena like extinction or invasion as discussed in [20]. Based on the dynamical system , we identify the equilibrium points by solving

\[\frac{dM}{dt} = 0, \quad \frac{dP}{dt} = 0, \quad \frac{dE}{dt} = 0.\]

The trivial equilibrium point is the state of extinction. This corresponds to the absence of population in all compartments or \((M^*, P^*, E^*) = (0, 0, 0)\) This equilibrium represents complete extinction and often serves as a reference point in stability analysis. The non-trivial equilibrium point assumes all populations are nonzero at equilibrium, we denote \((M^*, P^*, E^*) \neq (0, 0, 0)\). Now Solving the steady-state equations, from Equation 3), we have

\[0 = \beta P^* - (\gamma + \mu)E^* \quad \Rightarrow \quad E^* = \frac{\beta}{\gamma + \mu}P^*.\]

From Equation (2), we get

\[0 = \alpha M^* - (\beta + \mu)P^* \quad \Rightarrow \quad M^* = \frac{(c + \beta + \mu)}{\alpha}P^*.\]

By substituting \(E^*\) and \(M^*\) into Equation (1), we get

\[0 = rM^* \left( 1 - \frac{M^*}{K} \right) + \gamma E^* - \delta M^*.\]

By inserting expressions, we get

\[0 = r \left( \frac{(c+\beta+\mu)}{\alpha} P^* \right) \left( 1 - \frac{(c+\beta+\mu)}{\alpha K} P^* \right) + \gamma \cdot \frac{\beta}{\gamma+\mu} P^* - \delta \cdot \frac{(c+\beta+\mu)}{\alpha} P^*.\]

This is a nonlinear algebraic equation in \(P^*\), which can be solved numerically (or analyzed for stability), and then substituted back to find \(M^*\) and \(E^*\).

3.2. Existence and Stability Conditions

The stability of population dynamics is a crucial area of research in ecology. Several studies have explored various aspects of this topic, including structural stability in nonlinear systems [7]. In this paper, we follow the stability criteria for equilibria of nonlinear matrix population models as discussed in [43].

Theorem 3.1. (Existence of non-trivial equilibrium) Consider the dynamical system (1)-(3) with parameters \(r, K, \delta, \alpha, \beta, \gamma, \mu > 0\). Then, there exists a unique non-trivial equilibrium point \((M^*, P^*, E^*)\) with all components strictly positive, provided

\[\frac{r}{\delta} > \frac{\alpha(\beta + \mu)(\gamma + \mu)}{\gamma\beta}.\] (16)

Proof: To prove this theorem, we start by finding the equilibrium and setting the right-hand sides of the system to zero. Firstly we solve Equation (2) of the algebraic system, giving

\[0 = \alpha M^* - (\beta + \mu)P^* \quad \Rightarrow \quad P^* = \frac{\alpha}{\beta + \mu}M^*. \tag{17}\]

From Equation (3), we get

\[0 = \beta P^* - (\gamma + \mu)E^* \quad \Rightarrow \quad E^* = \frac{\beta}{\gamma + \mu}P^*. \tag{18}\]

Substituting (17) into (18), we get

\[E^* = \frac{\beta}{\gamma + \mu} \cdot \frac{\alpha}{\beta + \mu} M^* = \frac{\alpha \beta}{(\beta + \mu)(\gamma + \mu)} M^*. \tag{19}\]

Next, by substituting (19) into Equation (1), we get

\[0 = rM^* \left( 1 - \frac{M^*}{K} \right) + \gamma E^* - \delta M^*. \tag{20}\]

Substituting (19) into Equation (20), we get

\[0 = rM^* \left( 1 - \frac{M^*}{K} \right) + \gamma \cdot \frac{\alpha \beta}{(\beta + \mu)(\gamma + \mu)} M^* - \delta M^*. \tag{21}\]

By factoring out \(M^* \neq 0\) on Equation (21), we get

\[0 = M^* \left[ r \left( 1 - \frac{M^*}{K} \right) + \frac{\gamma \alpha \beta}{(\beta + \mu)(\gamma + \mu)} - \delta \right]. \tag{22}\]

Dividing both sides by \(M^*\) (assuming \(M^* \neq 0\)), we have

\[r\left(1 - \frac{M^*}{K}\right) = \delta - \frac{\gamma \alpha \beta}{(\beta + \mu)(\gamma + \mu)}.\] (23)

Now, we denote

\[A := \frac{\gamma \alpha \beta}{(\beta + \mu)(\gamma + \mu)}.\]

Therefore, Equation (23) becomes

\[1 - \frac{M^*}{K} = \frac{\delta - A}{r} \quad \Rightarrow \quad M^* = K \left( 1 - \frac{\delta - A}{r} \right) = K \left( \frac{r - \delta + A}{r} \right).\]

Hence,

\[M^* = K \cdot \frac{r - \delta + A}{r}. (24)\]

Now for \(M^* > 0\), Equation (24) implies \(r > \delta - A\). By replacing A, we get

\[r > \delta - \frac{\gamma \alpha \beta}{(\beta + \mu)(\gamma + \mu)}. (25)\]

Rearranging Equation (25), we get

\[\frac{r}{\delta} > \frac{\alpha(\beta + \mu)(\gamma + \mu)}{\gamma\beta},\tag{26}\] then the right-hand side of (24) is positive, and thus \(M^* > 0\), and by substitution into (17) and (18), both \(P^* > 0\), \(E^* > 0\) as well.

The existence of a non-trivial equilibrium guarantees that the system of differential equations admits a steady-state solution where all subpopulations persist over time. In other words, it confirms the model describes a biologically meaningful and sustainable outcome. It supports the validity of results on stability and optimal control, and provides a clear mathematical foundation for long-term conservation planning.

Theorem 3.2 (Local asymptotic stability of the nontrivial equilibrium). Let \((M^*, P^*, E^*)\) be the nontrivial equilibrium point of system (5)-(7). Then, \((M^*, P^*, E^*)\) is locally asymptotically stable if the Jacobian matrix evaluated at \((M^*, P^*, E^*)\) has all eigenvalues with negative real parts. A sufficient condition for this to hold is

• 1) \(0 < M^* < K\), ensuring the logistic term satisfies \(1 \frac{2M^*}{K} < 1\), and
• 2) The following inequality holds:

\[\alpha \gamma \beta > (\beta + \mu)(\gamma + \mu) \left(\delta + r \left(1 - \frac{2M^*}{K}\right)\right).\]

This guarantees that the effective reproductive rate exceeds the system's dampening effects due to mortality and density dependence, thus yielding local asymptotic stability.

Proof: Let the system be written as \(\dot{\mathbf{x}} = \mathbf{f}(\mathbf{x})\), with \(\mathbf{x} = [M, P, E]^T\). The Jacobian matrix is

\[J = \begin{bmatrix} \frac{\partial f_1}{\partial M} & \frac{\partial f_1}{\partial P} & \frac{\partial f_1}{\partial E} \\ \frac{\partial f_2}{\partial M} & \frac{\partial f_2}{\partial P} & \frac{\partial f_2}{\partial E} \\ \frac{\partial f_3}{\partial M} & \frac{\partial f_3}{\partial P} & \frac{\partial f_3}{\partial E} \end{bmatrix}.\]

Computing the partial derivatives, then the Jacobian matrix of the system is given by

\[J(M^*, P^*, E^*) = \begin{bmatrix} r\left(1 - \frac{2M^*}{K}\right) - \delta & 0 & \gamma \\ \alpha & -(c + \beta + \mu) & 0 \\ 0 & (\beta + c) & -(\gamma + \mu) \end{bmatrix}.\] (27)

For simplicity, let us define

\[a_{11} = r \left( 1 - \frac{2M^*}{K} \right) - \delta,\]
\[a_{22} = -(c + \beta + \mu),\]
\[a_{33} = -(\gamma + \mu).\]

The characteristic polynomial is given by

\[\det(J - \xi I) = (\xi - a_{11})(\xi - a_{22})(\xi - a_{33}) + \alpha \beta \gamma = 0.\] (28)

Expanding Equation (28) gives

\[\xi^3 - (a_{11} + a_{22} + a_{33})\xi^2 + (a_{11}a_{22} + a_{11}a_{33} + a_{22}a_{33})\xi - (a_{11}a_{22}a_{33} + \alpha\beta\gamma) = 0.\]

Let the polynomial be written as

\[\xi^3 + A_1 \xi^2 + A_2 \xi + A_3 = 0,\] where A¹ = −(a¹¹ + a²² + a33), A² = a11a²² + a11a³³ + a22a33, A³ = −(a11a22a³³ + αβγ). Routh–Hurwitz stability criteria state that the equilibrium is locally asymptotically stable if

\[A_1 > 0\], \(A_3 > 0\), \(A_1 A_2 > A_3\).

Since a22, a³³ < 0, their product is positive. If M^∗ < K, then a¹¹ is also negative. Hence

\[A_1 > 0\] (sum of negatives with a minus sign), \(A_3 > 0 \iff \alpha \beta \gamma > -a_{11}a_{22}a_{33}\).

This gives the sufficient condition

\[\alpha\beta\gamma > (c + \beta + \mu) (\gamma + \mu) \left(\delta + r \left(1 - \frac{2M^*}{K}\right)\right).\]

Therefore, if M^∗ < K and the above inequality is satisfied, then (M^∗ , P^∗ , E^∗ ) is locally asymptotically stable.

3.3. Existence and Uniqueness of the Optimal Control

The existence and uniqueness of the optimal control for the Bali Starling population management problem refers to the results from optimal control theory. The existence is often guaranteed under certain conditions on the dynamics of the system, the cost function, and the control set. Several papers explore conditions guaranteeing these properties. Pedregal & Tiago [29] present a general approach for proving existence without relying on typical convexity conditions. Clarke [9] discusses existence theorems requiring convex velocity sets and growth restrictions. Rosenblueth [33] investigates links between uniqueness of multipliers and secondorder conditions in optimal control. Song et al. [41] provides an overview of optimal control methods, including Pontryagin's maximum principle and the Hamilton-Jacobi-Bellman equation. By Referring to the papers, we state the existence of an optimal control, as follows.

Theorem 3.3. (Existence of Optimal Control). Consider the optimal control problem governed by the system of differential equations (5)-(7) and the objective function given by Equation (4). Assume the following conditions hold:

• 1) The admissible control set given by U = {c(t) | cmin ≤ c(t) ≤ cmax, for all t ∈ [0, T]}, where cmin and cmax are finite bounds.
• 2) The function f(M, P, E, c) defining the system dynamics is continuous and satisfies the linear growth condition:

\[||f(M, P, E, c)|| \le C(1 + ||M|| + ||P|| + ||E||), \quad \forall (M, P, E, c) \in \mathbb{R}^3 \times \mathcal{U}.\] 3) The function L(P, c) = bP + qc² is convex in c and satisfies the coercivity condition

\[L(P,c) \ge \psi(P)\] where ψ(P) is a non-negative function that grows unbounded as |P| → ∞. Then, there exists an optimal control function c ∗ (t) ∈ U such that:

\[J(c^*) = \min_{c \in \mathcal{U}} J(c).\]

Proof: The proof follows standard arguments in optimal control theory, relying on the Filippov-Cesari Theorem [8] and Weierstrass' Theorem [15]. The control set U is compact by assumption since it is a closed and bounded subset of R. Since the system is defined over the finite time interval [0, T], we consider the set of admissible controls:

\[\mathcal{U} = \{c : [0, T] \to [c_{\min}, c_{\max}] \text{ measurable}\}.\]

By Filippov–Cesari's Theorem [8], since f(M, P, E, c) is continuous and satisfies a linear growth condition, the system has a well-posed solution for any admissible control. Now, we define the set of all state trajectories:

\[\mathcal{X} = \{(M(t), P(t), E(t)) \text{ satisfying the system dynamics for some } c(t) \in \mathcal{U}\}.\]

Since the control set is bounded, the right-hand side of the system satisfies the linear growth condition:

\[||f(M, P, E, c)|| \le C(1 + ||M|| + ||P|| + ||E||),\] which implies that the set of solutions is precompact in C([0, T], R 3 ) (the space of continuous functions). By the Arzela–Ascoli Theorem [ ` 2], this ensures compactness. Now, we show that the cost function is lower semicontinuous. Consider the objective function (4), since L(P, c) = bP +qc² is convex in c and continuous in both P and c, so it is lower semicontinuous. By Fatou's Lemma, the integral functional J(c) is also lower semicontinuous in c. Thus, by Weierstrass' Theorem [15], the minimum of J(c) exists in the compact control set U. So, the control set U is compact, the system is well-posed with continuous and bounded dynamics, and the cost function is convex and lower semicontinuous, then there exists an optimal control c ∗ (t) that minimizes the cost functional

\[J(c^*) = \min_{c \in \mathcal{U}} J(c).\]

Thus, the theorem is proved.

Based on Theorem 3.3, we claim the following results: (a) the cost function is convex with respect to the control variable, since the term qc² is quadratic in the control variable c, which is convex. Therefore, the cost function is convex with respect to c satisfying this condition. (b) The dynamics of the system (5)-(7) are Lipschitz continuous in both the state and control variables [11]. Our dynamics are linear with respect to the control variable c and nonlinear with respect to the state variables M, P and E. However, these equations are Lipschitz continuous in both the state and control variables, as the nonlinear terms are smooth functions of the state variables. (c) The set of admissible controls is compact. Compactness of the admissible control set is often a key condition for ensuring existence [26]. In our case, the control c belongs to a finite-dimensional control space, and if we restrict to be bounded, c ∈ [cmin, cmax], so the existence is satisfied.

Theorem 3.4. (The uniqueness of optimal control solution). Consider the optimal control problem governed by the system of differential equations, (5)-(7) and the objective function given by Equation (4). Using result in Theorem 3.3, the control c is selected to minimize the objective function (4).

Proof: We will prove uniqueness by showing that the optimal control satisfies the necessary conditions of Pontryagin's Maximum Principle and that these conditions yield a unique solution. Our proof will go through 4 steps. Firstly, by Pontryagin's Maximum Principle, we have the Hamiltonian as given by 9. The necessary condition for optimality requires

\[\frac{\partial \mathcal{H}}{\partial c} = 2qc - \lambda_2 P + \lambda_3 P = 0. \tag{29}\]

Solving (29) for c, we obtain

\[c^*(t) = \frac{(\lambda_2 - \lambda_3)P}{2q}. (30)\] with bounds

\[c^*(t) = \min\left(1, \max\left(0, \frac{(\lambda_2 - \lambda_3)P}{2q}\right)\right).\]

Since q>0, this provides an explicit expression for the control function. Secondly, the objective function (4) contains the term \(qc^2\), which is strictly convex in c because q>0. This ensures that any minimizer of \(J_H\) is unique. Thirdly, since the costate equations (10)-(12) and the system dynamics (5)-(7) satisfy the Lipschitz condition, the Picard-Lindelöf theorem (also known as the Cauchy-Lipschitz theorem), see [13] or [32], guarantees that the system of costate equations has a unique solution for \(\lambda_1, \lambda_2, \lambda_3\). Lastly, since the Hamiltonian is strictly convex in c, its first derivative, (13), is strictly increasing in c, this implies that the solution \(c^*(t) = \frac{(\lambda_2 - \lambda_3)P}{2q}\) is unique. Thus, the theorem is proved.

In short, the uniqueness holds due to the facts that (a) the convexity of the cost function, which ensures a unique minimizer. (b) The Lipschitz continuity of the system, ensuring uniqueness of the state and costate solutions. (c) The strict convexity of the Hamiltonian, which guarantees a unique extremal control function. (d) the objective function is strictly convex with respect to the control variable c, and the system dynamics are linear with respect to c, the Pontryagin Maximum Principle ensures that there is a unique solution to the optimal control problem [34].

4. SIMULATION

4.1. Numerical Result

To solve the system of differential equations representing the population dynamics of the Bali Starling and the optimal control strategy, we employ numerical methods due to the complexity of the equations. The specific method used in this study is the Runge-Kutta method, implemented via the 'solve_ivp' function from the SciPy library in Python. We primarily use the RK45, a Runge-Kutta method of order 5(4) method due to its balance between accuracy and computational efficiency for non-stiff systems. In this paper, we work with real (unnormalized) values, keeping K=800, and use unnormalized initial conditions M(0)=520, P(0)=92, E(0)=40. This is preferred when interpreting real population sizes and making management decisions in practice.

Moreover, in this research we used the following parameter values: r=0.25 (Growth rate of the Bali Starling population), \(\beta=0.50\) (Percentage of chicks in the exposed population), \(\mu=0.20\) (Natural death rate), \(\alpha=0.17\) (Rate of transfer of starlings from wild to captivity), \(\gamma=0.40\) (Exposed population movement coefficient), \(\delta=0.1\) (Mortality rate in the wild), K=800 (Carrying capacity), b=0.30 (Weight of feed in captivity), q=0.25 (Weight of additional control feed for chicks). The initial parameter values used in the model were derived from secondary data provided by the West Bali National Park (WBNP) office in the form of internal population monitoring spreadsheets (WBNP, 2023, unpublished). These include observed transition rates, mortality data, and carrying capacity limits specific to the Bali Starling breeding program at US Tegal Bunder.

The next step is to work out on the numerical experiments. This is aimed to compare the population dynamics of the Bali Starling under two scenarios: with control, \(c^*(t) = (\lambda_2 - \lambda_3)P/2q\) and without control, \((c^*(t) = 0)\). The degree of control depends on q, lower q values lead to stronger intervention, resulting in more regulated population growth, while higher q values allow for less control, causing growth to approach the uncontrolled case. Our goal is to observe the impact of control measures on population growth and to evaluate the cost-effectiveness of the control strategy.

Then, we use the Phyton Library, 'solve_ivp' function to integrate the system of ODEs over the specified time span (0 to 30 months) and extract and plot the population trajectories for total population (M), population in captivity (P), and juvenile population (E) for both scenarios. Calculating the objective function \(J_H = \int_0^t (bP + qc^2) \, dt\) gives the controlled scenario to evaluate the cost-effectiveness. Figure 2 highlights the positive impact of applying an optimal control strategy to the Bali Starling population.

Figure 2 highlights the positive impact of applying an optimal control strategy to the Bali Starling population. With control in place, the wild population grows to around 795 birds, compared to just 735 without control. This shows that targeted interventions, such as carefully timed feeding, structured transitions, and improved juvenile management, help more birds successfully return to the wild. Notably, the juvenile population also increases from approximately 150 to nearly 190 under control, suggesting that the transition from captivity to the wild becomes more effective when guided by a well-planned strategy

Interestingly, the captive population is much lower with control, that is about 80 birds versus 175 without control. This is a good sign because one of the goals of the management plan is to reduce the number of birds kept in captivity, which can be costly and resource-intensive. Instead of holding too many birds in cages, the system focuses on preparing and releasing them efficiently. These results highlight how strategic

Comparison of M(t), P(t), and E(t) With and Without Control

Figure 2: Time evolution of the Bali Starling population under control using unnormalized initial conditions. The system starts from M(0)=520 birds, P(0)=92 birds, and E(0)=40 birds, with the carrying capacity set to K=800 birds. When the optimal control (c(t)) is applied, the trajectories for \(M^*\approx795\), \(P^*\approx80\), and \(E^*\approx190\). When the system is uncontrolled \(c^*(t)=0\), the trajectory for \(M^*\approx735\), \(P^*\approx175\), and \(E^*\approx150\).

efforts not only help the wild population grow but also make the overall breeding and release program more balanced and sustainable in the long run.

Figure 3: The control strategy changes over time for managing the Bali Starling population in a captive breeding system. The vertical axis represents the intensity of the control action, noted as c(t), and the horizontal axis shows time in months. The value of c(t) ranges from 0 to 1, where 0 means no intervention or control is applied, 1 means full intervention or maximum control effort is being used.

Figure 3 shows that at the beginning, the control value is flat at zero. This phase may represent a period of observation or natural adjustment, where the system is allowed to evolve on its own without interference. After that, there is a sudden spike in the control profile, where the intensity jumps from 0 to 1, and stays at maximum intensity for the rest of the simulation. This switch suggests that the optimal strategy recommends full intervention after a certain point.

4.2. Discussion

The simulation results show two very different paths for managing the Bali Starling population: one using a carefully planned control strategy and the other with no intervention at all. In the scenario without control, the wild population grows quickly, which might seem like a good thing at first. However, this rapid growth comes with a big increase in the number of birds in captivity, which could be difficult and costly to manage. On the other hand, the scenario with control (Figure 2) follows a more balanced path, where management efforts (like feeding or conditioning) are applied strategically over time to guide the population more efficiently.

What's especially interesting is that the controlled system actually ends up with a larger wild population (about 795 birds) than the uncontrolled one (around 735 birds). At the same time, it keeps the number of birds in captivity much lower, which is a win from both a practical and financial point of view. The juvenile group, which is essential for future reintroduction to the wild, is also larger when control is applied. These results show that the goal of control is not just to increase numbers but to create a stable, sustainable balance across all parts of the population. The control effort starts off strong, then gradually decreases as the system stabilizes, reducing long-term costs and effort.

By contrast, when no control is applied, the captive population grows quite large. While that might sound like success, overcrowding can cause real problems, like disease, stress, and lower breeding rates, which have been documented in many conservation programs [22]. Plus, managing a large group of birds in captivity requires more space, food, and staff—resources that may not always be available.

The controlled model, Figure 3 is a kind of bang-bang control profile (no control, then full control). This is common in optimal control theory when trying to balance ecological outcomes with limited resources [44]. It helps ensure the wild population thrives, keeps the captive group manageable, and prepares juveniles for successful release. This mirrors real-life conservation efforts, where decisions are often shaped by limited resources and the need for long-term success. So even though the uncontrolled model may produce bigger

numbers at first (captive population, P(t)), the controlled strategy is more thoughtful and realistic, it helps protect the species in a way that's effective, affordable, and sustainable over time.

5. CONCLUSION

This study presents a dynamic population model and an optimal control strategy designed to support the sustainable management of the Bali Starling population in a captive breeding setting. The model captures key biological processes, such as natural growth, transfer to captivity, and juvenile habituation, through a system of nonlinear differential equations. Using Pontryagin's Maximum Principle, we derive an optimal feeding strategy that not only promotes population health but also minimizes the operational costs of captive management.

The research addresses two main objectives: (1) to develop a model that accurately represents the Bali Starling breeding system at the USS Tegal Bunder facility in West Bali National Park, and (2) to determine the optimal carrying capacity and feeding intervention required to achieve both ecological and financial sustainability. Analytical results confirm the existence and stability of a meaningful non-trivial equilibrium, while numerical simulations demonstrate that the controlled approach leads to better outcomes than an unmanaged one. Specifically, the optimal strategy results in a higher wild population, a larger releaseready juvenile population, and a smaller captive population, achieving a more balanced and cost-effective distribution. These findings provide valuable insights for developing more targeted, data-informed conservation policies and can serve as a framework for managing other endangered species facing similar habitat and resource constraints.

ACKNOWLEDGEMENT

We would like to thank Udayana University, in particular the Institute of Research and Community Service, and the Faculty of Mathematics and Natural Science which has funded this research under the Letter of Agreement No.: B/1.116/UN14.4.A/PT.01.03/2023, 02 May 2023.