---
title: "Multi-objective forest restoration planning in a highly productive landscape"
output:
  rmarkdown::html_vignette:
    toc: false
    number_sections: false
vignette: >
  %\VignetteIndexEntry{Multi-objective forest restoration planning in a highly productive landscape}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include=FALSE}
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.

```{r load-data}
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.

```{r plot-services, fig.width=8, fig.height=6}
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.

```{r action-eligibility}
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.

```{r create-problem, message=TRUE}
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.

```{r constraints}
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.

```{r deficit-effects}
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.

```{r objectives}
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.

```{r method}
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.

```{r load-solutions}
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.

```{r regenerate-solutions, eval=FALSE}
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:

```{r compact-regenerated-solutions, eval=FALSE}
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.

```{r inspect-results}
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.

```{r prepare-objective-comparison, include=FALSE}
  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.

```{r cost-service-improvement, echo=FALSE, fig.width=9, fig.height=6.5}
  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.

```{r objective-heatmap, echo=FALSE, fig.width=9, fig.height=7}
  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.

```{r objective-extremes}
  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.

```{r compromise-solution}
  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.

```{r cost-ecological-distance, echo=FALSE, fig.width=7.5, fig.height=5.5}
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.

```{r representative-plans, echo=FALSE}
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.

```{r map-selected-solutions, echo=FALSE, eval=FALSE, fig.keep='last', fig.width=11, fig.height=8}
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.

```{r map-selected-solutions-static, echo=FALSE, fig.width=11, fig.height=8}
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.

```{r selection-similarity, echo=FALSE, fig.width=8, fig.height=7}
  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.

```{r selection-frequency}
  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.