I have a problem with the calibration of a simple crop model, when using a function (composed by different elements) and the optim package for optimizing 2 parameters of the crop model ( parameters vector). For each year, the function adds a row in the data.frame when an accumulation index exceeds a threshold ( parameters[2] ).
The function model_func (attached below) works and gives the results adequately (when a set of parameters previously calibrated are used), however when I want to optimize the parameters within a range, R gives the error:
Error in optim(par = parameters, x2 = x2, fn = model_func, upper = c(parameters[1], : (list) object cannot be coerced to type 'double'
At the beginning, I inizialized model_db as a vector() which created a matrix when the results were consecutively added with the use of rbind() . In that case, I thought that the error was done because I was working with matrices. However, it also appeared when I have worked directly with dataframes.
In theory, both columns of the resulting summarizing table are numeric. Is it wrong?
Thaks,
Sergi
Simplifying, the code is:
model_func <- function(x2, parameters)
{
for (m in unique(x2[,1])) ## Where x2 is the daily meteo and the first column defines locations
{
x <- subset(x2, x2[,1] == m)
##### Initializing data.frame and indices
model_db <- data.frame(matrix(ncol = 2, nrow = 0))
colnames(model_db) <- c("Year", "Simul")
accum_model <- 0
cc_model <- 0
###### Starting model
for(d in 1:(length(x[,1]) - 1))
{
###### Extracting of the weather database
temp_year <- x[d, 2]
temp_doy <- x[d, 3]
temp_tmean <- x[d, 4]
# Daily model accumulation
accum_d_model <- max(temp_tmean - parameters[3], 0)
if (temp_doy <= parameters[1])
{
accum_model <- accum_d_model
cc_model <- 0
} else {
accum_model <- accum_model + accum_d_model
}
if(accum_model >= parameters[2] & cc_model == 0)
{
cc_model <- 1
model_db <- model_db %>% add_row(Year = temp_year, Simul = temp_doy)
}
}
tab_model_db <- if(m == unique(x2[,1])[1])
{cbind(m, model_db[-1,])} else {rbind(tab_model_db, cbind(m, model_db[-1,]))}
}
tab_model_db <- data.frame(tab_model_db)
names(tab_model_db) <- c("Site", "Year", "Simul")
return(tab_model_db)
}
### Optimizer
optim(
par = parameters,
x2 = x2,
fn = model_func,
upper = c(parameters[1], parameters[2:3] + abs(parameters[2:3] * 0.3)),
lower = c(parameters[1], parameters[2:3] - abs(parameters[2:3] * 0.3)),
method = 'L-BFGS-B'
)$par```