1 Teaching goal

This RPubs article is written for students who may know physics or astronomy, but may not yet know operations research or game theory. The aim is to build the method slowly:

  1. A gravitational-wave sky map is a probability distribution on the sky.
  2. A telescope observation is an action that covers part of that probability distribution.
  3. Scheduling becomes an operations-research problem: choose actions under budgets.
  4. Several telescopes acting independently form a game.
  5. A carefully chosen local utility turns the game into an exact potential game.
  6. A post-campaign Shapley value gives an interpretable contribution score.
  7. The same logic explains the paper-level public-event examples such as GW170817 and GW190814.

The code is intentionally small and transparent. It is not meant to replace professional astronomical schedulers; it is a conceptual bridge from mathematical theory to multimessenger-astronomy decision making.

2 Physical motivation: why this is not just astronomy plotting

When a gravitational-wave alert arrives, the source position is uncertain. The event is not given as one exact point in the sky. Instead, the gravitational-wave pipelines release a sky-localization map: a grid of sky pixels, each with posterior probability.

A telescope network must answer:

Which telescope should observe which sky tile, at which time, and how should we avoid wasting scarce telescope time on redundant observations?

This is a resource-allocation problem under uncertainty. That is exactly the language of operations research.

If each telescope is autonomous, then another question appears:

How can local telescope decisions align with global scientific value?

That is the language of game theory.

3 Minimum mathematical language

3.1 Sky map

Let

\[ \mathcal H = \{1,\ldots,n\} \]

be a set of sky pixels. Pixel \(h\) has posterior probability

\[ \pi_h \geq 0, \qquad \sum_{h=1}^n \pi_h = 1. \]

Interpretation: before observing, \(\pi_h\) is the probability that the source lies in pixel \(h\).

3.2 Telescope action

An action \(a\) means:

\[ a = \text{``telescope } k \text{ observes tile } j \text{ at time slot } \tau \text{''}. \]

Let

\[ q_{ah} \in [0,1] \]

be the conditional probability that action \(a\) detects the counterpart if the source is in pixel \(h\). A simple model is

\[ q_{ah} = \mathbf 1\{h \in C_a\} \times \text{visibility quality} \times \text{depth quality} \times \text{time-decay factor}. \]

Here \(C_a\) is the footprint of the tile observed by action \(a\).

3.3 Posterior-weighted detection value

For a schedule \(S\), meaning a set of actions, the probability of at least one detection if the true source is in pixel \(h\) is

\[ 1 - \prod_{a\in S}(1-q_{ah}). \]

Averaging this over the sky posterior gives

\[ F(S) = \sum_{h\in\mathcal H} \pi_h \left[ 1-\prod_{a\in S}(1-q_{ah}) \right]. \]

This is the core objective.

If action \(a\) costs \(c_a\), define

\[ C(S)=\sum_{a\in S}c_a, \qquad \Phi(S)=F(S)-C(S). \]

Here \(F(S)\) is expected detection value, while \(\Phi(S)\) is net welfare.

4 A toy gravitational-wave sky map

We now make a small two-dimensional toy sky. This is not a true HEALPix sphere. It is a teaching analogue: a sky patch with posterior probability concentrated in two regions.

n_pix <- 450

sky <- data.frame(
  pixel = seq_len(n_pix),
  x = runif(n_pix, -1, 1),
  y = runif(n_pix, -1, 1)
)

gaussian_bump <- function(x, y, mx, my, sx, sy, weight = 1) {
  weight * exp(-0.5 * (((x - mx) / sx)^2 + ((y - my) / sy)^2))
}

raw_prob <-
  gaussian_bump(sky$x, sky$y, -0.35,  0.25, 0.22, 0.16, 1.00) +
  gaussian_bump(sky$x, sky$y,  0.42, -0.18, 0.30, 0.20, 0.70) +
  0.02

sky$pi <- raw_prob / sum(raw_prob)

ggplot(sky, aes(x, y)) +
  geom_point(aes(size = pi, colour = pi), alpha = 0.75) +
  coord_equal() +
  labs(
    title = "Toy gravitational-wave sky map",
    subtitle = "Each point is a sky pixel; point size and colour represent posterior probability.",
    x = "Toy sky coordinate x",
    y = "Toy sky coordinate y",
    colour = "Posterior",
    size = "Posterior"
  )

The brightest pixels represent the most probable sky regions. A good schedule should cover these regions quickly, but it should not blindly duplicate coverage if another telescope is already observing the same high-probability part.

5 Candidate tiles and telescope network

We create 12 possible sky tiles. Each tile has a centre and a circular footprint in this toy coordinate system.

tile_grid <- expand.grid(
  tx = seq(-0.75, 0.75, length.out = 4),
  ty = seq(-0.60, 0.60, length.out = 3)
)

tiles <- data.frame(
  tile = seq_len(nrow(tile_grid)),
  tx = tile_grid$tx,
  ty = tile_grid$ty,
  radius = 0.38
)

telescopes <- data.frame(
  telescope = c(
    "ZTF-like",
    "DECam-like",
    "MeerLICHT-like",
    "SkyMapper-like",
    "BlackGEM-like"
  ),
  hemisphere = c("north", "south", "south", "south", "south"),
  fov_factor = c(1.15, 1.00, 0.72, 0.85, 0.78),
  depth_factor = c(0.58, 0.92, 0.82, 0.68, 0.80),
  cost_per_action = c(0.0019, 0.0022, 0.0018, 0.0017, 0.0019),
  budget = c(3, 3, 3, 2, 3)
)

slots <- c(1, 2, 3, 4)

kable(telescopes, caption = "Toy heterogeneous telescope network.")
Toy heterogeneous telescope network.
telescope hemisphere fov_factor depth_factor cost_per_action budget
ZTF-like north 1.15 0.58 0.0019 3
DECam-like south 1.00 0.92 0.0022 3
MeerLICHT-like south 0.72 0.82 0.0018 3
SkyMapper-like south 0.85 0.68 0.0017 2
BlackGEM-like south 0.78 0.80 0.0019 3

The names are chosen to be astronomy-flavoured, not exact engineering models. The important feature is heterogeneity: different facilities have different depth, cost, and visibility.

ggplot() +
  geom_point(data = sky, aes(x, y, size = pi), alpha = 0.25) +
  geom_point(data = tiles, aes(tx, ty), size = 3, shape = 4, stroke = 1.2) +
  geom_text(data = tiles, aes(tx, ty, label = tile), nudge_y = 0.06) +
  coord_equal() +
  labs(
    title = "Candidate observation tiles",
    subtitle = "Crosses indicate tile centres. Each tile covers a circular neighbourhood.",
    x = "Toy sky coordinate x",
    y = "Toy sky coordinate y",
    size = "Posterior"
  )

6 Building telescope actions

An action combines telescope, tile, and time slot. We also build a matrix \(Q\), where row \(a\) and column \(h\) stores \(q_{ah}\).

make_actions <- function(sky, tiles, telescopes, slots) {
  action_list <- list()
  Q_list <- list()
  id <- 1

  for (k in seq_len(nrow(telescopes))) {
    tel <- telescopes[k, ]

    for (j in seq_len(nrow(tiles))) {
      tile <- tiles[j, ]

      distance <- sqrt((sky$x - tile$tx)^2 + (sky$y - tile$ty)^2)
      in_footprint <- as.numeric(distance <= tile$radius * tel$fov_factor)

      # Toy visibility:
      # northern facility prefers positive y; southern facilities prefer negative y,
      # but neither is completely blind in this simplified example.
      if (tel$hemisphere == "north") {
        visibility <- 0.45 + 0.55 * pmin(pmax((sky$y + 1) / 2, 0), 1)
      } else {
        visibility <- 0.45 + 0.55 * pmin(pmax((1 - sky$y) / 2, 0), 1)
      }

      for (s in slots) {
        # Counterpart fading: later slots have smaller detection probability.
        time_decay <- exp(-0.09 * (s - 1))

        q <- in_footprint * visibility * tel$depth_factor * time_decay

        # Keep q in [0, 1].
        q <- pmin(pmax(q, 0), 1)

        action_id <- paste0("a", id)
        action_list[[id]] <- data.frame(
          action = action_id,
          telescope = tel$telescope,
          tile = tile$tile,
          slot = s,
          cost = tel$cost_per_action,
          local_score = sum(sky$pi * q),
          stringsAsFactors = FALSE
        )
        Q_list[[id]] <- q
        id <- id + 1
      }
    }
  }

  actions <- do.call(rbind, action_list)
  Q <- do.call(rbind, Q_list)
  rownames(Q) <- actions$action

  list(actions = actions, Q = Q)
}

built <- make_actions(sky, tiles, telescopes, slots)
actions <- built$actions
Q <- built$Q

actions %>%
  arrange(desc(local_score)) %>%
  head(12) %>%
  mutate(across(where(is.numeric), ~ round(.x, 4))) %>%
  kable(caption = "Top 12 individual actions by posterior-weighted local detection score.")
Top 12 individual actions by posterior-weighted local detection score.
action telescope tile slot cost local_score
a69 DECam-like 6 1 0.0022 0.2151
a73 DECam-like 7 1 0.0022 0.2077
a70 DECam-like 6 2 0.0022 0.1966
a21 ZTF-like 6 1 0.0019 0.1959
a74 DECam-like 7 2 0.0022 0.1898
a71 DECam-like 6 3 0.0022 0.1797
a22 ZTF-like 6 2 0.0019 0.1790
a75 DECam-like 7 3 0.0022 0.1735
a72 DECam-like 6 4 0.0022 0.1642
a23 ZTF-like 6 3 0.0019 0.1636
a76 DECam-like 7 4 0.0022 0.1586
a25 ZTF-like 7 1 0.0019 0.1508

7 The detection and welfare functions in R

We now implement

\[ F(S) = \sum_h \pi_h \left[ 1-\prod_{a\in S}(1-q_{ah}) \right] \]

and

\[ \Phi(S)=F(S)-C(S). \]

F_value <- function(action_set, Q, pi_vec) {
  action_set <- unique(action_set)
  if (length(action_set) == 0) return(0)
  Qs <- Q[action_set, , drop = FALSE]
  detect_by_pixel <- 1 - apply(1 - Qs, 2, prod)
  sum(pi_vec * detect_by_pixel)
}

cost_value <- function(action_set, actions) {
  action_set <- unique(action_set)
  if (length(action_set) == 0) return(0)
  sum(actions$cost[match(action_set, actions$action)])
}

Phi_value <- function(action_set, Q, pi_vec, actions) {
  F_value(action_set, Q, pi_vec) - cost_value(action_set, actions)
}

example_set <- actions %>%
  arrange(desc(local_score)) %>%
  head(5) %>%
  pull(action)

data.frame(
  schedule_size = length(example_set),
  expected_detection = F_value(example_set, Q, sky$pi),
  cost = cost_value(example_set, actions),
  welfare = Phi_value(example_set, Q, sky$pi, actions)
) %>%
  mutate(across(where(is.numeric), ~ round(.x, 5))) %>%
  kable(caption = "Detection value and welfare for a simple top-five schedule.")
Detection value and welfare for a simple top-five schedule.
schedule_size expected_detection cost welfare
5 0.55353 0.0107 0.54283

8 Operations research idea: greedy scheduling

A central planner can greedily add actions that increase welfare the most, subject to each telescope’s budget and the restriction that a telescope cannot observe two tiles in the same slot.

This is an operations-research heuristic. It is simple, interpretable, and connected to the theory of submodular maximization.

is_feasible_global <- function(candidate_set, actions, telescopes) {
  if (length(candidate_set) == 0) return(TRUE)
  sel <- actions[match(candidate_set, actions$action), ]

  # No telescope uses the same time slot twice.
  duplicated_pair <- any(duplicated(paste(sel$telescope, sel$slot)))
  if (duplicated_pair) return(FALSE)

  # Telescope budget.
  tab <- table(sel$telescope)
  for (nm in names(tab)) {
    b <- telescopes$budget[match(nm, telescopes$telescope)]
    if (tab[[nm]] > b) return(FALSE)
  }
  TRUE
}

centralized_greedy <- function(actions, Q, pi_vec, telescopes, max_total = 12) {
  chosen <- character(0)
  history <- data.frame(step = integer(), added = character(), welfare = numeric())

  for (step in seq_len(max_total)) {
    base <- Phi_value(chosen, Q, pi_vec, actions)
    best_gain <- 0
    best_action <- NA_character_

    for (a in actions$action) {
      cand <- unique(c(chosen, a))
      if (length(cand) == length(chosen)) next
      if (!is_feasible_global(cand, actions, telescopes)) next

      gain <- Phi_value(cand, Q, pi_vec, actions) - base
      if (gain > best_gain) {
        best_gain <- gain
        best_action <- a
      }
    }

    if (is.na(best_action) || best_gain <= 1e-12) break

    chosen <- unique(c(chosen, best_action))
    history <- rbind(
      history,
      data.frame(
        step = step,
        added = best_action,
        welfare = Phi_value(chosen, Q, pi_vec, actions)
      )
    )
  }

  list(schedule = chosen, history = history)
}

greedy_out <- centralized_greedy(actions, Q, sky$pi, telescopes, max_total = 12)

greedy_table <- actions[match(greedy_out$schedule, actions$action), ] %>%
  mutate(across(where(is.numeric), ~ round(.x, 4)))

kable(greedy_table, caption = "Centralized greedy toy schedule.")
Centralized greedy toy schedule.
action telescope tile slot cost local_score
69 a69 DECam-like 6 1 0.0022 0.2151
74 a74 DECam-like 7 2 0.0022 0.1898
21 a21 ZTF-like 6 1 0.0019 0.1959
59 a59 DECam-like 3 3 0.0022 0.0815
38 a38 ZTF-like 10 2 0.0019 0.1050
173 a173 SkyMapper-like 8 1 0.0017 0.0682
19 a19 ZTF-like 5 3 0.0019 0.0808
217 a217 BlackGEM-like 7 1 0.0019 0.1147
158 a158 SkyMapper-like 4 2 0.0017 0.0280
133 a133 MeerLICHT-like 10 1 0.0018 0.0419
222 a222 BlackGEM-like 8 2 0.0019 0.0497
118 a118 MeerLICHT-like 6 2 0.0018 0.1056 
ggplot(greedy_out$history, aes(step, welfare)) +
  geom_line(linewidth = 1) +
  geom_point(size = 2.5) +
  labs(
    title = "Centralized greedy improvement path",
    subtitle = "The planner repeatedly adds the feasible action with the largest welfare gain.",
    x = "Greedy step",
    y = "Net welfare"
  )

9 Submodularity: diminishing returns

The function \(F(S)\) is submodular. Intuitively:

If a sky region has already been covered by several observations, one more observation of the same region adds less new information than it would have added at the beginning.

Formally, for \(S \subseteq T\) and \(b \notin T\),

\[ F(S\cup\{b\})-F(S) \geq F(T\cup\{b\})-F(T). \]

The next small experiment demonstrates this effect.

# Choose one candidate action b and compare its gain when added to a small
# schedule versus a larger schedule.
b <- actions %>%
  arrange(desc(local_score)) %>%
  slice(8) %>%
  pull(action)

S_small <- greedy_out$schedule[1:2]
T_large <- greedy_out$schedule[1:8]

gain_small <- F_value(c(S_small, b), Q, sky$pi) - F_value(S_small, Q, sky$pi)
gain_large <- F_value(c(T_large, b), Q, sky$pi) - F_value(T_large, Q, sky$pi)

data.frame(
  candidate_action = b,
  marginal_gain_after_small_schedule = gain_small,
  marginal_gain_after_large_schedule = gain_large
) %>%
  mutate(across(where(is.numeric), ~ round(.x, 6))) %>%
  kable(caption = "A numerical illustration of diminishing returns.")
A numerical illustration of diminishing returns.
candidate_action marginal_gain_after_small_schedule marginal_gain_after_large_schedule
a75 0.053571 0.018928

10 Game theory basics

Now suppose each telescope is a player.

  • Player \(k\): telescope \(k\).
  • Strategy \(S_k\): the set of actions chosen by telescope \(k\).
  • Strategy profile:

\[ S=(S_1,\ldots,S_K). \]

  • Aggregate schedule:

\[ \bar S = \bigcup_{k=1}^K S_k. \]

If every telescope greedily maximizes only its own local posterior coverage, it may duplicate other facilities. A better idea is to reward a telescope for its marginal contribution to the global welfare.

Define

\[ u_k(S_k,S_{-k}) = \Phi\left(S_k\cup \bar S_{-k}\right) - \Phi\left(\bar S_{-k}\right). \]

This means:

Telescope \(k\)’s utility is exactly the extra welfare added by its chosen actions beyond what the other telescopes already provide.

This discourages redundant observations.

11 Exact potential game

Define the potential function

\[ \Psi(S_1,\ldots,S_K)=\Phi\left(\bigcup_{k=1}^K S_k\right). \]

If telescope \(k\) changes from \(S_k\) to \(S'_k\), then

\[ u_k(S'_k,S_{-k})-u_k(S_k,S_{-k}) = \Psi(S'_k,S_{-k})-\Psi(S_k,S_{-k}). \]

Therefore, the game is an exact potential game.

This is powerful because in a finite game, repeated strict better responses must eventually stop. The stopping point is a pure-strategy Nash equilibrium.

12 Fast local strategy library

The exact local best-response problem can become combinatorial. In the manuscript-level problem, the best-response subproblem is an operations-research problem that may be solved by local search, dynamic programming, integer programming, or greedy approximation.

For RPubs teaching, we use a fast approximation:

  1. Score every local action cheaply.
  2. Keep only the best few actions per time slot.
  3. Enumerate feasible combinations only among the pruned actions.

This preserves the concept while avoiding long knitting times.

all_local_strategies <- function(actions, telescope_name, budget, top_per_slot = 3) {
  tel_actions <- actions[actions$telescope == telescope_name, ]
  tel_actions$net_score <- tel_actions$local_score - tel_actions$cost

  # Candidate pruning: keep only the top few actions in each time slot.
  pruned <- do.call(
    rbind,
    lapply(split(tel_actions, tel_actions$slot), function(d) {
      d <- d[order(-d$net_score), ]
      d[seq_len(min(nrow(d), top_per_slot)), ]
    })
  )

  ids <- pruned$action
  max_b <- min(budget, length(ids), length(unique(pruned$slot)))

  strategies <- list(character(0))

  if (max_b == 0) return(strategies)

  for (b in seq_len(max_b)) {
    cmb <- combn(ids, b, simplify = FALSE)
    for (s in cmb) {
      sel <- actions[match(s, actions$action), ]
      if (!any(duplicated(sel$slot))) {
        strategies[[length(strategies) + 1]] <- s
      }
    }
  }
  strategies
}

strategy_library <- list()

for (i in seq_len(nrow(telescopes))) {
  nm <- telescopes$telescope[i]
  strategy_library[[nm]] <- all_local_strategies(
    actions = actions,
    telescope_name = nm,
    budget = telescopes$budget[i],
    top_per_slot = 3
  )
}

data.frame(
  telescope = names(strategy_library),
  number_of_candidate_strategies = vapply(strategy_library, length, integer(1))
) %>%
  kable(caption = "Pruned local strategy-library sizes used for fast RPubs knitting.")
Pruned local strategy-library sizes used for fast RPubs knitting.
telescope number_of_candidate_strategies
ZTF-like ZTF-like 175
DECam-like DECam-like 175
MeerLICHT-like MeerLICHT-like 175
SkyMapper-like SkyMapper-like 67
BlackGEM-like BlackGEM-like 175

13 Best-response dynamics

A best-response update means that one telescope changes its schedule to maximize its marginal contribution, keeping all other telescopes fixed.

union_profile <- function(profile) {
  unique(unlist(profile, use.names = FALSE))
}

best_response_for_tel <- function(tel_name, profile, strategy_library, Q, pi_vec, actions) {
  others <- profile
  others[[tel_name]] <- character(0)
  other_actions <- union_profile(others)

  best_strategy <- profile[[tel_name]]
  best_utility <- -Inf

  for (candidate in strategy_library[[tel_name]]) {
    utility <- Phi_value(c(other_actions, candidate), Q, pi_vec, actions) -
      Phi_value(other_actions, Q, pi_vec, actions)

    if (utility > best_utility + 1e-12) {
      best_utility <- utility
      best_strategy <- candidate
    }
  }

  list(strategy = best_strategy, utility = best_utility)
}

potential_game_scheduler <- function(telescopes, strategy_library, Q, pi_vec, actions,
                                     max_rounds = 8, tol = 1e-10) {
  tel_names <- telescopes$telescope
  profile <- setNames(vector("list", length(tel_names)), tel_names)
  for (nm in tel_names) profile[[nm]] <- character(0)

  history <- data.frame(
    round = integer(),
    telescope = character(),
    welfare = numeric(),
    changed = logical(),
    stringsAsFactors = FALSE
  )

  for (r in seq_len(max_rounds)) {
    any_changed <- FALSE

    for (nm in tel_names) {
      old <- profile[[nm]]
      br <- best_response_for_tel(nm, profile, strategy_library, Q, pi_vec, actions)
      new <- br$strategy

      changed <- !setequal(old, new)
      if (changed) any_changed <- TRUE
      profile[[nm]] <- new

      history <- rbind(
        history,
        data.frame(
          round = r,
          telescope = nm,
          welfare = Phi_value(union_profile(profile), Q, pi_vec, actions),
          changed = changed,
          stringsAsFactors = FALSE
        )
      )
    }

    if (!any_changed) break
  }

  list(profile = profile, schedule = union_profile(profile), history = history)
}

pg_out <- potential_game_scheduler(
  telescopes = telescopes,
  strategy_library = strategy_library,
  Q = Q,
  pi_vec = sky$pi,
  actions = actions,
  max_rounds = 8
)

pg_schedule_table <- actions[match(pg_out$schedule, actions$action), ] %>%
  arrange(telescope, slot) %>%
  mutate(across(where(is.numeric), ~ round(.x, 4)))

kable(pg_schedule_table, caption = "Potential-game toy schedule after best-response convergence.")
Potential-game toy schedule after best-response convergence.
action telescope tile slot cost local_score
a221 BlackGEM-like 8 1 0.0019 0.0544
a214 BlackGEM-like 6 2 0.0019 0.1166
a223 BlackGEM-like 8 3 0.0019 0.0454
a73 DECam-like 7 1 0.0022 0.2077
a70 DECam-like 6 2 0.0022 0.1966
a79 DECam-like 8 3 0.0022 0.0881
a121 MeerLICHT-like 7 1 0.0018 0.1049
a118 MeerLICHT-like 6 2 0.0018 0.1056
a119 MeerLICHT-like 6 3 0.0018 0.0965
a169 SkyMapper-like 7 1 0.0017 0.1212
a170 SkyMapper-like 7 2 0.0017 0.1108
a21 ZTF-like 6 1 0.0019 0.1959
a38 ZTF-like 10 2 0.0019 0.1050
a39 ZTF-like 10 3 0.0019 0.0959 
ggplot(pg_out$history, aes(seq_along(welfare), welfare)) +
  geom_line(linewidth = 1) +
  geom_point(aes(shape = changed), size = 2.4) +
  labs(
    title = "Best-response dynamics in the toy potential game",
    subtitle = "Each point is one telescope update. The potential is the net welfare of the aggregate schedule.",
    x = "Update index",
    y = "Potential / net welfare",
    shape = "Changed?"
  )

14 Comparing scheduling rules

We now compare four methods:

  1. Centralized greedy: a global planner greedily maximizes net welfare.
  2. Potential game: telescopes best-respond using marginal-contribution utilities.
  3. Independent greedy: each telescope chooses its own best actions ignoring overlap.
  4. Random feasible: a naive feasible baseline.
independent_greedy <- function(actions, Q, pi_vec, telescopes) {
  chosen <- character(0)

  for (i in seq_len(nrow(telescopes))) {
    nm <- telescopes$telescope[i]
    b <- telescopes$budget[i]

    tel_actions <- actions[actions$telescope == nm, ]
    tel_actions$net <- tel_actions$local_score - tel_actions$cost
    tel_actions <- tel_actions[order(-tel_actions$net), ]

    local <- character(0)
    used_slots <- integer(0)

    for (ii in seq_len(nrow(tel_actions))) {
      if (length(local) >= b) break
      if (!(tel_actions$slot[ii] %in% used_slots)) {
        local <- c(local, tel_actions$action[ii])
        used_slots <- c(used_slots, tel_actions$slot[ii])
      }
    }
    chosen <- c(chosen, local)
  }
  unique(chosen)
}

random_feasible <- function(actions, telescopes) {
  chosen <- character(0)
  for (i in seq_len(nrow(telescopes))) {
    nm <- telescopes$telescope[i]
    b <- telescopes$budget[i]
    tel_actions <- actions[actions$telescope == nm, ]
    tel_actions <- tel_actions[sample(seq_len(nrow(tel_actions))), ]

    local <- character(0)
    used_slots <- integer(0)
    for (ii in seq_len(nrow(tel_actions))) {
      if (length(local) >= b) break
      if (!(tel_actions$slot[ii] %in% used_slots)) {
        local <- c(local, tel_actions$action[ii])
        used_slots <- c(used_slots, tel_actions$slot[ii])
      }
    }
    chosen <- c(chosen, local)
  }
  unique(chosen)
}

schedules <- list(
  "centralized_greedy" = greedy_out$schedule,
  "potential_game" = pg_out$schedule,
  "independent_greedy" = independent_greedy(actions, Q, sky$pi, telescopes),
  "random_feasible" = random_feasible(actions, telescopes)
)

method_summary <- data.frame(
  method = names(schedules),
  n_actions = vapply(schedules, length, integer(1)),
  expected_detection = vapply(schedules, F_value, numeric(1), Q = Q, pi_vec = sky$pi),
  cost = vapply(schedules, cost_value, numeric(1), actions = actions),
  welfare = vapply(schedules, Phi_value, numeric(1), Q = Q, pi_vec = sky$pi, actions = actions),
  unique_tiles = vapply(schedules, function(s) length(unique(actions$tile[match(s, actions$action)])), integer(1)),
  telescopes_used = vapply(schedules, function(s) length(unique(actions$telescope[match(s, actions$action)])), integer(1))
) %>%
  arrange(desc(welfare))

method_summary %>%
  mutate(across(where(is.numeric), ~ round(.x, 5))) %>%
  kable(caption = "Toy benchmark comparison.")
Toy benchmark comparison.
method n_actions expected_detection cost welfare unique_tiles telescopes_used
centralized_greedy centralized_greedy 12 0.75324 0.0231 0.73014 7 5
potential_game potential_game 14 0.69443 0.0268 0.66763 4 5
independent_greedy independent_greedy 14 0.54888 0.0268 0.52208 2 5
random_feasible random_feasible 14 0.48288 0.0268 0.45608 9 5 
method_long <- method_summary %>%
  select(method, expected_detection, welfare) %>%
  pivot_longer(cols = c(expected_detection, welfare), names_to = "metric", values_to = "value")

ggplot(method_long, aes(reorder(method, value), value, fill = metric)) +
  geom_col(position = "dodge") +
  coord_flip() +
  labs(
    title = "Toy scheduler comparison",
    subtitle = "The potential-game schedule is a decentralized alternative to a centralized planner.",
    x = "Method",
    y = "Score",
    fill = "Metric"
  )

15 Shapley value: fair scientific credit

The schedule tells us what was observed. But in multi-observatory campaigns, we also want to know:

Which telescope made how much marginal contribution to the final scientific opportunity?

The Shapley value answers this by averaging marginal contribution over all possible arrival orders of the players.

For the toy teaching example, we compute a realized Shapley value: the action set selected by each telescope is fixed by the potential-game schedule, and we allocate the final welfare among telescopes.

all_permutations <- function(x) {
  if (length(x) == 1) return(list(x))
  out <- list()
  for (i in seq_along(x)) {
    rest <- x[-i]
    sub <- all_permutations(rest)
    for (s in sub) out[[length(out) + 1]] <- c(x[i], s)
  }
  out
}

realized_shapley <- function(profile, Q, pi_vec, actions) {
  players <- names(profile)
  perms <- all_permutations(players)
  credit <- setNames(rep(0, length(players)), players)

  for (perm in perms) {
    coalition_actions <- character(0)
    current_value <- Phi_value(coalition_actions, Q, pi_vec, actions)

    for (p in perm) {
      new_actions <- unique(c(coalition_actions, profile[[p]]))
      new_value <- Phi_value(new_actions, Q, pi_vec, actions)
      credit[p] <- credit[p] + (new_value - current_value)
      coalition_actions <- new_actions
      current_value <- new_value
    }
  }

  credit / length(perms)
}

toy_shapley <- realized_shapley(pg_out$profile, Q, sky$pi, actions)

shapley_table <- data.frame(
  telescope = names(toy_shapley),
  realized_shapley_value = as.numeric(toy_shapley),
  n_actions = vapply(pg_out$profile, length, integer(1))
) %>%
  arrange(desc(realized_shapley_value))

shapley_table %>%
  mutate(realized_shapley_value = round(realized_shapley_value, 5)) %>%
  kable(caption = "Realized Shapley credit for the toy potential-game schedule.")
Realized Shapley credit for the toy potential-game schedule.
telescope realized_shapley_value n_actions
DECam-like DECam-like 0.21409 3
ZTF-like ZTF-like 0.19720 3
MeerLICHT-like MeerLICHT-like 0.09935 3
BlackGEM-like BlackGEM-like 0.08066 3
SkyMapper-like SkyMapper-like 0.07632 2 
ggplot(shapley_table, aes(reorder(telescope, realized_shapley_value), realized_shapley_value)) +
  geom_col() +
  coord_flip() +
  labs(
    title = "Realized Shapley credit in the toy telescope network",
    subtitle = "Credit is assigned by average marginal contribution to final welfare.",
    x = "Telescope",
    y = "Realized Shapley value"
  )

16 Physics/astro toy interpretation

The toy model can be read physically as follows:

  • The bright probability islands represent possible source localizations.
  • Tile footprints are telescope pointings.
  • Visibility factors mimic observability from different hemispheres.
  • Depth factors mimic limiting magnitude.
  • Time decay mimics fading of an electromagnetic counterpart.
  • Cost mimics exposure time, overhead, or opportunity cost.
  • \(F(S)\) measures expected chance of detection.
  • \(\Phi(S)\) measures expected detection after subtracting operational burden.
  • The potential game makes each telescope care about the value it adds beyond others.
  • Shapley credit tells us whose observations were quantitatively useful.

This conceptual structure remains the same when the toy sky is replaced by a public gravitational-wave HEALPix sky map.

17 Bridge to public-event-style results

The following table reproduces the manuscript-style summary for two public gravitational-wave events. Here the data are embedded directly in the R Markdown so that the document is standalone.

real_metrics <- data.frame(
  event = c(
    "GW170817", "GW170817", "GW170817", "GW170817",
    "GW190814", "GW190814", "GW190814", "GW190814"
  ),
  method = c(
    "proposed_potential_game", "centralized_greedy", "independent_greedy", "random_feasible",
    "proposed_potential_game", "centralized_greedy", "independent_greedy", "random_feasible"
  ),
  n_actions = c(9, 9, 9, 9, 12, 12, 12, 12),
  expected_detection = c(0.121277, 0.120291, 0.120580, 0.094848,
                         0.654073, 0.662937, 0.641847, 0.624945),
  unique_covered_mass = c(0.708058, 0.722862, 0.722862, 0.620815,
                          0.934254, 0.921806, 0.874235, 0.935519),
  overlap_mass = c(0.448968, 0.278122, 0.340782, 0.215602,
                   0.861279, 0.859900, 0.861279, 0.864565),
  total_cost = c(0.01925, 0.01925, 0.01925, 0.01925,
                 0.02500, 0.02500, 0.02500, 0.02500),
  welfare = c(0.102027, 0.101041, 0.101330, 0.075598,
              0.629073, 0.637937, 0.616847, 0.599945),
  telescopes_used = c(3, 3, 3, 3, 4, 4, 4, 4),
  unique_tiles = c(5, 6, 6, 5, 7, 6, 6, 8)
)

real_metrics %>%
  mutate(across(where(is.numeric), ~ round(.x, 6))) %>%
  kable(caption = "Public-event-style method comparison used in the manuscript.")
Public-event-style method comparison used in the manuscript.
event method n_actions expected_detection unique_covered_mass overlap_mass total_cost welfare telescopes_used unique_tiles
GW170817 proposed_potential_game 9 0.121277 0.708058 0.448968 0.01925 0.102027 3 5
GW170817 centralized_greedy 9 0.120291 0.722862 0.278122 0.01925 0.101041 3 6
GW170817 independent_greedy 9 0.120580 0.722862 0.340782 0.01925 0.101330 3 6
GW170817 random_feasible 9 0.094848 0.620815 0.215602 0.01925 0.075598 3 5
GW190814 proposed_potential_game 12 0.654073 0.934254 0.861279 0.02500 0.629073 4 7
GW190814 centralized_greedy 12 0.662937 0.921806 0.859900 0.02500 0.637937 4 6
GW190814 independent_greedy 12 0.641847 0.874235 0.861279 0.02500 0.616847 4 6
GW190814 random_feasible 12 0.624945 0.935519 0.864565 0.02500 0.599945 4 8 
ggplot(real_metrics, aes(method, welfare, fill = method)) +
  geom_col(show.legend = FALSE) +
  facet_wrap(~ event, scales = "free_y") +
  coord_flip() +
  labs(
    title = "Welfare comparison for public-event-style results",
    subtitle = "The proposed potential-game scheduler is decentralized and remains competitive with centralized greedy.",
    x = "Method",
    y = "Net welfare"
  )

ggplot(real_metrics, aes(method, expected_detection, fill = method)) +
  geom_col(show.legend = FALSE) +
  facet_wrap(~ event, scales = "free_y") +
  coord_flip() +
  labs(
    title = "Expected detection comparison for public-event-style results",
    subtitle = "Expected detection is posterior-weighted and accounts for repeated probabilistic coverage.",
    x = "Method",
    y = "Expected detection"
  )

18 Proposed schedules in the public-event-style summary

proposed_schedule <- data.frame(
  event = c(
    "GW170817", "GW170817", "GW170817",
    "GW190814", "GW190814", "GW190814", "GW190814"
  ),
  telescope = c(
    "Network_B_DECamLike",
    "Network_C_MeerLICHTLike",
    "Network_E_BlackGEMLike",
    "Network_A_ZTFLike",
    "Network_B_DECamLike",
    "Network_C_MeerLICHTLike",
    "Network_E_BlackGEMLike"
  ),
  selected_slot_hours = c(
    "10.5, 11.0, 11.5",
    "4.5, 5.0, 5.5",
    "10.5, 11.0, 11.5",
    "13.0, 13.5, 14.0",
    "10.5, 11.0, 11.5",
    "4.0, 4.5, 5.0",
    "10.0, 10.5, 11.0"
  ),
  selected_tiles = c(
    "3, 0, 2",
    "1, 0, 4",
    "3, 0, 2",
    "8, 1, 7",
    "0, 1, 3",
    "1, 0, 2",
    "4, 0, 1"
  ),
  mean_q = c(
    0.075302, 0.123426, 0.071225,
    0.034507, 0.380429, 0.556140, 0.364986
  )
)

proposed_schedule %>%
  mutate(mean_q = round(mean_q, 6)) %>%
  kable(caption = "Condensed proposed potential-game schedule for the two public-event-style examples.")
Condensed proposed potential-game schedule for the two public-event-style examples.
event telescope selected_slot_hours selected_tiles mean_q
GW170817 Network_B_DECamLike 10.5, 11.0, 11.5 3, 0, 2 0.075302
GW170817 Network_C_MeerLICHTLike 4.5, 5.0, 5.5 1, 0, 4 0.123426
GW170817 Network_E_BlackGEMLike 10.5, 11.0, 11.5 3, 0, 2 0.071225
GW190814 Network_A_ZTFLike 13.0, 13.5, 14.0 8, 1, 7 0.034507
GW190814 Network_B_DECamLike 10.5, 11.0, 11.5 0, 1, 3 0.380429
GW190814 Network_C_MeerLICHTLike 4.0, 4.5, 5.0 1, 0, 2 0.556140
GW190814 Network_E_BlackGEMLike 10.0, 10.5, 11.0 4, 0, 1 0.364986 
ggplot(proposed_schedule, aes(telescope, mean_q, fill = event)) +
  geom_col(position = "dodge") +
  coord_flip() +
  labs(
    title = "Mean action quality in proposed schedules",
    subtitle = "The schedule structure changes from event to event because visibility and posterior geometry change.",
    x = "Telescope",
    y = "Mean q",
    fill = "Event"
  )

19 Shapley credit in the public-event-style summary

real_shapley <- data.frame(
  event = c(
    "GW170817", "GW170817", "GW170817", "GW170817", "GW170817",
    "GW190814", "GW190814", "GW190814", "GW190814", "GW190814"
  ),
  telescope = c(
    "Network_C_MeerLICHTLike",
    "Network_E_BlackGEMLike",
    "Network_B_DECamLike",
    "Network_A_ZTFLike",
    "Network_D_SkyMapperLike",
    "Network_C_MeerLICHTLike",
    "Network_B_DECamLike",
    "Network_E_BlackGEMLike",
    "Network_A_ZTFLike",
    "Network_D_SkyMapperLike"
  ),
  realized_shapley_value = c(
    0.050335, 0.026128, 0.025564, 0, 0,
    0.270491, 0.160863, 0.157275, 0.040445, 0
  ),
  n_actions = c(3, 3, 3, 0, 0, 3, 3, 3, 3, 0)
)

real_shapley %>%
  mutate(realized_shapley_value = round(realized_shapley_value, 6)) %>%
  kable(caption = "Realized Shapley-value decomposition for the public-event-style examples.")
Realized Shapley-value decomposition for the public-event-style examples.
event telescope realized_shapley_value n_actions
GW170817 Network_C_MeerLICHTLike 0.050335 3
GW170817 Network_E_BlackGEMLike 0.026128 3
GW170817 Network_B_DECamLike 0.025564 3
GW170817 Network_A_ZTFLike 0.000000 0
GW170817 Network_D_SkyMapperLike 0.000000 0
GW190814 Network_C_MeerLICHTLike 0.270491 3
GW190814 Network_B_DECamLike 0.160863 3
GW190814 Network_E_BlackGEMLike 0.157275 3
GW190814 Network_A_ZTFLike 0.040445 3
GW190814 Network_D_SkyMapperLike 0.000000 0 
ggplot(real_shapley, aes(reorder(telescope, realized_shapley_value), realized_shapley_value, fill = event)) +
  geom_col(position = "dodge") +
  coord_flip() +
  labs(
    title = "Realized Shapley credit for public-event-style schedules",
    subtitle = "Credit reflects marginal scientific contribution, not merely nominal participation.",
    x = "Telescope",
    y = "Realized Shapley value",
    fill = "Event"
  )

20 What students should notice

20.1 Conceptual message

The main conceptual movement is:

\[ \text{sky map} \longrightarrow \text{probabilistic objective} \longrightarrow \text{constrained scheduling} \longrightarrow \text{decentralized game} \longrightarrow \text{credit allocation}. \]

This is not merely a computational pipeline. It is a modelling philosophy.

20.2 Mathematical message

The objective

\[ F(S) = \sum_h \pi_h \left[ 1-\prod_{a\in S}(1-q_{ah}) \right] \]

has diminishing returns. That is why greedy operations-research methods are sensible.

The utility

\[ u_k(S_k,S_{-k}) = \Phi(S_k\cup \bar S_{-k})-\Phi(\bar S_{-k}) \]

turns local telescope updates into an exact potential game. That is why decentralized best-response dynamics have a clean equilibrium interpretation.

The Shapley value then supplies a principled contribution score.

20.3 Scientific message

The public-event-style results should be interpreted carefully:

  • The decentralized potential-game scheduler is not guaranteed to beat a strong centralized scheduler.
  • Its value is that it gives a mathematically interpretable decentralized mechanism.
  • It improves over naive feasible scheduling and can improve over uncoordinated selfish behaviour.
  • It also gives a transparent credit-allocation layer.

21 Classroom exercises

  1. Increase the number of sky pixels. Does the algorithm still knit quickly?
  2. Change top_per_slot from 3 to 5 in the local strategy library. How does speed change?
  3. Increase the time-decay rate. Which telescopes become more important?
  4. Increase the cost of the deepest telescope. Does the schedule diversify?
  5. Create a third posterior probability island in the sky map. Does the Shapley allocation become more balanced?
  6. Replace the toy circular tiles with rectangular fields of view.
  7. Implement weather uncertainty by randomly reducing \(q_{ah}\).
  8. Compare best-response dynamics from an empty profile versus a centralized warm start.
  9. Interpret the difference between expected_detection and unique_covered_mass.
  10. Explain why overlap mass is not always bad when individual detection probabilities are less than one.

22 A short glossary

Action: a telescope-tile-time choice.

Budget: a constraint on the number of observations, exposure time, or cost.

Posterior sky probability: the probability assigned to each sky pixel after gravitational-wave localization.

Coverage: the amount of posterior probability mass reached by a schedule.

Submodularity: a diminishing-returns property of set functions.

Operations research: mathematical optimization of decisions under constraints.

Game theory: mathematical study of interacting decision makers.

Potential game: a game where unilateral utility improvements match improvements in a global potential function.

Nash equilibrium: a profile where no single player can improve by changing alone.

Shapley value: a cooperative-game-theoretic allocation of total value based on average marginal contribution.

23 Conclusion

This standalone tutorial has built the paper’s theory from the ground up. We began with a probability map on the sky, defined telescope actions, constructed a posterior-weighted detection objective, introduced cost-adjusted welfare, and showed how scheduling becomes an operations-research problem. We then moved from centralized greedy scheduling to decentralized telescope decision making through an exact potential game. Finally, we used the Shapley value to turn the final schedule into a scientifically interpretable credit allocation.

The same backbone scales conceptually from the toy examples shown here to real public gravitational-wave sky maps and multi-observatory follow-up campaigns.

24 Further reading

For students, the following topics are natural next steps:

  • HEALPix sky maps and Bayesian localization in gravitational-wave astronomy.
  • Greedy algorithms for monotone submodular maximization.
  • Exact potential games and best-response dynamics.
  • Cooperative game theory and Shapley-value credit allocation.
  • Multimessenger astronomy and electromagnetic follow-up scheduling.