Linear Optimization for Ecological Indices Based on Aggregation Functions

Abstract

We consider an optimization problem in ecology where our objective is to maximize biodiversity with respect to different land-use allocations. As it turns out, the main problem can be framed as learning the weights of a weighted arithmetic mean where the objective is the geometric mean of its outputs. We propose methods for approximating solutions to this and similar problems, which are non-linear by nature, using linear and bilevel techniques.

Appendix: Implementations

We have implemented all three approaches to optimization as functions in an R library available at our website⁴.

Optimizing total abundance: eco.opti().

Description of inputs

• species.data - matrix of species abundances per unit of land area, i.e. with $x_{ij}$ denoting the i-th species and its abundance for land type j;

• densities - vector of densities per unit of land area, i.e. with $d_j$ denoting the human population density for a given land-type j;

• tot.pop - the total population required to be fit into the given land area;

• tot.land - the total area over which we need to distribute the population.

Additional optional constraints

• w.min / w.max - vectors denoting minimum or maximum bounds on the land-types, e.g. if we want to ensure that at least 20 % of the land is populated at minimum density we incorporate the constraint $w_1 \ge 0.2$ (assuming $w_1$ is the land-type with minimum density), or alternatively we may wish to limit high density housing to at most 40 % of the land area etc.;

• spec.min / spec.max - vector placing minimum or maximum bounds for a particular species, for example, if we want to make sure that a rare species is above a given threshold $\gamma _i$, the linear constraint $\sum \limits _{j=1}^n w_j x_{ij} \ge \gamma _i$ is added for that species.

The function also gives as output a number of ecological indices such as the individual species abundances, and the Simpson and Shannon diversity indices.

Maximizing the geometric mean of species abundances: eco.opti.gm()

Description of inputs

In addition to all inputs and constraints used with eco.opti(), this function has two additional optional parameters to control the precision.

• fprec - a positive integer giving the number of tangent functions to be defined. The default setting is 100 linear segments, and so gains in accuracy can be achieved with settings of 500, 1000 etc., however obviously at the cost of computation time;

• max.x - a real number giving the maximum value for the domain over which the tangent functions are calculated, the default setting is 10000, and so depending on the scale given it could be necessary to increase this value (or decrease it for finer accuracy) or the optimization will be the same as it would be for maximizing abundance.

For the number of tangent lines K, optimizing over 5 species with 100 linear segments will require $5\times 100 = 500$ additional constraints, use of 1000 linear segments will require 5000 additional constraints and so on. We need to be careful when reducing the precision, since the log function’s gradient changes more drastically for values closer to zero than it does for large values. Rather than taking equal step sizes in calculating our tangent lines, they were distributed using $t_k = \exp (-k\cdot \max (s)/K)$ where $\max (s)$ is the max.x parameter above.

We note also that by default the setting for spec.min will be 1 for all species. It could be adjusted to a fractional value if desired.

Maximizing Shannon’s diversity index: eco.opti.sh()

This function uses the same inputs as with the previous two. The program first solves a maximum and minimum problem using eco.opti() in order to find the feasible bounds to search for M. Note that $-t \ln t$ is concave for $t \in (0,1]$.

Another biodiversity index used as an objective and included in the code made available online is Simpson’s diversity index $1/(\sum _{i=1}^m p_i^2)$. This is performed in a similar manner, however now we are minimizing for a convex function $y = t^2$ rather than maximizing for a concave function and so we need to make the appropriate changes when using the linear framework above.