This example is using NetLogo Flocking model (Wilensky, 1998) to demonstrate parameter fitting with simulated annealing.
There are many R packages for solving optimization problems (see CRAN Task View). This example uses GenSA function from GenSA
package (Xiang et al., 2013). See also Thiele, Kurth & Grimm (2014) chapter 2.31 Simulated annealing.
Experiment definition is the same as in L-BFGS-B Optimization example:
experiment <- nl_experiment(
model_file = "models/Sample Models/Biology/Flocking.nlogo",
setup_commands = c("setup", "repeat 100 [go]"),
iterations = 5,
param_values = list(
world_size = 50,
population = 80,
vision = 6,
min_separation = seq(from = 0, to = 4, by = 0.5),
max_align_turn = seq(from = 0, to = 20, by = 2.5)
),
mapping = c(
min_separation = "minimum-separation",
max_align_turn = "max-align-turn"),
step_measures = measures(
converged = "1 -
(standard-deviation [dx] of turtles +
standard-deviation [dy] of turtles) / 2",
mean_crowding =
"mean [count flockmates + 1] of turtles"
),
eval_criteria = criteria(
c_converged = mean(step$converged),
c_mcrowding = mean(step$mean_crowding)
),
repetitions = 10, # repeat simulations 10 times
eval_aggregate_fun = mean, # aggregate over repetitions
eval_mutate = criteria( # evaluation criterium
eval_value =
sqrt((c_mcrowding - 8)^2 + 400*(c_converged - 1)^2)
)
)
In this example the nl_eval_run
function is used as evaluation function instead of evaluating all parameter sets with nl_run
.
Note that this kind of evaluation requires started NetLogo instance. User have to take care to initialize NetLogo and load the model before optimization begins and close NetLogo when it is no longer needed (see nl_eval_run
in package documentation).
library(GenSA)
cl <- nl_eval_init(experiment, parallel = TRUE)
trace <- nl_eval_tracer(verbose = FALSE)
param_range <- nl_get_param_range(experiment)
set.seed(42)
result <- GenSA(
par = param_range$upper,
fn = nl_eval_run,
experiment = experiment,
criteria = "eval_value",
call_back = trace$add,
parallel = TRUE, cluster = cl,
param_names = names(param_range$lower),
control=list(max.time = 60*5), # let it run for 5 minutes
lower = param_range$lower,
upper = param_range$upper
)
nl_eval_close(parallel = TRUE, cl)
final <- setNames(data.frame(t(result$par)), names(param_range$lower))
final$result <- result$value
htmlTable::htmlTable(round(final,1))
min_separation | max_align_turn | result | |
---|---|---|---|
1 | 1.9 | 14.4 | 4.4 |
library(ggplot2)
dat <- as.data.frame(result$trace.mat)
ggplot(dat, aes(x = nb.steps, y =current.minimum)) +
geom_point(data = dat,
aes(x = nb.steps, y = function.value),
position = position_jitter(height = 0, width = 0.05),
alpha = 0.5)+
geom_step(alpha = 0.3) +
theme_minimal() +
ylab("Values, current minimum")
dat <- trace$get()
ggplot(dat, aes(x = min_separation, y = max_align_turn, color = result)) +
geom_point(size = 3, alpha = 0.8) +
geom_point(data = final , size = 10, shape = 8, color = "red") +
theme_minimal() +
#theme(legend.position = "none") +
coord_fixed(4/20)