Динамическая оптимизация (модель промысла)

Question

Я пытаюсь создать модель, в которой рыбаки выбирают уровень промыслового усилия, которое максимизирует их сумму прибыли с течением времени. Я использую простое уравнение логистического роста, и все работает нормально, но я просто не могу понять, как запустить optim, чтобы получить решение. Может ли optim найти вектор e[i], который максимизирует прибыль? Вот код, я использую:

# Optimal fishery management by choosing effort 

# Set parameters 

r = 0.1  # intrinsic growth rate 
K = 1   # carrying capacity 
q = 0.01  # catchability 
eta = 0.3  # stiffness parameter 
p = 200  # price of fish 
c = 1   # unit cost of effort 

eo = 1  # initial effort 
xo = 1  # initial biomass 
Yo = 0.01  # initial growth (meaningless) 
Ho = 0.01  # initial harvest (meaningless) 

# set time periods 

time <- seq(0,50) 

# Logitstic growth 

x <- rep(xo,length(time))     # vector for stock biomass 
Y <- rep(Yo,length(time))     # vector for growth in the stock 
H <- rep(Ho, length(time))     # vector for harvest 
e <- rep(eo, length(time))     # vector for effort 
profit <- rep(0, length(time))    # vector for profit 

for (i in 2:length(time)){ 
    x[i]=x[i-1]+Y[i-1]-H[i-1]    # stock equation 
    H[i]=q*x[i]*e[i]     # harvest equation 
    Y[i]=r*x[i]*(1-x[i]/K)    # growth equation 
    profit[i] = p*H[i]-c*e[i]    # profit equation 
    } 

totalprofit <- function(e, x){-sum(p*q*x[i]*e[i]-c*(e[i]))} 

optim(par = eo, totalprofit, x=x[i], method = "Brent", lower = 0, upper = 20) 

Answer 1

я продолжал копать и был в состоянии найти решение с помощью optimx. Самая большая проблема заключалась в том, что оптимизируемая функция должна включать все системные уравнения. Прошу прощения за ответ на мой вопрос.

# Optimal fishery management by choosing effort 

# Set parameters 

r = 0.1  # intrinsic growth rate 
K = 1   # carrying capacity 
q = 0.01  # catchability 
eta = 0.3  # stiffness parameter 
p = 200  # price of fish 
c = 1   # unit cost of effort 

eo = 1  # initial effort 
xo = 1  # initial biomass 
Yo = 0.01  # initial growth (meaningless) 
Ho = 0.01  # initial harvest (meaningless) 

# set time periods 

time <- seq(0,100) 

# Logitstic growth 

x <- rep(xo,length(time))     # vector for stock biomass 
Y <- rep(Yo,length(time))     # vector for growth in the stock 
H <- rep(Ho, length(time))     # vector for harvest 
e <- rep(eo, length(time))     # vector for effort 
profit <- rep(1, length(time))    # vector for profit 

# Create function to be optimized 
# function contains all equations 

totalprofit <- function (e, npar = TRUE, print = TRUE) { 
    for (i in 2:length(time)){ 
    x[i]=x[i-1]+Y[i-1]-H[i-1]    # stock equation 
    H[i]=q*x[i]*e[i]     # harvest equation 
    Y[i]=r*x[i]*(1-x[i]/K)    # growth equation 
    profit[i] = (p*q*x[i]*e[i])-(c*e[i]) # profit equation 
    } 
    result <- sum(profit) 
    return(result) 
    } 

# Optimize system by choosing effort 

maxfish <- optimx(par = e, fn = totalprofit, lower = 0, upper = 100, 
    method = c("L-BFGS-B"), 
    control = list(maximize = TRUE)) 

emax <- coef(maxfish) 

# Run system again to see performance 

for (i in 2:length(time)){ 
    x[i]=x[i-1]+Y[i-1]-H[i-1]    # stock equation 
    H[i]=q*x[i]*emax[i]     # harvest equation 
    Y[i]=r*x[i]*(1-x[i]/K)    # growth equation 
    profit[i] = (p*q*x[i]*emax[i])-(c*emax[i]) # profit equation 
    } 

plot (time, x, type='b') 

plot (x,emax, type='b') 

plot (time, emax, type='b') 

plot (time, H, type='b') 

plot (time, profit, type='b')