1 MARL for Action based PES design

Action based PES is when incentives are rewarded for actions that landowners or herders take to reduce land use pressures on ecosystems. In this model, the action is the chosen herd-size category,and the rangeland condition categories are the states.

This is a simplified version with 6 agents and 25cells(5x5). Each agent has 2 designated pastures that they switch between each iteration. In an emergency they are allowed to choose a third pasture,that is not a designated pasture of someone else’s. Co-location has a 0.7 penalty, so there is incentive to avoid overlapping with other agents in the same cell(pasture). Because simulation would be expensive on a continuous herd size, the model for now has four stocking categories as actions.Agents can’t change herd size dramatically between two iterations,and hence there is an adjacent-category constraint:an agent may stay in the same category or move by one category per iteration, but cannot jump directly from<200to>600or vice versa. There is a 5 year ecosystem memory with an exponential decay weight(last year’s condition is the most influential on this year’s). Information sharing is only for the purpose of estimating neighbourhood density, which then enables agents to avoid a co-location penalty.

Independent Q-learning algorithm assumptions: each agent has their own independent policy. Agents interact indirectly i.e. only through off-policy tables to estimate neighbourhood density. There is no centralised learning.

1.1 Model definition

The learning state for agent \(i\) at time \(t\) is

\[ s_{i,t} = ( \text{current cell}, \text{rangeland condition}, \text{weather}, \text{neighbour density}, \text{previous action} ). \]

Each agent learns with

\[ Q(s,a) \leftarrow Q(s,a) + \alpha \left[ r + \gamma \max_{a'}Q(s',a') - Q(s,a) \right]. \]

The behavioural policy is softmax over the actions that are feasible from the previous stocking category:

\[ \pi(a|s) = \operatorname{softmax} \left( Q(s,a)/\tau \right), \qquad a \in \mathcal{A}(a_{t-1}). \]

Direct livestock income is adjusted by expected stock survival indexed by current rangeland condition and stocking action. Stock survival affects current reward only in this baseline; surviving herd size is logged but not carried forward.

1.2 Spatial rules

  • There are six agents and 25 cells.
  • Each agent has two unique designated cells.
  • Under ordinary conditions, the agent alternates deterministically between its two designated cells.
  • Emergency movement is activated only when:
    1. current weather is drought; and
    2. both designated cells are degraded.
  • An emergency destination:
    • cannot be any agent’s designated cell;
    • must be ecologically better than both designated cells;
    • incurs travel and social costs.
  • Travel and social costs receive equal weight: 0.5 + 0.5.
  • Social cost measures passage through other agents’ designated cells on the route to the emergency destination.
  • If multiple agents occupy the same destination cell, their rewards receive a co-location penalty.

1.3 Information sharing

All agents observe a shared public ledger containing the current locations and current herd sizes of all other agents. This information is used to calculate neighbour density. The shared ledger is deliberately simple in this version.

1.4 Ecological memory

Every cell carries its own persistent 5-step exponential moving average (EMA) of ecological degradation severity. The stochastic ecological signal is drawn from the CSV-derived transition kernel in Step 5, which now permits constrained recovery as well as deterioration. This includes unassigned cells: if an emergency herder grazed a cell several iterations earlier, that legacy remains in the cell EMA and fades exponentially rather than being reset when the cell is ungrazed. The newest transition outcome receives the largest weight. The EMA defines the categorical rangeland_condition observed by agents. If one simulation step represents one year, this is approximately a five-year ecological memory; if one step represents a season, it is a five-season memory.

2 . Setup

library(tidyverse)
## Warning: package 'tidyverse' was built under R version 4.3.3
## Warning: package 'tibble' was built under R version 4.3.2
## Warning: package 'tidyr' was built under R version 4.3.3
## Warning: package 'dplyr' was built under R version 4.3.3
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.4     ✔ readr     2.2.0
## ✔ forcats   1.0.1     ✔ stringr   1.6.0
## ✔ ggplot2   4.0.0     ✔ tibble    3.2.1
## ✔ lubridate 1.9.5     ✔ tidyr     1.3.1
## ✔ purrr     1.2.0     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(purrr)
library(ggplot2)
set.seed(128)
#Grid
grid_nx<-5L
grid_ny<-5L
#Rangeland degradation severity from best to worst
r_states<-c("intact","good","poor","degraded")
degradation_score<-c(intact=0,good=1,poor=2,degraded=3)
#Four discrete stocking-rate actions
actions<-c("<200","200-400","400-600",">600")
#"none" is when there is no grazing in the cell
ecological_footprints<-c("none",actions)
#Representative herd size used for information sharing and combined pressure
action_herd_size<-c("<200"=100,"200-400"=300,"400-600"=500,">600"=700)
#Starting stocking categories as placeholders
initial_actions<-c(A1="<200",A2="200-400",A3="200-400",A4="400-600",A5="400-600",A6=">600")
#Direct herd income before survival adjustment
herd_income<-c("<200"=0.8,"200-400"=2.0,"400-600"=2.5,">600"=3.5)
#Action-based PES
pes_payment<-c("<200"=1.5,"200-400"=1.0,"400-600"=0.4,">600"=0.0)
#Expected stock survival by current rangeland condition and stocking action
stock_survival<-matrix(
  c(
    0.98,0.95,0.94,0.93,#intact
    0.96,0.94,0.93,0.92,#good
    0.84,0.83,0.8,0.7,#poor
    0.83,0.8,0.7,0.6#degraded
  ),
  nrow=4,
  byrow=TRUE,
  dimnames=list(r_states,actions)
)
#Weather
weather_prob<-c(drought=0.30,normal=0.50,wet=0.20)
#Q-learning paras
alpha<-0.20
gamma<-0.95
tau<-1.00
#Co-location penalty
colocation_multiplier<-0.70
#Emergency movement costs
travel_cost_weight<-0.50
social_cost_weight<-0.50
emergency_cost_scale<-1.00
#Ecosystem memory
memory_window<-5L
ema_alpha<-2/(memory_window+1)

3 . Create the 5 × 5 landscape

Cells are indexed from 1 to 25. All cells begin intact.

4 . Assign two designated cells to each of six agents

The 12 designated cells are mutually exclusive for each agent.

designated_cells<-tribble(
  ~agent_id, ~slot, ~cell_id,
  "A1", 1L,  1L,
  "A1", 2L,  7L,
  "A2", 1L,  5L,
  "A2", 2L,  9L,
  "A3", 1L, 21L,
  "A3", 2L, 17L,
  "A4", 1L, 25L,
  "A4", 2L, 19L,
  "A5", 1L, 11L,
  "A5", 2L, 12L,
  "A6", 1L, 14L,
  "A6", 2L, 15L
)
designated_lookup<-designated_cells %>%
  arrange(agent_id, slot) %>%
  group_by(agent_id) %>%
  summarise(
    cell_1 = first(cell_id),
    cell_2 = last(cell_id),
    .groups = "drop"
  )
designated_owner<-designated_cells %>%
  distinct(cell_id, agent_id)
unassigned_cells<-setdiff(
  land_grid_initial$cell_id,
  designated_cells$cell_id
)

5 . Weather

Weather is exogenous and shared by the whole landscape at each time step.

draw_weather<-function(n = 1L) {
  sample(
    names(weather_prob),
    size = n,
    replace = TRUE,
    prob = weather_prob
  )
}

6 . Expert derived ecological transition model

The expert derived transition table contains joint relationship between climate, stocking category, and three pasture outcomes; intact,good,degraded. There is no directly elicited poor state. We derive poor state by applying the exponential moving average (EMA) function, which produces a continuous number between 0 and 3. poor is assigned for values between 1.5 and 2.5.

Due to time limitations on experts’ time and cognitive load, not all starting conditions for the full transition table was possible to be elicited in the field. Only the starting condition intact is fully elicited. The remaining starting conditions (good, poor and degraded) are therefore are derived from the intact starting condition values i.e. purely assumption based at this stage.

The missing current-condition rows are constructed as follows for each action x weather combination:

\[ \mathbf p_G=\frac{2}{3}\mathbf p_I+\frac{1}{3}\mathbf p_D, \]

and

\[ \mathbf p_P^{raw}=\frac{1}{3}\mathbf p_G+\frac{2}{3}\mathbf p_D. \]

For current poor and degraded conditions, a direct transition signal to intact is prohibited. Any raw intact mass in the poor interpolation is transferred to good. Thus recovery from poor or degraded occurs through good rather than jumping directly to intact.

For a currently degraded cell:

For none, which is also absent from the elicitation, we have also made assumptions for now as a placeholder; intact-start persistence anchors are 0.99 under wet and normal and 0.01 under drought. The remaining probability is split between good and degraded using the same-weather <200 empirical good:degraded ratio.

#Empirical joint weather x stocking transition probabilities. These are the mean of 140 individual estimates. 
climate_transition_source<-tribble(
  ~weather,~action,~intact,~good,~degraded,
  "wet","<200",0.8665,0.0963,0.0372,
  "wet","200-400",0.7425,0.1576,0.0999,
  "wet","400-600",0.5945,0.2278,0.1777,
  "wet",">600",0.4349,0.2689,0.2962,
  "normal","<200",0.7438,0.1837,0.0725,
  "normal","200-400",0.5856,0.2404,0.1740,
  "normal","400-600",0.3926,0.2906,0.3168,
  "normal",">600",0.2581,0.2811,0.4608,
  "drought","<200",0.4906,0.3077,0.2017,
  "drought","200-400",0.3481,0.3067,0.3452,
  "drought","400-600",0.2327,0.2730,0.4943,
  "drought",">600",0.1261,0.2255,0.6484
)
transition_signal_states<-c("intact","good","degraded")
make_prob_vector<-function(intact,good,degraded){
  p<-c(intact=intact,good=good,poor=0,degraded=degraded)
  if(any(!is.finite(p))||any(p<0)||any(p>1)||abs(sum(p)-1)>1e-10){
    stop("Invalid ecological transition probability vector.")
  }
  p
}
empirical_intact_prob<-function(action,weather){
  stopifnot(action%in%actions,weather%in%names(weather_prob))
  row<-climate_transition_source[
    climate_transition_source$action==action&
      climate_transition_source$weather==weather,
    ,drop=FALSE
  ]
  if(nrow(row)!=1L)stop("Expected one empirical action x weather row.")
  make_prob_vector(row$intact[[1]],row$good[[1]],row$degraded[[1]])
}
#Explicit no-grazing anchors for an intact current condition.
none_intact_persistence<-c(wet=0.99,normal=0.99,drought=0.01)
make_none_intact_prob<-function(weather){
  stopifnot(weather%in%names(weather_prob))
  p_intact<-unname(none_intact_persistence[[weather]])
  reference<-empirical_intact_prob("<200",weather)
  residual<-1-p_intact
  gd_total<-reference[["good"]]+reference[["degraded"]]
  make_prob_vector(
    intact=p_intact,
    good=residual*reference[["good"]]/gd_total,
    degraded=residual*reference[["degraded"]]/gd_total
  )
}
intact_start_prob<-function(action,weather){
  if(action=="none")make_none_intact_prob(weather) else empirical_intact_prob(action,weather)
}
#Wet recovery curve for a currently degraded grazed cell.
wet_recovery_shape<-climate_transition_source%>%
  filter(weather=="wet")%>%
  mutate(recovery_favourability=1-degraded)%>%
  arrange(match(action,actions))
wet_fav_low<-wet_recovery_shape$recovery_favourability[
  wet_recovery_shape$action==">600"
]
wet_fav_high<-wet_recovery_shape$recovery_favourability[
  wet_recovery_shape$action=="<200"
]
wet_recovery_shape<-wet_recovery_shape%>%
  mutate(
    shape=(recovery_favourability-wet_fav_low)/(wet_fav_high-wet_fav_low),
    degraded_to_good=0.01+shape*(0.50-0.01)
  )
degraded_none_recovery<-c(wet=0.80,normal=0.66,drought=0.01)
degraded_recovery_prob<-function(action,weather){
  stopifnot(action%in%ecological_footprints,weather%in%names(weather_prob))
  if(action=="none")return(unname(degraded_none_recovery[[weather]]))
  if(weather!="wet")return(0)
  row<-wet_recovery_shape[wet_recovery_shape$action==action,,drop=FALSE]
  if(nrow(row)!=1L)stop("Expected one wet recovery row.")
  row$degraded_to_good[[1]]
}
degraded_start_prob<-function(action,weather){
  recovery<-degraded_recovery_prob(action,weather)
  make_prob_vector(intact=0,good=recovery,degraded=1-recovery)
}
apply_recovery_constraint<-function(prob,current_condition){
  stopifnot(current_condition%in%r_states)
  if(current_condition%in%c("poor","degraded")){
    prob[["good"]]<-prob[["good"]]+prob[["intact"]]
    prob[["intact"]]<-0
  }
  if(any(prob<0)||any(prob>1)||abs(sum(prob)-1)>1e-10){
    stop("Recovery constraint produced an invalid probability vector.")
  }
  prob
}
make_transition_matrix<-function(action,weather){
  stopifnot(action%in%ecological_footprints,weather%in%names(weather_prob))
  p_intact<-intact_start_prob(action,weather)
  p_degraded<-apply_recovery_constraint(
    degraded_start_prob(action,weather),
    "degraded"
  )
  p_good<-(2/3)*p_intact+(1/3)*p_degraded
  p_poor_raw<-(1/3)*p_good+(2/3)*p_degraded
  p_poor<-apply_recovery_constraint(p_poor_raw,"poor")
  M<-rbind(
    intact=p_intact,
    good=p_good,
    poor=p_poor,
    degraded=p_degraded
  )
  M<-M[r_states,r_states,drop=FALSE]
  if(any(!is.finite(M))||any(M<0)||any(M>1)||any(abs(rowSums(M)-1)>1e-10)){
    stop("Invalid ecological transition matrix.")
  }
  M
}
base_transition<-setNames(
  lapply(ecological_footprints,function(a){
    setNames(
      lapply(names(weather_prob),function(w){make_transition_matrix(a,w)}),
      names(weather_prob)
    )
  }),
  ecological_footprints
)
sample_next_condition<-function(current_condition,action,weather){
  transition_matrix<-base_transition[[action]][[weather]]
  probabilities<-transition_matrix[current_condition,,drop=TRUE]
  sample(names(probabilities),size=1,prob=probabilities)
}

7 . Five-step ecosystem memory

This is cell-level ecological memory. Every one of the 25 cells stores its own persistent EMA, including cells that are currently ungrazed or unassigned.

For the ecological update to time \(t+1\), first draw a stochastic transition outcome \(z_{t+1}\) from the base transition using \(C_t\), \(A_t\), and \(W_{t+1}\). Then update the cell EMA:

\[ E_{t+1} = \alpha_E z_{t+1} + (1-\alpha_E)E_t. \]

The newest outcome receives weight \(\alpha_E\); older ecological history decays geometrically. The categorical rangeland_condition observed by the agents is then defined from the EMA thresholds. This is the mechanism that defines the rangeland-condition component already present in \(s\).

condition_from_ema<-function(x){
  case_when(
    x<0.5~"intact",
    x<1.5~"good",
    x<2.5~"poor",
    TRUE~"degraded"
  )
}
update_cell_memory<-function(old_ema,transition_outcome,alpha=ema_alpha){
  outcome_score<-unname(degradation_score[[transition_outcome]])
  new_ema<-alpha*outcome_score+(1-alpha)*old_ema
  list(ema_score=new_ema,condition=condition_from_ema(new_ema))
}

8 . Agent objects and Q tables

Each agent has its own Q table. All agents start with the same learning parameters in this basic version.

empty_q_table<-function() {
  matrix(
    numeric(0),
    nrow = 0,
    ncol = length(actions),
    dimnames = list(
      character(0),
      actions
    )
  )
}
create_agent<-function(
    agent_id,
    initial_action = initial_actions[[agent_id]],
    alpha_value = alpha,
    gamma_value = gamma,
    tau_value = tau
) {
  stopifnot(initial_action %in% actions)
  first_cell<-designated_lookup %>%
    filter(.data$agent_id==.env$agent_id) %>%
    pull(cell_1)
  list(
    id = agent_id,
    alpha = alpha_value,
    gamma = gamma_value,
    tau = tau_value,
    current_cell = first_cell,
    previous_action = initial_action,
    herd_size = unname(action_herd_size[[initial_action]]),
    Q = empty_q_table()
  )
}
agents_initial<-setNames(
  lapply(
    names(initial_actions),
    function(id) {
      create_agent(
        agent_id = id,
        initial_action = initial_actions[[id]]
      )
    }
  ),
  names(initial_actions)
)
ensure_q_state<-function(agent, state_key) {
  if (!state_key %in% rownames(agent$Q)) {
    new_row<-matrix(
      0,
      nrow = 1,
      ncol = length(actions),
      dimnames = list(
        state_key,
        actions
      )
    )
    agent$Q<-rbind(agent$Q, new_row)
  }
  agent
}

9 . Shared information for neighbourhood density estimation

The public information ledger contains:

Neighbour density is the number of other agents in the neighbourhood of the focal agent (one-cell radius in x and y).

share_agent_information<-function(
    agents,
    land_grid
) {
  map_dfr(
    agents,
    function(agent) {
      loc<-land_grid %>%
        filter(cell_id == agent$current_cell)
      tibble(
        agent_id = agent$id,
        current_cell = agent$current_cell,
        x = loc$x,
        y = loc$y,
        previous_action = agent$previous_action,
        herd_size = agent$herd_size
      )
    }
  )
}
compute_neighbour_density<-function(
    agent_id,
    shared_info,
    radius = 1L
) {
  focal<-shared_info %>%
    filter(.data$agent_id==.env$agent_id)
  others<-shared_info %>%
    filter(.data$agent_id!=.env$agent_id)
  sum(
    abs(others$x - focal$x) <= radius &
      abs(others$y - focal$y) <= radius
  )
}

10 . Q-learning state encoding

The exact state requested for policy learning is

\[ ( \text{current cell}, \text{rangeland condition}, \text{weather}, \text{neighbour density}, \text{previous action} ). \]

encode_state<-function(
    current_cell,
    rangeland_condition,
    weather,
    neighbour_density,
    previous_action
) {
  paste(
    paste0("cell=", current_cell),
    paste0("land=", rangeland_condition),
    paste0("weather=", weather),
    paste0("neigh=", neighbour_density),
    paste0("prev=", previous_action),
    sep = "|"
  )
}
observe_agent_state<-function(
    agent,
    land_grid,
    weather,
    shared_info
) {
  current_condition<-land_grid %>%
    filter(cell_id == agent$current_cell) %>%
    pull(rangeland_condition)
  neighbour_density<-compute_neighbour_density(
    agent_id = agent$id,
    shared_info = shared_info
  )
  encode_state(
    current_cell = agent$current_cell,
    rangeland_condition = current_condition,
    weather = weather,
    neighbour_density = neighbour_density,
    previous_action = agent$previous_action
  )
}

11 . Softmax policy with stocking-category anchoring

Action choice is based only on Q values for the currently observed state, but the feasible action set is constrained by the previous stocking category. The agent may stay in the same category or move one category up/down in one iteration. This is the discrete approximation to the stocking-change constraint.

#Discrete analogue of a <=10% stocking-change anchor:
#an agent may stay in the same category or move only one category up/down. For example, an agent can not go from 600 goats down to 100 in one move. They always have to move through the middle actions i.e. 200-400 and 400-600. 
allowed_actions<-function(previous_action) {
  switch(
    previous_action,
    "<200" = c("<200", "200-400"),
    "200-400" = c(
      "<200",
      "200-400",
      "400-600"
    ),
    "400-600" = c(
      "200-400",
      "400-600",
      ">600"
    ),
    ">600" = c("400-600", ">600"),
    stop("Unknown previous action: ", previous_action)
  )
}
softmax<-function(q_values, tau = 1) {
  stopifnot(tau > 0)
  z<-(q_values - max(q_values)) / tau
  exp_z<-exp(z)
  exp_z / sum(exp_z)
}
select_action_q<-function(agent, state_key) {
  agent<-ensure_q_state(agent, state_key)
  feasible_actions<-allowed_actions(
    agent$previous_action
  )
  q_values<-agent$Q[
    state_key,
    feasible_actions,
    drop = TRUE
  ]
  probabilities<-softmax(
    q_values,
    tau = agent$tau
  )
  chosen_action<-sample(
    feasible_actions,
    size = 1,
    prob = probabilities
  )
  #Keep a full action probability vector for diagnostics.
  full_probabilities<-setNames(
    rep(0, length(actions)),
    actions
  )
  full_probabilities[feasible_actions]<-probabilities
  list(
    agent = agent,
    action = chosen_action,
    feasible_actions = feasible_actions,
    probabilities = full_probabilities,
    q_values = q_values
  )
}

12 . Normal movement between two designated cells

Normal location is deterministic: agents alternate between their two designated cells.

normal_cell_for_step<-function(
    agent_id,
    step
) {
  row<-designated_lookup %>%
    filter(.data$agent_id==.env$agent_id)
  if (step %% 2L == 1L) {
    row$cell_1
  } else {
    row$cell_2
  }
}

13 . Emergency trigger and emergency destination

The emergency gate is intentionally simple:

\[ \text{Emergency}_{i,t} = \mathbf{1} ( W_t = \text{drought} \land C_{i,1} = \text{degraded} \land C_{i,2} = \text{degraded} ). \]

There is no stochastic mobility propensity in this basic version.

13.1 .1 Straight-line route through the grid

The route approximation is used only to calculate social trespass cost.

cells_on_straight_line<-function(
    origin_cell,
    destination_cell,
    land_grid,
    n_points = 101L
) {
  origin<-land_grid %>%
    filter(cell_id == origin_cell)
  destination<-land_grid %>%
    filter(cell_id == destination_cell)
  xs<-seq(origin$x, destination$x, length.out = n_points)
  ys<-seq(origin$y, destination$y, length.out = n_points)
  tibble(
    x = round(xs),
    y = round(ys)
  ) %>%
    distinct(x, y) %>%
    left_join(
      land_grid %>% select(cell_id, x, y),
      by = c("x", "y")
    ) %>%
    filter(!is.na(cell_id)) %>%
    pull(cell_id)
}

13.2 .2 Travel and social costs

Travel cost is normalized Euclidean distance.Social cost is the fraction of intermediate route cells that are designated to another agent. The destination itself must be unassigned, so social cost comes from crossing other people’s pastures rather than occupying their designated destination.

emergency_movement_cost<-function(
    agent_id,
    origin_cell,
    destination_cell,
    land_grid
) {
  origin<-land_grid %>%
    filter(cell_id == origin_cell)
  destination<-land_grid %>%
    filter(cell_id == destination_cell)
  distance<-sqrt(
    (origin$x - destination$x)^2 +
      (origin$y - destination$y)^2
  )
  max_distance<-sqrt(
    (grid_nx - 1)^2 +
      (grid_ny - 1)^2
  )
  travel_cost<-distance / max_distance
  path_cells<-cells_on_straight_line(
    origin_cell,
    destination_cell,
    land_grid
  )
  intermediate_cells<-setdiff(
    path_cells,
    c(origin_cell, destination_cell)
  )
  other_designated<-designated_cells %>%
    filter(.data$agent_id!=.env$agent_id) %>%
    pull(cell_id) %>%
    unique()
  social_cost<-if (
    length(intermediate_cells) == 0
  ) {
    0
  } else {
    mean(
      intermediate_cells %in% other_designated
    )
  }
  total_cost<-travel_cost_weight*travel_cost+social_cost_weight*social_cost
  list(
    travel_cost = travel_cost,
    social_cost = social_cost,
    total_cost = total_cost,
    path_cells = path_cells
  )
}

13.3 .3 Emergency candidate set

An emergency destination must be:

  • be unassigned cell that is not a designated cell for another agent;
  • be better than both of the agent’s own designated cells;
  • therefore be non-degraded when both designated cells are degraded.

Among valid candidates, the basic version chooses the cell with the lowest combined emergency cost.

emergency_active<-function(
    agent_id,
    weather,
    land_grid
) {
  if (weather != "drought") {
    return(FALSE)
  }
  designated<-designated_cells %>%
    filter(.data$agent_id==.env$agent_id) %>%
    pull(cell_id)
  conditions<-land_grid %>%
    filter(cell_id %in% designated) %>%
    arrange(match(cell_id, designated)) %>%
    pull(rangeland_condition)
  length(conditions) == 2 &&
    all(conditions == "degraded")
}
select_emergency_cell<-function(
    agent_id,
    origin_cell,
    land_grid
) {
  designated<-designated_cells %>%
    filter(.data$agent_id==.env$agent_id) %>%
    pull(cell_id)
  designated_scores<-land_grid %>%
    filter(cell_id %in% designated) %>%
    pull(rangeland_condition) %>%
    degradation_score[.]
  #Candidate must be strictly better than BOTH designated cells.
  max_allowed_score<-min(designated_scores)
  candidates<-land_grid %>%
    filter(
      cell_id %in% unassigned_cells,
      degradation_score[rangeland_condition] < max_allowed_score
    )
  if (nrow(candidates) == 0) {
    return(NULL)
  }
  candidate_costs<-map_dfr(
    candidates$cell_id,
    function(destination_cell) {
      cost<-emergency_movement_cost(
        agent_id = agent_id,
        origin_cell = origin_cell,
        destination_cell = destination_cell,
        land_grid = land_grid
      )
      tibble(
        destination_cell = destination_cell,
        travel_cost = cost$travel_cost,
        social_cost = cost$social_cost,
        total_cost = cost$total_cost
      )
    }
  )
  best<-candidate_costs %>%
    filter(total_cost == min(total_cost)) %>%
    slice_sample(n = 1)
  list(
    cell_id = best$destination_cell,
    travel_cost = best$travel_cost,
    social_cost = best$social_cost,
    total_cost = best$total_cost,
    all_candidates = candidate_costs
  )
}
resolve_location<-function(
    agent_id,
    step,
    weather,
    land_grid
) {
  normal_cell<-normal_cell_for_step(
    agent_id,
    step
  )
  if (!emergency_active(
    agent_id,
    weather,
    land_grid
  )) {
    return(list(
      cell_id = normal_cell,
      movement_mode = "designated",
      emergency_cost = 0,
      travel_cost = 0,
      social_cost = 0
    ))
  }
  emergency<-select_emergency_cell(
    agent_id = agent_id,
    origin_cell = normal_cell,
    land_grid = land_grid
  )
  if (is.null(emergency)) {
    return(list(
      cell_id = normal_cell,
      movement_mode = "emergency_no_destination",
      emergency_cost = 0,
      travel_cost = 0,
      social_cost = 0
    ))
  }
  list(
    cell_id = emergency$cell_id,
    movement_mode = "emergency",
    emergency_cost = emergency$total_cost,
    travel_cost = emergency$travel_cost,
    social_cost = emergency$social_cost
  )
}

14 . Action-based PES reward

Immediate reward contains:

  1. herd income multiplied by expected stock survival for the current rangeland condition and stocking action;
  2. action-based PES payment;
  3. co-location penalty;
  4. emergency movement cost when applicable.

Stock survival is the only multiplier on direct livestock income in this baseline. Stock survival affects current-period reward, while the intended herd size selected by the action still determines current grazing pressure. Surviving herd size is recorded diagnostically but is not yet carried forward as an endogenous herd-state variable.

Importantly, this reward enters the Q update directly, so it influences future policy.

get_stock_survival<-function(action,rangeland_condition){
  unname(stock_survival[rangeland_condition,action])
}
compute_reward<-function(action,rangeland_condition,n_agents_same_cell=1L,emergency_cost=0){
  survival_rate<-get_stock_survival(action=action,rangeland_condition=rangeland_condition)
  intended_herd_size<-unname(action_herd_size[[action]])
  expected_surviving_herd<-intended_herd_size*survival_rate
  expected_stock_loss<-intended_herd_size-expected_surviving_herd
  livestock_income<-herd_income[[action]]*survival_rate
  pes_reward<-pes_payment[[action]]
  gross_reward<-livestock_income+pes_reward
  interaction_multiplier<-colocation_multiplier^max(n_agents_same_cell-1L,0L)
  reward_after_interaction<-gross_reward*interaction_multiplier
  net_reward<-reward_after_interaction-emergency_cost_scale*emergency_cost
  list(
    net_reward=net_reward,
    livestock_income=livestock_income,
    pes_reward=pes_reward,
    gross_reward=gross_reward,
    interaction_multiplier=interaction_multiplier,
    emergency_cost=emergency_cost,
    stock_survival=survival_rate,
    intended_herd_size=intended_herd_size,
    expected_surviving_herd=expected_surviving_herd,
    expected_stock_loss=expected_stock_loss
  )
}

15 . Q-learning update

The policy is updated through the standard one-step temporal-difference rule.

update_q<-function(
    agent,
    state_key,
    action,
    reward,
    next_state_key
) {
  agent<-ensure_q_state(agent, state_key)
  agent<-ensure_q_state(agent, next_state_key)
  q_old<-agent$Q[
    state_key,
    action
  ]
#At s_(t+1), the action just taken at t is the previous stocking category.Therefore only the same/adjacent actions are feasible in the Bellman max.
  next_feasible_actions<-allowed_actions(action)
  q_next_max<-max(
    agent$Q[
      next_state_key,
      next_feasible_actions,
      drop = TRUE
    ]
  )
  td_target<-reward+agent$gamma*q_next_max
  td_error<-td_target - q_old
  q_new<-q_old+agent$alpha*td_error
  agent$Q[
    state_key,
    action
  ]<-q_new
  list(
    agent = agent,
    q_old = q_old,
    q_new = q_new,
    td_target = td_target,
    td_error = td_error,
    next_feasible_actions = next_feasible_actions
  )
}

16 . Convert combined herd pressure to an ecological footprint

If several agents occupy the same cell, their herd sizes are added before the ecological transition is sampled to account for the total grazing pressure in each cell.

herd_size_to_action<-function(total_herd_size) {
  case_when(
    total_herd_size <= 0 ~ "none",
    total_herd_size < 200 ~ "<200",
    total_herd_size < 400 ~ "200-400",
    total_herd_size < 600 ~ "400-600",
    TRUE ~ ">600"
  )
}

17 . Update the landscape once per step

Each of the 25 cells is updated once per iteration and keeps its own persistent ecological EMA. For grazed cells, the stochastic ecological transition uses the combined herd pressure of all agents present. For ungrazed cells, ecological footprint is none: grazing pressure is zero, but weather still affects the stochastic transition outcome. A previously grazed emergency cell therefore retains its ecological history after the herder leaves. Every update is written to the cell-level ledger.

update_landscape<-function(land_grid,action_log,transition_weather,step){
  pressure_by_cell<-action_log%>%
    group_by(cell_id)%>%
    summarise(total_herd_size=sum(herd_size),n_agents=n(),.groups="drop")%>%
    mutate(ecological_footprint=map_chr(total_herd_size,herd_size_to_action))
  new_grid<-land_grid
  ledger_list<-vector("list",nrow(land_grid))
  for(j in seq_len(nrow(land_grid))){
    cell<-land_grid$cell_id[j]
    current_condition<-land_grid$rangeland_condition[j]
    old_ema<-land_grid$condition_ema_score[j]
    pressure<-pressure_by_cell%>%filter(.data$cell_id==cell)
    if(nrow(pressure)==0){
      ecological_footprint<-"none"
      total_herd_size<-0
      n_agents<-0
    }else{
      ecological_footprint<-pressure$ecological_footprint[[1]]
      total_herd_size<-pressure$total_herd_size[[1]]
      n_agents<-pressure$n_agents[[1]]
    }
    transition_outcome<-sample_next_condition(
      current_condition=current_condition,
      action=ecological_footprint,
      weather=transition_weather
    )
    memory_update<-update_cell_memory(old_ema=old_ema,transition_outcome=transition_outcome)
    new_grid$condition_ema_score[j]<-memory_update$ema_score
    new_grid$rangeland_condition[j]<-memory_update$condition
    ledger_list[[j]]<-tibble(
      step=step,
      from_time=step,
      to_time=step+1L,
      cell_id=cell,
      transition_weather=transition_weather,
      condition_before=current_condition,
      ecological_footprint=ecological_footprint,
      transition_outcome=transition_outcome,
      ema_before=old_ema,
      ema_after=memory_update$ema_score,
      condition_after=memory_update$condition,
      total_herd_size=total_herd_size,
      n_agents=n_agents
    )
  }
  list(land_grid=new_grid,ledger=bind_rows(ledger_list))
}

18 . One simultaneous MARL step

This is the core multi-agent step.

Order of operations:

  1. Resolve each agent’s current location using the current EMA-defined landscape and current weather.
  2. Share current locations and previous stocking-category herd sizes.
  3. Construct \(s_t\).
  4. Restrict each agent to the same or adjacent stocking category relative to \(A_{t-1}\), then sample \(A_t\) from softmax Q over that feasible set.
  5. Calculate shared-cell interaction counts and immediate rewards from the current landscape.
  6. Update every cell from \(t\) to \(t+1\): combined \(A_t\) (or none) + \(W_{t+1}\) generates a stochastic transition outcome, which updates that cell’s persistent EMA and defines \(C_{t+1}\).
  7. Resolve next locations using \(W_{t+1}\) and the updated EMA-defined landscape.
  8. Construct \(s_{t+1}\).
  9. Update Q, maximizing only over actions feasible from \(A_t\).
  10. Apply one Q-learning update to every agent.
  11. Update each agent’s previous action and herd size.
  12. Write agent and cell ledgers.
marl_step<-function(
    agents,
    land_grid,
    weather,
    next_weather,
    step
) {
#------------------------------------------------------------
#A. Resolve locations in current state s_t
#------------------------------------------------------------
  location_results<-map(
    names(agents),
    ~ resolve_location(
      agent_id = .x,
      step = step,
      weather = weather,
      land_grid = land_grid
    )
  )
  names(location_results)<-names(agents)
  for (id in names(agents)) {
    agents[[id]]$current_cell<-location_results[[id]]$cell_id
  }
#------------------------------------------------------------
#B. Information sharing before action selection
#------------------------------------------------------------
  shared_info<-share_agent_information(
    agents,
    land_grid
  )
#------------------------------------------------------------
#C. Observe s_t
#s_t = (cell, EMA-defined condition, weather, neighbour density, previous stocking category)
#------------------------------------------------------------
  state_keys<-setNames(
    map_chr(
      agents,
      observe_agent_state,
      land_grid = land_grid,
      weather = weather,
      shared_info = shared_info
    ),
    names(agents)
  )
#------------------------------------------------------------
#D. Choose A_t simultaneously, constrained by A_(t-1)
#------------------------------------------------------------
  action_results<-map2(
    agents,
    state_keys,
    select_action_q
  )
#Preserve any newly created Q rows.
  for (id in names(agents)) {
    agents[[id]]<-action_results[[id]]$agent
  }
  chosen_actions<-setNames(
    map_chr(action_results, "action"),
    names(agents)
  )
  herd_sizes<-action_herd_size[chosen_actions]
#------------------------------------------------------------
#E. Build simultaneous action log
#------------------------------------------------------------
  action_log<-tibble(
    agent_id = names(agents),
    cell_id = map_int(agents, "current_cell"),
    previous_action = map_chr(agents, "previous_action"),
    action = unname(chosen_actions),
    herd_size = as.numeric(herd_sizes),
    state_key = unname(state_keys),
    feasible_actions = map_chr(
      action_results,
      ~ paste(.x$feasible_actions, collapse = " | ")
    ),
    chosen_probability = map2_dbl(
      action_results,
      chosen_actions,
      ~ unname(.x$probabilities[[.y]])
    ),
    movement_mode = map_chr(location_results, "movement_mode"),
    emergency_cost = map_dbl(location_results, "emergency_cost"),
    travel_cost = map_dbl(location_results, "travel_cost"),
    social_cost = map_dbl(location_results, "social_cost")
  ) %>%
    group_by(cell_id) %>%
    mutate(n_agents_same_cell = n()) %>%
    ungroup()
#------------------------------------------------------------
#F. Immediate reward from current ecological state C_t
#------------------------------------------------------------
  action_log<-action_log %>%
    left_join(
      land_grid %>%
        select(cell_id, rangeland_condition),
      by = "cell_id"
    )
  reward_components<-pmap(
    list(
      action_log$action,
      action_log$rangeland_condition,
      action_log$n_agents_same_cell,
      action_log$emergency_cost
    ),
    function(
        action,
        rangeland_condition,
        n_agents_same_cell,
        emergency_cost
    ) {
      compute_reward(
        action = action,
        rangeland_condition = rangeland_condition,
        n_agents_same_cell = n_agents_same_cell,
        emergency_cost = emergency_cost
      )
    }
  )
  action_log$reward<-map_dbl(reward_components,"net_reward")
  action_log$livestock_income<-map_dbl(reward_components,"livestock_income")
  action_log$pes_reward<-map_dbl(reward_components,"pes_reward")
  action_log$stock_survival<-map_dbl(reward_components,"stock_survival")
  action_log$expected_surviving_herd<-map_dbl(reward_components,"expected_surviving_herd")
  action_log$expected_stock_loss<-map_dbl(reward_components,"expected_stock_loss")
  action_log$interaction_multiplier<-map_dbl(reward_components,"interaction_multiplier")
#------------------------------------------------------------
#G. Rangeland evolves from t -> t+1

#A_t is the previous stocking category from the perspective of s_(t+1).
#W_(t+1) is the current weather that helps generate C_(t+1).
#Every cell is updated, including cells with ecological_footprint = "none".
  #------------------------------------------------------------
  landscape_update<-update_landscape(
    land_grid = land_grid,
    action_log = action_log,
    transition_weather = next_weather,
    step = step
  )
  next_land_grid<-landscape_update$land_grid
#------------------------------------------------------------
#H. Construct agent history/location for s_(t+1)
#------------------------------------------------------------
  next_agents<-agents
  for (id in names(next_agents)) {
    next_agents[[id]]$previous_action<-chosen_actions[[id]]
    next_agents[[id]]$herd_size<-unname(
      action_herd_size[[chosen_actions[[id]]]]
    )
    next_location<-resolve_location(
      agent_id = id,
      step = step + 1L,
      weather = next_weather,
      land_grid = next_land_grid
    )
    next_agents[[id]]$current_cell<-next_location$cell_id
  }
  next_shared_info<-share_agent_information(
    next_agents,
    next_land_grid
  )
  next_state_keys<-setNames(
    map_chr(
      next_agents,
      observe_agent_state,
      land_grid = next_land_grid,
      weather = next_weather,
      shared_info = next_shared_info
    ),
    names(next_agents)
  )
#------------------------------------------------------------
#I. Q-learning update
#------------------------------------------------------------
  q_diagnostics<-vector("list", length(agents))
  names(q_diagnostics)<-names(agents)
  for (id in names(agents)) {
    reward_i<-action_log$reward[
      action_log$agent_id == id
    ]
    update<-update_q(
      agent = agents[[id]],
      state_key = state_keys[[id]],
      action = chosen_actions[[id]],
      reward = reward_i,
      next_state_key = next_state_keys[[id]]
    )
    updated_agent<-update$agent
    updated_agent$previous_action<-chosen_actions[[id]]
    updated_agent$herd_size<-unname(
      action_herd_size[[chosen_actions[[id]]]]
    )
    updated_agent$current_cell<-next_agents[[id]]$current_cell
    agents[[id]]<-updated_agent
    q_diagnostics[[id]]<-tibble(
      agent_id = id,
      state_key = state_keys[[id]],
      previous_action = action_log$previous_action[
        action_log$agent_id == id
      ],
      action = chosen_actions[[id]],
      reward = reward_i,
      next_state_key = next_state_keys[[id]],
      q_old = update$q_old,
      q_new = update$q_new,
      td_target = update$td_target,
      td_error = update$td_error,
      next_feasible_actions = paste(
        update$next_feasible_actions,
        collapse = " | "
      )
    )
  }
  q_log<-bind_rows(q_diagnostics) %>%
    mutate(step = step)
  action_log<-action_log %>%
    mutate(
      step = step,
      weather = weather,
      next_weather = next_weather
    )
  list(
    agents = agents,
    land_grid = next_land_grid,
    action_log = action_log,
    cell_ledger = landscape_update$ledger,
    q_log = q_log,
    shared_info = shared_info,
    next_shared_info = next_shared_info
  )
}

19 . Full MARL simulation

The weather sequence has length steps + 1 so that \(s'\) includes the next observed weather condition when Q is updated.

simulate_marl<-function(
    steps = 500L,
    seed = NULL,
    initial_agents = NULL,
    initial_land_grid = NULL
) {
  if (!is.null(seed)) {
    set.seed(seed)
  }
  agents<-if (is.null(initial_agents)) {
    setNames(
      lapply(
        names(initial_actions),
        function(id) {
          create_agent(
            agent_id = id,
            initial_action = initial_actions[[id]]
          )
        }
      ),
      names(initial_actions)
    )
  } else {
    initial_agents
  }
  land_grid<-if (is.null(initial_land_grid)) {
    make_land_grid()
  } else {
    initial_land_grid
  }
#W_t is observed in s_t. W_(t+1) is used with A_t to generate C_(t+1).
  weather_sequence<-draw_weather(steps + 1L)
  action_logs<-vector("list", steps)
  cell_ledgers<-vector("list", steps)
  q_logs<-vector("list", steps)
  shared_ledgers<-vector("list", steps)
  landscape_counts<-vector("list", steps)
  for (t in seq_len(steps)) {
    result<-marl_step(
      agents = agents,
      land_grid = land_grid,
      weather = weather_sequence[t],
      next_weather = weather_sequence[t + 1L],
      step = t
    )
    agents<-result$agents
    land_grid<-result$land_grid
    action_logs[[t]]<-result$action_log
    cell_ledgers[[t]]<-result$cell_ledger
    q_logs[[t]]<-result$q_log
    shared_ledgers[[t]]<-result$shared_info %>%
      mutate(step = t)
    condition_count<-table(
      factor(
        land_grid$rangeland_condition,
        levels = r_states
      )
    )
    landscape_counts[[t]]<-tibble(
      step = t,
      intact = unname(condition_count["intact"]),
      good = unname(condition_count["good"]),
      poor = unname(condition_count["poor"]),
      degraded = unname(condition_count["degraded"]),
      pct_degraded =
        100 * unname(condition_count["degraded"]) / nrow(land_grid),
      mean_ema = mean(land_grid$condition_ema_score),
      max_ema = max(land_grid$condition_ema_score)
    )
  }
  list(
    agents = agents,
    land_grid = land_grid,
    weather = weather_sequence,
    action_log = bind_rows(action_logs),
    cell_ledger = bind_rows(cell_ledgers),
    q_log = bind_rows(q_logs),
    shared_ledger = bind_rows(shared_ledgers),
    landscape_counts = bind_rows(landscape_counts)
  )
}

20 . Train one MARL system

This is one continuous learning trajectory. Q tables carry forward throughout the run.

results_train<-simulate_marl(
  steps = 500,
  seed = 128
)

21 . Core diagnostics

21.1 .1 Landscape degradation

p_degraded_pct<-results_train$landscape_counts %>%
  ggplot(aes(step, pct_degraded)) +
  geom_line() +
  labs(
    title = "Percentage of degraded cells",
    x = "Step",
    y = "Degraded cells (%)"
  ) +
  theme_minimal()
p_degraded_pct

21.2 .2 Action mix through time

results_train$action_log %>%
  mutate(
    period = cut(
      step,
      breaks = seq(0, 500, by = 50),
      include.lowest = TRUE
    )
  ) %>%
  count(period, action) %>%
  group_by(period) %>%
  mutate(prop = n / sum(n)) %>%
  ungroup() %>%
  ggplot(aes(period, prop, fill = action)) +
  geom_col() +
  labs(
    title = "Learned stocking-action mix through time",
    x = "Training period",
    y = "Proportion of choices",
    fill = "Action"
  ) +
  theme_minimal() +
  theme(
    axis.text.x = element_text(
      angle = 45,
      hjust = 1
    )
  )

21.3 .3 Emergency movements

results_train$action_log %>%
  count(step, movement_mode) %>%
  ggplot(aes(step, n, fill = movement_mode)) +
  geom_col() +
  labs(
    title = "Emergency and designated-pasture movement",
    x = "Step",
    y = "Number of agents",
    fill = "Movement mode"
  ) +
  theme_minimal()

21.4 .4 Q-learning changes

results_train$q_log %>%
  mutate(
    period = ceiling(step / 25)
  ) %>%
  group_by(
    period,
    agent_id,
    action
  ) %>%
  summarise(
    mean_abs_td_error =
      mean(abs(td_error)),
    .groups = "drop"
  ) %>%
  ggplot(
    aes(
      period,
      mean_abs_td_error,
      colour = action
    )
  ) +
  geom_line() +
  facet_wrap(~ agent_id) +
  labs(
    title = "Temporal-difference error through learning",
    x = "25-step period",
    y = "Mean absolute TD error",
    colour = "Action"
  ) +
  theme_minimal()

21.5 .5 Cell-level ecological ledger

This shows more details of the ecological dynamics, an insight that illustrates what is going on inside the aggregated degradation percentage.

condition_numeric<-c(
  intact = 0,
  good = 1,
  poor = 2,
  degraded = 3
)
results_train$cell_ledger %>%
  mutate(
    condition_value =
      condition_numeric[condition_after]
  ) %>%
  ggplot(
    aes(
      step,
      factor(cell_id),
      fill = condition_value
    )
  ) +
  geom_tile() +
  scale_fill_viridis_c(
    breaks = 0:3,
    labels = r_states
  ) +
  labs(
    title = "Cell-level ecological ledger",
    x = "Step",
    y = "Cell ID",
    fill = "Condition"
  ) +
  theme_minimal()

22 . Independent replicate workflow

We need a continuous training trajectory with independent simulation replicates.This function creates independent runs by resetting both:

run_independent_replicates<-function(
    n_runs = 20L,
    steps = 100L,
    seed_start = 1000L
) {
  map_dfr(
    seq_len(n_runs),
    function(run_id) {
      res<-simulate_marl(
        steps = steps,
        seed = seed_start + run_id
      )
      res$landscape_counts %>%
        mutate(run = run_id)
    }
  )
}
replicate_degradation<-run_independent_replicates(
  n_runs = 20,
  steps = 100,
  seed_start = 2000
)

22.1 Complexity backlog — deliberately excluded for now

The following components are intentionally not implemented in this file so that they can be reintroduced one at a time later:

  1. Direct empirical estimation of current-condition-specific recovery; current non-intact rows use the documented interpolation and recovery constraints in Step 5.
  2. Agent action bias / \(\eta\).
  3. Agent reward memory, spatial memory, or learned location preferences.
  4. Historical emergency-movement propensity and stochastic mobility eligibility.
  5. Heterogeneous \(\alpha\), \(\gamma\), and \(\tau\) between agents.
  6. Outcome-based PES.
  7. Explicit social networks beyond shared location/herd-size information.
  8. Grazing effects on cells crossed during emergency travel.
  9. Joint-action Q tables, centralized critics, or deep MARL.
  10. Full landscape rollouts / Monte Carlo tree search.
  11. Seasonal pasture polygons larger than two designated cells.
  12. Endogenous herd demography and persistent carry-forward of livestock losses. Current-period stock survival is implemented, but surviving herd size is not yet a state variable.