This week I will be continuing developing a simple bioeconomic model, as a follow-up to previous blog posts: [1] and [2]. Last time, I introduced more economics into the model, namely the profit of the fishery. How much effort is exerted into fishing is still exogenous to the model, but we can see how the stock is affected by having different effort strategies.

A natural step is to say that effort should be endogenous. Instead of defining a fishing strategy (i.e. effort level) and “see” what happens, we can start from the premise that whoever decides the effort level should aim at maximizing their benefits, whether these are profits, the social value of the fishery or net benefits. Different agents will be interest in maximizing different values: the Fisher wants to maximize their own profit, while the social planner usually aims at maximizing the social value. They then set the effort level to whatever maximizes their objective. I start by looking at profits.

I am using the same letters as last blog posts. As a summary, here are the parameters used and what they represent as a comment:

```
# Parameters
numberOfYears <- 5 # Number of time periods
X <- c(rep(0, times = numberOfYears+1)) # Stock size vector
X[1] <- 521832 # Stock size in the first time period
r <- 0.302 # stock growth parameter
q <- 0.1 # Catchability parameter
p <- 1293.362 # Price per unit of harvest
c <- 1.5892 # Cost per unit of effort
i <- 0.04 # Discount rate
m <- 521832 # Carrying capacity of the ecosystem (i.e. largest biomass of the species it can sustainably maintain)
```

Given these parameters, I can write the profit of a fishery as the amount harvested time the price per unit harvested, minus the cost per effort times the units of effort exerted. One can write profit at a generic time *t* as:

`Profit <- p*Y[t]-c*E[t]`

Effort is now endogenous, meaning my goal is to find the optimal effort that maximizes the profit expression. Harvest Y is also a function of effort (as well as stock X). For a generic period t, harvest can be written as:

`Y[t] <- q*E[t]*X[t]`

Finally, the stock level at time t is a function of multiple parameters: not only the stock at the previous time period, but also effort at the previous time period, carrying capacity, the catchability parameter, etc.

`X[t] <- X[t-1]+r*X[t-1]*(1-X[t-1]/m)-q*E[t-1]*X[t-1]`

It is time to join all of these conditions together in a single optimization problem. There are several ways to optimize functions in R. Both the functions “optimize” and “optim” are available for optimizaton exercises. Unfortunately neither of these can handle constraints. Instead, I used the function *nloptr*. Bhadani (2021) explains very clearly how to set up the model using *nloptr*. I followed these steps to a large extent.

`library(nloptr)`

*nloptr* needs several different arguments which we set up in an intuitive but specific way. One is the core function, where I write the profit function to be maximized. I call this function *profit_max*. To minimize the amounts of lines in my function, I wrote the full expression of harvest rather than *Y(t)*. There are two loops: one calculates the stock after time period 1, and once calculates the discounted profit for each time period. I then have the sum of discounted profits as the number of be maximized:

```
# Specifying harvest and profit vectors
Y <- c(rep(0, times = numberOfYears))
Profit <- c(rep(0, times = numberOfYears))
# Profit maximization function
profit_max <- function(par){
for(t in 2:(numberOfYears+1)){X[t] <- X[t-1]+r*X[t-1]*(1-X[t-1]/m)-q*par[t-1]*X[t-1] }
for(t in 1:numberOfYears){Profit[t] <- (p*q*par[t]*X[t]-c*par[t])/((1+i)^t) }
return(-sum(Profit))
}
```

The function *nloptr* is meant to find the minimum of a function, as the default. In order to find the maximum of a function, I turned the problem around by asking nloptr to find the minimum of the negative of profits: *return(-sum(Profit))*.

However, this problem has constraints. As one would expect, neither harvest, effort nor the stock can be negative at any point. This already implies many constraints: 3 times the number of time periods in the problem. Moreover, it also does not make sense that the harvest would ever exceed the stock, resulting in additional constraints, one for each time period.

The second argument I need are the many constraints of this model. For now, I wrote down one constraint per line, resulting in 14 constraints. For *nloptr* to work, the left-hand side of the constraint has to be less or equal to zero. Meaning all of the 14 constraints specified in the *constraint_1* function have to take values lower or equal to zero.

```
# Constraints function
constraint_1 <- function(par) {
for(t in 2:(numberOfYears+1)){X[t] <- X[t-1]+r*X[t-1]*(1-X[t-1]/m)-q*par[t-1]*X[t-1] }
constr_1 <- q*par[1]*X[1] - X[1] # Harvest has to be lower than stock
constr_2 <- q*par[2]*X[2] - X[2]
constr_3 <- q*par[3]*X[3] - X[3]
constr_4 <- q*par[4]*X[4] - X[4]
constr_5 <- q*par[5]*X[5] - X[5]
constr_6 <- -q*par[1]*X[1] # Harvest has to be non-negative
constr_7 <- -q*par[2]*X[2]
constr_8 <- -q*par[3]*X[3]
constr_9 <- -q*par[4]*X[4]
constr_10<- -q*par[5]*X[5]
constr_11<- -X[2] # Stock has to be non-negative
constr_12<- -X[3]
constr_13<- -X[4]
constr_14<- -X[5]
constr_15<- -X[6]
return (c(constr_1 ,
constr_2 ,
constr_3 ,
constr_4 ,
constr_5 ,
constr_6 ,
constr_7 ,
constr_8 ,
constr_9 ,
constr_10,
constr_11,
constr_12,
constr_13,
constr_14,
constr_15) )
}
```

Five additional constraints are not in the constraint list. These are the fact that effort cannot be negative. Instead, I specify that these decision variables have to be non-negative by specifying it in the optimizer line the lower bound (*lb*) of the decision variables: *lb = c(rep(0, times = numberOfYears))*.

With all the arguments set up, I can now call the nloptr function to maximize profits:

```
optimal <- nloptr( x0=c(rep(10, times = numberOfYears)), # Starting values
eval_f=profit_max, # Specifying the minimization problem function
eval_g_ineq = constraint_1 , # Specifying the constraints function
lb = c(rep(0, times = numberOfYears)), # Effort cannot be lower than zero
opts=list("algorithm"="NLOPT_LN_COBYLA", # which algorithm to use
"maxeval"= 160000, # How many iterations
"xtol_rel"=1.0e-8) # Tolerance margin
)
optimal
```

The output is as follows:

```
> optimal
Call:
nloptr(x0 = c(rep(10, times = numberOfYears)), eval_f = profit_max,
lb = c(rep(0, times = numberOfYears)), eval_g_ineq = constraint_1,
opts = list(algorithm = "NLOPT_LN_COBYLA", maxeval = 160000, xtol_rel = 1e-08))
Minimization using NLopt version 2.7.1
NLopt solver status: 4 ( NLOPT_XTOL_REACHED: Optimization stopped because xtol_rel or xtol_abs (above) was reached. )
Number of Iterations....: 481
Termination conditions: maxeval: 160000 xtol_rel: 1e-08
Number of inequality constraints: 15
Number of equality constraints: 0
Optimal value of objective function: -749650988.005694
Optimal value of controls: 5.662252 1.709997 1.771068 10 10
```

The optimizer was successful at finding the effort values that maximize profit, as I wanted. The maximum value of the fishery is 749 million euros. Effort in the first time period would be 5.66, decrease to 1.7 in the second and third time period, and effort would be 10 in the last time periods.

This makes intuitive sense: the discount rate of 4% is high enough to incentivize more fishing in the first time period. Then, effort is quite low to let the stock size grow. In the last two time periods, effort is so high that the fishery will be depleted at the end of the final time period.

One can now see how different parameters change optimal effort levels, such as price, cost, or initial stock. This model can also be adapted for a larger number of time periods.

**References: **

Rahul Bhadani. Nonlinear Optimization in R using nlopt. arXiv preprint arXiv:2101.02912, 2021.

## 1 Comment