Multi-objective forest restoration planning in a highly productive landscape"

In multiscape: Multi-Action, Multi-Objective Spatial Planning

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.path = "integrated-ecosystem-services-",
  fig.width = 7,
  fig.height = 5,
  fig.retina = 2,
  out.width = "100%",
  message = FALSE,
  warning = FALSE,
  error = FALSE
)

objective_aliases <- c(
  "cost",
  "rodent_control",
  "carbon_stock",
  "water_yield",
  "sediment_retention"
)

Overview

Restoration planning is especially difficult in highly productive landscapes. Most land is already committed to agriculture, forestry, or other economically valuable uses, so restoration cannot be treated as though it were being placed on vacant land. Every new restoration site competes with an existing use, requires public or private expenditure, and may alter areas that already provide important ecosystem services.

At the same time, ecosystem-service provision is often uneven across the landscape. Some locations already store substantial carbon, regulate water flows, retain sediment, or support biological control, whereas other locations provide comparatively little of these functions. In this example, the planning goal is therefore not to restore the sites with the highest current service values. Instead, it is to identify eligible sites with relatively low baseline levels of ecosystem services. These low-service locations are interpreted as areas of current ecological deficit where restoration may be most justified, while avoiding the displacement of functions that are already being supplied elsewhere.

This creates a genuine many-objective problem. A site with low carbon stock may still have high seasonal water yield, and a site with low sediment retention may be important for rodent and lagomorph control. Moreover, the least expensive sites are not necessarily those with the largest ecosystem-service deficits. Collapsing the four services into a single index would hide these conflicts and would require fixed exchange weights before decision makers have seen their consequences.

This vignette uses multiscape to extend an existing conservation network with a fixed amount of forest restoration. Existing conservation areas are treated as prior commitments and remain unchanged. The optimizer decides where to add restoration within the surrounding productive matrix while balancing five separate objectives: implementation cost and the deficits to common service ceilings in the selected restoration portfolio.

The example addresses the following planning question:

Which additional sites should be restored to complete the landscape commitment at low cost, while directing restoration towards locations with comparatively low current provision of multiple ecosystem services?

The analysis produces a set of alternative spatial plans rather than one predefined answer. These alternatives show how the restoration pattern changes when greater emphasis is placed on different ecosystem-service deficits.

Planning formulation

Let (i \in \mathcal{I}_R) index planning units that are eligible for forest restoration. The binary decision variable (r_i) equals one when unit (i) is selected for restoration and zero otherwise. The existing conservation network is fixed in every solution. It provides the spatial foundation of the plan, but the decisions compared in this vignette concern only the additional restoration units.

The total restoration cost is

[ C(r)=\sum_{i \in \mathcal{I}_R}c_i^R r_i, ]

where (c_i^R) is the local cost of restoring unit (i).

For each ecosystem service (f), define a common ceiling (M_f) as the largest reference amount observed in the eligible restoration units. The effect assigned to restoration is the deficit to that ceiling:

[ \Delta_{if}=M_f-q_{if},\qquad B_f(r)=\sum_{i\in\mathcal{I}R}(M_f-q{if})r_i. ]

Here (q_{if}) is the reference amount currently provided by unit (i). The four service objectives maximize (B_f), while implementation cost is minimized. A larger benefit means that restoration is allocated to a stronger current deficit relative to the common ceiling.

The complete objective vector is

[ \left{\min C(r),\ \max B_{\mathrm{control}}(r), \max B_{\mathrm{carbon}}(r),\ \max B_{\mathrm{water}}(r), \max B_{\mathrm{sediment}}(r)\right}. ]

Every plan restores exactly (K) equal-area units. Consequently,

[ B_f(r)=K M_f-\sum_i q_{if}r_i. ]

The first term is constant across plans. Maximizing this deficit therefore preserves the former criterion of minimizing current service provision in the restoration portfolio. This equivalence requires fixed effort, a common ceiling within each feature, and one restoration action per unit. It does not apply automatically to varying restoration effort or concurrent restoration actions.

The ceiling is an explicit accounting assumption for deficit-based prioritization. It is not an estimated ecological response to restoration. A low current value does not establish that a site can reach the ceiling or will respond more strongly than another site. Applications with estimated outcomes should supply those estimates to add_effects() instead.

The four services remain separate objectives because their spatial patterns differ; no exchange weights between services are imposed in advance.

Landscape and ecosystem-service data

The example represents a highly productive landscape divided into 30,496 planning units. Within this matrix, only a subset of units is available for new forest restoration, while other units are already conserved or excluded from intervention.

Four normalized raster layers describe current ecosystem-service provision: rodent and lagomorph control, carbon stock, seasonal water yield, and sediment retention. Raster-cell values range from 0 to 1 within each layer. create_problem() sums these values inside each planning unit, so the stored reference amounts can exceed 1. Higher values indicate stronger current provision of the corresponding service; lower values indicate areas where that service is comparatively weak and where restoration may help address an existing ecological deficit.

library(multiscape)
library(dplyr)
library(ggplot2)
library(terra)

service_data <- load_ecosystem_services()
planning_units <- service_data$planning_units
ecosystem_services <- service_data$feature_raster

service_labels <- c(
  L1_control_inv = "Rodent and lagomorph control",
  L1_stock_inv = "Carbon stock",
  L1_rend_inv = "Seasonal water yield",
  L1_reten_inv = "Sediment retention"
)

nrow(planning_units)
names(ecosystem_services)

The maps show that the four deficits do not occur in the same places. An area with low carbon stock may not have low water yield, and an area with weak sediment retention may still contribute substantially to biological control. This lack of spatial congruence is the central reason for retaining four separate objectives. A restoration plan that performs well for one service can still redirect intervention towards high-value areas for another.

names(ecosystem_services) <- unname(service_labels[names(ecosystem_services)])
terra::plot(ecosystem_services)
names(ecosystem_services) <- names(service_labels)

Define the restoration decision context

The planning-unit attributes divide the landscape into three management roles. Units marked locked_in constitute the existing conservation network and are included in every plan. Units that are neither locked in nor locked out form the candidate pool for additional forest restoration. Units marked locked_out are unavailable because of land-use, institutional, or other planning restrictions represented in the example.

The action table includes both conservation and restoration so that the final maps display the complete intervention system. However, only restoration is a variable decision. The conservation network is inherited from previous policy commitments and remains spatially fixed across all solutions.

planning_attributes <- planning_units |>
  sf::st_drop_geometry() |>
  select(id, cost, area, locked_in, locked_out)

actions <- data.frame(
  id = c("conservation", "restoration"),
  name = c("Conservation", "Restoration")
)

conservation_units <- planning_attributes |>
  filter(locked_in) |>
  transmute(pu = id, action = "conservation")

restoration_units <- planning_attributes |>
  filter(!locked_in, !locked_out) |>
  transmute(pu = id, action = "restoration")

feasible_actions <- bind_rows(
  conservation_units,
  restoration_units
)

action_costs <- feasible_actions |>
  left_join(
    planning_attributes |> select(id, cost),
    by = c("pu" = "id")
  )

feasible_actions |> count(action)

The problem is therefore not a competition between conservation and restoration inside the same planning unit. It is a network-completion problem: given an existing conservation foundation, which additional units should be restored to meet the remaining area commitment?

This distinction also clarifies the role of the many-objective analysis. AUGMECON does not decide whether existing conservation should be replaced. It compares alternative ways of placing the same amount of new restoration across the eligible productive matrix.

Build the planning problem

create_problem() links the planning-unit geometries to the four ecosystem- service layers and stores the baseline amount of each service in every unit. add_actions() then registers the feasible conservation and restoration assignments and their associated local costs.

The planning-unit cost field remains available in the problem object for inspection. The optimization objective, however, uses the explicit action costs defined for restoration. This avoids counting the same cost component twice.

problem <- create_problem(
  pu = planning_units,
  features = ecosystem_services,
  cost = "cost"
) |>
  add_actions(
    actions = actions,
    include_pairs = feasible_actions,
    cost = action_costs
  )

The landscape-level commitment is set at 20% of all 30,496 planning units. The existing conservation network already contributes 3,988 units towards this commitment. The remaining gap is therefore 2,112 units, which must be supplied through new forest restoration.

Every solution selects exactly these 2,112 restoration units. Holding restoration area constant is essential for interpretation: differences in cost or ecosystem-service objectives arise from where restoration is allocated, not from one solution restoring more land than another. Because all planning units have the same area, an equality constraint with a tolerance smaller than one unit enforces the same restoration effort in every plan.

No feature-representation target is imposed. The policy requirement concerns the amount of land restored, while the four ecosystem-service layers are used to compare alternative locations for that fixed effort.

unit_area <- stats::median(planning_units$area)
stopifnot(all(is.finite(planning_units$area)),
          all(abs(planning_units$area - unit_area) < 1e-8))
total_commitment_units <- ceiling(0.20 * nrow(planning_units))
restoration_target_units <- total_commitment_units - nrow(conservation_units)
restoration_target_area <- restoration_target_units * unit_area

problem <- problem |>
  add_constraint_area(
    area = restoration_target_area,
    sense = "equal",
    tolerance = 0.25 * unit_area,
    actions = "restoration"
  ) |>
  add_constraint_locked_planning_units(
    locked_in = "locked_in",
    locked_out = "locked_out"
  )

Express restoration effects as deficits

For each feature, use the maximum reference amount among eligible restoration units as a common ceiling. Supplying this ceiling as outcome makes add_effects() compute the signed effect as ceiling minus local reference. The ceilings share the units of the aggregated reference amounts.

restoration_ceilings <- problem$data$dist_features |>
  filter(pu %in% restoration_units$pu) |>
  group_by(feature) |>
  summarise(outcome = max(amount), .groups = "drop")

restoration_effects <- restoration_ceilings |>
  mutate(action = "restoration")

problem <- problem |>
  add_effects(restoration_effects)

restoration_ceilings

Only restoration contributes effects in this accounting scope. The fixed conservation network remains visible in the maps but contributes no restoration benefit. The common ceiling makes every restoration effect non-negative here; max_benefit() also supports signed effects when other outcomes imply losses.

Register five restoration objectives

The model registers five atomic objectives, each of which retains a distinct interpretation.

The first objective minimizes the explicit cost of the selected restoration actions. The other four objectives maximize the deficit to the common ceiling for each ecosystem service in the restoration portfolio. In practical terms, they encourage the optimizer to direct restoration towards sites with low baseline carbon stock, low seasonal water yield, low sediment retention, or low rodent and lagomorph control.

Each ecosystem-service objective is restricted to the restoration action. The fixed conservation network is therefore visible in the maps but does not add a constant value to every solution. This makes the reported objective values directly comparable as properties of the additional restoration portfolio.

Keeping the four services separate is deliberate. A site selected because it has low carbon stock may still have high water yield, so no single restoration portfolio is expected to maximize all four deficit totals simultaneously.

problem <- problem |>
  add_objective_min_cost(
    alias = "cost",
    include_pu_cost = FALSE,
    include_action_cost = TRUE,
    actions = "restoration"
  ) |>
  add_objective_max_effect(
    features = "L1_control_inv",
    actions = "restoration",
    alias = "rodent_control"
  ) |>
  add_objective_max_effect(
    features = "L1_stock_inv",
    actions = "restoration",
    alias = "carbon_stock"
  ) |>
  add_objective_max_effect(
    features = "L1_rend_inv",
    actions = "restoration",
    alias = "water_yield"
  ) |>
  add_objective_max_effect(
    features = "L1_reten_inv",
    actions = "restoration",
    alias = "sediment_retention"
  )

print(problem)

Generate alternatives with AUGMECON

The purpose of the analysis is to reveal the consequences of alternative restoration priorities before decision makers specify a preferred trade-off. A weighted-sum model would require weights for cost and the four services at the outset. Those weights would embed an exchange rule before stakeholders had seen how it changes the spatial plan.

AUGMECON is used instead as an a posteriori method. Cost is treated as the primary objective, while the four ecosystem-service objectives become explicit epsilon constraints. Each run asks a different version of the same management question: what is the least-cost restoration portfolio that also satisfies a particular combination of minimum deficit benefits at the selected sites?

set_runs_grid(3) constructs three epsilon levels for each of the four secondary objectives from their payoff ranges. Their Cartesian product produces (3^4=81) requested configurations. These configurations sample contrasting combinations of minimum deficit benefits, from relatively permissive requirements to combinations that force restoration strongly towards low-service areas.

The grid is intentionally small for a vignette involving more than 30,000 planning units. It is sufficient to demonstrate the many-objective workflow, but it should not be interpreted as an exhaustive representation of the five-dimensional Pareto frontier. A decision-support application could refine the grid or introduce policy-defined thresholds with set_runs_manual().

Lexicographic anchoring is used to construct the payoff table, and the AUGMECON augmentation rewards unused slack after cost has been optimized. This helps avoid weakly efficient plans that improve one criterion only trivially while performing unnecessarily poorly on others.

Some combinations of stringent service limits may be infeasible. In this context, infeasibility has a substantive interpretation: the fixed restoration area cannot always be concentrated simultaneously in locations with very low levels of all four services. Such runs reveal conflicts among the spatial deficits rather than representing a failure of the analysis.

problem <- problem |>
  set_method_augmecon(
    primary = "cost",
    aliases = objective_aliases,
    runs = set_runs_grid(3),
    lexicographic = TRUE,
    augmentation = 1e-3,
    control = set_runs_control(
      stop_on_infeasible = FALSE,
      stop_on_no_solution = FALSE
    )
  ) |>
  set_solver_gurobi(gap_limit = 0)

The complete optimization requires Gurobi and a valid licence. It took about 18 minutes on the system used to prepare this vignette. Re-solving all 81 configurations whenever the vignette or package website is built would add substantial and unnecessary computational cost.

The vignette therefore loads a precomputed SolutionSet. Its plans were originally optimized by minimizing reference amounts. The cached object has been migrated to deficit benefits by evaluating the same selected actions, using the exact fixed-effort identity above. Objective values, optimization senses, epsilon bounds, feature summaries, and model vectors have been updated; the spatial allocations and original solver diagnostics are retained. No new optimization is implied by this migration.

The affine identity establishes equivalence for fixed restoration effort; it does not guarantee bit-for-bit reproduction after reoptimization. A fresh solver run can select different plans at numerical epsilon boundaries or among plans of equal cost. Current MO solves use precise solver defaults, scale the scalar criterion without changing its meaning, and check secondary bounds against reevaluated decisions. Inspect each stored solution's diagnostics$epsilon_checks and diagnostics$solver_args when regenerating the analysis. The historical cached plans retain their original diagnostics.

The bundled cache preserves the action assignments, objective values, feature summaries, solver diagnostics, and run metadata needed for the remaining analysis. It stores one shared planning-data snapshot and omits compiled model matrices and repeated per-run problem snapshots. These omissions reduce the package size without changing the plans, trade-offs, or maps. To optimize again, use the problem constructed above rather than the analysis cache.

solutions <- readRDS(file.path(
  "data",
  "integrated-ecosystem-services-solutions.rds"
))

The code below shows how to regenerate and save the complete solution set after changing the formulation, run design, or solver configuration. It is displayed for reproducibility but is not evaluated during ordinary vignette or pkgdown builds.

solutions <- solve(problem)

saveRDS(
  solutions,
  file.path("..", "data-raw", "restoration-full-solutions.rds"),
  compress = "xz"
)

The complete result is saved outside the bundled analysis cache. When preparing a package release from a repository checkout, use data-raw/compact-restoration-cache.R to create the smaller cache used above:

source(file.path("..", "data-raw", "compact-restoration-cache.R"))
saveRDS(
  compact_restoration_solution_set(solutions),
  file.path("data", "integrated-ecosystem-services-solutions.rds"),
  compress = "xz"
)

Examine the many-objective solution set

Run design and objective values

The run table records the 81 requested combinations of ecosystem-service thresholds. A run describes an attempted model configuration; it does not necessarily have an associated solution. Some configurations may be infeasible, while several feasible configurations may lead to the same spatial allocation because the decision variables are discrete.

For interpretation, the complete result is filtered to feasible, non-dominated solutions and repeated decision vectors are removed. The resulting SolutionSet contains the distinct efficient restoration portfolios generated by this experiment.

efficient_solutions <- solutions |>
  solution_filter(feasible_only = TRUE, nondominated = TRUE) |>
  solution_unique(by = "decisions")

run_summary <- data.frame(
  requested_configurations = nrow(get_runs(solutions)),
  retained_efficient_plans = nrow(get_runs(efficient_solutions))
)
run_summary

get_runs(efficient_solutions)
get_objectives(efficient_solutions, format = "wide")

The reduction from 81 requested configurations to a smaller set of unique efficient plans is expected. Nearby epsilon combinations can select the same 2,112 restoration units, dominated solutions add no new trade-off information, and some combinations of very restrictive service limits cannot be achieved simultaneously.

The retained plans should therefore be read as a set of substantively different restoration strategies, not as 81 equally distinct alternatives.

Read the solution set as a sequence of decisions

Five-dimensional trade-offs are difficult to understand from a complete matrix of pairwise plots alone. The following diagnostics organize the results around a simpler management question:

How much additional expenditure is required to move restoration away from sites with high current ecosystem-service provision and towards sites with lower baseline service levels?

Cost is minimized and the four deficit benefits are maximized. An improvement in a service objective means a larger deficit benefit than in the least-cost portfolio and, at fixed effort, a lower selected reference amount. In the logic of this example, the plan is targeting a stronger current deficit for that service. It does not mean that the service has already increased after restoration.

  objective_values <- get_objectives(
    efficient_solutions,
    format = "wide"
  )

  service_aliases <- objective_aliases[-1]
  objective_labels <- c(
    cost = "Cost",
    rodent_control = "Rodent control",
    carbon_stock = "Carbon stock",
    water_yield = "Water yield",
    sediment_retention = "Sediment retention"
  )

  minimum_cost_row <- which.min(objective_values$cost)
  minimum_cost_id <- objective_values$solution_id[minimum_cost_row]
  minimum_cost <- objective_values$cost[minimum_cost_row]

  compromise <- frontier_knee(
    efficient_solutions,
    objectives = objective_aliases,
    method = "ideal"
  )
  compromise_id <- compromise$solution_id

  service_best_ids <- setNames(
    vapply(
      service_aliases,
      function(service) {
        objective_values$solution_id[
          which.max(objective_values[[service]])
        ]
      },
      numeric(1)
    ),
    service_aliases
  )

  improvement_data <- do.call(
    rbind,
    lapply(service_aliases, function(service) {
      baseline <- objective_values[[service]][minimum_cost_row]
      ids <- objective_values$solution_id
      roles <- rep("Other efficient plan", length(ids))
      roles[ids == service_best_ids[[service]]] <- "Best for this service"
      roles[ids == compromise_id] <- "Five-objective compromise"
      roles[ids == minimum_cost_id] <- "Least cost"

      data.frame(
        solution_id = ids,
        service_alias = service,
        service = factor(
          objective_labels[[service]],
          levels = unname(objective_labels[service_aliases])
        ),
        additional_cost = 100 * (
          objective_values$cost / minimum_cost - 1
        ),
        benefit_increase = 100 * (
          objective_values[[service]] - baseline
        ) / baseline,
        role = roles,
        label = ifelse(
          roles == "Other efficient plan",
          "",
          paste("Solution", ids)
        )
      )
    })
  )

How much deficit targeting accompanies additional cost?

Each panel compares an efficient plan with the least-cost restoration plan, which is located at zero on both axes. Moving to the right means accepting a higher implementation cost. Moving upwards means increasing the deficit benefit of the corresponding service relative to the least-cost restoration plan.

A higher point therefore represents a portfolio that directs restoration more strongly towards locations where that service is currently low. Separate vertical scales are retained because the observed percentage range differs among services and because equal percentage changes need not have equal ecological or policy significance.

A point below zero would indicate that a more expensive plan selects sites with a smaller deficit benefit for that service than the least-cost plan.

  ggplot(
    improvement_data,
    aes(x = additional_cost, y = benefit_increase)
  ) +
    geom_hline(yintercept = 0, colour = "grey80", linewidth = 0.4) +
    geom_point(
      aes(colour = role),
      size = 2.6,
      alpha = 0.9
    ) +
    geom_text(
      aes(
        label = label,
        colour = role,
        hjust = dplyr::case_when(
          role == "Least cost" ~ 0,
          role == "Best for this service" ~ 1,
          TRUE ~ 0.5
        )
      ),
      vjust = -0.65,
      size = 3,
      show.legend = FALSE,
      check_overlap = FALSE
    ) +
    facet_wrap(~service, scales = "free_y", ncol = 2) +
    scale_colour_manual(
      values = c(
        "Other efficient plan" = "grey70",
        "Least cost" = "#D95F02",
        "Five-objective compromise" = "#7570B3",
        "Best for this service" = "#1B9E77"
      )
    ) +
    scale_x_continuous(
      labels = function(x) paste0(round(x), "%"),
      expand = expansion(mult = c(0.03, 0.16))
    ) +
    scale_y_continuous(
      labels = function(x) paste0(round(x, 1), "%"),
      expand = expansion(mult = c(0.05, 0.18))
    ) +
    labs(
      x = "Additional cost relative to the least-cost plan",
      y = "Increase in deficit benefit relative to the least-cost plan",
      colour = NULL
    ) +
    theme_minimal(base_size = 11) +
    theme(
      panel.grid.minor = element_blank(),
      strip.text = element_text(face = "bold"),
      legend.position = "bottom"
    )

The points are not connected because they represent alternative spatial portfolios, not successive stages of one restoration trajectory. Labels identify the least-cost plan, the five-objective compromise, and the best observed plan for the service shown in each panel. The same solution identifiers are used in the later heatmap and maps.

The figure should first be read within each panel. A steep increase indicates that a relatively small additional expenditure can redirect restoration towards a substantially stronger deficit for that service. A long movement to the right with little vertical improvement indicates declining returns. Comparing panels then reveals whether the same costly plan addresses several deficits or whether each service requires a different spatial allocation.

Because the services are not aggregated, the graph preserves the central planning conflict: restoration sites that are attractive for one ecological reason may be unattractive for another.

Compare the complete performance profile of each plan

The heatmap summarizes the full five-objective profile of every retained plan. Solutions are ordered from lowest to highest cost. Within each objective column, values are normalized across the observed efficient set. Blue indicates the best observed value for each objective (low cost or high benefit), red the worst, and intermediate colours show relative performance between those extremes.

This normalization supports visual comparison but does not make the objectives commensurable. It does not imply that cost and ecosystem services share units, or that the four services should receive equal policy weight.

  objective_matrix <- as.matrix(
    objective_values[, objective_aliases, drop = FALSE]
  )
  objective_ranges <- apply(objective_matrix, 2, range)
  normalized_matrix <- sweep(
    objective_matrix,
    2,
    objective_ranges[1, ],
    "-"
  )
  normalized_matrix <- sweep(
    normalized_matrix,
    2,
    objective_ranges[2, ] - objective_ranges[1, ],
    "/"
  )
  # Express every column as regret: low cost is good; high benefits are good.
  normalized_matrix[, service_aliases] <- 1 - normalized_matrix[, service_aliases]
  normalized_matrix[!is.finite(normalized_matrix)] <- 0

  heatmap_data <- do.call(
    rbind,
    lapply(seq_along(objective_aliases), function(j) {
      data.frame(
        solution_id = objective_values$solution_id,
        objective = factor(
          objective_labels[[objective_aliases[j]]],
          levels = unname(objective_labels[objective_aliases])
        ),
        relative_performance = normalized_matrix[, j]
      )
    })
  )
  heatmap_data$is_best <- heatmap_data$relative_performance < 1e-10

  solution_order <- objective_values$solution_id[
    order(objective_values$cost)
  ]
  solution_labels <- as.character(solution_order)
  solution_labels[solution_order == minimum_cost_id] <- paste0(
    minimum_cost_id,
    " - least cost"
  )
  solution_labels[solution_order == compromise_id] <- paste0(
    compromise_id,
    " - compromise"
  )
  heatmap_data$solution <- factor(
    heatmap_data$solution_id,
    levels = solution_order,
    labels = solution_labels
  )

  ggplot(
    heatmap_data,
    aes(x = objective, y = solution, fill = relative_performance)
  ) +
    geom_tile(colour = "white", linewidth = 0.5) +
    geom_point(
      data = subset(heatmap_data, is_best),
      shape = 21,
      fill = "white",
      colour = "black",
      size = 2.2,
      stroke = 0.6
    ) +
    scale_fill_gradient2(
      low = "#2166AC",
      mid = "#F7F7F7",
      high = "#B2182B",
      midpoint = 0.5,
      limits = c(0, 1),
      breaks = c(0, 0.5, 1),
      labels = c("Best", "Intermediate", "Worst")
    ) +
    labs(
      x = NULL,
      y = "Solution (increasing cost from bottom to top)",
      fill = "Observed\nperformance"
    ) +
    theme_minimal(base_size = 11) +
    theme(
      panel.grid = element_blank(),
      axis.text.x = element_text(angle = 25, hjust = 1),
      legend.position = "right"
    )

White circles identify the best observed plan for each objective. Rows with blue cells across several service columns represent restoration portfolios that simultaneously target multiple ecosystem-service deficits. Rows that combine blue and red cells expose specialization: the plan performs strongly for some services by accepting weaker performance for others.

The cost column completes the story. It shows whether a portfolio that targets several ecological deficits can be obtained with a moderate budget increase or whether it requires a substantially more expensive spatial reallocation.

Identify extremes and representative compromises

frontier_extremes() identifies the best and worst observed solution for each objective. These extremes are useful reference points because they show the maximum contrast generated by the experiment. However, the plan that best targets one ecosystem-service deficit may be costly or may select sites with high baseline values of another service.

  extremes <- frontier_extremes(
    efficient_solutions,
    objectives = objective_aliases,
    ties = "first"
  )
  extremes

With five objectives, a two-dimensional geometric knee is not appropriate. frontier_knee(method = "ideal") instead identifies the observed plan closest to the normalized five-objective ideal. The ideal combines the best observed value of every criterion, even though no single feasible plan may attain all of them simultaneously.

The resulting plan is a descriptive compromise within the generated solution set. It is not a universally optimal restoration strategy and does not replace stakeholder preferences.

For ecological interpretation, the analysis also calculates distance to an ideal defined only by the four ecosystem-service objectives. We refer to this as joint ecosystem-service regret. A value of zero would mean that one plan simultaneously attains the largest observed deficit benefit for all four services. Larger values indicate increasing distance from that joint deficit- targeting profile.

Plotting this distance against additional cost keeps the financial and ecological dimensions visible. Because each service is normalized to its observed range, the diagnostic gives the four services equal influence. This equal weighting is used only to describe the generated alternatives; it does not alter the AUGMECON formulation.

  service_distances <- frontier_distances(
    efficient_solutions,
    objectives = service_aliases,
    reference = "ideal"
  )
  service_distances$cost <- objective_values$cost[
    match(service_distances$solution_id, objective_values$solution_id)
  ]
  service_distances$additional_cost <- 100 * (
    service_distances$cost / minimum_cost - 1
  )

  ecological_id <- service_distances$solution_id[
    which.min(service_distances$distance_to_ideal)
  ]
  specialist_id <- unname(service_best_ids[["sediment_retention"]])

  compromise[, c(
    "solution_id",
    objective_aliases,
    "knee_score",
    "knee_rank"
  )]

  data.frame(
    five_objective_compromise = compromise_id,
    closest_to_ecological_ideal = ecological_id,
    same_solution = compromise_id == ecological_id
  )

In this experiment, the five-objective compromise and the plan closest to the four-service ideal are both solution r compromise_id. This agreement is an empirical feature of the generated set, not a property guaranteed by AUGMECON.

The figure also highlights the best observed plan for each individual service. Because solution r compromise_id is simultaneously the water-yield extreme, the highlighted roles correspond to the same five unique plans mapped below. The colour scheme is retained across figures: orange represents the least-cost plan, magenta the joint compromise, and green a service-specific extreme.

highlighted_ids <- unique(c(
  minimum_cost_id,
  compromise_id,
  unname(service_best_ids)
))

distance_plot_data <- service_distances
distance_plot_data$role <- "Other efficient plan"
distance_plot_data$role[
  distance_plot_data$solution_id %in% unname(service_best_ids)
] <- "Service-specific extreme"
distance_plot_data$role[
  distance_plot_data$solution_id == compromise_id
] <- "Five-objective compromise"
distance_plot_data$role[
  distance_plot_data$solution_id == minimum_cost_id
] <- "Least cost"

ggplot(
  distance_plot_data,
  aes(x = additional_cost, y = distance_to_ideal)
) +
  geom_point(
    aes(colour = role),
    size = 2.7,
    alpha = 0.9
  ) +
  ggrepel::geom_text_repel(
    data = subset(
      distance_plot_data,
      solution_id %in% highlighted_ids
    ),
    aes(
      label = paste("Solution", solution_id),
      colour = role
    ),
    seed = 123,
    box.padding = 0.35,
    point.padding = 0.25,
    min.segment.length = 0,
    max.overlaps = Inf,
    show.legend = FALSE,
    size = 3.1
  ) +
  scale_colour_manual(
    values = c(
      "Other efficient plan" = "grey70",
      "Least cost" = "#D95F02",
      "Five-objective compromise" = "#7570B3",
      "Service-specific extreme" = "#1B9E77"
    )
  ) +
  scale_x_continuous(
    labels = function(x) paste0(round(x), "%"),
    expand = expansion(mult = c(0.06, 0.18))
  ) +
  scale_y_continuous(
    expand = expansion(mult = c(0.08, 0.16))
  ) +
  labs(
    x = "Additional cost relative to the least-cost plan",
    y = "Joint ecosystem-service regret\n(0 = best observed joint service profile)",
    colour = NULL
  ) +
  theme_minimal(base_size = 11) +
  theme(
    panel.grid.minor = element_blank(),
    legend.position = "bottom"
  )

Connect objective performance to spatial decisions

Objective values describe how a plan performs, but they do not show which parts of the productive landscape would actually change. The maps therefore compare the least-cost plan, the best plan for each ecosystem-service objective, and the five-objective compromise.

Because the compromise is also the best observed plan for water yield, the set contains five unique spatial allocations rather than six. Together they represent three distinct planning logics: minimizing expenditure, targeting a particular ecosystem-service deficit, and balancing cost with all four deficits.

representative_plans <- data.frame(
  role = c(
    "Least-cost plan",
    "Five-objective compromise and water-yield best",
    "Rodent-control best",
    "Carbon-stock best",
    "Sediment-retention best"
  ),
  solution_id = c(
    minimum_cost_id,
    compromise_id,
    unname(service_best_ids[["rodent_control"]]),
    unname(service_best_ids[["carbon_stock"]]),
    unname(service_best_ids[["sediment_retention"]])
  )
)
representative_plans <- representative_plans[
  !duplicated(representative_plans$solution_id),
]

representative_summary <- representative_plans
representative_summary$additional_cost_percent <- round(
  100 * (
    objective_values$cost[
      match(
        representative_summary$solution_id,
        objective_values$solution_id
      )
    ] / minimum_cost - 1
  ),
  1
)

for (service in service_aliases) {
  baseline <- objective_values[[service]][minimum_cost_row]
  values <- objective_values[[service]][
    match(
      representative_summary$solution_id,
      objective_values$solution_id
    )
  ]
  representative_summary[[paste0(service, "_benefit_increase_percent")]] <- round(
    100 * (values - baseline) / baseline,
    1
  )
}

representative_summary

Positive percentages in the service columns represent increases in deficit benefits relative to the least-cost plan. For example, a positive carbon-stock value means that the plan restores sites with lower current carbon stock than those selected by the least-cost portfolio. It should not be interpreted as a measured increase in carbon after restoration.

The table provides the numerical context for the maps by showing the additional cost required to achieve each spatial redirection.

facet_labels <- setNames(
  paste0(
    "Solution ",
    representative_plans$solution_id,
    "\n",
    representative_plans$role
  ),
  representative_plans$solution_id
)

map_plot <- plot_spatial_actions(
  efficient_solutions,
  solutions = representative_plans$solution_id,
  fill_values = c(
    conservation = "#2E7D32",
    restoration = "#E69F00"
  ),
  base_alpha = 0.08
) +
  ggplot2::facet_wrap(
    ~solution_id,
    labeller = ggplot2::as_labeller(facet_labels)
  )

solution_colours <- setNames(
  rep("#1B9E77", nrow(representative_plans)),
  representative_plans$solution_id
)
solution_colours[as.character(minimum_cost_id)] <- "#D95F02"
solution_colours[as.character(compromise_id)] <- "#7570B3"

facet_order <- sort(representative_plans$solution_id)
facet_colours <- unname(
  solution_colours[as.character(facet_order)]
)

style_facet_strip <- function(grob, fill) {
  if (inherits(grob, "rect")) {
    grob$gp$fill <- fill
    grob$gp$col <- fill
  }
  if (inherits(grob, "text")) {
    grob$gp$col <- "white"
  }
  if (!is.null(grob$children)) {
    for (i in seq_along(grob$children)) {
      grob$children[[i]] <- style_facet_strip(
        grob$children[[i]],
        fill
      )
    }
  }
  if (!is.null(grob$grobs)) {
    for (i in seq_along(grob$grobs)) {
      grob$grobs[[i]] <- style_facet_strip(
        grob$grobs[[i]],
        fill
      )
    }
  }
  grob
}

map_grob <- ggplot2::ggplotGrob(map_plot)
strip_index <- grep("^strip-t", map_grob$layout$name)
strip_index <- strip_index[
  order(
    map_grob$layout$t[strip_index],
    map_grob$layout$l[strip_index]
  )
]

for (i in seq_along(strip_index)) {
  map_grob$grobs[[strip_index[i]]] <- style_facet_strip(
    map_grob$grobs[[strip_index[i]]],
    facet_colours[i]
  )
}

grid::grid.newpage()
grid::grid.draw(map_grob)

The maps are bundled as a pre-rendered image to keep package and website builds fast. The map-generation code is retained in the vignette source and is not evaluated during ordinary builds.

knitr::include_graphics(
  "figures/ecosystem-services/integrated-ecosystem-services-map-1.png"
)

The maps reveal the spatial mechanism behind the objective trade-offs. They show whether improved targeting of a service deficit is achieved through a small number of local substitutions or through a broader relocation of the restoration portfolio across the landscape.

Because the conservation network is fixed, all differences among panels arise from the 2,112 restoration decisions. Comparing the service-specific extremes shows whether different deficits point to common restoration areas or to contrasting parts of the productive matrix. The compromise map identifies the allocation that balances these competing spatial signals rather than optimizing one of them in isolation.

Compare decision-space similarity

Two plans can have similar objective values while selecting different places. Jaccard similarity complements the performance diagnostics by measuring overlap in selected planning-unit--action pairs.

The heatmap follows the same cost ordering used above. Higher similarity means that two plans prescribe more of the same spatial actions. Exact duplicate decision vectors have already been removed, so values close to one indicate very similar but not necessarily identical portfolios.

The fixed conservation network contributes the same assignments to every solution and can therefore raise overall similarity. Differences among plans are generated only by the variable restoration component. Consequently, the heatmap should be interpreted as similarity of the complete intervention system in a setting where one large component is intentionally held constant.

  similarity_matrix <- selection_similarity(
    efficient_solutions,
    metric = "jaccard",
    format = "matrix"
  )

  similarity_data <- expand.grid(
    solution_x = rownames(similarity_matrix),
    solution_y = colnames(similarity_matrix),
    stringsAsFactors = FALSE
  )
  similarity_data$similarity <- as.vector(similarity_matrix)
  similarity_levels <- as.character(solution_order)
  similarity_data$solution_x <- factor(
    similarity_data$solution_x,
    levels = similarity_levels
  )
  similarity_data$solution_y <- factor(
    similarity_data$solution_y,
    levels = rev(similarity_levels)
  )

  ggplot(
    similarity_data,
    aes(x = solution_x, y = solution_y, fill = similarity)
  ) +
    geom_tile(colour = "white", linewidth = 0.25) +
    scale_fill_gradientn(
      colours = c(
        "#440154",
        "#31688E",
        "#35B779",
        "#FDE725"
      ),
      limits = range(similarity_data$similarity, na.rm = TRUE),
      breaks = c(0.6, 0.8, 1.0)
    ) +
    coord_fixed() +
    labs(
      x = "Solution",
      y = "Solution",
      fill = "Jaccard\nsimilarity"
    ) +
    theme_minimal(base_size = 10) +
    theme(
      panel.grid = element_blank(),
      axis.text.x = element_text(angle = 45, hjust = 1)
    )

Selection frequency provides a more local view of stability. It reports how often each planning-unit--action assignment appears across the retained non-dominated solutions.

For restoration units, a high frequency indicates that a location is selected under many different combinations of cost and ecosystem-service limits. Such a site is comparatively robust to the trade-off configuration explored here. For conservation units, frequency one is expected because those assignments are locked into every solution.

Frequency is therefore a measure of stability within this generated solution set. It should not be interpreted as ecological irreplaceability, restoration success, or the probability that a site is objectively optimal.

  assignment_frequency <- selection_frequency(efficient_solutions)

  head(
    assignment_frequency[
      order(-assignment_frequency$frequency),
    ],
    15
  )

Interpretation and limitations

This vignette demonstrates how multiscape can prioritize additional forest restoration in a highly productive landscape where land is scarce, existing uses matter, and ecosystem-service deficits are spatially heterogeneous. The existing conservation network is treated as fixed, and every alternative adds the same 2,112 restoration units. The decision problem is therefore entirely about spatial allocation: which parts of the productive matrix should be restored?

The formulation gives a specific answer to that question. It favours locations with low implementation cost and low current levels of four ecosystem services. These low-service locations are treated as areas of ecological deficit, making them plausible candidates for restoration while avoiding sites where important functions are already strongly provided. By retaining the services as separate objectives, the model shows that there is no single universal map of ecological need: carbon, water, sediment retention, and biological control can point to different restoration priorities.

This interpretation must remain within the limits of the data and formulation. The deficit effects assume a common service ceiling. Low baseline service provision is not evidence of high realized restoration potential. A degraded site may respond strongly to intervention, weakly, or only after a long time. Conversely, a site with high current service provision could still benefit from restoration. The model identifies where current service deficits and implementation costs favour intervention; it does not predict the ecological response produced by that intervention.

A more complete application could combine the present deficit-based objectives with empirically estimated action outcomes, feasibility, time to recovery, land-use opportunity cost, connectivity, or future climate risk. Those additions would distinguish more explicitly between where services are currently low and where restoration is expected to improve them most.

Finally, the generated alternatives depend on the epsilon design. The 81 requested boundary configurations provide a transparent first exploration of the five-dimensional trade-off space, not an exhaustive approximation of the Pareto frontier. A real decision process should refine the thresholds around policy-relevant regions, test the stability of the spatial recommendations, and use stakeholder input to identify acceptable compromises after the consequences of the alternatives have been made visible.

multiscape documentation built on Oct. 12, 2026, 5:08 p.m.