Walters and Martell (2004, Section 3.2) present an interesting example of harvest feedback policy exploration using basic spreadsheet techniques and optimization procedures (e.g. Solver in Excel, or Optimizer in Quattro Pro). For a given species, parameters describing growth, mortality, and reproduction/recruitment are used to model the population's dynamics. Given a random series of environmental effects to recruitment parameters (survival and carrying capacity) the optimization procedure can be used to determine the best series of fishing mortality (F) values given a cost function (e.g. maximize the sum of: harvest yield, log-transformed yield, discounted yield):
In this example, for a hypothetical population of tilapia, yearly recruitment survival varies by ca. +/-30%, while recruitment carrying capacity has two periods of poor recruitment (-50%, e.g. due to a regime shift) during the simulation years 20-35 and 60-75, as indicated by the shaded areas. The optimal F series are similar for maximizing yield or discounted yield (discount rate = 5%) likely due to the fact that tilapia is a fast-growing, high fecundity species and can rebound quickly following intensive fishing pressure. Maximizing log-transformed biomass is the "risk-averse" strategy, whereby peaks in harvest are less highly valued, and fluctuations in fishing pressure (and yield) are dampened:
Interesting patterns emerge, e.g. one should increase fishing pressure at the beginning of a "bad" recruitment regime (grey regions) and reduce fishing at the beginning a "good" recruitment regime. In both cases, this speeds the trajectory of the stock to it's new equilibrium and takes advantage of soon to be lost or gained harvestable biomass. In other words, fish hard before bad times so as not to "waste" a large spawning stock, and lay off of fishing before good times so as to make full use of the spawning potential. Martell and Walters (2004) state,
These anticipatory behaviors are a bizarre reversal of our usual prescriptions about how to hedge against uncertain future changes in productivity and would almost certainly never be acceptable as management recommendations in practical decision settings.
Finally, the authors also demonstrate the situation where all sizes are vulnerable to fishing (as opposed to a knife-edge selection of individuals above a given minimum size, as in the above example). Here the sum(Yt) optimization results in a "pulse-fishing" strategy whereby the stock is fished hard (usually E > 0.5) for 1-2 years, followed by 1-2 years of recovery where no fishing is allowed. This pattern results from "growth overfishing", whereby unselective harvesting wastes a large part of the biomass that is still growing rapidly:
The examples use the stockSim() and optim.stockSim() functions of the R package fishdynr. Instructions for direct installation from GitHub can be found here: https://github.com/marchtaylor/fishdynr. Vector optimization of the F series is computed with the optim() function (stats package).
References:
- Walters, C. J., Martell, S. J., 2004. Fisheries ecology and management. Princeton University Press.
- Walters, C., Parma, A. M., 1996. Fixed exploitation rate strategies for coping with effects of climate change. Canadian Journal of Fisheries and Aquatic Sciences, 53(1), 148-158.
Script to reproduce the example:
library(fishdynr) set.seed(1) data(tilapia) params <- tilapia params$knife_edge_size <- 20 params$N0 <- 1e9 nyears <- 100 Ft <- rep(0.5, nyears) env_at <- runif(nyears, min=0.7, max=1.3) env_bt <- rep(1, nyears); env_bt[c(20:35, 60:75)] <- 0.5 # Optimization Ft for sum(yield) out1 <- optim( par = Ft, # initial guess fn = optim.stockSim, params = params, nyears = nyears, env_at = env_at, env_bt = env_bt, opt = "sum", method = "L-BFGS-B", lower = 0, upper = 2, control = list(fnscale=-1, trace=4) ) res1 <- stockSim(params, nyears=nyears, Ft=out1$par, env_at=env_at, env_bt=env_bt) res1$Et <- res1$Yt/res1$Bt plot(Ft ~ t, res1, t="l") # Optimization Ft for sum(log(yield)) out2 <- optim( par = Ft, # initial guess fn = optim.stockSim, params = params, nyears = nyears, env_at = env_at, env_bt = env_bt, opt = "sum.log", method = "L-BFGS-B", lower = 0, upper = 2, control = list(fnscale=-1, trace=4) ) res2 <- stockSim(params, nyears=nyears, Ft=out2$par, env_at=env_at, env_bt=env_bt) res2$Et <- res2$Yt/res2$Bt # Optimization Ft for sum(discounted yield) out3 <- optim( par = Ft, # initial guess fn = optim.stockSim, params = params, nyears = nyears, env_at = env_at, env_bt = env_bt, opt = "sum.disc", disc.rate = 0.05, method = "L-BFGS-B", lower = 0, upper = 2, control = list(fnscale=-1, trace=4) ) res3 <- stockSim(params, nyears=nyears, Ft=out3$par, env_at=env_at, env_bt=env_bt) res3$Et <- res3$Yt/res3$Bt # plots and results ==================================== # optimum fishing mortality series png("optim_Ft.png", width=5, height=4, units="in", res=300, family="Arial", type="cairo") op <- par(mar=c(3,3,1,1), mgp=c(2,0.5,0), tcl=-0.25) plot(Ft ~ t, res1, xlab="t", ylab="Fishing mortality", t="n", log="") rect(20, -100, 35, 1000, col="grey90", border=NA) rect(60, -100, 75, 1000, col="grey90", border=NA) lines(res1$Ft, col=1) lines(res2$Ft, col=2) lines(res3$Ft, col=3) box() legend("topleft", legend=c("Y", "log(Y)", "disc. Y"), title="Optimize:", ncol=3, col=1:3, lt=1) par(op) dev.off() # optimum exploitation rate series png("optim_Et.png", width=5, height=4, units="in", res=300, family="Arial", type="cairo") op <- par(mar=c(3,3,1,1), mgp=c(2,0.5,0), tcl=-0.25) plot(Et ~ t, res1, ylim=c(0,1), xlab="t", ylab="Exploitation rate", t="n", log="") rect(20, -100, 35, 1000, col="grey90", border=NA) rect(60, -100, 75, 1000, col="grey90", border=NA) lines(res1$Et, col=1) lines(res2$Et, col=2) lines(res3$Et, col=3) box() legend("topleft", legend=c("Y", "log(Y)", "disc. Y"), title="Optimize:", ncol=3, col=1:3, lt=1) par(op) dev.off() # optimum yield series png("optim_Yt.png", width=5, height=4, units="in", res=300, family="Arial", type="cairo") op <- par(mar=c(3,3,1,1), mgp=c(2,0.5,0), tcl=-0.25) plot(Yt ~ t, res1, t="n", ylab="Yield") rect(20, -1e20, 35, 1e20, col="grey90", border=NA) rect(60, -1e20, 75, 1e20, col="grey90", border=NA) lines(Yt ~ t, res1, col=1) lines(Yt ~ t, res2, col=2) lines(Yt ~ t, res3, col=3) sum(res2$Yt)/sum(res1$Yt) sum(res3$Yt)/sum(res1$Yt) box() legend("topright", legend=c("Y", "log(Y)", "disc. Y"), title="Optimize:", ncol=3, col=1:3, lt=1) par(op) dev.off() # optimum yield versus stock biomass fit1 <- lm(Yt ~ Bt, res1) fit2 <- lm(Yt ~ Bt, res2) fit3 <- lm(Yt ~ Bt, res3) png("Yt~Bt.png", width=5, height=5, units="in", res=300, family="Arial", type="cairo") op <- par(mar=c(3,3,1,1), mgp=c(2,0.5,0), tcl=-0.25) plot(Yt ~ Bt, res1, t="o", lty=3, pch=1, xlab="Biomass", ylab="Yield") lines(Yt ~ Bt, res3, t="o", lty=3, pch=1, col=3) lines(Yt ~ Bt, res2, t="o", lty=3, pch=1, col=2) abline(fit1, lwd=2) abline(fit2, col=2, lwd=2) abline(fit3, col=3, lwd=2) usr <- par()$usr text(x=usr[1]+1*(usr[2]-usr[1]),y=usr[3]+0.15*(usr[4]-usr[3]), labels=bquote(Yt == .(sprintf("%.2e", coef(fit1)[1]))+.(sprintf("%.2f", coef(fit1)[2]))*Bt), pos=2, col=1 ) text(x=usr[1]+1*(usr[2]-usr[1]),y=usr[3]+0.1*(usr[4]-usr[3]), labels=bquote(Yt == .(sprintf("%.2e", coef(fit2)[1]))+.(sprintf("%.2f", coef(fit2)[2]))*Bt), pos=2, col=2 ) text(x=usr[1]+1*(usr[2]-usr[1]),y=usr[3]+0.05*(usr[4]-usr[3]), labels=bquote(Yt == .(sprintf("%.2e", coef(fit3)[1]))+.(sprintf("%.2f", coef(fit3)[2]))*Bt), pos=2, col=3 ) legend("topleft", legend=c("Y", "log(Y)", "disc. Y"), title="Optimize:", ncol=1, col=1:3, lt=1) par(op) dev.off() # Minimum stock size before allowing harvest (x intercept) sprintf("%.2e", -coef(fit1)[1]/coef(fit1)[2]) sprintf("%.2e", -coef(fit2)[1]/coef(fit2)[2]) sprintf("%.2e", -coef(fit3)[1]/coef(fit3)[2]) # Compare constant F to optimal solution c1 <- stockSim(params, nyears=nyears, Ft=-log(1-coef(fit1)[2]), env_at=env_at, env_bt=env_bt) c2 <- stockSim(params, nyears=nyears, Ft=-log(1-coef(fit2)[2]), env_at=env_at, env_bt=env_bt) c3 <- stockSim(params, nyears=nyears, Ft=-log(1-coef(fit3)[2]), env_at=env_at, env_bt=env_bt) sum(c1$Yt)/sum(res1$Yt) sum(c2$Yt)/sum(res2$Yt) sum(c3$Yt)/sum(res3$Yt) # Pulse fishing (all sizes vulnerable to fishing) ===================== params$knife_edge_size <- 0 # Optimization Ft for sum(yield) out4 <- optim( par = Ft, # initial guess fn = optim.stockSim, params = params, nyears = nyears, env_at = env_at, env_bt = env_bt, opt = "sum", method = "L-BFGS-B", lower = 0, upper = 2, control = list(fnscale=-1, trace=4) ) res4 <- stockSim(params, nyears=nyears, Ft=out4$par, env_at=env_at, env_bt=env_bt) res4$Et <- res4$Yt/res4$Bt png("Et_Bt_pulse.png", width=5, height=5, units="in", res=300, family="Arial", type="cairo") layout(matrix(1:2, nrow=2, ncol=1), widths=5, heights=c(2,3), respect=TRUE) layout.show(2) op <- par(mgp=c(2,0.5,0), tcl=-0.25) par(mar=c(0.5,3,1,1)) plot(Et ~ t, res4, ylim=c(0,1), t="n", xaxt="n", ylab="Exploitation rate") rect(20, -1e20, 35, 1e20, col="grey90", border=NA) rect(60, -1e20, 75, 1e20, col="grey90", border=NA) lines(Et ~ t, res4, col=1) axis(1, labels=FALSE) box() par(mar=c(3,3,0.5,1)) plot(Bt ~ t, res4, t="n", ylab="Biomass", ylim=c(1,0.7)*range(res4$Bt)) rect(20, -1e20, 35, 1e20, col="grey90", border=NA) rect(60, -1e20, 75, 1e20, col="grey90", border=NA) lines(Bt ~ t, res4, col=1) box() dev.off()
No comments:
Post a Comment