Here we present two simple examples for running an optimization algorithm on a NetLogo model with nlrx. In our example, we use the Simulated Annealing simdesign (simdesign_GenSA()). However, except for the parameter definitions in the simdesign function and the output of the function, the genetic algorithm optimization (simdesign_GenAlg()) works in the same way.
We use the Wolf Sheep Predation model from the models library to show a basic example of the optimization workflow. Example 1 shows, how a NetLogo reporter can be used as a fitness criterion for optimization. Example 2 uses a self-defined evaluation function that calculates landscape metrics that are then used as fitness criterion.
library(nlrx)
# Windows default NetLogo installation path (adjust to your needs!):
netlogopath <- file.path("C:/Program Files/NetLogo 6.0.3")
modelpath <- file.path(netlogopath, "app/models/Sample Models/Biology/Wolf Sheep Predation.nlogo")
outpath <- file.path("C:/out")
# Unix default NetLogo installation path (adjust to your needs!):
netlogopath <- file.path("/home/NetLogo 6.0.3")
modelpath <- file.path(netlogopath, "app/models/Sample Models/Biology/Wolf Sheep Predation.nlogo")
outpath <- file.path("/home/out")
nl <- nl(nlversion = "6.0.3",
nlpath = netlogopath,
modelpath = modelpath,
jvmmem = 1024)Because we want to apply an optimization algorithm, we need to define proper variable ranges. The algorithm is allowed to change the values of these parameters within these ranges in order to minimize our fitness criterion. In this example we want to use a reporter from the metrics slot for evaluating our model runs. Here we want to find a parameterization that leads to the maximum number of wolfs after 50 ticks. Because the algorithm automatically searches for minimum values, we add "1 / count wolves" to the metrics vector in order to find the maximum number of wolves.
It is also important to think about the settings for tickmetrics, runtime and evalticks. Because we only want to consider the last tick of the simulation, we set tickmetrics to “false” and runtime to 50. If more than one tick would be measured, the algorithm automatically calculates the mean value of the selected reporter. If you wish to apply other functions to aggregate temporal information into one value, you can use a self-defined evaluation function (see Example 2).
nl@experiment <- experiment(expname="wolf-sheep-GenSA1",
outpath=outpath,
repetition=1,
tickmetrics="false",
idsetup="setup",
idgo="go",
runtime=50,
metrics=c("(1 / count wolves)"),
variables = list('initial-number-sheep' = list(min=50, max=150, qfun="qunif"),
'initial-number-wolves' = list(min=50, max=150, qfun="qunif")),
constants = list("model-version" = "\"sheep-wolves-grass\"",
"grass-regrowth-time" = 30,
"sheep-gain-from-food" = 4,
"wolf-gain-from-food" = 20,
"sheep-reproduce" = 4,
"wolf-reproduce" = 5,
"show-energy?" = "false"))We use the simdesgin_GenSA() function to attach a Simulated Annealing simdesign. We select the evaluation criterion (evalcrit) by assigning the position of the reporter that we want to evaluate within the metrics vector of the experiment. In our case, there is only one reporter in the metrics vector thus we set evalcrit to use the first reporter (evalcrit = 1). The control parameter allows us to provide additional parameters for the GenSA function (see ?GenSA for details). For demonstration purposes, we set the maximum number of iterations to 20.
nl@simdesign <- simdesign_GenSA(nl,
evalcrit = 1,
nseeds = 1,
control=list(maxit = 20))For optimization simdesign, the run_nl_dyn() function lets you execute the simulations. There are some notable differences between run_nl_all() and run_nl_dyn(). First, because parameterizations depend of results from previous runs, run_nl_dyn() can not be parallelized. Second, the procedure does not automatically loop over created random seeds of the simdesign. If you want to repeat the same algorithm several times, just embed the run_nl_dyn() function in any kind of loop and iterate through the nl@simdesign@simseeds vector. Third, the output of run_nl_dyn() is reported as objects from the specific optimization procedures and not in tibble format. In order to attach these results to the nl object, the output needs to be converted to tibble format first. However, attaching optimization results to the nl does not enable any further post-processing functions of the nlrx package and is only relevant for storing results together with the nl object. This design decision was made in order to allow application of the method specific summary functions to the results of the optimization.
results <- run_nl_dyn(nl, seed = nl@simdesign@simseeds[1])The output list of the Simulated Annealing procedure contains four elements: value reports the minimum final value of the evaluation criterion. par reports the parameter settings of the final parameterisation in the same order as defined in the experiment of the nl object. trace.mat gives you detailed information on the optimization process over all iterations. counts indicates how often the optimization procedure was executed in total.
resultsIn order to store our results together with the nl object we need to attach the results to the nl object first. As explained above, we need to enframe the results as a tibble.
setsim(nl, "simoutput") <- tibble::enframe(results)
saveRDS(nl, file.path(nl@experiment@outpath, "genSA_1.rds"))library(nlrx)
# Windows default NetLogo installation path (adjust to your needs!):
netlogopath <- file.path("C:/Program Files/NetLogo 6.0.3")
modelpath <- file.path(netlogopath, "app/models/Sample Models/Biology/Wolf Sheep Predation.nlogo")
outpath <- file.path("C:/out")
# Unix default NetLogo installation path (adjust to your needs!):
netlogopath <- file.path("/home/NetLogo 6.0.3")
modelpath <- file.path(netlogopath, "app/models/Sample Models/Biology/Wolf Sheep Predation.nlogo")
outpath <- file.path("/home/out")
nl <- nl(nlversion = "6.0.3",
nlpath = netlogopath,
modelpath = modelpath,
jvmmem = 1024)Because we want to apply an optimization algorithm, we need to define proper variable ranges. The algorithm is allowed to change the values of these parameters within these ranges in order to minimize our fitness criterion. In this example we want to use a self-defined evaluation function to calculate a fitness criterion. Thus, we add the patch coordinates and patch color (as a patch class indicator) to the metrics.patches vector. We want to use spatial data to calculate the landscape edge density index of the final tick and find a parameterization that leads to the edge density. Because we only want to consider the last tick of the simulation, we set tickmetrics to “false” and runtime to 50.
nl@experiment <- experiment(expname="wolf-sheep-GenSA2",
outpath=outpath,
repetition=1,
tickmetrics="false",
idsetup="setup",
idgo="go",
runtime=50,
metrics.patches = c("pxcor", "pycor", "pcolor"),
variables = list('initial-number-sheep' = list(min=50, max=150),
'initial-number-wolves' = list(min=50, max=150)),
constants = list("model-version" = "\"sheep-wolves-grass\"",
"grass-regrowth-time" = 30,
"sheep-gain-from-food" = 4,
"wolf-gain-from-food" = 20,
"sheep-reproduce" = 4,
"wolf-reproduce" = 5,
"show-energy?" = "false"))We use the simdesgin_GenSA() function to attach a Simulated Annealing simdesign. Because we want to post-process our simulation results, we need to define an evaluation function. The evaluation function needs to accept the nl object as input and must return a single numeric value. First we load the package landscapemetrics. We then convert the spatial data to a raster format and calculate the landscape edge density index. Finally, we report only the index value of the resulting tibble.
critfun <- function(nl) {
library(landscapemetrics)
res_spat <- nl_to_raster(nl)
res_spat_raster <- res_spat$spatial.raster[[1]]
lsm <- lsm_l_ed(res_spat_raster)
crit <- lsm$value
return(crit)
}In the simdesign_GenSA() function we now provide our evaluation function (critfun) as evaluation criterion (evalcrit). The control parameter allows us to provide additional parameters for the GenSA function (see ?GenSA for details). For demonstration purposes, we set the maximum number of iterations to 20.
nl@simdesign <- simdesign_GenSA(nl,
evalcrit = critfun,
nseeds = 1,
control=list(maxit = 20))For optimization simdesign, the run_nl_dyn() function lets you execute the simulations. There are some notable differences between run_nl_all() and run_nl_dyn(). First, because parameterizations depend of results from previous runs, run_nl_dyn() can not be parallelized. Second, the procedure does not automatically loop over created random seeds of the simdesign. If you want to repeat the same algorithm several times, just embed the run_nl_dyn() function in any kind of loop and iterate through the nl@simdesign@simseeds vector. Third, the output of run_nl_dyn() is reported as objects from the specific optimization procedures and not in tibble format. In order to attach these results to the nl object, the output needs to be converted to tibble format first. However, attaching optimization results to the nl does not enable any further post-processing functions of the nlrx package and is only relevant for storing results together with the nl object. This design decision was made in order to allow application of the method specific summary functions to the results of the optimization.
results <- run_nl_dyn(nl, seed = nl@simdesign@simseeds[1])The output list of the Simulated Annealing procedure contains four elements: value reports the minimum final value of the evaluation criterion. par reports the parameter settings of the final parameterisation in the same order as defined in the experiment of the nl object. trace.mat gives you detailed information on the optimization process over all iterations. counts indicates how often the optimization procedure was executed in total.
resultsIn order to store our results together with the nl object we need to attach the results to the nl object first. As explained above, we need to enframe the results as a tibble.
setsim(nl, "simoutput") <- tibble::enframe(results)
saveRDS(nl, file.path(nl@experiment@outpath, "genSA_2.rds"))