Subsections

• Formulation
• Performance

6 Marine Population Dynamics

Given estimates of the abundance of the population of a marine species at each stage (for example, nauplius, juvenile, adult) as a function of time, determine stage specific growth and mortality rates. The model for the population dynamics of the $n_s$-stage population is

$$y_{j}' = g_{j-1} y_{j-1} - (m_{j} + g_{j}) y_{j} , \qquad 1 \le j \le n_s ,$$ (3)

where $m_i$ and $g_i$ are the unknown mortality and growth rates at stage $i$ with $g_0 = g_{n_s} = 0$. This model assumes that the species eventually dies or grows into the next stage, with the implicit assumption that the species cannot skip a stage. Initial conditions for the differential equations are unknown, since the stage abundance measurements at the initial time might also be contaminated with experimental error. We minimize the error between computed and observed data,

$$\sum_{j=1}^{n_m} \Vert y(\tau_j;m,g) - z_j\Vert^2 ,$$

where $z_j$ are the stage abundance measurements. This problem is based on the work of Rothschild, Sharov, Kearsley, and Bondarenko [19].

Formulation

Our formulation of the marine population dynamics uses a $k$-stage collocation method, a uniform partition with $n_h$ subintervals of $[0,\tau_{n_m}]$, and the standard [2, pages 247-249] basis representation,

$$\ v_{i} + \sum_{j=1}^{k} \frac{(t - t_i)^j}{j!\ h^{j-1}} w_{ij} , \qquad t \in [t_i,t_{i+1}] ,$$

for the components of the solution $y$ of . The constraints in the optimization problem are the continuity conditions and the collocation equations. The continuity equations are a set of $n_s(n_h -1 )$ linear equations. The collocation equations are a set of $k \, n_s \, n_h$ nonlinear equations obtained by requiring that the collocation approximation satisfy at the collocation points $\xi_{ij} = t_i + h \rho_j$ for $i=1,\ldots,n_h$ and $j=1,\ldots, k$.

Table 6.1: Marine population dynamics problem data

Variables	$(k+1) n_s n_h + 2n_s-1$
Constraints	$(k+1) n_s n_h - n_s$
Bounds	$2n_s-1$
Linear equality constraints	$n_s (n_h-1)$
Linear inequality constraints	0
Nonlinear equality constraints	$k n_s n_h$
Nonlinear inequality constraints	0
Nonzeros in $\nabla ^2 f(x)$	$(k+1)^2 n_s n_m$
Nonzeros in $c'(x)$	$(2k+1)(k+2) n_s n_h$

The parameters in the problem are the $n_s n_h$ initial conditions, the $n_s$ mortality rates, the $n_s -1$ growth rates, and the $(k+1) n_s n_h$ basis parameters in the representation of the collocation approximation. Data for this problem appears in Table 6.1.

We do not impose any initial conditions on the differential equations, since initial measurements are usually contaminated with experimental error. Introducing these extra degrees of freedom into the problem formulation should allow solvers to find a better fit to the data. A significant difference between this problem and other parameter estimation problems is that the population dynamics data usually contains large observation errors.

Performance

We provide results for the AMPL formulation with $k = 2$ in Table 6.2. We use a simulated dataset with $n_s = 8$ stages. The initial basis parameters are chosen so that the collocation approximation is piecewise constant and interpolates the data.

Table 6.2: Performance on marine population dynamics problem

Solver	$n_h=25$	$n_h=50$	$n_h=100$	$n_h=200$
LANCELOT	623.75 s	1084.33 s	3170.26 s	$\ddagger$
$f$	1.97522e+07$\dagger$	1.97465e+07$\dagger$	1.97465e+07$\dagger$	$\ddagger$
$c$ violation	1.80930e-06$\dagger$	3.27200e-06$\dagger$	6.49480e-06$\dagger$	$\ddagger$
iterations	245	281	243	$\ddagger$
LOQO	2.07 s	4.64 s	12.53 s	38.4 s
$f$	1.97522e+07	1.97465e+07	1.97465e+07	1.97465e+07
$c$ violation	5.3e-10	7.0e-10	5.8e-11	2.8e-10
iterations	25	25	27	27
MINOS	6.58 s	15.75 s	167.83 s	245.1 s
$f$	1.97522e+07	1.97465e+07	2.17862e+07	0.00000e+00$\dagger$
$c$ violation	4.5e-12	5.3e-11	4.2e-08	3.4e+05$\dagger$
iterations	12	11	72	31
SNOPT	85.37 s	184.4 s	477.59 s	1502.26 s
$f$	1.97522e+07	1.97465e+07	1.97465e+07	1.97465e+07
$c$ violation	4.5e-12	1.1e-11	7.3e-12	2.0e-11
iterations	411	379	416	519
$\dagger$ Errors or warnings. $\ddagger$ Timed out.

LANCELOT returns the message step got too small for the values of $n_h$ for which it terminates within 3,600 wall-clock seconds. The intermediate solution returned by LANCELOT upon termination is in close agreement with the optimal solutions returned by the other solvers. MINOS makes no progress with $n_h = 200$, returning with the error unbounded (or badly scaled) problem.

The graph on the left of Figure 6.1 shows the populations for stages 1, 2, 5, and 6, while the graph on the right shows the populations for stages 3, 4, 7, and 8. In both cases, the fit between the model and the data is not always tight.

For this problem we are using a relatively small number of collocation points ($k=2$), since in this case the number of parameters grows quickly with the number of stages. The quality of the solution does not seem to be affected, at least as measured by the population curves and the mortality and growth parameters.

Figure 6.1: Marine populations for stages 1, 2, 5, 6 (left) and 3, 4, 7, 8 (right)

Figure 6.2 plots the mortality and growth parameters for the eight stages. Mortality parameters are marked $*$, while growth parameters are marked $\circ$. The mortality parameters for stages 5 and 6 are not zero, but they are on the order of $10^{-3}$ and $10^{-9}$, respectively.

Figure 6.2: Mortality ($*$) and growth ($\circ$) parameters for the marine populations stages