0

我试图建立一个模型，让渔民选择随着时间的推移最大化他们的利润总和的捕鱼努力水平。我正在使用一个简单的逻辑增长方程，一切似乎都很好，除了我不知道如何运行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 )
4

1 回答 1

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')
于 2013-09-30T17:48:33.960 回答