This summary is part of the Brave New Green Algorithm series, that investigates energy consumption of the Brave New Algorithm and, by extension, by other population-based optimization algorithms.

Introduction

The LION paper analyzes the Brave New Algorithm (BNA), a caste-stratified population-based stochastic optimization algorithm, measuring fitness and energy jointly (Merelo-Guervós et al. 2026; Merelo, Merelo, and Garcı́a-Valdez 2022). The optimization task in this study is the BBOB Sphere function (Hansen et al. 2010).

In this context, baseline refers to runs that stop after generating the initial population (same code path, no evolutionary loop), computed per (dimension, population) configuration and subtracted from workload runs (Merelo-Guervós et al. 2026). This page condenses the core paper findings into four visual messages:

  1. Baseline energy is noisy and must be discounted carefully.
  2. Within the same Julia version, max_gens = 10 usually saves energy.
  3. Lower energy and better fitness do not always happen together.
  4. Upgrading Julia can change the entire energy profile, even with the same algorithm.

Data preparation with energyR

# LION (Julia 1.11.8)
lion_baseline <- load_bna_csv("data/lion-1.11.8-baseline.csv")
lion_summary <- summarize_baseline(lion_baseline)

lion_workload <- load_bna_csv("data/lion-1.11.8-bna-fix-rand.csv", drop_baseline_cols = FALSE)
lion_workload <- prepare_workload(lion_summary, lion_workload, label_col = "dim_pop")
lion_workload$julia <- "1.11.8"

# EvoApps baseline/workload used in the paper as the previous platform (Julia 1.11.7)
prev_baseline <- load_bna_csv("data/evoapps-1.11.7-baseline-bna-baseline-16-Oct-11-08-20.csv")
prev_summary <- summarize_baseline(prev_baseline)

prev_workload <- load_bna_csv("data/evoapps-1.11.7-fix-rand-bna-fix-rand-25-Oct-11-06-07.csv", drop_baseline_cols = FALSE)
prev_workload <- prepare_workload(prev_summary, prev_workload, label_col = "dim_pop")
prev_workload$julia <- "1.11.7"

1) Baseline is variable, so subtraction matters

Before looking at algorithmic effects, we need to separate them from runtime and platform overhead. The baseline plot shows non-negligible dispersion, and that dispersion depends on dimension and population size.

This is one of the main methodological conclusions of the paper: comparing raw Joules between configurations is misleading, so baseline-corrected deltas are the defensible unit for comparisons (Merelo-Guervós et al. 2026; Cotta and Martı́nez-Cruz 2024).

plot_data <- lion_baseline %>%
  mutate(
    dimension = as.factor(dimension),
    population_size = as.factor(population_size)
  )

ggplot(plot_data, aes(x = dimension, y = PKG, fill = population_size)) +
  geom_violin(alpha = 0.6, trim = FALSE) +
  theme_minimal(base_size = 13) +
  labs(
    title = "Baseline runs have visible spread",
    subtitle = "This is why LION computes delta energy against per-configuration baseline summaries",
    x = "Dimension",
    y = "PKG energy (J)",
    fill = "Population"
  )
Baseline PKG distribution by dimension and population size (LION, Julia 1.11.8).

Baseline PKG distribution by dimension and population size (LION, Julia 1.11.8).

2) Main LION configuration effect: stopping earlier usually saves energy

With baseline correction in place, the effect of max_gens becomes clearer: in most parameter combinations, max_gens = 10 shifts the distribution downward relative to larger values.

The practical conclusion is not “always use fewer generations.” The paper shows that this setting is an energy-lean default, but it must be interpreted jointly with solution quality and with problem-specific accuracy requirements.

lion_for_plot <- lion_workload %>%
  mutate(
    work = paste0("max_gens=", max_gens),
    dim_pop = factor(dim_pop)
  )

plot_delta_energy(
  lion_for_plot,
  geom = "boxplot",
  facet_col = "dim_pop",
  title = "Delta energy by stopping criterion (LION / Julia 1.11.8)"
)
Delta energy in LION workload runs, grouped by max_gens and faceted by parameter combination.

Delta energy in LION workload runs, grouped by max_gens and faceted by parameter combination.

3) Energy vs fitness is a trade-off, not a single objective

This plot makes explicit what the paper discusses in detail: there is no universal setting that simultaneously optimizes energy and fitness in every case. Points do not collapse into a single dominant frontier.

Operationally, that means configuration should be policy-driven. If your context prioritizes lower energy, choose from the lower-energy clusters; if it prioritizes better fitness improvements, accept the additional energy budget and select accordingly (Merelo-Guervós et al. 2026).

lion_tradeoff <- lion_for_plot %>%
  filter(is.finite(log_diff), diff_fitness > 0) %>%
  mutate(
    work = factor(paste0("max_gens=", max_gens)),
    dimension = factor(dimension),
    population_size = factor(population_size)
  )

ggplot(lion_tradeoff, aes(x = delta_PKG, y = log_diff, color = work, shape = population_size)) +
  geom_point(alpha = 0.55, size = 2) +
  facet_wrap(~dimension, nrow = 1, scales = "free_x") +
  theme_minimal(base_size = 13) +
  labs(
    title = "LION: energy versus fitness by dimension",
    subtitle = "Faceting by dimension clarifies the per-configuration max_gens effect",
    x = "Delta PKG energy (J)",
    y = "log10 fitness improvement",
    color = "max_gens",
    shape = "Population"
  )
Energy-versus-fitness view for LION runs, faceted by problem dimension and grouped by max_gens/population.

Energy-versus-fitness view for LION runs, faceted by problem dimension and grouped by max_gens/population.

4) Platform drift: Julia version alone changes the energy baseline-corrected profile

One of the strongest discussion points in the paper is reproducibility under software platform evolution. The comparison across Julia versions shows that runtime/platform updates can materially shift measured energy, even when the algorithmic procedure is unchanged (Merelo-Guervós et al. 2026).

So, any claim of “configuration A is greener than B” should be version-scoped and re-validated after upgrades. In practice, rerunning baseline + workload after runtime changes is part of the protocol, not a nice-to-have (Cotta and Martı́nez-Cruz 2024; Merelo-Guervós et al. 2026).

version_compare <- bind_rows(
  prev_workload %>% mutate(source = "EvoApps setup"),
  lion_workload %>% mutate(source = "LION setup")
) %>%
  mutate(
    population_size = as.factor(population_size),
    dimension = as.factor(dimension)
  )

ggplot(version_compare, aes(x = julia, y = delta_PKG, fill = julia)) +
  geom_boxplot(alpha = 0.8, outlier.alpha = 0.2) +
  facet_grid(dimension ~ population_size) +
  theme_minimal(base_size = 13) +
  labs(
    title = "Julia version impacts observed energy",
    subtitle = "Consistent with the LION paper: software platform changes can dominate small tuning effects",
    x = "Julia version",
    y = "Delta PKG energy (J)",
    fill = "Julia"
  )
Comparing baseline-corrected workload energy for similar BNA settings across Julia versions used in the paper.

Comparing baseline-corrected workload energy for similar BNA settings across Julia versions used in the paper.

Compact configuration summary

# 1. Create a table of pairwise Wilcoxon tests
pairwise_comparisons <- lion_workload %>%
  filter(is.finite(delta_PKG)) %>%
  group_by(dimension, population_size) %>%
  group_modify(~{
    gens <- sort(unique(.x$max_gens))
    
    # If there's nothing to compare, return an empty template
    if (length(gens) < 2) {
      return(data.frame(
        max_gens_A = numeric(0), 
        max_gens_B = numeric(0), 
        p_value = numeric(0), 
        sig_mark = character(0)
      ))
    }
    
    # Generate all unique pairs of max_gens
    pair_tests <- utils::combn(gens, 2, simplify = FALSE)
    
    # Calculate the test for each pair and keep them in separate columns
    res <- lapply(pair_tests, function(pair) {
      x <- .x$delta_PKG[.x$max_gens == pair[1]]
      y <- .x$delta_PKG[.x$max_gens == pair[2]]
      p_value <- stats::wilcox.test(x, y, exact = FALSE)$p.value
      
      data.frame(
        max_gens_A = pair[1],
        max_gens_B = pair[2],
        p_value = p_value,
        sig_mark = if (is.finite(p_value) && p_value < 0.05) "*" else ""
      )
    })
    
    dplyr::bind_rows(res)
  }) %>%
  ungroup()

# 2. Compute the base summary metrics (assuming create_summary returns one row per config)
base_summary <- create_summary(lion_workload)

# 3. Join the summary metrics for both A and B to the pairwise table
lion_summary_table <- pairwise_comparisons %>%
  # Join metrics for max_gens_A
  left_join(
    base_summary %>% select(dimension, population_size, max_gens, trimmed_mean_delta_PKG, trimmed_mean_energy_per_evaluation),
    by = c("dimension", "population_size", "max_gens_A" = "max_gens")
  ) %>%
  rename(
    delta_PKG_A = trimmed_mean_delta_PKG,
    energy_eval_A = trimmed_mean_energy_per_evaluation
  ) %>%
  # Join metrics for max_gens_B
  left_join(
    base_summary %>% select(dimension, population_size, max_gens, trimmed_mean_delta_PKG, trimmed_mean_energy_per_evaluation),
    by = c("dimension", "population_size", "max_gens_B" = "max_gens")
  ) %>%
  rename(
    delta_PKG_B = trimmed_mean_delta_PKG,
    energy_eval_B = trimmed_mean_energy_per_evaluation
  ) %>%
  # Clean up and arrange
  mutate(
    delta_PKG_A = round(delta_PKG_A, 2),
    delta_PKG_B = round(delta_PKG_B, 2),
    energy_eval_A = signif(energy_eval_A, 3),
    energy_eval_B = signif(energy_eval_B, 3),
    p_value = signif(p_value, 3)
  ) %>%
  # Select final columns in a logical reading order
  select(
    dimension, population_size, 
    max_gens_A, delta_PKG_A, energy_eval_A,
    max_gens_B, delta_PKG_B, energy_eval_B, 
    p_value, sig_mark
  ) %>%
  arrange(dimension, population_size, p_value)

# 4. Render the table
knitr::kable(
  lion_summary_table,
  col.names = c(
    "Dimension", "Pop. Size", 
    "Max Gens A", "Delta PKG A (J)", "J/eval A",
    "Max Gens B", "Delta PKG B (J)", "J/eval B", 
    "p-value", "Sig."
  ),
  caption = "Pairwise comparisons of energy-only efficiency between max_gens settings. Lower values are better. * marks a statistically significant difference (p < 0.05) in delta-PKG between the specific pair in the row."
)
Pairwise comparisons of energy-only efficiency between max_gens settings. Lower values are better. * marks a statistically significant difference (p < 0.05) in delta-PKG between the specific pair in the row.
Dimension Pop. Size Max Gens A Delta PKG A (J) J/eval A Max Gens B Delta PKG B (J) J/eval B p-value Sig.
3 200 10 141.36 0.00357 25 137.55 0.001620 0.1580
3 400 10 140.42 0.00178 25 137.38 0.000800 0.1470
5 200 10 137.03 0.00362 25 134.23 0.001720 0.2970
5 400 10 128.97 0.00165 25 133.64 0.000844 0.0292 *

The table is useful as an at-a-glance ranking, but it should be read as a screening tool. Final decisions should combine this with fitness outcomes and with the platform/version context highlighted above. We can see, however, that the only significant differences in energy expenses for different values of max_gens are in the top rows, with dimension 5 and population 400.

Expanded analysis: multi-stage model on the final LION dataset

This section applies the later-paper multi-stage modeling workflow directly to the final LION workload file data/lion-1.11.8-bna-fix-rand.csv (Merelo-Guervós et al. 2026; Cotta and Martı́nez-Cruz 2024):

  1. model runtime (seconds) from algorithm covariates (dimension, population_size, generations, evaluations);
  2. residualize runtime;
  3. model energy (PKG) on the same covariates plus residualized runtime.
lion_model_data <- read.csv("data/lion-1.11.8-bna-fix-rand.csv") %>%
  mutate(
    max_gens = factor(max_gens),
    alpha = factor(alpha),
    population_size = factor(population_size),
    dimension = factor(dimension)
  )

time_model <- glm(
  seconds ~ dimension*population_size*generations*evaluations,
  data = lion_model_data
)

lion_model_data$residual_seconds <- residuals(time_model)

These residual_seconds (click on “Show”) isolate operating-context and implementation influences, not affected by the algorithm parameters.

energy_model <- glm(
  PKG ~ dimension + population_size + generations + evaluations + residual_seconds,
  data = lion_model_data
)

energy_coef <- summary(energy_model)$coefficients
energy_coef_table <- data.frame(
  term = rownames(energy_coef),
  estimate = energy_coef[, "Estimate"],
  std_error = energy_coef[, "Std. Error"],
  statistic = energy_coef[, "t value"],
  p_value = energy_coef[, "Pr(>|t|)"],
  row.names = NULL
) %>%
  filter(term != "(Intercept)") %>%
  arrange(p_value)

energy_sig_table <- energy_coef_table %>%
  filter(p_value < 0.10) %>%
  mutate(
    estimate = signif(estimate, 4),
    std_error = signif(std_error, 3),
    statistic = signif(statistic, 3),
    p_value = signif(p_value, 3)
  )

if (nrow(energy_sig_table) == 0) {
  energy_sig_table <- data.frame(
    term = "No coefficient with p < 0.05",
    estimate = NA, std_error = NA, statistic = NA, p_value = NA
  )
}

knitr::kable(
  energy_sig_table,
  col.names = c("Term", "Estimate", "Std. error", "t", "p-value"),
  caption = "Stage 2 (energy model): statistically significant (90% level) coefficients influencing PKG."
)
Stage 2 (energy model): statistically significant (90% level) coefficients influencing PKG.
Term Estimate Std. error t p-value
evaluations 0.000055 3.06e-05 1.80 0.0739
generations -0.157000 9.13e-02 -1.72 0.0868

What you see is what you get: only two of the coefficients are significant, and at the 10% level; most changes that can be attributed to algorithm parameters are drowned by the operating context. This is what led us to improve the experimental protocol, something you can check in this paper published at CEC.

energy_anova <- anova(energy_model)
energy_variance <- data.frame(
  term = rownames(energy_anova),
  deviance = energy_anova$"Deviance",
  p_value = energy_anova$"Pr(>F)",
  row.names = NULL
) %>%
  filter(!is.na(deviance), term != "Residuals") %>%
  mutate(
    variance_pct = 100 * deviance / sum(deviance),
    term = factor(term, levels = term[order(variance_pct)])
  )

ggplot(energy_variance, aes(x = term, y = variance_pct, fill = term)) +
  geom_col(width = 0.72, alpha = 0.9) +
  geom_text(
    aes(label = paste0(round(variance_pct, 1), "%")),
    hjust = -0.1,
    size = 3.8
  ) +
  coord_flip() +
  theme_minimal(base_size = 13) +
  theme(legend.position = "none") +
  labs(
    title = "Energy variance explained by each modeled term",
    subtitle = "Computed from ANOVA sum-of-squares in the Stage 2 model",
    x = "Model term",
    y = "Share of explained variance (%)"
  )
Variance decomposition for the multi-stage energy model (percent of model sum of squares by term).

Variance decomposition for the multi-stage energy model (percent of model sum of squares by term). 

energy_variance_table <- energy_variance %>%
  mutate(
    variance_pct = round(variance_pct, 2),
    p_value = signif(p_value, 3)
  ) %>%
  arrange(desc(variance_pct)) %>%
  select(term, variance_pct, p_value)

knitr::kable(
  energy_variance_table,
  col.names = c("Term", "Explained variance (%)", "p-value"),
  caption = "ANOVA-based variance share by term for the multi-stage energy model."
)
ANOVA-based variance share by term for the multi-stage energy model.
Term Explained variance (%) p-value
population_size 47.02 0.00669
evaluations 20.24 0.07390
residual_seconds 16.93 0.10200
dimension 15.81 0.11400
generations 0.01 0.97100

Again, as you see in the table, only one of the coefficients (population size) explains variance at a 95% significance level, and one (evaluations) at a 90% level.

Discussion highlights from the paper

  • Measurement methodology matters as much as algorithm design: baseline correction and careful sampling choices directly influence the conclusions you can trust (Merelo-Guervós et al. 2026; Cotta and Martı́nez-Cruz 2024).

  • Energy is a first-class optimization objective: if it is not explicitly included in analysis, tuning decisions can look good on fitness while being unnecessarily expensive energetically (Merelo-Guervós et al. 2026).

  • Results are contextual, not absolute: hardware state, runtime version, and execution environment can dominate small parameter effects (Merelo-Guervós et al. 2026; Cotta and Martı́nez-Cruz 2024).

  • Reproducibility requires protocol, not just code: to compare studies over time, rerun baseline and workload suites under each relevant software/hardware change (Merelo-Guervós et al. 2026).

  • Actionable takeaway: report both fitness and energy, and present trade-off-aware recommendations instead of one-number “best setting” claims (Merelo-Guervós et al. 2026; Merelo, Merelo, and Garcı́a-Valdez 2022).

See also

You can go ahead to this expanded version of the paper, using the same data, which applies statistical models to understand what these measurements reveal about energy consumption (and what they don’t).

Questions? Suggestions?

You can use GitHub issues to make them. PRs encouraged and welcome!

References

Cotta, Carlos, and Jesús Martı́nez-Cruz. 2024. “Energy Consumption Analysis of Batch Runs of Evolutionary Algorithms.” In Proceedings of the Genetic and Evolutionary Computation Conference Companion, 87–88. GECCO ’24 Companion. New York, NY, USA: Association for Computing Machinery. https://doi.org/10.1145/3638530.3664093.
Hansen, Nikolaus, Anne Auger, Raymond Ros, Steffen Finck, and Petr Pošı́k. 2010. “Comparing Results of 31 Algorithms from the Black-Box Optimization Benchmarking BBOB-2009.” In Proceedings of the 12th Annual Conference Companion on Genetic and Evolutionary Computation, 1689–96.
Merelo, Cecilia, Juan J Merelo, and Mario Garcı́a-Valdez. 2022. “A Brave New Algorithm to Maintain the Exploration/Exploitation Balance.” In New Perspectives on Hybrid Intelligent System Design Based on Fuzzy Logic, Neural Networks and Metaheuristics, 305–16. Springer.
Merelo-Guervós, Juan J., Cecilia Merelo-Molina, Pablo García-Sánchez, and Mario García-Valdez. 2026. “Is There a (Carbon-) Free Lunch? Energy/Performance Tradeoffs in Population-Based Metaheuristics.”