subgame-cyerwar-example.R

library(ompr)
library(ompr.roi)
library(ROI.plugin.glpk)
library(dplyr)

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

# Payoff matrices (lower bound for robust case)
Amin <- matrix(c(1, 2, 2, 1), nrow = 2, byrow = TRUE)
Bmin <- matrix(c(3, 0, 4, 1), nrow = 2, byrow = TRUE)

# Store found equilibria as (row, col) pairs
equilibria <- list()
iteration <- 1

repeat {
  cat("Iteration", iteration, ":\n")
  
  model <- MIPModel() %>%
    add_variable(x[i], i = 1:2, type = "binary") %>%
    add_variable(y[j], j = 1:2, type = "binary") %>%
    add_variable(z[i, j], i = 1:2, j = 1:2, lb = 0, ub = 1) %>%
    add_variable(w, lb = -1e6, ub = 1e6) %>%
    
    # Strategy selection
    add_constraint(sum_expr(x[i], i = 1:2) == 1) %>%
    add_constraint(sum_expr(y[j], j = 1:2) == 1) %>%
    
    # Product linearization
    add_constraint(z[i, j] <= x[i], i = 1:2, j = 1:2) %>%
    add_constraint(z[i, j] <= y[j], i = 1:2, j = 1:2) %>%
    add_constraint(z[i, j] >= x[i] + y[j] - 1, i = 1:2, j = 1:2) %>%
    
    # Player 1 robust incentive
    add_constraint(
      sum_expr(z[ii, jj] * Amin[ii, jj], ii = 1:2, jj = 1:2) >=
        sum_expr(y[jj] * Amin[ip, jj], jj = 1:2),
      ip = 1:2
    ) %>%
    
    # Player 2 robust incentive
    add_constraint(
      sum_expr(z[ii, jj] * Bmin[ii, jj], ii = 1:2, jj = 1:2) >=
        sum_expr(x[ii] * Bmin[ii, jp], ii = 1:2),
      jp = 1:2
    ) %>%
    
    set_objective(w, "max")
  
  # Exclude previously found equilibria
  if (length(equilibria) > 0) {
    for (eq in equilibria) {
      row_idx <- eq$row
      col_idx <- eq$col
      # Exclude this exact (row, col) combination
      model <- model %>%
        add_constraint(x[row_idx] + y[col_idx] <= 1)
    }
  }
  
  # Solve
  result <- solve_model(model, with_ROI(solver = "glpk"))
  
  if (result$status == "success") {
    # Extract solution
    x_sol <- get_solution(result, x[i])$value
    y_sol <- get_solution(result, y[j])$value
    
    row_chosen <- which(x_sol == 1)
    col_chosen <- which(y_sol == 1)
    
    # Store this equilibrium
    equilibria[[length(equilibria) + 1]] <- list(row = row_chosen, col = col_chosen)
    
    cat("  Found pure robust Nash equilibrium:\n")
    cat("    Player 1 chooses row:", row_chosen, "\n")
    cat("    Player 2 chooses col:", col_chosen, "\n\n")
    
    iteration <- iteration + 1
  } else {
    cat("No more equilibria found.\n")
    break
  }
}

## Iteration 1 :
##   Found pure robust Nash equilibrium:
##     Player 1 chooses row: 2 
##     Player 2 chooses col: 1 
## 
## Iteration 2 :
## No more equilibria found.

# Summary
cat("\n=== All Pure Robust Nash Equilibria ===\n")

## 
## === All Pure Robust Nash Equilibria ===

for (i in seq_along(equilibria)) {
  cat("Equilibrium", i, ": (row", equilibria[[i]]$row, ", col", equilibria[[i]]$col, ")\n")
}

## Equilibrium 1 : (row 2 , col 1 )