The hardware and bandwidth for this mirror is donated by dogado GmbH, the Webhosting and Full Service-Cloud Provider. Check out our Wordpress Tutorial.
If you wish to report a bug, or if you are interested in having us mirror your free-software or open-source project, please feel free to contact us at mirror[@]dogado.de.

Multi-action, multi-objective spatial planning in R

CRAN status CRAN downloads lifecycle R-CMD-check Codecov test coverage

multiscape is an exact optimisation framework for multi-action, multi-objective spatial planning in R. It addresses planning problems in which management actions are explicit decision alternatives whose spatial allocation must satisfy planning requirements while accounting for multiple competing objectives. Building on the transition from place-based prioritisation towards spatially explicit action planning (Tallis et al., 2021; Salgado-Rojas et al., 2023), multiscape jointly represents where management actions can be implemented, which actions are selected, and the consequences expected from those decisions.

Planning units, features, feasible actions, action-specific costs and outcomes, and spatial or management constraints are represented within a mixed-integer linear programming (MILP) formulation. Multiple objectives are defined independently, allowing ecological, economic, and spatial criteria to remain explicit rather than being collapsed a priori into a single criterion. Alternative spatial plans can then be generated and compared to characterise trade-offs among objectives and the corresponding spatial allocation of management actions. This formulation supports conservation and restoration applications as well as broader multiple-use spatial planning problems involving ecological, economic, and social objectives (Neubert et al., 2025).

Installation

Install the released version of multiscape from Comprehensive R Archive Network (CRAN):

install.packages("multiscape")

Alternatively, install the latest development version from GitHub:

if (!requireNamespace("remotes", quietly = TRUE)) {
  install.packages("remotes")
}
remotes::install_github("josesalgr/multiscape")

Getting started

Multi-action planning in the Meseta Ibérica

We illustrate the multiscape workflow using a spatial planning problem from Spain’s Meseta Ibérica, adapted from Cánibe Iglesias et al. (2025). The example comprises 11,109 planning units, 151 species, four ecosystem-service features, and four alternative management actions: Afforestation, BAU, FarmReturn, and Firesmart.

Understand the planning inputs

Before constructing the optimisation problem, we examine its spatial units, features, management alternatives, implementation costs, and expected outcomes. The example data are included in multiscape and can be loaded directly.

library(multiscape)
meseta <- load_meseta()

Planning units

The study area is divided into 11,109 planning units of approximately 1 km². Each unit represents a location where a management action may be implemented. The supplied spatial polygons and identifiers are stored in planning_units.

head(meseta$planning_units)
#> Simple feature collection with 6 features and 2 fields
#> Geometry type: POLYGON
#> Dimension:     XY
#> Bounding box:  xmin: 2936400 ymin: 2279200 xmax: 2945400 ymax: 2281200
#> Projected CRS: GRS_1980_IUGG_1980_Lambert_Azimuthal_Equal_Area
#>   id                       geometry cost
#> 1  1 POLYGON ((2944400 2280200, ...  0.1
#> 2  2 POLYGON ((2945400 2280200, ...  0.1
#> 3  3 POLYGON ((2937400 2279200, ...  0.1
#> 4  4 POLYGON ((2938400 2279200, ...  0.1
#> 5  5 POLYGON ((2939400 2279200, ...  0.1
#> 6  6 POLYGON ((2944400 2279200, ...  0.1

Features and reference amounts

The planning problem considers 155 features: 151 species and four ecosystem services. A feature is a measurable quantity represented in the planning model. The features table identifies these quantities, while dist_features provides their reference amounts by planning unit and feature. These amounts serve as the comparison values for action-specific outcomes; they do not define management decisions or representation targets.

head(meseta$features)
#>   id    name
#> 1  1 ACCGENT
#> 2  2 ACCNISU
#> 3  3 ACRARUN
#> 4  4 ACRSCIR
#> 5  5 AEGCAUD
#> 6  6     AER

In this example, dist_features assigns zero credited representation to the reference. Representation is instead credited through the outcomes associated with selected actions.

Management actions and implementation costs

Four alternative actions can be implemented throughout the study area: Afforestation, BAU, FarmReturn, and Firesmart. Their implementation costs vary by planning unit, so a plan’s cost depends on which actions are chosen and where they are implemented.

meseta$actions
#>              id          name
#> 1 Afforestation Afforestation
#> 2           BAU           BAU
#> 3    FarmReturn    FarmReturn
#> 4     Firesmart     Firesmart
plot(
  meseta$maps$costs,
  main = "Action-specific implementation costs",
  border = NA,
  pal = function(n) hcl.colors(n, "YlGnBu", rev = TRUE)
)

Action-specific outcomes

The same feature may respond differently to different management actions and locations. The outcomes table supplies its expected representation for each recorded planning-unit, feature, and action combination. In this dataset, outcomes are binary: a value of 1 credits one unit of representation, while unrecorded combinations contribute zero. This coding also applies to the included ecosystem-service features.

head(meseta$outcomes)
#>   pu feature outcome        action
#> 1  1       7       1 Afforestation
#> 2  1       8       1 Afforestation
#> 3  1      12       1 Afforestation
#> 4  1      13       1 Afforestation
#> 5  1      14       1 Afforestation
#> 6  1      15       1 Afforestation

The maps below illustrate the action-specific representation of DENMINO. They show outcomes, rather than differences from the reference.

plot(
  meseta$maps$DENMINO,
  main = "Action-specific outcomes for DENMINO",
  border = NA,
  breaks = c(-0.1, 0.5, 1.1),
  pal = function(n) c("grey95", "#238b45")
)

Build the optimisation problem

We now use these inputs to build a Problem object. The stages below first establish its spatial and feature structure, then add management decisions and their consequences, feasibility requirements, and objectives. No optimisation is performed until the method and solver have been configured.

Stage 1: Create the base problem

create_problem() defines the planning units, the feature catalogue, and the reference amounts stored in dist_features. We pass the supplied tables directly. In Stage 4, include_pu_cost = FALSE excludes the separate planning-unit cost from the objective, so this example counts only the costs of implementing specific actions, which we introduce in Stage 2.

problem <- create_problem(
  pu = meseta$planning_units,
  features = meseta$features,
  dist_features = meseta$dist_features,
  cost = "cost"
)

The initial Problem contains the spatial structure and 155 features, but no registered actions, effects, constraints, or objectives. This incomplete state is expected: those components are added in subsequent stages. We can inspect the model as it is being constructed.

print(problem)
#> A multiscape object (<Problem>)
#> ├─data
#> │├─planning units: <tbl_df> (11109 total)
#> │├─costs: min: 0.1, max: 0.1
#> │└─features: 155 total ("ACCGENT", "ACCNISU", "ACRARUN", ...)
#> └─actions and effects
#> │├─actions: none
#> │├─feasible action pairs: none
#> │├─effect data: none
#> │└─profit data: none
#> └─spatial
#> │├─geometry: sf (11109 rows)
#> │├─coordinates: 11109 rows (x: 2868900..3007900, y: 2110700..2280700)
#> │└─relations: none
#> └─targets and constraints
#> │├─targets: none
#> │├─area constraints: none
#> │├─budget constraints: none
#> │├─planning-unit locks: none
#> │└─action locks: none
#> └─model
#> │├─status: not built yet (will build in solve())
#> │├─objectives: none
#> │├─method: single-objective
#> │├─solver: not set (auto)
#> │└─checks: incomplete (no objective registered)
#> # ℹ Use `x$data` to inspect stored tables and model snapshots.

Stage 2: Add actions and their effects

add_actions() registers the available management decisions and their implementation costs; it does not choose where to implement them. The loader already uses readable action names as identifiers consistently across actions, action_costs, and outcomes, so no preparation is needed.

problem <- problem |> 
  add_actions(
    actions = meseta$actions,
    cost = meseta$action_costs
  )

An outcome describes the amount of a feature under a particular action; an effect describes its change relative to the reference. If \(r\) denotes the reference amount and \(y\) the action-specific outcome, the signed effect is \(\Delta = y-r\). add_effects() accepts three equivalent ways of specifying a given action’s consequences, using one of the following columns:

Column Supplied value Outcome \(y\) Effect \(\Delta\)
outcome Amount under the action outcome \(y-r\)
effect Signed absolute change \(r+\text{effect}\) effect
relative_change Proportional change \(q\) \(r(1+q)\) \(rq\)

For example, with a reference of 100, an outcome of 130, an effect of 30, and a relative change of 0.30 all describe the same resulting amount. They are alternative inputs, not values to supply together. When the reference is zero, proportional changes cannot produce a positive outcome; this example therefore uses the supplied outcome values directly.

problem <- problem |> 
  add_effects(
    effects = meseta$outcomes
  )

At this point, the problem describes both the available spatial decisions and the consequences expected if they are implemented. The next stage determines which combinations of decisions are admissible.

Stage 3: Define the constraints

Constraints define the feasible set of spatial plans. We require at most one action per planning unit and impose an absolute representation target for each feature. sense = "max" makes the action-cardinality constraint an upper bound, so units may remain without a selected action. The feature targets require the sum of credited representation to reach the supplied minimum for each feature.

head(meseta$targets)
#>   feature target
#> 1       1   92.8
#> 2       2  109.6
#> 3       3   26.4
#> 4       4  180.0
#> 5       5  102.4
#> 6       6  269.6

problem <- problem |>
  add_constraint_action_cardinality(count = 1, sense = "max") |>
  add_constraint_targets_absolute(meseta$targets)

These targets are expressed in the outcome units rather than as percentages of the reference. Given the binary coding of this example, a fractional target must still be met or exceeded by the resulting total representation.

Stage 4: Define the objectives

Constraints determine whether a plan is feasible; objectives determine how feasible plans are evaluated. Here we register two quantities to minimise: the cost of implementing the chosen actions and their spatial fragmentation. They remain separate objectives until an optimisation method specifies how to combine them.

problem <- problem |>
  add_spatial_relations(meseta$boundary, name = "boundary") |>
  add_objective_min_cost(include_pu_cost = FALSE, alias = "cost") |>
  add_objective_min_fragmentation_action(
    relation_name = "boundary", alias = "spatial"
  )

The fragmentation objective sums weighted boundaries where the selection of a given action differs between neighbouring units. Matching action assignments contribute no boundary for that action; different actions contribute to their respective boundaries. Minimising this criterion favours more cohesive action patches without imposing strict connectivity. The supplied boundary table provides the spatial relations and their weights.

Stage 5: Configure and solve

The optimisation method specifies how registered objectives are handled. For a first illustrative run, we use a weighted sum with coefficients 1 for cost and 0.5 for fragmentation, minimising \(\text{cost}+0.5\,\text{spatial}\). Because weight normalisation and objective scaling are disabled, this choice is sensitive to the input units and should be understood as an illustrative preference, not a universally meaningful trade-off rate.

We save the formulation before configuring the method so that it can also be used for the multi-run exploration below.

problem <- problem |> 
  set_method_weighted_sum(
  aliases = c("cost", "spatial"),
  runs = set_runs_manual(data.frame(
    weight_cost = 1,
    weight_spatial = 0.5
  )),
  normalize_weights = FALSE,
  objective_scaling = FALSE
)

The solver then searches for a feasible assignment that minimises the weighted criterion. We use Gurobi, which requires a valid installation and licence, with a 300-second solver limit and a 5% relative optimality-gap limit. Building the optimisation model can require additional time and memory. The optimisation is executed when knitting this README and its result is cached, so subsequent renders need not repeat the solve unless the code changes. Rendering this chunk requires an installed and licensed Gurobi solver.

problem <- set_solver_gurobi(
  problem,
  gap_limit = 0.05,
  time_limit = 300,
  verbose = TRUE
)

solutions <- solve(problem)
#> Set parameter Username
#> Set parameter LicenseID to value 2844238
#> Set parameter TimeLimit to value 300
#> Set parameter FeasibilityTol to value 1e-09
#> Set parameter IntFeasTol to value 1e-09
#> Set parameter MIPGap to value 0.05
#> Set parameter MIPGapAbs to value 0
#> Set parameter OptimalityTol to value 1e-09
#> Set parameter NodefileStart to value 0.5
#> Academic license - for non-commercial use only - expires 2027-07-14
#> Gurobi Optimizer version 12.0.2 build v12.0.2rc0 (win64 - Windows 10.0 (19045.2))
#> 
#> CPU model: 12th Gen Intel(R) Core(TM) i7-12700H, instruction set [SSE2|AVX|AVX2]
#> Thread count: 14 physical cores, 20 logical processors, using up to 20 threads
#> 
#> Non-default parameters:
#> TimeLimit  300
#> FeasibilityTol  1e-09
#> IntFeasTol  1e-09
#> MIPGap  0.05
#> MIPGapAbs  0
#> OptimalityTol  1e-09
#> NodefileStart  0.5
#> LogToConsole  0
#> 
#> Optimize a model with 327122 rows, 142317 columns and 2120102 nonzeros
#> Model fingerprint: 0x18f5272c
#> Variable types: 86772 continuous, 55545 integer (55545 binary)
#> Coefficient statistics:
#>   Matrix range     [1e+00, 1e+00]
#>   Objective range  [1e+00, 7e+00]
#>   Bounds range     [1e+00, 1e+00]
#>   RHS range        [8e-01, 9e+03]
#> Presolve removed 142334 rows and 11125 columns
#> Presolve time: 4.78s
#> Presolved: 184788 rows, 131192 columns, 1681755 nonzeros
#> Variable types: 0 continuous, 131192 integer (131192 binary)
#> Deterministic concurrent LP optimizer: primal simplex, dual simplex, and barrier
#> Showing barrier log only...
#> 
#> Root barrier log...
#> 
#> Ordering time: 1.96s
#> 
#> Barrier statistics:
#>  AA' NZ     : 5.943e+06
#>  Factor NZ  : 2.145e+07 (roughly 300 MB of memory)
#>  Factor Ops : 4.108e+09 (less than 1 second per iteration)
#>  Threads    : 11
#> 
#>                   Objective                Residual
#> Iter       Primal          Dual         Primal    Dual     Compl     Time
#>    0   9.08549505e+07 -4.09748130e+05  2.40e+03 0.00e+00  2.36e+03     9s
#>    1   3.47762789e+07 -3.33400318e+05  8.66e+02 1.58e+00  7.33e+02    10s
#>    2   2.81132534e+06 -3.21967467e+05  8.05e+01 8.14e-01  6.08e+01    10s
#>    3   9.02184971e+05 -3.10629553e+05  2.45e+01 4.03e-01  1.91e+01    10s
#>    4   3.87322807e+05 -3.01982316e+05  9.29e+00 2.60e-01  7.71e+00    11s
#>    5   2.24942854e+05 -2.87930717e+05  4.49e+00 1.95e-01  4.08e+00    11s
#>    6   1.21722132e+05 -1.57332447e+05  1.92e+00 7.83e-05  1.43e+00    11s
#>    7   6.70445867e+04 -6.98603867e+04  5.85e-01 5.77e-15  4.66e-01    11s
#>    8   5.56541059e+04 -4.10644675e+04  3.67e-01 6.00e-15  2.90e-01    12s
#>    9   4.25341831e+04 -1.91292600e+04  1.67e-01 4.00e-15  1.60e-01    12s
#>   10   4.16902051e+04 -4.30829840e+03  1.59e-01 6.88e-15  1.17e-01    12s
#>   11   3.81011615e+04  9.81693227e+03  1.10e-01 3.66e-15  6.78e-02    13s
#>   12   3.33141639e+04  1.82504425e+04  5.66e-02 5.66e-15  3.44e-02    13s
#>   13   2.98016261e+04  2.19762987e+04  2.47e-02 7.69e-15  1.75e-02    13s
#>   14   2.88989373e+04  2.25255123e+04  1.84e-02 5.19e-15  1.42e-02    14s
#>   15   2.87336352e+04  2.36594964e+04  1.73e-02 4.00e-15  1.12e-02    14s
#>   16   2.76331951e+04  2.42516480e+04  1.04e-02 4.88e-15  7.48e-03    14s
#>   17   2.70453799e+04  2.45937589e+04  7.05e-03 5.50e-15  5.42e-03    14s
#>   18   2.65020957e+04  2.48596222e+04  4.22e-03 6.95e-15  3.63e-03    15s
#>   19   2.60274969e+04  2.52286068e+04  2.04e-03 4.52e-15  1.76e-03    15s
#>   20   2.58833211e+04  2.53155146e+04  1.48e-03 2.53e-15  1.25e-03    15s
#>   21   2.58179479e+04  2.53643921e+04  1.23e-03 4.88e-15  9.99e-04    16s
#>   22   2.57076401e+04  2.53843108e+04  8.25e-04 5.25e-15  7.13e-04    16s
#>   23   2.56171574e+04  2.54040800e+04  5.03e-04 8.76e-15  4.70e-04    16s
#>   24   2.55623126e+04  2.54257984e+04  3.11e-04 3.61e-15  3.01e-04    16s
#>   25   2.55186149e+04  2.54374382e+04  1.62e-04 3.58e-15  1.80e-04    17s
#>   26   2.54990910e+04  2.54513899e+04  9.99e-05 5.54e-15  1.05e-04    17s
#>   27   2.54875587e+04  2.54560979e+04  6.55e-05 4.63e-15  6.95e-05    17s
#>   28   2.54839890e+04  2.54580283e+04  5.50e-05 2.09e-15  5.74e-05    17s
#>   29   2.54810141e+04  2.54597586e+04  4.59e-05 5.38e-15  4.70e-05    18s
#>   30   2.54729208e+04  2.54612117e+04  2.15e-05 2.33e-15  2.59e-05    18s
#>   31   2.54703240e+04  2.54632141e+04  1.42e-05 3.87e-15  1.57e-05    18s
#>   32   2.54691539e+04  2.54640654e+04  1.09e-05 3.95e-15  1.12e-05    18s
#>   33   2.54659252e+04  2.54641982e+04  1.95e-06 6.77e-15  3.84e-06    19s
#>   34   2.54651950e+04  2.54648110e+04  1.66e-07 5.91e-15  8.56e-07    19s
#>   35   2.54651237e+04  2.54650654e+04  5.57e-08 4.41e-15  1.30e-07    19s
#>   36   2.54650851e+04  2.54650770e+04  3.69e-09 5.87e-15  1.79e-08    20s
#>   37   2.54650823e+04  2.54650821e+04  1.51e-10 1.23e-12  4.88e-10    20s
#> 
#> Barrier solved model in 37 iterations and 19.82 seconds (18.67 work units)
#> Optimal objective 2.54650823e+04
#> 
#> 
#> Root crossover log...
#> 
#>   165733 DPushes remaining with DInf 0.0000000e+00                20s
#>        0 DPushes remaining with DInf 0.0000000e+00                23s
#> 
#>    22178 PPushes remaining with PInf 8.2450785e-06                23s
#>        0 PPushes remaining with PInf 0.0000000e+00                24s
#> 
#>   Push phase complete: Pinf 0.0000000e+00, Dinf 2.7618265e+03     24s
#> 
#> 
#> Root simplex log...
#> 
#> Iteration    Objective       Primal Inf.    Dual Inf.      Time
#>   174690    2.5465082e+04   0.000000e+00   2.761826e+03     24s
#>   175374    2.5465082e+04   0.000000e+00   0.000000e+00     25s
#> Concurrent spin time: 0.00s
#> 
#> Solved with barrier
#> 
#> Root relaxation: objective 2.546508e+04, 175374 iterations, 19.20 seconds (13.70 work units)
#> Total elapsed time = 25.66s (DegenMoves)
#> 
#>     Nodes    |    Current Node    |     Objective Bounds      |     Work
#>  Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time
#> 
#>      0     0 25465.0822    0 2133          - 25465.0822      -     -   26s
#> H    0     0                    25633.000000 25465.0822  0.66%     -   27s
#> 
#> Explored 1 nodes (182306 simplex iterations) in 27.04 seconds (26.90 work units)
#> Thread count was 20 (of 20 available processors)
#> 
#> Solution count 1: 25633 
#> 
#> Optimal solution found (tolerance 5.00e-02)
#> Best objective 2.563300000000e+04, best bound 2.546550000000e+04, gap 0.6535%

We inspect the returned SolutionSet before analysing its objective values. The print method summarises the method, attempted runs, stored solutions, and solver diagnostics.

print(solutions)
#> A multiscape solution set (<SolutionSet>)
#> ├─method
#> │├─name: `weighted`
#> │├─objectives: 2 (cost, spatial)
#> │└─run design: manual
#> │ └─design detail: 1 requested runs
#> └─content
#> │├─design rows: 1
#> │├─attempted runs: 1
#> │├─stored solutions: 1
#> │└─without solution: 0
#> └─run summary
#> │├─statuses: optimal: 1
#> │├─runtime: 27.217
#> │├─gap: 0.0065
#> │├─design columns: none
#> │└─objective columns: value_cost, value_spatial
#> └─objective ranges
#> │├─cost: 23980
#> │└─spatial: 3306
#> # ℹ Use get_runs(), get_objectives(), get_pu(), and get_actions() to inspect
#> results.

A 5% gap is a stopping tolerance, not a guarantee that the exact optimum or a unique spatial plan has been found. Check the reported solver status and gap before interpreting any returned solutions.

Interpret the solutions

The following analyses use the solutions object from Stage 5. The first results concern a single weighted run; we then vary the weight on spatial fragmentation to examine alternative plans.

get_runs() links the requested run to the retained solution and reports its solver status, runtime, and gap. get_objectives() provides the component objective values in their original units. The weighted value of this run can be reconstructed using the same coefficients as above.

get_runs(solutions)
#>   run_id solution_id  status runtime         gap
#> 1      1           1 optimal  27.217 0.006534545

performance <- get_objectives(solutions, format = "wide")

performance
#>   solution_id  cost spatial
#> 1           1 23980    3306

Solver versions and stopping conditions can lead to different feasible plans. Interpret the values returned by your run rather than expecting a particular spatial solution.

Check target achievement

Targets express the required feature representation and should be checked independently of cost and fragmentation. get_targets() reports the required and achieved amounts for each feature, allowing all targets to be checked.

target_achievement <- get_targets(solutions)

head(target_achievement)
#>   solution_id feature feature_name target_level total_available target achieved
#> 1           1       1      ACCGENT         92.8              NA   92.8     2384
#> 2           1       2      ACCNISU        109.6              NA  109.6     2636
#> 3           1       3      ACRARUN         26.4              NA   26.4      572
#> 4           1       4      ACRSCIR        180.0              NA  180.0     1068
#> 5           1       5      AEGCAUD        102.4              NA  102.4      103
#> 6           1       6          AER        269.6              NA  269.6      270
#>      gap  met
#> 1 2291.2 TRUE
#> 2 2526.4 TRUE
#> 3  545.6 TRUE
#> 4  888.0 TRUE
#> 5    0.6 TRUE
#> 6    0.4 TRUE

all(target_achievement$met)
#> [1] TRUE

A value of TRUE in the final check indicates that all reported targets were met. Confirm that the solver returned a feasible solution before interpreting target achievement.

Decision space: inspect the selected actions

Objective values describe a plan’s performance; its decision space representation records which action was assigned to each planning unit. get_actions() returns the action-level results, including a selected indicator. Rather than plotting all 11,109 units, we display the first ten selected assignments and the number of units assigned to each action.

selected_actions <- get_actions(solutions, solution = 1)

selected_actions <- selected_actions[
  selected_actions$selected == 1,
  c("pu", "action", "cost"),
  drop = FALSE
]

head(selected_actions, 10)
#>    pu action cost
#> 2   1    BAU    1
#> 6   2    BAU    1
#> 10  3    BAU    1
#> 14  4    BAU    1
#> 18  5    BAU    1
#> 22  6    BAU    1
#> 26  7    BAU    1
#> 30  8    BAU    1
#> 34  9    BAU    1
#> 38 10    BAU    1

table(selected_actions$action)
#> 
#> Afforestation           BAU    FarmReturn     Firesmart 
#>          2363          2719           423          3185

The selected actions can also be inspected spatially. plot_spatial_actions() shows the complete spatial prescription for the first solution. Each colour identifies the selected management action, while grey planning units have no selected action. At most one action can be selected per planning unit.

plot_spatial_actions(
    solutions,
    solutions = 1,
    layout = "single",
    fill_values = c(
      Afforestation = "#E6C34A",
      BAU = "#41AB5D",
      FarmReturn = "#2B8CBE",
      Firesmart = "#9B7653"
  )
)

Objective space: explore the cost-fragmentation trade-off

To examine the compromise between cost and spatial cohesion, we keep the planning data and constraints fixed and vary the weight on fragmentation. We rebuild the same formulation in problemMO to keep this exploration independent of the preceding single-run problem. Each row in set_runs_manual() represents one run.

problemMO <- create_problem(
    pu = meseta$planning_units,
    features = meseta$features,
    dist_features = meseta$dist_features,
    cost = "cost"
  ) |> 
    add_actions(
    actions = meseta$actions,
    cost = meseta$action_costs
  ) |> 
    add_effects(
    effects = meseta$outcomes
  ) |> 
  add_constraint_action_cardinality(count = 1, sense = "max") |>
  add_constraint_targets_absolute(meseta$targets) |> 
  add_spatial_relations(meseta$boundary, name = "boundary") |>
  add_objective_min_cost(include_pu_cost = FALSE, alias = "cost") |>
  add_objective_min_fragmentation_action(
    relation_name = "boundary", alias = "spatial"
  ) |> 
  set_solver_gurobi(
    gap_limit = 0.05,
    time_limit = 300,
    verbose = FALSE
  )

spatial_weights <- c(0, 0.1, 0.25, 0.5, 1, 2)

spatial_exploration <- problemMO |> 
  set_method_weighted_sum(
  aliases = c("cost", "spatial"),
  runs = set_runs_manual(data.frame(
    weight_cost = 1,
    weight_spatial = spatial_weights
  )),
  normalize_weights = FALSE,
  objective_scaling = FALSE
)

alternatives <- solve(spatial_exploration)

The objective values of the six solutions can be inspected directly with get_objectives() and visualised using plot_tradeoff().

get_objectives(alternatives, format = "wide")
#>   solution_id  cost spatial
#> 1           1 22637   22485
#> 2           2 22765    7865
#> 3           3 23173    5231
#> 4           4 23980    3306
#> 5           5 24944    1822
#> 6           6 25996    1090

plot_tradeoff(
    alternatives,
    label_runs = TRUE,
    objectives = c("cost", "spatial"),
    connect = FALSE
)

Both objectives are minimised, so smaller values are preferred on either axis. A larger spatial coefficient gives the solver a stronger incentive to accept higher costs if this reduces fragmentation. These six weighted runs sample the trade-off: they do not establish its complete Pareto frontier.

The returned objective values show the observed compromises directly: a lower fragmentation score may require a higher implementation cost, while a cost-focused run may produce a more fragmented assignment. Because the solver stops at a finite optimality gap, the sample can include dominated solutions. Subsequent frontier diagnostics describe the observed set, not a complete Pareto frontier.

Decision space: map alternative action allocations

The trade-off plot summarises how the objectives change, but not which management actions are selected in each location. We therefore map the spatial assignments in the six retained alternatives, using consistent colours across panels. Grey planning units have no selected action.

plot_spatial_actions(
    alternatives,
    solutions = get_objectives(alternatives, format = "wide")$solution_id,
    layout = "single",
    ncol = 2,
    fill_values = c(
      Afforestation = "#E6C34A",
      BAU = "#41AB5D",
      FarmReturn = "#2B8CBE",
      Firesmart = "#9B7653"
    )
)

These panels illustrate how increasing the fragmentation weight changes the location and composition of selected actions. They show solutions for the bundled inputs, rather than the ten-replicate scenario analysis in the original study.

Compare spatial similarity and recurring assignments

Objective values reveal performance, but do not show how similar two spatial prescriptions are. In decision space, Jaccard similarity measures overlap in selected planning unit-action pairs. A value of 1 means identical assignments; 0 means no assignments in common.

selection_similarity(alternatives, metric = "jaccard", format = "matrix")
#>           1         2         3         4         5         6
#> 1 1.0000000 0.3256152 0.2932796 0.2495680 0.2002491 0.1745742
#> 2 0.3256152 1.0000000 0.6390935 0.4364974 0.2988033 0.2468110
#> 3 0.2932796 0.6390935 1.0000000 0.5766297 0.4054251 0.3430646
#> 4 0.2495680 0.4364974 0.5766297 1.0000000 0.6582665 0.5065351
#> 5 0.2002491 0.2988033 0.4054251 0.6582665 1.0000000 0.7065921
#> 6 0.1745742 0.2468110 0.3430646 0.5065351 0.7065921 1.0000000
#> attr(,"metric")
#> [1] "jaccard"

action_frequency <- selection_frequency(alternatives)
head(action_frequency[order(-action_frequency$frequency), ], 10)
#>        pu        action n_selected n_solutions frequency
#> 80  10014     Firesmart          6           6         1
#> 82  10015           BAU          6           6         1
#> 145  1003 Afforestation          6           6         1
#> 185 10039 Afforestation          6           6         1
#> 193 10040 Afforestation          6           6         1
#> 197 10041 Afforestation          6           6         1
#> 204 10042     Firesmart          6           6         1
#> 214 10045           BAU          6           6         1
#> 317 10069 Afforestation          6           6         1
#> 325 10070 Afforestation          6           6         1

Assignment frequency measures recurrence within the six retained plans, not an ecological selection probability or irreplaceability measure.

Identify an empirical compromise

An observed solution is dominated when another observed solution is no worse on either objective and better on at least one. frontier_knee() identifies a potential compromise among the observed non-dominated solutions. Its result depends on the sampled objective values and their normalisation, so it should be treated as a decision aid rather than a uniquely best plan. frontier_extremes() and frontier_distances() provide additional diagnostics relative to the observed extreme, ideal, and nadir values.

# Which observed solution offers an empirical compromise?
knee <- frontier_knee(alternatives, objectives = c("cost", "spatial"))
knee[, c("solution_id", "cost", "spatial", "knee_score"), drop = FALSE]
#>   solution_id  cost spatial knee_score
#> 1           3 23173    5231  0.4574124

# Best and worst observed values for each objective.
frontier_extremes(alternatives, objectives = c("cost", "spatial"))
#>   solution_id objective sense bound  role value
#> 1           1      cost   min   min  best 22637
#> 2           6      cost   min   max worst 25996
#> 3           6   spatial   min   min  best  1090
#> 4           1   spatial   min   max worst 22485

# Distance and ranking relative to the empirical ideal and nadir.
distances <- frontier_distances(
  alternatives,
  objectives = c("cost", "spatial"),
  reference = c("ideal", "nadir")
)
distances[, c(
  "solution_id", "distance_to_ideal", "rank_to_ideal",
  "distance_to_nadir", "rank_from_nadir"
), drop = FALSE]
#>   solution_id distance_to_ideal rank_to_ideal distance_to_nadir rank_from_nadir
#> 1           1         1.0000000             5          1.000000               5
#> 2           2         0.3189474             2          1.179910               1
#> 3           3         0.2508477             1          1.164767               2
#> 4           4         0.4130194             3          1.078791               3
#> 5           5         0.6876632             4          1.015298               4
#> 6           6         1.0000000             5          1.000000               5

Linkage analysis asks whether changes in objective performance correspond to small adjustments or substantial reconfigurations of selected actions. frontier_neighbors() identifies neighbouring solutions in objective space; linkage_distances() and linkage_turnover() relate performance differences to distances or changes in their spatial assignments.

neighbors <- frontier_neighbors(
  alternatives, objectives = c("cost", "spatial")
)

linkage <- linkage_distances(
  alternatives,
  objectives = c("cost", "spatial"),
  pairs = neighbors,
  decision_metric = "jaccard"
)
linkage[, c(
  "from_solution", "to_solution",
  "objective_distance", "decision_distance"
), drop = FALSE]
#>   from_solution to_solution objective_distance decision_distance
#> 1             6           5          0.3150517         0.2934079
#> 2             5           4          0.2952532         0.3417335
#> 3             4           3          0.2565453         0.4233703
#> 4             3           2          0.1729464         0.3609065
#> 5             2           1          0.6843989         0.6743848

turnover <- linkage_turnover(
  alternatives,
  objectives = c("cost", "spatial"),
  pairs = neighbors,
  decision_metric = "jaccard"
)
turnover[, c(
  "from_solution", "to_solution", "objective_distance",
  "decision_distance", "reconfiguration_rate"
), drop = FALSE]
#>   from_solution to_solution objective_distance decision_distance
#> 1             6           5          0.3150517         0.2934079
#> 2             5           4          0.2952532         0.3417335
#> 3             4           3          0.2565453         0.4233703
#> 4             3           2          0.1729464         0.3609065
#> 5             2           1          0.6843989         0.6743848
#>   reconfiguration_rate
#> 1            0.9313008
#> 2            1.1574252
#> 3            1.6502751
#> 4            2.0868113
#> 5            0.9853681

linkage_contrasts() identifies the pairs with the greatest reconfiguration relative to their objective-space separation. linkage_transition() shows which assignments changed between a selected pair of solutions, and selection_consistency() summarises recurring planning-unit states. To keep the README compact, we display the contrast table and only a few rows of planning-unit consistency; detailed transition records remain available from the function itself.

linkage_contrasts(turnover, type = "high_reconfiguration", n = 2)
#>   contrast_rank        contrast_type from_solution to_solution
#> 1             1 high_reconfiguration             3           2
#> 2             2 high_reconfiguration             4           3
#>   objective_distance decision_similarity decision_distance objective_tie
#> 1          0.1729464           0.6390935         0.3609065         FALSE
#> 2          0.2565453           0.5766297         0.4233703         FALSE
#>   reconfiguration_rate changed_assignments changed_planning_units additions
#> 1             2.086811                3822                   2175      1915
#> 2             1.650275                4663                   2655      2324
#>   removals activated_planning_units deactivated_planning_units action_switches
#> 1     1907                      268                        260            1647
#> 2     2339                      316                        331            2008
#>   composition_changes from_cost to_cost delta_cost improvement_cost
#> 1                   0     23173   22765       -408              408
#> 2                   0     23980   23173       -807              807
#>   from_spatial to_spatial delta_spatial improvement_spatial
#> 1         5231       7865          2634               -2634
#> 2         3306       5231          1925               -1925

transition <- linkage_transition(
  alternatives, from = 1, to = 2, objectives = c("cost", "spatial")
)
print(transition)
#> Spatial solution transition
#> From solution: 1 
#> To solution:   2 
#> 
#> Planning units changed: 4571 of 11109 (41.1%)
#> Activated:            168 
#> Deactivated:          149 
#> Action switches:      4254 
#> Composition changes:  0 
#> 
#> Objective changes:
#> - cost: 128 (improvement -128)
#> - spatial: -14620 (improvement 14620)
#> 
#> Use `$summary`, `$objectives`, `$transitions`, `$actions`, or `$state_matrix` for details.

head(selection_consistency(alternatives))
#>      pu solution_group n_solutions selected_frequency managed_frequency
#> 1     1            all           6                1.0               1.0
#> 2    10            all           6                1.0               1.0
#> 3   100            all           6                1.0               1.0
#> 4  1000            all           6                1.0               1.0
#> 5 10000            all           6                0.0               0.0
#> 6 10001            all           6                0.5               0.5
#>   unmanaged_frequency dominant_state dominant_frequency dominant_tie n_states
#> 1                 0.0            BAU          0.8333333        FALSE        2
#> 2                 0.0            BAU          0.8333333        FALSE        2
#> 3                 0.0            BAU          0.8333333        FALSE        2
#> 4                 0.0  Afforestation          0.8333333        FALSE        2
#> 5                 1.0      unmanaged          1.0000000        FALSE        1
#> 6                 0.5      unmanaged          0.5000000        FALSE        3
#>   variable   entropy normalized_entropy
#> 1     TRUE 0.4505612          0.6500224
#> 2     TRUE 0.4505612          0.6500224
#> 3     TRUE 0.4505612          0.6500224
#> 4     TRUE 0.4505612          0.6500224
#> 5    FALSE 0.0000000          0.0000000
#> 6     TRUE 1.0114043          0.9206198

A complete executable example is provided in examples/meseta_iberica.R. The simulated workflow introduces additional optimisation methods and diagnostics on a smaller dataset.

Learn more

The package provides several prepared examples:

Loader Example
load_sim_multiaction() 64 spatial cells, two features and two actions for quick examples
load_meseta() The Meseta Ibérica inputs used in this README
load_ecosystem_services() Planning units and four raster layers for the integrated-planning vignette

Browse the function reference for details of create_problem(), add_actions(), add_effects(), add_constraint_targets_absolute(), the set_method_*() functions, and solve(). Tools for inspecting alternatives are organised as frontier_*(), selection_*(), and linkage_*().

If you find a bug or would like to suggest an improvement, please open an issue.

References

Cánibe Iglesias, M., Hermoso, V., Azevedo, J. C., Campos, J. C., Salgado-Rojas, J., Sil, A., & Regos, A. (2025). Integrating multiple landscape management strategies to optimise conservation under climate and planning scenarios: a case study in the Iberian Peninsula. Ecosystem Services, 74, 101742. doi:10.1016/j.ecoser.2025.101742.

Tallis, H., Fargione, J., Game, E., et al. (2021). Prioritizing actions: spatial action maps for conservation. Annals of the New York Academy of Sciences, 1505(1), 118-141. doi:10.1111/nyas.14651.

Salgado-Rojas, J., Hermoso, V., & Alvarez-Miranda, E. (2023). prioriactions: Multi-action management planning in R. Methods in Ecology and Evolution. doi:10.1111/2041-210X.14220.

Neubert, S., McGowan, J., Metcalfe, K., et al. (2025). Multiple-use spatial planning for sustainable development and conservation. Trends in Ecology & Evolution, 40(11), 1126-1142. doi:10.1016/j.tree.2025.09.007.

These binaries (installable software) and packages are in development.
They may not be fully stable and should be used with caution. We make no claims about them.
Health stats visible at Monitor.