Sdp Reed Controls Comparisons

library(knitr)
opts_chunk$set(dev='png', fig.width=3, fig.height=3, results='hide', fig.show='hold')

library(pdgControl)
library(ggplot2)

Reed solution; harvest control

f <- function(x, h, p){ 
  x <- sapply(x, function(x) max(x-h,0))
  p[[1]] * x / {1 + p[[2]] * x} 
}
p <- c(a = 6, b = 0.05)
profit <- function(x, h) sapply(x, min, h) # vectorized for x
x_grid <- seq(0, 150, length = 100)
h_grid <- seq(0,150,length = 100)
delta <- 0.05
OptTime <- 1000
xT <- 0
sigma_g <- 0.1
z_g <- function() 1 + rlnorm(1, 0, sigma_g)
pdfn <- function(P, s) dlnorm(P, 0, s)

Note that in \(f\) we must define the harvesting take place before the population recruits.

SDP_Mat <- determine_SDP_matrix(f, p, x_grid, h_grid, sigma_g, pdfn)

opt <- value_iteration(SDP_Mat, x_grid, h_grid, OptTime, xT, profit, delta)

qplot(x_grid, h_grid[opt$D]) + xlab("Stock size") + ylab("Harvest")
qplot(x_grid, x_grid - h_grid[opt$D]) + xlab("Stock size") + ylab("Escapement")

Reed solution with effort-based control

f <- function(x, h, p){
  x <- (1 - h) * x
  p[[1]] * x / {1 + p[[2]] * x}
  }
p <- c(a = 6, b = 0.05)
h_grid <- seq(0, 1, length = 100)
profit <- function(x, h) sapply(x, function(x) x * h) # vectorized for x

SDP_Mat <- determine_SDP_matrix(f, p, x_grid, h_grid, sigma_g, pdfn)

opt <- value_iteration(SDP_Mat, x_grid, h_grid, OptTime, xT, profit, delta)

qplot(x_grid, h_grid[opt$D]) + xlab("Stock size") + ylab("Fishing Effort")
qplot(x_grid, x_grid - x_grid * h_grid[opt$D]) + xlab("Stock size") + ylab("Escapement")