动态优化(渔业模型)

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')